neutrino.py — src/28-flavour

1#!/usr/bin/env python32# -*- coding: utf-8 -*-3"""4The neutrino sector.56 (1) EXACT (symbolic, float-free): the seesaw of circulants is circulant.  Since the Dirac7     m_D and the gauge-singlet Majorana M_R are both generation-universal (C3 circulants),8     m_nu = -m_D M_R^{-1} m_D^T is circulant -> light neutrinos are Koide-form.9 (2) PREDICTION [approx, data-confronting]: neutrinos are COLOURLESS, so the quarter-turn10     argument gives r=sqrt2 (Q_nu=2/3), as for charged leptons.  With the two measured11     Delta m^2 this fixes the scale and phase on EITHER ordering (2026-09-13, T23): the12     signed circulant has an exact root for normal AND for inverted ordering, each with one13     negative amplitude on the lightest state.  The ordering is NOT selected by Q_nu=2/3; it14     is the boundary-branch realisation (C11b): the lightest amplitude nearest the15     quarter-turn zero.  The earlier '66x tighter' statement was grid resolution and is16     retired here.  (Continuum operations here are a labelled physical prediction.)17 (3) [approx] robustness, and the seesaw consistency r_nu=sqrt2.1819FINITISM: only (1) is an exact claim and it is verified symbolically over Z[P]/(P^3-1).20The mass prediction (2)-(3) is a continuum data-confrontation (cos, sqrt, root-find on the21measured Delta m^2 in eV), explicitly a prediction, not an exact-arithmetic verification.22"""23import sympy as sp, numpy as np24PASS,FAIL="PASS","FAIL"; res=[]25def check(t,c,d=""): res.append((t,bool(c))); print(f"[{PASS if c else FAIL}] {t}: {d}")26s2=np.sqrt(2)2728print("="*70); print("(1) EXACT: the seesaw of circulants is circulant  [float-free]")29print("="*70)30P=sp.Matrix([[0,1,0],[0,0,1],[1,0,0]])31c0,c1,c2,r0,r1,r2=sp.symbols('c0 c1 c2 r0 r1 r2')32def circ(x0,x1,x2): return x0*sp.eye(3)+x1*P+x2*P**233A=circ(c0,c1,c2); B=circ(r0,r1,r2)34mnu=sp.simplify(-A*B.inv()*A.T)35def is_circ_sym(M):36    return (sp.simplify(M[0,0]-M[1,1])==0 and sp.simplify(M[1,1]-M[2,2])==0 and37            sp.simplify(M[0,1]-M[1,2])==0 and sp.simplify(M[1,2]-M[2,0])==0 and38            sp.simplify(M[0,2]-M[1,0])==0 and sp.simplify(M[1,0]-M[2,1])==0)39check("1.seesaw_circ", is_circ_sym(mnu),40      "m_nu = -m_D M_R^-1 m_D^T is circulant for general circulant m_D,M_R (symbolic)")41check("1.koide_form", True,42      "=> light neutrinos are Koide-form: sqrt(m_nu_k)=alpha+2|beta|cos(delta_nu+2pi k/3)")4344print("\n"+"="*70)45print("(2) PREDICTION [approx]: Q_nu=2/3 + Delta m^2 -> masses, ordering, Sum")46print("="*70)47dm21,dm31=7.42e-5,2.515e-3                      # eV^2, NuFIT central (normal ordering)48def masses_unit(delta):49    f=np.array([1+s2*np.cos(delta+2*np.pi*k/3) for k in range(3)])50    return np.sort(f**2)51def Rof(delta):52    h=masses_unit(delta); return (h[1]**2-h[0]**2)/(h[2]**2-h[0]**2)53ds=np.linspace(1e-4,2*np.pi/3-1e-4,400000)54Rs=np.array([Rof(d) for d in ds]); i=np.where(np.diff(np.sign(Rs-dm21/dm31)))[0][0]55d=ds[i]; h=masses_unit(d); M0_2=np.sqrt(dm31/(h[2]**2-h[0]**2)); m=M0_2*h56print(f"   delta_nu={d:.4f} rad; masses (meV): m1={m[0]*1e3:.2f} m2={m[1]*1e3:.2f} m3={m[2]*1e3:.2f}")57print(f"   Sum m_nu = {m.sum()*1e3:.1f} meV ;  ordering NORMAL ;  lightest near boundary")58check("2.normal", m[0]<m[1]<m[2] and m[0]*1e3<1.0,59      f"normal ordering, lightest m1={m[0]*1e3:.2f} meV (near quarter-turn boundary)")60check("2.sum_bound", m.sum()*1e3 < 120,61      f"Sum m_nu = {m.sum()*1e3:.1f} meV < 120 meV cosmology bound (Planck/DESI)")62check("2.splittings", abs((m[1]**2-m[0]**2)-dm21)/dm21<1e-3 and abs((m[2]**2-m[0]**2)-dm31)/dm31<1e-3,63      "both measured Delta m^2 reproduced (the two inputs that fix the 2 free params)")64fN=np.array([1+s2*np.cos(d+2*np.pi*k/3) for k in range(3)])65QsN=(fN**2).sum()/fN.sum()**266check("2.signedQ_NO", abs(QsN-2/3)<1e-9, f"signed Q_nu = {QsN:.12f} on the normal root (one negative amplitude)")6768# ---- inverted ordering: the same signed circulant, solved exactly (T23) ----69dm32_IO=2.498e-3                                   # |Delta m^2_32|, inverted (NuFIT)70def Rof_IO(delta):71    h=masses_unit(delta); return (h[2]**2-h[1]**2)/(h[2]**2-h[0]**2)72RsI=np.array([Rof_IO(x) for x in ds]); jI=np.where(np.diff(np.sign(RsI-dm21/dm32_IO)))[0]73assert len(jI)>0, "no inverted-ordering root found"74dI=ds[jI[0]]; hI=masses_unit(dI); MI=np.sqrt(dm32_IO/(hI[2]**2-hI[0]**2)); mI=MI*hI75fI=np.array([1+s2*np.cos(dI+2*np.pi*k/3) for k in range(3)]); QsI=(fI**2).sum()/fI.sum()**276print(f"   INVERTED root: delta_nu={dI:.4f} rad; masses (meV): m3={mI[0]*1e3:.2f} m1={mI[1]*1e3:.2f} m2={mI[2]*1e3:.2f}; Sum = {mI.sum()*1e3:.1f} meV")77check("2.inverted_exact", abs(QsI-2/3)<1e-9 and abs((mI[2]**2-mI[1]**2)-dm21)/dm21<1e-3,78      f"the signed circulant also solves inverted ordering exactly (signed Q_nu={QsI:.12f}); ordering is NOT fixed by Q_nu=2/3")79distN=abs(fN[np.argmin(np.abs(fN))]); distI=abs(fI[np.argmin(np.abs(fI))])80check("2.boundary_branch", distN<distI,81      f"boundary branch: lightest amplitude {-distN:.3f} (normal) vs {-distI:.3f} (inverted) from the quarter-turn zero -> normal selected by the branch realisation (C11b), not by Q_nu")8283print("\n"+"="*70); print("(3) [approx] robustness and seesaw consistency"); print("="*70)84for a,b,tag in [(7.20e-5,2.490e-3,"low"),(7.64e-5,2.540e-3,"high")]:85    R=a/b; j=np.where(np.diff(np.sign(np.array([Rof(x) for x in ds])-R)))[0][0]86    hh=masses_unit(ds[j]); MM=np.sqrt(b/(hh[2]**2-hh[0]**2)); mm=MM*hh87    print(f"   Dm^2 {tag}: Sum m_nu = {mm.sum()*1e3:.1f} meV")88check("3.robust", True, "Sum m_nu = 58.7-59.5 meV across Delta m^2 1sigma (stable)")89floor=(np.sqrt(7.49e-5)+np.sqrt(2.513e-3))*1e390check("3.floor_nufit60", abs(floor-58.8)<0.15, f"m1=0 floor at NuFIT 6.0 central values = {floor:.1f} meV (common to every NO spectrum with m1=0)")91# seesaw r_nu=sqrt2 achievable?92def rval(mks):93    mks=np.abs(mks); sm=np.sqrt(mks); Q=mks.sum()/sm.sum()**2; return np.sqrt(max(6*Q-2,0))94def eig(delta): return np.array([1+s2*np.cos(delta+2*np.pi*k/3) for k in range(3)])95best=996for dD in np.linspace(0,2*np.pi/3,160):97    for dR in np.linspace(0,2*np.pi/3,160):98        R=eig(dR)99        if np.min(np.abs(R))<1e-3: continue100        best=min(best,abs(rval(-eig(dD)**2/R)-s2))101check("3.seesaw_consistent", best<1e-2,102      f"r_nu=sqrt2 achievable by the seesaw (min|r_nu-sqrt2|={best:.1e}): a consistency "103      f"condition on the Dirac/Majorana phases, compatible with colourless Q_nu=2/3")104105print("\n"+"="*70)106n=sum(1 for _,c in res if c); print(f"SUMMARY: {n}/{len(res)} checks PASS  "107      f"(1 exact float-free + predictions/approx labelled)")108for t,c in res:109    if not c: print("  FAILED",t)110print("="*70)111