#!/usr/bin/env python3
"""
SCOPE NOTE, 2026-10-11 (comp 202).  Result 1 below attributes the 0.0004 between its
0.09228 and the anchor 0.0927 to three-loop running and threshold matching.  comp 202 codes
both: the third loop is worth +0.00007 at M_* and the anchor itself is not reproducible.
The stability results (2 to 4) stand; the loop-order sequence the book quotes is now
comp 202's.

PROVENANCE: NUMERICAL

Computation 198 -- stability of the regulated vacuum at two loops: the running
quartic is positive on the whole range the theory contains
==============================================================================
WHAT THIS DOES.  The second structural item (sec:open-structural) is a rigorous
treatment of the interacting substrate vacuum beyond the free-field and one-loop
results of sec:renorm, and the book states the open part as "stability is settled at
leading loop order only".  Below M_* the theory is the Standard Model with no tower
and no threshold (sec:renorm, "Renormalisation group"), and above M_* it has no modes
at all: the cutoff is structural, the grain l_* = hbar c / M_* being a property of
the projection and not a regulator to be removed.  So the vacuum's stability beyond
one loop is a question about the Standard Model's effective potential on the field
range |phi| <= M_*, and that question has an answer at two loops in the
renormalisation-group-improved sense: the potential is V(phi) = lambda(mu = phi)
phi^4 / 4 up to the resummed leading logarithms, and it is stable on the range iff
lambda(mu) > 0 for mu <= M_*.

WHAT IS COMPUTED.  comp 123's two-loop running of {g1, g2, g3, y_t, lambda} is
reproduced here (copied, not imported; the book's scripts do not import siblings
across a cut) and completed with the two-loop beta function of lambda itself, which
comp 123 anchored to the literature value rather than coding.  The two-loop
beta_lambda is the standard one (Ford, Jack and Jones 1992, as reproduced in
Degrassi et al. 2012 and Buttazzo et al. 2013), with the GUT-normalised g1.

RESULTS.
  1. VALIDATION.  lambda(M_*) at one loop is comp 123's 0.09145 and at two loops
     without the lambda piece its 0.09234, both to five decimals; the two-loop
     lambda terms then move it by -0.00006 to 0.09228, which is within the 0.0007
     error of the literature value 0.0927 comp 123 anchored to (the remaining 0.0004
     is the three-loop running and the threshold matching Buttazzo et al. include and
     this file does not).  And the scale at which lambda crosses zero, the Standard
     Model's instability scale, lands at 10^9.7 GeV, inside the known window 10^9.5
     to 10^11 GeV for the MS-bar quartic.  So the added terms are the right ones.
  2. STABILITY ON THE RANGE.  lambda(mu) decreases monotonically from m_t to M_* and
     its minimum on the range is lambda(M_*) = 0.0927 > 0, so the RG-improved
     potential is positive and convex at every field value the regulated theory
     contains.  Margin: the instability scale is six and a half decades above M_*,
     and the top mass would have to exceed 220 GeV for lambda to cross zero below
     M_* (the file solves for the threshold).  Varying m_t, m_h and alpha_s across
     their 3 sigma boxes moves lambda(M_*) by less than 0.01, and its lower edge
     stays above 0.088.
  3. THE ONE-LOOP SIGN SURVIVES.  The book's one-loop result that the quantum
     correction deepens the sombrero (delta mu^2 > 0) is a statement about
     mu^2(mu), whose running is governed by gamma_m = 6 lambda + 3 y_t^2 - ... at one
     loop; the two-loop correction from the self-interaction is the 0.08% shift the
     book already quotes (sec:renorm).  This file does not redo that; it checks that
     the running mass parameter keeps its sign over the range, which it does since
     its anomalous dimension is a multiplicative factor.

WHAT THIS DOES TO THE ITEM.  "Stability beyond leading loop order" is established
at two loops in the RG-improved sense on the whole field range the theory contains,
with a margin of six and a half decades in scale.  What the item still carries is the other
clause, the rigorous CONSTRUCTION of the interacting vacuum at all orders (a Fock
space and a Hamiltonian bounded below, nonperturbatively), which no perturbative
statement supplies.  The item stays open, narrowed to that construction.

CHECKS.  (1) lambda(M_*) at one loop, at two loops without the lambda piece, and at
full two loops, against comp 123's three anchors 0.09145, 0.09234, 0.0927.  (2) The
instability scale at full two loops in [10^9, 10^12] GeV.  (3) lambda > 0 on
[m_t, M_*], monotone, minimum at M_*.  (4) The m_t threshold at which lambda(M_*) = 0,
above 200 GeV.  (5) The 3 sigma input box.  (6) Tamper: dropping the y_t^4 term of
the one-loop beta removes the instability altogether (lambda never crosses zero below
10^19 GeV), which locates the Standard Model's instability in the top Yukawa, as is
known, and shows the check is sensitive to the terms it carries.
"""

