DGLAP Evolution and Scaling Violation

DGLAP evolution predicts how renormalized parton distributions change with the factorization scale. It resums collinear logarithms and explains logarithmic violations of Bjorken scaling. At leading order, a normalized nonsinglet distribution can be evolved explicitly for a short scale interval in momentum fraction and exactly in Mellin moments. The calculation below checks conserved quark number, the sign of momentum redistribution, and cancellation of factorization-scale dependence against the hard coefficient. It uses leading-twist, unpolarized, spacelike evolution in a fixed-flavor perturbative interval; heavy thresholds and endpoint resummation require additional treatment.

Required background. Collinear factorization and operator-defined PDFs supplies the distributions, convolution, and scheme cancellation. Coefficient evolution supplies the compensating evolution of short-distance coefficients.

Helpful background. Mellin transforms and scaling asymptotics supplies the moment-space diagonalization used below.

Spacelike evolution in momentum fraction

Take massless active quarks in $SU(3)$ QCD, with $\alpha_s=g_s^2/(4\pi)$ , $C_F=4/3$ , $C_A=3$ , and $T_F=1/2$ . Set the renormalization scale equal to the factorization scale, denoting both by $\mu$ , and use the $\overline{\mathrm{MS}}$ scheme. Both endpoint scales lie well above $\Lambda_{\mathrm{QCD}}$ , and no heavy threshold is crossed in the worked evolution. The nonperturbative boundary PDFs are inputs; perturbation theory determines their scale dependence.

The lightlike operators underlying the operator-defined PDFs renormalize by a matrix of convolution kernels. Let $Z$ relate bare and renormalized forward matrix elements by $f^{\mathrm{bare}}=Z\otimes f$ . At fixed bare data, with $\ell=\ln\mu^2$ ,

0=\frac{\mathrm d Z}{\mathrm d\ell}\otimes f +Z\otimes\frac{\mathrm df}{\mathrm d\ell}, \qquad P=-Z^{-1}\otimes\frac{\mathrm d Z}{\mathrm d\ell}.

Here $Z^{-1}$ is the perturbative convolution inverse and the identity kernel is $\delta_{ij}\delta(1-z)$ . Thus operator renormalization gives a state-independent kernel multiplying the hadron-dependent PDF vector:

\frac{\partial f_i(x,\mu)}{\partial\ln\mu^2} =\sum_j\left[P_{ij}(\alpha_s(\mu))\otimes f_j(\mu)\right](x),

[P\otimes f](x)=\int_x^1\frac{dz}{z}\,P(z)f\!\left(\frac{x}{z}\right), \qquad P_{ij}=\frac{\alpha_s}{2\pi}P_{ij}^{(0)}+ \left(\frac{\alpha_s}{2\pi}\right)^2P_{ij}^{(1)}+\cdots.

$P_{ij}$ describes parton $j$ feeding parton $i$ in this index convention. Reversing that convention transposes the singlet matrix, so the convention must travel with any numerical kernel. The full kernel $P$ includes its coupling; $P^{(0)}$ does not. The relation of PDF moments to renormalized twist-two operators and their compensating Wilson coefficients is developed in Schwartz 2014, §§32.4.3–32.4.5, pp. 693–695. The specific finite definition used in that DIS discussion differs from $\overline{\mathrm{MS}}$ ; the leading-order kernels used here agree, while finite coefficient terms and higher-order kernels must be kept in the chosen scheme.

For a quark nonsinglet combination such as $q_f-\bar q_f$ , gluon mixing cancels. The real collinear splitting $q\to qg$ , with the daughter quark carrying fraction $z$ , has weight $C_F(1+z^2)/(1-z)$ for $0<z<1$ and a transverse-momentum integral $\alpha_s\,\mathrm dk_T^2/(2\pi k_T^2)$ . This logarithmic measure explains why repeated splittings generate evolution in $\ln\mu^2$ . The spin/color sum and collinear limit giving this weight are shown in Schwartz 2014, §32.3, p. 682, Eqs. (32.55)–(32.58).

