"""Bounded coupled AP1-M1 background pilot with a moving physical split.

This module is a new, hash-bound reference path.  It advances the immutable
128-shell junction state in relative e-fold time, transports the inherited
modes in the validated log-amplitude/Riccati chart, transfers newly resolved
modes from the adiabatic order-0/2/4 tail at q_cut = 0.6 Lambda, and feeds the
renormalized source back into the frozen background equations by a relaxed
Picard iteration.

The pilot is deliberately short and non-production.  It neither reaches the
present anchor nor constructs a two-time response kernel, observable, curve,
fit, or significance.  Only aggregate diagnostics may be serialized; no
trajectory rows or new seed are released.
"""
from __future__ import annotations

import argparse
from dataclasses import asdict, dataclass
from datetime import datetime, timezone
from decimal import Decimal, localcontext
from hashlib import sha256
import json
import math
from pathlib import Path
from typing import Any

import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import CubicSpline

from ap1_m1_background_multiscale_preflight import load_seed
from ap1_m1_log_amplitude_transport_reference import seed_mode_state
from ap1_m1_uv_tail_handoff_reference import direct_resolved_flux, expanded_tail_flux
from ap1_r2c_self_consistent_candidate import classical_terms, standard_stress
from chi_background_closure import (
    ChiParameters,
    ChiTrajectory,
    PV_C,
    PV_J,
    frozen_renormalization_constants,
    portal_terms,
)
from planck2018_neutrino_closure import Planck2018Pilot


EXPECTED_AUTHORITIES = {
    "AP1/APEIRON_AP1_M1_UV_TAIL_HANDOFF_REFERENCE_CHECKPOINT_LATEST.json": (
        "e97b499ddd8ab6c714abc57ccb27a210a5ea6b61bf30a1c184d476291477904c"
    ),
    "AP1/APEIRON_AP1_M1_UV_TAIL_HANDOFF_HELDOUT_LATEST.json": (
        "9d27fc6c19fccc59e63cb76e3f00c67012809ff070554dece3562fcb539fa7bc"
    ),
    "AP1/APEIRON_AP1_M1_BACKGROUND_MULTISCALE_PREFLIGHT_LATEST.json": (
        "c58ee0221e284dc9a5a892cae68b31138b4f3b9c7634f2570880e4f1086d0d60"
    ),
    "AP1/APEIRON_AP1_M1_BACKGROUND_SPEC_LATEST.md": (
        "9a565aef648f4a7229941bf9efa5c51a754fd4763dd11e59e8032701d6ed0259"
    ),
    "AP1/APEIRON_AP1_M1_JUNCTION_128_SEED_CHECKPOINT_LATEST.json": (
        "cf946bfcde4c10a36d5eee99a18fe1a8745483e111304d52f4cc17b3fdf103b4"
    ),
    "AP1/APEIRON_AP1_M1_JUNCTION_128_SEED_LATEST.npz": (
        "644126925fcff32df6f881ee604787751dbe48829072671bf541552bd65dc729"
    ),
    "AP1/CODE/ap1_m1_background_multiscale_preflight.py": (
        "1482b548881af15983b3887434b99713fe0547f3a1e23536e001b3e679b43488"
    ),
    "AP1/CODE/ap1_m1_uv_tail_handoff_reference.py": (
        "567c261cb5938f554984ca0bd695b1f551608d962556e16293fcc02bfa91459a"
    ),
    "AP1/CODE/ap1_m1_log_amplitude_transport_reference.py": (
        "129cf197f08a178eab43c52a2fd9946d161945d16e19ff7f49cdf1d1b12ea523"
    ),
    "AP1/CODE/ap1_r2c_self_consistent_candidate.py": (
        "2c22454f69817ecdc40d38dacad0efaa72fdddcfa50c9c00a48478d610ed8cb0"
    ),
    "AP1/CODE/chi_background_closure.py": (
        "1097888c72fbe9152e1905a2a3ab37c15df9ca681639571b87d76de6ddf072a7"
    ),
    "AP1/CODE/planck2018_neutrino_closure.py": (
        "ccd3158930dcb9ab74d7634ed06f2eaa00745ac95e2883f9fe02d9039a056acc"
    ),
}

SOURCE_KEYS = ("rho", "pressure", "chi2")
WRONSKIAN_LIMIT = 1.0e-11

# New G29 numerical-reference policies.  They are ceilings for deriving new
# bounded-pilot tolerances, not modifications of any pre-existing physical or
# gate threshold.  The held-out plan is embedded in CoupledConfig and therefore
# becomes hash-bound before the first pilot evaluation.
FREEZE_POLICY: dict[str, Any] = {
    "metrics": {
        "max_fixed_point_source_gross_relative_change": {
            "multiplier": 32.0, "floor": 1.0e-9, "ceiling": 2.0e-3,
        },
        "max_background_time_relative_delta": {
            "multiplier": 64.0, "floor": 1.0e-10, "ceiling": 2.0e-4,
        },
        "max_source_time_gross_relative_delta": {
            "multiplier": 64.0, "floor": 1.0e-9, "ceiling": 2.0e-3,
        },
        "max_entry_quadrature_gross_relative_delta": {
            "multiplier": 64.0, "floor": 1.0e-9, "ceiling": 2.0e-3,
        },
        "max_tail_quadrature_gross_relative_delta": {
            "multiplier": 64.0, "floor": 1.0e-9, "ceiling": 2.0e-3,
        },
        "max_uv_cutoff_gross_relative_delta": {
            "multiplier": 64.0, "floor": 1.0e-9, "ceiling": 2.0e-3,
        },
        "max_friedmann_normalized": {
            "multiplier": 32.0, "floor": 1.0e-10, "ceiling": 2.0e-4,
        },
        "max_raychaudhuri_normalized": {
            "multiplier": 32.0, "floor": 1.0e-9, "ceiling": 2.0e-3,
        },
        "ward_normalized": {
            "multiplier": 4.0, "floor": 1.0e-6, "ceiling": 5.0e-2,
        },
        "standard_continuity_normalized": {
            "multiplier": 32.0, "floor": 1.0e-9, "ceiling": 2.0e-3,
        },
    },
    "junction_source_relative_cap": 2.0e-3,
    "junction_raychaudhuri_relative_cap": 1.0e-6,
    "decimal_60_to_80_gross_relative_cap": 1.0e-28,
    "float64_to_decimal80_gross_relative_cap": 1.0e-10,
    "wronskian_relative_cap": WRONSKIAN_LIMIT,
    "handoff_state_relative_cap": 1.0e-12,
    "local_fixed_k_state_relative_cap": 1.0e-8,
    "local_fixed_k_source_gross_relative_cap": 1.0e-8,
    "boundary_direct_expanded_gross_relative_cap": 1.0e-10,
    "last_uv_octave_gross_relative_cap": 2.0e-3,
}


class CoupledBackgroundError(RuntimeError):
    """A chronology, authority, numerical, or physical pilot gate failed."""


@dataclass(frozen=True)
class CoupledConfig:
    qcut_over_Lambda: float = 0.6
    pilot_span_N: float = 0.01
    pilot_fine_nodes: int = 513
    heldout_span_N: float = 0.02
    heldout_fine_nodes: int = 1025
    heldout_evaluation_start_N: float = 0.0125
    adiabatic_support_stride: int = 1
    adiabatic_momentum_support_nodes: int = 257
    fixed_point_steps: int = 6
    relaxation: float = 0.9
    fine_entry_nodes_per_panel: int = 2
    coarse_entry_nodes_per_panel: int = 1
    fine_tail_nodes_per_octave: int = 4
    coarse_tail_nodes_per_octave: int = 2
    fine_K_over_Lambda: float = 64.0
    coarse_K_over_Lambda: float = 32.0
    fine_mode_step_N: float = 1.25e-5
    coarse_mode_step_N: float = 2.5e-5
    fine_background_rtol: float = 2.0e-10
    coarse_background_rtol: float = 2.0e-9
    decimal_low_digits: int = 60
    decimal_high_digits: int = 80

    def validate(self) -> None:
        if self.qcut_over_Lambda != 0.6:
            raise ValueError("the frozen moving split is q_cut/Lambda = 0.6")
        for span, nodes in (
            (self.pilot_span_N, self.pilot_fine_nodes),
            (self.heldout_span_N, self.heldout_fine_nodes),
        ):
            if not 0.0 < span <= 0.025:
                raise ValueError("bounded short background span required")
            if nodes < 17 or nodes % 2 != 1:
                raise ValueError("odd nested fine background grid required")
            if (nodes - 1) & (nodes - 2):
                raise ValueError("2^p+1 fine background nodes required")
        if not self.pilot_span_N < self.heldout_evaluation_start_N < self.heldout_span_N:
            raise ValueError("held-out evaluation must be strictly beyond pilot span")
        if self.fixed_point_steps < 2 or not 0.0 < self.relaxation <= 1.0:
            raise ValueError("bounded relaxed fixed-point iteration required")
        if self.adiabatic_support_stride < 1:
            raise ValueError("separate dense-to-adiabatic support stride required")
        for nodes in (self.pilot_fine_nodes, self.heldout_fine_nodes):
            if (nodes - 1) % self.adiabatic_support_stride:
                raise ValueError("adiabatic support must be nested in fine grid")
        if self.adiabatic_momentum_support_nodes < 65 or self.adiabatic_momentum_support_nodes % 2 != 1:
            raise ValueError("odd dense adiabatic momentum support required")
        if not 1 <= self.coarse_entry_nodes_per_panel < self.fine_entry_nodes_per_panel:
            raise ValueError("nested entry-time quadratures required")
        if not 2 <= self.coarse_tail_nodes_per_octave < self.fine_tail_nodes_per_octave:
            raise ValueError("nested UV-tail quadratures required")
        if not self.qcut_over_Lambda < self.coarse_K_over_Lambda < self.fine_K_over_Lambda:
            raise ValueError("ordered physical UV cutoffs required")
        if not 0.0 < self.fine_mode_step_N < self.coarse_mode_step_N <= 5.0e-5:
            raise ValueError("bounded nested mode steps required")
        ratio = self.coarse_mode_step_N / self.fine_mode_step_N
        if abs(ratio - round(ratio)) > 1.0e-12:
            raise ValueError("mode steps must be nested")
        if not 0.0 < self.fine_background_rtol < self.coarse_background_rtol <= 1.0e-8:
            raise ValueError("ordered bounded background tolerances required")
        if not 32 <= self.decimal_low_digits < self.decimal_high_digits:
            raise ValueError("ordered Decimal precisions required")


def file_sha256(path: Path) -> str:
    digest = sha256()
    with path.open("rb") as stream:
        for block in iter(lambda: stream.read(1024 * 1024), b""):
            digest.update(block)
    return digest.hexdigest()


def _utc_now() -> str:
    return datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")


def _write_pass_json(report: dict[str, Any], output: Path, pass_key: str) -> None:
    if not report.get(pass_key):
        raise CoupledBackgroundError("NONPASS report was not written")
    output.parent.mkdir(parents=True, exist_ok=True)
    output.write_text(json.dumps(report, indent=2) + "\n", encoding="utf-8")


def load_inputs(
    apeiron_root: Path,
) -> tuple[dict[str, Any], dict[str, np.ndarray], dict[str, str]]:
    observed = {
        name: file_sha256(apeiron_root / name) for name in EXPECTED_AUTHORITIES
    }
    if observed != EXPECTED_AUTHORITIES:
        changed = [
            name for name, expected in EXPECTED_AUTHORITIES.items()
            if observed.get(name) != expected
        ]
        raise CoupledBackgroundError(f"authority drift: {changed}")
    checkpoint = json.loads(
        (apeiron_root / "AP1/APEIRON_AP1_M1_UV_TAIL_HANDOFF_REFERENCE_CHECKPOINT_LATEST.json")
        .read_text(encoding="utf-8")
    )
    heldout = json.loads(
        (apeiron_root / "AP1/APEIRON_AP1_M1_UV_TAIL_HANDOFF_HELDOUT_LATEST.json")
        .read_text(encoding="utf-8")
    )
    multiscale = json.loads(
        (apeiron_root / "AP1/APEIRON_AP1_M1_BACKGROUND_MULTISCALE_PREFLIGHT_LATEST.json")
        .read_text(encoding="utf-8")
    )
    if not checkpoint.get("all_checkpoint_gates_pass"):
        raise CoupledBackgroundError("G28 checkpoint is not PASS")
    if not heldout.get("all_heldout_gates_pass"):
        raise CoupledBackgroundError("G28 held-out reference is not PASS")
    if not multiscale.get("all_preflight_gates_pass"):
        raise CoupledBackgroundError("multiscale preflight is not PASS")
    if any(
        report.get("physical_response_kernel_started")
        for report in (checkpoint, heldout, multiscale)
    ):
        raise CoupledBackgroundError("an authority already claims a physical kernel")
    seed_report, arrays, _ = load_seed(apeiron_root)
    if not seed_report.get("all_seed_gates_pass") or not seed_report.get("seed_released"):
        raise CoupledBackgroundError("chronology-safe 128-shell seed is unavailable")
    return seed_report, arrays, observed


def _trajectory_value(tr: ChiTrajectory, name: str, n_value: float) -> float:
    return float(np.interp(n_value, tr.N, np.asarray(getattr(tr, name), dtype=float)))


def _cosmic_derivative(values: np.ndarray, tr: ChiTrajectory) -> np.ndarray:
    array = np.asarray(values, dtype=float)
    derivative = CubicSpline(tr.N, array, axis=0)(tr.N, 1)
    shape = (len(tr.N),) + (1,) * (array.ndim - 1)
    return tr.H.reshape(shape) * derivative


