"""Published star-formation histories of the solar neighbourhood, transcribed from the papers' tables (values checked against the arXiv/ar5iv text on 2026-10-11), and their implied head-start statistics. Each SFH is a list of (age_lo, age_hi, mass formed in the bin) in arbitrary units; a uniform rate is assumed inside each bin. `kind` says whether the SFH is a column (surface-density) quantity for the solar cylinder or a mid-plane volume quantity. Sources ------- Mor+ 2019, A&A 624, L1, Table 1 (fiducial G12NP-S): surface mass density of stars formed per age interval, thin disc 0-10 Gyr. The thick disc is not part of the fitted SFH (only its local volume densities at ~10 and ~12 Gyr are fitted), so two variants are given: thin disc only, and thin + an approximate thick disc (see `mor_thick`). Mazzi+ 2024, MNRAS 527, 583, Table 4: vertically integrated SFR per unit area (|z| < 1.3 kpc) in bins of log(age). del Alcazar-Julia+ 2025, A&A (arXiv:2501.17236), Table 2: surface SFR per age bin for the thin disc (0-10 Gyr) and thick disc (8-12 Gyr), four runs (PARSEC or STAREVOL stellar models; thick disc fixed or fitted). Alzate+ 2021 (MNRAS 501, 302) as parametrised by Madau 2023 (ApJ 957, 66, Fig. 2): four Gaussians, stars within 100 pc (mid-plane volume sample). Gallart+ 2024 (A&A 687, A168): only quantiles are quoted in the text (50% of the mass formed within 100 pc was in place 4.5 Gyr ago), so it enters as a single number (SINGLE_NUMBERS: mid-plane fraction older than the Sun about 0.50), not as an SFH. """ import numpy as np from common import AGE_EDGES, headstart_stats # ---------------------------------------------------------------- Mor+ 2019 MOR19_THIN = [(0, 0.1, 0.16), (0.1, 1, 2.2), (1, 2, 6.5), (2, 3, 8.7), (3, 5, 12.1), (5, 7, 7.7), (7, 8, 5.8), (8, 9, 7.4), (9, 10, 9.8)] # Msun/pc^2 formed per bin # 16th/84th percentile half-widths from the same table (minus, plus) MOR19_THIN_ERR = [(0.04, 0.07), (0.5, 0.4), (1.4, 1.4), (2.1, 2.7), (4.5, 4.1), (4.1, 5.7), (2.9, 7.2), (4.4, 5.7), (5.3, 5.9)] def mor_thick(): """Approximate formed surface density of the Mor+19 thick disc. Fitted local densities: 2.6e-3 (age ~10 Gyr) and 0.6e-3 Msun/pc^3 (~12 Gyr). ASSUMPTIONS (ours, not Mor+19's): sech^2 vertical profiles with scale heights 400 and 800 pc (Robin+ 2014 give 345-535 and ~800 pc), Sigma = 2 h rho0, and a present-to-formed mass ratio of 0.6.""" young = 2 * 400.0 * 2.6e-3 / 0.6 old = 2 * 800.0 * 0.6e-3 / 0.6 return [(9.5, 10.5, young), (11.5, 12.5, old)] # -------------------------------------------------------------- Mazzi+ 2024 # (log t1, log t2, Psi [1e-9 Msun/yr/pc^2], minus, plus) MAZZI24 = [(6.60, 7.10, 6.20, 2.60, 3.00), (7.10, 7.30, 5.70, 2.70, 4.00), (7.30, 7.50, 3.10, 1.50, 2.60), (7.50, 7.70, 2.80, 1.40, 2.20), (7.70, 7.90, 2.80, 1.20, 2.30), (7.90, 8.10, 2.11, 0.82, 1.00), (8.10, 8.30, 2.44, 0.83, 1.09), (8.30, 8.50, 2.31, 0.72, 0.76), (8.50, 8.70, 4.21, 0.54, 0.59), (8.70, 8.90, 4.97, 0.47, 0.54), (8.90, 9.10, 6.81, 0.47, 0.50), (9.10, 9.30, 7.79, 0.43, 0.45), (9.30, 9.50, 12.99, 0.54, 0.53), (9.50, 9.70, 7.09, 0.59, 0.59), (9.70, 9.90, 8.60, 0.48, 0.49), (9.90, 10.10, 11.47, 0.38, 0.39)] def mazzi_bins(values=None): out = [] for k, (l1, l2, psi, _, _) in enumerate(MAZZI24): a1, a2 = 10 ** (l1 - 9), 10 ** (l2 - 9) v = psi if values is None else values[k] out.append((a1, a2, v * (a2 - a1))) return out # ------------------------------------------------ del Alcazar-Julia+ 2025 DAJ_THIN_EDGES = [0, 0.15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] DAJ_THICK_EDGES = [8, 9, 10, 11, 12] DAJ25 = { # run: (thin-disc rates, thick-disc rates), Msun/pc^2/Gyr "G13P-13": ([1.10, 0.55, 2.53, 3.42, 4.70, 4.42, 3.28, 0.74, 3.28, 8.88, 14.08], [1.5, 1.9, 1.9, 1.5]), "G13S-13": ([2.26, 1.69, 7.95, 7.05, 7.61, 5.13, 3.94, 8.85, 10.49, 2.13, 4.16], [1.5, 1.9, 1.9, 1.5]), "G13P-17": ([1.04, 0.51, 2.35, 3.13, 4.72, 3.68, 3.84, 0.97, 4.39, 9.07, 13.09], [1.48, 0.69, 0.85, 1.48]), "G13S-17": ([2.06, 1.35, 6.00, 5.01, 5.55, 4.48, 2.30, 6.25, 5.74, 1.50, 4.14], [14.56, 0.63, 0.81, 2.68]), } def daj_bins(run): thin, thick = DAJ25[run] out = [(DAJ_THIN_EDGES[k], DAJ_THIN_EDGES[k + 1], r * (DAJ_THIN_EDGES[k + 1] - DAJ_THIN_EDGES[k])) for k, r in enumerate(thin)] out += [(DAJ_THICK_EDGES[k], DAJ_THICK_EDGES[k + 1], r) for k, r in enumerate(thick)] return out # ------------------------------------------- Alzate+ 2021 / Madau 2023 fit def alzate_madau_phi(edges=AGE_EDGES): cen = np.array([10.0, 5.5, 2.0, 0.7]) sig = np.array([1.2, 0.9, 0.35, 0.3]) peak = np.array([1.0, 1.36, 0.71, 0.86]) from scipy.stats import norm cdf = (peak * sig)[None, :] * norm.cdf(edges[:, None], cen[None, :], sig[None, :]) phi = np.diff(cdf.sum(axis=1)) return phi / phi.sum() # ------------------------------------------------------------------ helpers def bins_to_phi(bins, edges=AGE_EDGES): """Rebin a piecewise-uniform SFH onto the common age grid.""" phi = np.zeros(len(edges) - 1) lo, hi = edges[:-1], edges[1:] for a1, a2, m in bins: overlap = np.clip(np.minimum(hi, a2) - np.maximum(lo, a1), 0, None) phi += m * overlap / (a2 - a1) return phi / phi.sum() def all_published(): """dict name -> (kind, phi on the common grid).""" out = { "Mor+19 thin disc only": ("column", bins_to_phi(MOR19_THIN)), "Mor+19 thin + approx. thick": ("column", bins_to_phi(MOR19_THIN + mor_thick())), "Mazzi+24": ("column", bins_to_phi(mazzi_bins())), } for run in DAJ25: out[f"del Alcazar-Julia+25 {run}"] = ("column", bins_to_phi(daj_bins(run))) out["Alzate+21 (Madau 23 fit), 100 pc"] = ("midplane", alzate_madau_phi()) return out # Results that are published only as single numbers (no tabulated SFH), and # the model predictions the measurement is compared with. value = fraction # older than the Sun (or the closest quoted equivalent). SINGLE_NUMBERS = [ dict(name="Gallart+24 (ChronoGal I), 100 pc", kind="midplane", f_gt_0=0.50, note="A&A 687, A168: '50% ... of the stellar mass was reached 4.5 " "... Gyr ago' (mass formed, bubble of 100 pc)"), ] MODELS = [ dict(name="Lineweaver 2001 (cosmic SFH + metallicity)", f=0.74, err=0.09, note="Icarus 151, 307: 74 +/- 9% of Earth-like planets older than " "Earth, on average 1.8 +/- 0.9 Gyr older"), dict(name="Lineweaver+ 2004 (Galactic habitable zone)", f=0.75, err=None, note="Science 303, 59"), dict(name="Legassick 2015 (chemical-evolution model)", f=0.77, err=None, note="arXiv:1509.02832: 77% of potentially habitable systems are on " "average 3.13 Gyr older than the Sun"), dict(name="Behroozi & Peeples 2015 (Milky Way-mass halo)", f=0.80, err=None, note="MNRAS 454, 1811"), dict(name="Madau 2023 (100 pc SFH, K-dwarf planets)", f=0.77, err=None, note="ApJ 957, 66: mid-plane SFH of Alzate+21"), ] def young_fraction(phi, t=2.0, edges=AGE_EDGES): lo, hi = edges[:-1], edges[1:] part = np.clip((t - lo) / (hi - lo), 0, 1) return float((phi * part).sum() / phi.sum()) def table(): import pandas as pd rows = [] for name, (kind, phi) in all_published().items(): st = headstart_stats(phi) st.update(name=name, kind=kind, young_frac_lt2=young_fraction(phi), frac_gt8=float(1 - young_fraction(phi, 8.0))) rows.append(st) for r in SINGLE_NUMBERS: # quantile-only results: one number each rows.append(dict(name=r["name"], kind=r["kind"] + " (single number)", f_gt_0=r["f_gt_0"])) return pd.DataFrame(rows).set_index("name") def mazzi_mc(n=2000, seed=3): """Propagate Mazzi+24's quoted bin errors (split normal) to f_older.""" rng = np.random.default_rng(seed) psi = np.array([m[2] for m in MAZZI24]) em = np.array([m[3] for m in MAZZI24]) ep = np.array([m[4] for m in MAZZI24]) out = [] for _ in range(n): z = rng.standard_normal(len(psi)) v = np.clip(psi + np.where(z > 0, ep, em) * z, 0, None) out.append(headstart_stats(bins_to_phi(mazzi_bins(v)))) return out if __name__ == "__main__": import pandas as pd pd.set_option("display.width", 250, "display.max_columns", 30) print(table().round(3).to_string()) mc = mazzi_mc() f = np.array([m["f_gt_0"] for m in mc]) print("Mazzi+24 f_older with quoted errors: %.3f +%.3f -%.3f" % (np.median(f), np.percentile(f, 84) - np.median(f), np.median(f) - np.percentile(f, 16)))