import math

PI2 = 16 * math.pi ** 2

# --- physical inputs (comp 123) --------------------------------------------
MH = 125.25            # GeV
V = 246.22             # GeV
MT_REF = 173.34        # GeV
MSTAR = 4 * math.pi * MH * math.sqrt(2.0 / 3.0)   # 1285 GeV
PST = 0.25 * math.exp(-1.0)                        # e^-1 / 4
BC0 = dict(g1=0.46270, g2=0.64779, g3=1.16655, yt=0.93690, lam=0.12604)
DROP_YT4 = False       # Part 4's tamper switch: remove the -6 y_t^4 term of the one-loop beta


def betas(y, order="full"):
    """SM beta functions; y = [g1(GUT), g2, g3, yt, lam].
    order: "1" one loop; "gy" two-loop gauge and top Yukawa only (comp 123); "full"
    adds the two-loop beta of lambda."""
    g1, g2, g3, yt, lam = y
    b1 = (41.0 / 10) * g1 ** 3
    b2 = (-19.0 / 6) * g2 ** 3
    b3 = (-7.0) * g3 ** 3
    byt = yt * (4.5 * yt ** 2 - (17.0 / 20) * g1 ** 2 - (9.0 / 4) * g2 ** 2 - 8 * g3 ** 2)
    blam = (24 * lam ** 2 + 12 * lam * yt ** 2 - (0 if DROP_YT4 else 6) * yt ** 4
            - 3 * lam * (3 * g2 ** 2 + (3.0 / 5) * g1 ** 2)
            + (3.0 / 8) * (2 * g2 ** 4 + ((3.0 / 5) * g1 ** 2 + g2 ** 2) ** 2))
    if order != "1":
        b1 += g1 ** 3 * ((199.0 / 50) * g1 ** 2 + (27.0 / 10) * g2 ** 2
                         + (44.0 / 5) * g3 ** 2 - (17.0 / 10) * yt ** 2) / PI2
        b2 += g2 ** 3 * ((9.0 / 10) * g1 ** 2 + (35.0 / 6) * g2 ** 2
                         + 12 * g3 ** 2 - (3.0 / 2) * yt ** 2) / PI2
        b3 += g3 ** 3 * ((11.0 / 10) * g1 ** 2 + (9.0 / 2) * g2 ** 2
                         - 26 * g3 ** 2 - 2 * yt ** 2) / PI2
        byt2 = (-12 * yt ** 4
                + yt ** 2 * ((393.0 / 80) * g1 ** 2 + (225.0 / 16) * g2 ** 2 + 36 * g3 ** 2)
                + (1187.0 / 600) * g1 ** 4 - (23.0 / 4) * g2 ** 4 - 108 * g3 ** 4
                - (9.0 / 20) * g1 ** 2 * g2 ** 2 + (19.0 / 15) * g1 ** 2 * g3 ** 2
                + 9 * g2 ** 2 * g3 ** 2)
        byt += yt * byt2 / PI2
    if order == "full":
        g, gp, gs = g2, math.sqrt(3.0 / 5) * g1, g3        # gp is the hypercharge coupling
        blam2 = (-312 * lam ** 3 - 144 * lam ** 2 * yt ** 2 + 36 * lam ** 2 * (3 * g ** 2 + gp ** 2)
                 - 3 * lam * yt ** 4 + lam * yt ** 2 * (80 * gs ** 2 + 22.5 * g ** 2 + (85.0 / 6) * gp ** 2)
                 - (73.0 / 8) * lam * g ** 4 + (39.0 / 4) * lam * g ** 2 * gp ** 2 + (629.0 / 24) * lam * gp ** 4
                 + 30 * yt ** 6 - 32 * yt ** 4 * gs ** 2 - (8.0 / 3) * yt ** 4 * gp ** 2
                 - (9.0 / 4) * yt ** 2 * g ** 4 + 10.5 * yt ** 2 * g ** 2 * gp ** 2 - (19.0 / 4) * yt ** 2 * gp ** 4
                 + (305.0 / 16) * g ** 6 - (289.0 / 48) * g ** 4 * gp ** 2 - (559.0 / 48) * g ** 2 * gp ** 4
                 - (379.0 / 48) * gp ** 6)
        blam += blam2 / PI2
    return [b1 / PI2, b2 / PI2, b3 / PI2, byt / PI2, blam / PI2]


