#!/usr/bin/env python3
"""
PROVENANCE: PROOF

SCOPE NOTE, 2026-10-10 (comp 201).  "a_* = d_0 at k T = 83 eV" below is evaluated at the
Casimir BOUND d_0 = 2.3 nm.  comp 201 shows the book's canonical normalisation fixes
d_0 = sqrt(2) hbar c / m_h = 2.2 am, nine orders lower, where a_* passes through d_0 at
k T ~ 40 GeV, the electroweak epoch.  No number here changes; the label does.

Computation 196 -- the truncation spacing a_* is comoving, and under the conformal
identification it is the thermal wavelength of the state's radiation
==================================================================================
WHAT THIS DOES.  comp 182 left the conformal identification of Pi_mn (the first item
of sec:open-structural) waiting on "one number": the physical spacing a_* at which
A1's idealised limit is truncated, with Pi_mn ~ 1/a_*^4, and comp 167 had shown that
neither d_0 nor l_* = hbar c / M_* is that spacing.  The item has carried a_* as a
constant of the theory to be derived.  This file shows it is not a constant, and says
what it is.

THE ARGUMENT HAS TWO STEPS, AND ONLY THE SECOND USES THE IDENTIFICATION.

  Step 1 (exact, no identification needed).  eq:trace-conservation says
  nabla^m Pi_mn = -(1/4) nabla_n T and comp 158 showed the trace vanishes on balanced
  clouds, which homogeneous radiation is, so Pi_mn is conserved exactly there (comp 199
  says what happens out of balance).  On the Friedmann-Robertson-Walker geometry eq:frw-metric that
  Chapter 16 derives from the Cosmological Principle, a homogeneous isotropic
  traceless tensor is Pi^0_0 = rho, Pi^i_j = -(rho/3) delta^i_j, and conservation
  reads rho-dot + 4 H rho = 0, so rho = C / a(t)^4.  A tensor of magnitude 1/a_*^4
  with a_* a FIXED physical length is the case rho = const, and its covariant
  divergence is 4 H rho, which is not zero in an expanding geometry.  So a_* cannot be
  a fixed length: it scales with a(t).  The truncation spacing is COMOVING.  That is
  a theorem of the conservation law alone, checked symbolically below for every
  spatial curvature k, and the exponent 4 is the tracelessness (dust, w = 0, gives
  rho ~ a^-3, the control).

  Step 2 (under the identification).  sec:einstein identifies Pi_mn with the
  conformal part of T^(fields)_mn.  On FRW the stress tensor of a perfect fluid has
  trace rho - 3p, which vanishes iff w = 1/3: the traceless homogeneous component of
  T^(fields) is the radiation, and nothing else.  So under the identification
  Pi^0_0 = rho_r, the radiation energy density of the state, and

        a_*  =  (hbar c / rho_r)^(1/4).

  Today, with T_CMB = 2.7255 K and three neutrino species (N_eff = 3.044),
  rho_r = 7.05e-14 J/m^3 and a_* = 0.82 mm.  That is the thermal wavelength of the
  radiation bath: hbar c / k T_CMB = 0.84 mm, and the two differ by the factor
  (15 / (pi^2 f))^(1/4) = 0.97 with f = 1.69 the neutrino multiplier.  It scales as
  a(t), which is what Step 1 demands: a_* ~ 1/T.

WHAT THIS DOES TO THE ITEM.  The "one number" is not a constant of the theory and is
not owed by it.  It is a state variable, the magnitude of Pi_mn being the radiation
content of the state exactly as the magnitude of T_mn in general relativity is
supplied by the state and not by the field equations.  Under the identification,
Pi_mn's magnitude is fixed by what the identification says it is, and the second gap
of the item is absorbed into the first.  What is left of the item is the first gap
alone: the identification holds on the semidefinite metric-traceless cone that
eq:Pimunu reaches, and not outside it (comp 164).  The item stays open on that gap
and is narrower by one number.

Three things this does NOT do.  It does not compute rho_r: the density is read from
the state, so comp 167's conclusion that no absolute scale is computed, and with it
the absence of a vacuum-energy catastrophe, is untouched.  It does not establish the
identification; it reads a_* through it.  And it does not make a_* equal to d_0 or
l_* at any privileged moment: a_* passes through d_0 = 2.3 nm at k T = 83 eV and
through l_* near k T = M_*, which is comp 167's negative result restated as a
function of time.

CHECKS.  (1) Conservation of a homogeneous traceless tensor on FRW with curvature k
as a symbol: the only nonzero component of nabla_m Pi^m_n is rho-dot + 4 H rho,
exactly, for every k.  (2) rho = C a^-4 solves it; rho = const leaves 4 H rho.
(3) Dust control: w = 0 gives rho-dot + 3 H rho.  (4) Trace of a perfect fluid is
zero iff w = 1/3.  (5) The numbers: rho_gamma from the Planck law, the neutrino
multiplier, a_* today, the ratio to hbar c / k T, and the temperatures at which a_*
equals d_0 and l_*.  (6) Tamper: a fixed a_* is rho = const and fails (2); a
non-traceless Pi (w = 0) is conserved with a^-3 and would contradict comp 158.
"""

