Evolution in practice
Mellin solutions, small-x growth, heavy-flavour schemes, and the EKO/FK pipeline.
Three chapters of theory now meet the machine room. DGLAP is a system of integro-differential equations; a Bayesian fit needs its solution — for every trial PDF, at every dataset's $Q^2$ — hundreds of thousands of times. This chapter shows how the system is actually solved: exactly in Mellin space (where the non-singlet solution is one line and the singlet is a matrix exponential), numerically on $x$-grids (what EKO does), and then — the decisive engineering move — precomputed once into FK tables so that a fit never evolves anything at all. Along the way: heavy-quark thresholds and what "perturbative charm" means in our runs, why evolving downward is treacherous, where DGLAP itself stops being the right equation, and how missing higher orders become a covariance matrix.
Solving DGLAP I — the non-singlet solution in closed form
In $N$-space (§2.3) each non-singlet moment obeys an ordinary differential equation, but with a running coefficient:
$$ \frac{d\, \tilde q(N, Q^2)}{d \ln Q^2} \;=\; \frac{\alpha_s(Q^2)}{2\pi}\, \gamma_{qq}(N)\; \tilde q(N, Q^2) $$
Trade the awkward variable $\ln Q^2$ for $\alpha_s$ itself using the one-loop RGE $d\alpha_s/d\ln Q^2 = -\beta_0 \alpha_s^2$ (§2.1):
$$ \frac{d \ln \tilde q}{d \alpha_s} \;=\; \frac{(\alpha_s/2\pi)\,\gamma_{qq}(N)}{-\beta_0\, \alpha_s^2} \;=\; -\,\frac{\gamma_{qq}(N)}{2\pi \beta_0}\, \frac{1}{\alpha_s} $$
which integrates immediately to the exact LO solution:
$$ \tilde q(N, Q^2) \;=\; \tilde q(N, Q_0^2) \left[ \frac{\alpha_s(Q_0^2)}{\alpha_s(Q^2)} \right]^{\gamma_{qq}(N)/(2\pi\beta_0)} $$
All the $Q^2$-dependence rides on the ratio of couplings raised to a fixed power per moment. Sanity checks: $\gamma_{qq}(1) = 0$ gives exponent zero — valence number frozen; $\gamma_{qq}(N) < 0$ for $N \ge 2$ with $\alpha_s(Q_0^2)/\alpha_s(Q^2) > 1$ gives shrinking higher moments — large-$x$ depletion. Worked example 1 runs the numbers. Beyond LO the same structure survives with corrections organized in powers of $\alpha_s(Q_0) - \alpha_s(Q)$; the principle — evolution is a function of the coupling ratio — stands.
Solving DGLAP II — the singlet is a matrix problem
The singlet pair $(\Sigma, g)$ evolves with the $2\times 2$ anomalous-dimension matrix $\Gamma(N) = \begin{pmatrix} \gamma_{qq} & 2n_f\,\gamma_{qg} \\ \gamma_{gq} & \gamma_{gg} \end{pmatrix}$. At LO, diagonalize once per moment: eigenvalues
and the solution is a sum of two non-singlet-style powers along the eigenvector projectors — the "+" eigenvalue owns the small-$x$ gluon growth of §2.3. Beyond LO a genuine subtlety appears: $\Gamma$ changes with scale (through $\alpha_s$) and matrices at different scales do not commute, so the formal solution is a path-ordered exponential $\mathcal{P}\exp\int \Gamma\, d\ln Q^2$ — no closed form. Every serious evolution code therefore goes numerical: discretize the scale integration (Runge–Kutta or a truncated expansion in $\alpha_s$), and handle the convolution either by staying in $N$-space and inverting at the end, or on an $x$-grid directly.
What EKO computes: evolution as a precomputed linear map
The key structural fact: DGLAP is linear in the PDFs. Whatever the numerical scheme, the map from input to output is therefore a linear operator, independent of which PDF it acts on. EKO exploits this completely: it solves the RGE in Mellin space (exact anomalous dimensions, path-ordering handled by controlled approximants), projects onto a Lagrange-interpolation $x$-grid, and outputs the evolution kernel operator
— a big matrix over (flavour $\times$ grid-node) indices, computed once per theory settings and stored. Composing $E$ with the process's perturbative coefficients (the PineAPPL grid of $C_i$ values, §2.2) and collapsing the scale sums yields the FK table: for DIS a single linear map
and for hadronic observables the quadratic analogue $T_d = \sum \mathrm{FK}_{d,\,ia,\,jb} f_i(x_a) f_j(x_b)$. Everything upstream is frozen into the table: perturbative order, $\alpha_s$, thresholds, evolution. In our pipeline that freezing is literal and versioned — theory id 40000000 pins NNLO, $\alpha_s(M_Z) = 0.118$, $Q_0 = 1.65$ GeV and the heavy-quark masses, so every fit from PPDF-1 to PPDF-15 confronts identical theory, and any difference between results is attributable to data, parametrization, or statistics — never to a drifting evolution code.
Nested sampling (§4.2) evaluates the likelihood $\sim 10^5$–$10^6$ times per run. Solving DGLAP from
scratch takes EKO minutes per theory; done inside the loop, a single fit would take years. As a precomputed
contraction, a likelihood call is a handful of small einsums — microseconds to milliseconds
(worked example 2 counts the flops). The factor $\sim 10^8$ between the two architectures is the difference
between "Bayesian PDF inference is impossible" and the routine of this lab. The price is rigidity: change
$\alpha_s$, the order, or a mass, and every FK table must be regenerated — which is exactly why theory ids
exist, and why our MSHT20-baseline comparisons insist on matching them
before comparing anything else.
Heavy flavours: schemes, matching, and what "perturbative charm" means
Charm ($m_c = 1.51$ GeV) and bottom ($m_b = 4.92$ GeV) sit inside our kinematic range, and there are two consistent ways to treat a quark of mass $m_h$. In a fixed flavour number scheme (FFNS) the heavy quark is never a parton: it appears only in coefficient functions, with full mass dependence — accurate near threshold ($Q^2 \sim m_h^2$), but at $Q^2 \gg m_h^2$ the uncollected logs $\alpha_s^n \ln^n(Q^2/m_h^2)$ wreck the expansion. In a variable flavour number scheme (VFNS) the heavy quark becomes a parton above its threshold, the logs are resummed by DGLAP, and the $n_f$- and $(n_f{+}1)$-flavour descriptions are stitched together by perturbative matching conditions,
where the $A_{ij}$ are computable; at our NNLO they are non-trivial, generating a small but non-zero charm distribution at threshold. General-mass VFNS variants (FONLL, TR′, ACOT) interpolate the two limits and are what production fits use. Now the physics distinction that matters for us. Perturbative charm: the charm PDF is not fitted — it is whatever matching plus evolution generate. Our fits implement this by parametrizing an $n_f{=}4$-flavour proton with $c(x, Q_0^2) = 0$ at $Q_0 = 1.65$ GeV, deliberately parked just above $m_c = 1.51$: charm then grows radiatively ($g \to c\bar c$) and by construction carries no free parameters. Fitted charm: parametrize $c(x, Q_0^2)$ like any other flavour and let data decide — the physical question being whether the proton has intrinsic charm, a non-perturbative $|uudc\bar c\rangle$ component (expected at $\mathcal{O}(\Lambda^2/m_c^2)$, peaked at large $x$), which NNPDF4.0-style analyses report as a $\sim\!3\sigma$-level valence-like bump at $x \sim 0.4$. That is the RUN-15 direction in this lab's ledger: same machinery, one more fitted flavour, and the Bayesian evidence (§4.3) as the referee for whether the data demand it.
Two honest caveats: backward evolution, and the edge of DGLAP
Backward evolution is unstable. Upward evolution is diffusive — it smooths: radiation spreads momentum, structure gets washed toward the universal small-$x$ rise. Running the same operator backward is anti-diffusion: features sharpen, and any small error — interpolation noise, a percent-level wiggle in the input — amplifies exponentially, often driving PDFs negative or oscillatory (backward matching across a threshold, where the operator must be inverted, is especially delicate). The practical consequences shape our whole setup: parametrize at the lowest scale you trust ($Q_0 = 1.65$ GeV) and only ever evolve up; and treat any "de-evolved" PDF with suspicion — it is the numerical analogue of unsmearing.
DGLAP is not the whole truth at small $x$. DGLAP resums $(\alpha_s \ln Q^2)^n$; at small enough $x$ the competing logarithm $\alpha_s \ln(1/x)$ is just as large, and its resummation is a different equation — BFKL — with its own prediction (a power-like $x^{-\omega}$ growth at fixed $Q^2$; the LL intercept $\omega \approx 0.5$ is far too steep, and NLL corrections to it are famously violent, which is the field's warning that this frontier is genuinely hard). State of the art combines both: small-$x$-resummed fits find that HERA data below $x \sim 10^{-4}$, $Q^2 \lesssim 10$ GeV$^2$ mildly prefer resummation. For this lab's current fixed-target-dominated datasets ($x \gtrsim 10^{-2}$) pure DGLAP is comfortably adequate; the caveat is filed for the day HERA's lowest-$x$ bins enter our likelihoods.
And the orders we did not compute remain an uncertainty even inside DGLAP's domain. The standard estimator is scale variation — vary $\mu_R$ and $\mu_F$ by factors of 2 around their central choices (the 7-point envelope), since the residual dependence is exactly the size of the first missing order (§2.2, Derivation 2). Promoted to a correlated theory covariance matrix added to the experimental one, this is the MHOU programme — and in the FPPDF direction this lab is joining, those theory-error degrees of freedom enter the Bayesian fit itself, alongside the PDF parameters.
Worked example 1 — evolving a moment by hand
Given: a non-singlet (valence-type) moment with $N = 2$: $\gamma_{qq}(2) = -\tfrac{16}{9}$ (§2.3). Work at $n_f = 4$ throughout for transparency ($\beta_0 = \tfrac{25}{12\pi}$, so $2\pi\beta_0 = \tfrac{25}{6}$; EKO would switch to $n_f = 5$ at $Q^2 = m_b^2 = 24.2$ GeV$^2$ — a percent-level refinement here). One-loop couplings from the §2.1 chain: $\alpha_s(2\ \text{GeV}^2) = 0.307$, $\alpha_s(100\ \text{GeV}^2) = 0.171$.
Exponent: $\;\gamma_{qq}(2)/(2\pi\beta_0) = \dfrac{-16/9}{25/6} = -\dfrac{32}{75} = -0.4267$.
Evolve: $$ \frac{\tilde q(2, 100)}{\tilde q(2, 2)} = \left[\frac{0.307}{0.171}\right]^{-0.4267} = (1.796)^{-0.4267} = e^{-0.4267 \times 0.5857} = e^{-0.2499} = 0.779 $$
Reading: the momentum carried by any valence combination falls by $22\%$ over this range — if $u_V$ carried $28\%$ of the proton's momentum at $2$ GeV$^2$, it carries $\approx 21.8\%$ at $100$ GeV$^2$, the difference radiated into glue and sea exactly as the paired second-moment sum rules of §2.3 dictate. Note what made this trivial: in $N$-space, evolution is one exponentiation. The hard part of real life is only that data live in $x$-space, at many $Q^2$, with thresholds — which is precisely the bookkeeping EKO's operator buries into the FK tables.
Worked example 2 — one parameter vector, two $Q^2$: tracing the tensors
Step 0 — parameters: a trial point from the sampler: $\theta \in \mathbb{R}^{36}$ (our MSHT20-style parametrization; PPDF-10 onward). Sum-rule normalizations are computed analytically from $\theta$ (§3.2).
Step 1 — PDF vector at the input scale: evaluate the parametrization on the interpolation grid: $f_i(x_a, Q_0^2)$ with $i = 1\ldots 14$ (evolution basis: $\Sigma, g, V, V_3, \ldots$) and $a = 1\ldots 50$ ($x$-nodes, log-spaced toward small $x$). Shape: $(14, 50)$ — the 700 numbers that are the only PDF information the likelihood ever sees.
Step 2 — dataset A (say BCDMS $F_2^p$ bins at $Q^2 \approx 12$ GeV$^2$, $N_A$ points): its FK
table has shape $(N_A, 14, 50)$, with evolution $Q_0^2 \to 12$ GeV$^2$ already folded in. Prediction:
$T^A_d = \sum_{ia} \mathrm{FK}^A_{d,ia} f_{ia}$ — einsum("dia,ia->d").
Step 3 — dataset B (bins at $Q^2 \approx 110$ GeV$^2$): a different table, shape $(N_B, 14, 50)$, containing the longer evolution — contracted with the same $(14, 50)$ vector. Nothing about $\theta$ or $f$ changes between Steps 2 and 3: all $Q^2$-dependence, thresholds included, lives in the tables.
Step 4 — cost: per data point, $14 \times 50 = 700$ multiply–adds; for $\sim 10^3$ DIS points, $\sim 10^6$ flops per likelihood call; across $5 \times 10^5$ nested-sampling calls, $\sim 5 \times 10^{11}$ flops — under a minute of a single modern core. The same fit with in-loop evolution ($\sim$ minutes per call) would need years. That ratio is the entire argument for the EKO $\to$ FK architecture — and the reason a Q for Maria about theory settings (order, $\alpha_s$, masses) is always also a question about which FK tables to regenerate.