#!/usr/bin/env python3 """Figures for the head-start study (reads only the CSVs in data/). fig1 age distribution and cumulative head-start distribution fig2 forest plot: every estimate of the fraction older than the Sun fig3 validation: selection shells, scale heights, layer test, dwarf ages fig4 the Fermi clock fig5 age scales and the systematic budget of the headline number """ import os import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import matplotlib.ticker import common as C import published_sfh as P import fermi_clock as F FIG = C.FIGS # Okabe-Ito colour-blind-safe palette (validated: worst adjacent CVD dE 11) BLUE, ORANGE, GREEN, RED, PURPLE, GREY, SKY = ( "#0072B2", "#E69F00", "#009E73", "#D55E00", "#CC79A7", "#7f7f7f", "#56B4E9") INK, INK2 = "#0b0b0b", "#52514e" plt.rcParams.update({"font.size": 10, "axes.spines.top": False, "axes.spines.right": False, "figure.dpi": 150, "savefig.bbox": "tight", "legend.frameon": False, "text.color": INK, "axes.labelcolor": INK, "axes.edgecolor": INK2, "xtick.color": INK2, "ytick.color": INK2, "savefig.facecolor": "white"}) FID = F.FID NAT = F.NAT + "formed: all [Fe/H]" NATU = F.NAT + "primary + annex" GXR = F.GSG + "ages on the X&R22 scale: formed, all [M/H]" def _csv(name): path = C.DATA + "/" + name return pd.read_csv(path) if os.path.exists(path) else None def step(ax, phi, rebin=1, **kw): """Plot a binned distribution as stars per Gyr (fraction / bin width). rebin: number of 0.25-Gyr bins merged for display (the deconvolved tracer histograms are noisy at full resolution).""" y = phi.reshape(-1, rebin).sum(axis=1) / phi.sum() e = C.AGE_EDGES[::rebin] y = y / np.diff(e) ax.step(e, np.r_[y, y[-1]], where="post", **kw) def survival(phi): """x (head start, Gyr), fraction of stars with age > T_SUN + x.""" x = np.linspace(-C.T_SUN, 9.4, 400) return x, np.array([C.frac_older(phi, xx) for xx in x]) def fig1(dists, adopted): fig, ax = plt.subplots(1, 2, figsize=(11.8, 4.5)) pub = {k: v[1] for k, v in dists.items() if k.startswith("published") and v[0] == "column"} a = ax[0] Y = np.array([p / p.sum() / 0.25 for p in pub.values()]) a.fill_between(C.AGE_EDGES, np.r_[Y.min(0), Y.min(0)[-1]], np.r_[Y.max(0), Y.max(0)[-1]], step="post", color="0.87", label="published CMD-fit column SFHs (7 solutions)") step(a, dists["published SFH | Mazzi+24"][1], color=GREY, lw=1.2, label="Mazzi+24 (Gaia CMD fit)") if GXR in dists: step(a, dists[GXR][1], rebin=4, color=ORANGE, lw=1.2, label="Gaia 200 pc subgiants, own ages put on the X&R scale") if NAT in dists: step(a, dists[NAT][1], rebin=4, color=GREEN, lw=1.6, label="Nataf+24 subgiants, primary sample") step(a, dists[FID][1], rebin=4, color=BLUE, lw=2.0, label="LAMOST subgiants (Xiang & Rix ages), column-weighted") a.axvline(C.T_SUN, color=INK, ls="--", lw=1) a.text(C.T_SUN - 0.12, 0.012, "Sun", rotation=90, va="bottom", ha="right") a.plot([2.0, 2.0], [0, 0.205], color=INK2, ls=":", lw=0.8) a.text(0.1, 0.255, "< 2 Gyr: not traced\nby subgiants; share\n" "(0.10 $\\pm$ 0.04) from\npublished SFHs", ha="left", va="top", fontsize=6.5, color=INK2) a.set_xlim(0, 14); a.set_ylim(0, 0.26) a.set_xlabel("stellar age [Gyr]") a.set_ylabel("fraction of Sun-like stars formed per Gyr") a.set_title("(a) Formation history of Sun-like stars now in the solar " "cylinder", fontsize=10, loc="left") a.legend(fontsize=7.5, loc="upper center", ncol=2, bbox_to_anchor=(0.5, -0.16)) b = ax[1] S = np.array([survival(p)[1] for p in pub.values()]) x = survival(dists[FID][1])[0] b.fill_between(x, S.min(0), S.max(0), color="0.87", label="published CMD-fit column SFHs") curves = [(FID, BLUE, "-", 2.0, "LAMOST, all metallicities"), (F.XR + "formed: [M/H]>-0.5", BLUE, "--", 1.3, "LAMOST, [M/H] > $-0.5$"), (F.XR + "formed: [Fe/H]>-0.5", BLUE, ":", 1.3, "LAMOST, [Fe/H] > $-0.5$"), (NAT, GREEN, "-", 1.6, "Nataf+24 primary sample"), (NATU, GREEN, "--", 1.2, "Nataf+24 primary + annex"), (GXR, ORANGE, "-", 1.2, "Gaia 200 pc, X&R age scale")] for key, col, ls, lw, lab in curves: if key not in dists: continue xx, s = survival(dists[key][1]) b.plot(xx, s, color=col, ls=ls, lw=lw, label=lab) b.errorbar([0], [0.74], yerr=[0.09], fmt="s", color=RED, ms=5, capsize=3, label="Lineweaver (2001) model: 74 $\\pm$ 9 %", zorder=5) b.errorbar([0.35], [adopted["f_adopted"]], yerr=[[adopted["f_adopted"] - adopted["f_range_lo"]], [adopted["f_range_hi"] - adopted["f_adopted"]]], fmt="o", color=INK, ms=5, capsize=3, zorder=6, label="this work: %.2f (range %.2f-%.2f)" % (adopted["f_adopted"], adopted["f_range_lo"], adopted["f_range_hi"])) b.axvline(0, color=INK, ls="--", lw=1) b.set_xlim(-4.6, 9.4); b.set_ylim(0, 1) b.set_xlabel("head start over the Sun [Gyr] (age $-$ 4.57 Gyr)") b.set_ylabel("fraction of Sun-like stars with a larger head start") b.set_title("(b) Cumulative head-start distribution", fontsize=10, loc="left") b.legend(fontsize=7.2, loc="upper right") fig.tight_layout() fig.savefig(FIG + "/fig1_age_distribution.png") plt.close(fig) def fig2(summ, est, adopted): """Forest plot: fraction of Sun-like stars older than the Sun.""" s = summ.set_index("distribution") ms = pd.read_csv(C.DATA + "/ms_census.csv").set_index("sample") tw = pd.read_csv(C.DATA + "/solar_twins.csv").set_index("sample") pub = pd.read_csv(C.DATA + "/published_sfh_table.csv").set_index("name") rows = [] # (label, value, lo, hi, colour, marker, filled) def tr(key, lab, col): if key not in s.index: return r = s.loc[key] rows.append((lab, r.f_gt_0, r.get("f_gt_0_lo", np.nan), r.get("f_gt_0_hi", np.nan), col, "o", True)) def head(t): rows.append((t, None, 0, 0, INK, "", True)) head("MODELS IN THE LITERATURE") for m in P.MODELS: e = m["err"] rows.append((m["name"], m["f"], m["f"] - e if e else np.nan, m["f"] + e if e else np.nan, RED, "s", True)) head("THIS WORK: SUBGIANT TRACERS, STARS FORMED, COLUMN (filled: " "core; open: bracketing sample; x: not used)") colmap = {"LAMOST": BLUE, "Nataf": GREEN, "Gaia": ORANGE} for r in est.itertuples(): col = colmap[r.estimate.split(",")[0].split("+")[0].split(" ")[0]] rows.append((r.estimate.replace(" (paired shift applied to the primary)", " (paired)"), r.f_older, r.stat_lo, r.stat_hi, col, "x" if r.tier == "excluded" else "o", r.tier == "core")) head("WITH A METALLICITY FLOOR (LOWER BOUNDS ON THE HEAD START)") tr(F.XR + "formed: [M/H]>-0.5", "LAMOST, [M/H] > -0.5", BLUE) tr(F.XR + "formed: [Fe/H]>-0.5", "LAMOST, [Fe/H] > -0.5", BLUE) tr(F.XR + "formed: [M/H]>-0.3", "LAMOST, [M/H] > -0.3", BLUE) tr(F.XR + "formed: [Fe/H]>-0.3", "LAMOST, [Fe/H] > -0.3", BLUE) tr(F.NAT + "formed: [Fe/H]>-0.5", "Nataf+24 primary, [Fe/H] > -0.5", GREEN) tr(F.NAT + "formed: [Fe/H]>-0.3", "Nataf+24 primary, [Fe/H] > -0.3", GREEN) head("PUBLISHED CMD-FIT SFHs, RE-INTEGRATED HERE (x = mid-plane)") for k, r in pub.iterrows(): mid = str(r.kind).startswith("midplane") rows.append((k + (" [mid-plane]" if mid else ""), r.f_gt_0, np.nan, np.nan, GREY, "x" if mid else "D", True)) head("LIVING MAIN-SEQUENCE STARS AND MID-PLANE SAMPLES") tr(F.XR + "living MS: all [Fe/H]", "LAMOST tracers: still on the main sequence, column", PURPLE) tr(F.XR + "no column weights (mid-plane slab)", "LAMOST tracers, formed, mid-plane (no column weights)", SKY) tr(F.NAT + "no column weights (mid-plane slab)", "Nataf+24 primary, formed, mid-plane", SKY) for key, lab, col in [ ("GCS d<40 pc, 0.9-1.2 Msun, Padova, column", "GCS 40 pc dwarfs, column (Padova ages)", PURPLE), ("GCS d<40 pc, 0.9-1.2 Msun, BaSTI, column", "GCS 40 pc dwarfs, column (BaSTI ages)", PURPLE), ("GCS d<40 pc, 0.9-1.2 Msun, Padova, volume", "GCS 40 pc dwarfs, mid-plane volume (Padova)", SKY), ("GCS d<40 pc, 0.9-1.2 Msun, BaSTI, volume", "GCS 40 pc dwarfs, mid-plane volume (BaSTI)", SKY), ("Gaia FLAME d<200 pc, 0.8-1.2 Msun, column", "Gaia FLAME 200 pc dwarfs, column", PURPLE), ("Gaia FLAME d<200 pc, 0.8-1.2 Msun, volume", "Gaia FLAME 200 pc dwarfs, mid-plane volume", SKY)]: r = ms.loc[key] rows.append((lab, r.deconv_f_gt_0, r.deconv_f_gt_0_lo, r.deconv_f_gt_0_hi, col, "o", True)) for key, lab in [("Spina+ 2018 solar twins (HARPS)", "Spina+18 solar twins (N=79), raw"), ("Nissen+ 2020 solar-type (HARPS)", "Nissen+20 solar-type (N=72), raw"), ("Taniguchi+ 2026 Gaia GSP-Spec twins, Age(GMAG)", "Taniguchi+26 Gaia twins (N=6594), raw")]: r = tw.loc[key] rows.append((lab, r.f_older, r.f_older - r.f_older_binom_err, r.f_older + r.f_older_binom_err, GREY, "^", True)) fig, ax = plt.subplots(figsize=(9.6, 0.245 * len(rows) + 1.0)) y = np.arange(len(rows))[::-1] ax.axvspan(adopted["f_range_lo"], adopted["f_range_hi"], color=BLUE, alpha=0.08, lw=0) ax.axvspan(adopted["f_core_lo"], adopted["f_core_hi"], color=BLUE, alpha=0.12, lw=0) ax.axvline(adopted["f_adopted"], color=BLUE, lw=0.9, alpha=0.7) ax.axvspan(0.65, 0.83, ymin=0, ymax=1, facecolor="none", edgecolor=RED, lw=0.8, ls="--") for yy, (lab, v, lo, hi, col, mk, filled) in zip(y, rows): if v is None: ax.text(0.185, yy, lab, fontsize=7.6, fontweight="bold", va="center") continue if np.isfinite(lo) and np.isfinite(hi): ax.plot([lo, hi], [yy, yy], color=col, lw=1.4) ax.plot([v], [yy], mk, color=col, ms=5, mfc=col if filled else "white", mew=1.2) ax.text(0.985, yy, f"{v:.2f}", va="center", ha="right", fontsize=7.6) ax.set_yticks([yy for yy, r in zip(y, rows) if r[1] is not None]) ax.set_yticklabels([r[0] for r in rows if r[1] is not None], fontsize=7.4) ax.set_xlim(0.18, 1.0) ax.set_ylim(-0.7, len(rows) - 0.3) ax.set_xlabel("fraction of Sun-like stars older than the Sun") ax.set_title("Who had a head start? Blue bands: core range %.2f-%.2f " "(dark) and full range %.2f-%.2f (light); line: %.2f.\n" "Dashed box: Lineweaver (2001), 74 $\\pm$ 9 %%" % (adopted["f_core_lo"], adopted["f_core_hi"], adopted["f_range_lo"], adopted["f_range_hi"], adopted["f_adopted"]), fontsize=9, loc="left") ax.tick_params(axis="y", length=0) fig.savefig(FIG + "/fig2_fraction_older.png") plt.close(fig) def fig3(): fig, axs = plt.subplots(2, 2, figsize=(11.8, 8.2)) a, b, c, d = axs.ravel() sh = pd.read_csv(C.DATA + "/xr_shell_test.csv") ns = _csv("nataf_shell_test.csv") a.axvspan(0.3, 0.8, color="0.90", lw=0, label="shell used (0.3-0.8 kpc)") mid = 0.5 * (sh.d_lo + sh.d_hi) a.errorbar(mid, sh.f_gt_sun, xerr=0.5 * (sh.d_hi - sh.d_lo), fmt="o", color=BLUE, ms=5, label="LAMOST: older than the Sun") a.errorbar(mid, sh.f_gt_8, xerr=0.5 * (sh.d_hi - sh.d_lo), fmt="s", color=BLUE, ms=5, mfc="white", label="LAMOST: older than 8 Gyr") if ns is not None: mid = 0.5 * (ns.d_lo + ns.d_hi) a.errorbar(mid, ns.f_gt_sun, xerr=0.5 * (ns.d_hi - ns.d_lo), fmt="o", color=GREEN, ms=5, label="Nataf+24: older than the Sun") a.errorbar(mid, ns.f_gt_8, xerr=0.5 * (ns.d_hi - ns.d_lo), fmt="s", color=GREEN, ms=5, mfc="white", label="Nataf+24: older than 8 Gyr") a.set_xlabel("distance from the Sun [kpc]") a.set_ylabel("tracer-weighted fraction ($|Z|<0.3$ kpc)") a.set_ylim(0, 1.22); a.set_xlim(0.1, 2.05) a.set_yticks([0, 0.2, 0.4, 0.6, 0.8]) a.set_title("(a) Selection: age mix against distance", fontsize=10, loc="left") a.legend(fontsize=7.5, loc="upper center", ncol=2) hz = pd.read_csv(C.DATA + "/xr_scale_heights.csv") xa = 10 ** (0.5 * (hz.logage_lo + hz.logage_hi) - 9) b.errorbar(xa * 0.97, hz.h_mazzi24_exp_pc, yerr=hz.h_mazzi24_err_pc, fmt="s", color=GREY, capsize=2, ms=5, label="Mazzi+24 star counts, exponential fit") b.plot(xa, hz.h_mazzi24_sech2_pc, "D", color=GREY, ms=4.5, mfc="white", label="Mazzi+24 star counts, sech$^2$ fit") b.errorbar(xa * 1.03, hz.h_eff_pc, yerr=hz.h_eff_err_pc, fmt="o", color=BLUE, capsize=2, ms=5, label="this work: vertical orbits of LAMOST subgiants") b.set_xlabel("age [Gyr]"); b.set_ylabel("effective scale height [pc]") b.set_yscale("log"); b.set_yticks([200, 300, 500, 1000]) b.set_yticklabels(["200", "300", "500", "1000"]) b.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter()) b.set_title("(b) Dynamical column weights against star-count scale " "heights:\n agree below 4 Gyr, dynamical heights larger " "above", fontsize=10, loc="left") b.legend(fontsize=7.5, loc="upper left") lt = pd.read_csv(C.DATA + "/xr_layer_test.csv") nl = _csv("nataf_layer_test.csv") lay = lt[(lt.model == "observed") & (lt.layer_lo > 0)] xs = np.arange(len(lay)) styles = [("observed", INK, "o", True, "LAMOST observed in the layer"), ("dynamical weights", BLUE, "s", True, "predicted from the slab: dynamical weights"), ("star counts (exp)", GREY, "D", False, "predicted: star-count heights (exp)"), ("star counts (exp_thick)", GREY, "^", False, "predicted: star-count heights, thick disc 600/900 pc")] for k, (model, col, mk, filled, lab) in enumerate(styles): t = lt[(lt.model == model) & (lt.layer_lo > 0)] c.errorbar(xs - 0.27 + 0.18 * k, t.f_older, yerr=t.f_older_err, fmt=mk, color=col, ms=6, capsize=2, mfc=col if filled else "white", label=lab) if nl is not None: for model, mk, filled, lab in [ ("observed", "o", True, "Nataf+24 observed in the layer"), ("dynamical weights", "s", False, "Nataf+24 predicted: dynamical weights")]: t = nl[(nl.model == model) & (nl.layer_lo > 0)] c.errorbar(np.arange(len(t)) + (0.36 if filled else 0.45), t.f_older, yerr=t.f_older_err, fmt=mk, color=GREEN, ms=5, capsize=2, mfc=GREEN if filled else "white", label=lab) c.set_xticks(xs + 0.05) c.set_xticklabels([f"{int(r.layer_lo)}-{int(r.layer_hi)} pc" for r in lay.itertuples()]) c.set_xlim(-0.6, len(lay) - 0.3) c.set_ylim(0.44, 1.0) c.set_xlabel("height above the plane, $|Z|$ (Nataf+24 top layer: " "500-750 pc)") c.set_ylabel("tracer fraction older than the Sun in the layer") c.set_title("(c) Layer test: which vertical model reproduces the age mix\n" " seen above the slab?", fontsize=10, loc="left") c.legend(fontsize=7.2, loc="lower right") fv = pd.read_csv(C.DATA + "/flame_vs_gcs.csv").iloc[1:] xs = np.arange(len(fv)) d.bar(xs, fv.median_diff, color=GREY, width=0.6) d.axhline(0, color=INK, lw=0.8) d.set_xticks(xs) d.set_xticklabels(["< -0.5", "-0.5 to -0.3", "-0.3 to -0.1", "-0.1 to 0.1", "> 0.1"], fontsize=8) for xx, (_, r) in zip(xs, fv.iterrows()): d.text(xx, r.median_diff + (0.15 if r.median_diff > 0 else -0.45), f"{r.median_diff:+.1f} Gyr\nN={int(r.n)}", ha="center", fontsize=7.5, va="bottom" if r.median_diff > 0 else "top") d.set_xlabel("[Fe/H] (Geneva-Copenhagen Survey)") d.set_ylabel("median age difference, FLAME $-$ GCS [Gyr]") d.set_ylim(-4.2, 5.8) d.set_title("(d) Catalogue ages of unevolved dwarfs disagree star by star", fontsize=10, loc="left") fig.tight_layout() fig.savefig(FIG + "/fig3_validation.png") plt.close(fig) def fig4(dists): frac = F.kroupa_number_fraction() n_sun = [n * frac for n in F.N_STARS] v = np.logspace(-5.1, -1, 90) ts = F.D_GAL_LY / v / 1e9 fig, ax = plt.subplots(1, 2, figsize=(11.8, 4.5)) a, b = ax eta, eta_lo, eta_hi = 0.37, 0.16, 1.50 # Bryson+21 full quoted span members = _csv("fermi_clock_range_members.csv") # band: all-metallicity column distributions on published age scales if members is not None: N = [] for k in members.distribution: _, phi, fp = dists[k] N.append(n_sun[0] * eta * fp * np.array([F.clock_factor(phi, t, "Earth clock") for t in ts])) N = np.array(N) lo, hi = N.min(0), N.max(0) ok = hi > 0 a.fill_between(v[ok], np.maximum(lo[ok], 1e6), hi[ok], color=BLUE, alpha=0.15, lw=0, label="range over %d age distributions (samples, age " "scales, published SFHs)" % len(N)) curves = [(FID, "Earth clock", BLUE, "-", "LAMOST, all metallicities"), (F.XR + "formed: [M/H]>-0.5", "Earth clock", BLUE, "--", "LAMOST, [M/H] > $-0.5$"), (F.XR + "formed: [Fe/H]>-0.3", "Earth clock", BLUE, ":", "LAMOST, [Fe/H] > $-0.3$"), (NAT, "Earth clock", GREEN, "-", "Nataf+24 primary sample"), (FID, "hard steps n=3", ORANGE, "-", "LAMOST, hard steps $n=3$")] base_p = None for key, model, col, ls, lab in curves: if key not in dists: continue _, phi, fp = dists[key] cf = np.array([F.clock_factor(phi, t, model) for t in ts]) neff = n_sun[0] * eta * fp * cf a.plot(v, np.where(neff > 0, neff, np.nan), color=col, ls=ls, lw=1.6, label=lab) p = np.where(neff > 0, F.POISSON95 / np.maximum(neff, 1e-30), np.nan) if key == FID and model == "Earth clock": base_p = p ok = cf > 0 lo = F.POISSON95 / np.maximum(n_sun[2] * eta_hi * fp * cf, 1e-30) hi = F.POISSON95 / np.maximum(n_sun[1] * eta_lo * fp * cf, 1e-30) b.fill_between(v[ok], lo[ok], hi[ok], color=BLUE, alpha=0.12, lw=0, label="conditional limit over the extremes of\n" "$N_{\\odot-like}$ and $\\eta_\\oplus$") b.plot(v, p, color=BLUE, lw=1.8, label="conditional on fiducial $N_{\\rm eff}$: " "$p < 3.0/N_{\\rm eff}$") elif model != "Earth clock" or key != NAT: b.plot(v, p, color=col, ls=ls, lw=1.3, label=lab) marg = _csv("fermi_marginal_limit.csv") if marg is not None and base_p is not None: for k, ls in ((marg.factor.min(), "-"), (marg.factor.max(), "--")): sig = float(marg.sigma_ln[marg.factor == k].iloc[0]) b.plot(v, base_p * k, color=RED, ls=ls, lw=1.4, label="$N_{\\rm eff}$ uncertainty marginalised, " "$\\sigma_{\\ln N}=%.1f$" % sig) for axx in (a, b): axx.set_xscale("log"); axx.set_yscale("log") axx.set_xlabel("expansion speed $v/c$") with np.errstate(divide="ignore"): top = axx.secondary_xaxis("top", functions=( lambda x: 1e-4 / np.maximum(x, 1e-300), lambda x: 1e-4 / np.maximum(x, 1e-300))) top.set_xlabel("Galaxy-settling time [Gyr] ($10^5$ ly / $v$)") axx.set_xlim(8e-6, 0.1) axx.grid(True, which="major", color="0.92", lw=0.6) axx.set_axisbelow(True) a.set_ylim(1e8, 1e10) a.set_ylabel("$N_{\\rm eff}$: habitable-zone rocky planets whose\n" "settlement wave would already be complete") a.set_title("(a) Size of the Fermi clock (1.07e10 Sun-like stars, " "$\\eta_\\oplus$ = 0.37)", fontsize=10, loc="left", pad=8) a.legend(fontsize=7.2, loc="lower right") b.set_ylim(1e-10, 1e-6) b.set_ylabel("95% upper limit on per-planet probability $p$") b.set_title("(b) Limit on $p$ if no settlement has occurred", fontsize=10, loc="left", pad=8) b.legend(fontsize=7.0, loc="upper right") fig.tight_layout() fig.savefig(FIG + "/fig4_fermi_clock.png") plt.close(fig) def fig5(adopted): """(a) age scales of the catalogues relative to Xiang & Rix; (b) how far each systematic moves the headline number.""" off = pd.read_csv(C.DATA + "/age_scale_offsets.csv") sysd = pd.read_csv(C.DATA + "/xr_systematics.csv") est = pd.read_csv(C.DATA + "/headline_estimates.csv") fig, ax = plt.subplots(1, 2, figsize=(12.4, 4.7), gridspec_kw=dict(width_ratios=[1, 1.25])) a, b = ax pairs = [("Nataf+24 vs X&R22", GREEN, "o", "Nataf+24 (MIST, photometry)"), ("StarHorse vs X&R22", PURPLE, "s", "StarHorse (PARSEC)"), ("Gaia FLAME vs X&R22", SKY, "D", "Gaia DR3 FLAME"), ("this study's Gaia-only ages vs X&R22", ORANGE, "^", "this study's Gaia-only ages (MIST, $G$, BP$-$RP)")] for pair, col, mk, lab in pairs: t = off[(off.pair == pair) & off.subset.str.startswith("X&R age")] if not len(t): continue lo = t.subset.str.extract(r"age ([\d.]+)-([\d.]+)").astype(float) xm = 0.5 * (lo[0] + lo[1]).values allr = off[(off.pair == pair) & (off.subset == "all")].iloc[0] a.plot(xm, t.median_ratio, marker=mk, color=col, lw=1.4, ms=5.5, label=f"{lab}: {allr.median_ratio:.2f} (N = {int(allr.n):,})") a.axhline(1.0, color=INK, lw=0.8) a.axvline(C.T_SUN, color=INK, ls="--", lw=0.8) a.text(C.T_SUN + 0.15, 0.815, "Sun", fontsize=8) a.set_xlabel("Xiang & Rix (2022) age [Gyr]") a.set_ylabel("median age ratio, catalogue / Xiang & Rix") a.set_ylim(0.8, 1.34) a.set_title("(a) Age scales, same stars (legend: overall median ratio)", fontsize=10, loc="left") a.legend(fontsize=7.5, loc="upper left") base = adopted["f_primary"] items = [] for r in sysd.itertuples(): if r.source.startswith("baseline"): continue items.append((r.source, r.f_min, r.f_max, BLUE, True)) e = est.set_index("estimate") def pick(label, names, col, filled=True): v = [e.loc[n, "f_older"] for n in names if n in e.index] if v: items.append((label, min(v + [base]) if len(v) == 1 else min(v), max(v + [base]) if len(v) == 1 else max(v), col, filled)) pick("age catalogue swapped star by star (X&R -> Nataf+24)", [n for n in e.index if "substituted" in n], BLUE) pick("other sample: Nataf+24 primary (four age scales; > 3 Gyr)", [n for n in e.index if n.startswith("Nataf+24 primary sample")], GREEN) pick("bracketing sample: Nataf+24 primary + annex (own / X&R scale)", [n for n in e.index if n.startswith("Nataf+24 primary + metal")], GREEN) pick("bracketing sample: Gaia 200 pc (three published scales)", [n for n in e.index if n.startswith("Gaia 200 pc, ages on")], ORANGE) y = np.arange(len(items))[::-1] b.axvspan(adopted["f_range_lo"], adopted["f_range_hi"], color=BLUE, alpha=0.08, lw=0, label="full range") b.axvspan(adopted["f_core_lo"], adopted["f_core_hi"], color=BLUE, alpha=0.12, lw=0, label="core range") b.axvline(base, color=INK, lw=0.9, label="primary estimate (%.3f)" % base) for yy, (lab, lo, hi, col, filled) in zip(y, items): disf = ("star-count" in lab) or ("Gaia-only ages (x0.87" in lab) b.plot([lo, hi], [yy, yy], color=col, lw=5 if not disf else 2, solid_capstyle="butt", alpha=1.0 if not disf else 0.6) b.text(max(hi, base) + 0.004, yy, f"{lo:.3f} - {hi:.3f}" if hi - lo > 5e-4 else f"{lo:.3f}", va="center", fontsize=7.2) b.set_yticks(y) b.set_yticklabels([i[0] for i in items], fontsize=7.2) b.set_xlim(0.62, 0.835) b.set_ylim(-0.7, len(items) - 0.3) b.set_xlabel("fraction of Sun-like stars older than the Sun") b.set_title("(b) What moves the headline number (thin = disfavoured or " "not a published age scale)", fontsize=10, loc="left") b.tick_params(axis="y", length=0) b.legend(fontsize=7.5, loc="upper center", ncol=3, bbox_to_anchor=(0.5, -0.14)) fig.tight_layout() fig.savefig(FIG + "/fig5_systematics.png") plt.close(fig) def main(): dists = F.distributions() summ = pd.read_csv(C.DATA + "/headstart_summary.csv") est = pd.read_csv(C.DATA + "/headline_estimates.csv") adopted = pd.read_csv(C.DATA + "/headline_adopted.csv").iloc[0].to_dict() fig1(dists, adopted) fig2(summ, est, adopted) fig3() fig4(dists) fig5(adopted) print("[figures] wrote 5 figures to", FIG) if __name__ == "__main__": main()