import math
import sympy as sp


# ---------------------------------------------------------------------------
# Part 1: conservation of a homogeneous isotropic traceless tensor on FRW
# ---------------------------------------------------------------------------
def frw_divergence(w):
    """Return nabla_m Pi^m_n for Pi = diag(rho, -w rho, -w rho, -w rho) (mixed
    components) on the FRW metric with curvature k, in coordinates (t, r, th, ph).
    The result is the list of the four covariant components."""
    t, r, th, ph, k = sp.symbols('t r theta phi k', real=True)
    a = sp.Function('a')(t)
    rho = sp.Function('rho')(t)
    x = [t, r, th, ph]
    g = sp.diag(1, -a ** 2 / (1 - k * r ** 2), -a ** 2 * r ** 2, -a ** 2 * r ** 2 * sp.sin(th) ** 2)
    ginv = g.inv()
    # Christoffel symbols Gamma^l_{mn}
    Gam = [[[sp.simplify(sum(ginv[l, s] * (sp.diff(g[s, m], x[n]) + sp.diff(g[s, n], x[m]) - sp.diff(g[m, n], x[s])) for s in range(4)) / 2)
             for n in range(4)] for m in range(4)] for l in range(4)]
    # mixed tensor Pi^m_n
    P = sp.diag(rho, -w * rho, -w * rho, -w * rho)
    div = []
    for n in range(4):
        expr = 0
        for m in range(4):
            expr += sp.diff(P[m, n], x[m])
            for l in range(4):
                expr += Gam[m][m][l] * P[l, n] - Gam[l][m][n] * P[m, l]
        div.append(sp.simplify(expr))
    return div, (a, rho, t)


def part1():
    print("  PART 1.  nabla_m Pi^m_n on FRW with curvature k, for a homogeneous isotropic tensor")
    ok = True
    div, (a, rho, t) = frw_divergence(sp.Rational(1, 3))
    H = sp.diff(a, t) / a
    target = sp.diff(rho, t) + 4 * H * rho
    c0 = sp.simplify(div[0] - target) == 0
    cs = all(sp.simplify(d) == 0 for d in div[1:])
    print(f"      traceless (w = 1/3): time component = rho-dot + 4 H rho, exactly, k symbolic : {c0}")
    print(f"      spatial components vanish identically                                     : {cs}")
    ok &= c0 and cs
    # solutions
    C = sp.symbols('C', positive=True)
    sol = target.subs(rho, C / a ** 4).doit()
    c1 = sp.simplify(sol) == 0
    print(f"      rho = C a^-4 solves it                                                    : {c1}")
    const = sp.simplify(target.subs(rho, C).doit())
    c2 = sp.simplify(const - 4 * H * C) == 0
    print(f"      rho = const (a FIXED a_*) leaves the residual 4 H rho, not zero             : {c2}  [{const}]")
    ok &= c1 and c2
    # dust control
    divd, (a, rho, t) = frw_divergence(0)
    H = sp.diff(a, t) / a
    c3 = sp.simplify(divd[0] - (sp.diff(rho, t) + 3 * H * rho)) == 0
    print(f"      control, dust (w = 0): rho-dot + 3 H rho, so the exponent 4 is the trace     : {c3}")
    ok &= c3
    # trace of a perfect fluid
    w = sp.symbols('w')
    rho_s = sp.symbols('rho', positive=True)
    trace = rho_s - 3 * w * rho_s
    c4 = sp.solve(sp.Eq(trace, 0), w) == [sp.Rational(1, 3)]
    print(f"      trace rho - 3p of a perfect fluid vanishes iff w = 1/3                      : {c4}")
    ok &= c4
    print("      => a tensor of magnitude 1/a_*^4 is conserved iff a_* scales with a(t): the spacing is comoving.")
    return ok


