#!/usr/bin/env python3 """Assemble report.md from the tables in data/ (last step of run_all.py). Every number in the report is read from a CSV written by the analysis scripts, so the text cannot drift away from the code. The prose around the numbers is fixed; if a rerun changes a number materially, re-read the prose. """ import os import datetime import numpy as np import pandas as pd import common as C import fermi_clock as F import published_sfh as P ROOT = os.path.dirname(C.DATA) OUT = os.path.join(ROOT, "report.md") # Raw catalogue sizes (rows downloaded). Counted when the files were # fetched and confirmed by the referee's rerun; "used" counts are read from # the cut-flow tables at run time. RAW_ROWS = { "gcs": 16682, "flame_sunlike": 227504, "flame_hip": 53443, "starhorse": 291807, } def rd(name, index=None): df = pd.read_csv(os.path.join(C.DATA, name)) return df.set_index(index) if index else df def table(headers, rows, align=None): align = align or ["l"] + ["r"] * (len(headers) - 1) out = ["| " + " | ".join(headers) + " |", "|" + "|".join(":--" if a == "l" else "--:" for a in align) + "|"] for r in rows: out.append("| " + " | ".join(str(c) for c in r) + " |") # a literal "|" inside a cell would end the cell esc = [out[0].replace("|Z|", "\\|Z\\|"), out[1]] for line in out[2:]: cells = line[2:-2].split(" | ") esc.append("| " + " | ".join(c.replace("|", "\\|") for c in cells) + " |") return "\n".join(esc) def f3(v): return "–" if v is None or not np.isfinite(v) else f"{v:.3f}" def f2(v): return "–" if v is None or not np.isfinite(v) else f"{v:.2f}" def rng(lo, hi, fmt="{:.2f}"): return fmt.format(lo) + "–" + fmt.format(hi) def sci(v, d=1): if not np.isfinite(v): return "–" e = int(np.floor(np.log10(abs(v)))) return f"{v / 10 ** e:.{d}f}e{e:d}" def n_(v): return f"{int(v):,}" def cut_rows(cf): return [(r.step, n_(r.n), r.reason if isinstance(r.reason, str) else "") for r in cf.itertuples()] def main(runtime_s=None, quick=False): xr = rd("xr_results.csv", "variant") nat = rd("nataf_results.csv", "variant") gs = rd("gaia_sg_results.csv", "variant") est = rd("headline_estimates.csv") ad = rd("headline_adopted.csv").iloc[0] sysd = rd("xr_systematics.csv") off = rd("age_scale_offsets.csv") swap = rd("age_swap_test.csv") lt = rd("xr_layer_test.csv") nl = rd("nataf_layer_test.csv") hz = rd("xr_scale_heights.csv") bl = rd("xr_bright_limit.csv") cfx = rd("xr_cutflow.csv") cfn = rd("nataf_cutflow.csv") cfg = rd("gaia_sg_cutflow.csv") gq = rd("gaia_sg_quality.csv") shc = rd("gaia_sg_vs_starhorse.csv", "subset") kin = rd("nataf_kinematics.csv") pub = rd("published_sfh_table.csv", "name") ms = rd("ms_census.csv", "sample") twins = rd("solar_twins.csv", "sample") fvg = rd("flame_vs_gcs.csv", "subset") rngc = rd("fermi_clock_range.csv", "v_over_c") marg = rd("fermi_marginal_limit.csv") clock = rd("fermi_clock.csv") rank = rd("sun_rank.csv", "distribution") gal = rd("galaxy_wide_bracket.csv") asm = rd("fermi_assumptions.csv", "quantity") dists = F.distributions() B = xr.loc["formed: all [Fe/H]"] # primary estimate NP = nat.loc["formed: all [Fe/H]"] NX = nat.loc["ages on the X&R22 scale"] NA = nat.loc["primary + annex"] N3 = nat.loc["tracer ages > 3 Gyr (young prior for 3 Gyr)"] G0 = gs.loc["formed: all [M/H]"] GX = gs.loc["ages on the X&R22 scale: formed, all [M/H]"] GXm = gs.loc["ages on the X&R22 scale: mid-plane"] GN = gs.loc["ages on the Nataf+24 scale: formed, all [M/H]"] def xcut(key): return int(cfx[cfx.step.str.startswith(key)].n.iloc[0]) def ncut(key): return int(cfn[cfn.step.str.startswith(key)].n.iloc[0]) def gcut(key): return int(cfg[cfg.step.str.contains(key, regex=False)].n.iloc[0]) def ratio(pair, col="median_ratio"): t = off[(off.pair == pair) & (off.subset == "all")] return float(t[col].iloc[0]), int(t.n.iloc[0]) r_nat, n_nat = ratio("Nataf+24 vs X&R22") r_sh, n_sh = ratio("StarHorse vs X&R22") r_fl, n_fl = ratio("Gaia FLAME vs X&R22") r_gx, n_gx = ratio("this study's Gaia-only ages vs X&R22") d_gx = ratio("this study's Gaia-only ages vs X&R22", "median_diff_gyr")[0] r_gn, n_gn = ratio("this study's Gaia-only ages vs Nataf+24") r_shn, n_shn = ratio("StarHorse vs Nataf+24") r_gs, n_gs = ratio("this study's Gaia-only ages vs StarHorse") d_gs = ratio("this study's Gaia-only ages vs StarHorse", "median_diff_gyr")[0] g_lo, g_hi = min(r_gx, r_gn, r_gs), max(r_gx, r_gn, r_gs) d_nat = ratio("Nataf+24 vs X&R22", "median_diff_gyr")[0] d_sh = ratio("StarHorse vs X&R22", "median_diff_gyr")[0] d_fl = ratio("Gaia FLAME vs X&R22", "median_diff_gyr")[0] sw = swap[swap.subset == "in common, Nataf quality flags"].set_index("ages") sw_x = sw.loc["X&R ages"] sw_n = sw.loc["Nataf+24 ages"] sw_r = sw[sw.index.str.startswith("Nataf+24 ages /")].iloc[0] d_swap = sw_n.f_older_column - sw_x.f_older_column e = est.set_index("estimate") inr = est[est.in_range] pub_col = pub[pub.kind == "column"] plo, phi_ = pub_col.f_gt_0.min(), pub_col.f_gt_0.max() mazzi = pub.loc["Mazzi+24"] alz = pub.loc["Alzate+21 (Madau 23 fit), 100 pc"] # ---------------------------------------------------------- Fermi clock c_slow, c_fast = rngc.loc[1e-4], rngc.loc[1e-1] kmin, kmax = marg.factor.min(), marg.factor.max() pm_lo, pm_hi = marg["p95_at_0.1c"].min(), marg["p95_at_0.1c"].max() pm_lo_s, pm_hi_s = marg["p95_at_1e-4c"].min(), marg["p95_at_1e-4c"].max() cf = clock[(clock.distribution == F.FID) & (clock.model == "Earth clock")] p_ext_lo = cf[cf.v_over_c == 1e-1].p95_low.min() p_ext_hi = cf[cf.v_over_c == 1e-4].p95_high.max() n_sun = float(asm.loc["N Sun-like stars ever formed", "value"]) kfrac = float(asm.loc["Kroupa IMF number fraction 0.8-1.2 Msun", "value"]) hard = clock[(clock.distribution == F.FID) & clock.eta_earth.str.contains("lower") & (clock.v_over_c == 1e-3)].set_index("model") rk = rank.loc[F.FID] phi_fid = dists[F.FID][1] t_voy = F.D_GAL_LY / (17.0 / 299792.458) / 1e9 f_voy = C.frac_older(phi_fid, t_voy) # ----------------------------------------- mid-plane to column corrections mid = xr.loc["no column weights (mid-plane slab)"] corr = [ ("LAMOST subgiants, dynamical weights", B.f_gt_0 - mid.f_gt_0), ("LAMOST subgiants, star-count heights (exp / sech² / thick)", None), ("Nataf+24 primary sample, dynamical weights", NP.f_gt_0 - nat.loc["no column weights (mid-plane slab)", "f_gt_0"]), ("Nataf+24 primary + annex, dynamical weights", NA.f_gt_0 - nat.loc["primary + annex, mid-plane", "f_gt_0"]), ("Gaia 200 pc subgiants, ages as computed", G0.f_gt_0 - gs.loc["no column weights (mid-plane sphere)", "f_gt_0"]), ("Gaia 200 pc subgiants, X&R age scale", GX.f_gt_0 - GXm.f_gt_0), ("Gaia FLAME dwarfs, 200 pc", ms.loc["Gaia FLAME d<200 pc, 0.8-1.2 Msun, column", "deconv_f_gt_0"] - ms.loc["Gaia FLAME d<200 pc, 0.8-1.2 Msun, volume", "deconv_f_gt_0"]), ("GCS dwarfs, 40 pc (Padova / BaSTI ages)", None), ] sc = [xr.loc[k, "f_gt_0"] - mid.f_gt_0 for k in ("star-count heights (Mazzi+24 exp)", "star-count heights (Mazzi+24 sech2)", "star-count heights (exp, thick disc 600/900 pc)")] gcs = [ms.loc[f"GCS d<40 pc, 0.9-1.2 Msun, {m}, column", "deconv_f_gt_0"] - ms.loc[f"GCS d<40 pc, 0.9-1.2 Msun, {m}, volume", "deconv_f_gt_0"] for m in ("Padova", "BaSTI")] corr_rows, corr_vals = [], [] for lab, v in corr: if "star-count" in lab: corr_rows.append((lab, " / ".join(f"+{x:.2f}" for x in sc))) corr_vals += sc elif lab.startswith("GCS"): corr_rows.append((lab, " / ".join(f"+{x:.2f}" for x in gcs))) else: corr_rows.append((lab, f"+{v:.2f}")) corr_vals.append(v) corr_lo, corr_hi = min(corr_vals), max(corr_vals) # ------------------------------------------------------------- layer test def layer_rows(t, models): rows = [] for (lo, hi), g in t[t.layer_lo > 0].groupby(["layer_lo", "layer_hi"]): g = g.set_index("model") rows.append([f"{int(lo)}–{int(hi)} pc"] + [f"{g.loc[m, 'f_older']:.3f} / {g.loc[m, 'f_gt8']:.3f}" if m in g.index else "–" for m in models]) return rows def layer_dev(t, model): g = t[t.layer_lo > 0].pivot(index="layer_lo", columns="model", values="f_older") d = g[model] - g["observed"] return float(d.min()), float(d.max()) dx = layer_dev(lt, "dynamical weights") dsc = layer_dev(lt, "star counts (exp)") dth = layer_dev(lt, "star counts (exp_thick)") dn = layer_dev(nl, "dynamical weights") zero = cfx[cfx.step.str.contains("zero vertical weight")].iloc[0] slab_mode = xr.loc["slab-fraction weights (uniform S(z))", "f_gt_0"] bl23 = bl.iloc[0] blw = xr.loc["2-3 Gyr tracers up-weighted for bright-limit loss", "f_gt_0"] qall = gq.iloc[-1] q_rv = gq[gq.age_lo >= 2].iloc[:-1] kp = kin[kin["sample"] == "primary"] ka = kin[kin["sample"] == "primary + annex"] def srow(key): return sysd[sysd.source.str.startswith(key)].iloc[0] s_age_pub = srow("age scale: other published") s_age_ref = srow("age scale: this study's Gaia-only") s_sc = srow("vertical model: star-count") s_col = srow("vertical model: column counted") s_young = srow("young fraction") K13 = "phi_young = 0.13 (tracers' own flat-rate value)" sc_all = [xr.loc[k, "f_gt_0"] for k in ("star-count heights (Mazzi+24 exp)", "star-count heights (Mazzi+24 sech2)", "star-count heights (exp, thick disc 600/900 pc)")] n_sc_in = sum(ad.f_range_lo <= x <= ad.f_range_hi for x in sc_all) sc_in = {3: "all three", 2: "two of the three", 1: "one of the three", 0: "none"}[n_sc_in] # LAMOST minus Gaia 200 pc on each of three common age scales dg = [B.f_gt_0 - GX.f_gt_0, xr.loc["ages on the Nataf+24 scale", "f_gt_0"] - GN.f_gt_0, xr.loc["ages x1.146 (referee: Gaia-only scale)", "f_gt_0"] - G0.f_gt_0] uv = max(abs(nat.loc[k, "f_gt_0"] - NA.f_gt_0) for k in ("primary + annex, UV ramp -0.6..-0.35", "primary + annex, UV ramp -0.4..-0.15")) when = datetime.date.today().isoformat() nb = "20 / 10 / 8" if quick else "100 / 50 / 30" L = [] A = L.append # =================================================================== text A("# Who had a head start? An empirical age census of Sun-like stars " "and the size of the Fermi clock") A("") A(f"*Study `headstart`, revised after referee review. This file is " f"written by `code/make_report.py`, the last step of " f"`code/run_all.py`; every number is read from the tables in `data/` " f"at run time (run of {when}" + (f", {runtime_s / 60:.0f} min" if runtime_s else "") + (", QUICK mode: fewer bootstrap resamples" if quick else "") + ").*") A("") A(f"**Result.** About {100 * ad.f_adopted:.0f}% of the Sun-like " f"(0.8–1.2 Msun) stars now residing in the solar cylinder formed " f"before the Sun. The primary estimate, from {n_(B.n)} LAMOST subgiants " f"with Xiang & Rix (2022) ages, is {B.f_gt_0:.3f} " f"({rng(B.f_gt_0_lo, B.f_gt_0_hi, '{:.3f}')}, 16–84%, statistical plus " f"the prior on the young fraction). That interval is not the " f"uncertainty. The two catalogue samples whose selection can be " f"modelled (LAMOST and the Nataf+24 primary sample), on four " f"published age scales, span {rng(ad.f_core_lo, ad.f_core_hi)} " f"({int(ad.n_core)} estimates); two weaker samples that cannot be " f"ruled out widen this to {rng(ad.f_range_lo, ad.f_range_hi)}. The " f"result is therefore **{ad.f_adopted:.2f}, systematic range " f"{rng(ad.f_range_lo, ad.f_range_hi)}**; published star-formation " f"histories re-integrated here give {rng(plo, phi_)}. " f"This is consistent with Lineweaver's (2001) model value of 74 ± 9% " f"and does not test it more sharply than that. Among the older stars " f"the mean head start is {B.mean_headstart_older:.1f} Gyr in the " f"primary sample ({rng(ad.mean_hs_lo, ad.mean_hs_hi, '{:.1f}')} Gyr " f"across all of these estimates), and " f"{100 * ad.f5_range_lo:.0f}–{100 * ad.f5_range_hi:.0f}% of all " f"Sun-like stars are more than 5 Gyr older than the Sun. With a floor " f"of [M/H] > −0.5 the fraction is " f"{xr.loc['formed: [M/H]>-0.5', 'f_gt_0']:.2f}, with [Fe/H] > −0.3 " f"{xr.loc['formed: [Fe/H]>-0.3', 'f_gt_0']:.2f}; these are lower " f"bounds, because small planets show no sharp metallicity floor. " f"Main caveats: the number depends on the absolute stellar age scale " f"(catalogues differ by {abs(1 - r_sh) * 100:.0f}–" f"{abs(1 - r_fl) * 100:.0f}% for the same stars) and on the weight the " f"thick disc receives in the vertical column; stars younger than 2 Gyr " f"are not traced and their share (0.10 ± 0.04) is imported; and the " f"statement is about the solar cylinder, not the Galaxy. For the Fermi " f"question none of this matters much. If no civilisation has settled " f"the Galaxy, and a settled Galaxy would be evident, then with a " f"fiducial {sci(c_slow.n_eff_fiducial)}–{sci(c_fast.n_eff_fiducial)} " f"habitable-zone rocky planets holding enough head start the " f"per-planet probability of producing a settling civilisation is below " f"{sci(c_fast.p95_fiducial)}–{sci(c_slow.p95_fiducial)} (95%) for that " f"fiducial count, and below {sci(pm_lo)}–{sci(pm_hi)} once the " f"uncertainty of the star count and of eta_Earth is marginalised.") A("") A("## 1. Question") A("") A("Fermi's argument needs older planets: other Sun-like stars must have " "had time to produce a civilisation and for it to cross the Galaxy " "before now. Lineweaver (2001) derived from a cosmic star-formation " "and metallicity model that 74 ± 9% of terrestrial planets are older " "than Earth, on average by 1.8 ± 0.9 Gyr, and that number has been " "repeated ever since. The question here is what measured stellar ages " "say: what fraction of Sun-like stars is older than the Sun, by how " "much, how sure is that, and how large a \"Fermi clock\" does it give " "(the number of planets with enough head start for a settlement wave " "of a given speed to have finished)?") A("") A("## 2. Data") A("") A(table(["Source", "Rows", "Used", "Role"], [ ("Xiang & Rix 2022, LAMOST subgiants (VizieR J/other/Nat/603.599)", n_(xcut("Xiang")), n_(xcut("J_z > 0")), "primary tracer sample"), ("Nataf+ 2024 subgiants (`TotalSample.fits`, authors' repository; " "ApJ 976, 87) with Gaia DR3 astrometry and radial velocities " "looked up by source_id", f"{n_(ncut('Nataf+ 2024 catalogue'))}; {n_(ncut('pre-selected'))} " "looked up", n_(ncut("age >= 2 Gyr")), "second tracer sample (added in revision)"), ("Gaia DR3 within 200 pc (gaia_source + astrophysical_parameters)", n_(gcut("parent download")), n_(gcut("sample used")), "third tracer sample, this study's own photometric ages " "(cross-check)"), ("Gaia DR3 FLAME ages, Sun-like dwarfs within 200 pc", n_(RAW_ROWS["flame_sunlike"]), n_(ms.loc["Gaia FLAME d<200 pc, 0.8-1.2 Msun, volume", "n"]), "dwarf cross-check"), ("Geneva-Copenhagen Survey, Casagrande+ 2011 (J/A+A/530/A138)", n_(RAW_ROWS["gcs"]), n_(ms.loc["GCS d<40 pc, 0.9-1.2 Msun, Padova, volume", "n"]), "dwarf cross-check (40 pc)"), ("StarHorse ages for the Gaia RVS sample, Queiroz+ 2023 " "(J/A+A/673/A155), local subset", n_(RAW_ROWS["starhorse"]), f"{n_(n_sh)} / {n_(n_shn)} / {n_(shc.loc['all', 'n'])} in common", "age-scale comparison"), ("Solar twins: Spina+ 2018, Nissen+ 2020, Taniguchi+ 2026", "79 / 72 / 6,594", "all", "raw fractions, cross-check"), ("Published star-formation histories: Mor+ 2019, Mazzi+ 2024, del " "Alcázar-Julià+ 2025 (four runs), Alzate+ 2021 as fitted by Madau " "2023; Gallart+ 2024 (one number)", "transcribed tables", "–", "comparison"), ("MIST v1.2 isochrones", "–", "–", "tracer weights, Gaia-only ages"), ("Bryson+ 2021 eta_Earth", "literature value", "–", "Fermi clock"), ], align=["l", "r", "r", "l"])) A("") A("**Selection cuts, LAMOST sample** (`data/xr_cutflow.csv`)") A("") A(table(["Step", "Stars left", "Reason"], cut_rows(cfx), align=["l", "r", "l"])) A("") A("**Selection cuts, Nataf+24 sample** (`data/nataf_cutflow.csv`)") A("") A(table(["Step", "Stars left", "Reason"], cut_rows(cfn), align=["l", "r", "l"])) A("") A("**Selection cuts, Gaia 200 pc sample** (`data/gaia_sg_cutflow.csv`)") A("") A(table(["Step", "Stars left", "Reason"], cut_rows(cfg), align=["l", "r", "l"])) A("") A(f"Two losses in these tables depend on age and are quantified. (i) " f"{n_(zero.n)} LAMOST stars ({100 * zero.n / B.n:.1f}%) receive zero " f"column weight: {zero.reason}. A weighting that has no zeros gives " f"{slab_mode:.3f} instead of {B.f_gt_0:.3f}. (ii) LAMOST's bright " f"limit removes the nearest 2–3 Gyr subgiants: " f"{100 * bl23.frac_beyond_0p55:.0f}% of them lie beyond 0.55 kpc " f"against {100 * bl.frac_beyond_0p55.iloc[1:].mean():.0f}% of older " f"tracers, an implied completeness of {bl23.implied_completeness:.2f}; " f"up-weighting them gives {blw:.3f} ({blw - B.f_gt_0:+.3f}). In the " f"Gaia sample {100 * (1 - qall.rv_fraction):.0f}% of the " f"membership-weighted candidates have no radial velocity; the loss is " f"flat with age above 2 Gyr ({rng(q_rv.rv_fraction.min(), q_rv.rv_fraction.max())} " f"retained in every bin). {100 * qall.ruwe_gt_1p4_fraction:.0f}% of " f"that sample has RUWE > 1.4 " f"({100 * q_rv.ruwe_gt_1p4_fraction.min():.0f}–" f"{100 * q_rv.ruwe_gt_1p4_fraction.max():.0f}% by age bin), a sign of " f"heavy contamination by unresolved binaries.") A("") A("## 3. Method") A("") A("**Subgiants as tracers.** A star's subgiant phase is brief and its " "age is well determined there. A subgiant observed inside a " "colour-magnitude selection box stands for I0/G(age, [M/H]) Sun-like " "stars born at the same time, where G is the IMF-weighted mass " "interval of a MIST isochrone that lies inside the box and I0 the same " "integral over 0.8–1.2 Msun (Kroupa slope 2.3). This counts stars " "ever formed; multiplying by the fraction still on the main sequence " "gives the living population. The " "survey's sampling rate cancels as long as it does not depend on age, " "which is tested by comparing distance shells (Fig. 3a): the age mix " "is flat between 0.3 and 0.8 kpc and that shell is used, within 0.5 kpc " "of the solar radius.") A("") A("**From mid-plane to column.** Old stars spend most of their time far " "from the plane, so a thin slab under-counts them. Each star in the " "slab |Z| < 200 pc is weighted by the inverse of the fraction of its " "vertical orbit spent inside the slab, |v_z| T(E) / (4 z0), in a " "one-dimensional Milky Way potential.") A("") A(f"**Age errors, young stars, uncertainties.** The weighted age " f"histogram (0.25 Gyr bins) is deconvolved for the catalogue age " f"errors by expectation-maximisation. No 0.8–1.2 Msun star younger " f"than 2 Gyr is a subgiant, so the share of stars formed in the last " f"2 Gyr is set to 0.10 ± 0.04, the mean and scatter of seven published " f"column star-formation histories. Statistical intervals come from " f"bootstrap resampling ({nb} resamples for the three samples) with " f"the young share redrawn each time.") A("") A(f"**Three samples.** The Nataf+24 primary sample is treated the same " f"way with its own selection box and Gaia kinematics. Its metal-poor " f"annex, selected with ultraviolet colours that the isochrones cannot " f"reproduce exactly, is modelled approximately and reported " f"separately. The Gaia 200 pc sample uses this study's own ages from " f"(M_G, BP−RP, [M/H]) and subgiant membership probabilities. Ages of " f"stars shared between catalogues are compared directly " f"(`code/age_scales.py`), and the LAMOST estimate is recomputed with " f"Nataf+24 ages substituted star by star.") A("") A(f"**Fermi clock.** N_eff(v) = N_Sun-like × eta_Earth × the fraction " f"of stars whose age minus the settling time (10^5 light-years / v) " f"exceeds the emergence time. Assumed, not measured: 2e11 stars in the " f"Galaxy (factor 2 either way), of which a Kroupa IMF puts " f"{100 * kfrac:.1f}% at 0.8–1.2 Msun ({sci(n_sun, 2)} stars); " f"eta_Earth = 0.37 (+0.48/−0.21) from Bryson+ 2021; emergence at " f"Earth's present age (\"Earth clock\"), or hard-step power laws; the " f"solar-cylinder age mix applied to the whole Galaxy.") A("") A("## 4. Results") A("") A("### 4.1 Fraction of Sun-like stars older than the Sun") A("") rows = [] for r in est.itertuples(): rows.append((r.estimate, f3(r.f_older), rng(r.stat_lo, r.stat_hi, "{:.3f}") if np.isfinite(r.stat_lo) else "–", f2(r.mean_headstart_older), f3(r.f_gt_5), {"core": "core", "bracket": "bracket", "excluded": "not used"}[r.tier], r.comment if isinstance(r.comment, str) else "")) A(table(["Estimate (stars formed, full column)", "f older", "16–84%", "mean head start of the older [Gyr]", "f > 5 Gyr older", "tier", "note"], rows, align=["l", "r", "r", "r", "r", "l", "l"])) A("") A(f"**Core** estimates use individual published subgiant ages and a " f"selection box that can be evaluated on isochrones: LAMOST and the " f"Nataf+24 primary sample, each on the four published age scales, the " f"star-by-star age swap, and Nataf+24 without its incompletely " f"sampled 2–3 Gyr tracers. They span " f"{rng(ad.f_core_lo, ad.f_core_hi, '{:.3f}')} and the adopted value, " f"{ad.f_adopted:.2f}, is the midpoint. On one age scale the two " f"samples agree ({B.f_gt_0:.3f} against {NX.f_gt_0:.3f}). " f"**Bracketing** estimates come from samples that are weaker but not " f"refuted. The Gaia 200 pc sample, with this study's own ages rescaled " f"by a constant onto published scales, gives " f"{rng(min(GX.f_gt_0, GN.f_gt_0), max(GX.f_gt_0, GN.f_gt_0, gs.loc['ages on the StarHorse scale: formed, all [M/H]', 'f_gt_0']), '{:.3f}')}, " f"lower than LAMOST by {rng(min(dg), max(dg))} on each common scale; " f"the cause is not established (home-made ages, a high binary " f"fraction and a smaller volume are candidates). The Nataf+24 sample " f"including its metal-poor annex gives " f"{rng(NA.f_gt_0, nat.loc['primary + annex, ages on the X&R22 scale', 'f_gt_0'], '{:.3f}')}; " f"the primary box under-samples metal-poor stars and the annex " f"appears contaminated (section 5), so the two bracket the truth. " f"With the bracketing estimates the **full range is " f"{rng(ad.f_range_lo, ad.f_range_hi, '{:.3f}')}**. Published column " f"star-formation histories give {rng(plo, phi_)} " f"(Mazzi+24: {mazzi.f_gt_0:.3f}); the two published mid-plane values " f"are 0.50 (Gallart+ 2024) and {alz.f_gt_0:.2f} (Alzate+ 2021 as " f"fitted by Madau 2023, who quotes 77%). See Fig. 2.") A("") A("### 4.2 How large is the head start?") A("") def th(col): return (f"{B[col]:.3f}", rng(B[col + '_lo'], B[col + '_hi'], "{:.3f}")) rows = [("older than the Sun",) + th("f_gt_0") + (rng(ad.f_range_lo, ad.f_range_hi),), ("older by more than 1 Gyr",) + th("f_gt_1") + (rng(ad.f1_range_lo, ad.f1_range_hi),), ("older by more than 3 Gyr",) + th("f_gt_3") + (rng(ad.f3_range_lo, ad.f3_range_hi),), ("older by more than 5 Gyr",) + th("f_gt_5") + (rng(ad.f5_range_lo, ad.f5_range_hi),), ("mean head start of the older stars [Gyr]", f2(B.mean_headstart_older), rng(B.mean_headstart_older_lo, B.mean_headstart_older_hi), rng(ad.mean_hs_lo, ad.mean_hs_hi)), ("median head start of the older stars [Gyr]", f2(B.median_headstart_older), rng(B.median_headstart_older_lo, B.median_headstart_older_hi), rng(ad.median_hs_lo, ad.median_hs_hi)), ("mean age minus the Sun's age, all stars [Gyr]", f2(B.mean_offset_all), rng(B.mean_offset_all_lo, B.mean_offset_all_hi), rng(ad.offset_all_lo, ad.offset_all_hi))] A(table(["Quantity (all metallicities)", "Primary (LAMOST)", "16–84% statistical", "Range, core and bracketing estimates"], rows)) A("") A(f"The last row is the quantity Lineweaver quoted as 1.8 ± 0.9 Gyr. " f"If the disfavoured variants are included (star-count scale heights, " f"ages × 0.87 and × 1.146) the mean head start of the older stars " f"ranges over {rng(min(s_sc.mean_hs_min, s_age_ref.mean_hs_min), s_age_ref.mean_hs_max, '{:.1f}')} Gyr. " f"Counting only stars still on the main sequence, the fraction older " f"than the Sun is {xr.loc['living MS: all [Fe/H]', 'f_gt_0']:.3f} " f"(LAMOST) and {nat.loc['living MS: all [Fe/H]', 'f_gt_0']:.3f} " f"(Nataf+24).") A("") A("### 4.3 Metallicity floors") A("") rows = [] for lab, tab, key in [ ("LAMOST, no floor", xr, "formed: all [Fe/H]"), ("LAMOST, [M/H] > −0.5", xr, "formed: [M/H]>-0.5"), ("LAMOST, [Fe/H] > −0.5", xr, "formed: [Fe/H]>-0.5"), ("LAMOST, [M/H] > −0.3", xr, "formed: [M/H]>-0.3"), ("LAMOST, [Fe/H] > −0.3", xr, "formed: [Fe/H]>-0.3"), ("Nataf+24 primary, [Fe/H] > −0.5", nat, "formed: [Fe/H]>-0.5"), ("Nataf+24 primary, [Fe/H] > −0.3", nat, "formed: [Fe/H]>-0.3"), ("Gaia 200 pc on the X&R scale, [M/H] > −0.5", gs, "ages on the X&R22 scale: [M/H]>-0.5")]: r = tab.loc[key] rows.append((lab, f3(r.f_gt_0), rng(r.f_gt_0_lo, r.f_gt_0_hi, "{:.3f}") if np.isfinite(r.get("f_gt_0_lo", np.nan)) else "–", f2(r.mean_headstart_older), f3(r.f_gt_5), f2(r.frac_pass))) A(table(["Sample and floor", "f older", "16–84%", "mean head start of the older [Gyr]", "f > 5 Gyr older", "share of stars passing"], rows)) A("") A("Old thick-disc stars are alpha-enhanced, so a floor on [Fe/H] " "prunes them harder than a floor on total metals; the [M/H] rows " "correct this. Small planets show no sharp floor at all: Buchhave et " "al. (2012, 2014) find planets below about 1.7 Earth radii around " "stars of widely different metallicity, and Kepler-444, a metal-poor " "thick-disc star 11.2 ± 1.0 Gyr old, has five planets between Mercury " "and Venus in size (Campante et al. 2015). The floor rows are " "therefore lower bounds on the head start.") A("") A("### 4.4 The age scale") A("") rows = [] for pair in ["Nataf+24 vs X&R22", "StarHorse vs X&R22", "Gaia FLAME vs X&R22", "StarHorse vs Nataf+24", "this study's Gaia-only ages vs X&R22", "this study's Gaia-only ages vs Nataf+24", "this study's Gaia-only ages vs StarHorse"]: t = off[(off.pair == pair) & (off.subset == "all")].iloc[0] rows.append((pair, n_(t.n), f"{t.median_diff_gyr:+.2f}", f"{t.median_ratio:.3f}")) A(table(["Ages compared for the same stars", "Stars", "Median difference [Gyr]", "Median ratio"], rows)) A("") A(f"Three published catalogues lie within " f"{100 * abs(1 - r_fl):.0f}% below to {100 * abs(1 - r_sh):.0f}% above " f"the Xiang & Rix scale (Fig. 5a). This study's own Gaia-only ages are " f"{100 * (g_lo - 1):.0f}–{100 * (g_hi - 1):.0f}% older than all three " f"published sets of ages for the same stars and " f"pile {100 * G0.f_last_bin:.0f}% of their weight into the last age " f"bins, so the Gaia sample \"as computed\" ({G0.f_gt_0:.3f}) is biased " f"and its agreement with LAMOST in the first pass was a coincidence of " f"age scale. Rescaling the LAMOST ages by the three published ratios " f"gives {rng(s_age_pub.f_min, s_age_pub.f_max, '{:.3f}')}. A harder " f"test swaps ages star by star: for the {n_(sw_x.n)} LAMOST shell " f"stars that also have Nataf+24 ages, the same estimator with the " f"same weights gives {sw_x.f_older_column:.3f} with Xiang & Rix ages " f"and {sw_n.f_older_column:.3f} with Nataf+24 ages " f"({d_swap:+.3f}; {sw_r.f_older_column:.3f} after dividing the " f"Nataf+24 ages by {r_nat:.2f}). The swap moves more than the constant " f"ratio predicts because {100 * sw_n.weight_share_xr_old_other_young:.0f}% " f"of the weight crosses the Sun's age downward and only " f"{100 * sw_n.weight_share_xr_young_other_old:.1f}% upward. The age " f"scale is the largest single uncertainty.") A("") A("### 4.5 Is the column weighting right? The layer test") A("") A(f"The orbit-based weights imply effective scale heights that match " f"Mazzi+24's star-count fits below 4 Gyr " f"({hz.h_eff_pc.iloc[0]:.0f}, {hz.h_eff_pc.iloc[1]:.0f}, " f"{hz.h_eff_pc.iloc[2]:.0f} pc against " f"{hz.h_mazzi24_exp_pc.iloc[0]:.0f}, {hz.h_mazzi24_exp_pc.iloc[1]:.0f}, " f"{hz.h_mazzi24_exp_pc.iloc[2]:.0f} pc) and exceed them above: by " f"{100 * (hz.h_eff_pc.iloc[3:7] / hz.h_mazzi24_exp_pc.iloc[3:7] - 1).min():.0f}–" f"{100 * (hz.h_eff_pc.iloc[3:7] / hz.h_mazzi24_exp_pc.iloc[3:7] - 1).max():.0f}% " f"between 4 and 10 Gyr and threefold in the oldest bin " f"({hz.h_eff_pc.iloc[7]:.0f} against {hz.h_mazzi24_exp_pc.iloc[7]:.0f} pc; " f"Fig. 3b). The weights also put " f"{100 * B.high_alpha_share_tracers:.0f}% of the column's Sun-like " f"stars in the high-alpha population, against " f"{100 * mid.high_alpha_share_tracers:.0f}% in the mid-plane. Using " f"star-count heights instead lowers the headline to " f"{rng(s_sc.f_min, s_sc.f_max, '{:.3f}')}.") A("") A("The two vertical models can be told apart with the data. If the " "weights are right, the slab stars and their orbits predict the age " "mix of the tracers found in layers above the slab; survey sampling " "varies with height but not, at fixed height, with age.") A("") A(table(["Layer |Z|", "LAMOST observed", "predicted: dynamical weights", "predicted: star counts (exp)", "predicted: star counts, thick disc 600/900 pc"], layer_rows(lt, ["observed", "dynamical weights", "star counts (exp)", "star counts (exp_thick)"]))) A("") A(table(["Layer |Z|", "Nataf+24 observed", "predicted: dynamical weights"], layer_rows(nl, ["observed", "dynamical weights"]))) A("") A(f"(Each cell: tracer fraction older than the Sun / older than 8 Gyr.) " f"For LAMOST the dynamical prediction is within " f"{dx[0]:+.3f} to {dx[1]:+.3f} of the observed fraction older than the " f"Sun in every layer up to 800 pc, while star-count heights " f"under-predict it by {-dsc[1]:.2f}–{-dsc[0]:.2f} (exponential) and " f"{-dth[1]:.2f}–{-dth[0]:.2f} (with canonical thick-disc heights). For " f"Nataf+24 the dynamical prediction is high by {dn[0]:.2f}–{dn[1]:.2f}. " f"The star-count variants are therefore kept in the tables but not " f"used for the range (which contains {sc_in} of them anyway), and the " f"test does not reach beyond 800 pc: counting the column only to " f"0.8, 1.3 or 2 kpc gives " f"{xr.loc['column to |z| < 0.8 kpc only (height verified by the layer test)', 'f_gt_0']:.3f}, " f"{xr.loc['column to |z| < 1.3 kpc only', 'f_gt_0']:.3f} and " f"{xr.loc['column to |z| < 2 kpc only', 'f_gt_0']:.3f} (Fig. 3c).") A("") A("### 4.6 Error budget of the primary estimate") A("") rows = [] for r in sysd.itertuples(): rows.append((r.source, f"{r.f_min:.3f}" if abs(r.f_max - r.f_min) < 5e-4 else rng(r.f_min, r.f_max, "{:.3f}"), f"{r.delta_min:+.3f}" if abs(r.f_max - r.f_min) < 5e-4 else f"{r.delta_min:+.3f} / {r.delta_max:+.3f}", f2(r.mean_hs_min) if abs(r.mean_hs_max - r.mean_hs_min) < 5e-3 else rng(r.mean_hs_min, r.mean_hs_max))) rows.append(("age catalogue swapped star by star (Nataf+24 for X&R)", f3(B.f_gt_0 + d_swap), f"{d_swap:+.3f}", "–")) for lab, t in (("other core sample: Nataf+24 primary (four age " "scales; > 3 Gyr)", est[(est.tier == "core") & est.estimate.str.startswith("Nataf")]), ("bracketing samples: Gaia 200 pc on published scales; " "Nataf+24 primary + annex", est[est.tier == "bracket"])): rows.append((lab, rng(t.f_older.min(), t.f_older.max(), "{:.3f}"), f"{t.f_older.min() - B.f_gt_0:+.3f} / " f"{t.f_older.max() - B.f_gt_0:+.3f}", rng(t.mean_headstart_older.min(), t.mean_headstart_older.max()))) A(table(["Source of uncertainty (LAMOST sample)", "f older", "shift from baseline", "mean head start of the older [Gyr]"], rows)) A("") A("The young-fraction prior is already inside the statistical interval. " "The × 0.87 and × 1.146 rows are the scale of this study's own " "Gaia-only ages, which no published catalogue supports; they are " "listed because the referee asked for them (Fig. 5b).") A("") A("### 4.7 Mid-plane samples understate the old fraction") A("") A(table(["Sample and method", "column minus mid-plane"], corr_rows)) A("") A(f"For subgiant and FLAME samples the correction is " f"+{corr_lo:.2f} to +{corr_hi:.2f}, depending on sample and method; " f"the 40 pc GCS dwarfs give more, with ages too uncertain to weigh " f"heavily (Gaia FLAME and GCS ages of the same dwarfs differ by " f"{fvg.loc['all Sun-like in common', 'median_diff']:+.1f} Gyr in the " f"median, Fig. 3d). This is why mid-plane literature values (0.50 to " f"{alz.f_gt_0:.2f}) cannot be compared directly with Lineweaver's " f"number.") A("") A("### 4.8 The Fermi clock") A("") rows = [] for v in rngc.index: r = rngc.loc[v] if r.n_eff_fiducial <= 0: rows.append((f"{v:g}", f"{F.D_GAL_LY / v / 1e9:.2f}", "0", "0", "no limit", "no limit", "–")) continue rows.append((f"{v:g}", f"{F.D_GAL_LY / v / 1e9:.3g}", sci(r.n_eff_fiducial, 2), f"{sci(r.n_eff_min)}–{sci(r.n_eff_max)}", sci(r.p95_fiducial, 2), f"{sci(r.p95_min)}–{sci(r.p95_max)}", f"{sci(r.p95_min_with_floor)}–{sci(r.p95_max_with_floor)}")) A(table(["v / c", "settling time [Gyr]", "N_eff, primary", f"N_eff over {int(rngc.n_distributions.iloc[0])} age distributions", "p95 = 3.0 / N_eff, primary", "p95 over the same distributions", "p95 with metallicity floors"], rows)) A("") A(f"These limits are conditional on the fiducial star count and " f"eta_Earth, and on the premise that no settlement has occurred and " f"would be evident if it had. The age-distribution uncertainty moves " f"N_eff by less than " f"{100 * max(1 - c_slow.n_eff_min / c_slow.n_eff_fiducial, c_slow.n_eff_max / c_slow.n_eff_fiducial - 1):.0f}% " f"at 1e-4 c; the assumed star count and eta_Earth move it by a factor " f"of several (limit {sci(p_ext_lo)} to {sci(p_ext_hi)} at their joint " f"extremes). Treating N_eff as log-normally uncertain " f"(sigma_ln = {marg.sigma_ln.min():.1f}–{marg.sigma_ln.max():.1f}) " f"with a flat prior on p loosens the 95% limit by a factor " f"{kmin:.1f}–{kmax:.1f}: **p < {sci(pm_lo)}–{sci(pm_hi)} at 0.1 c " f"and {sci(pm_lo_s)}–{sci(pm_hi_s)} at 1e-4 c**. The clock is nearly " f"independent of speed above about 3e-5 c " f"(half of the older stars have a head start exceeding the crossing " f"time at {sci(rk.v_half_over_c)} c) and stops below " f"{sci(rk.v_min_over_c)} c, where crossing takes longer than any head " f"start. Hard-step emergence models (n = 1, 3, 5; habitable lifetime 5.5 Gyr) " f"change N_eff at 1e-3 c by factors of " f"{hard.loc['hard steps n=1', 'n_eff'] / hard.loc['Earth clock', 'n_eff']:.2f}, " f"{hard.loc['hard steps n=3', 'n_eff'] / hard.loc['Earth clock', 'n_eff']:.2f} and " f"{hard.loc['hard steps n=5', 'n_eff'] / hard.loc['Earth clock', 'n_eff']:.2f} " f"(Fig. 4).") A("") A("### 4.9 The Sun's rank, and the Galaxy beyond the solar cylinder") A("") g_b = gal[gal.method == "bulge assumption"].iloc[1] g_a = gal[gal.method.str.startswith("Xiang")].iloc[0] A(f"If star formation continues at the recent rate for another 5 or 10 " f"Gyr, the Sun ends up with " f"{100 * rk.f_older_incl_future_const_5Gyr:.0f}% or " f"{100 * rk.f_older_incl_future_const_10Gyr:.0f}% of all Sun-like " f"stars older than itself: unremarkable. Two illustrative brackets " f"extend the local number to the Galaxy, both by assumption: a bulge " f"holding 15–35% of the stars and 80–100% old gives " f"{rng(g_b.f_galaxy_min, g_b.f_galaxy_max)} from the adopted local " f"value; taking Xiang+ 2025's present-day high-alpha disc mass " f"(2e10 Msun, formed about 11 Gyr ago) as 33–50% of the Galaxy's " f"stars and giving the rest the local low-alpha age mix " f"({g_a.f_local:.2f} older than the Sun) gives " f"{rng(g_a.f_galaxy_min, g_a.f_galaxy_max)}. Neither is a measurement.") A("") A("## 5. Robustness checks") A("") def v(tab, key): return f"{tab.loc[key, 'f_gt_0']:.3f}" A(f"- **LAMOST** (baseline {B.f_gt_0:.3f}): no error deconvolution " f"{v(xr, 'no error deconvolution')}; potentials " f"{rng(srow('vertical model: potential').f_min, srow('vertical model: potential').f_max, '{:.3f}')}; " f"shell, slab and annulus choices " f"{rng(srow('selection geometry').f_min, srow('selection geometry').f_max, '{:.3f}')}; " f"isochrone shifts and IMF slope " f"{rng(srow('tracer weights').f_min, srow('tracer weights').f_max, '{:.3f}')}; " f"age error below 15% only {v(xr, 'age error < 15%')}; tracers older " f"than 2.5 Gyr only {v(xr, 'tracer ages > 2.5 Gyr only')}; young " f"share 0.06 / 0.13 / 0.15: {v(xr, 'phi_young = 0.06')} / " f"{v(xr, K13)} / " f"{v(xr, 'phi_young = 0.15')}.") A(f"- **Nataf+24 primary** (baseline {NP.f_gt_0:.3f}): without the " f"catalogue quality flags {v(nat, 'without the catalogue quality flags')}; " f"RUWE < 1.2 {v(nat, 'RUWE < 1.2 only')}; |b| > 15° " f"{v(nat, '|b| > 15 deg only')}; A_V < 0.25 {v(nat, 'A_V < 0.25 only')}; " f"isochrone colour ± 0.02 mag " f"{v(nat, 'isochrone colours -0.02 mag')} / {v(nat, 'isochrone colours +0.02 mag')}; " f"tracers older than 3 or 4 Gyr (young share from published histories) " f"{v(nat, 'tracer ages > 3 Gyr (young prior for 3 Gyr)')} / " f"{v(nat, 'tracer ages > 4 Gyr (young prior for 4 Gyr)')}. The age mix " f"is flat with distance from 0.15 to 1.3 kpc (Fig. 3a).") A(f"- **Nataf+24 annex.** Adding the metal-poor annex raises the " f"estimate to {NA.f_gt_0:.3f} (annex weight share " f"{100 * NA.annex_weight_share:.0f}%); moving the assumed ultraviolet " f"completeness ramp by ± 0.1 dex changes it by at most {uv:.3f}. But " f"the combined sample is kinematically colder at old ages than the " f"primary sample (rms v_z {ka.rms_vz_kms.iloc[-2]:.0f} against " f"{kp.rms_vz_kms.iloc[-2]:.0f} km/s at 10–12 Gyr), which points to " f"contamination of the annex by younger stars with overestimated " f"ages. It defines the upper edge of the range and is not preferred.") A(f"- **Gaia 200 pc** (as computed {G0.f_gt_0:.3f}; on the X&R scale " f"{GX.f_gt_0:.3f}): RUWE < 1.4 {v(gs, 'RUWE < 1.4 only')}; d < 100 pc " f"{v(gs, 'd < 100 pc')}; p_sg > 0.5 {v(gs, 'p_sg > 0.5 only')}; " f"binary-safe window {v(gs, 'binary-safe window (above the twin sequence)')}; " f"[M/H] ± 0.1 dex {v(gs, '[M/H] shifted by -0.1 dex')} / " f"{v(gs, '[M/H] shifted by +0.1 dex')} (all on the ages as computed).") A(f"- **Unevolved dwarfs** (cross-check only): GCS 40 pc column " f"{ms.loc['GCS d<40 pc, 0.9-1.2 Msun, Padova, column', 'deconv_f_gt_0']:.2f} " f"(Padova ages) and " f"{ms.loc['GCS d<40 pc, 0.9-1.2 Msun, BaSTI, column', 'deconv_f_gt_0']:.2f} " f"(BaSTI); Gaia FLAME 200 pc column " f"{ms.loc['Gaia FLAME d<200 pc, 0.8-1.2 Msun, column', 'deconv_f_gt_0']:.2f}; " f"solar twins, raw, " f"{twins.loc['Spina+ 2018 solar twins (HARPS)', 'f_older']:.2f} (Spina+18), " f"{twins.loc['Nissen+ 2020 solar-type (HARPS)', 'f_older']:.2f} (Nissen+20), " f"{twins.loc['Taniguchi+ 2026 Gaia GSP-Spec twins, Age(GMAG)', 'f_older']:.2f} " f"(Taniguchi+26).") A("") A("## 6. Limitations") A("") A(f"- **Age scale.** The Sun's age falls where the cumulative age " f"distribution is steep. Published catalogues differ by up to " f"{100 * max(abs(1 - r_fl), abs(1 - r_sh), abs(1 - r_nat)):.0f}% for " f"the same stars, worth {s_age_pub.delta_min:+.3f} to " f"{s_age_pub.delta_max:+.3f}, and {d_swap:+.3f} in the star-by-star " f"swap. No catalogue is known to be right.") A(f"- **Vertical structure.** The weights are one-dimensional. The " f"layer test supports them to 800 pc for LAMOST and shows a " f"{dn[0]:.2f}–{dn[1]:.2f} over-prediction for Nataf+24; above 800 pc " f"they are untested, and the oldest bin's implied height " f"({hz.h_eff_pc.iloc[7]:.0f} pc) rests on {int(hz.n.iloc[7])} stars.") A(f"- **Three samples, two answers.** LAMOST and the Nataf+24 primary " f"sample agree to {abs(B.f_gt_0 - NX.f_gt_0):.2f} on a common age " f"scale; the Gaia 200 pc sample is lower by {rng(min(dg), max(dg))} " f"and the cause is not established.") A("- **Selection functions are not modelled star by star.** The method " "relies on age-independent sampling inside the shell, checked only " "through the distance and layer tests. The Nataf+24 primary box " "under-samples metal-poor stars and its annex cannot be modelled " "exactly; the two bracket the truth.") A("- **Young stars are imported.** The 0.10 ± 0.04 share of stars " "younger than 2 Gyr comes from published star-formation histories; " f"the tracers themselves suggest " f"{B.phi_young_if_flat:.2f} if the recent rate was flat.") A("- **Solar cylinder, not Galaxy; residence, not birthplace.** Radial " "migration means these stars were not all born at the Sun's radius. " "The inner Galaxy is older; the Galaxy-wide brackets are assumptions.") A("- **Stars, not planets.** Metallicity floors are applied to stars, " "every star younger than 2 Gyr is assumed to pass them, and eta_Earth " "is assumed independent of age and metallicity.") A("- **The Fermi limit is conditional** on an assumed star count, " "eta_Earth, an emergence-time model no slower than Earth's, the solar " "age mix holding Galaxy-wide, and the premise that settlement has not " "happened and would be evident. If settlement would have prevented " "our existence, the limit becomes an anthropic statement.") A("- **Not used:** GALAH and APOGEE ages; Kordopatis+ 2023 ages (service " "time-outs); the Gallart+ 2024 star-formation history beyond its one " "quoted number (not tabulated).") A("") A("## 7. Closest prior work and what is new") A("") A("The answer is not new. Lineweaver (2001) and Lineweaver et al. (2004) " "obtained about three quarters from models of star formation and " "metal build-up (cosmic and Galactic, respectively); " "Legassick (2015, arXiv:1509.02832) obtained 77% of potentially " "habitable systems older than the Sun, by 3.13 Gyr on average, from a " "Galactic model, which is the same pair of statistics; Behroozi & " "Peeples (2015) obtained 80% of present Earth-like planets already " "formed by the Sun's birth; Zackrisson et al. (2016) obtained a mean " "age of 7 ± 1 Gyr for terrestrial planets around FGK stars; Robles et " "al. (2008), ranking the Sun against a chromospheric star-formation " "history, found 53 ± 2% of stars older than the Sun. Madau " "(2023) is the closest empirical work: an analytic fit to a Gaia " "mid-plane star-formation history within 100 pc combined with planet " "occurrence, giving 77%. The data products used here are also " "published: Xiang & Rix (2022) and Nataf et al. (2024) for subgiant " "ages, Mazzi et al. (2024) and Gallart et al. (2024) for local " "star-formation histories, Xiang et al. (2025) for a selection-" f"function-corrected model of the same LAMOST subgiants. Integrating " f"Mazzi+24's published column history takes a few lines and gives " f"{mazzi.f_gt_0:.2f}.") A("") A("This study is best read as a reproduction of the known \"about three " "quarters\" with newer data and an honest systematic range. What it " "adds: (i) the statistic from individual subgiant ages in three " "samples rather than a model or a fitted history; (ii) an explicit, " "orbit-based mid-plane-to-column correction, with a layer test that " "discriminates it from star-count scale heights; (iii) an age-scale " "budget from stars shared between catalogues, including a star-by-star " "swap; (iv) the head-start distribution and metallicity-floor " "variants on [Fe/H] and [M/H]; (v) the Fermi clock computed from a " "measured distribution, with the limit stated both conditionally and " "marginalised. I found no publication that does (ii) or (iii) for " "this statistic, but I cannot exclude one.") A("") A("## 8. Implications for the Fermi Paradox") A("") A(f"The head start is real and large: " f"{100 * ad.f_range_lo:.0f}–{100 * ad.f_range_hi:.0f}% of Sun-like " f"stars near the Sun's orbit are older than the Sun, by a median of " f"{rng(ad.median_hs_lo, ad.median_hs_hi, '{:.1f}')} Gyr, and " f"{100 * ad.f5_range_lo:.0f}–{100 * ad.f5_range_hi:.0f}% of all " f"Sun-like stars by more than 5 Gyr. At 1e-3 c the Galaxy is crossed " f"in 0.1 Gyr, some {ad.median_hs_lo / 0.1:.0f} times less than that " f"median. Under an Earth-clock emergence model, **travel time alone " f"cannot explain the silence for front speeds above about 3e-5 c**: " f"the number of planets with enough head start is already " f"{100 * rngc.loc[3e-5, 'n_eff_fiducial'] / c_fast.n_eff_fiducial:.0f}% " f"of its maximum at 3e-5 c and " f"{100 * c_slow.n_eff_fiducial / c_fast.n_eff_fiducial:.0f}% at " f"1e-4 c. The 95% limit on the per-planet probability of " f"producing a Galaxy-settling civilisation, about 1e-9 for the " f"fiducial planet count and {sci(pm_lo, 0)} to {sci(pm_hi, 0)} with " f"the count's uncertainty marginalised, holds only if no settlement " f"has occurred and a settled Galaxy would be evident. The census " f"fixes the size of the clock; it says nothing about which " f"explanation is right (rarity, short lifetimes, no expansion, or " f"undetected presence), and nothing about emergence times longer than " f"Earth's. If intelligence typically needs much longer than 4.6 Gyr, " f"or more metals than the old disc had, the clock shrinks: with the " f"hardest floor considered ([Fe/H] > −0.3) the limit is " f"{sci(rngc.loc[1e-4, 'p95_max_with_floor'])} at 1e-4 c, still of " f"order 1e-9.") A("") A("## Figures") A("") A("- `figures/fig1_age_distribution.png`: formation history of Sun-like " "stars now in the solar cylinder and the cumulative head-start " "distribution, three samples and published histories.") A("- `figures/fig2_fraction_older.png`: every estimate of the fraction " "older than the Sun, with the adopted range and Lineweaver's interval.") A("- `figures/fig3_validation.png`: distance-shell test, scale heights, " "layer test, disagreement of dwarf ages.") A("- `figures/fig4_fermi_clock.png`: N_eff and the limit on p against " "expansion speed, conditional and marginalised.") A("- `figures/fig5_systematics.png`: catalogue age scales for stars in " "common, and what moves the headline number.") A("") A("## Response to referee") A("") resp = [ ("1. Headline wording and error bar", f"Accepted. The headline is now \"about {100 * ad.f_adopted:.0f}%\", " f"{ad.f_adopted:.2f} with a systematic range " f"{rng(ad.f_range_lo, ad.f_range_hi)} (core samples " f"{rng(ad.f_core_lo, ad.f_core_hi)}), with the primary sample's " f"{B.f_gt_0:.3f} ({rng(B.f_gt_0_lo, B.f_gt_0_hi, '{:.3f}')}) labelled " f"statistical plus young prior, and \"consistent with\" Lineweaver " f"instead of \"confirms\"."), ("2. Gaia sample is not an independent replication", f"Accepted. It is now a cross-check. Offsets reported: " f"{d_gx:+.2f} Gyr and ratio {r_gx:.2f} against Xiang & Rix, " f"{d_gs:+.2f} Gyr and ratio {r_gs:.2f} against StarHorse for the " f"{n_(n_gs)} stars with p_sg > 0.5 " f"({shc.loc['all', 'median_diff_gyr']:+.2f} Gyr in the first-pass " f"validation table), ratio {r_gn:.2f} against Nataf+24. On " f"the X&R scale it gives " f"{GX.f_gt_0:.3f} (column) and {GXm.f_gt_0:.3f} (mid-plane); the " f"{rng(min(dg), max(dg))} disagreement with LAMOST on a common " f"scale is stated. " f"The analysis was also redone with Nataf+24 ages as a new second " f"sample ({NP.f_gt_0:.3f}; {NX.f_gt_0:.3f} on the X&R scale)."), ("3. Systematic error budget", "Accepted: section 4.6, Fig. 5b, " "`data/xr_systematics.csv`, `data/headline_estimates.csv`."), ("4. Fig. 3b and scale heights", f"Accepted. Panel retitled, sech² heights added, the excess above " f"4 Gyr and the {100 * B.high_alpha_share_tracers:.0f}% high-alpha " f"column share stated (section 4.5)."), ("5. Mean head start with a systematic range", f"Accepted: {B.mean_headstart_older:.1f} Gyr, " f"{rng(ad.mean_hs_lo, ad.mean_hs_hi, '{:.1f}')} Gyr over core and " f"bracketing estimates, wider with the disfavoured variants " f"(section 4.2). " f"The × 1.146 values are now {s_age_ref.f_max:.3f} and " f"{s_age_ref.mean_hs_max:.1f} Gyr rather than 0.759 and 4.6 Gyr: " f"the age-rescaling option of the first-pass code dropped stars " f"pushed beyond 14 Gyr, which is fixed."), ("6. Mid-plane correction", f"Accepted: quoted as " f"+{corr_lo:.2f} to +{corr_hi:.2f} by sample and method " f"(section 4.7)."), ("7. Gallart+24", "Accepted: mid-plane 0.50 added to the tables and " "Fig. 2; the mid-plane literature range reads 0.50 to 0.76."), ("8. Metallicity floors", "Accepted: floors labelled [Fe/H], [M/H] " "variants added, Buchhave+ 2012/2014 and Kepler-444 cited, floors " "described as lower bounds (section 4.3)."), ("9. Fermi limit is conditional", f"Accepted: conditional and marginalised limits are given together " f"with the premise. My marginalised range ({sci(pm_lo)}–{sci(pm_hi)}) " f"reaches higher than the referee's 3e-9 to 6e-9 because reading " f"\"a factor 2 either way\" in the star count as 1 sigma gives " f"sigma_ln = {marg.sigma_ln.max():.2f}."), ("10. Implications", "Accepted: softened to the referee's wording " "(section 8)."), ("11. Prior work", "Accepted: Legassick 2015, Nataf+ 2024 and " "Xiang+ 2025 cited; the result is described as a reproduction; " "Xiang+25's high-alpha disc mass gives a second Galaxy-wide " "bracket beside the bulge assumption."), ("12. Document the cuts", "Accepted: cut-flow tables with reasons " "and counts for all three samples; zero-weight stars, bright-limit " "loss, radial-velocity loss and RUWE fraction quantified " "(section 2)."), ("13. Wording and report file", "Accepted: \"now residing in the " "solar cylinder\" throughout; this report replaces the " "auto-rendered file and is regenerated by the pipeline."), ("Not accepted in full: star-count scale heights as a −0.03 to " "−0.06 error term (referee 2.3)", f"The referee wrote that neither set of heights is clearly right. " f"The layer test (section 4.5) decides it for LAMOST: orbit weights " f"predict the observed age mix at 200–800 pc to within " f"{max(abs(dx[0]), abs(dx[1])):.2f}, star-count heights miss it by " f"{-dsc[1]:.2f}–{-dsc[0]:.2f}. I therefore do not treat the " f"star-count values as equally likely. What remains is smaller: " f"Nataf+24's layers are over-predicted by {dn[0]:.2f}–{dn[1]:.2f}, " f"and the column above 800 pc is untested " f"({s_col.delta_min:+.3f} if it is dropped altogether). The full " f"range contains {sc_in} of the star-count values regardless."), ("Not accepted in full: ages × 0.87 and × 1.146 as plausible age " "scales (referee 2.1, 2.2)", f"Both factors come from this study's Gaia-only ages, which are " f"{100 * (g_lo - 1):.0f}–{100 * (g_hi - 1):.0f}% older than three " f"published catalogues for the same stars, including Nataf+24, " f"which uses the same MIST models. Published scales span × " f"{r_fl:.2f} to × {r_sh:.2f}. The referee's size for the effect is " f"nevertheless confirmed by another route: the star-by-star swap " f"to Nataf+24 ages gives {d_swap:+.3f}."), ] A(table(["Referee point", "Response"], resp, align=["l", "l"])) A("") A(f"The referee's reproduction, arithmetic checks and transcription " f"checks needed no action. The referee proposed 0.64 to 0.76; the " f"range found here is {rng(ad.f_range_lo, ad.f_range_hi)}: the lower " f"edge is set, as in the referee's report, by the rescaled Gaia " f"sample, and the upper edge is higher because of the Nataf+24 annex, " f"which the referee did not have.") A("") A("## Reproducibility") A("") A("`python3 code/run_all.py` regenerates every table, figure and this " "report (about 15 minutes on one core once the raw data are " "downloaded; `--quick` uses fewer bootstrap resamples; `--no-fetch` " "skips the network). Raw downloads go to " "`fermi_raw/headstart/` and are re-fetched if missing: " "VizieR and Gaia TAP queries, MIST isochrones, and the Nataf+24 " "catalogue from the link in the authors' repository " "(github.com/DavidMoiseNataf/Subgiants). Random seeds are fixed. " "Python 3 with numpy, scipy, pandas, matplotlib and astropy.") A("") text = "\n".join(L) # typographic minus for signed numbers (not inside 1e-9 or 0.3-0.8) import re text = re.sub(r"(?