def run(mu0, mu1, bc, order="full", n=3000, trace=False):
    t0, t1 = math.log(mu0), math.log(mu1)
    h = (t1 - t0) / n
    y = [bc["g1"], bc["g2"], bc["g3"], bc["yt"], bc["lam"]]
    path = []
    for i in range(n):
        k1 = betas(y, order)
        k2 = betas([y[j] + 0.5 * h * k1[j] for j in range(5)], order)
        k3 = betas([y[j] + 0.5 * h * k2[j] for j in range(5)], order)
        k4 = betas([y[j] + h * k3[j] for j in range(5)], order)
        y = [y[j] + (h / 6) * (k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) for j in range(5)]
        if trace:
            path.append((math.exp(t0 + (i + 1) * h), y[4]))
    return (y, path) if trace else y


def instability_scale(bc=BC0, order="full"):
    lo, hi = 1e6, 1e19
    if run(MT_REF, hi, bc, order, n=2000)[4] > 0:
        return float('inf')
    for _ in range(28):
        mid = math.sqrt(lo * hi)
        if run(MT_REF, mid, bc, order, n=1500)[4] > 0:
            lo = mid
        else:
            hi = mid
    return lo


def bc_shift(dmt=0.0, das=0.0, dmh=0.0, bc=BC0):
    b = dict(bc)
    b["yt"] = bc["yt"] * (1 + dmt / MT_REF)
    as0 = bc["g3"] ** 2 / (4 * math.pi)
    b["g3"] = math.sqrt(4 * math.pi * (as0 + das))
    b["lam"] = bc["lam"] + 0.00206 * dmh
    return b


def part1():
    print("  PART 1.  Validation against comp 123's anchors and the known instability scale")
    ok = True
    l1 = run(MT_REF, MSTAR, BC0, "1")[4]
    lgy = run(MT_REF, MSTAR, BC0, "gy")[4]
    lf = run(MT_REF, MSTAR, BC0, "full")[4]
    print(f"      lambda(M_*): one loop {l1:.5f} (comp 123: 0.09145);  two-loop g + y_t {lgy:.5f} (comp 123: 0.09234);  full two loops {lf:.5f} (literature anchor 0.0927 +- 0.0007)")
    ok &= abs(l1 - 0.09145) < 2e-5 and abs(lgy - 0.09234) < 2e-5 and abs(lf - 0.0927) < 0.0007
    print(f"      the two-loop lambda terms move lambda(M_*) by {lf - lgy:+.5f}; the rest of the distance to 0.0927 is three-loop running and threshold matching, outside this file")
    ok &= abs(lf - lgy) < 0.0005
    mu_I = instability_scale()
    print(f"      instability scale (lambda -> 0) at full two loops: {mu_I:.2e} GeV = 10^{math.log10(mu_I):.2f}; known window 10^9.5 to 10^11")
    ok &= 1e9 < mu_I < 1e12
    print(f"      against the substrate value e^-1/4 = {PST:.5f}: {(lf - PST) / 0.0007:.1f} sigma (comp 123's test, unchanged in kind)")
    return ok