# ---------------------------------------------------------------------------
# Part 2: the numbers, under the identification Pi^0_0 = rho_r
# ---------------------------------------------------------------------------
HBAR = 1.054571817e-34      # J s
C_LIGHT = 2.99792458e8      # m/s
KB = 1.380649e-23           # J/K
EV = 1.602176634e-19        # J
T_CMB = 2.7255              # K
N_EFF = 3.044
D0 = 2.3e-9                 # m, the Casimir bound of Chapter 15
MSTAR_GEV = 4 * math.pi * 125.25 * math.sqrt(2.0 / 3.0)   # 1285 GeV, eq:mstar-pred
HBARC = HBAR * C_LIGHT      # J m


def rho_gamma(T):
    """Planck law energy density of photons at temperature T, J/m^3."""
    return (math.pi ** 2 / 15) * (KB * T) ** 4 / HBARC ** 3


def neutrino_multiplier(n_eff=N_EFF):
    return 1 + (7.0 / 8) * (4.0 / 11) ** (4.0 / 3) * n_eff


def a_star(rho):
    return (HBARC / rho) ** 0.25


def part2():
    print("  PART 2.  a_* today, under Pi^0_0 = rho_r")
    ok = True
    rg = rho_gamma(T_CMB)
    f = neutrino_multiplier()
    rr = rg * f
    print(f"      rho_gamma(2.7255 K) = {rg:.3e} J/m^3 = {rg / EV / 1e6:.3f} eV/cm^3;  neutrino multiplier f = {f:.3f};  rho_r = {rr:.3e} J/m^3")
    ok &= 4.1e-14 < rg < 4.3e-14 and 1.68 < f < 1.70
    ast = a_star(rr)
    lam_T = HBARC / (KB * T_CMB)
    n_gamma = (2 * 1.2020569 / math.pi ** 2) * (KB * T_CMB / HBARC) ** 3
    print(f"      a_* = (hbar c / rho_r)^(1/4) = {ast * 1e3:.3f} mm")
    print(f"      hbar c / k T_CMB = {lam_T * 1e3:.3f} mm;  a_* / (hbar c / k T) = {ast / lam_T:.4f} = (15/(pi^2 f))^(1/4) = {(15 / (math.pi ** 2 * f)) ** 0.25:.4f}")
    print(f"      n_gamma^(-1/3) = {n_gamma ** (-1 / 3) * 1e3:.3f} mm  (n_gamma = {n_gamma / 1e6:.1f} per cm^3)")
    ok &= 0.7e-3 < ast < 0.9e-3 and abs(ast / lam_T - (15 / (math.pi ** 2 * f)) ** 0.25) < 1e-9
    # scaling with temperature: a_* ~ 1/T, i.e. ~ a(t)
    ratios = []
    for T in (T_CMB, 10 * T_CMB, 1000 * T_CMB):
        ratios.append(a_star(rho_gamma(T) * f) * T)
    c_scale = max(ratios) / min(ratios) - 1 < 1e-12
    print(f"      a_* T is constant across T = T_CMB, 10 T_CMB, 1000 T_CMB (a_* ~ 1/T ~ a(t))      : {c_scale}")
    ok &= c_scale
    # when does a_* pass through d_0 and l_*?
    kT_d0 = (15 / (math.pi ** 2 * f)) ** 0.25 * HBARC / D0          # J, photon+neutrino era
    print(f"      a_* = d_0 = 2.3 nm at k T = {kT_d0 / EV:.1f} eV  (T = {kT_d0 / KB:.2e} K, z ~ {kT_d0 / KB / T_CMB:.1e}, radiation era)")
    g_star = 106.75
    l_star = HBARC / (MSTAR_GEV * 1e9 * EV)
    ast_Mstar = (30 / (math.pi ** 2 * g_star)) ** 0.25 * l_star     # a_* at k T = M_* with the SM's g_*
    print(f"      at k T = M_* = {MSTAR_GEV:.0f} GeV (g_* = 106.75): a_* = {ast_Mstar:.2e} m = {ast_Mstar / l_star:.2f} l_*, l_* = {l_star:.2e} m")
    ok &= 50 < kT_d0 / EV < 120 and 0.3 < ast_Mstar / l_star < 0.6
    print("      => a_* is neither d_0 nor l_* at any fixed value; it is a function of the state, as comp 167 found and this file explains.")
    return ok


