"""Limits, behaviour scenarios, robustness variants, Bayesian posteriors and Dyson-swarm numbers. Everything here is derived from coverage.build()/combine() (radio) and from the GCNS + WISE sensitivity model (waste heat). Outputs go to data/*.csv and data/key_numbers.json. Conventions of the revised analysis * Baseline record = every published null EXCEPT SETI@home (whose ~200 candidates await re-observation), transmitter drift class 'inertial', in-beam star counts calibrated to Mason+ 2025, eta = 0.5. SETI@home enters only in the scenarios labelled 'SETI@home counted as a null'. * Drift classes: 'inertial' (transmitter removes its own acceleration), 'orbit', 'surface' (Earth-like). * Limits use the exact binomial form: P(all nulls | N) = (1 - p)^N, N_max(95%) = ln(0.05) / ln(1 - p). * Values at a given EIRP are log-log interpolated on the 0.125-dex grid (no nearest-grid-point rounding). """ import json import os import numpy as np import pandas as pd from scipy.optimize import brentq from scipy.special import exp1 from scipy.stats import chi2 from common import DATA, RAW, LGRID, PC, FOURPI, L_ARECIBO, L_KI, L_KII, L_LTE, dmax_pc import coverage as C import galaxy as G ETA_FID = 0.5 # probability that a signal above threshold survives RFI excision / the pipeline ETA_RANGE = (0.25, 1.0) # (Price+ 2020 adopt 0.5 for their limits; SETI@home birdie recovery >= 0.8) CL = 0.95 UL95 = -np.log(1.0 - CL) # 95% Poisson upper limit for zero detections (2.996, the 'rule of three') DRIFTS = ("inertial", "orbit", "surface") BASE = "baseline (no SETI@home; in-beam counts calibrated to Mason+ 2025)" SAH = "SETI@home counted as a null (all pending candidates rejected)" MODERN = "modern surveys only (legacy tier-B surveys dropped; no SETI@home)" MODERN_SAH = "modern surveys + SETI@home (no legacy surveys)" NOSKY = "no sky surveys at all" MODERN_GC = "modern surveys only, Galactic-centre fields scaled to authors' star counts" def at(arr, L): """log-log interpolation of a function tabulated on LGRID.""" if L < LGRID[0]: return 0.0 v = 10 ** np.interp(np.log10(L), np.log10(LGRID), np.log10(np.maximum(arr, 1e-300))) return float(v) if v > 1e-250 else 0.0 def nmax_exact(p, cl=CL): """Largest N consistent with zero detections at confidence cl when each civilisation is detected with probability p independently: (1-p)^N = 1-cl.""" if p <= 0: return np.inf if p >= 1: return 0.0 return float(np.log(1.0 - cl) / np.log1p(-p)) BL_COMPS = ("bl_gbtL_targets", "bl_gbtL_bycatch", "bl_gbtS_targets", "bl_gbtS_bycatch", "bl_pks10_targets", "bl_pks10_bycatch", "bl_gbtC_archive", "bl_gbtX_archive", "bl_gc_parkes", "bl_gc_gbtC", "choza_gal", "uno_bg") # ------------------------------------------------------------------ robustness variants def variants(): """name -> dict(kw=build kwargs, sah=include SETI@home, inc=extra include filter, which=union, in_range=counted in quoted ranges, only_inertial=skip for drifting transmitters)""" A = lambda k, c: c["tier"] in ("A", "X") metas = lambda f: dict(smin_scale={"meta": f, "meta2": f, "meta_h2": f}) v = { BASE: dict(), SAH: dict(sah=True), "SETI@home counted, sky fraction 25% instead of 30%": dict(sah=True, kw=dict(n_scale={"setiathome": 25.0 / 30.0})), "SETI@home counted, union upper bound (plain sum)": dict(sah=True, which="hi"), "union lower bound (max over surveys)": dict(which="lo"), "union upper bound (plain sum)": dict(which="hi"), MODERN: dict(inc=A), MODERN_SAH: dict(inc=A, sah=True), MODERN_GC: dict(inc=A, kw=dict(n_scale={"bl_gc_parkes": 6e7 / 2.761e8, "bl_gc_gbtC": 6e5 / 3.651e6})), NOSKY: dict(inc=lambda k, c: c["group"] != "sky"), "named targets only (no in-beam stars, no fields)": dict(inc=lambda k, c: c["group"] in ("near", "cosmic", "fast_gc"), in_range=False), "tier B thresholds x10 (nominal -> confirmable)": dict(kw=dict(smin_scale={"tier_B": 10.0})), "all thresholds x2": dict(kw=dict(smin_scale={"all": 2.0})), "saturated scintillation weighting": dict(kw=dict(scint=True)), "Phoenix Arecibo at 1e-26 W/m2 (Technosearch value)": dict(kw=dict(smin_scale={"phoenix_arecibo": 1.0 / 16.0})), "Galactic-centre fields scaled to authors' star counts": dict(kw=dict(n_scale={"bl_gc_parkes": 6e7 / 2.761e8, "bl_gc_gbtC": 6e5 / 3.651e6})), "in-beam counts: raw analytic model (first-pass fiducial)": dict(kw=dict(bycatch_cal="none")), "in-beam counts: calibration omitting the two densest pointings": dict(kw=dict(bycatch_cal="mason_excl2")), "in-beam counts x0.5": dict(kw=dict(n_scale={"group_bl_beam": 0.5})), # Margot+ 2023 (AJ 166, 206, Sect. 5.3): a BL-like pipeline recovers only 12.7% of injected signals when # drift rates are spread uniformly over +-4 Hz/s; relative to the fiducial eta = 0.5 this is a factor 0.254 "BL pipeline efficiency 12.7% (Margot+ 2023)": dict(kw=dict(n_scale={k: 0.127 / 0.5 for k in BL_COMPS})), "META/META II thresholds x10 (drift smearing)": dict(kw=metas(10.0), only_inertial=True), "META/META II thresholds x30 (drift smearing)": dict(kw=metas(30.0), only_inertial=True), "META/META II thresholds x60 (drift smearing)": dict(kw=metas(60.0), only_inertial=True), "META/META II x30, SETI@home counted (non-barycentric column)": dict(kw=dict(smin_scale={"meta": 30.0, "meta2": 30.0, "meta_h2": 30.0}), sah=True, only_inertial=True, sah_nonbary=True, in_range=False), } return v def run_variants(): """V[drift][name] = dict(neff={prior: N_eff(L)}, in_range, ...). 'orbit' is run for the two main records only.""" out = {d: {} for d in DRIFTS} for drift in DRIFTS: for name, spec in variants().items(): if drift != "inertial" and spec.get("only_inertial"): continue if drift == "orbit" and name not in (BASE, SAH): continue kw = dict(spec.get("kw", {})) comps, meta = C.build(drift=drift, **kw) if spec.get("sah_nonbary"): # take the SETI@home component from the non-barycentric build comps["setiathome"] = C.build(drift="surface", **{k: v for k, v in kw.items() if k != "smin_scale"})[0]["setiathome"] extra, sah = spec.get("inc"), spec.get("sah", False) inc = lambda k, c, extra=extra, sah=sah: (sah or k != "setiathome") and (extra is None or extra(k, c)) which = spec.get("which", "fid") ne, res = C.neff_all(comps, include=inc, n_cap=meta["n_mw_total"], which=which + "_mw", priors=list(C.PRIORS) + list(C.AUX_PRIORS)) out[drift][name] = dict(neff=ne, in_range=spec.get("in_range", True), sah=sah, which=which, res=res if name in (BASE, SAH) else None, comps=comps if name == BASE else None, meta=meta if name == BASE else None) return out def envelope(V, p, drifts=("inertial", "surface")): arrs = [v["neff"][p] for d in drifts for v in V[d].values() if v["in_range"]] return np.min(arrs, axis=0), np.max(arrs, axis=0) def ncov_at(res, f_mhz, which="fid_mw"): for i, (a, b) in enumerate(res["intervals"]): if a <= f_mhz < b: return res[which][i] return np.zeros(len(LGRID)) # ------------------------------------------------------------------ Bayesian posterior def posterior_quantiles(p, n_lo, n_hi, qs=(0.05, 0.5, 0.95)): """Posterior of N for likelihood (1-p)^N = exp(-N lam), lam = -ln(1-p), and a log-uniform prior on [n_lo, n_hi]. CDF(n) = [E1(lam n_lo) - E1(lam n)] / [E1(lam n_lo) - E1(lam n_hi)].""" lam = -np.log1p(-min(p, 1 - 1e-12)) if p > 0 else 0.0 if lam * n_hi < 1e-9: # likelihood flat over the prior range: posterior = prior return [n_lo * (n_hi / n_lo) ** q for q in qs] a, b = exp1(lam * n_lo), exp1(lam * n_hi) cdf = lambda ln: (a - exp1(lam * np.exp(ln))) / (a - b) return [float(np.exp(brentq(lambda x: cdf(x) - q, np.log(n_lo), np.log(n_hi)))) for q in qs] def posterior_prob_above(p, n_lo, n_hi, n): lam = -np.log1p(-min(p, 1 - 1e-12)) if p > 0 else 0.0 if lam * n_hi < 1e-9: return float(np.log(n_hi / n) / np.log(n_hi / n_lo)) a, b = exp1(lam * n_lo), exp1(lam * n_hi) return float((exp1(lam * n) - b) / (a - b)) # ------------------------------------------------------------------ Dyson-swarm sensitivity H, KB, CLIGHT, SIG = 6.62607015e-34, 1.380649e-23, 2.99792458e8, 5.670374419e-8 LSUN = 3.828e26 W_LAM = {"W3": 11.56e-6, "W4": 22.09e-6} W_ZP = {"W3": 31.674, "W4": 8.363} # Jy, Jarrett et al. (2011) W_LIM = {"W3": 1.0e-3, "W4": 6.0e-3} # Jy, 5-sigma point-source sensitivity (Wright et al. 2010) # Suazo+ 2024 candidates: fitted covering fraction and temperature (their Table 5) and follow-up status SUAZO24 = {"A": (0.08, 138, "radio counterpart (Ren+ 2024)"), "B": (0.06, 275, "radio counterpart (Ren+ 2024); contaminated (Ren+ 2026)"), "C": (0.14, 187, "contaminated in archival imaging (Ren+ 2026), not spectroscopically refuted"), "D": (0.16, 178, "background galaxy, JWST/MIRI (Zackrisson+ 2026)"), "E": (0.08, 180, "background galaxy, JWST/MIRI (Zackrisson+ 2026)"), "F": (0.03, 137, "open"), "G": (0.13, 100, "radio counterpart (Ren+ 2024), not refuted")} K_SURV_10 = 2 # candidates with fitted covering fraction > 0.1 that are not yet refuted (C and G) def bb_fnu_per_lbol(lam, T): """F_nu / F_bol (1/Hz) for a blackbody of temperature T at wavelength lam.""" nu = CLIGHT / lam B = 2 * H * nu ** 3 / CLIGHT ** 2 / np.expm1(H * nu / (KB * T)) return np.pi * B / (SIG * T ** 4) def dyson_table(n_allsky_shells): """Number of stars for which a T_DS blackbody swarm of covering fraction gamma would be a >=5 sigma source in BOTH W3 and W4 (the Suazo+ 2024 requirement), with an excess >=30% of the W3 photosphere, and with the dimmed star still brighter than G=19. This is a SENSITIVITY count. Dividing 3 by it gives a floor that a search could reach if it were 100% complete, which Suazo+ 2024 do not claim (their pipeline rejects variable sources and 328 of 368 high-S/N sources as blends); the published upper limits remain the photometric-compatibility fractions of Suazo+ 2022 (columns Suazo2022_published_frac_*).""" g = pd.read_csv(os.path.join(RAW, "vizier", "gcns_table1c.csv")) rel = pd.read_csv(os.path.join(RAW, "gaia_mg_loglum_relation.csv")) rel = rel[rel["n"] >= 5] ok = (g["GCNSprob"] > 0.5) & ~(g["WDprob"] > 0.5) & g["Gmag"].notna() & (g["Dist50"] <= 0.1) g = g[ok].copy() d = g["Dist50"].values * 1000.0 MG = g["Gmag"].values - 5 * np.log10(d / 10.0) x, y = rel["mg_bin"].values + 0.25, rel["mean_loglum"].values slope_hi = (y[-1] - y[-5]) / (x[-1] - x[-5]) logL = np.interp(MG, x, y) logL = np.where(MG > x[-1], y[-1] + slope_hi * (MG - x[-1]), logL) # linear extrapolation for late M logL = np.where(MG < x[0], y[0], logL) Lw = 10 ** logL * LSUN phot_w3 = np.where(g["W3mag"].notna(), W_ZP["W3"] * 10 ** (-0.4 * g["W3mag"].values), 0.0) rows = [] # luminosity-function weights from the volume-complete 50 pc sample, for the extension beyond 100 pc in50 = d <= 50.0 mg_edges = np.arange(-2.0, 20.01, 0.5) lf, _ = np.histogram(MG[in50], bins=mg_edges) lf = lf / lf.sum() mg_mid = 0.5 * (mg_edges[:-1] + mg_edges[1:]) logL_mid = np.interp(mg_mid, x, y) logL_mid = np.where(mg_mid > x[-1], y[-1] + slope_hi * (mg_mid - x[-1]), logL_mid) cum_allsky = np.concatenate([[0.0], np.cumsum(n_allsky_shells)]) n_within = lambda dd: np.interp(np.log(np.maximum(dd, G.D_EDGES[0])), np.log(G.D_EDGES), cum_allsky) for T in (300.0, 100.0, 600.0): k3, k4 = bb_fnu_per_lbol(W_LAM["W3"], T), bb_fnu_per_lbol(W_LAM["W4"], T) for gam in (0.1, 0.5, 0.9): f3 = gam * Lw / (FOURPI * (d * PC) ** 2) * k3 * 1e26 # Jy f4 = gam * Lw / (FOURPI * (d * PC) ** 2) * k4 * 1e26 gdim = g["Gmag"].values - 2.5 * np.log10(1 - gam) det = (f3 >= np.maximum(W_LIM["W3"], 0.3 * phot_w3)) & (f4 >= W_LIM["W4"]) & (gdim < 19.0) n100 = int(det.sum()) fgk = MG < 7.5 # extension to 300 pc: luminosity function x model star counts, same flux criteria Lmid = 10 ** logL_mid * LSUN dmax3 = np.sqrt(gam * Lmid * k3 * 1e26 / (FOURPI * W_LIM["W3"])) / PC dmax4 = np.sqrt(gam * Lmid * k4 * 1e26 / (FOURPI * W_LIM["W4"])) / PC dmaxG = 10 ** ((19.0 + 2.5 * np.log10(1 - gam) - mg_mid + 5) / 5.0) dlim = np.minimum(np.minimum(dmax3, dmax4), dmaxG) n_sens = {D: float(np.sum(lf * n_within(np.minimum(dlim, D)))) for D in (100.0, 300.0, 1000.0)} n_tot = {D: float(n_within(D)) for D in (100.0, 300.0, 1000.0)} rows.append(dict(T_DS_K=T, covering_fraction=gam, N_stars_100pc_GCNS=int(len(g)), N_sensitive_100pc_GCNS=n100, N_sensitive_100pc_FGK=int((det & fgk).sum()), N_FGK_100pc=int(fgk.sum()), N_sensitive_100pc_model=n_sens[100.0], N_total_100pc_model=n_tot[100.0], N_sensitive_300pc_model=n_sens[300.0], N_total_300pc_model=n_tot[300.0], N_sensitive_1kpc_model=n_sens[1000.0], N_total_1kpc_model=n_tot[1000.0])) t = pd.DataFrame(rows) # published upper limits: 'fraction of stars whose photometry is compatible with a 300 K Dyson sphere' (Suazo+ 2022 Table 1) pub = {0.1: (6.6e-3, 4.8e-2, 1.3e-1, 1.6e-1), 0.5: (1.9e-4, 4.8e-4, 1.0e-3, 7.1e-3), 0.9: (1.8e-5, 8.6e-5, 2.1e-4, 8.3e-4)} for j, D in enumerate(("100pc", "300pc", "1kpc", "5kpc")): t[f"Suazo2022_published_frac_{D}"] = [pub[g_][j] if T_ == 300.0 else np.nan for T_, g_ in zip(t["T_DS_K"], t["covering_fraction"])] # sensitivity floors (NOT limits): 95% Poisson bound for k surviving candidates / number of sensitive stars. # k = 0 everywhere; k = 2 only for the >10% bin (candidates C and G; no candidate has gamma > 0.16). ul0, ul2 = UL95, 0.5 * chi2.ppf(CL, 2 * (K_SURV_10 + 1)) # 2.996, 6.296 for D, col in ((100, "N_sensitive_100pc_GCNS"), (300, "N_sensitive_300pc_model")): t[f"floor_frac_{D}pc_if_100pct_complete_k0"] = ul0 / t[col] t["floor_frac_300pc_if_100pct_complete_k2"] = np.where(np.isclose(t["covering_fraction"], 0.1), ul2 / t["N_sensitive_300pc_model"], np.nan) return t # ------------------------------------------------------------------ main def main(): V = run_variants() fid = V["inertial"][BASE] res, meta, comps = fid["res"], fid["meta"], fid["comps"] NMW = meta["n_mw_total"] P = C.prep() keyL = [1e10, 1e11, 1e12, 1e13, L_ARECIBO, 1e14, 1e15, 1e16, 1e17, L_KI, 1e18, 1e19, 1e20, 1e22, 1e24, 1e26, L_KII] nm = lambda ne, eta=ETA_FID: nmax_exact(eta * ne / NMW) cfg = {"base": ("inertial", BASE), "sah": ("inertial", SAH), "orbit": ("orbit", BASE), "orbit_sah": ("orbit", SAH), "surface": ("surface", BASE), "surface_sah": ("surface", SAH), "modern": ("inertial", MODERN)} NE = {k: V[d][n]["neff"] for k, (d, n) in cfg.items()} # ---- (a) limits vs EIRP (baseline record) lim = pd.DataFrame({"EIRP_W": LGRID}) spot = {"1420.4MHz_HI": 1420.4, "1500MHz_Lband": 1500.0, "2000MHz_Sband": 2000.0, "3000MHz": 3000.0, "6000MHz_Cband": 6000.0, "9000MHz_Xband": 9000.0, "430MHz": 430.0, "120MHz": 120.0} for name, f in spot.items(): lim[f"Ncov_{name}"] = ncov_at(res, f) lim["Ncov_1420.4MHz_HI_withSAH"] = ncov_at(V["inertial"][SAH]["res"], 1420.4) lim["Ncov_1420.4MHz_HI_surface_drift"] = ncov_at(V["surface"][BASE]["res"], 1420.4) for p in C.PRIORS: for k in cfg: lim[f"Neff_{p}__{k}"] = NE[k][p] lo, hi = envelope(V, p) lim[f"Neff_{p}__min_over_variants"] = lo lim[f"Neff_{p}__max_over_variants"] = hi lim[f"frac_UL95_{p}"] = [min(nm(x) / NMW, np.inf) for x in NE["base"][p]] lim[f"Nmax95_MW_{p}"] = [nm(x) for x in NE["base"][p]] lim["N_MW_model"] = NMW lim.to_csv(os.path.join(DATA, "limits_vs_eirp.csv"), index=False, float_format="%.5g") # ---- robustness table rob = [] for drift in ("inertial", "surface"): for name, v in V[drift].items(): for p in C.PRIORS: row = dict(variant=name, transmitter_drift=drift, prior=p, in_quoted_range=v["in_range"]) for L, tag in ((1e12, "1e12"), (L_ARECIBO, "2e13"), (1e15, "1e15"), (1e17, "1e17"), (1e20, "1e20"), (L_KII, "KII")): row[f"Neff_L{tag}"] = at(v["neff"][p], L) row["Nmax95_isotropic_1e17W"] = nm(at(v["neff"][p], 1e17)) row["Nmax95_isotropic_KII"] = nm(at(v["neff"][p], L_KII)) rob.append(row) rob = pd.DataFrame(rob) rob.to_csv(os.path.join(DATA, "robustness_variants.csv"), index=False, float_format="%.4g") # ---- scenario ladder for the HI line and the 1-10 GHz prior ladder_spec = [ ("A. baseline: legacy sky surveys at nominal thresholds, transmitter removes its own drift, no SETI@home", "inertial", BASE), ("B. A + SETI@home counted as a null (barycentric sensitivities)", "inertial", SAH), ("C. SETI@home + modern surveys only (legacy surveys not trusted)", "inertial", MODERN_SAH), ("D. A, but transmitter with Earth-like orbital drift (0.028 Hz/s at 1.42 GHz)", "orbit", BASE), ("E. A, but transmitter on an Earth-like rotating surface (0.16 Hz/s at 1.42 GHz)", "surface", BASE), ("F. E + SETI@home counted as a null (non-barycentric sensitivities)", "surface", SAH), ("G. META/META II thresholds x10, no SETI@home", "inertial", "META/META II thresholds x10 (drift smearing)"), ("H. META/META II thresholds x30, no SETI@home", "inertial", "META/META II thresholds x30 (drift smearing)"), ("I. META/META II thresholds x60, no SETI@home", "inertial", "META/META II thresholds x60 (drift smearing)"), ("J. modern surveys only, no SETI@home", "inertial", MODERN), ("K. no sky surveys at all", "inertial", NOSKY), ] lad = [] for desc, d, n in ladder_spec: ne = V[d][n]["neff"] row = dict(scenario=desc) for p in ("HI_line", "HI_rest", "log_1_10GHz"): for L, tag in ((1e17, "1e17W"), (1e20, "1e20W"), (L_KII, "KII")): x = at(ne[p], L) row[f"Neff_{p}_{tag}"] = x row[f"p_{p}_{tag}"] = ETA_FID * x / NMW row[f"Nmax95_{p}_{tag}"] = nm(x) lad.append(row) lad = pd.DataFrame(lad) lad.to_csv(os.path.join(DATA, "hi_line_scenarios.csv"), index=False, float_format="%.4g") # ---- (b) behaviours G_AR = 2.0e7 # Arecibo S-band gain: EIRP 2e13 W / 1e6 W transmitter (73 dBi) FB_AR = 1.0 / G_AR # fraction of the sky inside the beam DUTY = 0.05 # assumed fraction of time a given radar is transmitting # LTE-like leakage: Saide+ 2023 give ~4 GW peak summed over all LTE bands. ASSUMPTION: it is spread over only # 20 MHz (the observing bandwidth they consider); the true emission is spread over more spectrum, so the # power inside one search channel is even lower than computed here. LTE_BW = 20e6 chan = 2.79 L_lte_chan = L_LTE * chan / LTE_BW smin_best = 1.09e-26 # FAST, the most sensitive threshold in the record (W/m2) d_lte_pc = float(dmax_pc(L_lte_chan, smin_best)) beh = [ # id, description, L, duty, beam fraction, default drift class ("i_beacon_1e17W", "isotropic continuous narrowband beacon, EIRP 1e17 W", 1e17, 1.0, 1.0, "inertial"), ("i_beacon_KI", "isotropic continuous narrowband beacon, EIRP 1.74e17 W (Kardashev I)", L_KI, 1.0, 1.0, "inertial"), ("i_beacon_1e20W", "isotropic continuous narrowband beacon, EIRP 1e20 W", 1e20, 1.0, 1.0, "inertial"), ("i_beacon_KII", "isotropic continuous narrowband beacon, EIRP 3.8e26 W (Kardashev II)", L_KII, 1.0, 1.0, "inertial"), ("i_beacon_2e13W", "isotropic continuous narrowband beacon, EIRP 2e13 W", L_ARECIBO, 1.0, 1.0, "inertial"), ("ii_radar_random", "Arecibo-class beam (EIRP 2e13 W, gain 2e7) pointed at random, duty cycle 5%", L_ARECIBO, DUTY, FB_AR, "surface"), ("ii_radar_aimed", "Arecibo-class beam aimed at the Sun, duty cycle 5%", L_ARECIBO, DUTY, 1.0, "surface"), ("iii_leak_narrow", "Earth-like leakage, generous: all 4 GW in one isotropic continuous narrowband carrier", L_LTE, 1.0, 1.0, "surface"), ("iii_leak_literal", "Earth-like leakage, literal: 4 GW spread over an assumed 20 MHz (what a 3-Hz channel sees)", L_lte_chan, 1.0, 1.0, "surface"), ] rows, post = [], [] for bid, desc, L, duty, fb, ddef in beh: for p in C.PRIORS: for drift in DRIFTS: rlo, rhi = envelope(V, p, drifts=(drift,) if drift != "orbit" else ("inertial", "surface")) ne_lo, ne_hi = at(rlo, L), at(rhi, L) for sah in (False, True): ne = at(V[drift][SAH if sah else BASE]["neff"][p], L) pbar = ETA_FID * duty * fb * ne / NMW p_lo = ETA_RANGE[0] * duty * fb * ne_lo / NMW p_hi = ETA_RANGE[1] * duty * fb * ne_hi / NMW nmax = nmax_exact(pbar) is_def = (drift == ddef) and (not sah) rows.append(dict(behaviour=bid, description=desc, freq_prior=p, transmitter_drift=drift, setiathome_counted=sah, default_row=is_def, EIRP_W=L, duty_cycle=duty, beam_fraction=fb, eta=ETA_FID, Neff_stars=ne, Neff_min_variants=ne_lo, Neff_max_variants=ne_hi, p_detect_per_civ=pbar, p_detect_lo=p_lo, p_detect_hi=p_hi, Nmax95=nmax, Nmax95_lo=nmax_exact(p_hi), Nmax95_hi=nmax_exact(p_lo), Nmax_over_Nstars=nmax / NMW, expected_detections_if_every_star_hosts_one=pbar * NMW)) if drift == ddef: for prior_name, nlo in (("loguniform_1_to_Nstars", 1.0), ("loguniform_1e-3_to_Nstars", 1e-3)): q05, q50, q95 = posterior_quantiles(pbar, nlo, NMW) pr = [nlo * (NMW / nlo) ** q for q in (0.05, 0.5, 0.95)] post.append(dict(behaviour=bid, freq_prior=p, transmitter_drift=drift, setiathome_counted=sah, N_prior=prior_name, p_detect_per_civ=pbar, prior_median=pr[1], prior_q95=pr[2], post_median=q50, post_q95=q95, post_P_N_gt_100=posterior_prob_above(pbar, nlo, NMW, 100.0), post_P_N_gt_1e4=posterior_prob_above(pbar, nlo, NMW, 1e4), post_P_N_gt_1e6=posterior_prob_above(pbar, nlo, NMW, 1e6), prior_P_N_gt_1e6=float(np.log(NMW / 1e6) / np.log(NMW / nlo)), dex_removed_from_q95=float(np.log10(pr[2] / q95)))) beh_t = pd.DataFrame(rows) beh_t.to_csv(os.path.join(DATA, "behaviour_limits.csv"), index=False, float_format="%.4g") post_t = pd.DataFrame(post) post_t.to_csv(os.path.join(DATA, "posterior_summary.csv"), index=False, float_format="%.4g") # ---- (e) integrator table: fraction of Milky Way stars at which an isotropic, continuous, narrowband transmitter # would have been above threshold (eta = 1); multiply by eta * duty * beam fraction. Lt = LGRID[LGRID >= 0.999e9] det = pd.DataFrame({"EIRP_W": Lt}) for p in C.PRIORS: cols = {f"__{d}{'_sah' if s else ''}": V[d][SAH if s else BASE]["neff"][p] for d in DRIFTS for s in (False, True)} cols["__modern"] = V["inertial"][MODERN]["neff"][p] cols["__min"], cols["__max"] = envelope(V, p) for tag, arr in cols.items(): det[f"p_{p}{tag}"] = 10 ** np.interp(np.log10(Lt), np.log10(LGRID), np.log10(np.maximum(arr, 1e-300))) / NMW det["N_MW_model"] = NMW det.to_csv(os.path.join(DATA, "detection_probability_table.csv"), index=False, float_format="%.4g") # ---- per-survey contribution at reference EIRPs (baseline build; SETI@home listed but not in the baseline record) contrib = [] for k, c in comps.items(): row = dict(comp_id=k, label=c["label"], group=c["group"], tier=c["tier"], in_baseline_record=(k != "setiathome"), band_MHz=";".join(f"{a:g}-{b:g}" for a, b in c["bands"]), threshold_factor_surface_drift=float(C.drift_loss(P["comps"].set_index("comp_id").loc[k], "surface"))) for L, tag in ((1e12, "1e12"), (L_ARECIBO, "2e13"), (1e15, "1e15"), (1e17, "1e17"), (1e20, "1e20"), (L_KII, "KII")): row[f"N_L{tag}"] = at(c["N_mw"] + c["N_x"], L) for p in ("HI_line", "log_1_10GHz", "log_0p1_100GHz"): kind, f1, f2 = C.PRIORS[p] row[f"freq_weight_{p}"] = sum(C.prior_weight(kind, f1, f2, a, b) for a, b in c["bands"]) contrib.append(row) pd.DataFrame(contrib).to_csv(os.path.join(DATA, "survey_contributions.csv"), index=False, float_format="%.4g") # ---- shares of single components (leave-one-out on the baseline record) def share(drop, p, L, base_inc=C.NO_SAH): ne_all = at(C.neff_all(comps, include=base_inc, n_cap=NMW, priors=[p])[0][p], L) ne_wo = at(C.neff_all(comps, include=lambda k, c: base_inc(k, c) and not drop(k, c), n_cap=NMW, priors=[p])[0][p], L) return 1.0 - ne_wo / ne_all shares = dict( beta_of_KII_1_10GHz=share(lambda k, c: k == "beta", "log_1_10GHz", L_KII), serendip3_of_1e17_0p1_100GHz=share(lambda k, c: k == "serendip3", "log_0p1_100GHz", 1e17), inbeam_of_1e17_1_10GHz=share(lambda k, c: c["group"] == "bl_beam", "log_1_10GHz", 1e17), meta_meta2_of_1e17_HI=share(lambda k, c: k in ("meta", "meta2"), "HI_line", 1e17), setiathome_of_1e17_HI_when_counted=share(lambda k, c: k == "setiathome", "HI_line", 1e17, base_inc=lambda k, c: True), ) bp = C.bycatch_pointing_table() pt = {} for cid, g in bp.groupby("comp_id"): s = g["n_inbeam_detectable"] pt[cid] = dict(n_pointings=int(len(g)), top5_share=float(s.nlargest(5).sum() / s.sum()), top10_share=float(s.nlargest(10).sum() / s.sum()), share_absb_lt5=float(s[np.abs(g["b_deg"]) < 5].sum() / s.sum()), share_absb_lt10=float(s[np.abs(g["b_deg"]) < 10].sum() / s.sum()), median_d_pc=float(g["median_d_pc_of_set"].iloc[0]), top3=[(r.target, round(r.l_deg, 1), round(r.b_deg, 1)) for r in g.head(3).itertuples()]) # ---- continuum-survey caveat for Kardashev-II beacons: flux density of a 3.8e26 W narrowband line averaged # over a continuum band (NVSS: two 42-MHz IFs; 4.85 GHz single-dish surveys: 600 MHz) cont = {} for nm_, bw in (("NVSS_84MHz", 84e6), ("GB6_PMN_600MHz", 600e6), ("VLASS_spw_128MHz", 128e6)): cont[nm_] = {f"{dk:g}kpc": float(L_KII / (FOURPI * (dk * 1000 * PC) ** 2) / bw * 1e26) for dk in (8.2, 25.0)} # ---- comparison with Grimaldi & Marcy (2018): their idealised all-sky 'prior' search (S_min = 1e-23 W/m2 over # 1-2 GHz) gives pi_o = 2.6e-8 (L/L_Arecibo)^1.5 and R = 0.0422 sqrt(L/L_Arecibo) kly (valid for R <~ 1 kly) gm = [] for L in (L_ARECIBO, 1e15, 1e16): x = L / L_ARECIBO for k in ("base", "sah"): ne = at(NE[k]["uniform_1_2GHz"], L) gm.append(dict(EIRP_W=L, record="baseline" if k == "base" else "with SETI@home", GM2018_R_prior_kly=0.0422 * np.sqrt(x), GM2018_pi_prior=2.6e-8 * x ** 1.5, this_work_Neff_1_2GHz=ne, this_work_coverage_fraction=ne / NMW, ratio_this_over_GM=ne / NMW / (2.6e-8 * x ** 1.5), this_work_Nmax95_eta0p5=nm(ne), GM_style_bound_1_over_pi=1.0 / (2.6e-8 * x ** 1.5))) gm = pd.DataFrame(gm) gm.to_csv(os.path.join(DATA, "grimaldi_comparison.csv"), index=False, float_format="%.4g") # ---- (c) Dyson swarms dy = dyson_table(P["shells"]["__allsky__"]) dy.to_csv(os.path.join(DATA, "dyson_limits.csv"), index=False, float_format="%.4g") # ---- key numbers def g(b, p, col, drift=None, sah=False): m = (beh_t.behaviour == b) & (beh_t.freq_prior == p) & (beh_t.setiathome_counted == sah) m &= (beh_t.transmitter_drift == drift) if drift else (beh_t.transmitter_drift == dict((x[0], x[5]) for x in beh)[b]) return float(beh_t[m][col].iloc[0]) def pq(b, p, col, sah=False, pr="loguniform_1_to_Nstars"): m = (post_t.behaviour == b) & (post_t.freq_prior == p) & (post_t.N_prior == pr) & (post_t.setiathome_counted == sah) return float(post_t[m][col].iloc[0]) d300 = dy[dy.T_DS_K == 300.0].set_index("covering_fraction") kn = dict( N_MW_model=NMW, eta=ETA_FID, eta_range=list(ETA_RANGE), stellar_mass_equiv_Msun=NMW / P["stars_per_msun"], Neff_table={f"{L:.3g}": {k: {p: at(NE[k][p], L) for p in C.PRIORS} for k in cfg} for L in keyL}, Ncov_spot={f"{L:.3g}": {name: at(lim[f"Ncov_{name}"].values, L) for name in spot} for L in keyL}, behaviours={b: {p: {f"{d}{'_sah' if s else ''}": dict( p=g(b, p, "p_detect_per_civ", d, s), Nmax95=g(b, p, "Nmax95", d, s), Nmax95_lo=g(b, p, "Nmax95_lo", d, s), Nmax95_hi=g(b, p, "Nmax95_hi", d, s), Neff=g(b, p, "Neff_stars", d, s), mu_all=g(b, p, "expected_detections_if_every_star_hosts_one", d, s)) for d in DRIFTS for s in (False, True)} for p in C.PRIORS} for b in beh_t.behaviour.unique()}, default_drift={x[0]: x[5] for x in beh}, posteriors={b: {p: {("sah" if s else "base"): dict(post_q95=pq(b, p, "post_q95", s), post_median=pq(b, p, "post_median", s), prior_q95=pq(b, p, "prior_q95", s), prior_median=pq(b, p, "prior_median", s), post_q95_lowprior=pq(b, p, "post_q95", s, "loguniform_1e-3_to_Nstars")) for s in (False, True)} for p in C.PRIORS} for b in beh_t.behaviour.unique()}, beacon_overall_range={b: {p: dict(lo=min(g(b, p, "Nmax95_lo", d) for d in ("inertial", "surface")), hi=max(g(b, p, "Nmax95_hi", d) for d in ("inertial", "surface"))) for p in C.PRIORS} for b in ("i_beacon_1e17W", "i_beacon_KII", "i_beacon_1e20W")}, ladder=lad.set_index("scenario").to_dict(orient="index"), variants_1e17={d: {n: {p: nm(at(v["neff"][p], 1e17)) for p in C.PRIORS} for n, v in V[d].items()} for d in ("inertial", "surface")}, variants_KII={d: {n: {p: nm(at(v["neff"][p], L_KII)) for p in C.PRIORS} for n, v in V[d].items()} for d in ("inertial", "surface")}, lte=dict(L_per_channel_W=L_lte_chan, horizon_pc=d_lte_pc, horizon_AU=d_lte_pc * 206265.0, n_fast_stars_within_reach_4GW_narrow=at(comps["fast_nearby"]["N_mw"], L_LTE), dmax_pc_4GW_narrow_FAST=float(dmax_pc(L_LTE, smin_best)), assumed_bandwidth_Hz=LTE_BW), coverage_fraction={k: {p: at(NE[k][p], L_KII) / NMW for p in C.PRIORS} for k in cfg}, coverage_fraction_1e20={k: {p: at(NE[k][p], 1e20) / NMW for p in C.PRIORS} for k in cfg}, haystack_like_fraction_log0p1_100={f"{L:.3g}": at(NE["base"]["log_0p1_100GHz"], L) / NMW for L in keyL}, shares=shares, bycatch_pointings=pt, bycatch_dropped_pointings=[(str(a), float(b), float(c_)) for a, b, c_ in P["cal_diag"]["dropped_pointings"]], bycatch_shell_factors={tel: [float(x) for x in P["cal_diag"]["mason"][tel]["fk"]] for tel in ("GBT", "Parkes")}, bycatch_raw_over_mason={tel: [float(x) for x in P["cal_diag"]["mason"][tel]["raw"] / P["cal_diag"]["mason"][tel]["NB"]] for tel in ("GBT", "Parkes")}, bl_target_rows={k: [v["n_rows_flagged"], v["n_targets"]] for k, v in meta["diag"].items() if "n_targets" in v}, velocity_window_fraction={k: float(P["shells_v"][(v, "__allsky__")].sum() / NMW) for k, v in P["vwin"].items()}, velocity_windows_kms={k: list(v) for k, v in P["vwin"].items()}, drift_threshold_factors={d: {r["comp_id"]: C.drift_loss(r, d) for _, r in P["comps"].iterrows()} for d in ("orbit", "surface")}, drift_coverage_ratio={d: {p: at(NE[d][p], 1e17) / at(NE["base"][p], 1e17) for p in C.PRIORS} for d in ("orbit", "surface")}, continuum_flux_density_Jy_KII=cont, grimaldi=gm.to_dict(orient="records"), dyson_300K={str(k): dict(N_sens_100pc=float(v["N_sensitive_100pc_GCNS"]), N_100pc=float(v["N_stars_100pc_GCNS"]), N_sens_300pc=float(v["N_sensitive_300pc_model"]), N_300pc=float(v["N_total_300pc_model"]), floor_100pc_k0=float(v["floor_frac_100pc_if_100pct_complete_k0"]), floor_300pc_k0=float(v["floor_frac_300pc_if_100pct_complete_k0"]), floor_300pc_k2=(float(v["floor_frac_300pc_if_100pct_complete_k2"]) if np.isfinite(v["floor_frac_300pc_if_100pct_complete_k2"]) else None), published_300pc=float(v["Suazo2022_published_frac_300pc"]), published_100pc=float(v["Suazo2022_published_frac_100pc"]), published_5kpc=float(v["Suazo2022_published_frac_5kpc"]), galaxy_N_published_300pc=float(v["Suazo2022_published_frac_300pc"]) * NMW, galaxy_N_published_5kpc=float(v["Suazo2022_published_frac_5kpc"]) * NMW, galaxy_N_floor_k0=UL95 * NMW / float(v["N_sensitive_300pc_model"])) for k, v in d300.iterrows()}, dyson_candidates=SUAZO24, dyson_k_surviving_gt10=K_SURV_10, ) with open(os.path.join(DATA, "key_numbers.json"), "w") as f: json.dump(kn, f, indent=1) # ---- console summary pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 95) print("\n[analysis] N_MW (model) = %.3e ; eta = %.2f" % (NMW, ETA_FID)) show = beh_t[beh_t.freq_prior.isin(["HI_line", "HI_rest", "log_1_10GHz"]) & beh_t.default_row][ ["behaviour", "freq_prior", "transmitter_drift", "EIRP_W", "Neff_stars", "p_detect_per_civ", "Nmax95", "Nmax95_lo", "Nmax95_hi", "expected_detections_if_every_star_hosts_one"]] print(show.to_string(index=False, float_format=lambda x: f"{x:.3g}")) print(lad[["scenario", "Nmax95_HI_line_1e17W", "Nmax95_HI_rest_1e17W", "Nmax95_log_1_10GHz_1e17W", "Nmax95_HI_line_KII", "Nmax95_HI_rest_KII", "Nmax95_log_1_10GHz_KII"]].to_string(index=False, float_format=lambda x: f"{x:.3g}")) print(rob[rob.prior.isin(["HI_line", "log_1_10GHz"])][["variant", "transmitter_drift", "prior", "in_quoted_range", "Neff_L1e17", "Nmax95_isotropic_1e17W", "Nmax95_isotropic_KII"]].to_string( index=False, float_format=lambda x: f"{x:.3g}")) print("shares:", {k: round(v, 3) for k, v in shares.items()}) print("in-beam pointings:", pt) print("continuum:", cont) print(gm.to_string(index=False, float_format=lambda x: f"{x:.3g}")) print(dy[dy.T_DS_K == 300].T.to_string(float_format=lambda x: f"{x:.4g}")) print("LTE literal: %.3g W per %.2f-Hz channel -> FAST horizon %.4g pc (%.0f AU)" % (L_lte_chan, chan, d_lte_pc, d_lte_pc * 206265)) return V, beh_t, post_t, dy if __name__ == "__main__": main()