E-012-dbn-heat-flow-toy.py

{
  "experiment": "E-012",
  "task": "3.3: de Bruijn-Newman heat-flow toy (D-021 illustration)",
  "description": "Two-body Rodgers-Tao collision times from E-002 gaps (10^4 zeros) + full all-to-all backward ODE on first 300 zeros; extremal-vs-averaged mismatch",
  "script": "research/numerics/scripts/E-012-dbn-heat-flow-toy.py",
  "parameters": {
    "n_zeros": 10000,
    "n_ode": 300
  },
  "checks": {
    "V1_spread_ge_4_orders": true,
    "V2_ode_vs_two_body_factor4": true,
    "V3_smallest_gap_collides_first": true
  },
  "key_values": {
    "n_gaps": 9999,
    "tc_min": 0.00017764709135458984,
    "tc_median": 0.10824175495920002,
    "tc_max": 5.9293876226367646,
    "spread_orders": 4.523451741355605,
    "binding_pair_index": 6709,
    "binding_gap": 0.03769849772652378,
    "mean_gap": 0.9864634392302998,
    "min_gap": 0.03769849772652378,
    "max_gap": 6.887314497036861,
    "n_ode_zeros": 300,
    "t_collision_ode": -0.02290826416204641,
    "colliding_pair_index": 298,
    "two_body_pair_index": 298,
    "t_collision_two_body": -0.021863900791280276,
    "ratio_ode_over_two_body": 1.0477665619111591,
    "same_pair": true,
    "terminated_on_event": true
  },
  "caveats": [
    "toy conventions: gammas used directly as H_0 zero coordinates; absolute t_c convention-dependent — only the spread and ordering carry the argument",
    "two-body truncation ignores many-body forces; Part 2 ODE quantifies the correction on the binding window"
  ],
  "elapsed_seconds": 0.2,
  "timestamp_utc": "2026-08-18T07:25:35.418404+00:00",
  "row_count": 9999
}

"""E-012 — de Bruijn-Newman heat-flow toy: collision-time statistics.

Illustration for D-021 (task 3.3): the de Bruijn-Newman constant Lambda is
set by the FIRST zero collision under the backward heat flow — an EXTREMAL
statistic (smallest rescaled gap) — whereas the 2/3 certificate is an
AVERAGED (second-moment) statement. This script quantifies the mismatch on
the E-002 table (first 10^4 zeros, dps 25).

Two-body truncation of the Rodgers-Tao zero ODE
    dx_j/dt = 2 * sum_{k != j} 1/(x_j - x_k):
isolated pair with gap delta obeys d(delta^2)/dt = 8, hence backward
collision time t_c = -delta^2 / 8.

Part 1: t_c distribution over all 10^4 - 1 consecutive gaps (spread in
orders of magnitude; binding pair vs mean/median).
Part 2: full all-to-all backward ODE integration on the first N_ODE = 300
zeros (scipy solve_ivp, event at min-gap closure) to check the two-body
estimate against the many-body dynamics.

Checks:
  V1. two-body t_c computed for all gaps; min/max spread >= 4 orders of
      magnitude documented (extremal sensitivity).
  V2. full-ODE first-collision time agrees with two-body estimate for the
      binding pair within a factor of 4 (many-body correction modest).
  V3. Spearman-style sanity: the pair with the smallest gap is the pair
      that collides first in the full ODE.

NOTE (toy conventions): gamma's are used directly as H_0 zero coordinates;
the physical H_t flow uses rescaled variables. Absolute t_c values are
convention-dependent; the ORDER-OF-MAGNITUDE SPREAD and the
extremal-vs-averaged point are not.

Outputs: research/numerics/data/E-012-dbn-heat-flow-toy.csv (+ .json)
"""
import json
import time
from datetime import datetime, timezone
from pathlib import Path

import numpy as np
from scipy.integrate import solve_ivp

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
ZEROS_CSV = DATA / "E-002-zeros-first-10000.csv"
CSV = DATA / "E-012-dbn-heat-flow-toy.csv"
JSON = DATA / "E-012-dbn-heat-flow-toy.json"

N_ODE = 300
GAP_TOL = 1e-10


def load_zeros():
    with ZEROS_CSV.open() as f:
        f.readline()
        return np.array([float(line.split(",")[1]) for line in f])


def force(x):
    """dx_j/dt = 2 * sum_{k != j} 1/(x_j - x_k), pairwise-vectorized."""
    d = x[None, :] - x[:, None]          # d[j, k] = x_k - x_j
    np.fill_diagonal(d, np.inf)
    return -2.0 * np.sum(1.0 / d, axis=1)


