"""Figures for the seti_meta study (written to figures/*.png).""" import os import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import LogNorm from scipy.special import exp1 from common import DATA, FIGS, LGRID, L_ARECIBO, L_KI, L_KII, L_LTE import coverage as C import analysis as A # Okabe-Ito hues in an order that passes the adjacent-pair colour-vision checks; text stays in ink colours OI = {"black": "#000000", "orange": "#E69F00", "sky": "#56B4E9", "green": "#009E73", "yellow": "#F0E442", "blue": "#0072B2", "red": "#D55E00", "purple": "#CC79A7", "grey": "#7f7f7f"} INK, INK2 = "#1a1a1a", "#555555" plt.rcParams.update({"font.size": 10, "axes.grid": True, "grid.alpha": 0.22, "figure.dpi": 130, "axes.spines.top": False, "axes.spines.right": False, "legend.frameon": False, "text.color": INK, "axes.labelcolor": INK, "xtick.color": INK2, "ytick.color": INK2}) REFL = ((L_LTE, "Earth LTE\n4 GW"), (L_ARECIBO, "Arecibo\nradar"), (L_KI, "K-I"), (L_KII, "K-II")) def ref_lines(ax, y, va="bottom"): for L, name in REFL: ax.axvline(L, color=OI["grey"], lw=0.8, ls=":") ax.text(L * 1.25, y, name, fontsize=7.5, color=INK2, va=va) def fig1(V): fid = V["inertial"][A.BASE]; res = fid["res"]; NMW = fid["meta"]["n_mw_total"] eta = A.ETA_FID fig = plt.figure(figsize=(11.6, 6.4)) ax = fig.add_axes([0.075, 0.105, 0.60, 0.80]) kx = fig.add_axes([0.69, 0.105, 0.30, 0.80]); kx.axis("off") spots = [("1420.4 MHz (HI line)", 1420.4, OI["red"]), ("1.5 GHz (L band)", 1500.0, OI["orange"]), ("2.0 GHz (S band)", 2000.0, OI["green"]), ("6 GHz (C band)", 6000.0, OI["blue"]), ("9 GHz (X band)", 9000.0, OI["purple"]), ("120 MHz (MWA, LOFAR)", 120.0, OI["sky"])] with np.errstate(divide="ignore", invalid="ignore"): for name, f, col in spots: n = A.ncov_at(res, f) ax.plot(LGRID, np.where(n > 0, A.UL95 / (eta * n), np.nan), color=col, lw=1.7, label=name) n = A.ncov_at(V["inertial"][A.SAH]["res"], 1420.4) ax.plot(LGRID, np.where(n > 0, A.UL95 / (eta * n), np.nan), color=OI["red"], lw=1.3, ls="--", label="HI line, if SETI@home's candidates are all rejected") n = A.ncov_at(V["surface"][A.BASE]["res"], 1420.4) ax.plot(LGRID, np.where(n > 0, A.UL95 / (eta * n), np.nan), color=OI["red"], lw=1.3, ls=":", label="HI line, transmitter drifting like one on Earth's surface") n = fid["neff"]["log_1_10GHz"] ax.plot(LGRID, A.UL95 / (eta * n), color="k", lw=2.6, label="frequency unknown: log-uniform 1-10 GHz") # published results, plotted as published (3/N for a published star count N, i.e. efficiency 1) pub = [("Tao+22, FAST: 33 stars", 6.3e12, 3 / 33., "v"), ("Price+20, Parkes 10 cm: 189", 9.1e12, 3 / 189., "o"), ("Price+20, GBT L: 882", 2.1e12, 3 / 882., "o"), ("Price+20, GBT S: 1005", 2.1e12, 3 / 1005., "h"), ("Enriquez+17: 692", 1.0e13, 3 / 692., "s"), ("Wlodarczyk-Sroka+20, <50 pc: 1513", 1.0e13, 3 / 1513., "D"), ("Wlodarczyk-Sroka+20, <200 pc", 2.5e14, 3 * 3.9e-4, "D"), ("Margot+23, <100 pc: 47", 1.35e13, 0.066, "^"), ("Margot+23, <20 kly: 10230", 5.08e16, 3.0e-4, "^"), ("Mason+25, <2.5 kpc", 5.0e16, 9.95e-6, "P"), ("Tremblay+25, COSMIC: 511", 2.09e16, 3 / 511., "X"), ("Johnson+23, LOFAR: 4e5", 1e17, 3 / 4e5, "<"), ("Gajjar+21, GC with GBT C: 6e5", 5e17, 3 / 6e5, "*"), ("Gajjar+21, GC with Parkes: 6e7", 4e18, 3 / 6e7, "*")] # number offsets (points) chosen so that the numerals do not collide in the crowded 1e12-1e13 W region off = {1: (-4, 7), 2: (7, 2), 3: (-13, 3), 4: (-13, -9), 5: (8, -9), 6: (6, -12), 7: (7, 3), 8: (7, 5), 9: (7, 3), 10: (-15, -4), 11: (7, 3), 12: (8, -3), 13: (8, 0), 14: (8, 0)} for i, (name, L, f, m) in enumerate(pub, 1): ax.plot(L, f, m, ms=6.5, mfc="white", mec="k", mew=1.0, zorder=5) ax.annotate(str(i), (L, f), xytext=off[i], textcoords="offset points", fontsize=7.5, color=INK, zorder=6) ax.axhline(1.0, color="k", lw=0.8) ax.axhline(A.UL95 / (eta * NMW), color=OI["grey"], lw=1.0, ls="--") ax.text(2e28, A.UL95 / (eta * NMW) * 0.7, "every star in the Galaxy covered", fontsize=7.5, color=INK2, ha="right", va="top") ax.text(2e28, 1.4, "above this line: no constraint (limit > 100% of stars)", fontsize=7.5, color=INK2, ha="right") ref_lines(ax, 2.2e-12) ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlim(1e9, 3e28); ax.set_ylim(1.5e-12, 150) ax.set_xlabel("EIRP of a continuous narrowband transmitter (W)") ax.set_ylabel("95% upper limit on the fraction of stars with such a transmitter") fig.suptitle("What the null radio searches exclude, by transmitter power and frequency " "(curves: this work, 3 / (0.5 N$_{cov}$), SETI@home not counted unless stated)", fontsize=10, x=0.075, ha="left", y=0.965) h, l = ax.get_legend_handles_labels() kx.legend(h, l, loc="upper left", fontsize=8, title="Transmitter frequency (this work, pipeline efficiency 0.5)", title_fontsize=8, bbox_to_anchor=(0.0, 1.0), alignment="left") kx.text(0.02, 0.50, "Published single-survey limits (as published,\n3/N, efficiency 1)", fontsize=8, va="top", transform=kx.transAxes) for i, (name, L, f, m) in enumerate(pub, 1): kx.text(0.02, 0.50 - 0.033 * (i + 1.3), f"{i:>2} {name}", fontsize=7.6, va="top", transform=kx.transAxes, color=INK) fig.savefig(os.path.join(FIGS, "fig1_transmitter_fraction_limit.png"), dpi=160) plt.close(fig) def fig2(V): res = V["inertial"][A.BASE]["res"] iv = res["intervals"] fig, ax = plt.subplots(figsize=(9.6, 5.8)) Z = np.maximum(res["fid_mw"], 1e-3) norm = LogNorm(vmin=1, vmax=1e11) cmap = plt.get_cmap("viridis").copy() cmap.set_under("#cfc3dd") # searched, but fewer than one star within reach for i, (a, b) in enumerate(iv): if Z[i].max() <= 1e-3: continue # widen very narrow intervals so that they remain visible on a log-frequency axis a2, b2 = (a, b) if (b / a) > 1.006 else (np.sqrt(a * b) / 1.003, np.sqrt(a * b) * 1.003) m = ax.pcolormesh([a2, b2], LGRID, np.vstack([Z[i], Z[i]]).T[:-1, :1], norm=norm, cmap=cmap, shading="flat") cb = fig.colorbar(m, ax=ax, pad=0.02, extend="min") cb.set_label("N$_{cov}$: Milky Way stars for which a transmitter of this EIRP\nat this frequency would have been detected") for L, name in ((L_LTE, "Earth LTE 4 GW"), (L_ARECIBO, "Arecibo radar"), (L_KI, "Kardashev I"), (L_KII, "Kardashev II")): ax.axhline(L, color="w", lw=0.8, ls=":") ax.text(21000, L * 1.3, name, fontsize=7.5, color=INK, ha="right", va="bottom") for f, name in ((1420.4, "HI"), (430, "430 MHz\nSERENDIP III"), (150, "MWA / LOFAR"), (6000, "C"), (9500, "X"), (2400, "S")): ax.text(f, 2e28, name, fontsize=7.5, ha="center", va="bottom") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlim(90, 22000); ax.set_ylim(1e9, 1e28) ax.set_facecolor("#e9e9e9") ax.grid(False) ax.set_xlabel("Transmitter frequency (MHz)") ax.set_ylabel("EIRP (W)") ax.set_title("Frequency x power coverage of dedicated narrowband radio searches (SETI@home not counted)\n" "grey = never searched; lilac = searched, but less than one star within reach", fontsize=10, pad=26) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig2_coverage_map.png"), dpi=160) plt.close(fig) def fig3(V): lim = pd.read_csv(os.path.join(DATA, "limits_vs_eirp.csv")) beh = pd.read_csv(os.path.join(DATA, "behaviour_limits.csv")) lad = pd.read_csv(os.path.join(DATA, "hi_line_scenarios.csv")) NMW = lim["N_MW_model"].iloc[0] fig = plt.figure(figsize=(14.4, 6.3)) ax = fig.add_axes([0.06, 0.12, 0.42, 0.76]) bx = fig.add_axes([0.735, 0.17, 0.25, 0.71]) L = lim["EIRP_W"].values nmv = np.vectorize(lambda ne, eta: A.nmax_exact(eta * ne / NMW)) sty = [("HI_line", "21-cm line, tuned to our frame", OI["red"]), ("water_hole", "anywhere in the 'water hole' (1.42-1.67 GHz)", OI["orange"]), ("log_1_10GHz", "anywhere 1-10 GHz (log-uniform)", "k"), ("log_0p1_100GHz", "anywhere 0.1-100 GHz (log-uniform)", OI["blue"])] for p, name, col in sty: ax.plot(L, np.maximum(nmv(lim[f"Neff_{p}__base"].values, A.ETA_FID), 0.5), color=col, lw=2.4 if p == "log_1_10GHz" else 1.6, label=name) if p in ("log_1_10GHz", "HI_line"): lo = np.maximum(nmv(lim[f"Neff_{p}__max_over_variants"].values, A.ETA_RANGE[1]), 0.5) hi = nmv(lim[f"Neff_{p}__min_over_variants"].values, A.ETA_RANGE[0]) ax.fill_between(L, lo, hi, color=col, alpha=0.11, lw=0) ax.plot(L, np.maximum(nmv(lim["Neff_HI_line__sah"].values, A.ETA_FID), 0.5), color=OI["red"], lw=1.2, ls="--", label="21-cm line, if SETI@home's candidates are all rejected") ax.plot(L, np.maximum(nmv(lim["Neff_HI_line__surface"].values, A.ETA_FID), 0.5), color=OI["red"], lw=1.2, ls=":", label="21-cm line, transmitter drifting like one on Earth's surface") ax.axhline(NMW, color=OI["grey"], lw=1.2, ls="--") ax.fill_between([1e8, 1e29], NMW, 1e19, color="#bbbbbb", alpha=0.25, lw=0) ax.text(2e28, NMW * 2.2, "above this line: more civilizations than stars = no constraint", fontsize=8, color=INK2, ha="right") b = beh[(beh.freq_prior == "log_1_10GHz") & beh.default_row].set_index("behaviour") pts = [("ii_radar_random", "(ii) Arecibo-class beam,\nrandom pointing", OI["purple"], "v", (9, -4)), ("ii_radar_aimed", "(ii') same, aimed at the Sun", OI["purple"], "o", (9, -12)), ("iii_leak_narrow", "(iii) Earth-2025 leakage,\ngenerous narrowband bound", OI["green"], "s", (9, 2)), ("i_beacon_1e17W", "(i) isotropic beacon 1e17 W", "k", "D", (9, 4))] for bid, name, col, m, off in pts: ax.errorbar(b.loc[bid, "EIRP_W"], b.loc[bid, "Nmax95"], yerr=[[b.loc[bid, "Nmax95"] - b.loc[bid, "Nmax95_lo"]], [b.loc[bid, "Nmax95_hi"] - b.loc[bid, "Nmax95"]]], color=col, lw=1.2, capsize=3, zorder=5) ax.plot(b.loc[bid, "EIRP_W"], b.loc[bid, "Nmax95"], m, color=col, ms=8, mec="k", zorder=6) ax.annotate(name, (b.loc[bid, "EIRP_W"], b.loc[bid, "Nmax95"]), xytext=off, textcoords="offset points", fontsize=8, color=INK, ha="right" if off[0] < 0 else "left") ax.text(2e28, 3e9, "(iii) literal Earth leakage (4 GW over an assumed\n20 MHz): undetectable beyond ~430 AU by any\nsearch in the record, so N is unbounded", fontsize=8, color=INK, ha="right", va="top") ref_lines(ax, 0.62) ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlim(1e9, 3e28); ax.set_ylim(0.5, 3e18) ax.set_xlabel("EIRP per civilization (W)") ax.set_ylabel("95% upper limit on the number of transmitting\ncivilizations in the Milky Way") ax.set_title("(a) Maximum number of transmitting civilizations vs. power\n(isotropic, continuous, narrowband; bands: all in-range variants x efficiency 0.25-1)", fontsize=9.5, loc="left", pad=8) ax.legend(loc="upper right", fontsize=7.6, bbox_to_anchor=(1.0, 0.97)) # ---- (b) what the 1e17 W limit depends on short = ["A baseline: legacy sky surveys at nominal\n thresholds, drift-free transmitter, no SETI@home", "B A + SETI@home counted as a null", "C SETI@home + modern surveys only", "D A, transmitter with Earth-like orbital drift", "E A, transmitter on an Earth-like rotating surface", "F E + SETI@home counted as a null", "G META / META II thresholds x10", "H META / META II thresholds x30", "I META / META II thresholds x60", "J modern surveys only, no SETI@home", "K no sky surveys at all"] y = np.arange(len(lad))[::-1] ser = [("Nmax95_HI_line_1e17W", "21-cm line, tuned to our frame", OI["red"], "o"), ("Nmax95_HI_rest_1e17W", "21-cm line in the star's rest frame", OI["orange"], "s"), ("Nmax95_log_1_10GHz_1e17W", "anywhere 1-10 GHz", "k", "D")] for yy in y: bx.axhline(yy, color="#dddddd", lw=0.6, zorder=0) for (col, name, c, m), dy_ in zip(ser, (0.17, -0.17, 0.0)): bx.plot(lad[col].values, y + dy_, m, color=c, ms=7.0, mec="white", mew=0.8, label=name, zorder=3, ls="none") for yy, v in zip(y, lad["Nmax95_HI_line_1e17W"].values): bx.annotate(f"{v:,.0f}" if v < 1e4 else f"{v:.1e}".replace("e+0", "e"), (v, yy + 0.17), xytext=(0, 6), textcoords="offset points", fontsize=7.2, color=INK, ha="center") bx.set_yticks(y); bx.set_yticklabels(short, fontsize=7.8, ha="left") bx.tick_params(axis="y", pad=238, length=0) bx.set_xscale("log"); bx.set_xlim(1e3, 2e6); bx.set_ylim(-0.7, len(lad) - 0.3) bx.grid(axis="y", visible=False) bx.set_xlabel("95% upper limit on the number of 1e17 W beacons") fig.text(0.50, 0.905, "(b) What the limit on 1e17 W isotropic beacons depends on (efficiency 0.5; number shown: 21-cm line, our frame)", fontsize=9.5, ha="left") bx.legend(fontsize=7.8, loc="upper right", bbox_to_anchor=(1.0, -0.13), ncol=3, columnspacing=1.0, handletextpad=0.3) fig.savefig(os.path.join(FIGS, "fig3_max_civilizations.png"), dpi=160) plt.close(fig) def fig4(V): beh = pd.read_csv(os.path.join(DATA, "behaviour_limits.csv")) NMW = V["inertial"][A.BASE]["meta"]["n_mw_total"] n = np.logspace(0, np.log10(NMW), 400) fig, ax = plt.subplots(figsize=(9.0, 5.6)) ax.plot(n, np.log(NMW / n) / np.log(NMW), color=OI["grey"], lw=2.5, ls="--", label="prior (log-uniform, 1 to 7.5e10)") sel = [("i_beacon_KII", "HI_line", False, "(i) isotropic K-II beacon, 21-cm line", OI["red"], "-"), ("i_beacon_1e17W", "HI_line", True, "(i) 1e17 W beacon, 21-cm line, SETI@home counted", OI["orange"], "--"), ("i_beacon_1e17W", "HI_line", False, "(i) 1e17 W beacon, 21-cm line, baseline", OI["orange"], "-"), ("i_beacon_1e17W", "log_1_10GHz", False, "(i) 1e17 W beacon, anywhere 1-10 GHz", "k", "-"), ("i_beacon_2e13W", "log_1_10GHz", False, "isotropic 2e13 W beacon, 1-10 GHz", OI["blue"], "-"), ("ii_radar_aimed", "log_1_10GHz", False, "(ii') Arecibo-class beam aimed at us, 5% duty", OI["purple"], "-"), ("iii_leak_narrow", "log_1_10GHz", False, "(iii) Earth leakage, generous bound", OI["green"], "-"), ("ii_radar_random", "log_1_10GHz", False, "(ii) Arecibo-class beams, random pointing", OI["purple"], ":")] for bid, p, sah, name, col, ls in sel: m = (beh.behaviour == bid) & (beh.freq_prior == p) & (beh.setiathome_counted == sah) m &= beh.default_row if not sah else (beh.transmitter_drift == "inertial") pb = float(beh[m]["p_detect_per_civ"].iloc[0]) lam = -np.log1p(-min(pb, 1 - 1e-12)) if lam * NMW < 1e-9: sf = np.log(NMW / n) / np.log(NMW) else: a, b_ = exp1(lam * 1.0), exp1(lam * NMW) sf = (exp1(lam * n) - b_) / (a - b_) ax.plot(n, sf, color=col, lw=1.8, ls=ls, label=name) ax.axhline(0.05, color="k", lw=0.6, ls=":") ax.text(1.3, 0.065, "5%", fontsize=8) ax.set_xscale("log") ax.set_xlim(1, NMW); ax.set_ylim(0, 1.02) ax.set_xlabel("N (number of civilizations of this type currently transmitting in the Milky Way)") ax.set_ylabel("Probability that the number exceeds N") ax.set_title("Posterior after all nulls vs. prior: only loud, always-on, isotropic transmitters are constrained", fontsize=10) ax.legend(fontsize=7.6, loc="upper right") fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig4_posteriors.png"), dpi=160) plt.close(fig) def fig5(): dy = pd.read_csv(os.path.join(DATA, "dyson_limits.csv")) d = dy[dy.T_DS_K == 300.0] fig, axs = plt.subplots(1, 2, figsize=(11.2, 5.6), gridspec_kw={"width_ratios": [1.3, 1]}) ax = axs[0] x = np.arange(3) w = 0.15 series = [("Suazo2022_published_frac_5kpc", "published upper limit, < 5 kpc", "#c6dbef", None), ("Suazo2022_published_frac_1kpc", "published upper limit, < 1 kpc", "#9ecae1", None), ("Suazo2022_published_frac_300pc", "published upper limit, < 300 pc", "#6baed6", None), ("Suazo2022_published_frac_100pc", "published upper limit, < 100 pc", "#3182bd", None), ("floor_frac_300pc_if_100pct_complete_k0", "sensitivity floor, < 300 pc (NOT a limit:\nassumes a 100% complete search)", "white", "////")] for j, (col, name, c, hatch) in enumerate(series): ax.bar(x + (j - 2) * (w + 0.012), d[col].values, w, color=c, label=name, edgecolor="k", lw=0.5, hatch=hatch) for xi, v in zip(x, d["Suazo2022_published_frac_300pc"].values): ax.text(xi, v * 1.35, f"{v:.1e}".replace("e-0", "e-"), fontsize=7.5, ha="center", color=INK) ax.set_yscale("log"); ax.set_xticks(x); ax.set_xticklabels(["> 10%", "> 50%", "> 90%"]) ax.set_xlabel("covering fraction of a 300 K swarm") ax.set_ylabel("fraction of stars") ax.set_ylim(3e-7, 1.0) ax.grid(axis="x", visible=False) ax.legend(fontsize=7.4, loc="upper center", ncol=2, bbox_to_anchor=(0.5, -0.15)) ax.set_title("Dyson swarms: published limits (Suazo et al. 2022)\nand the floor a complete WISE search could reach", fontsize=9.5) ax = axs[1] for T, c in ((100.0, OI["blue"]), (300.0, "k"), (600.0, OI["red"])): dd = dy[dy.T_DS_K == T] ax.plot(dd.covering_fraction, dd.N_sensitive_100pc_GCNS / dd.N_stars_100pc_GCNS, "o-", color=c, label=f"< 100 pc (GCNS), T = {T:.0f} K") ax.plot(dd.covering_fraction, dd.N_sensitive_300pc_model / dd.N_total_300pc_model, "s--", color=c, mfc="white", label=f"< 300 pc (model), T = {T:.0f} K") ax.set_ylim(0, 1.45); ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0]); ax.set_xticks([0.1, 0.5, 0.9]) ax.set_xlabel("covering fraction") ax.set_ylabel("fraction of stars for which the swarm would be\na >5-sigma WISE W3+W4 source with Gaia G < 19") ax.legend(fontsize=7, loc="upper center", ncol=2) ax.set_title("Flux-limit sensitivity of the WISE data\n(not the completeness of any published pipeline)", fontsize=9.5) fig.tight_layout() fig.savefig(os.path.join(FIGS, "fig5_dyson_limits.png"), dpi=160) plt.close(fig) def main(V=None): if V is None: V = A.run_variants() fig1(V); fig2(V); fig3(V); fig4(V); fig5() print("[figures] wrote", sorted(f for f in os.listdir(FIGS) if f.endswith(".png"))) if __name__ == "__main__": main()