"""Figures for the mdwarf_early study (static PNG, matplotlib).""" import os import numpy as np import earliness as E import plotstyle as ps from common import FIG_DIR, V_ESC_EARTH, M_DWARF_MAX_MASS from curated_outcomes import SOLAR_SYSTEM plt = ps.apply() ORD4 = ["#86b6ef", "#3987e5", "#1c5cab", "#0d366b"] # validated ordinal ramp (4 steps) ORD5 = ["#86b6ef", "#5598e7", "#2a78d6", "#1c5cab", "#104281"] # validated ordinal ramp (5 steps) RING = dict(edgecolors=ps.SURFACE, linewidths=1.2) def _save(fig, name): fig.savefig(os.path.join(FIG_DIR, name)) plt.close(fig) # ------------------------------------------------------------------------------------------- def fig_shoreline(census, targets, samples): """Figure 1: the census in the (escape velocity, instellation) plane with the fitted shoreline.""" fig, ax = plt.subplots(figsize=(8.4, 6.2)) c = census cls = {"bare_rock_favoured": 0, "thick_atm_ruled_out": 0, "inconclusive": 1, "anomalous": 1, "unpublished": 1, "tentative_atm": 2, "atm_detected": 3} col = {0: ps.BLUE, 1: ps.NEUTRAL, 2: ps.ORANGE, 3: ps.ORANGE} # shoreline lines for three host luminosities (posterior median and 68% band, base fit) v = np.logspace(np.log10(1.0), np.log10(30.0), 80) for logL, lab in ((0.0, "Sun-like host"), (-2.0, "0.01 Lsun host (M3-M4)"), (-3.26, "TRAPPIST-1-like host")): y = samples[:, 0][:, None] + samples[:, 1][:, None] * np.log10(v / V_ESC_EARTH)[None, :] + samples[:, 2][:, None] * logL lo, med, hi = np.percentile(y, [16, 50, 84], axis=0) ax.fill_between(v, 10 ** lo, 10 ** hi, color=ps.INK2, alpha=0.06, lw=0) ax.plot(v, 10 ** med, color=ps.INK2, lw=1.6) ax.annotate(lab, (v[-1], 10 ** med[-1]), xytext=(6, 0), textcoords="offset points", fontsize=8.5, color=ps.INK2, ha="left", va="center", annotation_clip=False) # Solar System for b in SOLAR_SYSTEM: has = b["pi"] > 0.5 ax.scatter(b["v_esc"], b["S"], marker="s", s=70, color=ps.ORANGE if has else ps.BLUE, zorder=4, **RING) if b["name"] in ("Earth", "Venus", "Mars", "Mercury", "Titan", "Moon", "Pluto"): right = b["name"] in ("Earth", "Pluto") ax.annotate(b["name"], (b["v_esc"], b["S"]), xytext=(7, 0) if right else (-7, -3), textcoords="offset points", fontsize=8, color=ps.INK2, ha="left" if right else "right", va="center") # exoplanets label_set = {"TRAPPIST-1 b": "T-1 b", "TRAPPIST-1 c": "c", "TRAPPIST-1 d": "T-1 d", "TRAPPIST-1 e": "e", "TRAPPIST-1 f": "f", "TRAPPIST-1 g": "g", "TRAPPIST-1 h": "T-1 h", "LHS 1140 b": "LHS 1140 b", "LTT 1445 A b": "LTT 1445 A b", "55 Cnc e": "55 Cnc e", "K2-18 b": "K2-18 b", "L 98-59 b": "L 98-59 b", "TOI-561 b": "TOI-561 b", "GJ 367 b": "GJ 367 b"} for _, r in c.iterrows(): k = cls[r.outcome] if r.group == "sub_neptune": mk, size = "^", 70 elif r.host_class == "M": mk, size = "o", 62 else: mk, size = "D", 52 if k == 2: # tentative: hollow ax.scatter(r.v_esc_kms, r.S_earth, marker=mk, s=size - 12, facecolors=ps.SURFACE, edgecolors=col[k], linewidths=1.8, zorder=5) elif k == 1: ax.scatter(r.v_esc_kms, r.S_earth, marker=mk, s=size - 12, facecolors=ps.SURFACE, edgecolors=col[k], linewidths=1.5, zorder=4) else: ax.scatter(r.v_esc_kms, r.S_earth, marker=mk, s=size, color=col[k], zorder=5, **RING) if r.planet in label_set: dx, dy, ha = 6, 3, "left" if r.planet in ("TRAPPIST-1 c", "TRAPPIST-1 e", "TRAPPIST-1 f", "TRAPPIST-1 g"): dx, dy = 6, -1 if r.planet in ("TRAPPIST-1 b", "TRAPPIST-1 d", "TRAPPIST-1 h", "L 98-59 b", "GJ 367 b", "TOI-561 b"): dx, ha = -6, "right" if r.planet == "TRAPPIST-1 b": dx, dy, ha = -6, -6, "right" if r.planet == "TRAPPIST-1 e": dx, dy, ha = -6, -5, "right" if r.planet == "LTT 1445 A b": dx, dy = 5, 9 if r.planet == "55 Cnc e": dx, dy = 7, 4 ax.annotate(label_set[r.planet], (r.v_esc_kms, r.S_earth), xytext=(dx, dy), textcoords="offset points", fontsize=8, color=ps.INK2, ha=ha, va="center") ax.axhspan(0.25, 0.9, color=ps.AQUA, alpha=0.08, lw=0, zorder=0) ax.text(1.06, 0.47, "habitable-zone instellation", fontsize=8.5, color=ps.INK2, va="center") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlim(1.0, 30); ax.set_ylim(3e-4, 2e4) ax.set_xticks([1, 2, 3, 5, 10, 20, 30]); ax.set_xticklabels(["1", "2", "3", "5", "10", "20", "30"]) ax.set_xlabel("escape velocity (km/s)") ax.set_ylabel("bolometric instellation (Earth = 1)") ax.set_title("Atmosphere verdicts through Oct 2026 and the fitted cosmic shoreline") from matplotlib.lines import Line2D h = [Line2D([], [], marker="o", ls="", color=ps.BLUE, markersize=8, markeredgecolor=ps.SURFACE, label="no thick atmosphere / airless"), Line2D([], [], marker="o", ls="", markerfacecolor=ps.SURFACE, markeredgecolor=ps.NEUTRAL, markeredgewidth=1.5, markersize=7, label="inconclusive or unpublished"), Line2D([], [], marker="o", ls="", markerfacecolor=ps.SURFACE, markeredgecolor=ps.ORANGE, markeredgewidth=1.8, markersize=7, label="tentative atmosphere"), Line2D([], [], marker="o", ls="", color=ps.ORANGE, markersize=8, markeredgecolor=ps.SURFACE, label="atmosphere"), Line2D([], [], marker="o", ls="", color=ps.INK2, markersize=6, label="circle: M-dwarf planet"), Line2D([], [], marker="D", ls="", color=ps.INK2, markersize=5.5, label="diamond: G/K-dwarf planet"), Line2D([], [], marker="^", ls="", color=ps.INK2, markersize=6.5, label="triangle: sub-Neptune"), Line2D([], [], marker="s", ls="", color=ps.INK2, markersize=6, label="square: Solar System"), Line2D([], [], color=ps.INK2, lw=1.6, label="shoreline, median and 68% band")] ax.legend(handles=h, loc="upper center", bbox_to_anchor=(0.5, -0.11), ncol=3, fontsize=8.5, handletextpad=0.4, labelcolor=ps.INK2, borderaxespad=0, columnspacing=1.5) _save(fig, "fig1_census_shoreline.png") # ------------------------------------------------------------------------------------------- def fig_predictions(pred, fm): """Figure 2: predicted retention probabilities for temperate planets, and f_atm versus stellar mass.""" fig, (ax, ax2) = plt.subplots(1, 2, figsize=(11.5, 6.6), gridspec_kw=dict(width_ratios=[1.25, 1], wspace=0.34)) keep = ["TRAPPIST-1 d", "TRAPPIST-1 e", "TRAPPIST-1 f", "TRAPPIST-1 g", "Teegarden's Star b", "Proxima Cen b", "GJ 1002 b", "Ross 128 b", "Gliese 12 b", "Wolf 1069 b", "GJ 1061 d", "LP 890-9 c", "TOI-700 d", "LHS 1140 b", "TOI-715 b", "LTT 1445 A c", "TOI-771 b", "TOI-198 b", "HD 260655 c"] p = pred.set_index("planet").loc[keep].reset_index() p = p.sort_values("pooled_mean").reset_index(drop=True) y = np.arange(len(p)) # RV-only planets: m sin i is a lower limit on the mass, so the plotted probability is a lower limit too names = [f"{n} *" if lim else n for n, lim in zip(p.planet, p.mass_is_lower_limit)] specs = [("base", ps.BLUE, "o", "rocky planets + Solar System", 0.24), ("plus_subneptunes", ps.ORANGE, "s", "+ sub-Neptune anchors", 0.0), ("flat_slope_priors", ps.AQUA, "D", "flat priors on slopes", -0.24)] for key, colr, mk, lab, off in specs: ax.hlines(y + off, p[f"{key}_p05"], p[f"{key}_p95"], color=colr, lw=1.6, alpha=0.55) ax.scatter(p[f"{key}_mean"], y + off, marker=mk, s=34, color=colr, zorder=3, label=lab, **RING) # published atmosphere-interior evolution model for TRAPPIST-1 e and f (Krissansen-Totton 2023, # arXiv:2306.05397: "e and f retain atmospheres in ~98% of model runs") - a model, not a measurement kt = [i for i, n in enumerate(p.planet) if n in ("TRAPPIST-1 e", "TRAPPIST-1 f")] ax.scatter([0.98] * len(kt), y[kt], marker="X", s=70, color=ps.INK2, zorder=4, label="evolution model (K-T 2023)", **RING) ax.set_yticks(y); ax.set_yticklabels(names, fontsize=9) ax.set_xlim(-0.02, 1.02); ax.set_ylim(-0.7, len(p) - 0.3) ax.set_xlabel("P(thick atmosphere): posterior mean and 5-95% range\n* mass is m sin i (lower limit)") ax.set_title("Shoreline extrapolations: structural, not measured") ax.grid(axis="y", visible=False) ax.legend(loc="upper center", fontsize=8.5, labelcolor=ps.INK2, bbox_to_anchor=(0.45, -0.115), ncol=2, columnspacing=1.5, handletextpad=0.4) m = fm.stellar_mass_Msun ax2.fill_between(m, fm.pooled_p05, fm.pooled_p95, color=ps.NEUTRAL, alpha=0.16, lw=0, label="reference mixture, 5-95%") for key, colr, mk, lab, off in specs: ax2.plot(m, fm[f"{key}_mean"], color=colr, lw=2, label=lab) ax2.plot(m, fm.pooled_mean, color=ps.INK, lw=2.2, label="reference mixture, mean") ax2.plot(m, fm.solar_system_only_mean, color=ps.NEUTRAL, lw=1.6, label="Solar System only (pre-JWST)") ax2.text(m.iloc[-1] * 1.01, fm.solar_system_only_mean.iloc[-1], " pre-JWST", fontsize=8, color=ps.INK2, va="center") ax2.set_xscale("log"); ax2.set_xlim(0.08, 0.72); ax2.set_ylim(0, 1.0) ax2.set_xticks([0.08, 0.1, 0.15, 0.2, 0.3, 0.4, 0.6]); ax2.set_xticklabels(["0.08", "0.1", "0.15", "0.2", "0.3", "0.4", "0.6"]) ax2.set_xlabel("host-star mass (Msun)") ax2.set_ylabel("mean retention probability, HZ rocky planets") ax2.set_title("Extrapolated retention by host mass") ax2.legend(loc="upper left", fontsize=8.5, labelcolor=ps.INK2) _save(fig, "fig2_shoreline_predictions.png") # ------------------------------------------------------------------------------------------- def fig_hanson(n_h, L_grid, hanson_grid, rank_all, rank_noM): """Figure 3: rank of today's date versus hard-steps power n and maximum habitable lifetime. All three grids are evaluated on the same (n_h, L_grid) mesh.""" fig, axes = plt.subplots(1, 3, figsize=(13.2, 4.6), sharey=True, gridspec_kw=dict(wspace=0.08)) levels = [-20, -10, -5, -3, -2, -1, np.log10(0.25)] fmt = {-20: "1e-20", -10: "1e-10", -5: "1e-5", -3: "0.1%", -2: "1%", -1: "10%", np.log10(0.25): "25%"} cols = ["#0d366b", "#104281", "#1c5cab", "#2a78d6", "#3987e5", "#5598e7", "#86b6ef"] panels = [(hanson_grid, "Hanson et al. (2021) eq. 1, replicated"), (rank_all, "This work: all stars habitable (p_M = 1)"), (rank_noM, "This work: hosts above 0.6 Msun only")] for k, (ax, (grid, title)) in enumerate(zip(axes, panels)): z = np.log10(np.maximum(grid, 1e-300)) cs = ax.contour(L_grid, n_h, z, levels=levels, colors=cols, linewidths=1.8) if k < 2: ax.clabel(cs, fmt=fmt, fontsize=8.5, colors=ps.INK, inline=True, inline_spacing=2) else: # contours are horizontal beyond ~45 Gyr (longest-lived 0.6 Msun host): label them there jx = int(np.argmin(abs(np.log(L_grid) - np.log(600.0)))) for lev in levels: col_z = z[:, jx] if col_z.min() < lev < col_z.max(): yv = float(np.exp(np.interp(-lev, -col_z, np.log(n_h)))) ax.text(L_grid[jx], yv * 1.04, fmt[lev], fontsize=8.5, color=ps.INK, ha="center", va="bottom") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xticks([5, 10, 20, 50, 100, 500, 1000, 5000, 20000]) ax.set_xticklabels(["5", "10", "20", "50", "100", "500", "1k", "5k", "20k"], fontsize=8.5) ax.set_yticks([1, 2, 3, 5, 10, 20]); ax.set_yticklabels(["1", "2", "3", "5", "10", "20"]) ax.set_ylim(1, 20); ax.set_xlim(5, 2e4) ax.set_xlabel("maximum habitable planet lifetime (Gyr)") ax.set_title(title, fontsize=10.5) axes[0].set_ylabel("hard-steps power n") fig.suptitle("Percentile rank of today's date among civilisation birth dates (contours)", x=0.125, ha="left", fontsize=12, fontweight="bold", y=1.02) _save(fig, "fig3_rank_vs_n_lifetime.png") # ------------------------------------------------------------------------------------------- def fig_rank_vs_pm(pm_01): """Figure 4: rank versus the M-dwarf habitability weight p_M for several n.""" fig, axes = plt.subplots(1, 2, figsize=(11.5, 5.2), sharey=True, gridspec_kw=dict(wspace=0.06)) pgrid = np.logspace(-13, 0, 200) hw = E.host_weight(M_DWARF_MAX_MASS) med, up = np.percentile(pm_01, [50, 95]) for k, (ax, (L, title)) in enumerate(zip(axes, ((1e9, "Full stellar lifetimes (no planetary cap)"), (50.0, "Habitable lifetime capped at 50 Gyr")))): for yv in (0.10, 0.01): ax.axhline(yv, color=ps.AXIS, lw=1.0) ax.axvline(up, color=ps.INK2, lw=1.2) ax.axvline(med, color=ps.INK2, lw=1.2) if k == 0: ax.text(up * 0.6, 3e-10, "95th pct", fontsize=8.5, color=ps.INK2, ha="right", va="bottom") ax.text(med * 1.5, 3e-12, "median", fontsize=8.5, color=ps.INK2, ha="left", va="bottom") ax.text(med * 0.6, 1.5e-9, "p_M of 0.1 Msun hosts,\nshoreline x assumed penalties\n(median, 95th percentile)", fontsize=8.5, color=ps.INK2, ha="right", va="bottom") for n, colr in zip((1, 2, 3, 6), ORD4): num, den = E.birth_integrals(n, L) A, B, C, D = (num * hw).sum(), (num * (1 - hw)).sum(), (den * hw).sum(), (den * (1 - hw)).sum() r = (A + pgrid * B) / (C + pgrid * D) ax.plot(pgrid, r, color=colr, lw=2.2, label=f"n = {n}") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlim(1e-13, 1.0); ax.set_ylim(1e-13, 1.5) ax.set_xticks([1e-12, 1e-9, 1e-6, 1e-3, 1]) ax.set_title(title, fontsize=10.5) for yv, lab in ((0.10, "rank 10%"), (0.01, "rank 1%")): axes[1].annotate(lab, (1.0, yv), xytext=(6, 0), textcoords="offset points", fontsize=8.5, color=ps.INK2, va="center", ha="left", annotation_clip=False) axes[0].set_ylabel("humanity's percentile rank (fraction born earlier)") axes[1].legend(loc="lower left", title="hard steps", fontsize=9, title_fontsize=9, labelcolor=ps.INK2, bbox_to_anchor=(0.03, 0.05)) fig.supxlabel("p_M: habitability weight of M-dwarf planets relative to FGK planets", fontsize=10.5, color=ps.INK2, y=0.0) fig.suptitle("How airless must M dwarfs be? Rank of today's date versus the M-dwarf weight", x=0.125, ha="left", fontsize=12, fontweight="bold", y=0.99) _save(fig, "fig4_rank_vs_pM.png") # ------------------------------------------------------------------------------------------- def fig_early_by_n(cw, byn): """Figure 5. Left: rank of today's date when M dwarfs are excluded, versus the number of hard steps, for four host-mass cuts and two habitable-window fractions. Right: prior-weighted fraction of the Lmax range (log-uniform 5 Gyr - 20,000 Gyr) in which the rank is below 1%, by n, for four M-dwarf scenarios. Neither panel is a posterior probability.""" fig, (ax, ax2) = plt.subplots(1, 2, figsize=(13.2, 5.2), gridspec_kw=dict(wspace=0.62)) n = np.arange(1, 11) cuts = [0.6, 0.65, 0.7, 0.8] for yv in (0.10, 0.01): ax.axhline(yv, color=ps.AXIS, lw=1.0, zorder=1) for cut, colr in zip(cuts, ORD4[::-1]): for fwin, ls, lw in ((1.0, "-", 2.2), (0.6, (0, (4, 2)), 1.4)): r = cw[(cw.f_window == fwin) & (np.isclose(cw.host_mass_cut_Msun, cut))].iloc[0] v = np.array([r[f"rank_n{k}"] for k in n]) ax.plot(n, v, color=colr, lw=lw, ls=ls, marker="o" if fwin == 1.0 else None, ms=4.5, markeredgecolor=ps.SURFACE, markeredgewidth=0.8, zorder=3) if fwin == 1.0: ax.annotate(f"> {cut:g} Msun\n(to {r.habitable_lifetime_at_cut_Gyr:.0f} Gyr)", (10, v[-1]), xytext=(7, 0), textcoords="offset points", fontsize=8.5, color=ps.INK2, va="center", ha="left", annotation_clip=False) for yv, lab in ((0.10, "10%"), (0.01, "1%")): ax.text(1.05, yv * 1.12, lab, fontsize=8.5, color=ps.INK2, va="bottom", ha="left") ax.set_yscale("log"); ax.set_ylim(1e-6, 1.0); ax.set_xlim(0.8, 10.2) ax.set_xticks(n) ax.set_xlabel("number of hard steps n") ax.set_ylabel("rank of today's date (fraction born earlier)") ax.set_title("No M dwarfs: the rank depends on where the cut is", fontsize=10.5) from matplotlib.lines import Line2D ax.legend(handles=[Line2D([], [], color=ps.INK2, lw=2.2, marker="o", ms=4.5, label="whole main-sequence lifetime habitable"), Line2D([], [], color=ps.INK2, lw=1.4, ls=(0, (4, 2)), label="first 60% of it habitable")], loc="lower left", fontsize=8.5, labelcolor=ps.INK2, title="habitable hosts above the cut (longest lifetime)", title_fontsize=8.5) b = byn.set_index(["scenario", "n"]).frac_rank_below_1pct series = [("p_M = 1 (all stars equally habitable)", "all stars equally habitable", "o"), ("JWST shoreline (reference mixture), retention only", "shoreline retention only", "s"), ("JWST shoreline (reference mixture) x penalties", "retention x assumed penalties", "D"), ("M dwarfs excluded: p_M = 0 (hosts > 0.6 Msun)", "no M dwarfs (> 0.6 Msun)", "^")] fitn = [f"JWST shoreline fit '{k}', retention only" for k in ("base", "plus_subneptunes", "flat_slope_priors")] lo = np.min([[b[(f, k)] for k in n] for f in fitn], axis=0) hi = np.max([[b[(f, k)] for k in n] for f in fitn], axis=0) ax2.fill_between(n, lo, hi, color=ORD4[1], alpha=0.18, lw=0, zorder=1, label="range over three shoreline fits") ends = [] for (key, lab, mk), colr in zip(series, ORD4): v = np.array([b[(key, k)] for k in n]) ax2.plot(n, v, color=colr, lw=2.0, marker=mk, ms=5.5, markeredgecolor=ps.SURFACE, markeredgewidth=0.8, zorder=3, label=lab) ends.append(v[0]) for (key, lab, mk), v0 in zip(series, ends): if v0 > 0.005: ax2.annotate(f"{v0:.2f}", (1, v0), xytext=(-8, 0), textcoords="offset points", fontsize=8.5, color=ps.INK2, ha="right", va="center") ax2.set_xlim(0.15, 10.2); ax2.set_ylim(-0.03, 1.0) ax2.set_xticks(n) ax2.set_xlabel("number of hard steps n") ax2.set_ylabel("fraction of the Lmax prior with rank < 1%") ax2.set_title("M-dwarf weights change the answer mainly for n = 1-3", fontsize=10.5) ax2.legend(loc="lower right", fontsize=8.5, labelcolor=ps.INK2) fig.suptitle("Is humanity early without M dwarfs, and where can M-dwarf evidence change the answer?", x=0.07, ha="left", fontsize=12, fontweight="bold", y=1.0) _save(fig, "fig5_early_by_n.png")