def resample_trajectory(trajectory: ChiTrajectory, nodes: int) -> ChiTrajectory:
    """Nested smooth support for local adiabatic differential functionals."""

    tr = trajectory.validated()
    if nodes < 9 or nodes > len(tr.N):
        raise ValueError("bounded adiabatic support grid required")
    grid = np.linspace(float(tr.N[0]), float(tr.N[-1]), int(nodes))
    return ChiTrajectory(
        N=grid,
        H=CubicSpline(tr.N, tr.H)(grid),
        Hdot=CubicSpline(tr.N, tr.Hdot)(grid),
        sigma=CubicSpline(tr.N, tr.sigma)(grid),
        theta=CubicSpline(tr.N, tr.theta)(grid),
    ).validated()


def adiabatic_primitives_at(
    trajectory: ChiTrajectory,
    sample_N: float,
    q_now: np.ndarray,
    p: ChiParameters | None = None,
) -> dict[str, Any]:
    """Order-0/2/4 state for fixed-k modes crossing at an arbitrary N."""

    p = ChiParameters() if p is None else p
    tr = trajectory.validated()
    q = np.asarray(q_now, dtype=float)
    if q.ndim != 1 or len(q) == 0 or np.any(q <= 0.0) or not np.all(np.isfinite(q)):
        raise ValueError("positive finite one-dimensional physical q required")
    if not float(tr.N[0]) <= sample_N <= float(tr.N[-1]):
        raise ValueError("adiabatic sample must lie on the trajectory interval")
    physical_q = q[None, :] * np.exp(sample_N - tr.N[:, None])
    mass2 = np.asarray(portal_terms(tr.sigma, tr.theta, p)[0], dtype=float)
    muR2 = (p.muR_over_Lambda * p.Lambda) ** 2
    masses2 = np.broadcast_to(
        mass2[:, None, None] + PV_J[None, :, None] * muR2,
        (len(tr.N), 4, len(q)),
    ).copy()
    omega2 = physical_q[:, None, :] ** 2 + masses2
    if np.min(omega2) <= 0.0:
        raise CoupledBackgroundError("moving adiabatic sector entered omega^2 <= 0")
    omega = np.sqrt(omega2)
    omega_dot = _cosmic_derivative(omega, tr)
    omega_ddot = _cosmic_derivative(omega_dot, tr)
    logarithmic0 = omega_dot / omega
    s2 = -2.25 * tr.H[:, None, None] ** 2 - 1.5 * tr.Hdot[:, None, None]
    q2 = 0.75 * logarithmic0**2 - 0.5 * omega_ddot / omega
    W2 = (s2 + q2) / (2.0 * omega)
    W2_dot = _cosmic_derivative(W2, tr)
    W2_ddot = _cosmic_derivative(W2_dot, tr)
    logarithmic2 = W2_dot / omega - logarithmic0 * W2 / omega
    q4 = 1.5 * logarithmic0 * logarithmic2 - 0.5 * (
        W2_ddot / omega - omega_ddot * W2 / omega**2
    )
    W4 = (q4 - W2**2) / (2.0 * omega)
    W = omega + W2 + W4
    W_dot = _cosmic_derivative(W, tr)
    values = (omega, omega_dot, W2, W2_dot, W4, W, W_dot)
    if not all(np.all(np.isfinite(value)) for value in values) or np.min(W) <= 0.0:
        raise CoupledBackgroundError("nonfinite or nonpositive adiabatic state")

    def at(value: np.ndarray) -> np.ndarray:
        return np.asarray(CubicSpline(tr.N, value, axis=0)(sample_N), dtype=float)

    return {
        "q": q.copy(),
        "crossing_N": float(sample_N),
        "H": _trajectory_value(tr, "H", sample_N),
        "mass2": at(masses2),
        "omega": at(omega),
        "omega_dot": at(omega_dot),
        "W2": at(W2),
        "W2_dot": at(W2_dot),
        "W4": at(W4),
        "W": at(W),
        "W_dot": at(W_dot),
    }


def adiabatic_local_table(
    trajectory: ChiTrajectory,
    q_targets: np.ndarray,
    q_support_low: float,
    q_support_high: float,
    q_support_nodes: int,
    p: ChiParameters | None = None,
) -> dict[str, np.ndarray]:
    """Evaluate the fixed-k adiabatic hierarchy through D=H(d_N-d_lnq).

    The returned arrays have shape (time, PV sector, requested physical q).
    Unlike repeatedly constructing one fixed-comoving history per requested
    crossing, this local representation evaluates the same convective cosmic
    derivative on a common (N, ln q) grid.
    """

    p = ChiParameters() if p is None else p
    tr = trajectory.validated()
    targets = np.asarray(q_targets, dtype=float)
    if (
        targets.ndim != 1 or len(targets) == 0
        or np.any(targets <= 0.0) or not np.all(np.isfinite(targets))
    ):
        raise ValueError("positive finite q targets required")
    if not 0.0 < q_support_low < float(np.min(targets)) <= float(np.max(targets)) < q_support_high:
        raise ValueError("q support must strictly bracket every target")
    if q_support_nodes < 65 or q_support_nodes % 2 != 1:
        raise ValueError("odd dense q support required")
    logq = np.linspace(math.log(q_support_low), math.log(q_support_high), q_support_nodes)
    q = np.exp(logq)
    q2 = q[None, None, :] ** 2
    mass2 = np.asarray(portal_terms(tr.sigma, tr.theta, p)[0], dtype=float)
    mass2_N = CubicSpline(tr.N, mass2)(tr.N, 1)
    mass2_dot = tr.H * mass2_N
    mass2_ddot = tr.H * CubicSpline(tr.N, mass2_dot)(tr.N, 1)
    muR2 = (p.muR_over_Lambda * p.Lambda) ** 2
    masses2 = np.broadcast_to(
        mass2[:, None, None] + PV_J[None, :, None] * muR2,
        (len(tr.N), 4, len(q)),
    ).copy()
    omega2 = q2 + masses2
    if np.min(omega2) <= 0.0:
        raise CoupledBackgroundError("local adiabatic q support entered omega^2 <= 0")
    omega = np.sqrt(omega2)
    numerator = mass2_dot[:, None, None] - 2.0 * tr.H[:, None, None] * q2
    numerator_dot = (
        mass2_ddot[:, None, None]
        - 2.0 * tr.Hdot[:, None, None] * q2
        + 4.0 * tr.H[:, None, None] ** 2 * q2
    )
    omega_dot = numerator / (2.0 * omega)
    omega_ddot = numerator_dot / (2.0 * omega) - numerator**2 / (4.0 * omega**3)
    logarithmic0 = omega_dot / omega
    s2 = -2.25 * tr.H[:, None, None] ** 2 - 1.5 * tr.Hdot[:, None, None]
    q2_correction = 0.75 * logarithmic0**2 - 0.5 * omega_ddot / omega
    W2 = (s2 + q2_correction) / (2.0 * omega)

    def convective_cosmic_derivative(value: np.ndarray) -> np.ndarray:
        derivative_N = CubicSpline(tr.N, value, axis=0)(tr.N, 1)
        derivative_logq = CubicSpline(logq, value, axis=2)(logq, 1)
        return tr.H[:, None, None] * (derivative_N - derivative_logq)

    W2_dot = convective_cosmic_derivative(W2)
    W2_ddot = convective_cosmic_derivative(W2_dot)
    logarithmic2 = W2_dot / omega - logarithmic0 * W2 / omega
    q4 = 1.5 * logarithmic0 * logarithmic2 - 0.5 * (
        W2_ddot / omega - omega_ddot * W2 / omega**2
    )
    W4 = (q4 - W2**2) / (2.0 * omega)
    W = omega + W2 + W4
    W_dot = omega_dot + W2_dot + convective_cosmic_derivative(W4)
    fields = {
        "mass2": masses2,
        "omega": omega,
        "omega_dot": omega_dot,
        "W2": W2,
        "W2_dot": W2_dot,
        "W4": W4,
        "W": W,
        "W_dot": W_dot,
    }
    if np.min(W) <= 0.0 or not all(np.all(np.isfinite(value)) for value in fields.values()):
        raise CoupledBackgroundError("local adiabatic table became invalid")
    log_targets = np.log(targets)
    return {
        name: np.asarray(CubicSpline(logq, value, axis=2)(log_targets), dtype=float)
        for name, value in fields.items()
    }


