#!/usr/bin/env python3 """Step 4. Short-hop settlement in the real local star field. A. Static percolation of the actual Gaia (GCNS) star systems: for which hop length d_hop is the present-day star field connected? B. Pair-encounter statistics of the empirical Gaia velocity distribution and the kinetic growth rate, as a function of the rendezvous limit (the velocity change a ship can make to stop at a passing star) and of the settleable fraction f. C. Direct simulations in a periodic box populated at the observed density with velocities resampled from the observed 3D velocity distribution: early growth (doubling time) in a cube, front speed in a long slab. LEAD RULE: a 30 km/s ship may settle a star only if the two stars' relative speed is below 30 km/s (it has to stop there). The "any relative speed" rule is kept as an upper bound: it assumes ships that can match the velocity of any passing star (tens to hundreds of km/s). D. Cautious extrapolation to the Galaxy. E. Robustness: direction and slab width of the front, time step, origins that stall, settleable fraction f < 1, a 10 km/s (Voyager-class) variant. Seeds. Doubling times and fronts are measured from seeds placed on a typical disc star (slower than 30 km/s in the mean frame); section E3 starts from random stars instead and counts the origins that stall. Revision after refereeing: the capped rule now leads and covers the full delay grid; stalled seeds are counted instead of averaged away; fronts are measured along U, V and W; f < 1 and the dependence on the rendezvous limit are quantified; the simulator activates launchers within a time step. """ import os import sys import time import numpy as np import pandas as pd from scipy.sparse import coo_matrix from scipy.sparse.csgraph import connected_components from scipy.spatial import cKDTree from scipy.special import lambertw from common import DATA, KMS_TO_PCMYR, SEED, gcns_systems, load_gcns, load_json, save_json from settle_sim import arrival_time, make_box, simulate, time_to_reach D_HOPS = [0.5, 1.0, 2.0, 3.0, 5.0] T_DELAYS = [0.0, 0.01, 0.1, 1.0] # Myr V_SHIP = 30.0 # km/s DV_LIMITS = [30.0, np.inf] # rendezvous limits [km/s]; 30 = lead rule, inf = upper bound DV_TABLE = [5.0, 10.0, 15.0, 20.0, 30.0, 50.0, 100.0, np.inf] F_TABLE = [1.0, 0.3, 0.1, 0.03, 0.01] ETA_C = 0.34189 # continuum-percolation critical filling factor (3D spheres) N_GALAXY = 1e11 R0_KPC, R_OUT_KPC = 8.2, 15.0 AXIS_NAME = {0: "U (radial)", 1: "V (rotation)", 2: "W (vertical)"} AXES = {0: (0, 1, 2), 1: (1, 0, 2), 2: (2, 0, 1)} TRAVEL_FRAC_MIN = 0.4 # a front speed is quoted only above this (see below) QUICK = "--quick" in sys.argv def dv_key(dv): return "any" if np.isinf(dv) else f"dv{dv:.0f}" def dt_limit(d, td, v_ship=V_SHIP): """Upper limit on the time step. In the dense regime (d_hop >= 2 pc) most targets already lie inside d_hop when a launcher activates and are reached after sep/v_ship < d_hop/v_ship; the step must resolve that. In the dilute regime the default step is already converged (section E).""" if d >= 2.0: return 0.5 * d / (v_ship * KMS_TO_PCMYR) + td return None # ---------------------------------------------------------------------------- # A. static percolation # ---------------------------------------------------------------------------- def cluster_stats(xyz, d, sun_index=0): pairs = cKDTree(xyz).query_pairs(d, output_type="ndarray") n = len(xyz) m = coo_matrix((np.ones(len(pairs)), (pairs[:, 0], pairs[:, 1])), shape=(n, n)) _, lab = connected_components(m, directed=False) sizes = np.bincount(lab) sun_lab = lab[sun_index] r = np.linalg.norm(xyz - xyz[sun_index], axis=1) return sizes.max() / n, sizes[sun_lab], r[lab == sun_lab].max() def static_percolation(rng, n_missing_corr): g = gcns_systems(load_gcns()) prim = g[g["is_primary"]] out = {} rows = [] dgrid = np.round(np.arange(0.6, 4.01, 0.1), 2) for R in [25.0, 40.0]: s = prim[prim["dist"] < R] xyz = np.vstack([[0.0, 0.0, 0.0], s[["xcoord50", "ycoord50", "zcoord50"]].values]) sig = 0.5 * (s["Dist84"] - s["Dist16"]).values * 1000 out[f"R{R:.0f}_N"] = int(len(s)) out[f"R{R:.0f}_frac_sigma_d_gt_0p3pc"] = float((sig > 0.3).mean()) # Poisson comparison spheres with the same number of points npois = 10 pois = [] for _ in range(npois): u = rng.normal(size=(len(s), 3)) u /= np.linalg.norm(u, axis=1)[:, None] p = u * (R * rng.random(len(s)) ** (1 / 3))[:, None] pois.append(np.vstack([[0.0, 0.0, 0.0], p])) # Poisson spheres at the completeness-corrected density pois_c = [] nc = int(round(len(s) * n_missing_corr)) for _ in range(npois): u = rng.normal(size=(nc, 3)) u /= np.linalg.norm(u, axis=1)[:, None] p = u * (R * rng.random(nc) ** (1 / 3))[:, None] pois_c.append(np.vstack([[0.0, 0.0, 0.0], p])) for d in dgrid: f, ns, rmax = cluster_stats(xyz, d) fp = [cluster_stats(p, d)[0] for p in pois] fpc = [cluster_stats(p, d)[0] for p in pois_c] rows.append(dict(R_pc=R, d_hop=d, largest_frac=f, sun_cluster=int(ns), sun_reach_pc=rmax, poisson_largest_frac=np.mean(fp), poisson_sd=np.std(fp), poisson_corr_largest_frac=np.mean(fpc), poisson_corr_sd=np.std(fpc))) tab = pd.DataFrame(rows) tab.to_csv(os.path.join(DATA, "static_percolation.csv"), index=False) def d_half(sub, col): x, y = sub["d_hop"].values, sub[col].values k = np.flatnonzero(y >= 0.5) if len(k) == 0 or k[0] == 0: return np.nan k = k[0] return float(x[k - 1] + (0.5 - y[k - 1]) / (y[k] - y[k - 1]) * (x[k] - x[k - 1])) for R in [25.0, 40.0]: sub = tab[tab.R_pc == R] out[f"R{R:.0f}_d50_real"] = d_half(sub, "largest_frac") out[f"R{R:.0f}_d50_poisson"] = d_half(sub, "poisson_largest_frac") out[f"R{R:.0f}_d50_poisson_corrected_density"] = d_half(sub, "poisson_corr_largest_frac") # smallest hop for which the Sun's cluster reaches 90% of the sphere radius ok = sub[sub["sun_reach_pc"] >= 0.9 * R] out[f"R{R:.0f}_d_sun_spans"] = float(ok["d_hop"].min()) if len(ok) else np.nan print(tab[tab.R_pc == 40][["d_hop", "largest_frac", "sun_cluster", "sun_reach_pc", "poisson_largest_frac", "poisson_sd"]].round(3).to_string(index=False)) return out # ---------------------------------------------------------------------------- # B. pair statistics and kinetic growth # ---------------------------------------------------------------------------- def meanfield_rate(gamma, tau): """Growth rate s solving s = gamma * exp(-s tau) (Euler-Lotka).""" if tau <= 0: return gamma return float(np.real(lambertw(gamma * tau)) / tau) def pair_matrix(V, w, rng, ns): """Matrix of relative speeds |v_i - v_j| [km/s] for ns velocities drawn from the weighted sample (float32, built in chunks).""" idx = rng.choice(len(V), size=ns, p=w) v = V[idx].astype(np.float32) D = np.empty((ns, ns), np.float32) for a in range(0, ns, 500): diff = v[a:a + 500, None, :] - v[None, :, :] D[a:a + 500] = np.sqrt((diff * diff).sum(-1)) return D def leading_eigenvalue(D, dv_max, n_iter=400, tol=1e-6): """Leading eigenvalue [km/s] of the kernel |v_i - v_j| 1[|v_i - v_j| < dv_max] averaged over the velocity distribution. In the dilute limit the settled population grows as exp(s t) with s = pi d^2 n * (this speed) for zero lag; it exceeds the plain mean of the kernel because the settled population concentrates where the kernel is largest in velocity space.""" ns = len(D) K = D if np.isinf(dv_max) else np.where(D < dv_max, D, np.float32(0.0)) x = np.ones(ns, np.float32) / ns lam = 0.0 for _ in range(n_iter): y = K @ x / ns new = float(y.sum() / x.sum()) x = y / y.sum() if abs(new - lam) < tol * max(new, 1e-12): lam = new break lam = new return lam def first_hop_kernel(V, w, dv_max): """K_i = sum_j w_j |v_i - v_j| 1[|v_i - v_j| < dv_max] [km/s] for every star of the velocity sample: pi d^2 n K_i is the rate at which origin i meets its first reachable partner.""" v = V.astype(np.float32) w32 = w.astype(np.float32) K = np.empty(len(v), np.float64) for a in range(0, len(v), 400): diff = v[a:a + 400, None, :] - v[None, :, :] D = np.sqrt((diff * diff).sum(-1)) D[D >= dv_max] = 0.0 K[a:a + 400] = D @ w32 return K def wquantile(x, w, q): o = np.argsort(x) cw = np.cumsum(w[o]) / w.sum() return float(x[o][min(np.searchsorted(cw, q), len(x) - 1)]) # ---------------------------------------------------------------------------- # C. simulation helpers # ---------------------------------------------------------------------------- def front_from_run(r): """Arrival times of the half-settled front (mean of the +x and -x sides), the fraction of the time to 1 kpc spent travelling from 250 pc to 1 kpc, the front speed over 250-1000 pc [pc/Myr] (NaN if growth-limited or one-sided), and the same speed measured separately on the +x and -x sides. A side is valid only if it reached 1000 pc before the opposite side reached 1400 pc: after that the opposite front, wrapping around the periodic slab (half-length 1500 pc), fills this side from behind and the arrival time no longer measures this side's own front.""" ta, side = {}, {} for X in (100.0, 250.0, 500.0, 1000.0, 1400.0): tp = arrival_time(r["t"], r["f50p"], X) tm = arrival_time(r["t"], r["f50m"], X) side[X] = (tp, tm) ta[X] = 0.5 * (tp + tm) if np.isfinite(tp) and np.isfinite(tm) else np.nan valid = [np.isfinite(side[1000.0][k]) and not (side[1400.0][1 - k] < side[1000.0][k]) for k in (0, 1)] # a front exists on this scale only if a substantial part of the time to # reach 1 kpc is spent travelling rather than growing in place; otherwise # the slab is growth-limited (well mixed) and no front speed is quoted travel_frac = (ta[1000.0] - ta[250.0]) / ta[1000.0] ok = np.isfinite(travel_frac) and travel_frac >= TRAVEL_FRAC_MIN and all(valid) c_front = 750.0 / (ta[1000.0] - ta[250.0]) if ok else np.nan c_sides = [] for k in (0, 1): tf_k = (side[1000.0][k] - side[250.0][k]) / side[1000.0][k] if valid[k] else np.nan good = valid[k] and tf_k >= TRAVEL_FRAC_MIN c_sides.append(750.0 / (side[1000.0][k] - side[250.0][k]) if good else np.nan) return ta, travel_frac, c_front, c_sides def cube_runs(n_dens, V, w, rng, d, td, dv, n_real, t_end, v_ship=V_SHIP, side=80.0, dtm="auto", seeds="slow", keep=None, tag=None, n_lo=10, n_hi=1000): """Early growth from single seeds in periodic cubes. Returns one record per run: t10 and t1000 (times to n_lo and n_hi settled systems, NaN if not reached by t_end) and T2 = (t1000 - t10)/log2(n_hi/n_lo). seeds = "slow": the star nearest the centre among those moving slower than 30 km/s in the mean frame (a typical disc star), so that the doubling time can be measured in every run; seeds = "random": a random star, i.e. a random origin (used to count the origins that stall under a rendezvous limit).""" recs = [] box = (side, side, side) if isinstance(dtm, str): dtm = dt_limit(d, td, v_ship) for k in range(n_real): pos, vel = make_box(n_dens, box, V, w, rng) if seeds == "slow": slow = np.flatnonzero(np.linalg.norm(vel, axis=1) < 30.0 * KMS_TO_PCMYR) seed = int(slow[np.argmin(np.linalg.norm(pos[slow] - np.array(box) / 2, axis=1))]) else: seed = int(rng.integers(len(pos))) r = simulate(pos, vel, box, d, td, v_ship, seed, t_end=t_end, dv_max_kms=dv, v_cap_kms=min(200.0, dv), stop_n=int(1.2 * n_hi), dt_max=dtm) t10 = time_to_reach(r["t"], r["n"], n_lo) t1000 = time_to_reach(r["t"], r["n"], n_hi) recs.append(dict(t10=t10, t1000=t1000, T2=(t1000 - t10) / np.log2(n_hi / n_lo), seed_speed_kms=float(np.linalg.norm(vel[seed]) / KMS_TO_PCMYR), n_final=int(r["n"][-1]) if len(r["n"]) else 1, dt=r["dt"], t_end=t_end, N=len(pos))) if keep is not None: keep.append(pd.DataFrame(dict(kind="cube", t=r["t"], n=r["n"], run=k, **tag))) return recs def cube_summary(recs): T2 = np.array([r["T2"] for r in recs]) t10 = np.array([r["t10"] for r in recs]) ok = np.isfinite(T2) return dict(n_real=len(recs), n_started=int(np.isfinite(t10).sum()), n_reached=int(ok.sum()), T2_sim_Myr=float(T2[ok].mean()) if ok.any() else np.nan, T2_sim_min=float(T2[ok].min()) if ok.any() else np.nan, T2_sim_max=float(T2[ok].max()) if ok.any() else np.nan, t10_median_Myr=float(np.median(t10[np.isfinite(t10)])) if np.isfinite(t10).any() else np.nan, t1000_sim_Myr=float(np.mean([r["t1000"] for r in recs if np.isfinite(r["t1000"])])) if ok.any() else np.nan) def slab_run(n_dens, V, w, rng, d, td, dv, t_end, v_ship=V_SHIP, width=16.0, axis=0, length=3000.0, dtm="auto", v_cap=None): """Front propagation along the x axis of a long periodic slab, with the Galactic velocity component `axis` laid along x. The seed is the star nearest the centre among those moving slower than 30 km/s in the mean frame (a typical disc star), so that the seed phase does not stall; the front speed is measured between 250 and 1000 pc and does not depend on the seed.""" slab = (length, width, width) pos, vel = make_box(n_dens, slab, V, w, rng, axes=AXES[axis]) slow = np.flatnonzero(np.linalg.norm(vel, axis=1) < 30.0 * KMS_TO_PCMYR) seed = int(slow[np.argmin(np.linalg.norm(pos[slow] - np.array(slab) / 2, axis=1))]) if isinstance(dtm, str): dtm = dt_limit(d, td, v_ship) r = simulate(pos, vel, slab, d, td, v_ship, seed, t_end=t_end, dv_max_kms=dv, v_cap_kms=min(200.0, dv) if v_cap is None else v_cap, stop_front=1060.0, record_every=2, dt_max=dtm) ta, travel_frac, c_front, c_sides = front_from_run(r) vx = np.abs(vel[:, 0]) / KMS_TO_PCMYR info = dict(slab_N=r["N"], width_pc=width, axis=AXIS_NAME[axis], t_arr_100pc=ta[100.0], t_arr_250pc=ta[250.0], t_arr_500pc=ta[500.0], t_arr_1000pc=ta[1000.0], front_pcMyr=c_front, front_kms=c_front / KMS_TO_PCMYR, front_plus_kms=c_sides[0] / KMS_TO_PCMYR, front_minus_kms=c_sides[1] / KMS_TO_PCMYR, mean_speed_to_1kpc_kms=1000.0 / ta[1000.0] / KMS_TO_PCMYR, travel_frac=travel_frac, n_vx_gt_200=int((vx > 200).sum()), vx_max_kms=float(vx.max())) return r, info def main(): rng = np.random.default_rng(SEED + 4) s1 = load_json("s1_density_kinematics.json") n_sys = s1["n_sys_adopted"] out = dict(n_sys=n_sys, v_ship_kms=V_SHIP, settleable_fraction_main_grid=1.0) t_start = time.time() print("== A. static percolation of the real GCNS star field") out["static"] = static_percolation(rng, s1["missing_corr"]) out["static"]["d_c_infinite_poisson_adopted_density"] = float((ETA_C * 6 / np.pi / n_sys) ** (1 / 3)) out["static"]["d_c_infinite_poisson_gcns_density"] = float( (ETA_C * 6 / np.pi / s1["n_systems_gcns_50pc"]) ** (1 / 3)) print({k: (round(v, 3) if isinstance(v, float) else v) for k, v in out["static"].items()}) print("== B. pair statistics") vs = pd.read_csv(os.path.join(DATA, "velocity_sample_50pc.csv")) V = vs[["U", "V", "W"]].values w = vs["weight"].values i = rng.choice(len(V), 2_000_000, p=w) j = rng.choice(len(V), 2_000_000, p=w) vrel = np.linalg.norm(V[i] - V[j], axis=1) out["pair_mean_vrel_kms"] = float(vrel.mean()) out["pair_median_vrel_kms"] = float(np.median(vrel)) gam1 = np.pi * n_sys * vrel.mean() * KMS_TO_PCMYR out["gamma_pair_1pc_perMyr"] = float(gam1) sp = np.linalg.norm(V, axis=1) # speed relative to the Sun cuts = [5, 10, 20, 30, 50, 100, 200] rowsv = [] for c in cuts: rowsv.append(dict(dv_kms=c, pair_frac=float((vrel < c).mean()), pair_flux_frac=float(vrel[vrel < c].sum() / vrel.sum()), sun_frac=float(np.sum(w * (sp < c))), sun_flux_frac=float(np.sum(w * sp * (sp < c)) / np.sum(w * sp)))) dvt = pd.DataFrame(rowsv) dvt.to_csv(os.path.join(DATA, "deltav_flux_fractions.csv"), index=False) print(dvt.round(4).to_string(index=False)) flux_frac = {dv: (1.0 if np.isinf(dv) else float(vrel[vrel < dv].sum() / vrel.sum())) for dv in DV_TABLE} out["pair_flux_frac_dv30"] = flux_frac[30.0] # kinetic (leading-eigenvalue) speed for each rendezvous limit, from two # independent velocity draws of different size (finite-sample check) D_big = pair_matrix(V, w, rng, 2500 if QUICK else 5000) D_small = pair_matrix(V, w, rng, 2500) kin, kin_small = {}, {} for dv in DV_TABLE: kin[dv] = leading_eigenvalue(D_big, dv) kin_small[dv] = leading_eigenvalue(D_small, dv) del D_big, D_small rows = [] for dv in DV_TABLE: vship = min(V_SHIP, dv) # a ship limited to dv also cruises at <= dv row = dict(dv_max_kms=dv, pair_flux_frac=flux_frac[dv], mean_kernel_kms=float(vrel.mean() * flux_frac[dv]), kinetic_speed_kms=kin[dv], kinetic_speed_kms_small_sample=kin_small[dv], v_ship_kms=vship) for d in D_HOPS: g_ = np.pi * n_sys * d**2 * kin[dv] * KMS_TO_PCMYR s_ = meanfield_rate(g_, d / (vship * KMS_TO_PCMYR)) row[f"T2_Myr_d{d:g}"] = float(np.log(2) / s_) row[f"T_grow_1e11_Myr_d{d:g}"] = float(np.log(N_GALAXY) / s_) rows.append(row) kdv = pd.DataFrame(rows) kdv.to_csv(os.path.join(DATA, "kinetic_vs_dvcap.csv"), index=False) print(kdv.round(3).to_string(index=False)) out["kinetic_speed_any_kms"] = kin[np.inf] out["kinetic_speed_dv30_kms"] = kin[30.0] out["kinetic_speed_dv10_kms"] = kin[10.0] out["mean_kernel_any_kms"] = float(vrel.mean()) out["mean_kernel_dv30_kms"] = float(vrel.mean() * flux_frac[30.0]) # settleable fraction f: only a fraction f of systems can be settled (and # then launch). Kinetic growth: s = f * gamma * exp(-s tau). rows = [] for dv in DV_LIMITS: for d in D_HOPS: for td in [0.0, 0.1, 1.0]: for f in F_TABLE: g_ = f * np.pi * n_sys * d**2 * kin[dv] * KMS_TO_PCMYR s_ = meanfield_rate(g_, d / (V_SHIP * KMS_TO_PCMYR) + td) rows.append(dict(dv_max_kms=dv, d_hop_pc=d, t_delay_Myr=td, f=f, s_perMyr=s_, T2_Myr=np.log(2) / s_, T_grow_1e11_Myr=np.log(N_GALAXY) / s_)) pd.DataFrame(rows).to_csv(os.path.join(DATA, "f_scaling_kinetic.csv"), index=False) # origins: how long does a random origin star wait for its first reachable # partner under the 30 km/s limit? K30 = first_hop_kernel(V, w, 30.0) Kany = first_hop_kernel(V, w, np.inf) rate1 = np.pi * n_sys * K30 * KMS_TO_PCMYR # per Myr at d_hop = 1 pc wait1 = 1.0 / np.maximum(rate1, 1e-12) out["origins"] = dict( frac_no_partner_below_30=float(np.sum(w * (K30 == 0))), first_hop_wait_Myr_1pc_median=wquantile(wait1, w, 0.5), first_hop_wait_Myr_1pc_p90=wquantile(wait1, w, 0.9), first_hop_wait_Myr_1pc_p99=wquantile(wait1, w, 0.99), frac_wait_gt_10Myr_1pc=float(np.sum(w * (wait1 > 10.0))), frac_wait_gt_100Myr_1pc=float(np.sum(w * (wait1 > 100.0))), frac_wait_gt_100Myr_2pc=float(np.sum(w * (wait1 / 4.0 > 100.0))), frac_wait_gt_1Gyr_1pc=float(np.sum(w * (wait1 > 1000.0))), first_hop_wait_Myr_1pc_median_any_speed=wquantile( 1.0 / (np.pi * n_sys * Kany * KMS_TO_PCMYR), w, 0.5)) print("origins:", {k: round(v, 4) for k, v in out["origins"].items()}) # ---- validation of the encounter detection ---------------------------- # 40 permanently active seeds, nothing else ever activates, instantaneous # ships: the number of stars "settled" in T must equal the analytic number # of distinct stars that came within d of the seeds. box = (80.0, 80.0, 80.0) pos, vel = make_box(n_sys, box, V, w, rng) seeds = rng.choice(len(pos), 40, replace=False) Tval, dval = 2.0, 1.0 r = simulate(pos, vel, box, dval, 1e9, 1e12, seeds, t_end=Tval) n_contact = int(np.isfinite(r["t_arr"]).sum() - len(seeds)) n_box = len(pos) / np.prod(box) expect = 0.0 for sidx in seeds: rel = np.linalg.norm(vel - vel[sidx], axis=1) # pc/Myr expect += np.pi * dval**2 * n_box * rel.mean() * Tval + n_box * 4 / 3 * np.pi * dval**3 out["validation_contacts_sim"] = n_contact out["validation_contacts_expected"] = float(expect) out["validation_ratio"] = float(n_contact / expect) out["validation_ratio_poisson_err"] = float(np.sqrt(n_contact) / expect) print(f"validation: {n_contact} contacts in the simulator vs {expect:.1f} expected analytically " f"(ratio {n_contact / expect:.3f} +- {np.sqrt(n_contact) / expect:.3f})") # ---- C. simulations --------------------------------------------------- print("== C. simulations (settleable fraction f = 1)") rows = [] curves = [] n_real = 1 if QUICK else 3 for dv in DV_LIMITS: for d in D_HOPS: for td in T_DELAYS: t0 = time.time() gam = gam1 * d**2 * flux_frac[dv] tau = d / (V_SHIP * KMS_TO_PCMYR) + td s_mf = meanfield_rate(gam, tau) s_kin = meanfield_rate(np.pi * n_sys * d**2 * kin[dv] * KMS_TO_PCMYR, tau) eta = n_sys * d**3 # --- cube: early growth from a single seed t_end = 60.0 / max(s_mf, 1e-3) + 5 * tau + 5.0 recs = cube_runs(n_sys, V, w, rng, d, td, dv, n_real, t_end, keep=curves if d in (1.0, 2.0) else None, tag=dict(d_hop=d, t_delay=td, dv=dv)) cs = cube_summary(recs) # --- slab: arrival times of the half-settled front (along U) t_end = 25.0 / max(s_mf, 1e-3) + 1100.0 / 25.0 r, info = slab_run(n_sys, V, w, rng, d, td, dv, t_end) k = slice(None, None, max(1, len(r["t"]) // 300)) curves.append(pd.DataFrame(dict(kind="slab", d_hop=d, t_delay=td, dv=dv, run=0, t=r["t"][k], n=r["n"][k], f50p=r["f50p"][k], f50m=r["f50m"][k]))) rows.append(dict(dv_max_kms=dv, d_hop_pc=d, t_delay_Myr=td, eta=eta, gamma_perMyr=gam, tau_Myr=tau, s_meanfield_perMyr=s_mf, T2_meanfield_Myr=np.log(2) / s_mf, T2_kinetic_Myr=np.log(2) / s_kin, **cs, **info, wall_s=time.time() - t0)) print({k: (round(float(v), 4)) for k, v in rows[-1].items() if k in ("dv_max_kms", "d_hop_pc", "t_delay_Myr", "T2_kinetic_Myr", "T2_sim_Myr", "n_reached", "t_arr_250pc", "t_arr_1000pc", "front_kms", "travel_frac", "wall_s")}, flush=True) sim = pd.DataFrame(rows) print(f" grid done, {time.time() - t_start:.0f} s elapsed", flush=True) # ---- E. robustness and extensions -------------------------------------- rngE = np.random.default_rng(SEED + 5) # independent of the main sequence print("== E1. front speed: direction, slab width, realisation") cases = [ # (rule dv, d_hop, axis, width, v_cap or None, label) (30.0, 3.0, 0, 16.0, None, "second realisation"), (30.0, 3.0, 0, 8.0, None, "8 pc slab"), (30.0, 3.0, 0, 32.0, None, "32 pc slab"), (30.0, 3.0, 1, 16.0, None, "along V"), (30.0, 3.0, 2, 16.0, None, "along W"), (30.0, 2.0, 1, 16.0, None, "along V"), (30.0, 2.0, 2, 16.0, None, "along W"), (30.0, 5.0, 1, 16.0, None, "along V"), (30.0, 5.0, 2, 16.0, None, "along W"), (np.inf, 1.0, 0, 16.0, None, "second realisation"), (np.inf, 1.0, 0, 16.0, 400.0, "v_cap 400 km/s"), (np.inf, 1.0, 0, 8.0, None, "8 pc slab"), (np.inf, 1.0, 0, 32.0, None, "32 pc slab"), (np.inf, 1.0, 1, 16.0, None, "along V"), (np.inf, 1.0, 2, 16.0, None, "along W"), ] if QUICK: cases = [cases[0], cases[4], cases[9]] ft = [] for dv, d, axis, width, vcap, label in cases: t0 = time.time() gam = gam1 * d**2 * flux_frac[dv] s_mf = meanfield_rate(gam, d / (V_SHIP * KMS_TO_PCMYR)) r, info = slab_run(n_sys, V, w, rngE, d, 0.0, dv, 25.0 / s_mf + 1100.0 / 20.0, width=width, axis=axis, v_cap=vcap) ft.append(dict(rule=dv_key(dv), dv_max_kms=dv, d_hop_pc=d, case=label, v_cap_kms=min(200.0, dv) if vcap is None else vcap, **info, wall_s=time.time() - t0)) print({k: (round(v, 2) if isinstance(v, float) else v) for k, v in ft[-1].items() if k in ("rule", "d_hop_pc", "case", "axis", "width_pc", "slab_N", "front_kms", "front_plus_kms", "front_minus_kms", "travel_frac", "wall_s")}, flush=True) if axis != 0 and d in (1.0, 3.0): k = slice(None, None, max(1, len(r["t"]) // 300)) curves.append(pd.DataFrame(dict(kind=f"slab_axis{axis}", d_hop=d, t_delay=0.0, dv=dv, run=0, t=r["t"][k], n=r["n"][k], f50p=r["f50p"][k], f50m=r["f50m"][k]))) ft = pd.DataFrame(ft) ft.to_csv(os.path.join(DATA, "front_tests.csv"), index=False) print("== E2. time-step convergence of the doubling time (same boxes, shorter steps)") tsc = [] conv = [(30.0, 2.0), (30.0, 5.0), (30.0, 1.0), (np.inf, 2.0)] for dv, d in (conv[:1] if QUICK else conv): lag = d / (V_SHIP * KMS_TO_PCMYR) for label, dtm in (("as used", "auto"), ("5 x finer than the lag rule", 0.1 * lag)): if d < 2.0 and label != "as used": dtm, label = 0.5 * lag, "lag rule (finer than used)" rr = np.random.default_rng(SEED + 6) recs = cube_runs(n_sys, V, w, rr, d, 0.0, dv, n_real, 80.0, dtm=dtm) cs = cube_summary(recs) tsc.append(dict(rule=dv_key(dv), d_hop_pc=d, step=label, dt_kyr=1e3 * recs[0]["dt"], T2_sim_kyr=1e3 * cs["T2_sim_Myr"], n_reached=cs["n_reached"])) print(tsc[-1], flush=True) pd.DataFrame(tsc).to_csv(os.path.join(DATA, "timestep_convergence.csv"), index=False) print("== E3. random origins under the 30 km/s limit (stalled seeds)") orig = [] n_or = 4 if QUICK else 24 for d in (1.0, 2.0): recs = cube_runs(n_sys, V, w, rngE, d, 0.0, 30.0, n_or, 100.0, seeds="random") for k, rec in enumerate(recs): orig.append(dict(d_hop_pc=d, run=k, **rec)) t10 = np.array([r_["t10"] for r_ in recs]) t1k = np.array([r_["t1000"] for r_ in recs]) sp_ok = [r_["seed_speed_kms"] for r_ in recs if np.isfinite(r_["t1000"])] sp_no = [r_["seed_speed_kms"] for r_ in recs if not np.isfinite(r_["t1000"])] out[f"origins_sim_d{d:g}"] = dict( n=len(recs), n_first_10_within_100Myr=int(np.isfinite(t10).sum()), n_1000_within_30Myr=int((t1k < 30.0).sum()), n_1000_within_100Myr=int(np.isfinite(t1k).sum()), median_seed_speed_started=float(np.median(sp_ok)) if sp_ok else np.nan, median_seed_speed_stalled=float(np.median(sp_no)) if sp_no else np.nan) print(d, out[f"origins_sim_d{d:g}"], flush=True) pd.DataFrame(orig).to_csv(os.path.join(DATA, "origin_runs.csv"), index=False) print("== E4. settleable fraction f < 1 (30 km/s limit, no delay)") fr = [] fcases = [(2.0, 0.3), (2.0, 0.1), (2.0, 0.03), (1.0, 0.1)] for d, f in (fcases[:1] if QUICK else fcases): t0 = time.time() g_ = f * np.pi * n_sys * d**2 * kin[30.0] * KMS_TO_PCMYR s_kin = meanfield_rate(g_, d / (V_SHIP * KMS_TO_PCMYR)) recs = cube_runs(f * n_sys, V, w, rngE, d, 0.0, 30.0, n_real, 40.0 / s_kin + 20.0, side=80.0 / f ** (1 / 3), dtm=None) cs = cube_summary(recs) fr.append(dict(kind="cube", d_hop_pc=d, f=f, T2_kinetic_Myr=np.log(2) / s_kin, **cs, wall_s=time.time() - t0)) print({k: (round(float(v), 3) if not isinstance(v, str) else v) for k, v in fr[-1].items()}, flush=True) scases = [(3.0, 0.3), (3.0, 0.1), (3.0, 0.03)] for d, f in ([] if QUICK else scases): t0 = time.time() g_ = f * np.pi * n_sys * d**2 * kin[30.0] * KMS_TO_PCMYR s_kin = meanfield_rate(g_, d / (V_SHIP * KMS_TO_PCMYR)) r, info = slab_run(f * n_sys, V, w, rngE, d, 0.0, 30.0, 30.0 / s_kin + 1100.0 / 10.0, width=16.0 / np.sqrt(f), dtm=None) fr.append(dict(kind="slab", d_hop_pc=d, f=f, T2_kinetic_Myr=np.log(2) / s_kin, **info, wall_s=time.time() - t0)) print({k: (round(v, 3) if isinstance(v, float) else v) for k, v in fr[-1].items() if k in ("kind", "d_hop_pc", "f", "T2_kinetic_Myr", "slab_N", "front_kms", "front_plus_kms", "front_minus_kms", "travel_frac", "t_arr_1000pc", "wall_s")}, flush=True) pd.DataFrame(fr).to_csv(os.path.join(DATA, "f_scaling_runs.csv"), index=False) print("== E5. Voyager-class variant: 10 km/s ships, rendezvous limit 10 km/s, no delay") # Only ~1% of stars move within 10 km/s of a given disc star, so an 80 pc # cube holds too few partners (the periodic box cannot refresh its # slow-moving stars); use a 200 pc cube. lo = [] for d, nr in ([] if QUICK else [(5.0, 3), (3.0, 2)]): t0 = time.time() g_ = np.pi * n_sys * d**2 * kin[10.0] * KMS_TO_PCMYR s_kin = meanfield_rate(g_, d / (10.0 * KMS_TO_PCMYR)) recs = cube_runs(n_sys, V, w, rngE, d, 0.0, 10.0, nr, 40.0 / s_kin + 10.0, v_ship=10.0, side=200.0) cs = cube_summary(recs) lo.append(dict(dv_max_kms=10.0, v_ship_kms=10.0, d_hop_pc=d, box_pc=200.0, N=recs[0]["N"], T2_kinetic_Myr=np.log(2) / s_kin, **cs, wall_s=time.time() - t0)) print({k: (round(float(v), 3)) for k, v in lo[-1].items()}, flush=True) pd.DataFrame(lo, columns=["dv_max_kms", "v_ship_kms", "d_hop_pc", "box_pc", "N", "T2_kinetic_Myr", "n_real", "n_started", "n_reached", "T2_sim_Myr", "T2_sim_min", "T2_sim_max", "t10_median_Myr", "t1000_sim_Myr", "wall_s"]).to_csv(os.path.join(DATA, "voyager_class_runs.csv"), index=False) # ---- D. Galaxy extrapolation ----------------------------------------- # growth phase: log2(N_gal) doublings at the simulated local doubling time; # crossing phase: (R0 + R_out) / front speed. Fronts: all f = 1 slab runs # of the rule in which a front exists on the 1 kpc scale, along U, V and W. sim["T_grow_Myr"] = np.log2(N_GALAXY) * sim["T2_sim_Myr"] gal = {} L_PC = (R0_KPC + R_OUT_KPC) * 1000.0 cols = ["d_hop_pc", "axis", "front_kms", "front_plus_kms", "front_minus_kms"] for dv in DV_LIMITS: key = dv_key(dv) grid = sim[(sim.dv_max_kms == dv) & sim.front_kms.notna()][cols] extra = ft[(ft.dv_max_kms == dv) & ft.front_kms.notna() & (ft.v_cap_kms <= 200)][cols] allf = pd.concat([grid, extra], ignore_index=True) gal[f"front_n_runs_{key}"] = int(len(allf)) gal[f"front_min_hop_pc_{key}"] = float(allf.d_hop_pc.min()) for ax_i, ax in AXIS_NAME.items(): sub = allf[allf.axis == ax] tagx = "UVW"[ax_i] gal[f"front_kms_median_{tagx}_{key}"] = float(sub.front_kms.median()) if len(sub) else np.nan gal[f"front_kms_min_{tagx}_{key}"] = float(sub.front_kms.min()) if len(sub) else np.nan gal[f"front_kms_max_{tagx}_{key}"] = float(sub.front_kms.max()) if len(sub) else np.nan gal[f"front_n_runs_{tagx}_{key}"] = int(len(sub)) gal[f"front_kms_min_{key}"] = float(allf.front_kms.min()) gal[f"front_kms_max_{key}"] = float(allf.front_kms.max()) gal[f"T_cross_Myr_fastest_front_{key}"] = float(L_PC / (allf.front_kms.max() * KMS_TO_PCMYR)) gal[f"T_cross_Myr_slowest_front_{key}"] = float(L_PC / (allf.front_kms.min() * KMS_TO_PCMYR)) gal[f"T_cross_Myr_median_U_{key}"] = float(L_PC / (gal[f"front_kms_median_U_{key}"] * KMS_TO_PCMYR)) gal[f"front_side_diff_median_frac_{key}"] = float( (np.abs(allf.front_plus_kms - allf.front_minus_kms) / allf.front_kms).median()) allr = sim[sim.dv_max_kms == dv] gal[f"T_grow_Myr_min_{key}"] = float(allr.T_grow_Myr.min()) gal[f"T_grow_Myr_max_{key}"] = float(allr.T_grow_Myr.max()) out["galaxy"] = gal print(gal) sim.to_csv(os.path.join(DATA, "settlement_runs.csv"), index=False) pd.concat(curves, ignore_index=True).to_csv(os.path.join(DATA, "settlement_curves.csv"), index=False) # kinetic growth-limited fill time versus hop length (no simulation) mf = [] for d in [0.02, 0.05, 0.1, 0.2, 0.5, 1.0, 2.0]: for td in [0.0, 0.1, 1.0, 10.0, 100.0]: for dv in DV_LIMITS: gam = np.pi * n_sys * d**2 * kin[dv] * KMS_TO_PCMYR tau = d / (V_SHIP * KMS_TO_PCMYR) + td s = meanfield_rate(gam, tau) mf.append(dict(d_hop_pc=d, d_hop_au=d * 206264.8, t_delay_Myr=td, dv_max_kms=dv, s_perMyr=s, T2_Myr=np.log(2) / s, T_grow_1e11_Myr=np.log(N_GALAXY) / s)) pd.DataFrame(mf).to_csv(os.path.join(DATA, "meanfield_growth.csv"), index=False) # hop length below which even exponential growth at the local density # cannot produce 1e11 settled systems within 1 Gyr (zero delay, lag # ignored, f = 1); scales as f^-1/2 for dv in DV_LIMITS: s_need = np.log(N_GALAXY) / 1000.0 d_crit = np.sqrt(s_need / (np.pi * n_sys * kin[dv] * KMS_TO_PCMYR)) out[f"d_crit_1Gyr_pc_{dv_key(dv)}"] = float(d_crit) out[f"d_crit_1Gyr_au_{dv_key(dv)}"] = float(d_crit * 206264.8) out["n_runs"] = int(len(sim)) out["wall_s"] = time.time() - t_start save_json("s4_settlement.json", out) print(sim[["dv_max_kms", "d_hop_pc", "t_delay_Myr", "eta", "T2_meanfield_Myr", "T2_kinetic_Myr", "T2_sim_Myr", "T2_sim_min", "T2_sim_max", "n_reached", "t_arr_1000pc", "front_kms", "travel_frac", "T_grow_Myr"]].round(3).to_string(index=False)) return out if __name__ == "__main__": main()