#!/usr/bin/env python3 """Pillar 1b (added in revision): column SFH and head-start distribution from the Nataf et al. (2024, ApJ 976, 87) subgiants. Why this sample. It differs from the LAMOST / Xiang & Rix sample in the three respects that matter for the systematics the referee raised: * selection: an all-sky box in the Gaia colour - absolute-magnitude diagram (no spectroscopic targeting); * ages: MIST v1.2 fits to broadband photometry + Gaia parallaxes (Xiang & Rix use YY isochrones with spectroscopic Teff, [Fe/H], [alpha/Fe]); since the tracer weights G(tau, [Fe/H]) are also computed on MIST, selection box and ages are on the same models here; * vertical kinematics: Gaia DR3 radial velocities and astrometry give (z, v_z) directly, so the vertical energy needs no action-to-energy conversion and no star gets a clipped velocity. Selection handling (every cut is counted in data/nataf_cutflow.csv) ------------------------------------------------------------------ * "primary sample" (flag PS): colour-magnitude box 0.5 <= W_G <= 2.5, W_G <= 2.5 + 5 ((BP-RP)_0 - 0.9), (BP-RP)_0 <= 1.05 with W_G = M_G - 1.90 (BP-RP). This box can be evaluated exactly on the isochrones. It was designed for -0.5 < [Fe/H] < +0.5; more metal-poor subgiants fall mostly outside it (G/I0 drops from 0.04 at solar metallicity to 0.02 at [Fe/H] = -1), so they are traced with larger weights by fewer, more evolved stars. * "metal-poor annex" (flag MPannex): the region of the diagram on the blue side of the primary box, cleaned of metal-rich turn-off stars with ultraviolet colour cuts. The UV cuts cannot be evaluated on the isochrones used here; they are modelled as a completeness ramp in the catalogue metallicity, c = 1 for [Fe/H] <= -0.5 falling linearly to 0 at -0.25 (the annex's 84th/95th/99th percentiles of [Fe/H] are -0.35, -0.25, -0.16). Annex stars are kinematically cold for their catalogue ages (rms v_z 20-23 km/s at 8-12 Gyr against 27-32 km/s in the primary sample and 27-39 km/s in LAMOST), which points to contamination by UV-bright thin-disc turn-off stars; the annex is therefore used only in a sensitivity variant ("primary + annex"), not in the default. * the catalogue's own quality flags, as in Nataf+24's analysis sample: Blend = 0, VarExcess < 0.95, and no radial-velocity-noise p-value <= 0.001; * distance shell, solar annulus and slab as for the LAMOST sample; * a Gaia radial velocity (needed for v_z); * age >= 2 Gyr (default, as for LAMOST). The primary box only fills up between 2 and 3 Gyr (G/I0 rises from 0.012 to 0.038 at solar metallicity), so a 3 Gyr limit, with the young-fraction prior for 3 Gyr taken from the published SFHs, is run as a variant. """ import os import numpy as np import pandas as pd import common as C import sgweights as S import tracer as T import fetch_data as F import published_sfh as P R0 = 8.178 SHELL = (0.3, 0.8) Z0 = 200.0 DR = 0.5 G_FLOOR = 0.01 TAU_MIN_N = T.TAU_MIN AGE_ERR_FLOOR = 0.05 # fractional; catalogue percentiles can be very narrow UV_RAMP = (-0.50, -0.25) # annex completeness: 1 below, 0 above, linear between def _cw(iso, dc=0.0, dW=0.0): c = iso[:, 9].astype(float) - iso[:, 10].astype(float) + dc return c, iso[:, 8].astype(float) - 1.90 * c + dW def ps_box(iso, feh, dc=0.0, dW=0.0): """Nataf+24 primary-sample box on a MIST isochrone table.""" c, W = _cw(iso, dc, dW) return (W <= 2.5) & (W >= 0.5) & (W <= 2.5 + 5 * (c - 0.9)) & (c <= 1.05) def annex_box(iso, feh, dc=0.0, dW=0.0): """Colour-magnitude region of the metal-poor annex (without UV cuts).""" c, W = _cw(iso, dc, dW) return ((W > 2.5 + 5 * (c - 0.9)) & (W >= 1.5 - 2.5 * (c - 0.3)) & (W <= 1.5 + 3.625 * (c - 0.3)) & (W <= 2.95 - 3.75 * (c - 0.7))) def young_prior(tau_min): """Fraction of all stars formed in the last tau_min Gyr according to the published column SFHs (mean and scatter of the seven solutions). For tau_min = 2 Gyr this is the 0.10 +/- 0.04 used for the LAMOST sample.""" if abs(tau_min - T.TAU_MIN) < 1e-9: return T.PHI_YOUNG v = [P.young_fraction(phi, tau_min) for kind, phi in P.all_published().values() if kind == "column"] return float(np.mean(v)), float(np.std(v, ddof=1)) class Tracer: """G(tau, [Fe/H]) for the primary box, optionally plus the annex region times the UV completeness ramp.""" def __init__(self, union=False, ramp=UV_RAMP, dc=0.0, dW=0.0, slope=C.IMF_SLOPE): self.union, self.ramp, self.slope = union, ramp, slope self.ps = S.TracerWeights(lambda i, f: ps_box(i, f, dc, dW), min_eep=202, slope=slope) self.ann = (S.TracerWeights(lambda i, f: annex_box(i, f, dc, dW), min_eep=202, slope=slope) if union else None) def uv_completeness(self, feh): f1, f2 = self.ramp return np.clip((f2 - np.asarray(feh, float)) / (f2 - f1), 0.0, 1.0) def g(self, age, feh): g = self.ps.g(age, feh) if self.union: g = g + self.uv_completeness(feh) * self.ann.g(age, feh) return g def m_turnoff(self, age, feh): return self.ps.m_turnoff(age, feh) def available(): return os.path.exists(F.NATAF_FITS) and os.path.exists(F.NATAF_GAIA) def load(): from astropy.table import Table t = Table.read(F.NATAF_FITS).to_pandas() t["Sample"] = t.Sample.str.decode("ascii").str.strip() g = pd.read_csv(F.NATAF_GAIA) info = dict(n_catalogue=len(t), n_primary=int((t.Sample == "PS").sum()), n_gaia_lookup=len(g)) m = t.merge(g, left_on="dr3_source_id", right_on="source_id", how="inner", suffixes=("", "_g")) m["d"] = 1.0 / m.parallax # kpc l, b = np.radians(m.l.values), np.radians(m.b.values) m["Zpc"] = 1000.0 * m.d * np.sin(b) + C.Z_SUN m["R"] = np.hypot(R0 - m.d * np.cos(b) * np.cos(l), m.d * np.cos(b) * np.sin(l)) m["Age"] = 10 ** (m.age - 9.0) m["e_lo"] = m.Age - 10 ** (m.age_16 - 9.0) m["e_hi"] = 10 ** (m.age_84 - 9.0) - m.Age rv = m.radial_velocity.notna().values W = np.full(len(m), np.nan) _, _, w_, _ = C.galactic_uvw(m.l[rv].values, m.b[rv].values, m.parallax[rv].values, m.pmra[rv].values, m.pmdec[rv].values, m.radial_velocity[rv].values, m.ra_g[rv].values, m.dec_g[rv].values) W[rv] = w_ + C.W_SUN m["W"] = W m["quality"] = ((m.Blend == 0) & (m.VarExcess < 0.95) & ~(m.Chance_pvalue <= 0.001)) return m, info def select(m, tr=None, shell=SHELL, z0=Z0, dR=DR, tau_min=TAU_MIN_N, quality=True, need_rv=True, zmin=0.0): """tr: Tracer; with tr.union the annex stars whose UV completeness exceeds 0.1 are included, otherwise the primary sample only.""" k = ((m.Age >= tau_min) & (m.d > shell[0]) & (m.d < shell[1]) & (np.abs(m.R - R0) < dR) & (np.abs(m.Zpc) < z0) & (np.abs(m.Zpc) >= zmin)) ps = (m.Sample == "PS").values if tr is not None and tr.union: k &= ps | (tr.uv_completeness(m.feh.values) > 0.1) else: k &= ps if quality: k &= m.quality if need_rv: k &= m.W.notna() return m[k].reset_index(drop=True) def cutflow(m, info): rows = [dict(step="Nataf+ 2024 catalogue", n=info["n_catalogue"], reason="TotalSample.fits (authors' repository)"), dict(step="primary sample (PS)", n=info["n_primary"], reason="colour-magnitude box that can be evaluated on " "isochrones; the metal-poor annex needs UV colours"), dict(step="pre-selected for the Gaia look-up (PS + annex)", n=info["n_gaia_lookup"], reason="catalogue-based distance; slab surroundings and " "the layers above the slab")] k = (m.Sample == "PS").values steps = [("primary sample among the looked-up stars", k, ""), ("0.3 < d < 0.8 kpc (Gaia parallax)", (m.d > SHELL[0]) & (m.d < SHELL[1]), "same shell as for LAMOST; age mix flat with distance"), ("|R - R0| < 0.5 kpc", np.abs(m.R - R0) < DR, "solar annulus"), ("|Z| < 200 pc", np.abs(m.Zpc) < Z0, "slab"), ("catalogue quality flags", m.quality, "Blend = 0, VarExcess < 0.95, RV-noise p-value not <= 0.001"), ("Gaia radial velocity available", m.W.notna(), "needed for the vertical velocity"), ("age >= 2 Gyr (sample used)", m.Age >= TAU_MIN_N, "tracer limit"), ("of which: age >= 3 Gyr", m.Age >= 3.0, "variant; the primary box only fills up between 2 and 3 Gyr")] acc = np.ones(len(m), bool) for name, cut, why in steps: acc &= np.asarray(cut) rows.append(dict(step=name, n=int(acc.sum()), reason=why)) return pd.DataFrame(rows) class Sample: def __init__(self, s, tr, pot_name="MW14", z0=Z0, mode="vT", zcap=None): self.s, self.tr, self.z0, self.mode, self.zcap = s, tr, z0, mode, zcap self.pot = C.get_potential(pot_name) self.age = s.Age.values.astype(float) self.e_lo = np.maximum(s.e_lo.values, AGE_ERR_FLOOR * self.age) self.e_hi = np.maximum(s.e_hi.values, AGE_ERR_FLOOR * self.age) self.zpc = s.Zpc.values self.I0 = S.imf_integral(0.8, 1.2, tr.slope) self.g = np.maximum(tr.g(self.age, s.feh.values) / self.I0, G_FLOOR) self.E = self.pot.energy(self.zpc, s.W.values) self.fcap = (1.0 if zcap is None else self.pot.slab_fraction(self.E, zcap)) self.alive = S.alive_fraction(tr, self.age, s.feh.values) def weights(self, idx=None, living=False): idx = np.arange(len(self.age)) if idx is None else idx wd = T.vertical_weights(self.pot, self.E[idx], self.zpc[idx], self.z0, self.mode) w = wd / self.g[idx] if self.zcap is not None: w = w * self.fcap[idx] return w * self.alive[idx] if living else w def shell_test(m, tr): I0 = S.imf_integral(0.8, 1.2) rows = [] for lo, hi in [(0.15, 0.3), (0.3, 0.45), (0.45, 0.6), (0.6, 0.8), (0.8, 1.0), (1.0, 1.3)]: s = select(m, tr, (lo, hi), 300.0, need_rv=False) if len(s) < 50: continue w = 1 / np.maximum(tr.g(s.Age.values, s.feh.values) / I0, G_FLOOR) a = s.Age.values rows.append(dict(d_lo=lo, d_hi=hi, n=len(s), f_gt_sun=np.sum(w * (a > C.T_SUN)) / w.sum(), f_gt_8=np.sum(w * (a > 8)) / w.sum(), rv_fraction=float(s.W.notna().mean()))) return pd.DataFrame(rows) def kinematics_by_age(m, tr): """rms vertical velocity by catalogue age (a check on age quality: the age - velocity-dispersion relation flattens when ages are noisy).""" s = select(m, tr) rows = [] for lo, hi in [(2, 3), (3, 4.567), (4.567, 6), (6, 8), (8, 10), (10, 12), (12, 20)]: k = ((s.Age >= lo) & (s.Age < hi)).values if k.sum() < 20: continue rows.append(dict(age_lo=lo, age_hi=hi, n=int(k.sum()), rms_vz_kms=float(np.sqrt(np.mean(s.W.values[k] ** 2))), annex_fraction=float((s.Sample.values[k] != "PS").mean()))) return pd.DataFrame(rows) def layer_test(m, tr, layers=((200.0, 350.0), (350.0, 500.0), (500.0, 750.0)), z0=Z0): """As analysis_xr.layer_test: the orbits of the slab stars (here from Gaia velocities directly) predict the age mix in layers above the slab, which is compared with the mix the catalogue contains there.""" slab = select(m, tr, SHELL, z0) sm = Sample(slab, tr) w_dyn, inv_g = sm.weights(), 1.0 / sm.g rows = [] def frac(w, a, cut): return float(np.sum(w * (a > cut)) / w.sum()) def ferr(w, a, cut): f = frac(w, a, cut) return float(np.sqrt(np.sum(w ** 2 * ((a > cut) - f) ** 2)) / w.sum()) def add(lo, hi, model, n, w, a): rows.append(dict(layer_lo=lo, layer_hi=hi, model=model, n=n, f_older=frac(w, a, C.T_SUN), f_older_err=ferr(w, a, C.T_SUN), f_gt8=frac(w, a, 8.0), f_gt8_err=ferr(w, a, 8.0))) add(0.0, z0, "observed", len(slab), inv_g, sm.age) for a, b in layers: o = select(m, tr, SHELL, b, zmin=a, need_rv=False) go = np.maximum(tr.g(o.Age.values, o.feh.values) / sm.I0, G_FLOOR) add(a, b, "observed", len(o), 1 / go, o.Age.values.astype(float)) fl = sm.pot.frac_weighted( sm.E, lambda z: ((np.abs(z) >= a) & (np.abs(z) < b)).astype(float)) add(a, b, "dynamical weights", len(slab), w_dyn * fl, sm.age) return pd.DataFrame(rows) def one(m, tr, name, n_boot=0, pot="MW14", z0=Z0, mode="vT", shell=SHELL, dR=DR, tau_min=TAU_MIN_N, feh_floor=None, living=False, deconvolve=True, phi_young=None, zcap=None, quality=True, cut=None, age_scale=1.0): """phi_young: fraction of all stars formed in the last tau_min Gyr (default: mean of the published column SFHs for that tau_min).""" prior = young_prior(tau_min) phi_young = prior[0] if phi_young is None else phi_young if age_scale != 1.0: m = m.assign(Age=m.Age * age_scale, e_lo=m.e_lo * age_scale, e_hi=m.e_hi * age_scale) s = select(m, tr, shell, z0, dR, tau_min, quality=quality) if cut is not None: s = s[cut(s)].reset_index(drop=True) sm = Sample(s, tr, pot, z0, mode, zcap) keep = None if feh_floor is None else (s.feh.values > feh_floor) w = sm.weights(living=living) st, phi = T.run_variant(sm.age, sm.e_lo, sm.e_hi, w, keep, phi_young, deconvolve=deconvolve, tau_min=tau_min, rescale_young=False) row = dict(variant=name, n=len(s) if keep is None else int(keep.sum()), n_eff=float(w.sum() ** 2 / (w ** 2).sum()), **st) row["annex_weight_share"] = float( np.sum(w * (s.Sample.values != "PS")) / w.sum()) if n_boot: b = T.bootstrap_variant(sm.age, sm.e_lo, sm.e_hi, lambda idx: sm.weights(idx, living), keep, n_boot=n_boot, tau_min=tau_min, phi_young=prior, rescale_young=False) for k, (p16, p50, p84) in b.items(): row[k + "_lo"], row[k + "_hi"] = p16, p84 return row, phi def main(n_boot=40): if not available(): print("[nataf] catalogue or Gaia look-up missing: step skipped " "(run fetch_data.py)") return None import age_scales m, info = load() tr = Tracer() tru = Tracer(union=True) pd.set_option("display.width", 250, "display.max_columns", 30, "display.max_colwidth", 80) cf = cutflow(m, info) cf.to_csv(C.DATA + "/nataf_cutflow.csv", index=False) print(cf[["step", "n"]].to_string(index=False)) st = shell_test(m, tr) st.to_csv(C.DATA + "/nataf_shell_test.csv", index=False) print("[nataf] 1/G-weighted age mix by distance (|Z| < 300 pc):") print(st.round(3).to_string(index=False)) kin = pd.concat([kinematics_by_age(m, tr).assign(sample="primary"), kinematics_by_age(m, tru).assign(sample="primary + annex")]) kin.to_csv(C.DATA + "/nataf_kinematics.csv", index=False) print("[nataf] rms v_z by catalogue age:") print(kin.round(2).to_string(index=False)) lt = layer_test(m, tr) lt.to_csv(C.DATA + "/nataf_layer_test.csv", index=False) print("[nataf] age mix above the slab: observed vs predicted") print(lt.round(3).to_string(index=False)) rows, phis = [], {} def add(name, nb=0, tr_=None, **kw): r, phi = one(m, tr_ or tr, name, nb, **kw) rows.append(r); phis[name] = phi print(f"[nataf] {name:46s} N={r['n']:6d} f_older={r['f_gt_0']:.3f} " f"f>1={r['f_gt_1']:.3f} f>3={r['f_gt_3']:.3f} " f"f>5={r['f_gt_5']:.3f} ={r['mean_headstart_older']:.2f}") return r ratio = age_scales.measured_ratios({"Nataf": 0.95})["Nataf"] add("formed: all [Fe/H]", n_boot) add("formed: [Fe/H]>-0.5", n_boot, feh_floor=-0.5) add("formed: [Fe/H]>-0.3", n_boot, feh_floor=-0.3) add("living MS: all [Fe/H]", living=True) add("no column weights (mid-plane slab)", mode="none") add("ages on the X&R22 scale", n_boot, age_scale=1 / ratio)["age_scale"] = 1 / ratio # the other two published scales, via their ratios to X&R rat = age_scales.measured_ratios({"FLAME": 0.93, "StarHorse": 1.04}) for key, lab in (("FLAME", "ages on the Gaia FLAME scale"), ("StarHorse", "ages on the StarHorse scale")): sc = rat[key] / ratio add(lab, age_scale=sc)["age_scale"] = sc add("primary + annex", n_boot, tr_=tru) add("primary + annex, mid-plane", tr_=tru, mode="none") add("primary + annex, [Fe/H]>-0.5", tr_=tru, feh_floor=-0.5) add("primary + annex, ages on the X&R22 scale", tr_=tru, age_scale=1 / ratio)["age_scale"] = 1 / ratio add("primary + annex, UV ramp -0.6..-0.35", tr_=Tracer(union=True, ramp=(-0.60, -0.35))) add("primary + annex, UV ramp -0.4..-0.15", tr_=Tracer(union=True, ramp=(-0.40, -0.15))) add("tracer ages > 3 Gyr (young prior for 3 Gyr)", n_boot, tau_min=3.0) add("tracer ages > 3 Gyr, primary + annex", tr_=tru, tau_min=3.0) add("tracer ages > 4 Gyr (young prior for 4 Gyr)", tau_min=4.0) add("no error deconvolution", deconvolve=False) add("slab-fraction weights (uniform S(z))", mode="slab") add("column to |z| < 1.3 kpc only", zcap=1300.0) add("without the catalogue quality flags", quality=False) add("RUWE < 1.2 only", cut=lambda s: s.ruwe < 1.2) add("|b| > 15 deg only", cut=lambda s: np.abs(s.b) > 15) add("A_V < 0.25 only", cut=lambda s: s.AV_model < 0.25) add("z0 = 150 pc", z0=150.0) add("z0 = 300 pc", z0=300.0) add("shell 0.3-0.6 kpc", shell=(0.3, 0.6)) add("shell 0.45-0.8 kpc", shell=(0.45, 0.8)) add("|R-R0| < 0.25 kpc", dR=0.25) add("fiducial KG potential", pot="fiducial") add("light-disc potential", pot="light") add("heavy-disc potential", pot="heavy") add("phi_young = 0.06", phi_young=0.06) add("phi_young = 0.15", phi_young=0.15) for lab, kw in [("isochrone colours +0.02 mag", dict(dc=0.02)), ("isochrone colours -0.02 mag", dict(dc=-0.02)), ("isochrones +0.1 mag in W_G", dict(dW=0.1)), ("isochrones -0.1 mag in W_G", dict(dW=-0.1))]: add(lab, tr_=Tracer(**kw)) pd.DataFrame(rows).to_csv(C.DATA + "/nataf_results.csv", index=False) out = pd.DataFrame({"age_lo": C.AGE_EDGES[:-1], "age_hi": C.AGE_EDGES[1:]}) for k in ["formed: all [Fe/H]", "formed: [Fe/H]>-0.5", "formed: [Fe/H]>-0.3", "living MS: all [Fe/H]", "no column weights (mid-plane slab)", "ages on the X&R22 scale", "primary + annex", "primary + annex, [Fe/H]>-0.5"]: out[k] = phis[k] out.to_csv(C.DATA + "/nataf_age_distribution.csv", index=False) return rows, phis if __name__ == "__main__": main()