DGLAP: derivation and splitting functions
P_qq derived at LO, plus-distributions, the full LO kernel set, moments.
Factorization handed us an evolution equation as a consistency condition (§2.2, Derivation 2); this chapter earns it as physics and computes its kernels. The core idea fits in one sentence: how many partons you count depends on how closely you look, and the change per decade of resolution is perturbatively calculable. We derive $P_{qq}$ from collinear gluon emission, tame its soft divergence with the plus-distribution (and see that the delta term is conservation of quark number, not a trick), assemble the full coupled DGLAP system — meeting the reason the evolution basis of our pipeline exists — and turn the whole apparatus into algebra with Mellin moments.
Resolution-dependent counting
A quark is never bare: quantum mechanics dresses it in a cloud of virtual gluons and $q\bar q$ pairs. A probe of virtuality $Q^2$ resolves transverse structure down to $\sim 1/Q$: fluctuations with transverse momentum $k_T < Q$ are seen as separate partons (each carrying a share of the momentum), while those with $k_T > Q$ stay unresolved inside the effective quark. Raise $Q^2$ and previously hidden splittings come into view: the counted quark loses momentum to newly-visible gluons, and the small-$x$ population swells. Because the emission probability per logarithmic interval of $k_T^2$ is scale-invariant ($\alpha_s\, dk_T^2/k_T^2$ — no scale in the problem between $\Lambda^2$ and $Q^2$), the change is uniform per decade: PDFs evolve in $\ln Q^2$. DGLAP is the bookkeeping of this, and its kernels $P_{ji}(z)$ answer one question each: given a parton $i$, what is the probability of finding inside it a parton $j$ with momentum fraction $z$, per unit $\alpha_s/2\pi$ per unit $\ln k_T^2$?
Take a quark of momentum $p$ (call its large light-cone component the reference) emitting a gluon and retaining fraction $z$: $q(p) \to q(zp + k_T\text{-recoil}) + g\big((1-z)p, -k_T\big)$, with $k_T \ll$ all large scales. Kinematics first: the intermediate quark's virtuality is
$$ t \;=\; (p - k_g)^2 \;\simeq\; -\,\frac{k_T^2}{1-z} $$
— it vanishes as $k_T \to 0$: the propagator $1/t$ blows up, and the emission becomes long-distance physics. Now assemble the rate. The amplitude carries the vertex $g_s$, colour factor $C_F$ after averaging and summing, and the propagator $1/t$; the squared amplitude thus scales as $g_s^2\, C_F\, |M_{\rm split}|^2 / t^2$, where $M_{\rm split}$ is the spin-dependent numerator. The two-body phase space in the collinear region supplies $dz\, dk_T^2 \times$ constants, and one power of the numerator cancels one $1/t$: schematically,
$$ d\sigma \;=\; \sigma_{\rm hard}\;\times\; \frac{\alpha_s}{2\pi}\, \frac{dk_T^2}{k_T^2}\; P_{qq}(z)\, dz $$
The spin algebra (average over initial helicities, sum over final — half a page of Dirac traces we compress, honestly flagged) leaves the celebrated kernel:
$$ P_{qq}^{\rm real}(z) \;=\; C_F\, \frac{1 + z^2}{1 - z} $$
Read its anatomy. The $\frac{1}{1-z}$ is the eikonal factor: soft gluons ($z \to 1$, gluon energy $\to 0$) couple to the quark's charge, not its structure, and their emission diverges logarithmically — every gauge theory does this. The numerator $1 + z^2$ is helicity conservation: a vector coupling to (nearly) massless quarks preserves the quark helicity, and $1+z^2 = $ (no spin flip)$^2$ summed over the two transverse gluon polarizations. And the $\frac{dk_T^2}{k_T^2}$, integrated from a regulator $m^2$ up to the probe scale $Q^2$, gives $\ln(Q^2/m^2)$ — the very same collinear logarithm that mass factorization filed into the PDF in §2.2. DGLAP evolution is nothing but these logarithms, emission by emission, with $k_T$-ordering ($k_{T,1}^2 \ll k_{T,2}^2 \ll \ldots \ll Q^2$) organizing the leading tower $\big(\alpha_s \ln Q^2\big)^n$.
$P_{qq}^{\rm real}$ is not integrable at $z = 1$. Physics fixes this: real emission is not the whole $\mathcal{O}(\alpha_s)$ story — there are virtual corrections (gluon emitted and reabsorbed), which live entirely at $z = 1$ (the quark keeps all its momentum) and carry the opposite-sign soft divergence. To combine them, define the plus-distribution by its action under an integral:
$$ \int_0^1 dz\; \frac{f(z)}{(1-z)_+} \;\equiv\; \int_0^1 dz\; \frac{f(z) - f(1)}{1 - z} $$
— finite for any smooth $f$, since the numerator vanishes linearly at $z = 1$. The subtracted piece ($-f(1)$) is exactly the virtual contribution's job. Rather than compute the virtual diagrams, use unitarity: emitting a gluon does not change the number of quarks, so the total kernel must satisfy $\int_0^1 P_{qq}(z)\, dz = 0$ (what leaves the momentum bin $z=1$ must reappear at $z < 1$, with zero net quark creation). Writing $P_{qq}(z) = C_F\big[\tfrac{1+z^2}{(1-z)_+}\big] + A\,\delta(1-z)$ and imposing the sum rule fixes $A$ (worked example 1 does the integral): $A = \tfrac32 C_F$. Hence the complete LO kernel:
$$ P_{qq}(z) \;=\; C_F\left[ \frac{1 + z^2}{(1-z)_+} \;+\; \frac{3}{2}\,\delta(1-z) \right] $$
The delta term is not a regularization artefact: it is the depletion of unsplit quarks, sized by conservation of probability.
The four kernels
The same collinear analysis for each splitting channel gives the complete LO set:
Read them like a field guide. $P_{qg}$ ($g \to q\bar q$): finite everywhere — no soft divergence because a soft quark costs a power (only gauge bosons eikonalize); symmetric under $z \leftrightarrow 1-z$ since quark and antiquark are twins. $P_{gq}$ ($q \to g$): the $1/z$ is the eikonal factor seen from the gluon's side — soft gluons are abundant, and this term feeds the small-$x$ sea. $P_{gg}$: has both divergences, $1/(1-z)_+$ and $1/z$ (a gluon is soft-singular as parent and as child), with the largest colour charge $2C_A$ in front — the gluon is the most radiative object in the theory, and its $2C_A/z$ term drives everything at small $x$ (worked example 2). And a pleasing echo: the $\delta$-coefficient is $\tfrac{11C_A - 4n_f T_R}{6} = 2\pi\beta_0$ — the beta function of §2.1 reappearing as the gluon's probability bookkeeping, gluon splitting against $q\bar q$ dilution.
The DGLAP system — and why the evolution basis exists
Each quark flavour couples to the gluon, and the gluon couples back to everything:
— thirteen coupled equations for $n_f = 6$ quarks, antiquarks, and the gluon. But the coupling matrix has huge redundancy: at LO the gluon talks to every $q_i$ identically. Diagonalize by flavour symmetry. Any non-singlet combination — differences in which the gluon contribution cancels, such as valence $V = \sum_i (q_i - \bar q_i)$ or flavour differences like $T_3 = (u + \bar u) - (d + \bar d)$ — evolves alone, by a scalar equation:
Only the singlet $\Sigma = \sum_i (q_i + \bar q_i)$ mixes with the gluon, in a single $2\times 2$ system:
This block-diagonalization is why the evolution basis exists. The 14 flavour combinations ($\Sigma, g, V, V_3, V_8, \ldots, T_3, T_8, \ldots$) in which our pipeline stores PDFs — the "14 flavours" axis of every FK table and every parametrization in our fits — are precisely the combinations in which DGLAP is one $2\times2$ block plus twelve independent scalars. Evolution code (EKO, §2.4) would be needlessly expensive and ill-conditioned in the physical $u, \bar u, d, \ldots$ basis; in the evolution basis it falls apart into trivial pieces.
Sum rules ride through evolution
The sum rules imposed on our fits at $Q_0$ would be worthless if evolution broke them at other scales. It cannot — the kernels' moments protect them. Valence number: since $\int_0^1 P_{\rm NS}(z)\, dz = \int_0^1 P_{qq}(z)\, dz = 0$ (worked example 1), every non-singlet number integral, e.g. $\int_0^1 (u - \bar u)\, dx = 2$, is exactly $\mu^2$-independent. Momentum: the second moments obey two conditions,
The first says a splitting quark's momentum is exactly shared between daughter quark and gluon (check: $\int z P_{qq} = -\tfrac43 C_F$ and $\int z P_{gq} = +\tfrac43 C_F$ — computed via the same technique as worked example 1); the second is the identical statement for a splitting gluon. Together they force $\frac{d}{d\ln\mu^2}\int_0^1 x\big[\Sigma(x) + g(x)\big]\, dx = 0$: total momentum, once normalized to 1 at $Q_0$, stays 1 forever. Evolution redistributes — quarks bleed momentum into glue and sea — but never creates or destroys. Our fits exploit this by imposing sum rules once, analytically, at $Q_0$ (§3.2), and trusting DGLAP with them everywhere else.
Moments: convolution becomes multiplication
Define the Mellin transform $\tilde f(N) = \int_0^1 dx\, x^{N-1} f(x)$. Its magic property: convolutions factorize,
(insert the definitions and swap integration order — the nested $x \le \xi$ region splits into independent $z$ and $\xi$ integrals). DGLAP in $N$-space is therefore an ordinary differential equation per moment, with the kernel moments — the anomalous dimensions $\gamma(N) \equiv \tilde P(N)$ — as coefficients. For the quark kernel (using $\int_0^1 \frac{z^{a}-1}{1-z} dz = -S_1(a)$, $S_1$ the harmonic number):
Checks: $\gamma_{qq}(1) = C_F[\tfrac32 + \tfrac12 - 2] = 0$ (number conservation again); $\gamma_{qq}(2) = C_F[\tfrac32 + \tfrac16 - 3] = -\tfrac{16}{9}$ (momentum drain — the number §2.4 evolves); and $\gamma_{qq}(N) \sim -2C_F \ln N \to -\infty$ at large $N$, the moment-space image of soft suppression at $z \to 1$. Large-$x$ PDFs get steeper with $Q^2$; small-$x$ ones get taller. This machinery is not just formal: §2.4 solves DGLAP exactly in $N$-space, and EKO's numerics run on precisely these anomalous dimensions.
No time passes and nothing happens to the proton: a 1.65-GeV proton and a 100-GeV-probed proton are the same quantum state. What changes is the question — how finely the probe bins the proton's constituents. DGLAP transports the answer between resolutions, which is why it is first-order in $\ln Q^2$ (a flow, like a renormalization group, which is exactly what it is) and why "evolving down" is conceptually suspect: coarse answers do not determine fine ones (§2.4 on backward instability). The safest mental model: PDFs at different $Q^2$ are different coordinates for the same object, with DGLAP the change-of-coordinates.
Worked example 1 — the plus-prescription normalization, explicitly
Plus-distribution piece: by the definition (with $f(z) = 1 + z^2$, $f(1) = 2$):
$$ \int_0^1 \frac{1+z^2}{(1-z)_+}\, dz = \int_0^1 \frac{(1+z^2) - 2}{1-z}\, dz = \int_0^1 \frac{z^2 - 1}{1 - z}\, dz = -\int_0^1 (1+z)\, dz = -\left[z + \tfrac{z^2}{2}\right]_0^1 = -\frac32 $$
using $z^2 - 1 = -(1-z)(1+z)$ — the would-be logarithmic divergence cancels algebraically before any integration, which is the plus-prescription doing its job.
Delta piece: $\int_0^1 \tfrac32\,\delta(1-z)\, dz = +\tfrac32$.
Total: $C_F\left[-\tfrac32 + \tfrac32\right] = 0$ ✓. Reading: the $-\tfrac32$ is the rate at which quarks leave any given momentum fraction by radiating; the $+\tfrac32$ is the survival probability deficit at $z = 1$. They must cancel because a quark that radiates is still a quark: gluon emission moves quarks down in $x$, never off the books. This single zero is what holds $\int (u - \bar u)\,dx = 2$ fixed from $Q_0 = 1.65$ GeV to the highest HERA scales — the valence sum rules our parametrization imposes analytically survive evolution to every dataset's $Q^2$.
Worked example 2 — small-$x$ gluon growth, estimated honestly
Setup: at small $x$ the singlet system is gluon-dominated and the dominant kernel piece is $P_{gg} \approx 2C_A/z$, whose Mellin transform has a pole $\gamma_{gg}(N) \approx \frac{2C_A}{N-1}$. Each step of evolution then pairs one $\ln Q^2$ with one $\ln(1/x)$ — the double-leading-log (DLA) tower. Summing it (saddle point of the inverse Mellin integral) gives
$$ x g(x, Q^2) \;\sim\; \exp\left[ 2\sqrt{\,\xi(Q^2)\, \ln\tfrac1x\,} \right], \qquad \xi(Q^2) = \frac{C_A}{\pi \beta_0} \ln\!\left[\frac{\ln(Q^2/\Lambda^2)}{\ln(Q_0^2/\Lambda^2)}\right] $$
Numbers ($n_f = 4$: $C_A/\pi\beta_0 = 36/25 = 1.44$; $\Lambda_4 = 0.25$ GeV; start the growth at $Q_0^2 = 2$ GeV$^2$; $\ln(1/x) = 6.91$): $\;\xi(10) = 1.44 \ln(5.08/3.47) = 0.549$; $\xi(100) = 1.44\ln(7.38/3.47) = 1.088$. Exponents: $2\sqrt{0.549 \times 6.91} = 3.90$ and $2\sqrt{1.088 \times 6.91} = 5.48$. Predicted growth factor: $e^{5.48 - 3.90} = e^{1.59} \approx \mathbf{4.9}$.
Honesty ledger: full NNLO evolution of a realistic input (what EKO does) gives a factor $\approx 2.3$ for the same step — the DLA captures the mechanism but overshoots the size, because it ignores (i) momentum conservation (the $2\pi\beta_0\,\delta(1-z)$ and large-$z$ parts of $P_{gg}$ drain what the $1/z$ pole pumps), (ii) the quark back-reaction ($g \to q\bar q$ leakage into the sea), and (iii) the fact that a realistic $g(x, 10\,\text{GeV}^2)$ has already absorbed part of the growth the flat-input DLA re-counts. The qualitative claims survive quantitative scrutiny: growth is doubly logarithmic (faster than any power of $\ln Q^2$, slower than any power of $1/x$), it is gluon-driven, and it is why the gluon — unmeasured directly by $F_2$ at LO — is nonetheless pinned by scaling violations $\partial F_2/\partial \ln Q^2$ wherever data span a lever arm in $Q^2$.