Sum rules analytically: the Euler-Beta machinery

Every normalization integral derived — the exact algebra our MSHT20 port runs.

Four of MSHT20's normalizations are not fitted — they are computed, exactly, at every single sampler step, from the other 52 parameters. This chapter derives the machinery that makes that possible: one classical special function (the Euler Beta function), one master integral, and a pile of careful bookkeeping for the Chebyshev terms. Everything here is the mathematics behind msht20_jax.py, the port whose parity gate passed at the $10^{-9}$ level (PPDF-10).

Why analytic beats penalties

The sum rules of Chapter 1.4 — valence counting $\int_0^1 u_V\,dx = 2$, $\int_0^1 d_V\,dx = 1$, zero net strangeness $\int_0^1 (s - \bar s)\,dx = 0$, and momentum $\int_0^1 x \sum_i f_i(x)\,dx = 1$ — are exact consequences of the proton's quantum numbers. There are two ways to impose them on a parametrized model. The soft way adds a penalty $\lambda\,(\text{sum rule violation})^2$ to the $\chi^2$: now the rules hold only approximately, the answer depends on an arbitrary $\lambda$, and the sampler still wanders the violating directions. The hard way solves the constraints analytically: express $A_{u_V}, A_{d_V}, A_g, x_0$ in closed form in terms of the remaining parameters and substitute. The rules then hold to machine precision at every point ever evaluated; four dimensions vanish from the sampling space (52 instead of 56 — nested sampling cost scales roughly linearly in dimension, and flat normalization directions are exactly the kind that wreck slice-sampler efficiency); and no numerical quadrature runs inside the likelihood. The only requirement is that the parametrization's integrals be doable in closed form. For powers times polynomials in $\sqrt{x}$, they are — that is the Beta function's job.

One more reason, specific to this lab's evidence program: $\log Z$ comparisons require every model to be normalized on the same footing. A penalty term is an ad-hoc modification of the likelihood whose strength $\lambda$ shifts $\log Z$ by an arbitrary amount — two models with different penalty weights are not comparable. Exact constraints change nothing about the likelihood; they only restrict the hypothesis space, which is precisely what a model is. There is also a subtler statistical effect: once $A_{u_V}$ is a derived quantity, its posterior uncertainty comes entirely from the shape parameters it depends on — the sum rules induce exact correlations between normalizations and shapes, which is physics (a harder $u_V$ must be taller to hold two quarks), not an artifact.

The Euler Beta function and the master integral

The momentum sum rule assembles $A_g$

Momentum is carried by all quarks plus the gluon. In the MSHT basis (with the FPPDF convention $S = 2(\bar u + \bar d) + s + \bar s$, so $u_V + d_V + S + c_+$ exhausts the quarks) the pieces are the $\delta$-index momentum integrals, and the two-term gluon contributes two integrals of its own:

with $\langle x \rangle_q = A_q [ I(\delta_q, \eta_q) + \sum_i a_{q,i} Ic_i(\delta_q, \eta_q)]$ — where $A_{u_V}, A_{d_V}$ are themselves the counting-rule outputs, so the order of evaluation matters: valence normalizations first, then momentum. Only the first gluon term's normalization is solved for; the second term (that negative small-$x$ dip) is fitted freely and enters as the constant $g_2$.

Zero strangeness fixes the crossing point $x_0$

$x s_- = A_- x^{\delta_-}(1-x)^{\eta_-}(1 - x/x_0)$ must integrate (as a number density) to zero:

— the crossing point is a pure ratio of two Beta functions (Chebyshev-dressed in general; MSHT freeze those coefficients to zero). At the publication values $\delta_- = 0.0015$, $\eta_- = 7.46$ this gives $x_0 = 1.76 \times 10^{-4}$: with $\delta_-$ that close to zero, $\mathcal{I}_1 = I(-0.9985, 7.46)$ is nearly divergent ($\Gamma(0.0015) \approx 670$), so the number integral is dominated by tiny $x$ and the crossing must sit far down to balance it. The amplitude $A_-$ stays free — the sum rule fixes where $s - \bar s$ changes sign, not how big it is.

Numerical stability: the Stirling branch

One engineering wrinkle. Float64 $\Gamma$ overflows beyond $\Gamma(171.6) \approx 10^{308}$, and even before that, ratios of huge Gammas lose digits. Large $\eta$ is physical (sea $\eta \sim 7$ is fine, but samplers explore, and other FPPDF configurations run $\eta$ large), so the code switches at $b \ge 100$ to a Stirling-type closed form,

which needs only $\Gamma(a+1)$ (small argument) and powers — no overflow, relative error $\sim 10^{-5}$ at the $b = 100$ switch point and shrinking as $b$ grows. Our port reproduces the branch and its threshold verbatim: parity means matching the reference's choices, including its approximations.