"""AP1-R2b isolated chi-mode background closure preflight. This module is a new AP1 code path. It neither imports nor executes the frozen v7.13 solver or Physical Map. It implements the canonical portal, Pauli--Villars moment structure, and an N=ln(a) implicit-midpoint/Yoshida mode propagator. Production quantum sources remain fail-closed until the canonical UV tail and finite renormalization are ported and independently gated. """ from __future__ import annotations from dataclasses import asdict, dataclass from hashlib import sha256 import json import math from pathlib import Path from typing import Mapping import numpy as np from scipy.interpolate import CubicSpline PV_C = np.array([1.0, -3.0, 3.0, -1.0], dtype=float) PV_J = np.array([0.0, 1.0, 2.0, 3.0], dtype=float) LD = np.longdouble PI_LD = LD(np.pi) CANONICAL_SOURCE_HASHES = { "model.py": "db4c0fa9cf4f2013fef299f558ca7b6b509d17c1acbd76f84ece6bd2dd0c5a65", "frw_noncanonical_background_v6_6.py": "f96b13df83e6d6dce8f9bec073e2aab24c23660073cee994d25a327f2b61f7f1", "frw_pv_fixed_mpl_renorm_v6_17.py": "aecd4dc3202604228d1a46c1a4cb439b874e129314e55378e1fbb5b999df482d", "frw_pv_covariant_uv_complete_v6_22.py": "9c54b4095845687b9824dab66da91031436b67a5decd01055ef487b14b4c19b7", "frw_pv_uv_yoshida_v6_31.py": "8fa1b93a2f24377c751f3d9fdffe700394cdb2e00842fc4dc5f35d19ae3943e1", } class ChiClosureError(RuntimeError): """Base class for fail-closed R2b errors.""" class CanonicalSourceError(ChiClosureError): """Canonical source identity does not match the frozen specification.""" class TrajectoryError(ChiClosureError): """A proposed new N-trajectory is outside the mode-core domain.""" class IncompleteRenormalizationError(ChiClosureError): """Production sources were requested before UV/renormalization closure.""" @dataclass(frozen=True) class ChiParameters: Lambda: float = 2.3e-3 phase_scale: float = 0.8549468502440697 displacement_v: float = 0.03272387838931097 portal_lambda: float = 1.6740193604000272e-5 bare_mass: float = 9.817951142008693e-5 muR_over_Lambda: float = 1.5 def validate(self) -> None: values = np.array(list(asdict(self).values()), dtype=float) if not np.all(np.isfinite(values)) or np.any(values <= 0.0): raise ValueError("canonical chi parameters must be finite and positive") @dataclass(frozen=True) class ChiTrajectory: N: np.ndarray H: np.ndarray Hdot: np.ndarray sigma: np.ndarray theta: np.ndarray def validated(self) -> "ChiTrajectory": arrays = { name: np.asarray(getattr(self, name), dtype=float) for name in ("N", "H", "Hdot", "sigma", "theta") } shape = arrays["N"].shape if len(shape) != 1 or shape[0] < 2: raise TrajectoryError("one-dimensional N grid with at least two nodes required") if any(value.shape != shape for value in arrays.values()): raise TrajectoryError("all trajectory arrays must have the same shape") if not all(np.all(np.isfinite(value)) for value in arrays.values()): raise TrajectoryError("trajectory must be finite") if not np.all(np.diff(arrays["N"]) > 0.0): raise TrajectoryError("N grid must be strictly increasing") if not np.all(arrays["H"] > 0.0): raise TrajectoryError("expanding positive-H branch required") return ChiTrajectory(**arrays) 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 audit_canonical_sources(source_dir: Path) -> dict[str, str]: observed = {} for name, expected in CANONICAL_SOURCE_HASHES.items(): path = source_dir / name if not path.is_file(): raise CanonicalSourceError(f"missing canonical source: {name}") value = file_sha256(path) if value != expected: raise CanonicalSourceError(f"canonical source hash mismatch: {name}") observed[name] = value return observed def portal_terms( sigma: np.ndarray | float, theta: np.ndarray | float, p: ChiParameters, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: p.validate() sigma = np.asarray(sigma, dtype=float) theta = np.asarray(theta, dtype=float) phase = 2.0 * theta / p.phase_scale shifted = p.displacement_v + sigma mass2 = p.bare_mass**2 + 2.0 * p.portal_lambda * shifted**2 * np.cos(phase) d_sigma = 4.0 * p.portal_lambda * shifted * np.cos(phase) d_theta = ( -4.0 * p.portal_lambda * shifted**2 * np.sin(phase) / p.phase_scale ) return mass2, d_sigma, d_theta def pv_moment_residuals(max_order: int = 3) -> dict[int, float]: return { order: float(np.sum(PV_C * PV_J**order)) for order in range(max_order + 1) } def physical_shells( Lambda: float, nodes: int = 128, kmax_over_Lambda: float = 0.6, ) -> tuple[np.ndarray, np.ndarray]: if not math.isfinite(Lambda) or Lambda <= 0.0: raise ValueError("positive finite Lambda required") if nodes < 4 or not (0.0 < kmax_over_Lambda <= 1.0): raise ValueError("valid physical shell configuration required") x, weights = np.polynomial.legendre.leggauss(int(nodes)) high = kmax_over_Lambda * Lambda k = 0.5 * high * (x + 1.0) measure = 0.5 * high * weights * k * k / (2.0 * np.pi**2) return k, measure def raw_pv_stress_n( trajectory: ChiTrajectory, k: np.ndarray, weights: np.ndarray, modes: Mapping[str, np.ndarray | float], p: ChiParameters | None = None, ) -> dict[str, np.ndarray | float]: """Resolved physical-support PV stress from the isolated mode core.""" p = ChiParameters() if p is None else p tr = trajectory.validated() k = np.asarray(k, dtype=float) weights = np.asarray(weights, dtype=float) if k.shape != weights.shape or k.ndim != 1: raise ValueError("matching one-dimensional k and weight arrays required") u = np.asarray(modes["u"], dtype=np.complex128) v = np.asarray(modes["v"], dtype=np.complex128) expected = (len(tr.N), 4, len(k)) if u.shape != expected or v.shape != expected: raise ValueError(f"mode history shape must be {expected}") a = np.exp(tr.N) f = u / a[:, None, None] ** 1.5 fdot = (v - 1.5 * tr.H[:, None, None] * u) / a[:, None, None] ** 1.5 mass2 = np.asarray(portal_terms(tr.sigma, tr.theta, p)[0], dtype=float) muR2 = (p.muR_over_Lambda * p.Lambda) ** 2 physical_k2 = k[None, :] ** 2 / a[:, None] ** 2 masses2 = mass2[:, None, None] + PV_J[None, :, None] * muR2 f2 = np.abs(f) ** 2 fdot2 = np.abs(fdot) ** 2 rho_each = np.sum( weights[None, None, :] * 0.5 * (fdot2 + (physical_k2[:, None, :] + masses2) * f2), axis=2, ) pressure_each = np.sum( weights[None, None, :] * ( 0.5 * fdot2 - (physical_k2[:, None, :] / 6.0 + 0.5 * masses2) * f2 ), axis=2, ) chi2_each = np.sum(weights[None, None, :] * f2, axis=2) return { "rho": rho_each @ PV_C, "pressure": pressure_each @ PV_C, "chi2": chi2_each @ PV_C, "gross_rho": np.sum(np.abs(rho_each) * np.abs(PV_C)[None, :], axis=1), } def _binomial_series(alpha: np.longdouble, maximum: int) -> np.ndarray: output = np.empty(maximum + 1, dtype=np.longdouble) output[0] = 1.0 for order in range(1, maximum + 1): output[order] = output[order - 1] * (alpha - (order - 1)) / order return output _BINOMIAL_HALF = _binomial_series(LD("0.5"), 12) _BINOMIAL_MINUS_HALF = _binomial_series(LD("-0.5"), 12) def _pv_mass_moments(masses2: np.ndarray, maximum: int = 12) -> np.ndarray: masses2 = np.asarray(masses2, dtype=np.longdouble) coefficients = PV_C.astype(np.longdouble) return np.stack( [ np.sum(coefficients[None, :] * masses2**order, axis=1) for order in range(maximum + 1) ], axis=1, ) def order0_uv_tail_n( trajectory: ChiTrajectory, p: ChiParameters | None = None, kcut_over_Lambda: float = 0.6, qswitch_over_Lambda: float = 16.0, quadrature_nodes: int = 72, ) -> dict[str, np.ndarray | float]: """Canonical PV adiabatic-order-0 tail from kcut to infinity.""" p = ChiParameters() if p is None else p tr = trajectory.validated() if not (0.0 < kcut_over_Lambda < qswitch_over_Lambda): raise ValueError("tail switch must lie above the physical boundary") a = np.exp(tr.N) q0 = np.asarray(kcut_over_Lambda * p.Lambda / a, dtype=np.longdouble) mass2 = np.asarray(portal_terms(tr.sigma, tr.theta, p)[0], dtype=np.longdouble) muR2 = LD((p.muR_over_Lambda * p.Lambda) ** 2) masses2 = mass2[:, None] + PV_J.astype(np.longdouble)[None, :] * muR2 if np.min(q0[:, None] ** 2 + masses2) <= 0.0: raise TrajectoryError("physical tail boundary entered omega2 <= 0") Q = LD(qswitch_over_Lambda * p.Lambda) if np.any(q0 >= Q): raise TrajectoryError("tail switch below physical boundary") x, weights = np.polynomial.legendre.leggauss(int(quadrature_nodes)) x = x.astype(np.longdouble) weights = weights.astype(np.longdouble) midpoint = (Q + q0) / 2.0 halfwidth = (Q - q0) / 2.0 q = midpoint[:, None] + halfwidth[:, None] * x[None, :] integration_weights = halfwidth[:, None] * weights[None, :] omega = np.sqrt(q[:, None, :] ** 2 + masses2[:, :, None]) coefficients = PV_C.astype(np.longdouble)[None, :, None] signed_omega = np.sum(coefficients * omega, axis=1) signed_inverse_omega = np.sum(coefficients / omega, axis=1) rho_finite = np.sum( integration_weights * q**2 * signed_omega / (4.0 * PI_LD**2), axis=1 ) pressure_finite = np.sum( integration_weights * q**4 * signed_inverse_omega / (12.0 * PI_LD**2), axis=1, ) chi2_finite = np.sum( integration_weights * q**2 * signed_inverse_omega / (4.0 * PI_LD**2), axis=1, ) moments = _pv_mass_moments(masses2, 12) rho_asymptotic = np.zeros(len(tr.N), dtype=np.longdouble) pressure_asymptotic = np.zeros(len(tr.N), dtype=np.longdouble) chi2_asymptotic = np.zeros(len(tr.N), dtype=np.longdouble) for order in range(3, 13): rho_asymptotic += ( _BINOMIAL_HALF[order] * moments[:, order] * Q ** (4 - 2 * order) / (LD(2 * order - 4) * 4.0 * PI_LD**2) ) pressure_asymptotic += ( _BINOMIAL_MINUS_HALF[order] * moments[:, order] * Q ** (4 - 2 * order) / (LD(2 * order - 4) * 12.0 * PI_LD**2) ) chi2_asymptotic += ( _BINOMIAL_MINUS_HALF[order] * moments[:, order] * Q ** (2 - 2 * order) / (LD(2 * order - 2) * 4.0 * PI_LD**2) ) last_order = 12 remainder_proxy = np.max( np.abs( _BINOMIAL_HALF[last_order] * moments[:, last_order] * Q ** (4 - 2 * last_order) / (LD(2 * last_order - 4) * 4.0 * PI_LD**2) ) ) return { "rho": np.asarray(rho_finite + rho_asymptotic, dtype=float), "pressure": np.asarray(pressure_finite + pressure_asymptotic, dtype=float), "chi2": np.asarray(chi2_finite + chi2_asymptotic, dtype=float), "min_tail_omega2": float(np.min(q0[:, None] ** 2 + masses2)), "asymptotic_last_term_abs": float(remainder_proxy), "Q_over_Lambda": float(qswitch_over_Lambda), } def cosmic_derivative(values: np.ndarray, trajectory: ChiTrajectory) -> np.ndarray: """Differentiate along a prescribed trajectory using d/dt = H d/dN.""" tr = trajectory.validated() values = np.asarray(values, dtype=float) if values.shape[0] != len(tr.N): raise ValueError("time axis must match trajectory") if len(tr.N) < 5: raise TrajectoryError("at least five N nodes required for curvature derivatives") derivative_N = CubicSpline(tr.N, values, axis=0)(tr.N, 1) expansion = tr.H.reshape((len(tr.N),) + (1,) * (values.ndim - 1)) return expansion * derivative_N def _logarithmic_shells( low: float, high: float, nodes_per_octave: int = 20, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: if not high > low > 0.0 or nodes_per_octave < 4: raise ValueError("valid logarithmic shell bounds required") edges = [float(low)] while edges[-1] < high * (1.0 - 1.0e-14): edges.append(min(edges[-1] * 2.0, high)) x, weights = np.polynomial.legendre.leggauss(int(nodes_per_octave)) shells = [] measures = [] segments = [] for segment, (left, right) in enumerate(zip(edges[:-1], edges[1:])): k = 0.5 * (right - left) * x + 0.5 * (left + right) measure = 0.5 * (right - left) * weights * k * k / (2.0 * np.pi**2) shells.extend(k) measures.extend(measure) segments.extend([segment] * len(k)) return np.asarray(shells), np.asarray(measures), np.asarray(segments, dtype=int) def curvature_uv_tail_n( trajectory: ChiTrajectory, p: ChiParameters | None = None, kcut_over_Lambda: float = 0.6, K_over_Lambda: float = 64.0, nodes_per_octave: int = 20, ) -> dict[str, np.ndarray | float | list[dict]]: """Canonical local adiabatic curvature orders 2+4 on an N grid.""" p = ChiParameters() if p is None else p tr = trajectory.validated() if len(tr.N) < 5: raise TrajectoryError("curvature tail requires at least five N nodes") a = np.exp(tr.N) kcut = kcut_over_Lambda * p.Lambda K = K_over_Lambda * p.Lambda * float(np.max(a)) k, weights, segments = _logarithmic_shells(kcut, K, nodes_per_octave) physical_k2 = k[None, :] ** 2 / a[:, None] ** 2 mass2 = np.asarray(portal_terms(tr.sigma, tr.theta, p)[0], dtype=float) muR2 = (p.muR_over_Lambda * p.Lambda) ** 2 rho = np.zeros(len(tr.N)) pressure = np.zeros(len(tr.N)) chi2 = np.zeros(len(tr.N)) segmented_rho = np.zeros((int(np.max(segments)) + 1, len(tr.N))) hierarchy = [] for sector in range(4): masses2 = mass2[:, None] + PV_J[sector] * muR2 omega2 = physical_k2 + masses2 if np.min(omega2) <= 0.0: raise TrajectoryError(f"curvature tail sector {sector} entered omega2 <= 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] ** 2 - 1.5 * tr.Hdot[:, 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) prefactor = 1.0 / (2.0 * a[:, None] ** 3) A0 = prefactor / omega A2 = -prefactor * W2 / omega**2 A4 = prefactor * (W2**2 / omega**3 - W4 / omega**2) B1 = 1.5 * tr.H[:, None] + 0.5 * logarithmic0 B3 = 0.5 * logarithmic2 C0 = omega2 C2 = 2.0 * omega * W2 + B1**2 C4 = W2**2 + 2.0 * omega * W4 + 2.0 * B1 * B3 D2 = A0 * C2 + A2 * C0 D4 = A0 * C4 + A2 * C2 + A4 * C0 rho24 = 0.5 * (D2 + omega2 * A2 + D4 + omega2 * A4) pressure_weight = physical_k2 / 6.0 + 0.5 * masses2 pressure24 = 0.5 * (D2 + D4) - pressure_weight * (A2 + A4) chi24 = A2 + A4 coefficient = PV_C[sector] rho += coefficient * np.sum(weights[None, :] * rho24, axis=1) pressure += coefficient * np.sum(weights[None, :] * pressure24, axis=1) chi2 += coefficient * np.sum(weights[None, :] * chi24, axis=1) for segment in range(segmented_rho.shape[0]): mask = segments == segment segmented_rho[segment] += coefficient * np.sum( weights[mask][None, :] * rho24[:, mask], axis=1 ) hierarchy.append( { "sector": sector, "min_tail_omega2": float(np.min(omega2)), "max_abs_W2_over_W0": float(np.max(np.abs(W2 / omega))), "max_abs_W4_over_W0": float(np.max(np.abs(W4 / omega))), } ) last_octave = float(np.max(np.abs(segmented_rho[-1]))) total_scale = max(float(np.max(np.abs(rho))), 1.0e-300) return { "rho": rho, "pressure": pressure, "chi2": chi2, "hierarchy": hierarchy, "last_octave_rho_over_total": last_octave / total_scale, "K_comoving_over_Lambda": K / p.Lambda, "nodes": int(len(k)), } def completed_quantum_sources_n( trajectory: ChiTrajectory, p: ChiParameters | None = None, physical_nodes: int = 48, curvature_nodes_per_octave: int = 12, ) -> dict[str, np.ndarray | float | dict | list]: """Assemble the isolated renormalized source on a prescribed trajectory.""" p = ChiParameters() if p is None else p tr = trajectory.validated() k, weights = physical_shells(p.Lambda, nodes=physical_nodes) modes = propagate_modes_n(tr, k, p.Lambda, p) resolved = raw_pv_stress_n(tr, k, weights, modes, p) tail0 = order0_uv_tail_n(tr, p) tail24 = curvature_uv_tail_n( tr, p, nodes_per_octave=curvature_nodes_per_octave ) pre_rho = resolved["rho"] + tail0["rho"] + tail24["rho"] pre_pressure = ( resolved["pressure"] + tail0["pressure"] + tail24["pressure"] ) pre_chi2 = resolved["chi2"] + tail0["chi2"] + tail24["chi2"] constants = frozen_renormalization_constants(p) C = constants["C_chi2"] A = constants["A_g"] B = constants["B_G_rhs"] mass2 = np.asarray(portal_terms(tr.sigma, tr.theta, p)[0], dtype=float) rho = pre_rho - 0.5 * C * mass2 + A + 3.0 * B * tr.H**2 pressure = ( pre_pressure + 0.5 * C * mass2 - A - B * (2.0 * tr.Hdot + 3.0 * tr.H**2) ) chi2 = pre_chi2 - C mass2_dot = cosmic_derivative(mass2, tr) rho_dot = cosmic_derivative(rho, tr) ward = rho_dot + 3.0 * tr.H * (rho + pressure) - 0.5 * mass2_dot * chi2 ward_scale = max( float(np.max(np.abs(rho_dot))), float(np.max(np.abs(3.0 * tr.H * (rho + pressure)))), float(np.max(np.abs(0.5 * mass2_dot * chi2))), 1.0e-300, ) return { "rho": rho, "pressure": pressure, "chi2": chi2, "ward_residual": ward, "ward_normalized": float(np.max(np.abs(ward)) / ward_scale), "wronskian_relative_error": modes["wronskian_relative_error"], "tail0": tail0, "tail24": tail24, "renormalization": constants, } def _benchmark_trajectory(kind: str, nodes: int) -> ChiTrajectory: if kind == "de_sitter": N = np.linspace(0.0, 0.05, nodes) H = np.full_like(N, 1.0e-4) return ChiTrajectory( N=N, H=H, Hdot=np.zeros_like(N), sigma=np.full_like(N, 0.02), theta=np.full_like(N, -0.69), ) if kind == "smooth_portal": N = np.linspace(0.0, 0.02, nodes) H = 1.0e-4 * np.exp(-0.1 * N) return ChiTrajectory( N=N, H=H, Hdot=-0.1 * H**2, sigma=0.02 + 0.001 * N, theta=-0.69 + 0.002 * N, ) raise ValueError("unknown analytic benchmark") def analytic_canary_report(p: ChiParameters | None = None) -> dict: """Three-level de Sitter convergence plus an independent portal canary.""" p = ChiParameters() if p is None else p levels = [(33, 16, 6), (65, 24, 8), (129, 32, 10)] de_sitter = [] for N_nodes, physical_nodes, curvature_nodes in levels: result = completed_quantum_sources_n( _benchmark_trajectory("de_sitter", N_nodes), p, physical_nodes, curvature_nodes, ) de_sitter.append( { "N_nodes": N_nodes, "physical_nodes": physical_nodes, "curvature_nodes_per_octave": curvature_nodes, "ward_normalized": result["ward_normalized"], "wronskian_relative_error": result[ "wronskian_relative_error" ], "endpoint": { key: float(result[key][-1]) for key in ("rho", "pressure", "chi2") }, } ) endpoint_relative = {} for key in ("rho", "pressure", "chi2"): coarse = de_sitter[0]["endpoint"][key] fine = de_sitter[-1]["endpoint"][key] endpoint_relative[key] = abs(coarse - fine) / max(abs(fine), 1.0e-300) portal_ward = [] portal_finite = True for N_nodes, physical_nodes, curvature_nodes in levels: result = completed_quantum_sources_n( _benchmark_trajectory("smooth_portal", N_nodes), p, physical_nodes, curvature_nodes, ) portal_ward.append(result["ward_normalized"]) portal_finite = portal_finite and all( np.all(np.isfinite(result[key])) for key in ("rho", "pressure", "chi2", "ward_residual") ) de_sitter_wards = [item["ward_normalized"] for item in de_sitter] max_wronskian = max(item["wronskian_relative_error"] for item in de_sitter) gates = { "de_sitter_Ward_refines_monotonically": bool( de_sitter_wards[2] < de_sitter_wards[1] < de_sitter_wards[0] ), "de_sitter_fine_Ward_below_2e_6": bool(de_sitter_wards[-1] < 2.0e-6), "de_sitter_endpoint_convergence_below_2e_5": bool( max(endpoint_relative.values()) < 2.0e-5 ), "de_sitter_wronskian_below_1e_12": bool(max_wronskian < 1.0e-12), "smooth_portal_finite": bool(portal_finite), "smooth_portal_Ward_refines_monotonically": bool( portal_ward[2] < portal_ward[1] < portal_ward[0] ), } return { "de_sitter_levels": de_sitter, "de_sitter_endpoint_coarse_to_fine_relative": endpoint_relative, "smooth_portal_Ward_levels": portal_ward, "gates": gates, "pass": bool(all(gates.values())), "role": "analytic code canaries only; not an Apeiron M1 trajectory", } def flat_pv_exact(mass2: float, muR2: float) -> tuple[float, float, float]: """Canonical closed infinite-PV Minkowski stress and variance.""" masses2 = np.asarray(mass2 + PV_J * muR2, dtype=np.longdouble) coefficients = PV_C.astype(np.longdouble) if np.min(masses2) <= 0.0: raise ValueError("flat renormalization masses must be positive") pi = np.longdouble(np.pi) rho = np.sum(coefficients * masses2**2 * np.log(masses2)) / (64.0 * pi**2) chi2 = np.sum(coefficients * masses2 * np.log(masses2)) / (16.0 * pi**2) return float(rho), float(-rho), float(chi2) def newton_loop_response( mass2: float, muR2: float, Lambda: float, Qratio: float = 1024.0, nodes: int = 48, ) -> dict[str, float]: """Canonical H->0 order-2 PV response fixing the Einstein term.""" masses2 = np.asarray(mass2 + PV_J * muR2, dtype=np.longdouble) coefficients = PV_C.astype(np.longdouble) xg, wg = np.polynomial.legendre.leggauss(int(nodes)) xg = xg.astype(np.longdouble) wg = wg.astype(np.longdouble) scale = np.longdouble(Lambda) edges = [np.longdouble(0.0), scale / 16.0] q_edge = scale / 16.0 Q = np.longdouble(Qratio) * scale while q_edge < Q: q_edge = min(q_edge * 2.0, Q) edges.append(q_edge) rho_response = np.longdouble(0.0) pressure_response = np.longdouble(0.0) pi = np.longdouble(np.pi) for low, high in zip(edges[:-1], edges[1:]): q = 0.5 * (high - low) * xg + 0.5 * (high + low) weights = 0.5 * (high - low) * wg omega = np.sqrt(q[None, :] ** 2 + masses2[:, None]) fraction = q[None, :] ** 2 / (q[None, :] ** 2 + masses2[:, None]) b = 1.5 - 0.5 * fraction kernel_rho = np.sum(coefficients[:, None] * b**2 / (4.0 * omega), axis=0) F = -2.25 - fraction + 1.25 * fraction**2 kernel_pressure = np.sum( coefficients[:, None] * (b**2 / 4.0 + F * (0.25 - fraction / 12.0)) / omega, axis=0, ) measure = weights * q**2 / (2.0 * pi**2) rho_response += np.sum(measure * kernel_rho) pressure_response += np.sum(measure * kernel_pressure) loop_B = rho_response / 3.0 return { "rho_over_H2": float(rho_response), "pressure_over_H2": float(pressure_response), "enthalpy_over_H2": float(rho_response + pressure_response), "B_loop": float(loop_B), "B_rhs_counterterm": float(-loop_B), "Q_over_Lambda": float(Qratio), } def frozen_renormalization_constants( p: ChiParameters | None = None, sigma_reference: float = 0.015257914139527498, theta_reference: float = -0.69, ) -> dict: """Trajectory-independent finite constants at the canonical flat slice.""" p = ChiParameters() if p is None else p p.validate() mass2_reference = float(portal_terms(sigma_reference, theta_reference, p)[0]) muR2 = (p.muR_over_Lambda * p.Lambda) ** 2 raw_rho, raw_pressure, raw_chi2 = flat_pv_exact(mass2_reference, muR2) C_chi2 = raw_chi2 A_g = -(raw_rho - 0.5 * C_chi2 * mass2_reference) response = newton_loop_response(mass2_reference, muR2, p.Lambda) return { "mass2_reference": mass2_reference, "sigma_reference": sigma_reference, "theta_reference": theta_reference, "rho_flat_raw": raw_rho, "pressure_flat_raw": raw_pressure, "chi2_flat_raw": raw_chi2, "C_chi2": C_chi2, "A_g": A_g, "B_G_rhs": response["B_rhs_counterterm"], "B_response": response, "R2": 0.0, } def flat_calibration_audit(p: ChiParameters | None = None) -> dict: constants = frozen_renormalization_constants(p) mass2 = constants["mass2_reference"] C = constants["C_chi2"] A = constants["A_g"] rho = constants["rho_flat_raw"] - 0.5 * C * mass2 + A pressure = constants["pressure_flat_raw"] + 0.5 * C * mass2 - A chi2 = constants["chi2_flat_raw"] - C response = constants["B_response"] response_scale = max( abs(response["rho_over_H2"]), abs(response["pressure_over_H2"]), 1.0e-300, ) return { "rho_ren_flat": rho, "pressure_ren_flat": pressure, "chi2_ren_flat": chi2, "flat_enthalpy": rho + pressure, "curvature_enthalpy_relative": abs(response["enthalpy_over_H2"]) / response_scale, "constants": constants, "gates": { "minkowski_rho_zero": abs(rho) < 1.0e-24, "minkowski_pressure_zero": abs(pressure) < 1.0e-24, "minkowski_chi2_zero": abs(chi2) < 1.0e-20, "curvature_response_covariant": ( abs(response["enthalpy_over_H2"]) / response_scale < 1.0e-9 ), "newton_shift_small": abs(response["B_rhs_counterterm"]) < 1.0e-6, }, } def _interp_mid(x: np.ndarray, y: np.ndarray, value: float) -> float: return float(np.interp(value, x, y)) def _omega2_at( N_value: float, trajectory: ChiTrajectory, k: np.ndarray, regulator_mass2: float, p: ChiParameters, ) -> np.ndarray: sigma = _interp_mid(trajectory.N, trajectory.sigma, N_value) theta = _interp_mid(trajectory.N, trajectory.theta, N_value) H = _interp_mid(trajectory.N, trajectory.H, N_value) Hdot = _interp_mid(trajectory.N, trajectory.Hdot, N_value) mass2 = float(portal_terms(sigma, theta, p)[0]) + regulator_mass2 physical = k * k * math.exp(-2.0 * N_value) + mass2 return physical - 2.25 * H * H - 1.5 * Hdot def propagate_modes_n( trajectory: ChiTrajectory, k: np.ndarray, Lambda: float, p: ChiParameters | None = None, ) -> dict[str, np.ndarray | float]: """Propagate u=a^(3/2)chi on a prescribed expanding N trajectory. The first-order system is du/dN=v/H and dv/dN=-Omega^2 u/H. Each Yoshida substep uses an implicit midpoint update at its physical midpoint. This is a mode-core test only; it does not construct renormalized sources. """ p = ChiParameters() if p is None else p p.validate() tr = trajectory.validated() k = np.asarray(k, dtype=float) if k.ndim != 1 or len(k) == 0 or not np.all(np.isfinite(k)) or np.any(k < 0.0): raise ValueError("finite non-negative one-dimensional k grid required") if not math.isfinite(Lambda) or Lambda <= 0.0: raise ValueError("positive finite Lambda required") mass2_0 = float(portal_terms(tr.sigma[0], tr.theta[0], p)[0]) muR2 = (p.muR_over_Lambda * Lambda) ** 2 omega2_0 = ( k[None, :] ** 2 * math.exp(-2.0 * tr.N[0]) + mass2_0 + PV_J[:, None] * muR2 ) if np.any(omega2_0 <= 0.0): raise TrajectoryError("non-positive initial canonical frequency") omega0 = np.sqrt(omega2_0) u = (1.0 / np.sqrt(2.0 * omega0)).astype(np.complex128) v = (-1j * omega0 * u).astype(np.complex128) u_history = np.empty((len(tr.N), 4, len(k)), dtype=np.complex128) v_history = np.empty_like(u_history) u_history[0] = u v_history[0] = v cbrt2 = 2.0 ** (1.0 / 3.0) gamma1 = 1.0 / (2.0 - cbrt2) gamma0 = -cbrt2 / (2.0 - cbrt2) min_heavy_physical_omega2 = math.inf for index in range(len(tr.N) - 1): dN = float(tr.N[index + 1] - tr.N[index]) cursor = float(tr.N[index]) for gamma in (gamma1, gamma0, gamma1): substep = gamma * dN midpoint = cursor + 0.5 * substep Hmid = _interp_mid(tr.N, tr.H, midpoint) alpha = substep / (2.0 * Hmid) for sector in range(4): regulator = PV_J[sector] * muR2 Omega2 = _omega2_at(midpoint, tr, k, regulator, p) physical = ( k * k * math.exp(-2.0 * midpoint) + float(portal_terms( _interp_mid(tr.N, tr.sigma, midpoint), _interp_mid(tr.N, tr.theta, midpoint), p, )[0]) + regulator ) if sector > 0: min_heavy_physical_omega2 = min( min_heavy_physical_omega2, float(np.min(physical)) ) denominator = 1.0 + alpha * alpha * Omega2 if np.any(np.abs(denominator) < 1.0e-12): raise TrajectoryError("implicit midpoint denominator approached zero") factor = 1.0 - alpha * alpha * Omega2 old_u = u[sector].copy() old_v = v[sector].copy() u[sector] = (factor * old_u + 2.0 * alpha * old_v) / denominator v[sector] = ( -2.0 * alpha * Omega2 * old_u + factor * old_v ) / denominator cursor += substep u_history[index + 1] = u v_history[index + 1] = v wronskian = -2.0 * np.imag(np.conj(u_history) * v_history) wronskian_relative_error = float( np.max(np.abs(wronskian / wronskian[0:1] - 1.0)) ) return { "u": u_history, "v": v_history, "wronskian": wronskian, "wronskian_relative_error": wronskian_relative_error, "min_heavy_physical_omega2": float(min_heavy_physical_omega2), } def require_production_renormalization( uv_tail_provider: object | None, finite_counterterm_provider: object | None, ) -> None: missing = [] if uv_tail_provider is None: missing.append("canonical order-0 and adiabatic order-2/4 UV tail") if finite_counterterm_provider is None: missing.append("trajectory-independent frozen finite counterterms") if missing: raise IncompleteRenormalizationError("; ".join(missing)) def build_preflight_report(source_dir: Path) -> dict: hashes = audit_canonical_sources(source_dir) p = ChiParameters() # Synthetic prescribed trajectory: a mode-core reference, not Apeiron data. N = np.linspace(0.0, 2.0e-4, 17) H = np.full_like(N, 1.0e-4) Hdot = np.zeros_like(N) sigma = np.full_like(N, 0.02) theta = np.full_like(N, -0.69) tr = ChiTrajectory(N=N, H=H, Hdot=Hdot, sigma=sigma, theta=theta) modes = propagate_modes_n(tr, np.array([0.0, 2.0e-4]), Lambda=1.0e-3, p=p) k_shell, shell_weights = physical_shells(p.Lambda, nodes=24) shell_modes = propagate_modes_n(tr, k_shell, Lambda=p.Lambda, p=p) resolved = raw_pv_stress_n(tr, k_shell, shell_weights, shell_modes, p) tail0 = order0_uv_tail_n(tr, p, quadrature_nodes=48) N_curvature = np.linspace(0.0, 0.05, 33) curvature_reference = ChiTrajectory( N=N_curvature, H=np.full_like(N_curvature, 1.0e-4), Hdot=np.zeros_like(N_curvature), sigma=np.full_like(N_curvature, 0.02), theta=np.full_like(N_curvature, -0.69), ) tail24 = curvature_uv_tail_n(curvature_reference, p, nodes_per_octave=8) closure_error = None try: require_production_renormalization(None, frozen_renormalization_constants) except IncompleteRenormalizationError as exc: closure_error = str(exc) moments = pv_moment_residuals() flat_audit = flat_calibration_audit(p) canaries = analytic_canary_report(p) return { "schema": "apeiron-ap1-r2b-chi-mode-core-preflight-v1.0", "status": "LOCAL_CHI_CLOSURE_ANALYTIC_CANARIES_PASS_PRODUCTION_TRAJECTORY_NOT_RUN", "canonical_sources": hashes, "parameters": asdict(p), "pv_moment_residuals": {str(key): value for key, value in moments.items()}, "frozen_flat_renormalization": flat_audit, "N_mode_core": { "parameterization": "N_equals_ln_a", "old_solver_or_physical_map_called": False, "synthetic_reference_only": True, "finite": bool( np.all(np.isfinite(modes["u"])) and np.all(np.isfinite(modes["v"])) ), "wronskian_relative_error": modes["wronskian_relative_error"], "min_heavy_physical_omega2": modes["min_heavy_physical_omega2"], }, "resolved_and_order0_source_preflight": { "finite": bool( all(np.all(np.isfinite(value)) for value in resolved.values()) and all( np.all(np.isfinite(value)) for key, value in tail0.items() if isinstance(value, np.ndarray) ) ), "order0_tail_min_omega2": tail0["min_tail_omega2"], "order0_asymptotic_last_term_abs": tail0[ "asymptotic_last_term_abs" ], }, "curvature_tail_preflight": { "finite": bool( all( np.all(np.isfinite(tail24[key])) for key in ("rho", "pressure", "chi2") ) ), "last_octave_rho_over_total": tail24[ "last_octave_rho_over_total" ], "hierarchy": tail24["hierarchy"], "synthetic_reference_only": True, }, "analytic_canaries": canaries, "production_background_sources": { "runnable": False, "finite_counterterms_ported": bool(all(flat_audit["gates"].values())), "order0_uv_tail_ported": True, "curvature_order2_order4_code_ported": True, "analytic_canaries_pass": canaries["pass"], "blocking_closure": ( "candidate-trajectory Ward and multi-resolution convergence gates" ), "rho_q_released": False, "p_q_released": False, "chi2_q_released": False, }, "claim_boundary": "no Apeiron trajectory, observable curve, fit, or empirical claim", } def main() -> None: import argparse parser = argparse.ArgumentParser() parser.add_argument("canonical_source_dir", type=Path) parser.add_argument("--output", type=Path) args = parser.parse_args() rendered = json.dumps(build_preflight_report(args.canonical_source_dir), indent=2) + "\n" if args.output is None: print(rendered, end="") else: args.output.write_text(rendered, encoding="utf-8") if __name__ == "__main__": main()