Writing $(1+z^2)/(1-z)=2/(1-z)-1-z$ isolates the soft endpoint. After the plus subtraction, the regular terms integrate to $-3/2$ . Quark-number conservation therefore fixes the remaining virtual endpoint coefficient to $+3/2$ :

P_{\mathrm{NS}}^{(0)}(z)=C_F \left[\frac{2}{(1-z)_+}-1-z+\frac32\delta(1-z)\right].

The plus distribution is defined by

\int_0^1dz\,\frac{\phi(z)}{(1-z)_+} =\int_0^1dz\,\frac{\phi(z)-\phi(1)}{1-z}.

For smooth test functions $\phi$ , it represents the real-emission singularity together with the endpoint subtraction needed for a finite distribution. The $\delta(1-z)$ term contains the remaining virtual endpoint contribution in this decomposition. Dropping either piece breaks conserved moments even if the pointwise kernel looks plausible. The distribution and complete leading-order kernel appear in Schwartz 2014, §32.2, pp. 679–681, Eqs. (32.39), (32.43), and (32.49). His derivative $\mu\,\mathrm d/\mathrm d\mu$ has coefficient $\alpha_s/\pi$ ; our $\mathrm d/\mathrm d\ln\mu^2$ has half that coefficient.

For a convolution whose lower limit is $x$ , the endpoint subtraction must still refer to the distribution defined on the whole unit interval. Applying that definition to a test function supported on $[x,1]$ gives

\begin{aligned} I_+[f](x) &\equiv\int_x^1\mathrm dz\,\frac{f(x/z)/z}{(1-z)_+}\\ &=\int_x^1\mathrm dz\, \frac{f(x/z)/z-f(x)}{1-z} +f(x)\ln(1-x). \end{aligned}

The logarithm is the contribution of the omitted interval $0<z<x$ . It must not be discarded when implementing the finite integral.

Moment-space derivation and exact checks

Define Mellin moments

f_N(\mu)=\int_0^1dx\,x^{N-1}f(x,\mu), \qquad \gamma_N=\int_0^1dz\,z^{N-1}P^{(0)}(z).

Changing variables $x=zy$ in the convolution gives

\int_0^1dx\,x^{N-1}[P^{(0)}\otimes f](x)=\gamma_Nf_N.

For convergent moments, interchange the triangular integration region $0<x<z<1$ and substitute $x=zy$ ; its Jacobian cancels the convolution’s $1/z$ , leaving separate $z$ and $y$ moments. For the plus term, perform its subtraction before this interchange. The integro-differential equation therefore becomes an ordinary one. For a fixed coupling and $t=\ln(\mu^2/\mu_0^2)$ ,

f_N(t)=\exp\!\left[\frac{\alpha_s}{2\pi}\gamma_Nt\right]f_N(0).

At one loop, $\beta_0=11-2n_f/3$ and

\frac{\mathrm d\alpha_s}{\mathrm d\ln\mu^2} =-\frac{\beta_0}{4\pi}\alpha_s^2, \qquad \frac{1}{f_N}\frac{\mathrm df_N}{\mathrm d\alpha_s} =-\frac{2\gamma_N}{\beta_0\alpha_s}.

The second equality applies to a nonzero scalar nonsinglet moment. Integrating it in a fixed- $n_f$ interval gives

f_N(\mu)=f_N(\mu_0) \left[\frac{\alpha_s(\mu)}{\alpha_s(\mu_0)}\right]^{-2\gamma_N/\beta_0}.

The $\gamma_N$ here is a moment of the coupling-stripped splitting kernel. It is not the operator anomalous dimension with the sign convention of the coefficient-evolution prerequisite: if $\mu\,\mathrm dO_N/\mathrm d\mu=-\gamma_N^{\mathrm{op}}O_N$ , then at this order $\gamma_N^{\mathrm{op}}=-(\alpha_s/\pi)\gamma_N$ .

For positive integer $N$ , the plus term contributes $-2H_{N-1}$ , where $H_N=\sum_{k=1}^N1/k$ and $H_0=0$ . The regular terms contribute $-1/N-1/(N+1)$ , so

\gamma_N=C_F\left[\frac32-2H_N+\frac{1}{N(N+1)}\right].

In particular,