def main() -> int:
    if CSV.exists() and "--force" not in __import__("sys").argv:
        print(f"cached: {CSV.exists()} (use --force to recompute)")
        return 0
    t0 = time.time()
    gam = load_zeros()
    gaps = np.diff(gam)

    # ---- Part 1: two-body collision-time distribution
    tc = gaps ** 2 / 8.0                 # backward time |t_c| to collision
    order = np.argsort(tc)
    spread = tc.max() / tc.min()
    stats = {
        "n_gaps": int(len(gaps)),
        "tc_min": float(tc.min()),
        "tc_median": float(np.median(tc)),
        "tc_max": float(tc.max()),
        "spread_orders": float(np.log10(spread)),
        "binding_pair_index": int(order[0]) + 1,   # 1-based zero index j of pair (j, j+1)
        "binding_gap": float(gaps[order[0]]),
        "mean_gap": float(gaps.mean()),
        "min_gap": float(gaps.min()),
        "max_gap": float(gaps.max()),
    }

    # ---- Part 2: full backward ODE on first N_ODE zeros
    x0 = gam[:N_ODE].copy()

    def rhs(t, x):
        return force(x)

    def event(t, x):
        return np.diff(x).min() - GAP_TOL
    event.terminal = True
    event.direction = 1   # integrating backward (t decreasing): gap shrinks

    t_span = (0.0, -10.0 * tc[: N_ODE - 1].min())
    sol = solve_ivp(rhs, t_span, x0, method="DOP853", events=event,
                    rtol=1e-11, atol=1e-13, dense_output=False)
    t_coll_ode = float(sol.t[-1])
    x_end = sol.y[:, -1]
    end_gaps = np.diff(x_end)
    j_coll = int(np.argmin(end_gaps))
    # two-body prediction: binding pair within the ODE window
    j_tb = int(np.argmin(gaps[: N_ODE - 1]))
    t_coll_tb = -float(gaps[j_tb] ** 2 / 8.0)
    ode_stats = {
        "n_ode_zeros": N_ODE,
        "t_collision_ode": t_coll_ode,
        "colliding_pair_index": j_coll + 1,
        "two_body_pair_index": j_tb + 1,
        "t_collision_two_body": t_coll_tb,
        "ratio_ode_over_two_body": t_coll_ode / t_coll_tb if t_coll_tb else None,
        "same_pair": bool(j_coll == j_tb),
        "terminated_on_event": bool(sol.status == 1),
    }

    check_v1 = bool(stats["spread_orders"] >= 4.0)
    check_v2 = bool(ode_stats["terminated_on_event"] and
                    0.25 <= ode_stats["ratio_ode_over_two_body"] <= 4.0)
    check_v3 = bool(ode_stats["same_pair"])
    checks = {
        "V1_spread_ge_4_orders": check_v1,
        "V2_ode_vs_two_body_factor4": check_v2,
        "V3_smallest_gap_collides_first": check_v3,
    }

    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("pair_index_j,gap,two_body_collision_time_abs\n")
        for j in range(len(gaps)):
            f.write(f"{j + 1},{gaps[j]:.17e},{tc[j]:.17e}\n")

    meta = {
        "experiment": "E-012",
        "task": "3.3: de Bruijn-Newman heat-flow toy (D-021 illustration)",
        "description": "Two-body Rodgers-Tao collision times from E-002 "
                       "gaps (10^4 zeros) + full all-to-all backward ODE on "
                       "first 300 zeros; extremal-vs-averaged mismatch",
        "script": "research/numerics/scripts/E-012-dbn-heat-flow-toy.py",
        "parameters": {"n_zeros": len(gam), "n_ode": N_ODE},
        "checks": checks,
        "key_values": {**stats, **ode_stats},
        "caveats": [
            "toy conventions: gammas used directly as H_0 zero "
            "coordinates; absolute t_c convention-dependent — only the "
            "spread and ordering carry the argument",
            "two-body truncation ignores many-body forces; Part 2 ODE "
            "quantifies the correction on the binding window",
        ],
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": int(len(gaps)),
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    ok = all(checks.values())
    print(f"spread: {stats['spread_orders']:.1f} orders; "
          f"t_c min {stats['tc_min']:.3e} vs median {stats['tc_median']:.3e}")
    print(f"ODE first collision at t={t_coll_ode:.6e} (pair {j_coll + 1}); "
          f"two-body {t_coll_tb:.6e} (pair {j_tb + 1}); "
          f"ratio {ode_stats['ratio_ode_over_two_body']:.3f}")
    print(f"wrote {CSV.name}, {JSON.name} — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE: ' + str(checks)}")
    return 0 if ok else 1


if __name__ == "__main__":
    raise SystemExit(main())