def _log_rhs(
    n_value: float,
    L: np.ndarray,
    x: np.ndarray,
    log_abs_beta: np.ndarray,
    k: np.ndarray,
    trajectory: ChiTrajectory,
    p: ChiParameters,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    H = np.longdouble(_trajectory_value(trajectory, "H", n_value))
    Hdot = np.longdouble(_trajectory_value(trajectory, "Hdot", n_value))
    sigma = _trajectory_value(trajectory, "sigma", n_value)
    theta = _trajectory_value(trajectory, "theta", n_value)
    mass2 = np.longdouble(float(portal_terms(sigma, theta, p)[0]))
    muR2 = np.longdouble((p.muR_over_Lambda * p.Lambda) ** 2)
    omega2 = (
        k[None, :] ** 2 * np.exp(np.longdouble(-2.0 * n_value))
        + mass2
        + PV_J.astype(np.longdouble)[:, None] * muR2
        - 2.25 * H * H
        - 1.5 * Hdot
    )
    if not np.all(np.isfinite(omega2)):
        raise CoupledBackgroundError("nonfinite moving-split mode frequency")
    if float(np.max(log_abs_beta)) > 0.5 * math.log(np.finfo(float).max):
        raise CoupledBackgroundError("Riccati imaginary component overflow")
    beta2 = np.exp(
        np.maximum(
            np.longdouble(2.0) * log_abs_beta,
            np.longdouble(math.log(np.nextafter(0.0, 1.0))),
        )
    )
    dL = x / H
    dx = (-omega2 - x * x + beta2) / H
    dlogb = -2.0 * x / H
    return dL, dx, dlogb


def _advance_log_state(
    state: dict[str, np.ndarray],
    k: np.ndarray,
    left: float,
    right: float,
    trajectory: ChiTrajectory,
    max_step: float,
    p: ChiParameters,
) -> dict[str, np.ndarray]:
    if right < left or not max_step > 0.0:
        raise ValueError("ordered interval and positive step required")
    if right == left:
        return {key: value.copy() for key, value in state.items()}
    steps = max(1, int(math.ceil((right - left) / max_step)))
    h = (right - left) / steps
    L = np.asarray(state["log_amplitude"], dtype=np.longdouble).copy()
    x = np.asarray(state["riccati_real"], dtype=np.longdouble).copy()
    b = np.asarray(state["log_abs_riccati_imag"], dtype=np.longdouble).copy()
    h = np.longdouble(h)
    cursor = np.longdouble(left)
    for _ in range(steps):
        k1 = _log_rhs(cursor, L, x, b, k, trajectory, p)
        k2 = _log_rhs(
            cursor + 0.5 * h,
            L + 0.5 * h * k1[0], x + 0.5 * h * k1[1], b + 0.5 * h * k1[2],
            k, trajectory, p,
        )
        k3 = _log_rhs(
            cursor + 0.5 * h,
            L + 0.5 * h * k2[0], x + 0.5 * h * k2[1], b + 0.5 * h * k2[2],
            k, trajectory, p,
        )
        k4 = _log_rhs(
            cursor + h,
            L + h * k3[0], x + h * k3[1], b + h * k3[2],
            k, trajectory, p,
        )
        L += h * (k1[0] + 2.0 * k2[0] + 2.0 * k3[0] + k4[0]) / 6.0
        x += h * (k1[1] + 2.0 * k2[1] + 2.0 * k3[1] + k4[1]) / 6.0
        b += h * (k1[2] + 2.0 * k2[2] + 2.0 * k3[2] + k4[2]) / 6.0
        cursor += h
    if not all(np.all(np.isfinite(value)) for value in (L, x, b)):
        raise CoupledBackgroundError("log/Riccati transport became nonfinite")
    return {
        "log_amplitude": L,
        "riccati_real": x,
        "log_abs_riccati_imag": b,
    }


def _mode_sector_sources(
    n_value: float,
    H: float,
    sigma: float,
    theta: float,
    k: np.ndarray,
    weights: np.ndarray,
    state: dict[str, np.ndarray],
    p: ChiParameters,
) -> dict[str, np.ndarray]:
    L = np.asarray(state["log_amplitude"], dtype=np.longdouble)
    x = np.asarray(state["riccati_real"], dtype=np.longdouble)
    logb = np.asarray(state["log_abs_riccati_imag"], dtype=np.longdouble)
    n_ld = np.longdouble(n_value)
    H_ld = np.longdouble(H)
    amplitude2 = np.exp(np.longdouble(2.0) * L - np.longdouble(3.0) * n_ld)
    beta2 = np.exp(
        np.maximum(
            np.longdouble(2.0) * logb,
            np.longdouble(math.log(np.nextafter(0.0, 1.0))),
        )
    )
    kinetic = ((x - np.longdouble(1.5) * H_ld) ** 2 + beta2) * amplitude2
    mass2 = np.longdouble(float(portal_terms(sigma, theta, p)[0]))
    muR2 = np.longdouble((p.muR_over_Lambda * p.Lambda) ** 2)
    k_ld = np.asarray(k, dtype=np.longdouble)
    q2 = k_ld[None, :] ** 2 * np.exp(np.longdouble(-2.0) * n_ld)
    masses2 = mass2 + PV_J.astype(np.longdouble)[:, None] * muR2
    each = {
        "rho": np.longdouble(0.5) * (kinetic + (q2 + masses2) * amplitude2),
        "pressure": np.longdouble(0.5) * kinetic - (
            q2 / np.longdouble(6.0) + np.longdouble(0.5) * masses2
        ) * amplitude2,
        "chi2": amplitude2,
    }
    return {
        key: np.asarray(
            np.sum(
                np.asarray(value, dtype=np.longdouble)
                * np.asarray(weights[None, :], dtype=np.longdouble),
                axis=1,
                dtype=np.longdouble,
            ),
            dtype=np.longdouble,
        )
        for key, value in each.items()
    }


def boundary_handoff_table(
    trajectory: ChiTrajectory,
    support_nodes: int,
    qcut: float,
    p: ChiParameters,
    momentum_support_nodes: int = 257,
) -> dict[str, Any]:
    """Precompute the smooth q=q_cut adiabatic state as a function of crossing N."""

    support = resample_trajectory(trajectory, support_nodes)
    local = adiabatic_local_table(
        support, np.array([qcut]), 0.5 * qcut, 2.0 * qcut,
        momentum_support_nodes, p,
    )
    arrays = {
        name: np.asarray(local[name][:, :, 0], dtype=float)
        for name in ("W", "W_dot", "omega", "W2", "W4")
    }
    splines = {name: CubicSpline(support.N, value, axis=0) for name, value in arrays.items()}
    comparison = 0.0
    for index in sorted(set((len(support.N) // 4, len(support.N) // 2, 3 * len(support.N) // 4))):
        independent = adiabatic_primitives_at(
            support, float(support.N[index]), np.array([qcut]), p
        )
        W_reference = np.asarray(independent["W"], dtype=float)[:, 0]
        Wdot_reference = np.asarray(independent["W_dot"], dtype=float)[:, 0]
        u_reference = 1.0 / np.sqrt(2.0 * W_reference)
        v_reference = (
            -Wdot_reference / (2.0 * W_reference) - 1j * W_reference
        ) * u_reference
        W_local = arrays["W"][index]
        Wdot_local = arrays["W_dot"][index]
        u_local = 1.0 / np.sqrt(2.0 * W_local)
        v_local = (-Wdot_local / (2.0 * W_local) - 1j * W_local) * u_local
        comparison = max(
            comparison,
            float(np.max(np.abs(u_local - u_reference) / np.maximum(np.abs(u_reference), 1e-300))),
            float(np.max(np.abs(v_local - v_reference) / np.maximum(np.abs(v_reference), 1e-300))),
        )
    return {
        "support": support,
        "arrays": arrays,
        "splines": splines,
        "max_local_to_independent_fixed_k_state_relative_delta": comparison,
        "momentum_support_nodes": momentum_support_nodes,
    }


def _handoff_state_from_table(
    table: dict[str, Any], entry_N: float
) -> tuple[dict[str, np.ndarray], dict[str, float]]:
    value = {name: np.asarray(spline(entry_N), dtype=float)[:, None] for name, spline in table["splines"].items()}
    value = {name: np.asarray(item, dtype=np.longdouble) for name, item in value.items()}
    W = value["W"]
    Wdot = value["W_dot"]
    if np.min(W) <= 0.0 or not all(np.all(np.isfinite(item)) for item in value.values()):
        raise CoupledBackgroundError("interpolated handoff table became invalid")
    L = np.longdouble(-0.5) * np.log(np.longdouble(2.0) * W)
    x = -Wdot / (np.longdouble(2.0) * W)
    logb = np.log(W)
    tail_u = 1.0 / np.sqrt(2.0 * W)
    tail_v = (x - 1j * W) * tail_u
    reconstructed_u = np.asarray(np.exp(L), dtype=float).astype(complex)
    reconstructed_v = (
        np.asarray(x, dtype=float) - 1j * np.asarray(np.exp(logb), dtype=float)
    ) * reconstructed_u
    observed = -2.0 * np.imag(np.conj(reconstructed_u) * reconstructed_v)
    return (
        {
            "log_amplitude": L,
            "riccati_real": x,
            "log_abs_riccati_imag": logb,
        },
        {
            "min_omega2": float(np.min(value["omega"] ** 2)),
            "min_W": float(np.min(W)),
            "max_abs_W2_over_W0": float(np.max(np.abs(value["W2"] / value["omega"]))),
            "max_abs_W4_over_W0": float(np.max(np.abs(value["W4"] / value["omega"]))),
            "state_transfer_relative_error": max(
                float(np.max(np.abs(reconstructed_u - tail_u) / np.maximum(np.abs(tail_u), 1e-300))),
                float(np.max(np.abs(reconstructed_v - tail_v) / np.maximum(np.abs(tail_v), 1e-300))),
            ),
            "initial_wronskian_relative_error": float(np.max(np.abs(observed - 1.0))),
        },
    )


def transport_resolved_modes(
    trajectory: ChiTrajectory,
    arrays: dict[str, np.ndarray],
    entry_nodes_per_panel: int,
    max_step: float,
    qcut_over_Lambda: float = 0.6,
    adiabatic_support_nodes: int | None = None,
    adiabatic_momentum_support_nodes: int = 257,
) -> dict[str, Any]:
    """Transport inherited modes and all tail-to-resolved crossings."""

    tr = trajectory.validated()
    p = ChiParameters()
    qcut = qcut_over_Lambda * p.Lambda
    support_nodes = (
        min(len(tr.N), 65) if adiabatic_support_nodes is None
        else int(adiabatic_support_nodes)
    )
    handoff_table = boundary_handoff_table(
        tr, support_nodes, qcut, p, adiabatic_momentum_support_nodes
    )
    seed_indices = tuple(range(len(np.asarray(arrays["k_comoving_Mpl"]))))
    seed = seed_mode_state(arrays, seed_indices)
    state = {
        "log_amplitude": np.asarray(seed["log_amplitude"], dtype=np.longdouble).copy(),
        "riccati_real": np.asarray(seed["riccati_real"], dtype=np.longdouble).copy(),
        "log_abs_riccati_imag": np.asarray(seed["log_abs_riccati_imag"], dtype=np.longdouble).copy(),
    }
    log_wronskian_target = np.log(np.asarray(seed["wronskian"], dtype=np.longdouble))
    k = np.asarray(arrays["k_comoving_Mpl"], dtype=np.longdouble).copy()
    weights = np.asarray(arrays["shell_weights_Mpl3"], dtype=np.longdouble).copy()
    seed_count = len(k)
    sector_total = {key: np.zeros((len(tr.N), 4), dtype=np.longdouble) for key in SOURCE_KEYS}
    sector_seed = {key: np.zeros((len(tr.N), 4), dtype=np.longdouble) for key in SOURCE_KEYS}
    sector_new = {key: np.zeros((len(tr.N), 4), dtype=np.longdouble) for key in SOURCE_KEYS}
    entry_x, entry_w = np.polynomial.legendre.leggauss(int(entry_nodes_per_panel))
    max_identity = 0.0
    max_state_error = 0.0
    max_initial_wronskian = 0.0
    min_entry_omega2 = math.inf
    min_entry_W = math.inf
    max_W2 = 0.0
    max_W4 = 0.0
    crossings = 0

    def record(index: int) -> None:
        nonlocal max_identity
        local = _mode_sector_sources(
            float(tr.N[index]), float(tr.H[index]), float(tr.sigma[index]),
            float(tr.theta[index]), k, weights, state, p,
        )
        seed_local = _mode_sector_sources(
            float(tr.N[index]), float(tr.H[index]), float(tr.sigma[index]),
            float(tr.theta[index]), k[:seed_count], weights[:seed_count],
            {name: value[:, :seed_count] for name, value in state.items()}, p,
        )
        for key in SOURCE_KEYS:
            sector_total[key][index] = local[key]
            sector_seed[key][index] = seed_local[key]
            sector_new[key][index] = local[key] - seed_local[key]
        identity = (
            math.log(2.0) + state["log_abs_riccati_imag"]
            + 2.0 * state["log_amplitude"] - log_wronskian_target
        )
        max_identity = max(max_identity, float(np.max(np.abs(identity))))

    record(0)
    for panel in range(len(tr.N) - 1):
        left = float(tr.N[panel])
        right = float(tr.N[panel + 1])
        state = _advance_log_state(state, k, left, right, tr, max_step, p)
        midpoint = 0.5 * (left + right)
        halfwidth = 0.5 * (right - left)
        entries = midpoint + halfwidth * entry_x
        entry_weights = halfwidth * entry_w
        birth_states = []
        birth_k = []
        birth_weights = []
        for entry, quadrature_weight in zip(entries, entry_weights):
            initial, audit = _handoff_state_from_table(handoff_table, float(entry))
            k_entry = qcut * math.exp(float(entry))
            advanced = _advance_log_state(
                initial, np.array([k_entry], dtype=np.longdouble), float(entry), right, tr, max_step, p
            )
            birth_states.append(advanced)
            birth_k.append(k_entry)
            birth_weights.append(
                float(quadrature_weight) * k_entry**3 / (2.0 * math.pi**2)
            )
            min_entry_omega2 = min(min_entry_omega2, audit["min_omega2"])
            min_entry_W = min(min_entry_W, audit["min_W"])
            max_W2 = max(max_W2, audit["max_abs_W2_over_W0"])
            max_W4 = max(max_W4, audit["max_abs_W4_over_W0"])
            max_state_error = max(max_state_error, audit["state_transfer_relative_error"])
            max_initial_wronskian = max(
                max_initial_wronskian, audit["initial_wronskian_relative_error"]
            )
            crossings += 1
        for name in state:
            state[name] = np.concatenate(
                [state[name]] + [birth[name] for birth in birth_states], axis=1
            )
        k = np.concatenate([k, np.asarray(birth_k, dtype=np.longdouble)])
        weights = np.concatenate([weights, np.asarray(birth_weights, dtype=np.longdouble)])
        log_wronskian_target = np.concatenate(
            [log_wronskian_target, np.zeros((4, len(birth_k)))], axis=1
        )
        record(panel + 1)

    return {
        "sector_total": sector_total,
        "sector_seed": sector_seed,
        "sector_new": sector_new,
        "diagnostics": {
            "seed_modes": seed_count,
            "new_entry_quadrature_modes": crossings,
            "final_resolved_modes": int(len(k)),
            "max_log_wronskian_identity_error": max_identity,
            "max_wronskian_relative_error": float(abs(math.expm1(min(max_identity, 1.0)))),
            "min_entry_omega2_Mpl2": min_entry_omega2,
            "min_entry_W_Mpl": min_entry_W,
            "max_abs_entry_W2_over_W0": max_W2,
            "max_abs_entry_W4_over_W0": max_W4,
            "max_handoff_state_relative_error": max_state_error,
            "max_initial_handoff_wronskian_relative_error": max_initial_wronskian,
            "vacuum_resets": 0,
            "entry_measure": "q_cut^3*exp(3*N_entry)/(2*pi^2) dN_entry",
            "adiabatic_support_nodes": support_nodes,
            "adiabatic_momentum_support_nodes": adiabatic_momentum_support_nodes,
            "max_local_to_independent_fixed_k_state_relative_delta": handoff_table[
                "max_local_to_independent_fixed_k_state_relative_delta"
            ],
        },
    }


def _log_shells(
    qcut: float, K: float, nodes_per_octave: int
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    edges = [float(qcut)]
    while edges[-1] < K * (1.0 - 1.0e-14):
        edges.append(min(2.0 * edges[-1], K))
    x, w = np.polynomial.legendre.leggauss(int(nodes_per_octave))
    nodes: list[float] = []
    weights: list[float] = []
    segments: list[int] = []
    for segment, (left, right) in enumerate(zip(edges[:-1], edges[1:])):
        log_left, log_right = math.log(left), math.log(right)
        logq = 0.5 * (log_right - log_left) * x + 0.5 * (log_right + log_left)
        nodes.extend(np.exp(logq))
        weights.extend(0.5 * (log_right - log_left) * w)
        segments.extend([segment] * len(x))
    return np.asarray(nodes), np.asarray(weights), np.asarray(segments, dtype=int)


def _binomial_series(alpha: np.longdouble, maximum: int) -> np.ndarray:
    values = np.empty(maximum + 1, dtype=np.longdouble)
    values[0] = np.longdouble(1.0)
    for order in range(1, maximum + 1):
        values[order] = values[order - 1] * (alpha - (order - 1)) / order
    return values


_BINOMIAL_HALF = _binomial_series(np.longdouble("0.5"), 12)
_BINOMIAL_MINUS_HALF = _binomial_series(np.longdouble("-0.5"), 12)


def expanded_tail_flux_longdouble(
    primitives: dict[str, Any],
) -> dict[str, np.ndarray]:
    """Order-0/2/4 boundary flux with signed-sum-ready extended precision."""

    q = np.asarray(primitives["q"], dtype=np.longdouble)[None, :]
    mass2 = np.asarray(primitives["mass2"], dtype=np.longdouble)
    omega = np.asarray(primitives["omega"], dtype=np.longdouble)
    omega_dot = np.asarray(primitives["omega_dot"], dtype=np.longdouble)
    W2 = np.asarray(primitives["W2"], dtype=np.longdouble)
    W2_dot = np.asarray(primitives["W2_dot"], dtype=np.longdouble)
    W4 = np.asarray(primitives["W4"], dtype=np.longdouble)
    H = np.longdouble(primitives["H"])
    A0 = np.longdouble(1.0) / (np.longdouble(2.0) * omega)
    A2 = -W2 / (np.longdouble(2.0) * omega**2)
    A4 = np.longdouble(0.5) * (W2**2 / omega**3 - W4 / omega**2)
    logarithmic0 = omega_dot / omega
    logarithmic2 = W2_dot / omega - logarithmic0 * W2 / omega
    B1 = np.longdouble(1.5) * H + np.longdouble(0.5) * logarithmic0
    B3 = np.longdouble(0.5) * logarithmic2
    C0 = omega**2
    C2 = np.longdouble(2.0) * omega * W2 + B1**2
    C4 = W2**2 + np.longdouble(2.0) * omega * W4 + np.longdouble(2.0) * B1 * B3
    D0 = A0 * C0
    D2 = A0 * C2 + A2 * C0
    D4 = A0 * C4 + A2 * C2 + A4 * C0
    amplitude = A0 + A2 + A4
    kinetic = D0 + D2 + D4
    measure = q**3 / (np.longdouble(2.0) * np.longdouble(np.pi) ** 2)
    return {
        "rho": measure * np.longdouble(0.5) * (
            kinetic + (q**2 + mass2) * amplitude
        ),
        "pressure": measure * (
            np.longdouble(0.5) * kinetic
            - (q**2 / np.longdouble(6.0) + np.longdouble(0.5) * mass2) * amplitude
        ),
        "chi2": measure * amplitude,
    }


def expanded_tail_flux_table_longdouble(
    fields: dict[str, np.ndarray], q: np.ndarray, H: np.ndarray
) -> dict[str, np.ndarray]:
    """Vectorized time-by-sector-by-q form of ``expanded_tail_flux_longdouble``."""

    q3 = np.asarray(q, dtype=np.longdouble)[None, None, :]
    mass2 = np.asarray(fields["mass2"], dtype=np.longdouble)
    omega = np.asarray(fields["omega"], dtype=np.longdouble)
    omega_dot = np.asarray(fields["omega_dot"], dtype=np.longdouble)
    W2 = np.asarray(fields["W2"], dtype=np.longdouble)
    W2_dot = np.asarray(fields["W2_dot"], dtype=np.longdouble)
    W4 = np.asarray(fields["W4"], dtype=np.longdouble)
    H3 = np.asarray(H, dtype=np.longdouble)[:, None, None]
    A0 = 1.0 / (2.0 * omega)
    A2 = -W2 / (2.0 * omega**2)
    A4 = 0.5 * (W2**2 / omega**3 - W4 / omega**2)
    logarithmic0 = omega_dot / omega
    logarithmic2 = W2_dot / omega - logarithmic0 * W2 / omega
    B1 = 1.5 * H3 + 0.5 * logarithmic0
    B3 = 0.5 * logarithmic2
    C0 = omega**2
    C2 = 2.0 * omega * W2 + B1**2
    C4 = W2**2 + 2.0 * omega * W4 + 2.0 * B1 * B3
    D0 = A0 * C0
    D2 = A0 * C2 + A2 * C0
    D4 = A0 * C4 + A2 * C2 + A4 * C0
    amplitude = A0 + A2 + A4
    kinetic = D0 + D2 + D4
    measure = q3**3 / (2.0 * np.longdouble(np.pi) ** 2)
    return {
        "rho": measure * 0.5 * (kinetic + (q3**2 + mass2) * amplitude),
        "pressure": measure * (
            0.5 * kinetic - (q3**2 / 6.0 + 0.5 * mass2) * amplitude
        ),
        "chi2": measure * amplitude,
    }


def _order0_asymptotic_tail(
    trajectory: ChiTrajectory, K: float, p: ChiParameters
) -> dict[str, np.ndarray]:
    mass2 = np.asarray(portal_terms(trajectory.sigma, trajectory.theta, p)[0], dtype=np.longdouble)
    muR2 = np.longdouble((p.muR_over_Lambda * p.Lambda) ** 2)
    masses2 = mass2[:, None] + PV_J.astype(np.longdouble)[None, :] * muR2
    moments = np.stack(
        [
            np.sum(PV_C.astype(np.longdouble)[None, :] * masses2**order, axis=1)
            for order in range(13)
        ],
        axis=1,
    )
    Q = np.longdouble(K)
    pi = np.longdouble(np.pi)
    output = {key: np.zeros(len(trajectory.N), dtype=np.longdouble) for key in SOURCE_KEYS}
    last = np.zeros(len(trajectory.N), dtype=np.longdouble)
    for order in range(3, 13):
        rho_term = (
            _BINOMIAL_HALF[order] * moments[:, order] * Q ** (4 - 2 * order)
            / (np.longdouble(2 * order - 4) * 4.0 * pi**2)
        )
        output["rho"] += rho_term
        output["pressure"] += (
            _BINOMIAL_MINUS_HALF[order] * moments[:, order] * Q ** (4 - 2 * order)
            / (np.longdouble(2 * order - 4) * 12.0 * pi**2)
        )
        output["chi2"] += (
            _BINOMIAL_MINUS_HALF[order] * moments[:, order] * Q ** (2 - 2 * order)
            / (np.longdouble(2 * order - 2) * 4.0 * pi**2)
        )
        if order == 12:
            last = rho_term
    output["last_rho_term"] = last
    return output


def moving_adiabatic_tail(
    trajectory: ChiTrajectory,
    nodes_per_octave: int,
    K_over_Lambda: float,
    qcut_over_Lambda: float = 0.6,
    adiabatic_support_nodes: int | None = None,
    adiabatic_momentum_support_nodes: int = 257,
) -> dict[str, Any]:
    """Local order-0/2/4 moving physical tail from q_cut to infinity."""

    tr = trajectory.validated()
    p = ChiParameters()
    support_nodes = (
        min(len(tr.N), 65) if adiabatic_support_nodes is None
        else int(adiabatic_support_nodes)
    )
    support = resample_trajectory(tr, support_nodes)
    qcut = qcut_over_Lambda * p.Lambda
    K = K_over_Lambda * p.Lambda
    q, log_weights, segments = _log_shells(qcut, K, nodes_per_octave)
    last_segment = int(np.max(segments))
    last_mask = segments == last_segment
    targets = np.concatenate([q, np.array([qcut])])
    local = adiabatic_local_table(
        support, targets, 0.5 * qcut, 2.0 * K,
        adiabatic_momentum_support_nodes, p,
    )
    finite_fields = {name: value[:, :, :-1] for name, value in local.items()}
    flux = expanded_tail_flux_table_longdouble(finite_fields, q, support.H)
    support_sector = {
        key: np.sum(
            value * np.asarray(log_weights[None, None, :], dtype=np.longdouble),
            axis=2,
            dtype=np.longdouble,
        )
        for key, value in flux.items()
    }
    support_last_sector = {
        key: np.sum(
            value[:, :, last_mask]
            * np.asarray(log_weights[None, None, last_mask], dtype=np.longdouble),
            axis=2,
            dtype=np.longdouble,
        )
        for key, value in flux.items()
    }
    boundary_fields = {name: value[:, :, -1:] for name, value in local.items()}
    boundary_expanded_each = expanded_tail_flux_table_longdouble(
        boundary_fields, np.array([qcut]), support.H
    )["rho"][:, :, 0]
    boundary_W = np.asarray(boundary_fields["W"][:, :, 0], dtype=np.longdouble)
    boundary_Wdot = np.asarray(boundary_fields["W_dot"][:, :, 0], dtype=np.longdouble)
    boundary_mass2 = np.asarray(boundary_fields["mass2"][:, :, 0], dtype=np.longdouble)
    boundary_amplitude = 1.0 / (2.0 * boundary_W)
    boundary_B = boundary_Wdot / (2.0 * boundary_W) + 1.5 * support.H[:, None]
    boundary_kinetic = (boundary_W**2 + boundary_B**2) * boundary_amplitude
    boundary_measure = np.longdouble(qcut) ** 3 / (2.0 * np.longdouble(np.pi) ** 2)
    boundary_direct_each = boundary_measure * 0.5 * (
        boundary_kinetic
        + (np.longdouble(qcut) ** 2 + boundary_mass2) * boundary_amplitude
    )
    support_boundary_direct_rho = np.sum(
        boundary_direct_each * PV_C.astype(np.longdouble)[None, :],
        axis=1, dtype=np.longdouble,
    )
    support_boundary_expanded_rho = np.sum(
        boundary_expanded_each * PV_C.astype(np.longdouble)[None, :],
        axis=1, dtype=np.longdouble,
    )
    support_boundary_gross_rho = np.maximum(
        np.sum(
            np.abs(boundary_direct_each * PV_C.astype(np.longdouble)[None, :]),
            axis=1, dtype=np.longdouble,
        ),
        np.sum(
            np.abs(boundary_expanded_each * PV_C.astype(np.longdouble)[None, :]),
            axis=1, dtype=np.longdouble,
        ),
    )
    omega = np.asarray(finite_fields["omega"], dtype=float)
    W = np.asarray(finite_fields["W"], dtype=float)
    W2 = np.asarray(finite_fields["W2"], dtype=float)
    W4 = np.asarray(finite_fields["W4"], dtype=float)
    min_omega2 = float(np.min(omega**2))
    min_W = float(np.min(W))
    max_W2 = float(np.max(np.abs(W2 / omega)))
    max_W4 = float(np.max(np.abs(W4 / omega)))
    local_fixed_k_delta = 0.0
    for index in sorted(set((len(support.N) // 4, len(support.N) // 2, 3 * len(support.N) // 4))):
        independent = adiabatic_primitives_at(
            support, float(support.N[index]), np.array([qcut]), p
        )
        independent_flux = expanded_tail_flux_longdouble(independent)
        local_primitive = {
            "q": np.array([qcut]),
            "H": float(support.H[index]),
            **{
                name: np.asarray(value[index], dtype=float)
                for name, value in boundary_fields.items()
            },
        }
        local_flux = expanded_tail_flux_longdouble(local_primitive)
        for key in SOURCE_KEYS:
            reference_each = np.asarray(independent_flux[key][:, 0], dtype=np.longdouble)
            observed_each = np.asarray(local_flux[key][:, 0], dtype=np.longdouble)
            reference = np.sum(
                reference_each * PV_C.astype(np.longdouble), dtype=np.longdouble
            )
            observed = np.sum(
                observed_each * PV_C.astype(np.longdouble), dtype=np.longdouble
            )
            gross = max(
                np.sum(np.abs(reference_each * PV_C), dtype=np.longdouble),
                np.sum(np.abs(observed_each * PV_C), dtype=np.longdouble),
                np.longdouble(1e-300),
            )
            local_fixed_k_delta = max(
                local_fixed_k_delta,
                float(abs(observed - reference) / gross),
            )
    sector = {
        key: np.asarray(CubicSpline(support.N, value, axis=0)(tr.N), dtype=np.longdouble)
        for key, value in support_sector.items()
    }
    last_sector = {
        key: np.asarray(CubicSpline(support.N, value, axis=0)(tr.N), dtype=np.longdouble)
        for key, value in support_last_sector.items()
    }
    asymptotic = _order0_asymptotic_tail(tr, K, p)
    boundary_direct_rho = CubicSpline(
        support.N, np.asarray(support_boundary_direct_rho, dtype=float)
    )(tr.N)
    boundary_expanded_rho = CubicSpline(
        support.N, np.asarray(support_boundary_expanded_rho, dtype=float)
    )(tr.N)
    boundary_gross_rho = CubicSpline(
        support.N, np.asarray(support_boundary_gross_rho, dtype=float)
    )(tr.N)
    boundary_closure = np.abs(boundary_direct_rho - boundary_expanded_rho) / np.maximum(
        np.abs(boundary_gross_rho), 1.0e-300
    )
    signed = {}
    gross = {}
    last_ratio = 0.0
    for key in SOURCE_KEYS:
        finite_signed = np.sum(
            sector[key] * PV_C.astype(np.longdouble)[None, :], axis=1,
            dtype=np.longdouble,
        )
        signed[key] = finite_signed + asymptotic[key]
        gross[key] = (
            np.sum(
                np.abs(sector[key] * PV_C.astype(np.longdouble)[None, :]),
                axis=1,
                dtype=np.longdouble,
            )
            + np.abs(asymptotic[key])
        )
        last_signed = np.sum(
            last_sector[key] * PV_C.astype(np.longdouble)[None, :], axis=1,
            dtype=np.longdouble,
        )
        last_gross = np.sum(
            np.abs(last_sector[key] * PV_C.astype(np.longdouble)[None, :]),
            axis=1,
            dtype=np.longdouble,
        )
        last_ratio = max(
            last_ratio,
            float(np.max(np.abs(last_signed) / np.maximum(last_gross, np.longdouble(1e-300)))),
        )
    return {
        "signed": signed,
        "gross": gross,
        "sector_finite": sector,
        "asymptotic_signed": {key: asymptotic[key] for key in SOURCE_KEYS},
        "boundary_direct_rho_flux": boundary_direct_rho,
        "boundary_expanded_rho_flux": boundary_expanded_rho,
        "boundary_gross_rho_flux": boundary_gross_rho,
        "diagnostics": {
            "nodes": int(len(q)),
            "segments": int(last_segment + 1),
            "adiabatic_support_nodes": support_nodes,
            "adiabatic_momentum_support_nodes": adiabatic_momentum_support_nodes,
            "K_over_Lambda": float(K_over_Lambda),
            "min_tail_omega2_Mpl2": min_omega2,
            "min_tail_W_Mpl": min_W,
            "max_abs_tail_W2_over_W0": max_W2,
            "max_abs_tail_W4_over_W0": max_W4,
            "max_last_octave_signed_to_gross_ratio": last_ratio,
            "asymptotic_order12_rho_last_term_abs": float(np.max(np.abs(asymptotic["last_rho_term"]))),
            "max_boundary_direct_expanded_gross_relative_delta": float(
                np.max(boundary_closure)
            ),
            "max_local_to_independent_fixed_k_source_gross_relative_delta": local_fixed_k_delta,
        },
    }


def _signed_sector(sector: np.ndarray) -> np.ndarray:
    return np.asarray(
        np.sum(
            np.asarray(sector, dtype=np.longdouble)
            * PV_C.astype(np.longdouble)[None, :],
            axis=1,
            dtype=np.longdouble,
        ),
        dtype=np.longdouble,
    )


def assemble_quantum_sources(
    trajectory: ChiTrajectory,
    arrays: dict[str, np.ndarray],
    entry_nodes_per_panel: int,
    tail_nodes_per_octave: int,
    K_over_Lambda: float,
    max_mode_step: float,
    adiabatic_support_nodes: int,
    adiabatic_momentum_support_nodes: int,
) -> tuple[dict[str, np.ndarray], dict[str, Any], dict[str, Any]]:
    p = ChiParameters()
    resolved = transport_resolved_modes(
        trajectory, arrays, entry_nodes_per_panel, max_mode_step,
        adiabatic_support_nodes=adiabatic_support_nodes,
        adiabatic_momentum_support_nodes=adiabatic_momentum_support_nodes,
    )
    tail = moving_adiabatic_tail(
        trajectory, tail_nodes_per_octave, K_over_Lambda,
        adiabatic_support_nodes=adiabatic_support_nodes,
        adiabatic_momentum_support_nodes=adiabatic_momentum_support_nodes,
    )
    constants = frozen_renormalization_constants(p)
    mass2 = np.asarray(portal_terms(trajectory.sigma, trajectory.theta, p)[0], dtype=float)
    C = constants["C_chi2"]
    A = constants["A_g"]
    B = constants["B_G_rhs"]
    source: dict[str, np.ndarray] = {}
    pre_sector: dict[str, np.ndarray] = {}
    extra_signed: dict[str, np.ndarray] = {}
    for key in SOURCE_KEYS:
        pre_sector[key] = resolved["sector_total"][key] + tail["sector_finite"][key]
    extra_signed["rho"] = np.asarray(tail["asymptotic_signed"]["rho"], dtype=np.longdouble) + (
        -0.5 * C * mass2 + A + 3.0 * B * trajectory.H**2
    )
    extra_signed["pressure"] = np.asarray(tail["asymptotic_signed"]["pressure"], dtype=np.longdouble) + (
        0.5 * C * mass2 - A - B * (2.0 * trajectory.Hdot + 3.0 * trajectory.H**2)
    )
    extra_signed["chi2"] = np.asarray(tail["asymptotic_signed"]["chi2"], dtype=np.longdouble) - C
    for key in SOURCE_KEYS:
        source[key] = np.asarray(_signed_sector(pre_sector[key]) + extra_signed[key], dtype=float)

    stored = np.asarray(arrays["quantum_rho_pressure_chi2"], dtype=float)
    pre_enforcement = {key: float(source[key][0]) for key in SOURCE_KEYS}
    junction_relative: dict[str, float] = {}
    for index, key in enumerate(SOURCE_KEYS):
        gross0 = float(
            np.sum(
                np.abs(pre_sector[key][0] * PV_C.astype(np.longdouble)),
                dtype=np.longdouble,
            )
            + abs(extra_signed[key][0])
        )
        junction_relative[key] = abs(pre_enforcement[key] - stored[index]) / max(
            gross0, abs(stored[index]), 1.0e-300
        )
        source[key][0] = stored[index]

    diagnostics = {
        "resolved": resolved["diagnostics"],
        "tail": tail["diagnostics"],
        "junction_pre_enforcement": pre_enforcement,
        "junction_stored": dict(zip(SOURCE_KEYS, map(float, stored))),
        "junction_source_gross_relative_delta_by_source": junction_relative,
        "max_junction_source_gross_relative_delta": max(junction_relative.values()),
        "junction_values_enforced_exactly": all(
            source[key][0] == stored[index] for index, key in enumerate(SOURCE_KEYS)
        ),
        "renormalization_constants_unchanged": True,
        "source_values_real": all(np.isrealobj(source[key]) for key in SOURCE_KEYS),
        "source_values_finite": all(np.all(np.isfinite(source[key])) for key in SOURCE_KEYS),
    }
    internal = {
        "pre_sector": pre_sector,
        "extra_signed": extra_signed,
        "resolved": resolved,
        "tail": tail,
    }
    return source, diagnostics, internal


def integrate_background(
    arrays: dict[str, np.ndarray],
    relative_N: np.ndarray,
    quantum: dict[str, np.ndarray],
    rtol: float,
) -> tuple[ChiTrajectory, np.ndarray]:
    grid = np.asarray(relative_N, dtype=float)
    state6 = np.asarray(arrays["m1_state_physical_N"], dtype=float)
    initial = state6[:5].copy()
    N_seed = float(state6[5])
    p = ChiParameters()
    pilot = Planck2018Pilot()
    splines = {key: CubicSpline(grid, np.asarray(quantum[key], dtype=float)) for key in SOURCE_KEYS}

    def rhs(n_value: float, state: np.ndarray) -> np.ndarray:
        sigma, sigma_dot, theta, theta_dot, H = map(float, state)
        if not np.all(np.isfinite(state)) or H <= 0.0:
            raise CoupledBackgroundError("background left finite positive-H branch")
        local = classical_terms(state)
        if local["PX"] <= 0.0 or local["K"] <= 0.0:
            raise CoupledBackgroundError("background left hyperbolic branch")
        rho_standard, pressure_standard = standard_stress(
            np.array([N_seed + n_value]), pilot
        )
        _mass2, dm_sigma, dm_theta = portal_terms(sigma, theta, p)
        chi2 = float(splines["chi2"](n_value))
        sigma_ddot = (
            -3.0 * H * sigma_dot + local["P_sigma"]
            - 0.5 * float(dm_sigma) * chi2
        )
        theta_ddot = (
            local["P_theta"] - 0.5 * float(dm_theta) * chi2
            - 3.0 * H * local["PX"] * theta_dot
            - local["PX_sigma"] * sigma_dot * theta_dot
        ) / local["K"]
        Hdot = -0.5 * (
            local["rho"] + local["pressure"]
            + float(splines["rho"](n_value) + splines["pressure"](n_value))
            + float(rho_standard[0] + pressure_standard[0])
        )
        return np.array(
            [sigma_dot / H, sigma_ddot / H, theta_dot / H, theta_ddot / H, Hdot / H]
        )

    solution = solve_ivp(
        rhs,
        (float(grid[0]), float(grid[-1])),
        initial,
        t_eval=grid,
        method="DOP853",
        rtol=rtol,
        atol=np.array([2.0e-13, 2.0e-17, 2.0e-13, 2.0e-17, 2.0e-17]),
        max_step=min(2.5e-4, float(grid[-1] - grid[0]) / (len(grid) - 1)),
    )
    if not solution.success or solution.y.shape != (5, len(grid)):
        raise CoupledBackgroundError(f"coupled background integration failed: {solution.message}")
    rows = solution.y.T
    Hdot = np.empty(len(grid))
    for index, (n_value, row) in enumerate(zip(grid, rows)):
        Hdot[index] = rhs(float(n_value), row)[4] * row[4]
    Hdot[0] = float(np.asarray(arrays["Hdot_Mpl2"], dtype=float)[0])
    trajectory = ChiTrajectory(
        N=grid, H=rows[:, 4], Hdot=Hdot, sigma=rows[:, 0], theta=rows[:, 2]
    ).validated()
    return trajectory, rows


def _max_source_relative(
    left: dict[str, np.ndarray], right: dict[str, np.ndarray], mask: np.ndarray | None = None
) -> float:
    values = []
    for key in SOURCE_KEYS:
        a = np.asarray(left[key], dtype=float)
        b = np.asarray(right[key], dtype=float)
        if mask is not None:
            a, b = a[mask], b[mask]
        values.append(float(np.max(np.abs(a - b))) / max(float(np.max(np.abs(a))), float(np.max(np.abs(b))), 1e-300))
    return max(values)


def source_gross_scale(
    internal: dict[str, Any], key: str
) -> np.ndarray:
    """Absolute signed-basis scale before PV/counterterm cancellations."""

    return np.asarray(
        np.sum(
            np.abs(
                np.asarray(internal["pre_sector"][key], dtype=np.longdouble)
                * PV_C.astype(np.longdouble)[None, :]
            ),
            axis=1,
            dtype=np.longdouble,
        )
        + np.abs(np.asarray(internal["extra_signed"][key], dtype=np.longdouble)),
        dtype=float,
    )


def _max_source_gross_relative(
    left: dict[str, np.ndarray],
    right: dict[str, np.ndarray],
    internal: dict[str, Any],
    mask: np.ndarray | None = None,
) -> float:
    values = []
    for key in SOURCE_KEYS:
        a = np.asarray(left[key], dtype=float)
        b = np.asarray(right[key], dtype=float)
        gross = source_gross_scale(internal, key)
        if mask is not None:
            a, b, gross = a[mask], b[mask], gross[mask]
        values.append(
            float(np.max(np.abs(a - b)))
            / max(float(np.max(np.abs(gross))), 1.0e-300)
        )
    return max(values)


def run_coupled(
    arrays: dict[str, np.ndarray],
    span_N: float,
    nodes: int,
    config: CoupledConfig,
    *,
    background_rtol: float,
    max_mode_step: float,
    entry_nodes_per_panel: int,
    tail_nodes_per_octave: int,
    K_over_Lambda: float,
    adiabatic_support_nodes: int,
    adiabatic_momentum_support_nodes: int,
) -> dict[str, Any]:
    grid = np.linspace(0.0, span_N, nodes)
    stored = np.asarray(arrays["quantum_rho_pressure_chi2"], dtype=float)
    quantum = {
        key: np.full(nodes, stored[index], dtype=float)
        for index, key in enumerate(SOURCE_KEYS)
    }
    iteration_log = []
    for iteration in range(config.fixed_point_steps):
        used = {key: value.copy() for key, value in quantum.items()}
        trajectory, rows = integrate_background(arrays, grid, used, background_rtol)
        fresh, source_diagnostics, fresh_internal = assemble_quantum_sources(
            trajectory, arrays, entry_nodes_per_panel, tail_nodes_per_octave,
            K_over_Lambda, max_mode_step, adiabatic_support_nodes,
            adiabatic_momentum_support_nodes,
        )
        change = _max_source_gross_relative(fresh, used, fresh_internal)
        signed_change = _max_source_relative(fresh, used)
        quantum = {
            key: config.relaxation * fresh[key] + (1.0 - config.relaxation) * used[key]
            for key in SOURCE_KEYS
        }
        for index, key in enumerate(SOURCE_KEYS):
            quantum[key][0] = stored[index]
        iteration_log.append(
            {
                "iteration": iteration + 1,
                "source_gross_relative_change": change,
                "source_signed_relative_change_diagnostic_not_gated": signed_change,
                "H_endpoint_Mpl": float(trajectory.H[-1]),
                "max_wronskian_relative_error": source_diagnostics["resolved"]["max_wronskian_relative_error"],
                "min_entry_omega2_Mpl2": source_diagnostics["resolved"]["min_entry_omega2_Mpl2"],
            }
        )
    used = {key: value.copy() for key, value in quantum.items()}
    trajectory, rows = integrate_background(arrays, grid, used, background_rtol)
    fresh, source_diagnostics, internal = assemble_quantum_sources(
        trajectory, arrays, entry_nodes_per_panel, tail_nodes_per_octave,
        K_over_Lambda, max_mode_step, adiabatic_support_nodes,
        adiabatic_momentum_support_nodes,
    )
    validation_change = _max_source_gross_relative(fresh, used, internal)
    validation_signed_change = _max_source_relative(fresh, used)
    return {
        "trajectory": trajectory,
        "background_rows": rows,
        "quantum_used": used,
        "quantum_fresh": fresh,
        "source_diagnostics": source_diagnostics,
        "source_internal": internal,
        "iteration_log": iteration_log,
        "validation_source_gross_relative_change": validation_change,
        "validation_source_signed_relative_change_diagnostic_not_gated": validation_signed_change,
        "config": {
            "span_N": span_N,
            "nodes": nodes,
            "background_rtol": background_rtol,
            "max_mode_step_N": max_mode_step,
            "entry_nodes_per_panel": entry_nodes_per_panel,
            "tail_nodes_per_octave": tail_nodes_per_octave,
            "K_over_Lambda": K_over_Lambda,
            "adiabatic_support_nodes": adiabatic_support_nodes,
            "adiabatic_momentum_support_nodes": adiabatic_momentum_support_nodes,
        },
    }


def _evaluation_mask(grid: np.ndarray, start: float) -> np.ndarray:
    mask = np.asarray(grid >= start, dtype=bool)
    indices = np.flatnonzero(mask)
    if len(indices) < 5:
        raise CoupledBackgroundError("evaluation window has too few nodes")
    mask[indices[-2:]] = False
    if np.count_nonzero(mask) < 3:
        raise CoupledBackgroundError("evaluation derivative margin is empty")
    return mask


def physical_diagnostics(
    run: dict[str, Any], arrays: dict[str, np.ndarray], evaluation_start: float
) -> dict[str, Any]:
    tr: ChiTrajectory = run["trajectory"]
    q = run["quantum_fresh"]
    mask = _evaluation_mask(tr.N, evaluation_start)
    pilot = Planck2018Pilot()
    N_seed = float(np.asarray(arrays["m1_state_physical_N"], dtype=float)[5])
    rows = np.asarray(run["background_rows"], dtype=float)
    classical = [classical_terms(row) for row in rows]
    rho_standard, pressure_standard = standard_stress(N_seed + tr.N, pilot)
    total_density = (
        np.array([item["rho"] for item in classical]) + rho_standard + q["rho"]
    )
    friedmann_residual = 3.0 * tr.H**2 - total_density
    friedmann_scale = np.maximum(
        3.0 * tr.H**2 + np.abs(total_density), 1.0e-300
    )
    independent_Hdot = tr.H * CubicSpline(tr.N, tr.H)(tr.N, 1)
    rhs_Hdot = -0.5 * (
        np.array([item["rho"] + item["pressure"] for item in classical])
        + q["rho"] + q["pressure"] + rho_standard + pressure_standard
    )
    ray_scale = np.maximum(np.abs(independent_Hdot) + np.abs(rhs_Hdot), 1.0e-300)
    p = ChiParameters()
    mass2 = np.asarray(portal_terms(tr.sigma, tr.theta, p)[0], dtype=float)
    mass2_dot = tr.H * CubicSpline(tr.N, mass2)(tr.N, 1)
    internal = run["source_internal"]
    tail = internal["tail"]
    tail_source = {
        key: np.asarray(
            _signed_sector(tail["sector_finite"][key])
            + tail["asymptotic_signed"][key],
            dtype=float,
        )
        for key in SOURCE_KEYS
    }
    # The resolved bulk Ward identity follows independently from the transported
    # mode equation.  Its moving-domain term is +H times the direct boundary
    # flux.  The local tail is differentiated separately and supplies the
    # opposite expanded-tail flux.  Their difference is precisely the G28
    # handoff remainder, without differentiating a 15-digit cancelled total.
    tail_rho_dot = tr.H * CubicSpline(tr.N, tail_source["rho"])(tr.N, 1)
    tail_enthalpy = 3.0 * tr.H * (
        tail_source["rho"] + tail_source["pressure"]
    )
    tail_portal = 0.5 * mass2_dot * tail_source["chi2"]
    boundary_direct = tr.H * np.asarray(tail["boundary_direct_rho_flux"], dtype=float)
    boundary_expanded = tr.H * np.asarray(tail["boundary_expanded_rho_flux"], dtype=float)
    ward = tail_rho_dot + tail_enthalpy - tail_portal + boundary_direct
    ward_scale = np.maximum(
        np.abs(tail_rho_dot) + np.abs(tail_enthalpy) + np.abs(tail_portal)
        + np.abs(boundary_expanded),
        1.0e-300,
    )
    cancelled_rho_dot = tr.H * CubicSpline(tr.N, q["rho"])(tr.N, 1)
    cancelled_ward = (
        cancelled_rho_dot + 3.0 * tr.H * (q["rho"] + q["pressure"])
        - 0.5 * mass2_dot * q["chi2"]
    )
    cancelled_scale = np.maximum(
        np.abs(cancelled_rho_dot)
        + np.abs(3.0 * tr.H * (q["rho"] + q["pressure"]))
        + np.abs(0.5 * mass2_dot * q["chi2"]),
        1.0e-300,
    )
    standard_derivative = tr.H * CubicSpline(tr.N, rho_standard)(tr.N, 1)
    standard_residual = standard_derivative + 3.0 * tr.H * (
        rho_standard + pressure_standard
    )
    standard_scale = np.maximum(
        np.abs(standard_derivative)
        + np.abs(3.0 * tr.H * (rho_standard + pressure_standard)),
        1.0e-300,
    )
    local_PX = np.array([item["PX"] for item in classical])
    local_K = np.array([item["K"] for item in classical])
    stored_Hdot = float(np.asarray(arrays["Hdot_Mpl2"])[0])
    junction_rhs = float(rhs_Hdot[0])
    junction_ray = abs(junction_rhs - stored_Hdot) / max(
        3.0 * float(tr.H[0]) ** 2,
        abs(junction_rhs) + abs(stored_Hdot),
        1.0e-300,
    )
    return {
        "evaluation_start_N": float(evaluation_start),
        "evaluation_nodes": int(np.count_nonzero(mask)),
        "all_background_values_finite": bool(
            all(np.all(np.isfinite(value)) for value in (tr.H, tr.Hdot, tr.sigma, tr.theta))
        ),
        "positive_H": bool(np.all(tr.H > 0.0)),
        "min_total_density_Mpl4": float(np.min(total_density)),
        "min_PX": float(np.min(local_PX)),
        "min_K": float(np.min(local_K)),
        "max_friedmann_normalized": float(np.max(np.abs(friedmann_residual[mask]) / friedmann_scale[mask])),
        "max_raychaudhuri_normalized": float(np.max(np.abs(independent_Hdot[mask] - rhs_Hdot[mask]) / ray_scale[mask])),
        "ward_normalized": float(np.max(np.abs(ward[mask]) / ward_scale[mask])),
        "ward_evaluation": "resolved_bulk_mode_EOM_plus_direct_boundary_flux_plus_independently_differentiated_local_tail",
        "cancelled_total_finite_difference_ward_diagnostic_not_gated": float(
            np.max(np.abs(cancelled_ward[mask]) / cancelled_scale[mask])
        ),
        "max_boundary_direct_expanded_gross_relative_delta": source_boundary
        if (source_boundary := float(
            tail["diagnostics"]["max_boundary_direct_expanded_gross_relative_delta"]
        )) >= 0.0 else math.inf,
        "standard_continuity_normalized": float(np.max(np.abs(standard_residual[mask]) / standard_scale[mask])),
        "junction_raychaudhuri_normalized": float(junction_ray),
        "H_endpoint_Mpl": float(tr.H[-1]),
        "delta_H_over_H": float(tr.H[-1] / tr.H[0] - 1.0),
        "physical_N_endpoint": N_seed + float(tr.N[-1]),
        "present_anchor_reached": False,
    }


def _interpolate_run_values(
    fine: dict[str, Any], coarse_grid: np.ndarray
) -> tuple[np.ndarray, dict[str, np.ndarray], dict[str, np.ndarray]]:
    tr: ChiTrajectory = fine["trajectory"]
    rows = np.asarray(fine["background_rows"], dtype=float)
    interpolated_rows = CubicSpline(tr.N, rows, axis=0)(coarse_grid)
    sources = {
        key: CubicSpline(tr.N, fine["quantum_fresh"][key])(coarse_grid)
        for key in SOURCE_KEYS
    }
    gross = {
        key: CubicSpline(
            tr.N, source_gross_scale(fine["source_internal"], key)
        )(coarse_grid)
        for key in SOURCE_KEYS
    }
    return interpolated_rows, sources, gross


def time_convergence_metrics(
    fine: dict[str, Any], coarse: dict[str, Any], evaluation_start: float
) -> dict[str, Any]:
    coarse_tr: ChiTrajectory = coarse["trajectory"]
    fine_rows, fine_sources, fine_gross = _interpolate_run_values(fine, coarse_tr.N)
    coarse_rows = np.asarray(coarse["background_rows"], dtype=float)
    mask = _evaluation_mask(coarse_tr.N, evaluation_start)
    per_field = {}
    for index, name in enumerate(("sigma", "sigmadot", "theta", "thetadot", "H")):
        scale = max(
            float(np.max(np.abs(fine_rows[mask, index]))),
            float(np.max(np.abs(coarse_rows[mask, index]))), 1.0e-300,
        )
        per_field[name] = float(np.max(np.abs(fine_rows[mask, index] - coarse_rows[mask, index]))) / scale
    per_source = {}
    for key in SOURCE_KEYS:
        a = np.asarray(fine_sources[key])[mask]
        b = np.asarray(coarse["quantum_fresh"][key])[mask]
        coarse_gross = source_gross_scale(coarse["source_internal"], key)[mask]
        per_source[key] = float(np.max(np.abs(a - b))) / max(
            float(np.max(np.abs(fine_gross[key][mask]))),
            float(np.max(np.abs(coarse_gross))),
            1.0e-300,
        )
    return {
        "background_relative_delta_by_field": per_field,
        "source_gross_relative_delta_by_source": per_source,
        "max_background_time_relative_delta": max(per_field.values()),
        "max_source_time_gross_relative_delta": max(per_source.values()),
    }


def _replacement_metric(
    reference: dict[str, np.ndarray], replacement: dict[str, np.ndarray],
    gross_scale: dict[str, np.ndarray], mask: np.ndarray
) -> tuple[dict[str, float], float]:
    by_source = {}
    for key in SOURCE_KEYS:
        a = np.asarray(reference[key], dtype=float)[mask]
        b = np.asarray(replacement[key], dtype=float)[mask]
        gross = np.asarray(gross_scale[key], dtype=float)[mask]
        by_source[key] = float(np.max(np.abs(a - b))) / max(
            float(np.max(np.abs(gross))), 1.0e-300
        )
    return by_source, max(by_source.values())


def operator_convergence_metrics(
    run: dict[str, Any], arrays: dict[str, np.ndarray], config: CoupledConfig,
    evaluation_start: float,
) -> dict[str, Any]:
    tr: ChiTrajectory = run["trajectory"]
    support_nodes = int(run["config"]["adiabatic_support_nodes"])
    mask = _evaluation_mask(tr.N, evaluation_start)
    reference = run["quantum_fresh"]
    ref_internal = run["source_internal"]
    entry_coarse = transport_resolved_modes(
        tr, arrays, config.coarse_entry_nodes_per_panel, config.coarse_mode_step_N,
        adiabatic_support_nodes=support_nodes,
        adiabatic_momentum_support_nodes=config.adiabatic_momentum_support_nodes,
    )
    entry_source = {}
    for key in SOURCE_KEYS:
        replacement_sector = (
            entry_coarse["sector_total"][key]
            + ref_internal["tail"]["sector_finite"][key]
        )
        entry_source[key] = np.asarray(
            _signed_sector(replacement_sector)
            + ref_internal["tail"]["asymptotic_signed"][key]
            + (ref_internal["extra_signed"][key] - ref_internal["tail"]["asymptotic_signed"][key]),
            dtype=float,
        )
        entry_source[key][0] = reference[key][0]
    tail_coarse = moving_adiabatic_tail(
        tr, config.coarse_tail_nodes_per_octave, config.fine_K_over_Lambda,
        adiabatic_support_nodes=support_nodes,
        adiabatic_momentum_support_nodes=config.adiabatic_momentum_support_nodes,
    )
    tail_source = {}
    for key in SOURCE_KEYS:
        replacement_sector = (
            ref_internal["resolved"]["sector_total"][key]
            + tail_coarse["sector_finite"][key]
        )
        counterterm = (
            ref_internal["extra_signed"][key]
            - ref_internal["tail"]["asymptotic_signed"][key]
        )
        tail_source[key] = np.asarray(
            _signed_sector(replacement_sector) + tail_coarse["asymptotic_signed"][key] + counterterm,
            dtype=float,
        )
        tail_source[key][0] = reference[key][0]
    uv_coarse = moving_adiabatic_tail(
        tr, config.fine_tail_nodes_per_octave, config.coarse_K_over_Lambda,
        adiabatic_support_nodes=support_nodes,
        adiabatic_momentum_support_nodes=config.adiabatic_momentum_support_nodes,
    )
    uv_source = {}
    for key in SOURCE_KEYS:
        replacement_sector = (
            ref_internal["resolved"]["sector_total"][key]
            + uv_coarse["sector_finite"][key]
        )
        counterterm = (
            ref_internal["extra_signed"][key]
            - ref_internal["tail"]["asymptotic_signed"][key]
        )
        uv_source[key] = np.asarray(
            _signed_sector(replacement_sector) + uv_coarse["asymptotic_signed"][key] + counterterm,
            dtype=float,
        )
        uv_source[key][0] = reference[key][0]
    gross = {key: source_gross_scale(ref_internal, key) for key in SOURCE_KEYS}
    entry_by, entry_max = _replacement_metric(reference, entry_source, gross, mask)
    tail_by, tail_max = _replacement_metric(reference, tail_source, gross, mask)
    uv_by, uv_max = _replacement_metric(reference, uv_source, gross, mask)
    return {
        "entry_quadrature_gross_relative_delta_by_source": entry_by,
        "tail_quadrature_gross_relative_delta_by_source": tail_by,
        "uv_cutoff_gross_relative_delta_by_source": uv_by,
        "max_entry_quadrature_gross_relative_delta": entry_max,
        "max_tail_quadrature_gross_relative_delta": tail_max,
        "max_uv_cutoff_gross_relative_delta": uv_max,
        "coarse_entry_diagnostics": entry_coarse["diagnostics"],
        "coarse_tail_diagnostics": tail_coarse["diagnostics"],
        "coarse_uv_diagnostics": uv_coarse["diagnostics"],
    }


def decimal_precision_metrics(
    run: dict[str, Any], evaluation_start: float, low_digits: int, high_digits: int
) -> dict[str, Any]:
    tr: ChiTrajectory = run["trajectory"]
    mask = _evaluation_mask(tr.N, evaluation_start)
    candidates = np.flatnonzero(mask)
    sample_indices = tuple(sorted(set((int(candidates[0]), int(candidates[len(candidates)//2]), int(candidates[-1])))))
    internal = run["source_internal"]
    max_precision = 0.0
    max_float = 0.0
    min_cancel = math.inf
    max_cancel = 0.0
    points = []
    coefficients = tuple(Decimal(int(value)) for value in (1, -3, 3, -1))
    for index in sample_indices:
        point = {"N_relative": float(tr.N[index]), "sources": {}}
        for key in SOURCE_KEYS:
            sector = np.asarray(internal["pre_sector"][key][index], dtype=float)
            extra = float(internal["extra_signed"][key][index])
            results = {}
            for digits in (low_digits, high_digits):
                with localcontext() as context:
                    context.prec = digits
                    terms = [coefficient * Decimal(repr(float(value))) for coefficient, value in zip(coefficients, sector)]
                    signed = sum(terms, Decimal(0)) + Decimal(repr(extra))
                    gross = sum((abs(value) for value in terms), Decimal(0)) + abs(Decimal(repr(extra)))
                    results[digits] = (signed, gross)
            low, _ = results[low_digits]
            high, gross = results[high_digits]
            precision = float(abs(low - high) / max(gross, Decimal("1e-300")))
            float_signed = float(np.dot(PV_C, sector) + extra)
            float_delta = float(abs(Decimal(repr(float_signed)) - high) / max(gross, Decimal("1e-300")))
            cancellation = float(abs(high) / max(gross, Decimal("1e-300")))
            max_precision = max(max_precision, precision)
            max_float = max(max_float, float_delta)
            min_cancel = min(min_cancel, cancellation)
            max_cancel = max(max_cancel, cancellation)
            point["sources"][key] = {
                "signed_decimal80": str(high),
                "gross_abs_decimal80": str(gross),
                "decimal_60_to_80_gross_relative_delta": precision,
                "float64_to_decimal80_gross_relative_delta": float_delta,
                "cancellation_ratio": cancellation,
            }
        points.append(point)
    return {
        "sample_points": points,
        "max_decimal_60_to_80_gross_relative_delta": max_precision,
        "max_float64_to_decimal80_gross_relative_delta": max_float,
        "min_PV_cancellation_ratio": min_cancel,
        "max_PV_cancellation_ratio": max_cancel,
    }


def evaluate_ensemble(
    arrays: dict[str, np.ndarray], config: CoupledConfig, *, heldout: bool
) -> dict[str, Any]:
    span = config.heldout_span_N if heldout else config.pilot_span_N
    fine_nodes = config.heldout_fine_nodes if heldout else config.pilot_fine_nodes
    coarse_nodes = (fine_nodes + 1) // 2
    fine_support_nodes = (fine_nodes - 1) // config.adiabatic_support_stride + 1
    coarse_support_nodes = (coarse_nodes - 1) // config.adiabatic_support_stride + 1
    evaluation_start = config.heldout_evaluation_start_N if heldout else span / 4.0
    fine = run_coupled(
        arrays, span, fine_nodes, config,
        background_rtol=config.fine_background_rtol,
        max_mode_step=config.fine_mode_step_N,
        entry_nodes_per_panel=config.fine_entry_nodes_per_panel,
        tail_nodes_per_octave=config.fine_tail_nodes_per_octave,
        K_over_Lambda=config.fine_K_over_Lambda,
        adiabatic_support_nodes=fine_support_nodes,
        adiabatic_momentum_support_nodes=config.adiabatic_momentum_support_nodes,
    )
    coarse = run_coupled(
        arrays, span, coarse_nodes, config,
        background_rtol=config.coarse_background_rtol,
        max_mode_step=config.coarse_mode_step_N,
        entry_nodes_per_panel=config.fine_entry_nodes_per_panel,
        tail_nodes_per_octave=config.fine_tail_nodes_per_octave,
        K_over_Lambda=config.fine_K_over_Lambda,
        adiabatic_support_nodes=coarse_support_nodes,
        adiabatic_momentum_support_nodes=config.adiabatic_momentum_support_nodes,
    )
    physical = physical_diagnostics(fine, arrays, evaluation_start)
    time_metrics = time_convergence_metrics(fine, coarse, evaluation_start)
    operator = operator_convergence_metrics(fine, arrays, config, evaluation_start)
    precision = decimal_precision_metrics(
        fine, evaluation_start, config.decimal_low_digits, config.decimal_high_digits
    )
    source_diag = fine["source_diagnostics"]
    metrics = {
        "max_fixed_point_source_gross_relative_change": fine[
            "validation_source_gross_relative_change"
        ],
        **{key: time_metrics[key] for key in (
            "max_background_time_relative_delta", "max_source_time_gross_relative_delta"
        )},
        **{key: operator[key] for key in (
            "max_entry_quadrature_gross_relative_delta",
            "max_tail_quadrature_gross_relative_delta",
            "max_uv_cutoff_gross_relative_delta",
        )},
        **{key: physical[key] for key in (
            "max_friedmann_normalized", "max_raychaudhuri_normalized",
            "ward_normalized", "standard_continuity_normalized",
        )},
    }
    return {
        "phase": "heldout" if heldout else "pilot",
        "fine_run_config": fine["config"],
        "coarse_time_run_config": coarse["config"],
        "fine_iteration_log": fine["iteration_log"],
        "coarse_iteration_log": coarse["iteration_log"],
        "physical_diagnostics": physical,
        "time_convergence": time_metrics,
        "operator_convergence": operator,
        "precision": precision,
        "source_diagnostics": source_diag,
        "metrics_for_freeze_or_gate": metrics,
        "trajectory_rows_persisted": 0,
        "new_AP1_M1_background_runs": 2,
    }


def _absolute_gates(ensemble: dict[str, Any]) -> dict[str, bool]:
    physical = ensemble["physical_diagnostics"]
    source = ensemble["source_diagnostics"]
    precision = ensemble["precision"]
    tail = source["tail"]
    resolved = source["resolved"]
    return {
        "finite_real_quantum_and_background": bool(
            physical["all_background_values_finite"]
            and source["source_values_real"] and source["source_values_finite"]
        ),
        "positive_expanding_branch": bool(
            physical["positive_H"] and physical["min_total_density_Mpl4"] > 0.0
        ),
        "hyperbolicity_PX_and_K_positive": bool(
            physical["min_PX"] > 0.0 and physical["min_K"] > 0.0
        ),
        "entry_and_tail_frequencies_positive": bool(
            resolved["min_entry_omega2_Mpl2"] > 0.0
            and resolved["min_entry_W_Mpl"] > 0.0
            and tail["min_tail_omega2_Mpl2"] > 0.0
            and tail["min_tail_W_Mpl"] > 0.0
        ),
        "no_vacuum_reset_and_exact_handoff": bool(
            resolved["vacuum_resets"] == 0
            and resolved["max_handoff_state_relative_error"]
            <= FREEZE_POLICY["handoff_state_relative_cap"]
        ),
        "local_convective_derivative_matches_independent_fixed_k_reference": bool(
            resolved["max_local_to_independent_fixed_k_state_relative_delta"]
            <= FREEZE_POLICY["local_fixed_k_state_relative_cap"]
            and tail["max_local_to_independent_fixed_k_source_gross_relative_delta"]
            <= FREEZE_POLICY["local_fixed_k_source_gross_relative_cap"]
        ),
        "direct_and_expanded_boundary_flux_close": bool(
            tail["max_boundary_direct_expanded_gross_relative_delta"]
            <= FREEZE_POLICY["boundary_direct_expanded_gross_relative_cap"]
        ),
        "wronskian_below_existing_limit": bool(
            max(
                resolved["max_wronskian_relative_error"],
                resolved["max_initial_handoff_wronskian_relative_error"],
            ) <= WRONSKIAN_LIMIT
        ),
        "junction_quantum_triplet_enforced_exactly": bool(
            source["junction_values_enforced_exactly"]
            and source["max_junction_source_gross_relative_delta"]
            <= FREEZE_POLICY["junction_source_relative_cap"]
        ),
        "junction_raychaudhuri_within_predeclared_cap": bool(
            physical["junction_raychaudhuri_normalized"]
            <= FREEZE_POLICY["junction_raychaudhuri_relative_cap"]
        ),
        "decimal_precision_below_predeclared_cap": bool(
            precision["max_decimal_60_to_80_gross_relative_delta"]
            <= FREEZE_POLICY["decimal_60_to_80_gross_relative_cap"]
        ),
        "float64_decimal80_agreement_below_predeclared_cap": bool(
            precision["max_float64_to_decimal80_gross_relative_delta"]
            <= FREEZE_POLICY["float64_to_decimal80_gross_relative_cap"]
        ),
        "last_uv_octave_below_predeclared_cap": bool(
            tail["max_last_octave_signed_to_gross_ratio"]
            <= FREEZE_POLICY["last_uv_octave_gross_relative_cap"]
        ),
        "present_anchor_not_claimed": not physical["present_anchor_reached"],
    }


def build_pilot(
    apeiron_root: Path, code_path: Path, config: CoupledConfig | None = None
) -> dict[str, Any]:
    config = CoupledConfig() if config is None else config
    config.validate()
    seed_report, arrays, authorities = load_inputs(apeiron_root)
    ensemble = evaluate_ensemble(arrays, config, heldout=False)
    absolute = _absolute_gates(ensemble)
    numerical = {
        f"{name}_below_predeclared_ceiling": value <= policy["ceiling"]
        for name, value in ensemble["metrics_for_freeze_or_gate"].items()
        for policy in (FREEZE_POLICY["metrics"][name],)
    }
    gates = {
        "G28_checkpoint_hash_bound_and_PASS": True,
        "chronology_safe_128_seed_hash_bound_and_PASS": bool(seed_report["all_seed_gates_pass"]),
        "pilot_and_heldout_windows_disjoint": config.pilot_span_N < config.heldout_evaluation_start_N,
        "frozen_qcut_and_PV_basis_unchanged": bool(
            config.qcut_over_Lambda == 0.6
            and np.array_equal(PV_C, np.array([1.0, -3.0, 3.0, -1.0]))
            and np.array_equal(PV_J, np.array([0.0, 1.0, 2.0, 3.0]))
        ),
        **absolute,
        **numerical,
        "no_two_time_kernel_seed_observable_or_curve": True,
    }
    passed = bool(all(gates.values()))
    return {
        "schema": "apeiron-ap1-m1-coupled-moving-split-background-pilot-v1.0",
        "updated_utc": _utc_now(),
        "classification": (
            "COUPLED_MOVING_SPLIT_BACKGROUND_PILOT_COMPLETE_TOLERANCE_FREEZE_OPEN"
            if passed else "COUPLED_MOVING_SPLIT_BACKGROUND_PILOT_INCOMPLETE"
        ),
        "authority_sha256": {
            **authorities,
            "AP1/CODE/ap1_m1_coupled_moving_split_background.py": file_sha256(code_path),
        },
        "seed_classification": seed_report["classification"],
        "config": asdict(config),
        "predeclared_freeze_policy": FREEZE_POLICY,
        "pilot_ensemble": ensemble,
        "gates": gates,
        "all_pilot_gates_pass": passed,
        "heldout_evaluated": False,
        "checkpoint_eligible": False,
        "bounded_background_pilot_started": True,
        "new_AP1_M1_background_runs": ensemble["new_AP1_M1_background_runs"],
        "background_tolerances_frozen": False,
        "long_background_candidate_assessed": False,
        "trajectory_rows_persisted": 0,
        "physical_response_kernel_started": False,
        "seed_released": False,
        "nonpass_stored_or_used_as_seed": False,
        "equations_changed": False,
        "physics_changed": False,
        "parameters_changed": False,
        "existing_gate_thresholds_changed": False,
        "memory_convergence_applicable": False,
        "claim_boundary": (
            "bounded coupled background pilot only; no present-anchor reach, "
            "two-time memory kernel, observable, production curve, fit or significance"
        ),
        "next_required": (
            "freeze pilot-derived G29 tolerances before the registered disjoint "
            "held-out interval is evaluated"
        ),
    }


def build_freeze(pilot_path: Path, code_path: Path) -> dict[str, Any]:
    pilot = json.loads(pilot_path.read_text(encoding="utf-8"))
    if not pilot.get("all_pilot_gates_pass") or pilot.get("heldout_evaluated"):
        raise CoupledBackgroundError("complete pre-held-out pilot required")
    code_hash = file_sha256(code_path)
    if pilot["authority_sha256"].get(
        "AP1/CODE/ap1_m1_coupled_moving_split_background.py"
    ) != code_hash:
        raise CoupledBackgroundError("coupled background code changed after pilot")
    metrics = pilot["pilot_ensemble"]["metrics_for_freeze_or_gate"]
    thresholds = {
        name: max(float(policy["floor"]), float(policy["multiplier"]) * float(metrics[name]))
        for name, policy in FREEZE_POLICY["metrics"].items()
    }
    gates = {
        "pilot_hash_bound_before_heldout": True,
        **{
            f"{name}_threshold_below_predeclared_ceiling": thresholds[name] <= policy["ceiling"]
            for name, policy in FREEZE_POLICY["metrics"].items()
        },
        "heldout_interval_registered_disjoint_and_unseen": True,
        "physical_kernel_remains_locked": True,
    }
    passed = bool(all(gates.values()))
    config = pilot["config"]
    return {
        "schema": "apeiron-ap1-m1-coupled-moving-split-background-freeze-v1.0",
        "updated_utc": _utc_now(),
        "classification": (
            "COUPLED_MOVING_SPLIT_BACKGROUND_TOLERANCES_FROZEN_HELDOUT_UNSEEN"
            if passed else "COUPLED_MOVING_SPLIT_BACKGROUND_TOLERANCE_FREEZE_FAILED"
        ),
        "input_sha256": {
            pilot_path.name: file_sha256(pilot_path),
            "AP1/CODE/ap1_m1_coupled_moving_split_background.py": code_hash,
        },
        "predeclared_freeze_policy": FREEZE_POLICY,
        "thresholds": thresholds,
        "heldout_plan": {
            "span_N": config["heldout_span_N"],
            "fine_nodes": config["heldout_fine_nodes"],
            "coarse_nodes": (config["heldout_fine_nodes"] + 1) // 2,
            "evaluation_start_N": config["heldout_evaluation_start_N"],
            "fine_entry_nodes_per_panel": config["fine_entry_nodes_per_panel"],
            "fine_tail_nodes_per_octave": config["fine_tail_nodes_per_octave"],
            "fine_K_over_Lambda": config["fine_K_over_Lambda"],
            "decimal_high_digits": config["decimal_high_digits"],
        },
        "gates": gates,
        "all_tolerance_freeze_gates_pass": passed,
        "tolerances_frozen_before_heldout": passed,
        "heldout_evaluated": False,
        "new_AP1_M1_background_runs": pilot["new_AP1_M1_background_runs"],
        "trajectory_rows_persisted": 0,
        "physical_response_kernel_started": False,
        "claim_boundary": "new bounded-background numerical tolerances only; heldout, long background and physical kernel remain unseen",
        "next_required": "evaluate the registered disjoint held-out interval once and persist only a full PASS",
    }


def _config_from_pilot(pilot: dict[str, Any]) -> CoupledConfig:
    config = CoupledConfig(**pilot["config"])
    config.validate()
    return config


def build_heldout(
    apeiron_root: Path, pilot_path: Path, freeze_path: Path, code_path: Path
) -> dict[str, Any]:
    pilot = json.loads(pilot_path.read_text(encoding="utf-8"))
    freeze = json.loads(freeze_path.read_text(encoding="utf-8"))
    if not pilot.get("all_pilot_gates_pass"):
        raise CoupledBackgroundError("pilot is not PASS")
    if not freeze.get("all_tolerance_freeze_gates_pass"):
        raise CoupledBackgroundError("tolerance freeze is not PASS")
    code_hash = file_sha256(code_path)
    if freeze["input_sha256"].get(pilot_path.name) != file_sha256(pilot_path):
        raise CoupledBackgroundError("pilot changed after tolerance freeze")
    if freeze["input_sha256"].get(
        "AP1/CODE/ap1_m1_coupled_moving_split_background.py"
    ) != code_hash:
        raise CoupledBackgroundError("code changed after tolerance freeze")
    config = _config_from_pilot(pilot)
    plan = freeze["heldout_plan"]
    if not (
        plan["span_N"] == config.heldout_span_N
        and plan["fine_nodes"] == config.heldout_fine_nodes
        and plan["evaluation_start_N"] == config.heldout_evaluation_start_N
    ):
        raise CoupledBackgroundError("held-out plan changed after freeze")
    _, arrays, authorities = load_inputs(apeiron_root)
    ensemble = evaluate_ensemble(arrays, config, heldout=True)
    absolute = _absolute_gates(ensemble)
    metrics = ensemble["metrics_for_freeze_or_gate"]
    thresholds = freeze["thresholds"]
    numerical = {
        f"{name}_below_frozen_threshold": float(value) <= float(thresholds[name])
        for name, value in metrics.items()
    }
    gates = {
        "pilot_and_freeze_hash_bound_before_heldout": True,
        "registered_disjoint_heldout_interval_used_once": True,
        **absolute,
        **numerical,
        "no_two_time_kernel_seed_observable_or_curve": True,
    }
    passed = bool(all(gates.values()))
    return {
        "schema": "apeiron-ap1-m1-coupled-moving-split-background-heldout-v1.0",
        "updated_utc": _utc_now(),
        "classification": (
            "COUPLED_MOVING_SPLIT_BOUNDED_BACKGROUND_HELDOUT_PASS_LONG_EXTENSION_OPEN"
            if passed else "COUPLED_MOVING_SPLIT_BOUNDED_BACKGROUND_HELDOUT_NONPASS"
        ),
        "authority_sha256": {
            **authorities,
            pilot_path.name: file_sha256(pilot_path),
            freeze_path.name: file_sha256(freeze_path),
            "AP1/CODE/ap1_m1_coupled_moving_split_background.py": code_hash,
        },
        "frozen_thresholds": thresholds,
        "heldout_was_unseen_when_tolerances_frozen": True,
        "heldout_ensemble": ensemble,
        "gates": gates,
        "all_heldout_gates_pass": passed,
        "checkpoint_eligible": passed,
        "new_AP1_M1_background_runs": pilot["new_AP1_M1_background_runs"] + ensemble["new_AP1_M1_background_runs"],
        "bounded_background_candidate_assessed": True,
        "long_background_candidate_assessed": False,
        "bounded_background_tolerances_frozen": True,
        "production_background_tolerances_frozen": False,
        "trajectory_rows_persisted": 0,
        "physical_response_kernel_started": False,
        "seed_released": False,
        "nonpass_stored_or_used_as_seed": False,
        "equations_changed": False,
        "physics_changed": False,
        "parameters_changed": False,
        "existing_gate_thresholds_changed": False,
        "claim_boundary": "bounded coupled held-out background reference only; no present anchor, production trajectory, kernel, observable, curve, fit or significance",
        "next_required": (
            "promote a G29 checkpoint and then extend the same moving-split architecture through staged logarithmic background windows"
            if passed else "stop fail-closed; do not persist, checkpoint or seed this result"
        ),
    }


def build_checkpoint(
    pilot_path: Path, freeze_path: Path, heldout_path: Path, code_path: Path
) -> dict[str, Any]:
    pilot = json.loads(pilot_path.read_text(encoding="utf-8"))
    freeze = json.loads(freeze_path.read_text(encoding="utf-8"))
    heldout = json.loads(heldout_path.read_text(encoding="utf-8"))
    if not (
        pilot.get("all_pilot_gates_pass")
        and freeze.get("all_tolerance_freeze_gates_pass")
        and heldout.get("all_heldout_gates_pass")
    ):
        raise CoupledBackgroundError("full pilot/freeze/heldout PASS chain required")
    expected = {
        pilot_path.name: file_sha256(pilot_path),
        freeze_path.name: file_sha256(freeze_path),
        "AP1/CODE/ap1_m1_coupled_moving_split_background.py": file_sha256(code_path),
    }
    for name, digest in expected.items():
        if heldout["authority_sha256"].get(name) != digest:
            raise CoupledBackgroundError(f"held-out authority mismatch: {name}")
    gates = {
        "pilot_complete": True,
        "tolerances_frozen_before_heldout": True,
        "heldout_full_PASS": True,
        "only_PASS_promoted": True,
        "no_production_background_kernel_seed_or_observable": True,
    }
    return {
        "schema": "apeiron-ap1-m1-coupled-moving-split-background-checkpoint-v1.0",
        "updated_utc": _utc_now(),
        "classification": "M1_BOUNDED_COUPLED_MOVING_SPLIT_BACKGROUND_REFERENCE_PASS_STAGED_LONG_EXTENSION_OPEN_PHYSICAL_KERNEL_BLOCKED",
        "authority_sha256": {
            pilot_path.name: file_sha256(pilot_path),
            freeze_path.name: file_sha256(freeze_path),
            heldout_path.name: file_sha256(heldout_path),
            "AP1/CODE/ap1_m1_coupled_moving_split_background.py": file_sha256(code_path),
        },
        "frozen_thresholds": freeze["thresholds"],
        "heldout_metrics": heldout["heldout_ensemble"]["metrics_for_freeze_or_gate"],
        "heldout_physical_diagnostics": heldout["heldout_ensemble"]["physical_diagnostics"],
        "heldout_source_diagnostics": heldout["heldout_ensemble"]["source_diagnostics"],
        "heldout_precision": heldout["heldout_ensemble"]["precision"],
        "gates": gates,
        "all_checkpoint_gates_pass": True,
        "checkpoint_eligible": True,
        "new_AP1_M1_background_runs": heldout["new_AP1_M1_background_runs"],
        "bounded_background_candidate_assessed": True,
        "long_background_candidate_assessed": False,
        "bounded_background_tolerances_frozen": True,
        "production_background_tolerances_frozen": False,
        "trajectory_rows_persisted": 0,
        "physical_response_kernel_started": False,
        "seed_released": False,
        "nonpass_stored_or_used_as_seed": False,
        "equations_changed": False,
        "physics_changed": False,
        "parameters_changed": False,
        "existing_gate_thresholds_changed": False,
        "AP1_status": "ORANGE",
        "next_required": "extend the same coupled moving-split solver through staged logarithmic N windows; freeze production-background tolerances before any full candidate or physical response kernel",
        "claim_boundary": "bounded coupled background architecture reference only; no present anchor, production trajectory, response kernel, observable, curve, fit or significance",
    }


def checkpoint_markdown(report: dict[str, Any]) -> str:
    metrics = report["heldout_metrics"]
    physical = report["heldout_physical_diagnostics"]
    source = report["heldout_source_diagnostics"]
    return f"""# Apeiron AP1 – G29 bounded coupled moving-split background reference

**Updated UTC:** {report['updated_utc']}  
**Classification:** `{report['classification']}`  
**AP1 status:** ORANGE

The chronology-safe 128-shell state was advanced by a genuinely coupled,
N-parametrized background/mode iteration.  Existing modes remained in the
log-amplitude/Riccati chart.  Every newly resolved mode inherited the same
adiabatic order-0/2/4 state validated by G28 at the moving physical boundary
`q_cut = 0.6 Lambda`; no instantaneous order-0 reset occurred.

## Held-out aggregate

- coupled time-grid background delta: `{metrics['max_background_time_relative_delta']:.17g}`
- coupled time-grid source gross-relative delta: `{metrics['max_source_time_gross_relative_delta']:.17g}`
- entry-time quadrature gross-relative delta: `{metrics['max_entry_quadrature_gross_relative_delta']:.17g}`
- UV-tail quadrature gross-relative delta: `{metrics['max_tail_quadrature_gross_relative_delta']:.17g}`
- UV cutoff gross-relative delta: `{metrics['max_uv_cutoff_gross_relative_delta']:.17g}`
- fixed-point source gross-relative closure: `{metrics['max_fixed_point_source_gross_relative_change']:.17g}`
- Ward residual: `{physical['ward_normalized']:.17g}`
- Friedmann residual: `{physical['max_friedmann_normalized']:.17g}`
- independent Raychaudhuri residual: `{physical['max_raychaudhuri_normalized']:.17g}`
- Wronskian error: `{source['resolved']['max_wronskian_relative_error']:.17g}`

All frozen G29 held-out gates pass.  Four bounded background integrations were
performed across pilot and held-out time-refinement pairs; zero trajectory rows
were promoted to a production curve or seed.

## Boundary

This checkpoint validates only the bounded coupled architecture.  The staged
long-background extension to the present anchor, production tolerances, AP2/AP3,
production curves, observables, and the physical two-time response kernel remain
locked.  No fit or significance is claimed.
"""


def main() -> None:
    parser = argparse.ArgumentParser()
    sub = parser.add_subparsers(dest="command", required=True)
    pilot_parser = sub.add_parser("pilot")
    pilot_parser.add_argument("apeiron_root", type=Path)
    pilot_parser.add_argument("--output", type=Path, required=True)
    freeze_parser = sub.add_parser("freeze")
    freeze_parser.add_argument("pilot", type=Path)
    freeze_parser.add_argument("--output", type=Path, required=True)
    heldout_parser = sub.add_parser("heldout")
    heldout_parser.add_argument("apeiron_root", type=Path)
    heldout_parser.add_argument("pilot", type=Path)
    heldout_parser.add_argument("freeze", type=Path)
    heldout_parser.add_argument("--output", type=Path, required=True)
    checkpoint_parser = sub.add_parser("checkpoint")
    checkpoint_parser.add_argument("pilot", type=Path)
    checkpoint_parser.add_argument("freeze", type=Path)
    checkpoint_parser.add_argument("heldout", type=Path)
    checkpoint_parser.add_argument("--json-output", type=Path, required=True)
    checkpoint_parser.add_argument("--md-output", type=Path, required=True)
    args = parser.parse_args()
    code_path = Path(__file__).resolve()
    if args.command == "pilot":
        report = build_pilot(args.apeiron_root, code_path)
        _write_pass_json(report, args.output, "all_pilot_gates_pass")
    elif args.command == "freeze":
        report = build_freeze(args.pilot, code_path)
        _write_pass_json(report, args.output, "all_tolerance_freeze_gates_pass")
    elif args.command == "heldout":
        report = build_heldout(args.apeiron_root, args.pilot, args.freeze, code_path)
        _write_pass_json(report, args.output, "all_heldout_gates_pass")
    else:
        report = build_checkpoint(args.pilot, args.freeze, args.heldout, code_path)
        _write_pass_json(report, args.json_output, "all_checkpoint_gates_pass")
        args.md_output.write_text(checkpoint_markdown(report), encoding="utf-8")


if __name__ == "__main__":
    main()