\gamma_1=C_F\left[0-\frac32+\frac32\right]=0,

so the net nonsinglet quark number is scale independent. The next moment is an especially useful sign check:

\gamma_2=C_F\left[-2-\frac56+\frac32\right] =-\frac{16}{9}.

Thus a positive nonsinglet second moment decreases as the resolution scale grows: radiation redistributes longitudinal momentum toward smaller $x$ . For asymptotically free running, $\alpha_s(\mu)<\alpha_s(\mu_0)$ when $\mu>\mu_0$ , and the exponent $-2\gamma_2/\beta_0$ is positive. These are two independent ways to check the sign.

A normalized distribution evolved in momentum fraction

Choose the dimensionless toy boundary distribution

q_0(x)=6x(1-x),\quad 0<x<1, \qquad (q_0)_N=\frac{6}{(N+1)(N+2)}.

It has unit nonsinglet number and second moment $1/2$ . This is an analytic input for testing evolution, not a fit to a hadron. Define the integrated leading-order coupling by

\tau(\mu,\mu_0) =\int_{\ln\mu_0^2}^{\ln\mu^2} \frac{\alpha_s(\mu')}{2\pi}\,\mathrm d\ln\mu'^2 =\frac{2}{\beta_0} \ln\frac{\alpha_s(\mu_0)}{\alpha_s(\mu)}.

The final equality uses one-loop running. At fixed coupling, $\tau=\alpha_s t/(2\pi)$ . Iterating the evolution equation once gives $q(x,\mu)=q_0(x)+\tau r(x)+O(\tau^2)$ , where $r=P_{\mathrm{NS}}^{(0)}\otimes q_0$ . This is a short-step expansion of the leading-order evolution solution, not a calculation of the next-to-leading-order kernel.

The finite plus integral and the regular part can both be done explicitly:

\begin{aligned} I_+[q_0](x) &=3-12x+9x^2-6x(1-x)\ln x\\ &\quad+6x(1-x)\ln(1-x),\\ J(x) &\equiv\int_x^1\frac{\mathrm dz}{z}(1+z)q_0(x/z)\\ &=3-12x+9x^2-6x\ln x. \end{aligned}

For example, the subtracted integrand in $I_+$ reduces to $6x[(1-x)/z+(1-x)/z^2-x/z^3]$ , so no singular numerical cancellation is needed to obtain these expressions. Combining $2I_+-J+(3/2)q_0$ yields

\begin{aligned} r(x)=C_F\bigl[& 3(1-x)+6x(2x-1)\ln x\\ &+12x(1-x)\ln(1-x)\bigr]. \end{aligned}

The logarithms are locally finite for $0<x<1$ and have integrable endpoint limits in the moments used here. Direct integration gives

\int_0^1 r(x)\,\mathrm dx=0, \qquad \int_0^1 x\,r(x)\,\mathrm dx=-\frac{2C_F}{3}=-\frac89.

These agree with the independent Mellin calculation $\gamma_1(q_0)_1=0$ and $\gamma_2(q_0)_2=(-16/9)(1/2)$ . For instance, $\int_0^1x^a\ln x\,\mathrm dx=-1/(a+1)^2$ and $\int_0^1x^a\ln(1-x)\,\mathrm dx=-H_{a+1}/(a+1)$ for nonnegative integer $a$ suffice to check both integrals.

Within the leading-order evolution equation the corresponding moments are known to all orders in $\tau$ :

q_N(\mu)=\frac{6\,e^{\gamma_N\tau}}{(N+1)(N+2)}, \qquad q_1=1,\qquad q_2=\frac12e^{-16\tau/9}.

Expanding this result gives $q_2=1/2-(8/9)\tau+O(\tau^2)$ , exactly as in momentum fraction. An implementation can compare both moments without introducing a numerical Mellin inversion.

The short-step expansion is not a uniform relative approximation at the endpoints. In particular,

\frac{r(x)}{q_0(x)} =C_F\left[2\ln(1-x)-\frac12+O(1-x)\right], \qquad x\to1.

For any fixed positive $\tau$ , the truncated expression eventually becomes negative sufficiently close to $1$ ; this does not demonstrate a negative exact leading-order evolved distribution. At small $x$ , $r(x)\to3C_F$ while $q_0(x)\sim6x$ , so a small absolute correction need not be a small relative one there either. Use the expansion at fixed interior $x$ with $|\tau r(x)/q_0(x)|\ll1$ , or use its controlled finite-moment expansion. Large evolution intervals require solving or resumming the equation.

Cancellation against the coefficient function

Changing an auxiliary factorization scale at a fixed hard scale $Q$ must not change the factorized observable. For a scalar nonsinglet channel write its leading-power contribution as $F=C\otimes q$ , with electroweak weights absorbed into $C$ . Differentiating and using associativity of the same distributional convolution gives

\frac{\mathrm dC}{\mathrm d\ln\mu^2}=-C\otimes P, \qquad \frac{\mathrm d(C\otimes q)}{\mathrm d\ln\mu^2} =-C\otimes P\otimes q+C\otimes P\otimes q=0.

For the singlet system $C$ is a row vector and $q$ is replaced by the PDF column vector; that matrix order matters. The cancellation is exact within a consistently solved evolution system and holds only through the retained perturbative order in a truncated prediction.

To expose the cancellation with the toy input, choose a unit Born coefficient at the reference scale and retain only its scale-dependent evolution:

\begin{aligned} C(x,\mu)&=\delta(1-x)-\tau P_{\mathrm{NS}}^{(0)}(x) +O(\tau^2),\\ C(\mu)\otimes q(\mu) &=q_0+\tau r-\tau P_{\mathrm{NS}}^{(0)}\otimes q_0 +O(\tau^2)\\ &=q_0+O(\tau^2). \end{aligned}

In moments the same check is $C_N=e^{-\gamma_N\tau}$ and $C_Nq_N=(q_0)_N$ . Process-specific finite matching terms are not computed by this example; they must be supplied for an actual structure function in the same scheme. Nor does cancellation of $\mu$ imply Bjorken scaling: varying the physical hard scale $Q$ , commonly evaluating the factors at $\mu\sim Q$ , changes the observable. Evolution organizes that physical scaling violation while canceling dependence on the arbitrary separation scale. See the operator/coefficient argument in Schwartz 2014, §§32.4.3–32.4.5, pp. 693–695.

Singlet–gluon mixing and momentum conservation

Let $\Sigma=\sum_f(q_f+\bar q_f)$ . Its leading-order evolution mixes with the gluon:

\frac{d}{d\ln\mu^2} \begin{pmatrix}\Sigma\\ g\end{pmatrix} =\frac{\alpha_s}{2\pi} \begin{pmatrix}P_{qq}^{(0)}&P_{qg}^{(0)}\\P_{gq}^{(0)}&P_{gg}^{(0)}\end{pmatrix} \otimes \begin{pmatrix}\Sigma\\ g\end{pmatrix}.

In this singlet basis $P_{qg}^{(0)}$ includes the $2n_f$ quark/antiquark multiplicity. The off-diagonal kernels and their second moments are

\begin{aligned} P_{qg}^{(0)}(z)&=2n_fT_F[z^2+(1-z)^2], &\int_0^1\mathrm dz\,zP_{qg}^{(0)}(z)&=\frac{n_f}{3},\\ P_{gq}^{(0)}(z)&=C_F\frac{1+(1-z)^2}{z}, &\int_0^1\mathrm dz\,zP_{gq}^{(0)}(z)&=\frac{4C_F}{3}. \end{aligned}

The single-flavor kernels are given in Schwartz 2014, §32.2, pp. 681–682, Eqs. (32.50)–(32.54). At moment $N=2$ , momentum conservation requires the column sums of the moment matrix to vanish, fixing its diagonal entries to the negatives of these off-diagonal moments. For $n_f=3$ in the leading-order normalization used here, this check is

\Gamma_2= \begin{pmatrix} -16/9&1\\[2pt] 16/9&-1 \end{pmatrix}, \qquad (1,1)\Gamma_2=(0,0).

Therefore $\Sigma_2+g_2$ remains constant under this evolution. This condition is stronger than checking either component independently and catches an accidental transpose immediately.

The physical pattern is a redistribution, not creation of hadron momentum. Higher resolution typically resolves more gluons and sea quarks at small $x$ , while valence-like large- $x$ moments tend to decrease. Conservation of the total second moment follows from the conserved energy-momentum tensor, as explained in Schwartz 2014, §32.4.4, p. 694; it does not require each parton contribution to be conserved.

Thresholds, accuracy, and solution methods

A complete evolution record states:

Item	Required statement
boundary data	distributions and covariance at $\mu_0$ in a named scheme
kernels	spacelike order and index convention
coupling	running order, reference input, and active flavors
thresholds	matching scales, PDF matching matrices, and heavy-mass scheme
solver	$x$ -space grid, Mellin inversion, or another method with tolerances
validation	sum-rule residuals, convergence, and round-trip error

When the scale crosses a heavy flavor, coupling evolution, PDFs, and coefficient functions must all be matched. A variable-flavor description is not obtained by merely adding a new component to the PDF vector.

Logarithmic accuracy also must be consistent: an evolution kernel at one order, boundary distributions extracted with another, and coefficient functions at a third do not automatically combine into the highest of those labels.

Failure modes and limits

Space- versus timelike evolution. Fragmentation functions obey timelike equations. Leading kernels are related by a transpose under a suitable convention, but the equality does not persist unchanged beyond leading order.

Small $x$ . When $\alpha_s\ln(1/x)$ is order one, fixed-order DGLAP kernels may need high-energy resummation and matching.

Large $x$ . Threshold logarithms and power corrections are enhanced as $x\to1$ ; evolution alone is not an endpoint prediction.

Numerical drift. Positivity-looking curves can still violate exact moments. Monitor valence and momentum sum rules at every scale and make grid refinement part of the validation record.

Exercises

A missing endpoint term. Suppose a numerical implementation subtracts the $z=1$ value inside $I_+[q_0]$ but omits $q_0(x)\ln(1-x)$ . Find the error in the first moment after one short step. Can refining the integration grid repair it?

Solution

The implemented correction is $r_{\mathrm{wrong}}=r-2C_Fq_0\ln(1-x)$ . Since $\int_0^1q_0(x)\ln(1-x)\,\mathrm dx=-5/6$ ,

\int_0^1q_{\mathrm{wrong}}(x,\mu)\,\mathrm dx =1+\frac{5C_F}{3}\tau+O(\tau^2) =1+\frac{20}{9}\tau+O(\tau^2).

The spurious number increase is a distributional error, not a quadrature error. Grid refinement converges to the wrong answer unless the endpoint term is restored.

Running and compensation. In a three-flavor interval take $\alpha_s(\mu_0)=0.20$ and $\alpha_s(\mu)=0.16$ as toy perturbative inputs, with $\mu>\mu_0$ . Find $\tau$ , the exact leading-order $q_2$ , and its short-step approximation. Determine the coefficient moment $C_2$ that keeps the factorized second moment $C_2q_2$ fixed.

Solution

With $\beta_0=9$ ,

\tau=\frac29\ln(1.25)\simeq0.0495875,\qquad q_2=\frac12(0.8)^{32/81}\simeq0.457809.

The first-step result is $q_2\simeq1/2-(8/9)\tau\simeq0.455922$ . The difference is the omitted higher powers of $\tau$ in this leading-order equation; it is not an estimate of all higher-order QCD effects. The evolved coefficient is $C_2=e^{16\tau/9}\simeq1.09216$ , and $C_2q_2=1/2$ . Evolving both factors with the same sign would instead double the scale dependence.

Using the evolved distributions

The evolution map passed onward is

U^{\mathrm S}_{ij}(\mu,\mu_0; P^{\mathrm S},\alpha_s,n_f,\mathcal S,\text{thresholds}),

together with boundary distributions, their covariance, and numerical tolerances. Use the superscript $\mathrm S$ to prevent accidental reuse as the timelike map for fragmentation functions. Use high-energy QCD when the small- $x$ logarithm changes the counting.

References

Schwartz, Matthew D. Quantum Field Theory and the Standard Model. Cambridge University Press, 2014. DOI.