"""Interstellar objects: how many are known, the number density that implies for inert bodies, and the limit this places on free-floating artefacts. 1. Count: every object in JPL SBDB with eccentricity > 1.05 (downloaded in s00); those with hyperbolic excess speed > 20 km/s are unambiguously interstellar (the two others are comets with weakly hyperbolic or poorly determined orbits). 2. Effective survey volume: hyperbolic trajectories are injected into the same survey simulator that was calibrated on the near-Earth asteroid discovery history (s02). For speed-at-infinity v, impact parameter b and perihelion time Tp, an interstellar population of number density n delivers n * v * pi * b^2 objects per unit time, so the expected number of discoveries in 1998-2026.8 is N(H) = n * V_eff(H), V_eff = sum_i w_i * P_detect,i(H), w_i = v_i * pi * (b_hi^2 - b_lo^2) * W / N_stratum (W = sampled range of Tp). Inert (non-outgassing) objects only: 2I and 3I were found as active comets. 3. Limits: with k inert interstellar objects found at absolute magnitude <= H (k = 1, 'Oumuamua, for H >= 22.08; k = 0 for brighter), the 95% upper limit on the number density is 3.0/V_eff (k=0) or 4.74/V_eff (k=1). The density implied by the one detection carries the Poisson interval of a single event (95%: 0.025-5.57 times 1/V_eff). 4. Rubin/LSST: ten more years with the 2021-2026 geometry and nightly probability but limiting magnitude 24.0 (scaling exercise, as in s05). Outputs: data/iso_objects.csv, data/iso_detection_volume.csv, data/iso_summary.json """ import json import os import numpy as np import pandas as pd import simlib as S from common import AU_KM, DATA, K_GAUSS, RAW, SEED, D_from_H, H_from_D GM = K_GAUSS ** 2 # au^3 / day^2 KMS = 86400.0 / AU_KM # 1 km/s in au/day PAD_YR = 6.0 # perihelion times sampled this far outside the survey window Q_SPLIT, Q_MAX = 1.5, 10.0 # au; strata in perihelion distance H_GRID = np.arange(12.0, 33.01, 0.5) H_MC_MAX = 24.0 # fainter than this the Monte Carlo has too few close passages: extrapolate AU3_PER_PC3 = 206264.806 ** 3 STARS_PER_PC3 = 0.1 # local stellar number density used by Do et al. (2018) KINEMATICS = { # velocity dispersions (U,V,W) km/s of the population relative to the local standard of rest "disc": (35.0, 25.0, 17.0), # old thin-disc stars "cold": (12.0, 9.0, 7.0), # young, kinematically cold population } V_SUN = np.array([11.1, 12.24, 7.25]) # solar motion relative to the LSR (Schoenrich et al. 2010) def sample_speed(rng, n, model): v = rng.normal(0.0, KINEMATICS[model], (n, 3)) - V_SUN return np.maximum(np.linalg.norm(v, axis=1), 3.0) def hyperbolic_positions(q, v_au_day, Tp_days, inc, Om, w): """Heliocentric positions [n, NDAYS] on hyperbolic orbits.""" a = GM / v_au_day ** 2 # |a| e = 1.0 + q / a nh = np.sqrt(GM / a ** 3) M = nh[:, None] * (S.T_DAYS[None, :] - Tp_days[:, None]) e_ = e[:, None] F = np.where(np.abs(M) < 5 * e_, np.arcsinh(M / e_), np.sign(M) * np.log(2 * np.abs(M) / e_ + 1.8)) for _ in range(12): F = F - (e_ * np.sinh(F) - F - M) / (e_ * np.cosh(F) - 1.0) xp = a[:, None] * (e_ - np.cosh(F)) yp = a[:, None] * np.sqrt(e_ ** 2 - 1.0) * np.sinh(F) cw, sw = np.cos(w)[:, None], np.sin(w)[:, None] cO, sO = np.cos(Om)[:, None], np.sin(Om)[:, None] ci, si = np.cos(inc)[:, None], np.sin(inc)[:, None] x1 = cw * xp - sw * yp y1 = sw * xp + cw * yp x = cO * x1 - sO * ci * y1 y = sO * x1 + cO * ci * y1 z = si * y1 return x, y, z DISC_VARIANTS = ["base", "shallow", "deep", "vdeep", "hard22p5"] def inject(rng, n_per_stratum, model, variants, chunk=200): """Returns ({variant: hist [n, N_ERA, 2, N_G]}, weight [n] in au^3); each variant with its own measured elongation weight.""" W_yr = (S.T1_YEAR - S.T0_YEAR) + 2 * PAD_YR hs, ws = [], [] for (q_lo, q_hi) in [(0.0, Q_SPLIT), (Q_SPLIT, Q_MAX)]: n = n_per_stratum v_kms = sample_speed(rng, n, model) v = v_kms * KMS # au/day b2 = lambda q: q * q + 2 * GM * q / v ** 2 # noqa: E731 (gravitational focusing) b2_lo, b2_hi = b2(q_lo), b2(q_hi) bb = b2_lo + rng.uniform(0, 1, n) * (b2_hi - b2_lo) # invert b^2 = q^2 + 2 GM q / v^2 k = GM / v ** 2 q = -k + np.sqrt(k * k + bb) Tp = (rng.uniform(-PAD_YR, (S.T1_YEAR - S.T0_YEAR) + PAD_YR, n)) * S.DAYS_PER_YEAR inc = np.arccos(rng.uniform(-1, 1, n)) Om = rng.uniform(0, 2 * np.pi, n) w = rng.uniform(0, 2 * np.pi, n) weight = (v * S.DAYS_PER_YEAR) * np.pi * (b2_hi - b2_lo) * W_yr / n # au^3 h = {vv: np.zeros((n, S.N_ERA, 2, S.N_G), dtype=np.float32) for vv in variants} wf = {vv: S.make_wfun(variant=vv) for vv in variants} for s in range(0, n, chunk): sl = slice(s, min(s + chunk, n)) x, y, z = hyperbolic_positions(q[sl], v[sl], Tp[sl], inc[sl], Om[sl], w[sl]) r = np.sqrt(x * x + y * y + z * z) hh = S.hist_from_positions(x, y, z, tmask=(r < 30.0), wfun=wf) for vv in variants: h[vv][sl] = hh[vv] hs.append(h) ws.append(weight) return {vv: np.concatenate([h[vv] for h in hs]) for vv in variants}, np.concatenate(ws) def veff(h, wgt, f, H, lsst=False): vl, vh = np.array(f["vlim"], float), np.array(f["vhard"], float) L = S.survival_log(h, H, vl, f["p"], f["rho"], vh) P = 1 - np.exp(L.sum(axis=1)) out = float((wgt * P).sum()) err = float(np.sqrt(((wgt * P) ** 2).sum())) if not lsst: return out, err vl2, vh2 = vl.copy(), vh.copy() vl2[-1], vh2[-1] = 24.0, 24.7 era_len = S.T1_YEAR - S.ERA_EDGES[-2] L2 = S.survival_log(h, H, vl2, f["p"], f["rho"], vh2)[:, -1] * 10.0 / era_len # volume surveyed by the extra ten years alone (a fresh set of passing objects, statistically # the same as those passing in the last era) P_era_now = 1 - np.exp(L[:, -1]) P_new = 1 - np.exp(L2) return float((wgt * P_new).sum()) , float((wgt * P_era_now).sum()) * 10.0 / era_len def main(n_per_stratum=3000): # ---------------- 1. count ---------------- d = json.load(open(os.path.join(RAW, "sbdb_hyperbolic.json"))) tab = pd.DataFrame(d["data"], columns=d["fields"]) for c in ["e", "q", "i", "H"]: tab[c] = pd.to_numeric(tab[c], errors="coerce") tab["full_name"] = tab.full_name.str.strip() tab["v_inf_kms"] = np.sqrt(GM * (tab.e - 1) / tab.q) / KMS tab["interstellar"] = tab.v_inf_kms > 20.0 tab.to_csv(os.path.join(DATA, "iso_objects.csv"), index=False) pd.set_option("display.width", 220) print(tab[["full_name", "class", "H", "e", "q", "v_inf_kms", "first_obs", "last_obs", "interstellar"]].round(3).to_string(index=False)) n_iso = int(tab.interstellar.sum()) print(f"[iso] objects with e > 1.05 in JPL SBDB: {len(tab)}; unambiguously interstellar (speed at infinity > 20 km/s): {n_iso}") # ---------------- 2. effective survey volume ---------------- fits = json.load(open(os.path.join(DATA, "survey_fit.json"))) if not S.load_elong_table(): # measured elongation weight (s02) raise SystemExit("data/elong_weight_measured.csv missing: run s02 first") rng = np.random.default_rng(SEED + 9) rows = [] summary = dict(n_interstellar_known=n_iso, names=tab[tab.interstellar].full_name.tolist(), au3_per_star=AU3_PER_PC3 / STARS_PER_PC3) for model, n_s in [("disc", n_per_stratum), ("cold", n_per_stratum * 2 // 3)]: vlist = DISC_VARIANTS if model == "disc" else ["base"] h, wgt = inject(rng, n_s, model, vlist) for v in vlist: f = fits[v] for H in H_GRID: ve, err = veff(h[v], wgt, f, H) v_lsst, v_now10 = veff(h[v], wgt, f, H, lsst=True) rows.append(dict(kinematics=model, variant=v, H=H, D_m_p0p2=float(D_from_H(H, 0.2)), D_m_p0p05=float(D_from_H(H, 0.05)), Veff_au3_1998_2026=ve, Veff_mc_err=err, Veff_au3_per_year_2021_2026=v_now10 / 10.0, Veff_au3_per_year_lsst_depth=v_lsst / 10.0)) out = pd.DataFrame(rows) # Fainter than H_MC_MAX only a handful of injected objects pass close enough to Earth, so the # Monte Carlo sum is noisy or zero. The effective volume there is extrapolated with the # log-linear slope fitted over 21 <= H <= 24 (close to the Euclidean -0.6 per magnitude; extra # trailing losses on very close passages would make the true volume smaller still). out["extrapolated"] = out.H > H_MC_MAX slopes = {} for (model, v), g in out.groupby(["kinematics", "variant"]): sel = g[(g.H >= 21.0) & (g.H <= H_MC_MAX)] for col in ["Veff_au3_1998_2026", "Veff_au3_per_year_2021_2026", "Veff_au3_per_year_lsst_depth"]: sl, ic = np.polyfit(sel.H.values, np.log10(sel[col].values), 1) idx = g.index[g.H > H_MC_MAX] out.loc[idx, col] = 10 ** (ic + sl * out.loc[idx, "H"]) if col == "Veff_au3_1998_2026": slopes[f"{model}_{v}"] = float(sl) summary["faint_end_slope_log10Veff_per_mag"] = slopes # ---------------- 3. limits ---------------- out["k_inert_found"] = (out.H >= 22.0).astype(int) # 'Oumuamua, H = 22.08 out["n_UL95_au3"] = np.where(out.k_inert_found == 0, 2.996, 4.744) / out.Veff_au3_1998_2026 out["n_UL95_per_star"] = out.n_UL95_au3 * AU3_PER_PC3 / STARS_PER_PC3 out["n_UL95_after_lsst_au3"] = np.where(out.k_inert_found == 0, 2.996, 4.744) / (out.Veff_au3_1998_2026 + 10 * out.Veff_au3_per_year_lsst_depth) out.to_csv(os.path.join(DATA, "iso_detection_volume.csv"), index=False) base = out[(out.kinematics == "disc") & (out.variant == "base")].set_index("H") print(base[["D_m_p0p2", "Veff_au3_1998_2026", "Veff_mc_err", "Veff_au3_per_year_2021_2026", "Veff_au3_per_year_lsst_depth", "n_UL95_au3", "n_UL95_per_star"]].to_string(float_format=lambda x: f"{x:.3g}")) # 'Oumuamua-like density: one detection at H = 22 for (model, v), g in out.groupby(["kinematics", "variant"]): g = g.set_index("H") ve = g.loc[22.0, "Veff_au3_1998_2026"] summary[f"n_oumuamua_like_au3_{model}_{v}"] = 1.0 / ve summary[f"Veff_H22_au3_{model}_{v}"] = ve summary[f"lsst_gain_rate_H22_{model}_{v}"] = g.loc[22.0, "Veff_au3_per_year_lsst_depth"] / g.loc[22.0, "Veff_au3_per_year_2021_2026"] summary[f"Veff_rate_now_H22_au3_yr_{model}_{v}"] = g.loc[22.0, "Veff_au3_per_year_2021_2026"] summary[f"Veff_rate_now_H12p5_au3_yr_{model}_{v}"] = g.loc[12.5, "Veff_au3_per_year_2021_2026"] # (revision) Poisson uncertainty of a single detection: 95% interval of the mean 0.0253-5.572, 68% 0.173-3.30 ve0 = summary["Veff_H22_au3_disc_base"] summary["n_oumuamua_like_au3_poisson95"] = [0.0253 / ve0, 5.572 / ve0] summary["n_oumuamua_like_au3_poisson68"] = [0.173 / ve0, 3.30 / ve0] vv = [summary[k] for k in summary if k.startswith("Veff_H22_au3_disc_")] summary["Veff_H22_au3_range_over_variants"] = [float(min(vv)), float(max(vv))] summary["n_oumuamua_like_au3_range_over_variants"] = [1.0 / float(max(vv)), 1.0 / float(min(vv))] gg = [summary[k] for k in summary if k.startswith("lsst_gain_rate_H22_disc_")] summary["lsst_gain_rate_H22_range_over_variants"] = [float(min(gg)), float(max(gg))] summary["interstellar_candidates_unconfirmed"] = ("none with an MPC designation; a June 2026 Rubin community-forum post reported a " "candidate linkage of alert detections that the LSST solar-system team judged " "almost certainly spurious (see report)") for Dm in [1.0, 10.0, 100.0, 1000.0, 10000.0]: Hd = float(H_from_D(Dm, 0.2)) ve = float(np.exp(np.interp(Hd, base.index.values, np.log(base.Veff_au3_1998_2026.values)))) k = 1 if Hd >= 22.0 else 0 ul = (4.744 if k else 2.996) / ve summary[f"D{int(Dm)}m_p0.2"] = dict(H=Hd, Veff_au3=ve, n_UL95_au3=ul, n_UL95_per_star=ul * AU3_PER_PC3 / STARS_PER_PC3, extrapolated=bool(Hd > H_MC_MAX)) print(json.dumps(summary, indent=1)) json.dump(summary, open(os.path.join(DATA, "iso_summary.json"), "w"), indent=1) if __name__ == "__main__": main()