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

Computation 109 -- Higgs-quartic normalisation + O(4) one-loop self-energy:
                   explicit verification of two referee findings
=========================================================================
A referee note on v26.03 flagged two compounding issues in the Higgs-sector
derivation that together drive the M_* = 4 pi m_h claim:

  Finding A: postulate (1) has quartic (1/2) <psi,psi>^2; working potential
             (23) has (1/4) phi^4.  Are they reconcilable?

  Finding B: paper's one-loop self-energy uses "4 lambda" (N lambda for
             N=4 real components).  Standard O(N) calculation gives
             (N+2) lambda = 6 lambda.  Which is right?

This computation performs both verifications symbolically.  No code paths
are skipped; every contraction is enumerated explicitly.

FINDING A: V_7 reduction of postulate (1)
==========================================
Postulate (1) in the paper:

    F[psi, eps] = integral [ -eps <psi,psi> + (1/4) <psi,psi>^2
                              + c |grad psi|^2 ] dmu                  (1)

with psi V_7-valued.  V_7 = Im(O) carries the inner product
<v,w> = sum_a v_a w_a (orthonormal basis from the octonion norm; equivalent
to the Euclidean inner product on R^7 under the standard identification).

Radial reduction:
    Write psi = phi * n-hat with n-hat in S^6 (unit V_7 vector) and
    phi a real magnitude.  Then:
        <psi, psi>   = phi^2 * <n-hat, n-hat> = phi^2 * 1 = phi^2
        <psi, psi>^2 = phi^4

    Substituting into (1):
        F[phi, eps] = int[ -eps phi^2 + (1/4) phi^4 + c|grad phi|^2 + ... ]

In the committed doublet convention the postulate quartic is (1/4), so the
radial reduction gives (1/4) phi^4, matching the working potential (23).

The paper's working potential (Eq 23) reads:
    F(phi) = -eps phi^2 + (1/4) phi^4

with quadratic coefficient -eps (same as postulate) and quartic (1/4).

Reconciliation attempt: a field rescaling phi -> alpha phi gives
    -eps phi^2 -> -alpha^2 eps phi^2
    (1/2) phi^4 -> (alpha^4/2) phi^4

Going from (1/2) -> (1/4) requires alpha^4 = 1/2 i.e. alpha = 2^(-1/4).
Then the quadratic becomes -2^(-1/2) eps phi^2 = -eps/sqrt(2) * phi^2.
The paper's (23) has quadratic -eps, NOT -eps/sqrt(2).  So the rescaling
that fixes the quartic spoils the quadratic.  There is no field rescaling
that takes (1) to (23) simultaneously in both coefficients.

VERDICT (A): the referee is correct that the paper, as originally
written, stated the quartic in inconsistent normalisations.

RESOLUTION (adopted in v26.09): the paper commits to the DOUBLET
CONVENTION throughout.  In this convention:
  - <psi,psi> is identified with the doublet bilinear Phi^dag Phi
    (NOT the real-component norm; the two differ by a factor of 2).
  - lambda_PST = 1/4 is the coefficient of (Phi^dag Phi)^2, the SAME
    convention as the SM lambda_SM in V_SM = lambda_SM (Phi^dag Phi)^2.
  - The postulate (1) is corrected to (1/4)<psi,psi>^2, so that the
    radial reduction <psi,psi> -> psi^2 gives (1/4)psi^4, matching the
    working potential (23) with no factor mismatch.
  - The complex-component sombrero of sec:mstar (F = -eps|psi|^2 +
    (1/2)|psi|^4) is a different field normalisation used only for the
    m_h/v/r_0 relations; it is bridged explicitly to the doublet
    convention in that section.
  - Z^2 = lambda_SM(M_*)/lambda_PST is then a ratio of two same-
    convention couplings, = 4 lambda_SM(M_*) ~ 0.37 ~ e^(-1).