def part2():
    print("  PART 2.  lambda on the range the regulated theory contains, [m_t, M_*]")
    ok = True
    yM, path = run(MT_REF, MSTAR, BC0, "full", trace=True)
    lams = [l for _, l in path]
    mono = all(lams[i + 1] < lams[i] for i in range(len(lams) - 1))
    lmin = min(lams)
    print(f"      lambda decreases monotonically from {BC0['lam']:.5f} at m_t to {lams[-1]:.5f} at M_*            : {mono}")
    print(f"      minimum on the range {lmin:.5f} > 0, at M_* itself; V(phi) = lambda(phi) phi^4/4 is positive on |phi| <= M_*  : {lmin > 0}")
    ok &= mono and lmin > 0
    mu_I = instability_scale()
    print(f"      decades between M_* and the instability scale: {math.log10(mu_I / MSTAR):.1f}")
    ok &= math.log10(mu_I / MSTAR) > 5
    # the top-mass threshold at which lambda(M_*) = 0
    lo, hi = 0.0, 400.0
    for _ in range(40):
        mid = 0.5 * (lo + hi)
        if run(MT_REF, MSTAR, bc_shift(dmt=mid), "full")[4] > 0:
            lo = mid
        else:
            hi = mid
    mt_thr = MT_REF + lo
    print(f"      m_t at which lambda(M_*) would reach zero (y_t scaled with m_t, other inputs fixed): {mt_thr:.0f} GeV")
    ok &= mt_thr > 200
    # 3 sigma box
    vals = []
    for smt in (-1, 0, 1):
        for sas in (-1, 0, 1):
            for smh in (-1, 0, 1):
                vals.append(run(MT_REF, MSTAR, bc_shift(dmt=3 * 0.76 * smt, das=3 * 0.0009 * sas, dmh=3 * 0.17 * smh), "full")[4])
    print(f"      lambda(M_*) over the 3 sigma box of (m_t, alpha_s, m_h): [{min(vals):.5f}, {max(vals):.5f}], spread {max(vals) - min(vals):.5f}")
    ok &= min(vals) > 0.085 and max(vals) - min(vals) < 0.01
    return ok


def part3():
    print("  PART 3.  The mass parameter keeps its sign: its running is multiplicative")
    # d ln m^2 / d ln mu = gamma_m2 = (12 lam + 6 yt^2 - 9/2 g2^2 - 9/10 g1^2) / (16 pi^2) at one loop
    y = [BC0["g1"], BC0["g2"], BC0["g3"], BC0["yt"], BC0["lam"]]
    g1, g2, g3, yt, lam = y
    gam = (12 * lam + 6 * yt ** 2 - 4.5 * g2 ** 2 - 0.9 * g1 ** 2) / PI2
    factor = math.exp(gam * math.log(MSTAR / MT_REF))
    print(f"      one-loop anomalous dimension of m^2 at m_t: {gam:.4f};  m^2(M_*)/m^2(m_t) = {factor:.3f} > 0, so the sign of mu^2 that drives SSB is the same at every scale on the range")
    print("      the two-loop self-interaction shift of m_h is the 0.08% the book quotes (sec:renorm); it does not flip a sign")
    return factor > 0


def part4():
    print("  PART 4.  Tamper: drop the y_t^4 term of the one-loop beta")
    global DROP_YT4
    DROP_YT4 = True
    try:
        mu_I = instability_scale()
    finally:
        DROP_YT4 = False
    gone = mu_I == float('inf')
    print(f"      without -6 y_t^4 the quartic never crosses zero below 10^19 GeV                  : {gone}")
    print("      so the instability is the top Yukawa's, as is known, and the check is sensitive to the terms it carries")
    return gone


def main():
    print("=" * 78)
    print("  COMPUTATION 198.  Stability of the regulated vacuum at two loops")
    print("=" * 78)
    results = [
        ("full two-loop running reproduces comp 123's anchors and the known instability scale", part1()),
        ("lambda > 0 on [m_t, M_*] at two loops, six and a half decades of margin, m_t threshold far above its value", part2()),
        ("the mass parameter's running is multiplicative, the SSB sign survives", part3()),
        ("tamper: removing the top Yukawa's quartic term removes the instability", part4()),
    ]
    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("""
  Stability beyond leading loop order is established at two loops in the RG-improved
  sense on the whole field range the regulated theory contains, lambda(M_*) = 0.0923
  being the minimum of a positive monotone function and the instability scale six and
  a half decades above the cutoff.  What the second structural item still carries is the
  construction of the interacting vacuum at all orders; the item is narrowed to it.""")
    return allok


if __name__ == "__main__":
    main()
