"""Step 11: figures (figures/*.png) and the slim per-star table (data/stars_slim.csv.gz).""" import os import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd from common import RAW, DATA, FIG, T_GRID import importlib ex = importlib.import_module("03_excess") inj = importlib.import_module("04_injection") im = importlib.import_module("07_images") lm = importlib.import_module("09_limits") DERIVED = os.path.join(RAW, "derived") TCOL = {100: "#313695", 150: "#4575b4", 200: "#74add1", 300: "#f46d43", 400: "#d73027", 600: "#a50026"} # Suazo et al. 2022 (MNRAS 512, 2988), 100 pc sample, 300 K: fraction of stars compatible with # a Dyson sphere of covering fraction gamma (their Section 4 / abstract) SUAZO22_100PC_300K = {0.1: 6.6e-3, 0.5: 1.9e-4, 0.9: 1.8e-5} def fig_significance(d): fig, axs = plt.subplots(1, 2, figsize=(10, 4.2), sharey=True) bins = np.arange(0, 60.01, 0.5) for ax, b in zip(axs, ("3", "4")): q = d["q" + b].values pos = d.loc[q, "chi" + b].values neg = -d.loc[q, f"chi{b}_neg"].values ax.hist(np.clip(pos, 0, 59.9), bins=bins, histtype="step", color="#d73027", lw=1.4, label=r"excess side: $\chi=\min(\chi_{Ks},\chi_{W1})$") ax.hist(np.clip(neg, 0, 59.9), bins=bins, histtype="step", color="#4575b4", lw=1.4, label=r"deficit side (mirror): $-\max(\chi_{Ks},\chi_{W1})$") x = np.linspace(0, 8, 200) sig = 0.9 ax.plot(x, q.sum() * 0.5 * np.exp(-0.5 * (x / sig) ** 2) / (sig * np.sqrt(2 * np.pi)), "k:", lw=1, label=r"Gaussian, $\sigma=0.9$") ax.axvline(ex.THR, color="gray", ls="--", lw=1) ax.set_yscale("log") ax.set_ylim(0.5, 3e5) ax.set_xlim(0, 60) ax.set_xlabel(f"W{b} significance (last bin: $\\geq$ 60)") n_pos, n_neg = int((pos >= ex.THR).sum()), int((neg >= ex.THR).sum()) ax.set_title(f"W{b}: {q.sum():,} searchable stars; $\\geq${ex.THR:g}$\\sigma$: {n_pos} excess, {n_neg} deficit", fontsize=10) axs[0].set_ylabel("stars per 0.5 bin") axs[0].legend(fontsize=8, loc="upper right") fig.tight_layout() fig.savefig(os.path.join(FIG, "fig1_significance.png"), dpi=150) plt.close(fig) def fig_completeness(d): z = np.load(os.path.join(DERIVED, "injection.npz")) rec, f_fine = z["rec"], z["f_fine"] fig, axs = plt.subplots(1, 4, figsize=(15.5, 4), sharey=True) for ax, cls in zip(axs, ("A", "FGK", "M")): sel = (d.cls == cls).values for it, T in enumerate(T_GRID): ax.plot(f_fine, rec[sel, it, :].mean(axis=0), color=TCOL[T], lw=1.6, label=f"{T} K") srch = ((d.q3 | d.q4).values & sel).sum() / sel.sum() ax.axhline(srch, color="gray", ls=":", lw=1) ax.text(1.2e-5, srch + 0.012, "searchable fraction", fontsize=8, color="gray") ax.set_title(f"{cls} stars by $M_G$ (N = {sel.sum():,})", fontsize=10) axs[0].legend(title="swarm temperature", fontsize=8) # M dwarfs split by absolute magnitude, 300 K ax = axs[3] it = T_GRID.index(300) sub_col = ["#081d58", "#225ea8", "#41b6c4", "#fd8d3c", "#bd0026"] for name, col in zip(lm.M_SUB, sub_col): sel = lm.stratum_mask(d, name) ax.plot(f_fine, rec[sel, it, :].mean(axis=0), color=col, lw=1.6, label=f"{name} ({sel.sum():,})") ax.set_title("M dwarfs by sub-range, 300 K", fontsize=10) ax.legend(title="Mamajek $M_G$ range (N)", fontsize=8, loc="upper left") for ax in axs: ax.set_xscale("log") ax.set_xlim(1e-5, 0.9) ax.set_ylim(0, 0.66) ax.set_xlabel("covering fraction f") ax.grid(alpha=0.25) axs[0].set_ylabel("completeness C(f, T): recovered / all parent stars") fig.tight_layout() fig.savefig(os.path.join(FIG, "fig2_completeness.png"), dpi=150) plt.close(fig) def fig_limits(): lim = pd.read_csv(os.path.join(DATA, "limits.csv")) rng = pd.read_csv(os.path.join(DATA, "limit_ranges.csv")) fig, axs = plt.subplots(2, 2, figsize=(11.5, 8.6), sharex=True, sharey=True) titles = {"FGK": "FGK stars", "M": "M dwarfs (in practice M0-M4.5 for f < 0.05)", "all": "all classes pooled (valid only for a type-independent rate)"} for ax, cls in zip(axs.flat[:3], ("FGK", "M", "all")): s = lim[(lim.cls == cls) & (lim.T_K == 300)] r = rng[rng.cls == cls] ax.fill_between(r.f, r.UL95_all_min, r.UL95_all_max, color="#d73027", alpha=0.18, lw=0, label="every excess: range over 100-600 K, tolerance 0.2-0.5") ax.fill_between(r.f, r.UL95_unexplained_min, r.UL95_unexplained_max, color="#1a9850", alpha=0.18, lw=0, label="identifications accepted: same range") ax.plot(s.f, s.UL95_all, "-s", color="#d73027", ms=4, lw=2.0, label="every excess counted (label-free), 300 K") ax.plot(s.f, s.UL95_strict, "-.", color="#ff7f00", lw=1.4, label="unexplained + catalogue-only identifications, 300 K") ax.plot(s.f, s.UL95_unexplained, "--o", color="#1a9850", ms=3.5, lw=1.4, label="unexplained only (every identification accepted), 300 K") ax.plot(s.f, s.UL95_sensitivity, "k:", lw=1.3, label="sensitivity floor (zero candidates)") if cls == "all": g = np.array(list(SUAZO22_100PC_300K)) ax.plot(g, [SUAZO22_100PC_300K[k] for k in g], "*", color="purple", ms=13, label="Suazo et al. 2022, 100 pc, 300 K:\nfraction of sources compatible (not a 95% limit)") ax.set_title(titles[cls], fontsize=10) ax.set_xscale("log") ax.set_yscale("log") ax.grid(alpha=0.25, which="both") ax.legend(fontsize=6.8, loc="upper right") ax = axs[1, 1] sub_col = ["#081d58", "#225ea8", "#41b6c4", "#fd8d3c", "#bd0026"] for name, col in zip(lm.M_SUB, sub_col): s = lim[(lim.cls == name) & (lim.T_K == 300)] ok = s.completeness >= 0.01 ax.plot(s.f[ok], s.UL95_all[ok], "-", color=col, lw=1.6, label=f"{name} (N = {int(s.N_parent.iloc[0]):,})") ax.plot(s.f[~ok], s.UL95_all[~ok], ":", color=col, lw=1.0) ax.set_title("M sub-ranges, every excess counted, 300 K\n(dotted: completeness below 1%)", fontsize=10) ax.grid(alpha=0.25, which="both") ax.legend(fontsize=8, title="Mamajek $M_G$ range") ax.set_ylim(1e-5, 1) ax.set_xlim(8e-5, 1) for a in axs[1]: a.set_xlabel("covering fraction f") for a in axs[:, 0]: a.set_ylabel("95% upper limit on fraction of stars hosting the swarm") fig.tight_layout() fig.savefig(os.path.join(FIG, "fig3_limits.png"), dpi=150) plt.close(fig) def fig_accounting(): c = pd.read_csv(os.path.join(DATA, "candidates.csv")) null = pd.read_csv(os.path.join(DATA, "image_centroid_null.csv")) acc = pd.read_csv(os.path.join(DATA, "accounting.csv")).set_index("cls") fig, axs = plt.subplots(1, 3, figsize=(14, 4.4), gridspec_kw=dict(width_ratios=[1.1, 1, 1])) order = ["young_star", "catalogued_excess_only", "other_known_type", "comoving_companion", "background_catalogued", "background_image_offset", "background_visual", "unexplained"] cols = ["#1a9850", "#fdae61", "#66bd63", "#a6d96a", "#4575b4", "#74add1", "#abd9e9", "#d73027"] names = {"young_star": "independent youth evidence", "catalogued_excess_only": "published excess catalogue only"} ax = axs[0] bottom = np.zeros(3) xs = np.arange(3) for e, col in zip(order, cols): v = np.array([((c.cls == k) & (c.explanation == e)).sum() for k in ("A", "FGK", "M")]) ax.bar(xs, v, bottom=bottom, color=col, label=names.get(e, e.replace("_", " ")), width=0.6) bottom += v for i, k in enumerate(("A", "FGK", "M")): # range of the chance-alignment expectation: direct count with a 4 arcsec cap up to the # analytic trigger-area estimate (caps 6 / 8 arcsec) hi = acc.loc[k, "expected_chance_alignments"] lo = acc.loc[k, "direct_count_cap4"] if "direct_count_cap4" in acc else hi mid = acc.loc[k, "direct_count_cap6"] if "direct_count_cap6" in acc else hi ax.plot([i + 0.38, i + 0.38], [lo, hi], "-", color="k", lw=1.5) ax.plot([i + 0.38], [mid], "D", color="k", ms=5) ax.plot([], [], "-kD", ms=5, label="chance alignments expected from control\nfields (range over blending radius)") ax.set_xticks(xs) ax.set_xticklabels([f"{k}\n({int(acc.loc[k, 'n_searchable']):,} searchable)" for k in ("A", "FGK", "M")]) ax.set_ylabel(r"stars with a $\geq5\sigma$ W3 or W4 excess") ax.legend(fontsize=7.5, loc="upper left") ax.set_title(f"How the {len(c)} excess stars are identified", fontsize=10) for ax, b, k in ((axs[1], "W3", "det3"), (axs[2], "W4", "det4")): nb = null[null.band == b]["n" + b].dropna().values groups = [("non-excess comparison stars", nb, "k", "-"), ("young / catalogued excess", c[c[k] & c.explanation.isin(["young_star", "catalogued_excess_only"])] ["n" + b].dropna().values, "#1a9850", "-"), ("catalogued background", c[c[k] & (c.explanation == "background_catalogued")]["n" + b].dropna().values, "#4575b4", "-"), ("no catalogue explanation", c[c[k] & c.explanation.isin(["background_image_offset", "background_visual", "unexplained", "comoving_companion"])] ["n" + b].dropna().values, "#d73027", "-")] for lab, v, col, ls in groups: if len(v) < 3: continue v = np.sort(v) ax.step(v, np.arange(1, len(v) + 1) / len(v), where="post", color=col, ls=ls, lw=1.6, label=f"{lab} (n={len(v)})") import json thr = json.load(open(os.path.join(DATA, "image_centroid_calib.json")))[b]["threshold"] ax.axvline(thr, color="gray", ls="--", lw=1) ax.set_xscale("log") ax.set_xlim(0.1, 40) ax.set_xlabel(f"{b} centroid offset from the star / expected noise") ax.set_ylabel("cumulative fraction") ax.set_title(f"{b} excess stars: is the light on the star?", fontsize=10) ax.legend(fontsize=7.5, loc="lower right") ax.grid(alpha=0.25) fig.tight_layout() fig.savefig(os.path.join(FIG, "fig4_accounting.png"), dpi=150) plt.close(fig) def fig_cutouts(d): c = pd.read_csv(os.path.join(DATA, "candidates.csv")) cen = pd.read_csv(os.path.join(DATA, "image_centroids.csv")) examples = {3478940625208241920: "TWA 33: known young disc", 5072573534418641536: "LP 886-52: catalogued galaxy 4.3\" away", 5095258452165352704: "W3/W4 source 15\" south (image offset)", 5260623420435383168: "SCR J0615-7615: galaxy at WISE epoch", 4936999799344230400: "unexplained, on the star", 6869593898383897600: "V4333 Sgr: unexplained W4 excess"} order = [s for s in examples if s in set(c.sid)] if not order: return dets = d[d.det].merge(c[["sid", "explanation"]], on="sid", how="left") labels = {s: examples[s][:38] for s in order} im.sheets(dets, cen, order, "fig5_cutouts", per=len(order), outdir=FIG, labels=labels, dpi=140) def slim_table(d): """Per-star table of the quantities the limits rest on (searchable stars only).""" s = d[d.q3 | d.q4] out = pd.DataFrame({"gaia_dr3_source_id": s.sid.values, "cls": s.cls.values, "dist_pc": s.dist_pc.round(2).values, "Gmag": s.Gmag.round(3).values, "M_G": s.M_G.round(3).values, "G_RP": s.g_rp.round(3).values, "Lbol_Lsun": s.Lbol.map(lambda v: float(f"{v:.3g}")).values, "q3": s.q3.astype(int).values, "q4": s.q4.astype(int).values, "E3_mJy": (1e3 * s.E3K).round(3).values, "sig3_mJy": (1e3 * s.sig3K).round(3).values, "chi3": s.chi3.round(2).values, "E4_mJy": (1e3 * s.E4K).round(3).values, "sig4_mJy": (1e3 * s.sig4K).round(3).values, "chi4": s.chi4.round(2).values, "det": s.det.astype(int).values}) out.to_csv(os.path.join(DATA, "stars_slim.csv.gz"), index=False) def main(): d = inj.load_stars() fig_significance(d) fig_completeness(d) fig_limits() fig_accounting() fig_cutouts(d) slim_table(d) print("figures written to", FIG) if __name__ == "__main__": main()