The choice matters: b enters the bridge premise (B) prediction
lambda_SM(M_*) = b * e^(-1).  b = 1/4 gives 0.092; b = 1/2 gives 0.184.


FINDING B: O(N) one-loop scalar self-energy
============================================
Lagrangian (symmetric phase, all four Higgs components live in a
common O(4) multiplet):
    V = lambda (Phi^dagger Phi)^2
      = (lambda / 4) (phi . phi)^2
      = (lambda / 4) sum_{a,b} phi_a^2 phi_b^2

4-point vertex tensor:
    V_{ijkl} = partial^4 V / partial phi_i partial phi_j ...   evaluated at phi = 0

Direct computation (see Part 2 below):
    V_{ijkl} = 2 lambda (delta_{ij} delta_{kl}
                       + delta_{ik} delta_{jl}
                       + delta_{il} delta_{jk})

Feynman rule = -i V_{ijkl}.

One-loop seagull self-energy: close legs k, l into the loop with propagator
i/p^2.  Loop integral = i * Lambda^2 / (16 pi^2) (Euclidean cutoff Lambda).
Symmetry factor 1/2 (two internal legs are interchangeable).

  Sigma_{ij} = (1/2) * (-i V_{ijkl}) * (i / loop) * delta_{kl}
             = (1/2) * V_{ijkl} * delta_{kl} * Lambda^2 / (16 pi^2)

Contracting V_{ijkl} with delta_{kl}:

  V_{ijkl} delta_{kl}
    = 2 lambda [ delta_{ij} delta_{kl} delta_{kl}
               + delta_{ik} delta_{jl} delta_{kl}
               + delta_{il} delta_{jk} delta_{kl} ]
    = 2 lambda [ delta_{ij} * N  +  delta_{ij}  +  delta_{ij} ]
    = 2 lambda (N + 2) delta_{ij}

  Sigma_{ij} = (1/2) * 2 lambda (N + 2) delta_{ij} * Lambda^2 / (16 pi^2)
             = (N + 2) lambda * delta_{ij} * Lambda^2 / (16 pi^2)

For N = 4:  Sigma = 6 lambda * Lambda^2 / (16 pi^2).

The paper writes "4 lambda * M_*^2 / (16 pi^2)" (Eq 107).  That is N lambda,
keeping only the delta_{ij} delta_{kl} contraction and dropping the +2 from
the symmetric tensor structure of V_{ijkl}.

Cross-check via broken-phase contributions (SM textbook):
  - physical Higgs (h) self-loop: vertex 3! lambda v^2 -> 3 lambda
  - 3 Goldstone modes, each contributing 1 lambda
  - total: 3 lambda + 3 * 1 lambda = 6 lambda                  (agrees)

VERDICT (B): the referee is correct.  The one-loop self-energy
coefficient is (N+2) lambda = 6 lambda for N=4, not 4 lambda.


COMPOUND IMPACT ON M_* PREDICTION
==================================
With C_H = (loop coefficient) * lambda_PST and m_h^2 = C_H * M_*^2 / (16 pi^2):

    M_* = 4 pi m_h / sqrt(C_H)

For the four combinations of (count, lambda_PST):

  count=4, lambda=1/4 (original reading, superseded)    -> C_H = 1     -> M_* = 4 pi m_h ~= 1573 GeV
  count=6, lambda=1/4 (paper b, fix B)  -> C_H = 3/2   -> M_* ~= 1285 GeV
  count=4, lambda=1/2 (fix A, keep B)   -> C_H = 2     -> M_* ~= 1113 GeV
  count=6, lambda=1/2 (both fixes)      -> C_H = 3     -> M_*  ~= 909 GeV

The headline 1573 GeV is the unique corner where both errors hold.
Either correction alone shifts M_*.