# ---------------------------------------------------------------------------
# Part 3: tamper tests
# ---------------------------------------------------------------------------
def part3():
    print("  PART 3.  Tamper tests")
    ok = True
    # a fixed a_* in an expanding geometry: residual 4 H rho; quantified for H_0
    H0 = 67.4e3 / 3.0857e22      # 1/s
    rr = rho_gamma(T_CMB) * neutrino_multiplier()
    residual = 4 * H0 * rr       # J/m^3/s
    print(f"      fixed a_*: residual 4 H_0 rho_r = {residual:.2e} J m^-3 s^-1, i.e. a relative non-conservation of 4 H_0 = {4 * H0:.2e} per second ({4 * H0 * 3.156e16:.2f} per Gyr)")
    ok &= residual > 0
    # a non-traceless Pi would be conserved with a different exponent, contradicting comp 158's zero trace
    ok_dust = True
    try:
        div, (a, rho, t) = frw_divergence(sp.Rational(1, 2))   # an arbitrary w
        H = sp.diff(a, t) / a
        ok_dust = sp.simplify(div[0] - (sp.diff(rho, t) + sp.Rational(9, 2) * H * rho)) == 0
    except Exception:
        ok_dust = False
    print(f"      general w: rho-dot + 3(1 + w) H rho (checked at w = 1/2), so only the traceless case gives the fourth power : {ok_dust}")
    ok &= ok_dust
    return ok


def main():
    print("=" * 78)
    print("  COMPUTATION 196.  a_* is comoving; under the identification it is the radiation's thermal wavelength")
    print("=" * 78)
    results = [
        ("conservation on FRW forces rho ~ a^-4, so a fixed a_* is excluded and the spacing is comoving", part1()),
        ("under Pi^0_0 = rho_r, a_* = 0.82 mm today, the thermal wavelength of the radiation, scaling as a(t)", part2()),
        ("tamper tests: a fixed a_* fails conservation by 4 H rho; only the traceless case gives the fourth power", part3()),
    ]
    print("-" * 78)
    for i, (name, r) in enumerate(results, 1):
        print(f"  {i}. {name}: {'PASS' if r else 'FAIL'}")
    allok = all(r for _, r in results)
    print("-" * 78)
    print("  VERDICT:", "PASS" if allok else "FAIL")
    print("""
  The first structural item carried a_* as one number the theory owes.  It owes no
  such number: a_* is comoving, by the conservation law alone, and under the conformal
  identification it is the thermal wavelength of the state's radiation, 0.82 mm today.
  The magnitude of Pi_mn is supplied by the state, as T_mn's is in general relativity.
  What remains of the item is its first gap, the identification on the semidefinite
  metric-traceless cone (comp 164), and nothing else.""")
    return allok


if __name__ == "__main__":
    main()