Integrate $i\,\partial_t|\psi\rangle=H_I|\psi\rangle$ from $t_0$ to $t$ and iterate:

In the $n$-th term the times are already ordered $t_1>t_2>\cdots$; symmetrizing the nested limits over all $n!$ orderings and dividing by $n!$ replaces the ordering constraint by the time-ordering symbol $T$, resumming to $U=T\exp(-i\int H_I\,dt')$. The $T$ is essential because $[H_I(t_1),H_I(t_2)]\ne0$, so the factors must be stacked later-to-the-left. Writing $H_I=-\int d^3x\,\mathcal{L}_{\mathrm{int}}$ makes it covariant, $S=U(\infty,-\infty)=T\exp(i\int d^4x\,\mathcal{L}_{\mathrm{int}})$, and each $\mathcal{L}_{\mathrm{int}}\propto\lambda$ is one vertex, so the $n$-th term is order $\lambda^n$. The lowest order scattering $\phi\phi\to\phi\phi$ is a single vertex ($\lambda^1$), hence $|\mathcal{M}|^2=\lambda^2$: phi4_amplitude_squared(0.3)$=0.09$.

Wick's theorem gives $T\,\phi(x)\phi(y)={:}\phi(x)\phi(y){:}+\langle0|T\,\phi(x)\phi(y)|0\rangle$. The normal-ordered piece has all annihilators on the right, so $\langle0|{:}\phi\phi{:}|0\rangle=0$ and only the single contraction survives — and that contraction is the Feynman propagator $D_F(x-y)$. Fourier-transforming,

On shell $p^2=m^2$ the denominator would vanish; the $+i\epsilon$ holds it at $i\epsilon$, so $|\widetilde D_F|\to1/\epsilon$ stays finite, and the infinitesimally displaced pole fixes the causal boundary condition (positive energy forward, antiparticles backward). The code shows both limits: on shell propagator(np.array([np.sqrt(5),2,0,0]), 1.0, 1e-6) has $|D_F|\approx10^6=1/\epsilon$ (the regulated pole), while off shell at $p^2=9$, propagator([3,0,0,0],1.0)$=i/(9-1)=0.125\,i$.

Expand each invariant, e.g. $s=p_1^2+2\,p_1\!\cdot\!p_2+p_2^2$, $t=p_1^2-2\,p_1\!\cdot\!p_3+p_3^2$, $u=p_1^2-2\,p_1\!\cdot\!p_4+p_4^2$, and add:

Momentum conservation $p_1+p_2=p_3+p_4$ gives $p_2-p_3-p_4=-p_1$, so the cross term is $-2p_1^2$ and the $3p_1^2$ collapses to $p_1^2$, leaving

where $p_i^2=m_i^2$ is the on-shell condition. Only two of $s,t,u$ are independent. For $m=1$, mandelstam(5.0, np.pi/3, 1.0)$=(25,-5.25,-15.75)$ sums to $4.0=$mandelstam_sum(1.0)$=4m^2$.

In the CM frame the 3-momenta cancel, $\mathbf p_1+\mathbf p_2=0$, so $s=(p_1+p_2)^2=(E_1+E_2)^2=E_{\mathrm{cm}}^2$ — $\sqrt s$ is the total energy. With both legs on-shell and equal mass, $E_1=E_2=\sqrt{|\mathbf p|^2+m^2}=\tfrac12\sqrt s$, which inverts to

This is real only for $s\ge4m^2$, i.e. $\sqrt s\ge2m$ — the threshold to make two real particles; below it cm_momentum(1.5,1.0)$=0$ and at it cm_momentum(2.0,1.0)$=0$ (no phase space), while above it cm_momentum(5.0,1.0)$=\tfrac12\sqrt{21}=2.2913$ and cm_energy of two back-to-back on-shell momenta returns $2E$. Then $t=-2|\mathbf p|^2(1-\cos\theta)$ and $u=-2|\mathbf p|^2(1+\cos\theta)$ are both $\le0$, sweeping from $t=0$ (forward, $\cos\theta=1$) to $u=0$ (backward, $\cos\theta=-1$).

The first-order S-matrix term is $i\int d^4x\,\mathcal{L}_{\mathrm{int}}=-\tfrac{i\lambda}{4!}\int d^4x\,\phi^4$. Attaching the four external legs $\phi\phi\to\phi\phi$ to the four fields at the vertex can be done in $4!$ ways — the four fields pair with the four distinct legs in any order — and these $4!$ identical contractions cancel the $1/4!$:

There is a single tree diagram — a point vertex, no internal propagator, no momentum transfer — so $|\mathcal{M}|^2=\lambda^2$ is isotropic, not forward-peaked: no $\theta$ dependence at this order. Hence phi4_amplitude_squared(0.3)$=0.09=\lambda^2$ for every scattering angle.

Elastic scattering keeps the same CM momentum in and out, $|\mathbf p_f|=|\mathbf p_i|$, so the kinematic ratio is $1$ and with $|\mathcal{M}|^2=\lambda^2$ the rate is angle-independent, $d\sigma/d\Omega=\lambda^2/(64\pi^2 s)$. Integrating an isotropic $d\sigma/d\Omega$ over the sphere gives $4\pi$, but the two final-state $\phi$'s are identical, so the full-solid-angle integral counts each distinct configuration twice; dividing by $2$,

At $\sqrt s=5,\ \lambda=0.3$: $\sigma=0.09/(32\pi\cdot25)=3.581\times10^{-5}$, exactly phi4_cross_section(5.0,0.3,1.0) and equal to $\tfrac12\cdot4\pi\cdot$phi4_differential_cross_section(5.0,0.3,1.0).

Read the closed form $\sigma=\lambda^2/(32\pi s)$ directly. The numerator $\lambda^2>0$ and denominator $s=E_{\mathrm{cm}}^2>0$ force $\sigma>0$; it is exactly $\propto\lambda^2$, so doubling the coupling quadruples it — phi4_cross_section(5.0,0.4,1.0)/phi4_cross_section(5.0,0.2,1.0)$=(0.4/0.2)^2=4.0$. The explicit $1/s$ makes $\sigma$ fall with energy,

so phi4_cross_section(3.0,0.3,1.0)$>$phi4_cross_section(6.0,0.3,1.0). Below threshold the incoming pair cannot be on-shell ($|\mathbf p_i|=\tfrac12\sqrt{s-4m^2}$ imaginary), there is no phase space, and $\sigma=0$ — phi4_cross_section(1.99,0.3,1.0)$=0.0$. The shape mirrors the non-relativistic ~QM-18: tree level $\leftrightarrow$ Born approximation, each first order in the coupling, each giving $d\sigma/d\Omega\propto|\text{amplitude}|^2$ — here $|\mathcal{M}|^2$, there $|f(\theta)|^2$.

For a global phase ($\alpha$ constant) $\partial_\mu(e^{i\alpha}\psi)=e^{i\alpha}\partial_\mu\psi$, so $\bar\psi\gamma^\mu\partial_\mu\psi$ is unchanged. Localizing $\alpha\to\alpha(x)$, the product rule produces a stray gradient,

and the extra $i(\partial_\mu\alpha)\psi$ spoils the invariance of $\bar\psi\gamma^\mu\partial_\mu\psi$. Introduce $D_\mu=\partial_\mu+ieA_\mu$ with the paired shift $A_\mu\to A_\mu-\tfrac1e\partial_\mu\alpha$; then under the combined transformation

the $i(\partial_\mu\alpha)\psi$ cancelling exactly against the $-\tfrac1e\partial_\mu\alpha$ shift. Hence $|D_\mu\psi|$ is gauge invariant while $|\partial_\mu\psi|$ is not — matching the check: $|D_\mu\psi|$ changes by only $\approx5\times10^{-3}$ (the $O(dx^2)$ finite-difference error) while the bare $|\partial_\mu\psi|$ jumps by $\approx1.0$.

Apply the shift $A_\mu\to A_\mu+\partial_\mu\lambda$ directly to the field strength:

Mixed partials commute, $\partial_\mu\partial_\nu\lambda=\partial_\nu\partial_\mu\lambda$, so the added term vanishes and $F_{\mu\nu}\to F_{\mu\nu}$ exactly. A pure-gauge potential is the special case $A_\mu=\partial_\mu\lambda$ (built from $A=0$), for which

identically. So $A_\mu$ carries gauge-dependent (unphysical) information — the piece $\partial_\mu\lambda$ that any $\lambda$ can add — while $F_{\mu\nu}$, being invariant, is what every observer agrees on and what is measured ($F_{0i}=E_i$, $F_{ij}=-\varepsilon_{ijk}B_k$). Numerically the gauge shift moves $F$ by only $\max|\Delta F|\approx2\times10^{-15}$ and the pure gauge gives $\max|F|\approx1.7\times10^{-15}$ — both zero to machine precision, since the discrete cross-derivatives along distinct axes commute exactly.

With $A_0=0$ and $A_1=Bx^0$,

a constant, and antisymmetry fixes $F_{10}=-F_{01}=-B$ with $F_{00}=F_{11}=0$. In the full tensor the time–space entries are the electric field, $F_{0i}=E_i$, and the space–space entries the magnetic field, $F_{ij}=-\varepsilon_{ijk}B_k$ — exactly ~EM-18's field tensor; here the single nonzero component is a uniform field of strength $B$ along the $x^1$ axis. Because np.gradient is exact on the linear field $Bx^0$, the code returns field_strength(np.stack([0X0, BX0]), dx)[0,1] $=B=0.75$ to $\approx10^{-10}$, with $F[1,0]=-F[0,1]=-0.75$ — confirming the value $B$ and the antisymmetry.

Insert $D_\mu=\partial_\mu+ieA_\mu$ into the Dirac term:

because $i\gamma^\mu\cdot ieA_\mu=-e\gamma^\mu A_\mu$. The cross term is exactly $-e\,j^\mu A_\mu$ with the Noether current $j^\mu=\bar\psi\gamma^\mu\psi$; adding the photon kinetic term gives

The interaction was never put in by hand — it rode in entirely on the single $e$ inside $D_\mu$ (the same $e$ fixing the gauge shift $A_\mu\to A_\mu-\tfrac1e\partial_\mu\alpha$), so the photon couples to the current with strength exactly $e$, no other value allowed. This is bridge B8: local U(1) symmetry has generated the interaction. The check isolates the mechanism — the coupling sits in the $ieA_\mu\psi$ piece, and setting $e=0$ collapses $D_\mu\psi$ to the free $\partial_\mu\psi$: covariant_derivative(psi, A, 0.0, dx) equals np.gradient(psi, dx, axis=mu) to $\approx10^{-14}$.

(a) Under $A_\mu\to A_\mu-\tfrac1e\partial_\mu\alpha$ a photon mass term transforms as

since the cross and $(\partial\alpha)^2$ pieces do not cancel ($A_\mu$, unlike the gauge-invariant $F_{\mu\nu}F^{\mu\nu}$, is not invariant). So gauge symmetry forbids a photon mass: $m_\gamma=0$ exactly. (b) The static potential is the Fourier transform of the propagator $1/(\mathbf q^2+m^2)$, the Green's function of $(-\nabla^2+m^2)$:

A massive mediator screens the force over range $1/m$; the massless photon of (a) gives the infinite-range $1/r$ Coulomb law of ~EM-01, so finding electrostatics to be exactly long-range is a bound on $m_\gamma$. The check confirms $V/V_{\text{Coulomb}}=e^{-mr}$: at $m=r=1$ it is $e^{-1}\approx0.3679$, rising monotonically to $1.0$ as $m\to0$, and $V=e^{-mr}/4\pi r$ solves $(-\nabla^2+m^2)V=0$ for $r>0$.

An external photon enters an amplitude through its polarization, $\mathcal M=\epsilon_\mu(k)\mathcal M^\mu(k)$. A gauge transformation shifts the polarization along the photon momentum,

the momentum-space image of $A_\mu\to A_\mu+\partial_\mu\lambda$. Physics cannot depend on the gauge parameter $c$, so the extra piece must not contribute: $\mathcal M\to\mathcal M+c\,k_\mu\mathcal M^\mu$ for all $c$ forces

the Ward identity. The longitudinal ($k_\mu$-direction) polarization decouples, so of the four components of $\epsilon_\mu$ only the two transverse ones are physical — the massless photon has two helicities. This is the amplitude-level shadow of current conservation $\partial_\mu j^\mu=0$ (in momentum space $k_\mu j^\mu=0$). The check shows the gauge mode is inert already classically: a pure-gauge $A_\mu=\partial_\mu\lambda$ — the "$k_\mu$ direction" — carries zero field strength, field_strength(gauge_transform(0*A, lam, dx), dx) $\approx1.7\times10^{-15}$, numerically zero.

Power-count the integrand at large loop momentum $k$: the measure $d^4k\sim k^3\,dk$ scales as $k^4$ and the two propagators as $(k^2-m^2)^{-2}\sim k^{-4}$, so

a logarithmic divergence (superficial degree $D=4-4=0$). Carrying out the integral,

A coupling counterterm $\delta_\lambda=-\tfrac{3\lambda^2}{32\pi^2}\ln(\Lambda^2/\mu^2)$ cancels the cutoff dependence, leaving a finite renormalized vertex. The theory is renormalizable because $\lambda$ is dimensionless in $d=4$: only the 2- and 4-point functions have $D\ge0$, so finitely many counterterms suffice at every order. The cutoff-independent coefficient of $\ln(\Lambda/\mu)$ is exactly the beta-function slope $\beta(\lambda)=3\lambda^2/16\pi^2$, matching beta_phi4(0.1) $=3(0.1)^2/16\pi^2\approx1.90\times10^{-4}$.

Holding the bare coupling fixed, $\mu$-independence of physics forces the renormalized $\alpha(\mu)$ to run; the one-loop vacuum-polarization log gives

Separate variables, $d\alpha/\alpha^2=(2/3\pi)\,d\ln\mu$, and integrate from $\mu$ to $Q$:

The positive $\beta$ makes $1/\alpha$ decrease with energy, so $\alpha$ is largest in the UV and smallest in the IR: the $e^+e^-$ vacuum polarizes around a bare charge and screens it, weakening the effective charge seen from afar. The check confirms both: beta_qed(1/137.036) $\approx1.13\times10^{-5}>0$, and $\alpha$ runs up monotonically — $\alpha(m_e)=1/137.0\to\alpha(1\,\text{GeV})=1/135.4\to\alpha(M_Z)=1/134.5$ (electron-only).

With only the electron in the loop ($S=1$) the closed form gives

Each charged fermion $f$ contributes the same way once $Q>m_f$, weighted by $Q_f^2N_c^f$ (squared charge $\times$ colour), steepening the slope:

Summing $e,\mu,\tau$ and the five quarks ($u,d,s,c,b$ with $N_c=3$) brings $1/\alpha(M_Z)$ down to $\approx128$, the measured value. The hadronic part is not a clean free-quark calculation because near $m_u,m_d$ the strong coupling is large and the quarks are confined; the genuine low-energy hadronic vacuum polarization is extracted from $e^+e^-\to$ hadrons data (the $R$-ratio), the $\sim0.3$ GeV light-quark "mass" being only a stand-in. The check matches: 1/qed_running_alpha(M_Z) $\approx134.5$ (electron only) and 1/qed_alpha_sm(M_Z) $\approx128.4$ (all charged fermions).

The running coupling $\alpha(Q)=\alpha_0/[1-(\alpha_0/3\pi)\ln(Q^2/Q_0^2)]$ blows up when its denominator vanishes. Setting it to zero,

and $\ln(Q^2/Q_0^2)=2\ln(Q/Q_0)$ gives the Landau pole

With $\alpha_0\approx1/137$ and $Q_0=m_e$ this is $\approx1.45\times10^{277}$ GeV — vastly above the Planck scale $\sim10^{19}$ GeV. Physically it does not mean the charge truly becomes infinite: long before such energies gravity and new physics intervene and the one-loop formula has broken down. The pole only signals that perturbative QED cannot be complete up to arbitrarily high scales. The check gives qed_landau_pole() $\approx1.45\times10^{277}$ GeV, the scale at which the closed-form denominator vanishes.

Separate variables in $d\lambda/d\ln\mu=3\lambda^2/16\pi^2$:

Since $1/\lambda$ falls as $\mu$ grows, $\lambda(\mu)$ rises in the UV and diverges when $1/\lambda\to0$, at the Landau pole

A genuine continuum theory needs the cutoff $\Lambda\to\infty$ with the pole pushed to infinity, which from the formula requires $1/\lambda_0\to\infty$, i.e. $\lambda_0\to0$: the renormalized coupling is driven to zero and the interacting theory collapses to a free one — the triviality of $\phi^4$ in four dimensions. The check confirms the integration: phi4_running(1.0, 1.0, 1e10) $\approx1.7776$, matching the closed form, and phi4_landau_pole(1.0, 1.0) $\approx7.25\times10^{22}$ GeV.

The non-abelian one-loop beta function is

The first term $+\tfrac{11}{3}C_A$ comes from the gluon self-interaction (gluons carry colour charge, with no abelian analogue) and antiscreens; the second $-\tfrac43 T_Fn_f$ is the QED-like fermion screening. For QCD the gluon term wins whenever $n_f<33/2=16.5$, so the bracket is positive and $\beta<0$. Equivalently $b_0=11-\tfrac23 n_f>0$, and the closed form

shows $1/\alpha_s$ growing with $Q$, so $\alpha_s$ shrinks in the UV — asymptotic freedom, why quarks are nearly free at short distance. This is the exact opposite of QED (P2), where $\beta>0$ makes $\alpha$ grow with energy. The check confirms the sign and the running: beta_qcd(0.1181) $\approx-0.017<0$, and $\alpha_s(2\,\text{GeV})=0.263>\alpha_s(M_Z)=0.118>\alpha_s(1\,\text{TeV})=0.088$.

Factor the beta function, $\beta(g)=ag^2-bg^3=g^2(a-bg)$, so the zeros are $g=0$ and

Stability follows from the slope $\beta'(g)=2ag-3bg^2$. Near the origin $\beta\approx ag^2>0$ for small $g>0$, so the coupling flows away from $0$ as $\mu$ rises — UV-repulsive. At $g_*$,

a negative slope, so nearby couplings flow into $g_$ as $\mu\to\infty$ — UV-attractive. This linearized slope $\beta'(g_)$ is the RG eigenvalue at the fixed point; at the Wilson–Fisher fixed point of $\phi^4$ in $d=4-\epsilon$ it fixes the critical exponents (e.g. the correlation-length exponent $\nu$) of the Ising universality class — the bridge to ~SM-05. The check matches: fixed_point(2.0, 4.0) returns g_star $=0.5$ and beta_prime $=-1.0$ with stability $=$ "UV-attractive", and beta_toy(0.5, 2.0, 4.0) $\approx0$.

Expand the adjoint $a_{j'}^{\prime\dagger}=\sum_l(\alpha_{j'l}\,a_l^\dagger-\beta_{j'l}\,a_l)$ and commute it against $a'_j$ using $[a_k,a_l^\dagger]=\delta_{kl}$ (all other commutators vanish):

Setting this to $\delta_{jj'}$ and taking $j=j'$ gives the normalization $\sum_k(|\alpha_{jk}|^2-|\beta_{jk}|^2)=1$. A single squeezed mode $\alpha=\cosh r,\ \beta=\sinh r$ satisfies it automatically, $\cosh^2 r-\sinh^2 r=1$. Its rotated vacuum holds $\langle N\rangle=\lVert a'|0\rangle\rVert^2$; since $a'|0\rangle=(\alpha^a-\beta^a^\dagger)|0\rangle=-\beta^*|1\rangle$,

At $r=1.3$, $(\cosh r,\sinh r)=(1.971,1.698)$, so bogoliubov_check$=1.971^2-1.698^2=1.0$ and particle_number$=\sinh^2 1.3\approx2.884$.

In $T_U=\dfrac{\hbar a}{2\pi c\,k_B}$ every factor but $a$ is a constant, so $T_U\propto a$ exactly: $a\mapsto\lambda a$ sends $T_U\mapsto\lambda T_U$ (the check is $\lambda=2$). Invert for the acceleration that reaches $1\,$K,

The prefactor $2\pi c\,k_B/\hbar\approx2.47\times10^{20}$ is enormous (because $\hbar$ is tiny and $c$ huge), so even $1\,$K demands an acceleration $\sim10^{19}g$ — hopelessly large. Hence $T_U(10^{20})=10^{20}/2.47\times10^{20}\approx0.405\,$K, while everyday $g=9.81$ gives $T_U=9.81/2.47\times10^{20}\approx4.0\times10^{-20}\,$K, undetectable. These match unruh_temperature(1e20)$\approx0.4055$ and unruh_temperature(9.81)$\approx3.98\times10^{-20}$, with the ratio $2.0$.

Write $R\equiv|\beta_\omega|^2/|\alpha_\omega|^2=e^{-2\pi c\omega/a}$. The norm $|\alpha_\omega|^2-|\beta_\omega|^2=1$ becomes $|\alpha_\omega|^2(1-R)=1$, so

The exponent is exactly the Boltzmann factor at $T_U$, since $\dfrac{\hbar\omega}{k_BT_U}=\hbar\omega\cdot\dfrac{2\pi c\,k_B}{k_B\,\hbar a}=\dfrac{2\pi c\omega}{a}$; thus $\langle n_\omega\rangle=1/(e^{\hbar\omega/k_BT_U}-1)$, a Planck spectrum. And since $|\alpha|^2=|\beta|^2+1$, $R=|\beta|^2/|\alpha|^2=\langle n\rangle/(\langle n\rangle+1)$. At $\omega=a/2\pi c$ the exponent is $1$: $\langle n\rangle=1/(e-1)=0.582$ and $R=e^{-1}=0.368=0.582/1.582$. So unruh_occupation(w, a) equals bose_occupation(w, unruh_temperature(a)) and unruh_beta_ratio(w, a)$=\langle n\rangle/(\langle n\rangle+1)$.

From $r_s=2GM/c^2$, the surface gravity is $\kappa=\dfrac{c^2}{2r_s}=\dfrac{c^2}{2}\cdot\dfrac{c^2}{2GM}=\dfrac{c^4}{4GM}$. Hawking's temperature is the Unruh formula evaluated at the horizon's acceleration $a\to\kappa$:

For the Sun, $r_s\approx2953\,$m and $\kappa=c^2/2r_s\approx1.52\times10^{13}\,$m/s$^2$, giving $T_H\approx6.17\times10^{-8}\,$K. This is $\sim8$ orders of magnitude below the $2.7\,$K CMB, so the hole absorbs far more radiation than it emits — effectively absorbing, it grows rather than evaporates. Matches schwarzschild_radius(M_SUN)$\approx2953$, surface_gravity(M_SUN)$\approx1.52\times10^{13}$, and hawking_temperature(M_SUN)$\approx6.17\times10^{-8}=\hbar\cdot$surface_gravity$/2\pi c\,k_B$.

Since $T_H=\hbar c^3/8\pi GM\,k_B\propto1/M$, scale off the solar value:

Its horizon is $r_s=2G/c^2\approx1.49\times10^{-27}\,$m — about twelve orders below a proton radius. At $k_BT_H$ of order $10^{23}\,$K the thermal quanta exceed every Standard-Model rest mass, so the hole pours out radiation; and because $\dot M\propto-1/M^2$, shedding mass only raises $T_H$ — a runaway ending in an explosive final burst. Matches hawking_temperature(1.0)$\approx1.23\times10^{23}$ K and schwarzschild_radius(1.0)$\approx1.5\times10^{-27}$ m (sub-nuclear).

Both factors scale with mass: $A=4\pi r_s^2=4\pi(2GM/c^2)^2\propto M^2$ and $T_H^4\propto M^{-4}$, so $L=\sigma A T_H^4\propto M^{-2}$. With every constant collected, the balance $L=-c^2\dot M$ reads

Integrating from $M$ down to $0$ with $\int_0^M M'^2\,dM'=M^3/3$,

so $\tau(2M)=8\,\tau(M)$. For $M_\odot$, $\tau\approx6.6\times10^{74}\,$s$\,\approx2.1\times10^{67}\,$yr, dwarfing the $\sim1.4\times10^{10}\,$yr age of the universe; setting $\tau$ equal to that age inverts to $M\sim10^{11}\,$kg — the primordial holes expiring now. Matches evaporation_lifetime(M_SUN)$\approx6.6\times10^{74}$ s and the ratio $8.0$.

Let $x=\hbar\omega/k_BT>0$ for any real mode at $T>0$. The two builders are $n_B=1/(e^x-1)$ and $n_F=1/(e^x+1)$. Because $e^x-1<e^x+1$, the bosonic denominator is the smaller one, so

bosons pile up while Pauli exclusion throttles fermions. Moreover $e^x+1>2$ for all $x>0$, hence

rising toward $\tfrac12$ as $x\to0^+$ and falling to $0$ as $x\to\infty$, i.e. $n_F\in(0,\tfrac12]$. For instance at $x\approx0.076$, $n_B=12.6\gg n_F=0.481<\tfrac12$, and $n_F(x\to0)=\tfrac12$. So thermal_beta_squared(w, T, 'bose') $=$ bose_occupation(w, T) exceeds thermal_beta_squared(w, T, 'fermi') for every $w$, with fermi_occupation bounded in $(0,\tfrac12]$.