"""Stage functions used by run_all.py (abundance summary, Earth calibration, population runs, sensitivities, hyper-parameter Monte Carlo, ranks).""" import os import numpy as np import pandas as pd import abundances as ab import mdwarf_earliness_snapshot as E import population as pop import rank import thermal from common import DATA_DIR, RAW_DIR, AGE_SUN, AGE_EARTH, SEED N_VALUES = (1, 2, 3, 4, 6, 8) # end-of-habitability rules: ridge-law degassing threshold (baseline), plate-speed-law degassing threshold # (variant law; identical to the baseline for the stagnant lid) and the end of dry decompression melting CRITERIA = {"degassing>=10%": "lifetime", "platespeed>=10%": "lifetime_arc", "any_melting": "melt_end"} MOBILE = ("F1_classical", "F2h_hybrid", "F2_flat") # --------------------------------------------------------------------------------------------- # abundances # --------------------------------------------------------------------------------------------- def _ols(x, y): x, y = np.asarray(x, dtype=float), np.asarray(y, dtype=float) A = np.column_stack([np.ones_like(x), x]) c, *_ = np.linalg.lstsq(A, y, rcond=None) r = y - A @ c cov = np.linalg.inv(A.T @ A) * r.var(ddof=2) return float(c[1]), float(np.sqrt(cov[1, 1])), float(np.sqrt(np.mean(r ** 2))) def abundance_stage(n_boot=300): S, table, summary = ab.analyse("MgSi", verbose=False, n_boot=n_boot) emp = pop.hyper_from_fits(summary) tw, dm, bb = S["solar_twins"], S["harps_gto"], S["battistini_bensby"] # test of the one-zone decay model: [Th/Eu] at birth versus age (solar twins) ft = summary["fit_theu"] emp["theu_slope_obs"], emp["theu_slope_err"] = ft["b"], ft["e_b"] emp["theu_slope_model"] = summary["theu_model_slope"] emp["ism_strength_from_theu"] = ft["b"] / summary["theu_model_slope"] emp["ism_strength_from_theu_err"] = ft["e_b"] / summary["theu_model_slope"] # bivariate fits: y = a + b (age - 4.57) + c [Fe/H]; Eu rows are DATA ONLY (no model decay factor) rows = [] for label, t, use, y, e in (("Th direct, solar twins (Botelho+2019)", tw, "use_th", "logf_th_direct", "e_logf_th_direct"), ("Eu, solar twins (Spina+2018)", tw, "use_eu", "logf_eu_rock", "e_logf_eu_rock"), ("Eu, HARPS-GTO (Delgado Mena+2017)", dm, "use_eu", "logf_eu_rock", "e_logf_eu_rock"), ("Eu, F/G dwarfs (Battistini & Bensby 2016)", bb, "use_eu", "logf_eu_rock", "e_logf_eu_rock")): m = t[use] f = ab.fit_plane_scatter(t.loc[m, "age"], t.loc[m, "feh"], t.loc[m, y], t.loc[m, e], x10=AGE_SUN, n_boot=n_boot) pct = np.percentile(t.loc[m, y], [2.5, 16, 50, 84, 97.5]) rows.append({"sample": label, "n": f["n"], "zero_point": f["a"], "e_zero_point": f["e_a"], "age_coeff_dex_per_gyr": f["b"], "e_age_coeff": f["e_b"], "feh_coeff_dex_per_dex": f["c"], "e_feh_coeff": f["e_c"], "intrinsic_scatter_dex": f["s_int"], "e_scatter": f["e_s"], "age_min": t.loc[m, "age"].min(), "age_max": t.loc[m, "age"].max(), "feh_min": t.loc[m, "feh"].min(), "feh_max": t.loc[m, "feh"].max(), "p2.5": pct[0], "p16": pct[1], "p50": pct[2], "p84": pct[3], "p97.5": pct[4]}) biv = pd.DataFrame(rows) biv.to_csv(os.path.join(DATA_DIR, "abundance_bivariate_fits.csv"), index=False) # RAW age slopes by reference element (ordinary least squares; no model factor of any kind) rows = [] sel = [("Eu, solar twins", tw, tw["use_eu"]), ("Eu, HARPS-GTO", dm, dm["use_eu"]), ("Eu, HARPS-GTO thin disc only", dm, dm["use_eu"] & (dm["pop"] == "thin")), ("Eu, HARPS-GTO thick disc and high-alpha", dm, dm["use_eu"] & (dm["pop"] != "thin")), ("Eu, Battistini & Bensby", bb, bb["use_eu"])] for label, t, m in sel: if m.sum() < 10: continue r = {"sample": label, "n": int(m.sum())} for col, tag in (("logf_eu_rock", "eu_rock"), ("eu_mg", "eu_mg"), ("eu_si", "eu_si")): b, e, rms = _ols(t.loc[m, "age"] - AGE_SUN, t.loc[m, col]) r[f"slope_{tag}"], r[f"e_slope_{tag}"], r[f"rms_{tag}"] = b, e, rms b, e, _ = _ols(t.loc[m, "age"] - AGE_SUN, t.loc[m, "eu_h"] - t.loc[m, "feh"]) r["slope_eu_fe"], r["e_slope_eu_fe"] = b, e rows.append(r) m = tw["use_th"] r = {"sample": "Th direct (decay-corrected inside the star), solar twins", "n": int(m.sum())} for ref, tag in ((0.5 * (tw["mg_h"] + tw["si_h"]), "eu_rock"), (tw["mg_h"], "eu_mg"), (tw["si_h"], "eu_si")): y = tw["th_h"] - ref + ab.DEX_PER_GYR_TH * (tw["age"] - AGE_SUN) b, e, rms = _ols(tw.loc[m, "age"] - AGE_SUN, y[m]) r[f"slope_{tag}"], r[f"e_slope_{tag}"], r[f"rms_{tag}"] = b, e, rms rows.append(r) raw = pd.DataFrame(rows) raw.to_csv(os.path.join(DATA_DIR, "abundance_raw_slopes.csv"), index=False) # how robust is the "intrinsic" scatter to the adopted measurement errors? (errors scaled by 1, 0.8, 0.6) rows = [] for label, t, use, y, e, note in ( ("Th direct, solar twins", tw, "use_th", "logf_th_direct", "e_logf_th_direct", "per-star errors from Botelho+2019"), ("Eu, solar twins", tw, "use_eu", "logf_eu_rock", "e_logf_eu_rock", "per-star errors from Spina+2018, Bedell+2018"), ("Eu, HARPS-GTO", dm, "use_eu", "logf_eu_rock", "e_logf_eu_rock", "per-star errors from Delgado Mena+2017"), ("Eu, Battistini & Bensby", bb, "use_eu", "logf_eu_rock", "e_logf_eu_rock", "Eu error ASSUMED 0.08 dex for every star")): m = t[use] r = {"sample": label, "n": int(m.sum()), "error_model": note} for scale in (1.0, 0.8, 0.6): f = ab.fit_line_scatter(t.loc[m, "age"], t.loc[m, y], scale * t.loc[m, e], x0=AGE_SUN, n_boot=20) r[f"s_int_err_x{scale:g}"] = f["s_int"] if scale == 1.0: r["rms_about_trend"], r["median_quoted_error"] = f["rms"], f["med_err"] rows.append(r) scat = pd.DataFrame(rows) scat.to_csv(os.path.join(DATA_DIR, "abundance_scatter_robustness.csv"), index=False) # age-binned medians of the broad Eu samples (population-level trend, data only) broad = pd.concat([dm[dm["use_eu"]], bb[bb["use_eu"]]]) rows = [] for lo, hi in ((0, 3), (3, 6), (6, 9), (9, 14)): m = (broad["age"] >= lo) & (broad["age"] < hi) v = broad.loc[m, "logf_eu_rock"] rows.append({"age_lo": lo, "age_hi": hi, "n": int(m.sum()), "median_eu_rock": v.median(), "sem": v.std() / np.sqrt(m.sum()), "sd": v.std(), "median_eu_mg": broad.loc[m, "eu_mg"].median(), "median_eu_si": broad.loc[m, "eu_si"].median(), "median_feh": broad.loc[m, "feh"].median(), "median_eu_fe": (broad.loc[m, "eu_h"] - broad.loc[m, "feh"]).median()}) binned = pd.DataFrame(rows) binned.to_csv(os.path.join(DATA_DIR, "abundance_age_bins_broad.csv"), index=False) # observed (not deconvolved) [Eu/Mg] ranges, for comparison with the ranges quoted by Nimmo et al. (2020) rows = [] for label, t in (("Spina+2018 solar twins", tw), ("Delgado Mena+2017 HARPS-GTO", dm), ("Battistini & Bensby 2016", bb)): for sel_, m in (("all stars with Eu", t["eu_h"].notna() & t["mg_h"].notna()), ("analysis cuts", t["use_eu"])): v = (t.loc[m, "eu_h"] - t.loc[m, "mg_h"]).values q = np.percentile(v, [1, 10, 50, 90, 99]) rows.append({"sample": label, "selection": sel_, "n": int(m.sum()), "p1": q[0], "p10": q[1], "p50": q[2], "p90": q[3], "p99": q[4], "sd": float(np.std(v)), "median_quoted_error": float(np.nanmedian(np.hypot(t.loc[m, "e_eu_h"], t.loc[m, "e_mg_h"])))}) eumg = pd.DataFrame(rows) eumg.to_csv(os.path.join(DATA_DIR, "eu_mg_observed_ranges.csv"), index=False) # choice of rock reference (Mg, Si or their mean): zero point AND slope rows = [] for ref in ("MgSi", "Mg", "Si"): _, _, sm = (S, table, summary) if ref == "MgSi" else ab.analyse(ref, verbose=False, n_boot=60) e2 = pop.hyper_from_fits(sm) rows.append({"rock_reference": ref, "zero_point": e2["a"], "zero_point_sd_between_samples": e2["e_a"], "raw_age_slope_broad": e2["b"], "raw_age_slope_harps": sm["fits_eu_raw"]["harps_gto"]["b"], "raw_age_slope_bb": sm["fits_eu_raw"]["battistini_bensby"]["b"], "raw_age_slope_twins": sm["fits_eu_raw"]["solar_twins"]["b"], "slope_with_full_decay_model_broad": 0.5 * (sm["fits_eu"]["harps_gto"]["b"] + sm["fits_eu"]["battistini_bensby"]["b"]), "intrinsic_scatter": e2["s"], "th_direct_zero_point": sm["fit_th_age"]["a"], "th_direct_slope": sm["fit_th_age"]["b"], "th_direct_scatter": sm["fit_th_age"]["s_int"]}) refs = pd.DataFrame(rows) refs.to_csv(os.path.join(DATA_DIR, "abundance_rock_reference_check.csv"), index=False) # Mishenina et al. 2022 thorium (no ages, no per-star errors): scatter and metallicity trend only mi = pd.read_csv(os.path.join(RAW_DIR, "mishenina2022_tablea1.csv")) mi = mi[mi["[Th/Fe]"].notna() & (mi["[Fe/H]"] > -1.0) & mi["Teff"].between(5300, 6400)] sl, ic = np.polyfit(mi["[Fe/H]"], mi["[Th/Fe]"], 1) mish = {"n": int(len(mi)), "mean_th_fe": float(mi["[Th/Fe]"].mean()), "sd_th_fe": float(mi["[Th/Fe]"].std()), "slope_vs_feh": float(sl), "sd_about_trend": float((mi["[Th/Fe]"] - ic - sl * mi["[Fe/H]"]).std())} return {"S": S, "table": table, "summary": summary, "emp": emp, "bivariate": biv, "binned": binned, "rock_refs": refs, "mishenina": mish, "eu_mg": eumg, "raw_slopes": raw, "scatter": scat} def heat_budget_stage(emp, n=20000, seed=SEED + 5): """Initial radiogenic heat production per unit rock mass relative to the Earth's (94 TW at formation for the bulk silicate Earth), (a) for the population model at several formation epochs and three strengths of the one-zone decay model (0 = data only, 0.5 = baseline, 1 = full model), and (b) star by star for the Eu samples, binned by age and by metallicity. Only the r-process-per-rock level is measured; the nuclide-specific factors of Th, U and 40K are the decay model and the K prior.""" h0 = thermal._H0_ARR # initial Earth heat by nuclide (238U, 235U, 232Th, 40K) lam = thermal._LAM h45 = h0 * np.exp(-lam * AGE_EARTH) z = pop.draw_normals(n, seed) rows = [] for sc in (0.0, pop.ISM_SCALE_BASE, 1.0): h = base_hyper(emp, ism_scale=sc) for age in (0.0, 2.0, AGE_SUN, 6.0, 8.5, 10.0, 12.0): for label, vary in (("measured r-process scatter only (K fixed at median)", ("th",)), ("r-process scatter + K scatter and volatility prior", ("th", "k", "vol"))): fac = pop.hpe_factors(age, z, h, vary) tot0 = (fac * h0).sum(axis=1) / h0.sum() tot45 = (fac * h45).sum(axis=1) / h45.sum() thu0 = (fac[:, :3] * h0[:3]).sum(axis=1) / h0[:3].sum() q0, q45 = np.percentile(tot0, [2.5, 16, 50, 84, 97.5]), np.percentile(tot45, [2.5, 16, 50, 84, 97.5]) rows.append({"decay_model_strength": sc, "stellar_age_today_gyr": age, "scatter_included": label, "median_f_r_data_only": float(10.0 ** (h["a"] + h["b"] * (min(age, pop.AGE_CLIP) - AGE_SUN))), "median_f_Th": float(np.median(fac[:, 2])), "median_f_U238": float(np.median(fac[:, 0])), "median_f_U235": float(np.median(fac[:, 1])), "median_f_K40": float(np.median(fac[:, 3])), "ThU_heat_initial_rel_median": float(np.median(thu0)), **{f"heat_initial_rel_p{k}": v for k, v in zip(("2.5", "16", "50", "84", "97.5"), q0)}, **{f"heat_at_4.5Gyr_rel_p{k}": v for k, v in zip(("2.5", "16", "50", "84", "97.5"), q45)}}) model_tab = pd.DataFrame(rows) model_tab.to_csv(os.path.join(DATA_DIR, "initial_heat_budget_model.csv"), index=False) # star by star (measured Eu per rock; own age for the model factors at baseline strength; median K relation) h = base_hyper(emp) sc = h["ism_scale"] st = pd.read_csv(os.path.join(DATA_DIR, "hpe_stars_MgSi.csv")) st = st[st["use_eu"] & st["sample"].isin(["harps_gto", "battistini_bensby"])].copy() age_c = np.clip(st["age"].values, 0.0, pop.AGE_CLIP) lf = st["logf_eu_rock"].values f_th = 10.0 ** (lf + sc * ab.ism_factor("Th232", age_c)) f_u8 = 10.0 ** (lf + sc * ab.ism_factor("U238", age_c)) f_u5 = 10.0 ** (lf + sc * ab.ism_factor("U235", age_c)) f_k = 10.0 ** (h["k_a"] + h["k_b"] * (age_c - AGE_SUN) + sc * ab.ism_factor("K40", age_c)) fac = np.column_stack([f_u8, f_u5, f_th, f_k]) st["heat_initial_rel"] = (fac * h0).sum(axis=1) / h0.sum() st["heat_at_4.5Gyr_rel"] = (fac * h45).sum(axis=1) / h45.sum() st["f_eu_data_only"] = 10.0 ** lf st["f_th"] = f_th rows = [] bins = [("age", lo, hi) for lo, hi in ((0, 3), (3, 6), (6, 9), (9, 14))] + \ [("feh", lo, hi) for lo, hi in ((-1.0, -0.4), (-0.4, -0.15), (-0.15, 0.05), (0.05, 0.25), (0.25, 0.6))] for col, lo, hi in bins: m = (st[col] >= lo) & (st[col] < hi) if m.sum() < 5: continue r = {"binned_by": "stellar age [Gyr]" if col == "age" else "[Fe/H]", "lo": lo, "hi": hi, "n": int(m.sum()), "median_age": float(st.loc[m, "age"].median()), "median_feh": float(st.loc[m, "feh"].median())} for c in ("f_eu_data_only", "f_th", "heat_initial_rel", "heat_at_4.5Gyr_rel"): q = np.percentile(st.loc[m, c], [16, 50, 84]) r.update({f"{c}_p16": q[0], f"{c}_median": q[1], f"{c}_p84": q[2]}) rows.append(r) star_tab = pd.DataFrame(rows) star_tab.to_csv(os.path.join(DATA_DIR, "initial_heat_budget_stars.csv"), index=False) return model_tab, star_tab # --------------------------------------------------------------------------------------------- # Earth calibration and histories # --------------------------------------------------------------------------------------------- def earth_stage(): rows = [] for name in thermal.FORMULATIONS: m = thermal.make_model(name) r = m.earth_summary() rows.append({"formulation": name, **r}) for beta in (0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 1.0 / 3.0): m = thermal.make_model("power", beta=beta) rows.append({"formulation": f"power beta={beta:.2f}", **m.earth_summary()}) cal = pd.DataFrame(rows) cal.to_csv(os.path.join(DATA_DIR, "earth_calibration.csv"), index=False) hist = [] rec = np.round(np.arange(0.0, 30.001, 0.1), 3) one = lambda: {"fac": np.ones((1, 4)), "mass": np.ones(1), "eta_fac": np.ones(1)} # noqa: E731 for name in thermal.FORMULATIONS: m = thermal.make_model(name) res = thermal.integrate(m, np.array([m.T_init]), one(), 30.0, dt=0.005, record=rec) hist.append(pd.DataFrame({"formulation": name, "t_gyr": res["t_rec"], "T_p_K": res["T_rec"][:, 0], "degassing_rel": res["D_rec"][:, 0]})) hist = pd.concat(hist, ignore_index=True) hist.to_csv(os.path.join(DATA_DIR, "earth_twin_histories.csv"), index=False) # Earth twin lifetime versus Th (and U) budget fgrid = np.array([0.5, 0.7, 0.85, 1.0, 1.2, 1.4, 2.0, 3.0]) rows = [] for name in thermal.FORMULATIONS: m = thermal.make_model(name) fac = np.ones((len(fgrid), 4)) fac[:, :3] = fgrid[:, None] p = {"fac": fac, "mass": np.ones(len(fgrid)), "eta_fac": np.ones(len(fgrid))} res = thermal.integrate(m, np.full(len(fgrid), m.T_init), p, pop.T_MAX[name], dt=0.01, dt_early=0.0025) for f, a, c, b in zip(fgrid, res["lifetime"], res["lifetime_arc"], res["melt_end"]): rows.append({"formulation": name, "f_Th_U": f, "lifetime_degassing": a, "lifetime_platespeed_law": c, "lifetime_melting": b}) sens = pd.DataFrame(rows) sens.to_csv(os.path.join(DATA_DIR, "lifetime_vs_th_budget.csv"), index=False) return cal, hist, sens # --------------------------------------------------------------------------------------------- # population helpers # --------------------------------------------------------------------------------------------- def base_hyper(emp, **over): h = dict(pop.DEFAULT_HYPER) for k in ("a", "b", "s", "k_a", "k_b", "s_k"): h[k] = emp[k] h.update(over) return h def ranks_for(life_nodes, n_values=N_VALUES, **kw): return {n: rank.rank_from_lifetimes(life_nodes, pop.NODE_AGES, n, **kw) for n in n_values} def alive_today_fraction(life_nodes, sfh="MD14"): """Share of all rocky planets formed so far (weighted by the cosmic star-formation history) whose geological lifetime exceeds their present age.""" t0 = E.age_today() b, psi, past, _ = E._time_grids(sfh) age = t0 - b surv = np.array([(life_nodes[k][None, :] > age[:, None]).mean(axis=1) for k in range(len(pop.NODE_AGES))]) a_c = np.clip(age, pop.NODE_AGES[0], pop.NODE_AGES[-1]) s = np.array([np.interp(a_c[i], pop.NODE_AGES, surv[:, i]) for i in range(len(age))]) return float(np.trapezoid(psi * s, b) / past) def tail_shares(L, n): """Share of sum(L^n) carried by the longest-lived 1%, 5% and 10% of planets.""" w = np.sort(np.asarray(L, dtype=float)) ** n tot = w.sum() return {f"top{int(q * 100)}pct": float(w[int(np.ceil(len(w) * (1 - q))):].sum() / tot) for q in (0.01, 0.05, 0.10)} def baseline_stage(emp, n=5000, seed=SEED): z = pop.draw_normals(n, seed) h = base_hyper(emp) runs, surv_rows, pct_rows, rank_rows, curve_rows, peak_rows, tail_rows = {}, [], [], [], [], [], [] tgrid = np.round(np.arange(0.0, 60.01, 0.25), 2) for name in thermal.FORMULATIONS: r = pop.run_population(name, h, z) runs[name] = r for k, age in enumerate(pop.NODE_AGES): tp = r["T_peak"][k] peak_rows.append({"formulation": name, "thermostat_K": getattr(r["model"], "T_cap", None), "stellar_age_today_gyr": age, "T_peak_median_K": float(np.median(tp)), "T_peak_p95_K": float(np.percentile(tp, 95)), "T_peak_max_K": float(tp.max()), "frac_reaching_thermostat": float(r["at_cap"][k].mean())}) for crit, key in CRITERIA.items(): if name == "F3_stagnant" and crit == "platespeed>=10%": continue for k, age in enumerate(pop.NODE_AGES): L = r[key][k] sv = pop.survival(L) surv_rows.append({"formulation": name, "criterion": crit, "stellar_age_today_gyr": age, **{f"P_hab_{int(a)}Gyr": v for a, v in zip(pop.TARGET_AGES, sv)}}) pc = np.percentile(L, [1, 5, 16, 50, 84, 95, 99]) pct_rows.append({"formulation": name, "criterion": crit, "stellar_age_today_gyr": age, "p1": pc[0], "p5": pc[1], "p16": pc[2], "median": pc[3], "p84": pc[4], "p95": pc[5], "p99": pc[6], "max": L.max(), "frac_censored": float((L >= r["t_max"] - 0.06).mean())}) curve_rows.append(pd.DataFrame({"formulation": name, "criterion": crit, "stellar_age_today_gyr": age, "t_gyr": tgrid, "P_hab": [(L > t).mean() for t in tgrid]})) rk_all = ranks_for(r[key]) rk_k = ranks_for(r[key], m_min_host=0.6) rank_rows.append({"formulation": name, "criterion": crit, "hosts": "all stars", "lid_mix": "all planets this regime", **{f"rank_n{n_}": v for n_, v in rk_all.items()}, "alive_today_fraction": alive_today_fraction(r[key])}) rank_rows.append({"formulation": name, "criterion": crit, "hosts": "> 0.6 Msun", "lid_mix": "all planets this regime", **{f"rank_n{n_}": v for n_, v in rk_k.items()}}) # how tail-dominated is the rank? share of sum(L^n) in the longest-lived planets (solar-epoch # node) and the rank when the longest-lived 5% / 1% of planets are clipped to the 95th / 99th percentile k_sun = int(np.argmin(np.abs(pop.NODE_AGES - AGE_SUN))) row = {"formulation": name, "criterion": crit} for n_ in (2, 4, 6): ts = tail_shares(r[key][k_sun], n_) row.update({f"share_sumLn_n{n_}_{q}": v for q, v in ts.items()}) row[f"rank_n{n_}"] = rk_all[n_] for q in (95, 99): clip = np.minimum(r[key], np.percentile(r[key], q, axis=1)[:, None]) row[f"rank_n{n_}_clipped_p{q}"] = rank.rank_from_lifetimes(clip, pop.NODE_AGES, n_) tail_rows.append(row) # other star-formation histories (same variants as studies/mdwarf_early) for name in thermal.FORMULATIONS: for crit, key in CRITERIA.items(): if name == "F3_stagnant" and crit == "platespeed>=10%": continue for sfh in ("MD14_exp", "Sobral13", "BP15_long"): rk_s = ranks_for(runs[name][key], sfh=sfh) rank_rows.append({"formulation": name, "criterion": crit, "hosts": "all stars", "lid_mix": f"all planets this regime; SFH {sfh}", **{f"rank_n{n_}": v for n_, v in rk_s.items()}}) # mixed lid regimes: a fraction p_mobile of planets follows a mobile-lid formulation, the rest F3 for name in MOBILE: for crit, key in CRITERIA.items(): for p_mob in (0.5, 0.17): key3 = "lifetime" if key == "lifetime_arc" else key life = np.concatenate([runs[name][key], runs["F3_stagnant"][key3]], axis=1) w = np.concatenate([np.full(n, p_mob), np.full(n, 1.0 - p_mob)]) rk = ranks_for(life, weights=w) rank_rows.append({"formulation": name, "criterion": crit, "hosts": "all stars", "lid_mix": f"{p_mob:.2f} mobile + {1 - p_mob:.2f} stagnant (F3)", **{f"rank_n{n_}": v for n_, v in rk.items()}}) surv = pd.DataFrame(surv_rows) pct = pd.DataFrame(pct_rows) rk = pd.DataFrame(rank_rows) curves = pd.concat(curve_rows, ignore_index=True) peaks = pd.DataFrame(peak_rows) tails = pd.DataFrame(tail_rows) surv.to_csv(os.path.join(DATA_DIR, "survival_baseline.csv"), index=False) pct.to_csv(os.path.join(DATA_DIR, "lifetime_percentiles_baseline.csv"), index=False) rk.to_csv(os.path.join(DATA_DIR, "rank_baseline.csv"), index=False) curves.to_csv(os.path.join(DATA_DIR, "survival_curves_baseline.csv"), index=False) peaks.to_csv(os.path.join(DATA_DIR, "peak_temperature_baseline.csv"), index=False) tails.to_csv(os.path.join(DATA_DIR, "rank_tail_dependence.csv"), index=False) # reference ranks with fixed caps (reproduces studies/mdwarf_early) ref_rows = [] for cap in (10.0, 15.0, 20.0, 1.0e5): ref_rows.append({"case": f"fixed cap {cap:g} Gyr, all stars", **{f"rank_n{n_}": rank.rank_fixed_cap(n_, cap) for n_ in N_VALUES}}) ref_rows.append({"case": "no cap, hosts > 0.6 Msun (K dwarfs and up)", **{f"rank_n{n_}": rank.rank_fixed_cap(n_, 1.0e5, m_min_host=0.6) for n_ in N_VALUES}}) ref = pd.DataFrame(ref_rows) ref.to_csv(os.path.join(DATA_DIR, "rank_reference_fixed_caps.csv"), index=False) return {"runs": runs, "survival": surv, "percentiles": pct, "ranks": rk, "curves": curves, "reference": ref, "peaks": peaks, "tails": tails} def variance_stage(emp, n=2000, seed=SEED + 1): """One-at-a-time spread of the lifetime (solar formation epoch) from each group of inputs.""" z = pop.draw_normals(n, seed) h = base_hyper(emp) groups = {"Th+U budget (measured r-process scatter)": ("th",), "K budget and volatility (prior)": ("k", "vol"), "planet mass 0.5-2 (prior)": ("mass",), "mantle viscosity, 0.5 dex (prior)": ("eta",), "initial temperature, 100 K (prior)": ("Ti",), "all together": ("th", "k", "vol", "mass", "eta", "Ti")} rows = [] node = np.array([AGE_SUN]) for name in thermal.FORMULATIONS: model = pop.make_model(name, h) for label, vary in groups.items(): r = pop.run_population(name, h, z, node_ages=node, vary=vary, model=model) for crit, key in CRITERIA.items(): if crit == "platespeed>=10%": continue L = r[key][0] p5, p50, p95 = np.percentile(L, [5, 50, 95]) rows.append({"formulation": name, "criterion": crit, "input_varied": label, "sd_gyr": L.std(), "p5": p5, "median": p50, "p95": p95, "width_5_95_gyr": p95 - p5}) # zero-point shifts of the whole population (not scatter): median lifetime. (The decay-model strength # does not act at the solar epoch; its effect by epoch is in scenarios.csv, median_lifetime_by_epoch.) for label, over in (("zero point of r-process per rock -0.044 dex (1 sd between samples)", {"a": h["a"] - 0.044}), ("zero point of r-process per rock +0.044 dex", {"a": h["a"] + 0.044}), ("zero point +0.1 dex", {"a": h["a"] + 0.1})): h2 = dict(h); h2.update(over) r = pop.run_population(name, h2, z, node_ages=node, model=model) for crit, key in CRITERIA.items(): if crit == "platespeed>=10%": continue L = r[key][0] p5, p50, p95 = np.percentile(L, [5, 50, 95]) rows.append({"formulation": name, "criterion": crit, "input_varied": "SHIFT: " + label, "sd_gyr": L.std(), "p5": p5, "median": p50, "p95": p95, "width_5_95_gyr": p95 - p5}) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA_DIR, "lifetime_variance_decomposition.csv"), index=False) return out # --------------------------------------------------------------------------------------------- # scenarios: what the rank depends on. Each scenario changes ONE assumption (or a stated combination) # relative to the baseline. Keys of a scenario: # over overrides of the population hyper-parameters (mass range, decay strength, thermostat, ...) # f_crit degassing threshold: a number, or ("loguniform", lo, hi) for one threshold per planet # law "lifetime" (ridge law, default) or "lifetime_arc" (plate-speed law) # mix_melt share of planets judged by the melting criterion instead of the degassing threshold # tau "hz" = stellar limit from the habitable-zone exit times instead of the main-sequence lifetime # sfh star-formation history # group "tail" (assumptions about the long-lived tail and the criterion), "heat" (heat budget), # "thermostat", "parameter" (thermal-model parameters), "other" # --------------------------------------------------------------------------------------------- SCENARIOS = { "baseline (0.5-2 Earth masses, ridge law, 10% threshold, decay strength 0.5, thermostat 2200 K)": {"group": "baseline"}, "planet mass 0.5-3 Earth masses": {"over": {"m_hi": 3.0}, "group": "tail"}, "planet mass 0.5-5 Earth masses": {"over": {"m_hi": 5.0}, "group": "tail"}, "5% of planets on the melting criterion": {"mix_melt": 0.05, "group": "tail"}, "10% of planets on the melting criterion": {"mix_melt": 0.10, "group": "tail"}, "threshold per planet log-uniform 3-100% of modern Earth": {"f_crit": ("loguniform", 0.03, 1.0), "group": "tail"}, "threshold 3% of modern Earth for all planets": {"f_crit": 0.03, "group": "tail"}, "threshold 100% of modern Earth for all planets (outer habitable zone)": {"f_crit": 1.0, "group": "tail"}, "plate-speed degassing law instead of ridge law": {"law": "lifetime_arc", "group": "tail"}, "decay-model strength 0 (heat budget from data only)": {"over": {"ism_scale": 0.0}, "group": "heat"}, "decay-model strength 1 (full one-zone model, first-pass baseline)": {"over": {"ism_scale": 1.0}, "group": "heat"}, "no thermostat (first-pass behaviour)": {"over": {"T_cap": None}, "group": "thermostat"}, "thermostat 2100 K": {"over": {"T_cap": 2100.0}, "group": "thermostat"}, "thermostat 2300 K": {"over": {"T_cap": 2300.0}, "group": "thermostat"}, "first-pass configuration (decay strength 1, no thermostat)": {"over": {"ism_scale": 1.0, "T_cap": None}, "group": "other"}, "COMBINED: mass 0.5-5 + plate-speed degassing law": {"over": {"m_hi": 5.0}, "law": "lifetime_arc", "group": "tail"}, "COMBINED low: mass 0.5-5 + 10% on melting criterion": {"over": {"m_hi": 5.0}, "mix_melt": 0.10, "group": "tail"}, "COMBINED lowest: mass 0.5-5 + 10% on melting + decay strength 1": {"over": {"m_hi": 5.0, "ism_scale": 1.0}, "mix_melt": 0.10, "group": "tail"}, "COMBINED highest: mass 0.5-2 + decay strength 0 + thermostat 2100 K": {"over": {"ism_scale": 0.0, "T_cap": 2100.0}, "group": "tail"}, "referee's combination: mass 0.5-5 + thermostat 2300 K + decay strength 0": {"over": {"m_hi": 5.0, "ism_scale": 0.0, "T_cap": 2300.0}, "group": "tail"}, "referee's combination + 10% on melting criterion": {"over": {"m_hi": 5.0, "ism_scale": 0.0, "T_cap": 2300.0}, "mix_melt": 0.10, "group": "tail"}, "stellar limit = habitable-zone exit instead of main-sequence lifetime": {"tau": "hz", "group": "other"}, "tenfold future star formation (Behroozi & Peeples 2015 upper bound)": {"sfh": "BP15_long", "group": "other"}, "no viscosity scatter": {"over": {"s_eta": 0.0}, "group": "parameter"}, "viscosity scatter 1 dex": {"over": {"s_eta": 1.0}, "group": "parameter"}, "activation energy 200 kJ/mol": {"over": {"E": 200.0e3}, "group": "parameter"}, "activation energy 400 kJ/mol": {"over": {"E": 400.0e3}, "group": "parameter"}, "Katz et al. solidus (1359 K at zero pressure)": {"over": {"t_sol0": thermal.T_SOL0_LOW}, "group": "parameter"}, "Foley & Smye solidus (1423 K at zero pressure)": {"over": {"t_sol0": thermal.T_SOL0_HIGH}, "group": "parameter"}, "no continental sequestration of heat sources (f_cc = 0)": {"over": {"f_cc": 0.0}, "group": "parameter"}, "r-process zero point +0.1 dex for every planet": {"da": 0.1, "group": "heat"}, "r-process zero point -0.044 dex (1 sd between samples)": {"da": -0.044, "group": "heat"}, "r-process scatter 0.10 dex (rms upper bound of the broad samples)": {"over": {"s": 0.10}, "group": "heat"}, "no trend of r-process per rock with stellar age (b = 0)": {"over": {"b": 0.0}, "group": "heat"}, "trend of solar twins (direct Th slope -0.0095 dex/Gyr)": {"over": {"b": -0.0095}, "group": "heat"}, } def scenario_stage(emp, n=2000, seed=SEED + 2, formulations=None, scenarios=None): """Ranks and lifetime statistics for every scenario and formulation. Population runs are cached, so scenarios that differ only in post-processing (criterion mix, law, stellar limit, star-formation history) share the same planets.""" formulations = thermal.FORMULATIONS if formulations is None else formulations scenarios = SCENARIOS if scenarios is None else scenarios z = pop.draw_normals(n, seed) k_sun = int(np.argmin(np.abs(pop.NODE_AGES - AGE_SUN))) tau_hz = None cache, rows = {}, [] for label, spec in scenarios.items(): over = dict(spec.get("over", {})) fc = spec.get("f_crit", 0.1) h = base_hyper(emp, **over) h["a"] += spec.get("da", 0.0) ckey = (tuple(sorted((k, str(v)) for k, v in over.items())), str(fc), spec.get("da", 0.0)) fc_val = pop.threshold_draw(n, fc[1], fc[2], seed + 7) if isinstance(fc, tuple) else fc kw = {} if spec.get("tau") == "hz": tau_hz = rank.tau_hz() if tau_hz is None else tau_hz kw["tau"] = tau_hz if "sfh" in spec: kw["sfh"] = spec["sfh"] for name in formulations: if name == "F3_stagnant" and ("T_cap" in over or spec.get("law") == "lifetime_arc" or spec.get("mix_melt") or spec.get("tau") or spec.get("sfh")): continue # no thermostat / plate-speed law for the stagnant lid; the rest is not needed if (name, ckey) not in cache: cache[(name, ckey)] = pop.run_population(name, h, z, f_crit=fc_val) r = cache[(name, ckey)] law = spec.get("law", "lifetime") for crit, key in (("degassing-type", law), ("any_melting", "melt_end")): if crit == "any_melting" and (spec.get("f_crit") is not None or spec.get("law") or spec.get("mix_melt")): continue # these scenarios only change the degassing criterion L = r[key] wts = None if crit == "degassing-type" and spec.get("mix_melt"): p = spec["mix_melt"] L = np.concatenate([r[key], r["melt_end"]], axis=1) wts = np.concatenate([np.full(n, 1.0 - p), np.full(n, p)]) rk = ranks_for(L, n_values=(2, 4, 6), weights=wts, **kw) Ls = L[k_sun] ws = np.full(len(Ls), 1.0 / len(Ls)) if wts is None else wts / wts.sum() order = np.argsort(Ls) cdf = np.cumsum(ws[order]) q = lambda f: float(Ls[order][min(np.searchsorted(cdf, f), len(Ls) - 1)]) # noqa: E731 rows.append({"scenario": label, "group": spec["group"], "formulation": name, "criterion": crit, "median_lifetime_solar_epoch": q(0.5), "p95_lifetime_solar_epoch": q(0.95), "p99_lifetime_solar_epoch": q(0.99), **{f"P_hab_{int(a)}Gyr": float(np.sum(ws * (Ls > a))) for a in pop.TARGET_AGES}, "P_hab_20Gyr_oldest_epoch": float(np.sum(ws * (L[-1] > 20.0))), "frac_censored": float((r[key] >= r["t_max"] - 0.06).mean()), "frac_reaching_thermostat_oldest_epoch": float(r["at_cap"][-1].mean()), "T_peak_median_oldest_epoch_K": float(np.median(r["T_peak"][-1])), "T_peak_p95_oldest_epoch_K": float(np.percentile(r["T_peak"][-1], 95)), "frac_T_peak_gt_2500K_oldest_epoch": float((r["T_peak"][-1] > 2500.0).mean()), "median_lifetime_by_epoch": "/".join(f"{np.median(r[key][k]):.1f}" for k in range(len(pop.NODE_AGES))), **{f"rank_n{n_}": v for n_, v in rk.items()}}) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA_DIR, "scenarios.csv"), index=False) # range of the rank over the scenarios, by formulation: this is the summary quoted in the report rng_rows = [] deg = out[(out["criterion"] == "degassing-type") & (out["group"] != "other")] not100 = ~deg["scenario"].str.startswith("threshold 100%") sets = {"tail and criterion assumptions (incl. baseline)": deg[deg["group"].isin(["baseline", "tail"])], "tail and criterion assumptions without the uniform 100% threshold": deg[deg["group"].isin(["baseline", "tail"]) & not100], "all degassing-type scenarios without the uniform 100% threshold": deg[not100], "heat budget (zero point, trend, decay strength, scatter)": deg[deg["group"].isin(["baseline", "heat"])], "thermostat": deg[deg["group"].isin(["baseline", "thermostat"])], "thermal-model parameters": deg[deg["group"].isin(["baseline", "parameter"])], "all degassing-type scenarios": deg} for sname, d in sets.items(): for name in list(MOBILE) + ["mobile-lid formulations together"]: g = d[d["formulation"].isin(MOBILE)] if name.startswith("mobile") else d[d["formulation"] == name] if len(g) == 0: continue row = {"scenario_set": sname, "formulation": name, "n_scenarios": int(g["scenario"].nunique())} for n_ in (2, 4, 6): c = f"rank_n{n_}" row[f"{c}_min"], row[f"{c}_max"] = float(g[c].min()), float(g[c].max()) row[f"{c}_min_scenario"] = g.loc[g[c].idxmin(), "scenario"] + " | " + g.loc[g[c].idxmin(), "formulation"] row[f"{c}_max_scenario"] = g.loc[g[c].idxmax(), "scenario"] + " | " + g.loc[g[c].idxmax(), "formulation"] rng_rows.append(row) mel = out[(out["criterion"] == "any_melting") & (out["group"] != "other")] for name in MOBILE: g = mel[mel["formulation"] == name] row = {"scenario_set": "melting criterion for every planet, all population scenarios", "formulation": name, "n_scenarios": int(g["scenario"].nunique())} for n_ in (2, 4, 6): c = f"rank_n{n_}" row[f"{c}_min"], row[f"{c}_max"] = float(g[c].min()), float(g[c].max()) row[f"{c}_min_scenario"] = g.loc[g[c].idxmin(), "scenario"] + " | " + name row[f"{c}_max_scenario"] = g.loc[g[c].idxmax(), "scenario"] + " | " + name rng_rows.append(row) rng = pd.DataFrame(rng_rows) rng.to_csv(os.path.join(DATA_DIR, "rank_ranges.csv"), index=False) return out, rng def hyper_stage(emp, n_draws=40, n=600, seed=SEED + 3): """Propagate hyper-parameter uncertainty (abundance zero point, raw trend and scatter, decay-model strength, activation energy, solidus, continental fraction, Earth's heat loss, viscosity scatter, thermostat). NOT varied here: mass range, degassing threshold and law, criterion mix; those are the scenarios of scenario_stage and they move the rank more than these draws do.""" rng = np.random.default_rng(seed) z = pop.draw_normals(n, seed + 100) k_sun = int(np.argmin(np.abs(pop.NODE_AGES - AGE_SUN))) rows = [] for d in range(n_draws): h = pop.draw_hyper(rng, base_hyper(emp), emp) for name in thermal.FORMULATIONS: r = pop.run_population(name, h, z) for crit, key in CRITERIA.items(): if name == "F3_stagnant" and crit == "platespeed>=10%": continue L = r[key] sv = pop.survival(L[k_sun]) rk = ranks_for(L, n_values=(2, 4, 6)) rows.append({"draw": d, "formulation": name, "criterion": crit, "a": h["a"], "b": h["b"], "s": h["s"], "ism_scale": h["ism_scale"], "T_cap": h["T_cap"], "cap_raised": bool(getattr(r["model"], "cap_raised", False)), "E_kJ": h["E"] / 1e3, "t_sol0": h["t_sol0"], "f_cc": h["f_cc"], "Q0": h["Q0"], "s_eta": h["s_eta"], "median_lifetime_solar_epoch": float(np.median(L[k_sun])), **{f"P_hab_{int(a)}Gyr": v for a, v in zip(pop.TARGET_AGES, sv)}, **{f"rank_n{n_}": v for n_, v in rk.items()}}) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA_DIR, "hyper_draws.csv"), index=False) cols = ["median_lifetime_solar_epoch"] + [f"P_hab_{int(a)}Gyr" for a in pop.TARGET_AGES] + [f"rank_n{n_}" for n_ in (2, 4, 6)] srows = [] for (name, crit), g in out.groupby(["formulation", "criterion"], sort=False): row = {"formulation": name, "criterion": crit, "n_draws": len(g)} for c in cols: q = np.percentile(g[c], [2.5, 16, 50, 84, 97.5]) row[f"{c}_p2.5"], row[f"{c}_p16"], row[f"{c}_p50"], row[f"{c}_p84"], row[f"{c}_p97.5"] = q srows.append(row) summ = pd.DataFrame(srows) summ.to_csv(os.path.join(DATA_DIR, "hyper_summary.csv"), index=False) return out, summ def conditional_shares(hyper): """Share of hyper-parameter draws with rank >= 10% (and >= 1%), separately for every formulation and criterion. These are CONDITIONAL shares: conditional on the formulation, the criterion, the 0.5-2 Earth-mass prior and the 10% threshold. They are not probabilities that geophysics resolves the earliness puzzle, and no average over cells is formed (each cell is near 0 or 1, so an average would only reflect the weights).""" rows = [] for (name, crit), g in hyper.groupby(["formulation", "criterion"], sort=False): row = {"formulation": name, "criterion": crit, "n_draws": len(g)} for n_ in (2, 4, 6): row[f"share_rank_n{n_}_ge_10pct"] = float((g[f"rank_n{n_}"] >= 0.10).mean()) row[f"share_rank_n{n_}_ge_1pct"] = float((g[f"rank_n{n_}"] >= 0.01).mean()) rows.append(row) out = pd.DataFrame(rows) out.to_csv(os.path.join(DATA_DIR, "conditional_shares.csv"), index=False) return out