"""Compact likelihood of the null technosignature record, for use by an integrator. Stand-alone: needs only numpy/pandas and two small tables written by analysis.py (data/detection_probability_table.csv, data/dyson_limits.csv). Radio model ----------- A civilisation is assumed to sit at a random Milky Way star (civilisations distributed like stars). For a continuous narrowband (<~1 Hz) radio transmitter of isotropic-equivalent power L, the probability that at least one published search would have detected it is p = eta * duty * beam_fraction * C(L | frequency prior, drift class, record) C = N_eff / N_MW : tabulated coverage fraction - the fraction of Milky Way stars at which such a transmitter would have been above the detection threshold of some survey, averaged over the frequency prior and de-duplicated between overlapping surveys eta : probability that an above-threshold signal survives the pipeline / RFI rejection (0.5 fiducial; 0.13-1 in the literature) duty : fraction of time the transmitter is on (single-epoch detection; use duty**2 if a re-detection at a later epoch is required for a signal to count as confirmed) beam_fraction: fraction of the transmitter's sky that its beam covers (1/gain for a randomly pointed beam; 1 if isotropic or deliberately aimed at the Sun; then L is the EIRP in the beam) With N such civilisations the probability of the observed all-null record is exactly (1 - p)^N, so ln L(nulls | N_1..N_k) = sum_k N_k ln(1 - p_k) (the Poisson form -N p is its small-p limit and is NOT used: p reaches eta for powerful beacons). Choices the caller must make (they change the 21-cm answer by up to a factor of ~50) ------------------------------------------------------------------------------- drift 'inertial' (default): the received carrier does not drift in the barycentric frame, i.e. the transmitter is in free fall far from a star or corrects its own Doppler drift towards us. 'orbit' : drifts like a transmitter on an Earth-like orbit (a/c = 2.0e-11 /s). 'surface': drifts like a transmitter on Earth's rotating surface (a/c = 1.1e-10 /s). The 1980s-90s sky surveys (META, META II) integrated 20 s at 0.05 Hz and lose a factor 11 (orbit) to 64 (surface) in sensitivity; modern searches lose a factor 1-7. setiathome_null False (default): SETI@home is NOT counted. As of October 2026 its ~200 final candidates (92 re-observed with FAST) have no published disposition, so the survey is not a completed null. True: treat it as a completed null (all candidates rejected); with drift != 'inertial' its 'non-barycentric' sensitivity column is used. freq_prior 'HI_line' 21-cm line as received in our frame (the transmitter tunes to the Sun's rest frame, the window META was built for) 'HI_rest' 21-cm line in the transmitter's host-star rest frame 'water_hole' 1.42-1.67 GHz uniform; 'log_1_10GHz'; 'log_0p1_100GHz' (log-uniform) bound 'fid' the tabulated fiducial for the chosen drift / record; 'min', 'max' the envelope over all robustness variants and drift classes; 'modern' post-2015 searches only (no SETI@home, no legacy sky surveys), drift-free. Dyson swarms ------------ The only published limits are those of Suazo et al. (2022); dyson_published_limit() returns them. dyson_sensitive_fraction() is a flux-limit sensitivity count, not a limit: it becomes a detection probability only when multiplied by a pipeline completeness that no published search has measured, so p_detect_dyson() has no default completeness. """ import os import numpy as np import pandas as pd _HERE = os.path.dirname(os.path.abspath(__file__)) _DATA = os.path.join(os.path.dirname(_HERE), "data") _TAB = None _DY = None FREQ_PRIORS = ("HI_line", "HI_rest", "water_hole", "log_1_10GHz", "log_0p1_100GHz") DRIFTS = ("inertial", "orbit", "surface") SETIATHOME_COUNTED_BY_DEFAULT = False def _tab(): global _TAB if _TAB is None: _TAB = pd.read_csv(os.path.join(_DATA, "detection_probability_table.csv")) return _TAB def _dy(): global _DY if _DY is None: _DY = pd.read_csv(os.path.join(_DATA, "dyson_limits.csv")) return _DY def coverage_fraction(eirp_w, freq_prior="log_1_10GHz", drift="inertial", setiathome_null=SETIATHOME_COUNTED_BY_DEFAULT, bound="fid"): """C = N_eff(L)/N_MW (efficiency, duty cycle and beaming NOT included). See the module docstring.""" if freq_prior not in FREQ_PRIORS: raise ValueError(f"freq_prior must be one of {FREQ_PRIORS}") if drift not in DRIFTS: raise ValueError(f"drift must be one of {DRIFTS}") t = _tab() if bound == "fid": col = f"p_{freq_prior}__{drift}" + ("_sah" if setiathome_null else "") elif bound in ("min", "max", "modern"): col = f"p_{freq_prior}__{bound}" else: raise ValueError("bound must be 'fid', 'min', 'max' or 'modern'") L = np.atleast_1d(np.asarray(eirp_w, float)) x, y = np.log10(t["EIRP_W"].values), np.log10(np.maximum(t[col].values, 1e-300)) out = 10 ** np.interp(np.log10(np.maximum(L, 1.0)), x, y) out = np.where(L < t["EIRP_W"].values[0], 0.0, out) # below the table: nothing in the record reaches any star out = np.clip(out, 0.0, 1.0) return out if np.ndim(eirp_w) else float(out[0]) def p_detect(eirp_w, freq_prior="log_1_10GHz", duty=1.0, beam_fraction=1.0, eta=0.5, drift="inertial", setiathome_null=SETIATHOME_COUNTED_BY_DEFAULT, bound="fid"): """Probability that the published null record would have detected ONE civilisation of this type.""" return eta * duty * beam_fraction * coverage_fraction(eirp_w, freq_prior, drift, setiathome_null, bound) def log_likelihood(populations): """populations: iterable of dicts {N, eirp_w, freq_prior, duty, beam_fraction, eta, drift, setiathome_null, bound}. Returns the exact ln P(all nulls) = sum N ln(1 - p).""" ll = 0.0 for pop in populations: kw = {k: v for k, v in pop.items() if k != "N"} p = p_detect(**kw) if pop["N"] <= 0: continue if p >= 1.0: return -np.inf ll += pop["N"] * np.log1p(-p) return ll def n_max(eirp_w, cl=0.95, **kw): """Exact binomial upper limit on the number of such civilisations at confidence cl (may exceed the number of stars, in which case the record gives no constraint).""" p = p_detect(eirp_w, **kw) if p <= 0: return np.inf if p >= 1: return 0.0 return np.log(1.0 - cl) / np.log1p(-p) def n_stars_milky_way(): """Number of stellar objects in the Galaxy model used for the table (7.5e10, about 4.2e10 Msun).""" return float(_tab()["N_MW_model"].iloc[0]) def _dyson_row(covering_fraction, t_ds): d = _dy() d = d[np.isclose(d.T_DS_K, t_ds) & np.isclose(d.covering_fraction, covering_fraction)] if len(d) != 1: raise ValueError("tabulated: T = 100/300/600 K, covering fraction = 0.1/0.5/0.9") return d.iloc[0] def dyson_published_limit(covering_fraction_min, volume="300pc"): """Suazo et al. (2022) upper limit on the fraction of stars hosting a ~300 K Dyson swarm with covering fraction above 0.1 / 0.5 / 0.9, within '100pc', '300pc', '1kpc' or '5kpc'. These are the published limits.""" return float(_dyson_row(covering_fraction_min, 300.0)[f"Suazo2022_published_frac_{volume}"]) def dyson_sensitive_fraction(covering_fraction, t_ds=300.0, volume="300pc"): """Fraction of ALL Milky Way stars that lie inside the searched volume ('100pc' Gaia GCNS, '300pc' model) and around which a blackbody swarm (covering fraction 0.1/0.5/0.9, T = 100/300/600 K) would be a >5-sigma WISE W3+W4 source with Gaia G < 19. Sensitivity only: NOT a detection probability.""" r = _dyson_row(covering_fraction, t_ds) col = {"100pc": "N_sensitive_100pc_GCNS", "300pc": "N_sensitive_300pc_model"}[volume] return float(r[col]) / n_stars_milky_way() def p_detect_dyson(covering_fraction, completeness, t_ds=300.0, volume="300pc"): """Probability that a swarm at a random Milky Way star would have been found AND recognised: dyson_sensitive_fraction x completeness. completeness (0-1) is the caller's assumption: no published WISE/Gaia search has an injection-recovery completeness, and the Suazo+ 2024 pipeline kept 7 candidates out of ~5 million stars, of which at most 2 with covering fraction > 0.1 are not yet refuted - so for covering fractions below ~0.2 the record is not a clean null.""" if not 0.0 <= completeness <= 1.0: raise ValueError("completeness must be in [0, 1]") return completeness * dyson_sensitive_fraction(covering_fraction, t_ds, volume) if __name__ == "__main__": print("N_MW (model) = %.3e; SETI@home counted by default: %s" % (n_stars_milky_way(), SETIATHOME_COUNTED_BY_DEFAULT)) ex = [("isotropic 1e17 W beacon, 21-cm (our frame), drift-free", dict(eirp_w=1e17, freq_prior="HI_line")), (" same, SETI@home counted as a null", dict(eirp_w=1e17, freq_prior="HI_line", setiathome_null=True)), (" same, Earth-surface drift, no SETI@home", dict(eirp_w=1e17, freq_prior="HI_line", drift="surface")), (" same, Earth-surface drift, SETI@home counted", dict(eirp_w=1e17, freq_prior="HI_line", drift="surface", setiathome_null=True)), ("isotropic 1e17 W beacon, 21-cm in the star's frame", dict(eirp_w=1e17, freq_prior="HI_rest")), ("isotropic 1e17 W beacon, 1-10 GHz", dict(eirp_w=1e17)), ("isotropic K-II beacon, 21-cm (our frame)", dict(eirp_w=3.83e26, freq_prior="HI_line")), ("isotropic K-II beacon, 1-10 GHz", dict(eirp_w=3.83e26)), (" same, post-2015 searches only", dict(eirp_w=3.83e26, bound="modern")), ("Arecibo-class beam, random pointing, 5% duty", dict(eirp_w=2e13, duty=0.05, beam_fraction=5e-8, drift="surface")), ("Arecibo-class beam aimed at the Sun, 5% duty", dict(eirp_w=2e13, duty=0.05, drift="surface")), ("Earth leakage as a 4 GW narrowband carrier", dict(eirp_w=4e9, drift="surface")), ("Earth leakage, literal (558 W per 3-Hz channel)", dict(eirp_w=558.0, drift="surface"))] for name, kw in ex: print(f" {name:52s} p = {p_detect(**kw):9.3e} N_max(95%) = {n_max(**kw):9.3e}") pops = [dict(N=1000, eirp_w=1e17, freq_prior="HI_line"), dict(N=10, eirp_w=3.83e26)] print(" example ln L(nulls | 1000 x 1e17 W HI beacons + 10 K-II 1-10 GHz beacons) = %.3f" % log_likelihood(pops)) p_kii = p_detect(3.83e26, "HI_line") assert abs(n_max(3.83e26, freq_prior="HI_line") - np.log(0.05) / np.log(1 - p_kii)) < 1e-9 assert abs(log_likelihood([dict(N=3, eirp_w=3.83e26, freq_prior="HI_line")]) - 3 * np.log(1 - p_kii)) < 1e-9 print(" Dyson 300 K: published limits <300 pc (Suazo+ 2022): >10%% %.1e, >50%% %.1e, >90%% %.1e" % tuple(dyson_published_limit(g) for g in (0.1, 0.5, 0.9))) print(" Dyson 300 K, covering 0.5: sensitive fraction of all Galactic stars (300 pc volume) = %.3e (not a limit)" % dyson_sensitive_fraction(0.5))