PYTHON VERIFICATION
====================
Symbolic verification of the V_{ijkl} tensor and the loop contraction
for arbitrary N, then evaluation at N = 4.
"""

import math
from itertools import product


def vijkl(i, j, k, l, lam):
    """
    Vertex tensor V_{ijkl} = 2*lambda * (delta_{ij} delta_{kl}
                                        + delta_{ik} delta_{jl}
                                        + delta_{il} delta_{jk}).

    Derived from V = (lambda/4) (sum_a phi_a^2)^2 via four derivatives.
    """
    def kron(a, b): return 1 if a == b else 0
    return 2 * lam * (kron(i, j) * kron(k, l)
                      + kron(i, k) * kron(j, l)
                      + kron(i, l) * kron(j, k))


def loop_contraction_coefficient(N, lam):
    """
    Compute Sigma_{ii} after closing the seagull loop:

        Sigma_{ij} = (1/2) * V_{ijkl} * delta_{kl} * (loop integral)

    Return the coefficient in front of Lambda^2 / (16 pi^2).
    Verify that Sigma_{ij} is isotropic (proportional to delta_{ij}).
    """
    sigma_diag = 0
    sigma_off = 0
    for k in range(N):
        sigma_diag += (1 / 2) * vijkl(0, 0, k, k, lam)
        sigma_off += (1 / 2) * vijkl(0, 1, k, k, lam)
    return sigma_diag, sigma_off


def main():
    print("=" * 72)
    print("Computation 109: Higgs-quartic normalisation + O(4) self-energy")
    print("=" * 72)
    print()
    print("FINDING A: V_7 radial reduction of postulate (1)")
    print("-" * 72)
    print("  psi V_7-valued, <psi,psi> = sum_a psi_a^2 = phi^2 (radial)")
    print("  <psi,psi>^2 = phi^4")
    print()
    print("  Postulate (1) coefficient:  (1/2) <psi,psi>^2  ->  (1/2) phi^4")
    print("  Working potential (23):     (1/4) phi^4")
    print()
    print("  Rescaling test: phi -> alpha*phi to send (1/2) -> (1/4)")
    print("    requires alpha^4 = 1/2, so alpha = 2^(-1/4) ~= 0.8409")
    print(f"    alpha = {2 ** (-1/4):.6f}")
    print(f"    new quadratic coeff = -alpha^2 * eps = -{2 ** (-1/2):.6f} eps")
    print("    paper (23) quadratic = -eps  (no rescaling factor)")
    print("    => no single rescaling reconciles BOTH coefficients.")
    print()
    print("  VERDICT: postulate (1) and working (23) are in different")
    print("           normalisations.  Cannot be unified by field rescaling.")
    print()

    print("FINDING B: O(N) one-loop self-energy")
    print("-" * 72)
    print("  Vertex V_{ijkl} = 2*lambda * (d_ij d_kl + d_ik d_jl + d_il d_jk)")
    print()
    print("  Sample V_{ijkl} values (lambda = 1):")
    for (i, j, k, l) in [(0, 0, 0, 0), (0, 0, 1, 1), (0, 1, 0, 1), (0, 1, 1, 0)]:
        v = vijkl(i, j, k, l, 1.0)
        print(f"    V_{{{i}{j}{k}{l}}} = {v}")
    print()
    print("  One-loop contraction (1/2) * V_{ijkl} * delta_{kl} for")
    print("  N in {1, 2, 3, 4, 5} with lambda = 1:")
    print()
    print(f"    {'N':>3} | {'diag (Sigma_00)':>15} | {'off (Sigma_01)':>15} | expected (N+2)*lambda")
    print(f"    {'-'*3}-+-{'-'*15}-+-{'-'*15}-+-{'-'*21}")
    for N in [1, 2, 3, 4, 5]:
        diag, off = loop_contraction_coefficient(N, 1.0)
        expected = N + 2
        match = "OK" if abs(diag - expected) < 1e-12 else "MISMATCH"
        print(f"    {N:>3} | {diag:>15.6f} | {off:>15.6f} | "
              f"{expected:>21} ({match})")
    print()
    print("  Off-diagonal (Sigma_01) is zero -> Sigma is isotropic, as expected.")
    print("  Diagonal coefficient matches (N+2) * lambda exactly for all N.")
    print()
    print("  VERDICT: one-loop self-energy = (N+2) lambda * Lambda^2/(16 pi^2)")
    print("           For N = 4:  6 lambda, NOT 4 lambda.")
    print()
    print("  Cross-check: broken-phase total")
    print("    physical Higgs h-loop (vertex 6 lambda v^2): contributes 3 lambda")
    print("    3 Goldstones, each 1 lambda: contributes 3 lambda")
    print("    SUM: 3 + 3 = 6 lambda.  Agrees with O(4) symmetric-phase result.")
    print()

    print("COMPOUND IMPACT ON M_* PREDICTION")
    print("-" * 72)
    print("  M_* = 4 pi m_h / sqrt(C_H),  with C_H = (count) * lambda_PST")
    print(f"  m_h = 125.25 GeV  =>  4 pi m_h = {4 * math.pi * 125.25:.2f} GeV")
    print()
    print(f"    {'count':>5} | {'lambda_PST':>10} | {'C_H':>6} | "
          f"{'M_* (GeV)':>12} | note")
    print(f"    {'-'*5}-+-{'-'*10}-+-{'-'*6}-+-{'-'*12}-+-{'-'*30}")
    for count, lam, note in [(4, 0.25, "original reading, superseded"),
                              (6, 0.25, "fix B only (paper's b)"),
                              (4, 0.50, "fix A only (paper's count)"),
                              (6, 0.50, "both fixes")]:
        c_h = count * lam
        m_star = 4 * math.pi * 125.25 / math.sqrt(c_h)
        print(f"    {count:>5} | {lam:>10.3f} | {c_h:>6.2f} | "
              f"{m_star:>12.2f} | {note}")
    print()
    print("  The 1573 GeV headline is the unique corner where both undercounts")
    print("  simultaneously hold.")
    print()

    print("IMPACT ON BRIDGE PREMISE (B) NUMERICAL CLAIM")
    print("-" * 72)
    print("  lambda_SM(M_*) = b * Z_H(beta_KO) = b * e^(-1)")
    print(f"  e^(-1) = {math.exp(-1):.6f}")
    print()
    print(f"    {'b':>6} | {'b * e^(-1)':>11} | {'vs lambda_SM ~= 0.0927':>22} | match?")
    print(f"    {'-'*6}-+-{'-'*11}-+-{'-'*22}-+-{'-'*10}")
    for b, label in [(0.25, "paper's b"), (0.50, "postulate's b (Finding A)")]:
        pred = b * math.exp(-1)
        ratio_pct = 100 * abs(pred - 0.0927) / 0.0927
        print(f"    {b:>6.3f} | {pred:>11.6f} | "
              f"{ratio_pct:>21.1f}% | {label}")
    print()
    print("  If Finding A holds (b = 1/2), the 0.8% match becomes ~100% off.")
    print()

    print("=" * 72)
    print("Summary")
    print("=" * 72)
    print()
    print("  Finding A (1/2 vs 1/4 in postulate vs working potential):")
    print("    confirmed by direct V_7 radial reduction; no rescaling")
    print("    reconciles both coefficients.")
    print()
    print("  Finding B (4 lambda vs 6 lambda in one-loop self-energy):")
    print("    confirmed by explicit O(N) vertex contraction for all N tested;")
    print("    agrees with the broken-phase textbook calculation 3*lambda")
    print("    (physical Higgs) + 3*lambda (Goldstones) = 6*lambda for N=4.")
    print()
    print("  Combined: M_* = 4 pi m_h ~= 1573 GeV is the unique corner where")
    print("    both undercounts hold.  Any correction shifts M_* into")
    print("    the 900-1285 GeV range and breaks the clean 4 pi factor.")
    print()
    print("  This computation does not propose a fix.  It confirms the referee's")
    print("  technical claims at proof-detail level.")


if __name__ == "__main__":
    main()
