#!/usr/bin/env python3 """Step 3. Earth's exposure to passing star systems and the local Fermi constraint. N(d, era) = Gamma_1pc * d^2 * T_era * F(T_era), where F(T) is the ratio of the encounter rate averaged over the last T Myr to the present-day rate. F is estimated by integrating the Sun's orbit in a simple axisymmetric Milky Way potential and averaging an exponential-disc density model along it; the uncertainty range additionally allows for spiral structure, radial migration and the growth and heating of the disc (compare Martinez-Barbosa et al. 2017, MNRAS 464, 2290, whose three solar-orbit scenarios span 0.54-1.6 of their no-migration rate). Revision (referee point 2.5): F is evaluated over the window that is actually used -- the last min(L, T) Myr -- and set to exactly 1 for windows of 10 Myr or less, where the present-day rate applies. """ import os import numpy as np import pandas as pd from common import DATA, KMS_TO_PCMYR, load_json, save_json, tenpc_systems G = 4.30091e-6 # kpc (km/s)^2 / Msun KPC_MYR = KMS_TO_PCMYR * 1e-3 # 1 km/s in kpc/Myr D_LIST = [0.1, 0.5, 1.0, 2.0, 5.0, 10.0] # Eras in which Earth showed a remotely detectable biosphere (durations in Myr, # as (low, central, high)); sources in the report. ERAS = [ ("O2 atmosphere (since GOE)", 2220.0, 2330.0, 2430.0, "wide"), ("High O2/O3 (since NOE)", 540.0, 670.0, 800.0, "narrow"), ("Land vegetation (red edge)", 430.0, 470.0, 510.0, "narrow"), ] # Multiplicative uncertainty of the time-averaged rate relative to today's: F_RANGE = {"narrow": (0.6, 1.3), "wide": (0.5, 2.0)} # ---------------------------------------------------------------------------- # Sun's orbit # ---------------------------------------------------------------------------- def accel(x, y, z, scale): """Acceleration [(km/s)^2/kpc] in a Miyamoto-Nagai disc + Hernquist bulge + NFW halo (structural parameters as in Bovy 2015's MWPotential2014, bulge replaced by a Hernquist sphere); `scale` rescales all masses.""" R2 = x * x + y * y r = np.sqrt(R2 + z * z) # disc Md, a, b = 6.8e10 * scale, 3.0, 0.28 zb = np.sqrt(z * z + b * b) den = (R2 + (a + zb) ** 2) ** 1.5 ax = -G * Md * x / den ay = -G * Md * y / den az = -G * Md * z * (a + zb) / (zb * den) # bulge (Hernquist) Mb, ab = 0.5e10 * scale, 0.6 fb = -G * Mb / (r * (r + ab) ** 2) ax += fb * x; ay += fb * y; az += fb * z # halo (NFW) Mh, ah = 4.37e11 * scale, 16.0 mr = Mh * (np.log(1 + r / ah) - (r / ah) / (1 + r / ah)) fh = -G * mr / r**3 ax += fh * x; ay += fh * y; az += fh * z return ax, ay, az def vcirc(R, scale): ax, _, _ = accel(R, 0.0, 0.0, scale) return np.sqrt(-ax * R) def sun_orbit(R0=8.2, Vc=233.0, uvw=(11.1, 12.24, 7.25), z0=0.0208, t_total=2400.0, dt=0.25): """Integrate the Sun backwards for t_total Myr (leapfrog). Returns R, z [kpc].""" scale = (Vc / vcirc(R0, 1.0)) ** 2 # Galactocentric frame: Sun at (-R0, 0, z0); rotation toward +y; U toward GC (+x) x, y, z = -R0, 0.0, z0 vx, vy, vz = uvw[0], Vc + uvw[1], uvw[2] h = -dt n = int(t_total / dt) Rs = np.empty(n); zs = np.empty(n) ax, ay, az = accel(x, y, z, scale) for i in range(n): vx += 0.5 * h * ax * KPC_MYR; vy += 0.5 * h * ay * KPC_MYR; vz += 0.5 * h * az * KPC_MYR x += h * vx * KPC_MYR; y += h * vy * KPC_MYR; z += h * vz * KPC_MYR ax, ay, az = accel(x, y, z, scale) vx += 0.5 * h * ax * KPC_MYR; vy += 0.5 * h * ay * KPC_MYR; vz += 0.5 * h * az * KPC_MYR Rs[i] = np.hypot(x, y); zs[i] = z return Rs, zs def density_factor(R, z, R0, z0, hR=2.6, hz=0.30, f_thick=0.10, hR_thick=2.0, hz_thick=0.9, profile="sech2"): def vert(zz, h): if profile == "exp": return np.exp(-np.abs(zz) / h) return 1.0 / np.cosh(zz / (2 * h)) ** 2 # ~ 4 exp(-|z|/h) at large |z| def rho(RR, zz): return ((1 - f_thick) * np.exp(-(RR - R0) / hR) * vert(zz, hz) + f_thick * np.exp(-(RR - R0) / hR_thick) * vert(zz, hz_thick)) return float(np.mean(rho(R, z)) / rho(R0, z0)) F_PRESENT_BELOW_MYR = 10.0 # windows this short use the present-day rate (F = 1) def window_factor(Rs, zs, T, dt=0.25, R0=8.2, z0=0.0208, **kw): """F(T): density factor averaged over the last T Myr of the orbit.""" if T <= F_PRESENT_BELOW_MYR: return 1.0 n = max(1, min(len(Rs), int(round(T / dt)))) return density_factor(Rs[:n], zs[:n], R0, z0, **kw) def main(): out = {} s1 = load_json("s1_density_kinematics.json") gam = s1["gamma1pc_adopted"] gam_err = s1["gamma1pc_adopted_err"] n_sys = s1["n_sys_adopted"] # ---- orbit-average factor --------------------------------------------- rows = [] for (R0, Vc) in [(8.2, 233.0), (8.0, 220.0), (8.3, 240.0)]: Rs, zs = sun_orbit(R0=R0, Vc=Vc) for hR in [2.1, 2.6, 3.5]: for prof in ["sech2", "exp"]: f = density_factor(Rs, zs, R0, 0.0208, hR=hR, profile=prof) rows.append(dict(R0=R0, Vc=Vc, hR=hR, profile=prof, R_min=Rs.min(), R_max=Rs.max(), z_max_pc=np.abs(zs).max() * 1000, F_orbit=f)) orb = pd.DataFrame(rows) orb.to_csv(os.path.join(DATA, "orbit_factor.csv"), index=False) print(orb.round(3).to_string(index=False)) base = orb[(orb.R0 == 8.2) & (orb.hR == 2.6) & (orb.profile == "sech2")].iloc[0] out["F_orbit_fiducial"] = float(base["F_orbit"]) out["F_orbit_min"] = float(orb["F_orbit"].min()) out["F_orbit_max"] = float(orb["F_orbit"].max()) out["sun_R_min_kpc"] = float(base["R_min"]); out["sun_R_max_kpc"] = float(base["R_max"]) out["sun_z_max_pc"] = float(base["z_max_pc"]) F0 = out["F_orbit_fiducial"] Rs0, zs0 = sun_orbit(R0=8.2, Vc=233.0) F_of = lambda T: window_factor(Rs0, zs0, T) out["F_window"] = {f"{T:g}": F_of(T) for T in [1.0, 10.0, 30.0, 100.0, 470.0, 670.0, 1000.0, 2330.0]} print("F(T) over the last T Myr:", {k: round(v, 3) for k, v in out["F_window"].items()}) # ---- exposure table --------------------------------------------------- rows = [] for name, tlo, tmid, thi, kind in ERAS: flo, fhi = F_RANGE[kind] Fm = F_of(tmid) for d in D_LIST: n_mid = gam * d**2 * tmid * Fm n_lo = (gam - gam_err) * d**2 * tlo * flo n_hi = (gam + gam_err) * d**2 * thi * fhi rows.append(dict(era=name, T_Myr=tmid, d_pc=d, F_central=Fm, N_central=n_mid, N_low=n_lo, N_high=n_hi, x95_central=-np.log(0.05) / n_mid, x95_conservative=-np.log(0.05) / n_lo)) # radio era: stars currently inside d (the count barely changes in 100 yr) _, pr = tenpc_systems() pr = pr[(pr["rk"] < 3)] dist = 1000.0 / pr["plx"].values for d in D_LIST: n_now = int((dist < d).sum()) if d <= 10 else np.nan n_pass = gam * d**2 * 1e-4 rows.append(dict(era="Radio technology (100 yr)", T_Myr=1e-4, d_pc=d, F_central=1.0, N_central=n_now + n_pass, N_low=n_now, N_high=n_now + 2 * n_pass, x95_central=np.nan, x95_conservative=np.nan)) exp = pd.DataFrame(rows) exp.to_csv(os.path.join(DATA, "exposure_table.csv"), index=False) pd.set_option("display.float_format", lambda v: f"{v:.4g}") print(exp.to_string(index=False)) out["radio_sphere_systems_30p7pc"] = float(n_sys * 4 / 3 * np.pi * 30.66**3) # ---- constraint versus persistence time and delta-v -------------------- # x95 = 3.0 / (Gamma d^2 F(L) L); optionally only slow passers count. # "v<30 km/s" is the self-consistent case for a civilisation whose ships # can change velocity by no more than 30 km/s: to stay in the Solar System # a ship must shed the full relative speed of its home star. dv = pd.read_csv(os.path.join(DATA, "deltav_flux_fractions.csv")) f30 = float(dv.loc[dv.dv_kms == 30, "sun_flux_frac"].iloc[0]) f10 = float(dv.loc[dv.dv_kms == 10, "sun_flux_frac"].iloc[0]) out["sun_flux_frac_v_lt_30"] = f30 out["sun_flux_frac_v_lt_10"] = f10 rows = [] for L in [1.0, 10.0, 100.0, 470.0, 670.0, 2330.0]: FL = F_of(L) for d in D_LIST: for label, ff in (("any speed", 1.0), ("v<30 km/s", f30), ("v<10 km/s", f10)): kind = "wide" if L > 800 else "narrow" flo = 1.0 if L <= F_PRESENT_BELOW_MYR else F_RANGE[kind][0] n = gam * d**2 * L * FL * ff n_cons = (gam - gam_err) * d**2 * L * flo * ff rows.append(dict(persistence_Myr=L, d_pc=d, passers=label, F=FL, N=n, x95=min(1.0, -np.log(0.05) / n), x95_conservative=min(1.0, -np.log(0.05) / n_cons))) con = pd.DataFrame(rows) con.to_csv(os.path.join(DATA, "fermi_constraint.csv"), index=False) for lab in ("any speed", "v<30 km/s", "v<10 km/s"): print(lab) print(con[con.passers == lab].pivot(index="persistence_Myr", columns="d_pc", values="x95").to_string()) print("conservative, any speed") print(con[con.passers == "any speed"].pivot(index="persistence_Myr", columns="d_pc", values="x95_conservative").to_string()) # ---- waiting times and closest expected approaches ---------------------- AU = 206264.8 out["wait_time_Myr"] = {f"{d:g}pc": 1.0 / (gam * d**2) for d in [0.01, 0.05, 0.1, 0.5, 1.0, 2.0]} out["closest_expected_approach_au_1Gyr"] = float(np.sqrt(1.0 / (gam * 1000.0)) * AU) out["closest_expected_approach_au_4p5Gyr_F"] = float(np.sqrt(1.0 / (gam * 4500.0 * F0)) * AU) out["passages_within_1500au_per_Gyr"] = float(gam * (1500.0 / AU) ** 2 * 1000.0) # alternative like-for-like rates for the literature comparison v = s1["mean_speed_adopted_kms"] * KMS_TO_PCMYR out["gamma1pc_all_stellar_objects_10pc"] = float(np.pi * s1["tenpc"]["n_stellar_objects"] / (4 / 3 * np.pi * 1000.0) * v) out["gamma1pc_all_systems_incl_BD_10pc"] = float(np.pi * s1["tenpc"]["dens_systems_all"] * v) out["gamma_slow_passers_v_lt_30_1pc"] = float(gam * f30) out["gamma_slow_passers_v_lt_10_1pc"] = float(gam * f10) key = exp[(exp.era.str.startswith("O2")) & (exp.d_pc == 1.0)].iloc[0] out["N_GOE_1pc_central"] = float(key["N_central"]) out["N_GOE_1pc_low"] = float(key["N_low"]); out["N_GOE_1pc_high"] = float(key["N_high"]) key = exp[(exp.era.str.startswith("Land")) & (exp.d_pc == 1.0)].iloc[0] out["N_land_1pc_central"] = float(key["N_central"]) out["N_land_1pc_low"] = float(key["N_low"]); out["N_land_1pc_high"] = float(key["N_high"]) key = exp[(exp.era.str.startswith("High")) & (exp.d_pc == 1.0)].iloc[0] out["N_NOE_1pc_central"] = float(key["N_central"]) out["N_NOE_1pc_low"] = float(key["N_low"]); out["N_NOE_1pc_high"] = float(key["N_high"]) # Martinez-Barbosa et al. 2017 (abstract): about 21, 39 and 63 encounters # per Myr within 4e5 au for three solar orbits (migrating from the inner # disc, no substantial migration, migrating from the outer disc, in the # order of the abstract); scaled by d^2 to 1 pc (4e5 au = 1.939 pc) mb = np.array([21.0, 39.0, 63.0]) / (4.0e5 / AU) ** 2 out["martinez_barbosa_2017_gamma1pc"] = [float(x) for x in mb] out["martinez_barbosa_2017_ratio_to_no_migration"] = [float(x) for x in mb / mb[1]] save_json("s3_exposure.json", out) return out if __name__ == "__main__": main()