""" reproduce_letter.py =================== Reproduces every number in the Letter from scratch, with no input other than measured constants. Run it and compare the output with the tables and figures: python reproduce_letter.py # all sections python reproduce_letter.py --section 3 # one section python reproduce_letter.py --table2 # Table 2 (needs the Gaia sample) Nothing here is fitted. The only dimensional input is the interface tension, obtained by inversion from the vacuum energy density and Newton's constant; the local baryonic slope is a measured input, declared, and varied in Sect. 5. The unit of x ------------- The control variable is referred to the FRONT acceleration, x = g_N / (c^2 k), c^2 k = 12 pi G sigma / c^2 = 1.1548e-9 m s^-2, so that the threshold equals unity at the front and falls at r_t = sqrt(A). The same normalisation gives a_0 = (5/48) c^2 k = 1.203e-10 m s^-2 in the galactic sector: one equation, one unit, from stars to galaxies. A second acceleration appears in the Letter and must not be confused with it: 6 pi G sigma / c^2 = 5.77e-10 m s^-2, half the front value, is the acceleration at the CROSSING radius sqrt(2A). It also reads sqrt(3 pi G rho_Lambda) and it coincides with the scale where the published analyses locate the onset of the anomaly — but it does not mark the turnover. Both are printed in Section 1. Sections -------- 1 The tension and the derived scales -> sigma, c^2 k, r_t, a_0 2 The transition function -> gamma(x), the peak, the dip 3 Projection into the observable plane -> Table 1 of the Letter 4 The bump and its projected fine structure -> the peak, the dip, the step 5 Sensitivity to the local baryonic slope -> plateau range 6 Table 2: predictions on a real Gaia sample -> needs the catalogue (see below) Every printed line carries the value the Letter states, as `(Letter: ...)`, so that a disagreement is visible without cross-referencing the text. The Gaia sample for Section 6 ----------------------------- Not redistributed here. Rebuild it in one VizieR query, catalogue J/MNRAS/506/2269 (El-Badry, Rix & Heintz 2021), with the cuts of Chae: BinType = MSMS, R < 0.01, 200 < sepAU < 40000, Plx1 > 5 mas, RPlx1 > 100, RPlx2 > 100, RUWE1 < 1.2, RUWE2 < 1.2 Masses follow the magnitude-mass relation of Banik et al. (2024), inverted on its monotone branch; the inversion is included below so that nothing has to be taken on trust. Requires only numpy and scipy. """ from __future__ import annotations import argparse import numpy as np # -------------------------------------------------------------------------- # Measured constants. These are the only inputs. # -------------------------------------------------------------------------- G_NEWTON = 6.674e-11 # m^3 kg^-1 s^-2 C_LIGHT = 2.99792458e8 # m s^-1 RHO_LAMBDA = 5.3e-10 # J m^-3, cosmological vacuum energy density AU = 1.495979e11 # m M_SUN = 1.98892e30 # kg YEAR = 3.156e7 # s # Local baryonic slope at the position of the pairs, in m s^-2. A measured # input, not a parameter: Sect. 5 varies it over the published range. G_BAR_DEFAULT = 1.80e-10 G_BAR_RANGE = (1.5e-10, 1.9e-10) # MOND, for comparison only: a0 and the Galactic external field in units of a0. A0_MOND = 1.2e-10 G_EXT_OVER_A0 = 1.8 # -------------------------------------------------------------------------- # Section 1 — the tension and the scales derived from it # -------------------------------------------------------------------------- def interface_tension(): """sigma = sqrt(c^4 rho_Lambda / 12 pi G), by inversion of G = c^4 rho_L/(12 pi sigma^2). No acceleration scale enters. This is the single dimensional input. """ return np.sqrt(C_LIGHT**4 * RHO_LAMBDA / (12.0 * np.pi * G_NEWTON)) def front_acceleration(sigma=None): """The front acceleration c^2 k = 12 pi G sigma / c^2. This is the unit of x and the place where the threshold equals unity. The mass cancels against the energy, so it is independent of the source. """ if sigma is None: sigma = interface_tension() return 12.0 * np.pi * G_NEWTON * sigma / C_LIGHT**2 def crossing_acceleration(sigma=None): """The acceleration at the crossing radius sqrt(2A): 6 pi G sigma / c^2. Half the front value. It also reads sqrt(3 pi G rho_Lambda). Reported for comparison with the published onset scale; it does NOT set the threshold. """ return 0.5 * front_acceleration(sigma) def threshold_radius(m_total_kg, sigma=None): """r_t = sqrt(A) = sqrt(E / 12 pi sigma) with E = M c^2. Equivalently sqrt(G M / c^2 k), so r_t grows as sqrt(M). The crossing radius sqrt(2A) is a factor sqrt(2) larger and is a different place. """ if sigma is None: sigma = interface_tension() return np.sqrt(m_total_kg * C_LIGHT**2 / (12.0 * np.pi * sigma)) def section1(verbose=True): sigma = interface_tension() a_front = front_acceleration(sigma) a_cross = crossing_acceleration(sigma) if verbose: print("Section 1 — tension and derived scales") print(f" sigma = {sigma:.4e} J m^-2 (Letter: 4.13e16)") print(f" c^2 k (front) = {a_front:.4e} m s^-2 (Letter: 12 pi G sigma / c^2)") print(f" a_0 = (5/48) c^2 k = {5/48*a_front:.4e} m s^-2 (Letter: 1.203e-10)") print(f" crossing accel. = {a_cross:.4e} m s^-2 (Letter: 5.8e-10, half the front)") print(f" sqrt(3 pi G rho) = {np.sqrt(3*np.pi*G_NEWTON*RHO_LAMBDA):.4e} m s^-2" f" (must equal the line above)") letter_rt = {0.5: 1603, 1.5: 2776, 2.5: 3583} for m in (0.5, 1.5, 2.5): rt = threshold_radius(m * M_SUN, sigma) / AU print(f" r_t({m} Msun) = {rt:8.1f} AU " f"(Letter: {letter_rt[m]})") r_mond = np.sqrt(G_NEWTON * 1.5 * M_SUN / A0_MOND) rt15 = threshold_radius(1.5 * M_SUN, sigma) print(f" r_t / r_MOND = {rt15/r_mond:.4f} (Letter: 0.322)") print(f" crossing radius = {np.sqrt(2)*rt15/AU:8.1f} AU " f"(sqrt(2A) at 1.5 Msun; not the turnover)") return sigma, a_front # -------------------------------------------------------------------------- # Section 2 — the transition function # -------------------------------------------------------------------------- def presence_fraction(x, n_dir=4): """Total presence fraction. Below the threshold (x < 1) equipartition fixes it at (1/2)^n_dir = 1/16. Above, the crossing coordinate becomes evanescent and the presence of each of the n_dir directions rises as f1 = 1 - (1/2) exp(-2 pi sqrt(x-1)); the total is their product. At x = 1, f1 = 1/2 and the product is again 1/16, so the two branches match exactly with no adjustment. """ x = np.asarray(x, float) f1 = 1.0 - 0.5 * np.exp(-2.0 * np.pi * np.sqrt(np.clip(x - 1.0, 0.0, None))) return np.where(x >= 1.0, f1 ** n_dir, (0.5) ** n_dir) def gamma_of_x(x, coeff=4.0 / 3.0): """g_obs / g_N as a function of the control variable x = g_N / (c^2 k). gamma = 1 + (4/3) sqrt(f (1-f)) / sqrt(x). The coefficient 4/3 is a factor 2 from the interference of the two amplitudes times the 2/3 of the shape equation; it has no free parameter. """ x = np.maximum(np.asarray(x, float), 1e-300) f = presence_fraction(x) return 1.0 + coeff * np.sqrt(f * (1.0 - f)) / np.sqrt(x) def x_of_separation(r_m, m_total_kg, g_bar=G_BAR_DEFAULT, sigma=None): """The control variable, composed with the background field. The field equation is linear in the displacement, so the displacements of the pair and of the background add, and so do their slopes: x = x_proper + x_loc. The composition carries no free parameter, and x >= x_loc > 0 by construction. Both terms are referred to the front acceleration c^2 k. """ a_front = front_acceleration(sigma) x_proper = (G_NEWTON * m_total_kg / r_m**2) / a_front x_loc = g_bar / a_front return x_proper + x_loc def gamma_mond(r_m, m_total_kg, nu="simple"): """MOND with an interpolating function and the Galactic external field. For comparison only; not part of the construction. """ g_n = G_NEWTON * m_total_kg / r_m**2 y = (g_n + G_EXT_OVER_A0 * A0_MOND) / A0_MOND if nu == "simple": return 0.5 * (1.0 + np.sqrt(1.0 + 4.0 / y)) return np.sqrt(0.5 * (1.0 + np.sqrt(1.0 + 4.0 / y**2))) def section2(verbose=True): """The peak, the dip and the plateau, all closed-form.""" from scipy.optimize import minimize_scalar res = minimize_scalar(lambda x: -gamma_of_x(x), bounds=(1.0, 1.5), method="bounded", options={"xatol": 1e-10}) x_peak, g_peak = float(res.x), float(-res.fun) g_thr = float(gamma_of_x(1.0)) a_front = front_acceleration() x_loc = G_BAR_DEFAULT / a_front rt15 = threshold_radius(1.5 * M_SUN) if verbose: print("\nSection 2 — the transition function") print(f" below threshold, f = {presence_fraction(0.5):.6f}" f" (must be 1/16 = 0.062500)") print(f" at threshold from above, f = {presence_fraction(1.0):.6f}" f" (must equal the line above)") print(f" gamma at x = 1 - 1e-8 : {float(gamma_of_x(1-1e-8)):.6f}") print(f" gamma at x = 1 + 1e-8 : {float(gamma_of_x(1+1e-8)):.6f}" f" (continuous; the derivative is not)") print(f" peak: x = {x_peak:.4f}, gamma = {g_peak:.4f} (Letter: 1.03, 1.656)") print(f" dip at the threshold: gamma = {g_thr:.4f} (Letter: 1.323)") print(f" height of the dip = {g_peak-g_thr:.4f} (Letter: 0.33)") print(f" peak in separation = {1.07*rt15/AU/1000:.2f} kAU at 1.5 Msun" f" (Letter: 1.07 r_t, 2.97 kAU)") print(f" x_loc (g_bar = {G_BAR_DEFAULT:.2e}) = {x_loc:.4f}") print(f" plateau gamma(x -> x_loc) = {float(gamma_of_x(x_loc)):.4f}" f" (Letter: within 1.80-1.90)") return x_peak, g_peak # -------------------------------------------------------------------------- # Section 3 — projection into the observable plane (Table 1) # -------------------------------------------------------------------------- def project_population(m_total_kg=1.5 * M_SUN, n=600_000, g_bar=G_BAR_DEFAULT, a_min_kau=0.8, a_max_kau=40.0, seed=20260911): """Draw a population and project it, exactly as stated in the Letter. Log-uniform semi-major axes, thermal eccentricities p(e) = 2e, sampled orbital phases through Kepler's equation, isotropic line of sight. This is a first-order smoothing: no selection cuts, no measurement noise, no orbit integration under the modified force. """ rng = np.random.default_rng(seed) a = 10 ** rng.uniform(np.log10(a_min_kau), np.log10(a_max_kau), n) * 1000 * AU e = np.sqrt(rng.random(n)) # p(e) = 2e mean_anom = rng.uniform(0.0, 2.0 * np.pi, n) ecc_anom = mean_anom.copy() for _ in range(60): # Newton-Raphson on Kepler ecc_anom = ecc_anom - (ecc_anom - e * np.sin(ecc_anom) - mean_anom) / ( 1.0 - e * np.cos(ecc_anom)) true_anom = 2.0 * np.arctan2(np.sqrt(1 + e) * np.sin(ecc_anom / 2), np.sqrt(1 - e) * np.cos(ecc_anom / 2)) r = a * (1.0 - e * np.cos(ecc_anom)) incl = np.arccos(rng.uniform(-1.0, 1.0, n)) # isotropic node = rng.uniform(0.0, 2.0 * np.pi, n) phase = node + true_anom s_proj = r * np.sqrt(np.cos(phase) ** 2 + np.cos(incl) ** 2 * np.sin(phase) ** 2) x = x_of_separation(r, m_total_kg, g_bar) return s_proj, r, x BINS_KAU = [(1, 1.5), (2, 2.5), (2.5, 3), (3, 3.5), (3.5, 4), (4, 5), (5, 7), (7, 10), (15, 25)] # Table 1 of the Letter, framework and MOND+EFE rows, in the order of BINS_KAU. LETTER_TABLE1_FRAME = [1.08, 1.30, 1.45, 1.43, 1.47, 1.52, 1.60, 1.68, 1.78] LETTER_TABLE1_MOND = [1.05, 1.10, 1.13, 1.15, 1.18, 1.21, 1.26, 1.31, 1.37] def section3(verbose=True, **kw): """Table 1 of the Letter.""" s_proj, r, x = project_population(**kw) m = kw.get("m_total_kg", 1.5 * M_SUN) g_frame = gamma_of_x(x) g_mond = gamma_mond(r, m) rows = [] for lo, hi in BINS_KAU: sel = (s_proj >= lo * 1000 * AU) & (s_proj < hi * 1000 * AU) rows.append((f"{lo:g}-{hi:g}", int(sel.sum()), float(np.mean(g_frame[sel])), float(np.mean(g_mond[sel])))) if verbose: print("\nSection 3 — Table 1, projected means") print(f" {'s_p (kAU)':>12}{'N':>9}{'framework':>12}{'(Letter)':>10}" f"{'MOND+EFE':>11}{'(Letter)':>10}") for (lab, n_sel, gf, gm), lf, lm in zip(rows, LETTER_TABLE1_FRAME, LETTER_TABLE1_MOND): flag = "" if abs(gf - lf) < 0.005 else " <-- differs" print(f" {lab:>12}{n_sel:>9}{gf:>12.3f}{lf:>10.2f}" f"{gm:>11.3f}{lm:>10.2f}{flag}") return rows # -------------------------------------------------------------------------- # Section 4 — the bump and the fine structure that survives projection # -------------------------------------------------------------------------- def gamma_without_bump(x): """The same chain with the presence held at 1/16 throughout, i.e. the prediction one would make without the rise of the presence above the threshold. The difference with gamma_of_x is the bump.""" x = np.maximum(np.asarray(x, float), 1e-300) f = 1.0 / 16.0 return np.where(x >= 1.0, 1.0, 1.0 + (4.0 / 3.0) * np.sqrt(f * (1 - f)) / np.sqrt(x)) # The Letter's fine-structure claim, in the observable plane, at 1.5 Msun. LETTER_FINE = [(2.75, 2.88, 1.480), (3.00, 3.13, 1.414)] LETTER_FINE_STEP = 0.067 def section4(verbose=True, **kw): s_proj, r, x = project_population(**kw) with_bump = gamma_of_x(x) without = gamma_without_bump(x) rows = [] for lo, hi in [(2, 2.5), (2.5, 3), (3, 3.5), (3.5, 4), (4, 5), (5, 7)]: sel = (s_proj >= lo * 1000 * AU) & (s_proj < hi * 1000 * AU) rows.append((f"{lo:g}-{hi:g}", float(np.mean(with_bump[sel]) - np.mean(without[sel])))) fine = [] for lo, hi, target in LETTER_FINE: sel = (s_proj >= lo * 1000 * AU) & (s_proj < hi * 1000 * AU) fine.append((lo, hi, float(np.mean(with_bump[sel])), target)) if verbose: print("\nSection 4 — the projected shoulder (difference between the two traces)") for lab, d in rows: print(f" {lab:>10} kAU : {d:+.2f}") print("\n the fine structure the Letter quotes, at 1.5 Msun:") for lo, hi, got, target in fine: print(f" {lo:.2f}-{hi:.2f} kAU : mean gamma = {got:.3f}" f" (Letter: {target:.3f})") step = fine[0][2] - fine[1][2] print(f" step between the two = {step:.3f}" f" (Letter: {LETTER_FINE_STEP:.3f})") return rows, fine # -------------------------------------------------------------------------- # Section 5 — sensitivity to the local baryonic slope # -------------------------------------------------------------------------- def section5(verbose=True): a_front = front_acceleration() out = [] if verbose: print("\nSection 5 — plateau against the local baryonic slope") print(" (a measured input, not a parameter; the plateau is not universal)") for g_bar in (G_BAR_RANGE[0], G_BAR_DEFAULT, G_BAR_RANGE[1]): plateau = float(gamma_of_x(g_bar / a_front)) out.append((g_bar, plateau)) if verbose: print(f" g_bar = {g_bar:.2e} m s^-2 -> plateau gamma = {plateau:.3f}") if verbose: print(f" asymptotic plateau: {out[-1][1]:.2f}-{out[0][1]:.2f}" f" (Letter: 1.80-1.90)") print(f" in the observable plane the widest bin reads " f"{LETTER_TABLE1_FRAME[-1]:.2f} at 15-25 kAU (Table 1),") print(" the projection not having reached the asymptote at 25 kAU.") return out # -------------------------------------------------------------------------- # Section 6 — Table 2, predictions on a real Gaia sample # -------------------------------------------------------------------------- def mass_from_absolute_magnitude(m_g): """Invert the magnitude-mass relation of Banik et al. (2024), eq. 6: M_G = 4.887 - 5.693 x + 0.4164 x^2 + 0.9611 x^3, x = ln(M / Msun) The cubic is monotone decreasing only between M ~ 0.21 and ~3.55 Msun; we invert on that branch and extend linearly outside, as they do. """ def mg(x): return 4.887 - 5.693 * x + 0.4164 * x**2 + 0.9611 * x**3 def dmg(x): return -5.693 + 0.8328 * x + 2.8833 * x**2 x_lo, x_hi = np.log(0.21), np.log(2.0) grid_x = np.linspace(x_lo, x_hi, 6000) grid_m = mg(grid_x) m_g = np.asarray(m_g, float) out = np.full_like(m_g, np.nan) inside = (m_g >= grid_m.min()) & (m_g <= grid_m.max()) out[inside] = np.exp(np.interp(m_g[inside], grid_m[::-1], grid_x[::-1])) faint = m_g > grid_m.max() out[faint] = np.exp(x_lo + (m_g[faint] - mg(x_lo)) / dmg(x_lo)) bright = m_g < grid_m.min() out[bright] = np.exp(x_hi + (m_g[bright] - mg(x_hi)) / dmg(x_hi)) return out BINS_TABLE2_AU = [(200, 1000), (1000, 2000), (2000, 3000), (3000, 5000), (5000, 8000), (8000, 15000), (15000, 40000)] def section6(tsv_path, verbose=True, g_bar=G_BAR_DEFAULT): """Table 2: median predicted gamma per separation bin on the Gaia sample. `tsv_path` is the VizieR export described in the module docstring. Required columns: sepAU, Plx1, Gmag1, Gmag2. """ rows_raw, header = [], None with open(tsv_path, encoding="utf-8", errors="replace") as fh: for line in fh: if line.startswith("#") or not line.strip(): continue parts = line.rstrip("\n").split("\t") if header is None: header = [p.strip() for p in parts] continue if parts[0].startswith("---") or not parts[0].strip(): continue if len(parts) == len(header): rows_raw.append(parts) idx = {name: i for i, name in enumerate(header)} def col(name): return np.array([(r[idx[name]].strip() or "nan") for r in rows_raw], float) sep_au, plx = col("sepAU"), col("Plx1") abs_g1 = col("Gmag1") + 5 * np.log10(plx) - 10 abs_g2 = col("Gmag2") + 5 * np.log10(plx) - 10 m_tot = (mass_from_absolute_magnitude(abs_g1) + mass_from_absolute_magnitude(abs_g2)) ok = np.isfinite(m_tot) & (m_tot > 0.2) & (m_tot < 4.0) & np.isfinite(sep_au) m_kg, r_m = m_tot[ok] * M_SUN, sep_au[ok] * AU g_frame = gamma_of_x(x_of_separation(r_m, m_kg, g_bar)) g_mond = gamma_mond(r_m, m_kg) rows = [] for lo, hi in BINS_TABLE2_AU: sel = (sep_au[ok] >= lo) & (sep_au[ok] < hi) if sel.sum() < 20: continue rows.append((f"{lo/1000:g}-{hi/1000:g}", int(sel.sum()), float(np.median(g_frame[sel])), float(np.median(g_mond[sel])))) if verbose: print(f"\nSection 6 — Table 2, on {ok.sum()} pairs from {tsv_path}") print(f" {'kAU':>10}{'pairs':>8}{'framework':>12}{'MOND+EFE':>11}{'diff':>8}") for lab, n_sel, gf, gm in rows: print(f" {lab:>10}{n_sel:>8}{gf:>12.2f}{gm:>11.2f}{gf-gm:>+8.2f}") print(" (Letter: the two plateaux settle near 1.58 and 1.40, gap ~0.18") print(" beyond 8 kAU; inside the threshold the prediction is exactly 1.00)") return rows # -------------------------------------------------------------------------- def main(): ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--section", type=int, choices=[1, 2, 3, 4, 5], help="run one section only") ap.add_argument("--table2", metavar="TSV", help="path to the VizieR export, for Section 6") args = ap.parse_args() if args.table2: section6(args.table2) return todo = [args.section] if args.section else [1, 2, 3, 4, 5] for s in todo: {1: section1, 2: section2, 3: section3, 4: section4, 5: section5}[s]() if __name__ == "__main__": main()