Four-Loop Flavour-Singlet Splitting in QCD

• The paper presents the analytic reconstruction of four-loop QCD flavour-singlet splitting functions by computing Mellin moments up to N=22.
• It employs advanced reduction techniques and Diophantine reconstruction to derive precise all‑x approximations with validated uncertainties.
• The results enable percent-level precision in N³LO QCD evolution, improving parton distribution analyses and collider phenomenology.

Four-loop flavour-singlet splitting functions $P_{ik}^{(3)}(x)$ are universal quantities in perturbative QCD that govern the scale evolution of parton distributions in the singlet sector at next-to-next-to-next-to-leading order (N³LO) in the strong coupling. The singlet splitting-function matrix appears in the coupled DGLAP equations, controlling the evolution of the singlet quark and gluon densities. Recent computational advances have now determined exact analytic moments up to $N = 22$ for all singlet channels, and approximate all-$x$ expressions are available with validated uncertainties, rendering the four-loop singlet sector ready for high-precision collider phenomenology (Falcioni et al., 11 Dec 2025).

1. Formalism and Definitions

The singlet evolution equations relate the singlet quark ($q_s$) and gluon ($g$) densities to their scale derivatives through the splitting-function matrix,

$$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$

where $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$, $a_s = \alpha_s/(4\pi)$, and

$$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$

The four-loop (N³LO) splitting functions $P_{ik}^{(3)}(x)$ determine the $N = 22$0 evolution. Their Mellin moments define the anomalous dimensions,

$N = 22$1

with the sign convention ensuring consistency with renormalization-group equations (Falcioni et al., 11 Dec 2025).

2. Computation of Four-Loop Moments

The recent determination of singlet splitting functions exploits advanced reduction techniques for four-loop operator matrix elements (OMEs), especially the Forcer program within Form for massless propagator-type integrals. All diagrams contributing to OMEs for even $N = 22$2 are computed for general gauge groups, and the color decomposition is made explicit. The endpoint constraints—large-$N = 22$3 threshold resummation (plus-distributions and $N = 22$4 for off-diagonal channels) and small-$N = 22$5 (high-energy) $N = 22$6 behavior—are incorporated to fully constrain the approximations (Falcioni et al., 11 Dec 2025, Falcioni et al., 2024, Falcioni et al., 2023).

Diophantine reconstruction methods are used to infer closed analytic forms for all non-rational contributions (i.e., $N = 22$7-functions), such that for $N = 22$8 the analytic expressions for all four channels are given in explicit numerical form, e.g.,

$N = 22$9

All $x$0 and $x$1 contributions have been reconstructed for general $x$2; only the $x$3 part of $x$4 and the $x$5 part of $x$6 remain incomplete (Falcioni et al., 11 Dec 2025).

3. All-$x$7 Approximations and Endpoint Constraints

To obtain $x$8-space approximations valid for all $x$9, the set of Mellin moments is supplemented by theoretical endpoint behaviors. For large $q_s$0, threshold resummation fixes plus-distribution and logarithmic terms up to $q_s$1 for off-diagonal channels and up to $q_s$2 for diagonal ones (Vogt et al., 2010). For small $q_s$3, the leading small-$q_s$4 singularities $q_s$5 are determined by BFKL resummation up to $q_s$6 at N³LO (Bonvini et al., 2018), and the coefficients are fixed from resummation and large-$q_s$7 limits (Davies, 2017, Davies et al., 2016).

The functional ansatz is constructed as a linear combination of interpolating basis functions in $q_s$8, powers of $q_s$9 ($g$0), $g$1 ($g$2), and plus-distributions, with coefficients fitted to moments and constrained by endpoints. Two representative approximations (labelled A, B) are selected to bracket uncertainties, e.g., for $g$3 in the pure-singlet case,

$g$4

where $g$5 and $g$6, $g$7 (Falcioni et al., 11 Dec 2025, Falcioni et al., 2023, Falcioni et al., 2024).

4. Color Structure and Moments: Status of Analytic Results

The singlet splitting functions admit a full color decomposition in terms of $g$8, $g$9, $$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$0, and quartic group invariants ($$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$1, $$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$2), with Mellin moments expressed as rational functions times harmonic sums and $$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$3-values. In the large-$$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$4 and $$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$5 sectors, analytic all-$$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$6 results have been established by explicit computation and Diophantine algebraic reconstruction (Davies, 2017, Gehrmann et al., 2023, Ruijl et al., 2016).

The highest power ($$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$7) terms, corresponding to maximal fermionic contributions, are given analytically for all channels in both Mellin space and $$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$8-space, exposing the hierarchy of logarithmic and polylogarithmic terms and confirming all theoretical endpoint predictions (Davies et al., 2016). The $$\frac{d}{d\ln \mu^2} \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} = \mathbf{P}(x,\alpha_s) \otimes \begin{pmatrix} q_s(x,\mu^2) \ g(x,\mu^2) \end{pmatrix} ,$$9 contributions to the pure-singlet component are also analytically available with explicit harmonic polylogarithm (HPL) representation up to weight 6 (Gehrmann et al., 2023).

Quartic invariants are essential for the correct scaling of the cusp anomalous dimensions and are included in all lowest-moment analytic results (Ruijl et al., 2016, Moch et al., 2021).

5. Threshold and Small-$$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$0 Resummation Properties

The large-$$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$1 asymptotics for all four-loop splitting functions are dictated by threshold resummation, with off-diagonal channels suppressed by explicit powers of $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$2 and double logarithms: $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$3 The exact coefficients for the leading three double-logarithmic terms are predicted and verified by explicit calculations (Vogt et al., 2010, Grunberg, 2011).

At small-$$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$4, the leading powers $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$5 in $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$6 and $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$7 are identified through BFKL-inspired resummations, and their impact on evolution is determined. The two-loop and three-loop boundary terms for the resummation are incorporated, and the hierarchy of the ensuing tower of logarithms is preserved (Bonvini et al., 2018).

6. Validation and Phenomenological Applications

The A/B envelope approximations constructed from $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$8 Mellin moments and endpoint constraints have been numerically validated against exact four-loop computations (Falcioni et al., 11 Dec 2025), with deviations typically within a few units in the final digit—for all singlet channels and $$\mathbf{P}(x,\alpha_s)=\sum_{n=0}^\infty a_s^{n+1} \mathbf{P}^{(n)}(x)$$9.

These functions can be used directly in N³LO QCD analyses such as DGLAP evolution for parton distribution functions and predictions for hard processes (Higgs, Drell–Yan, etc.) at LHC energies. Residual theoretical uncertainties from small-$a_s = \alpha_s/(4\pi)$0 are conservatively estimated. The four-loop corrections to the scale derivatives of the singlet quark and gluon distributions are generally below $a_s = \alpha_s/(4\pi)$1 for $a_s = \alpha_s/(4\pi)$2 at $a_s = \alpha_s/(4\pi)$3 (Falcioni et al., 11 Dec 2025, Falcioni et al., 2024, Moch et al., 2023).

All extensions relevant for ultra-high-scale applications, including $a_s = \alpha_s/(4\pi)$4 results, are available. The update of $a_s = \alpha_s/(4\pi)$5 at lower $a_s = \alpha_s/(4\pi)$6 further refines the estimate for the $a_s = \alpha_s/(4\pi)$7 coefficient, closing the remaining sources of uncertainty (Falcioni et al., 11 Dec 2025).

7. Outlook and Remaining Challenges

The analytic reconstruction of the four-loop singlet splitting functions is now complete for all $a_s = \alpha_s/(4\pi)$8 and $a_s = \alpha_s/(4\pi)$9 terms, and only limited $$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$0 color structures remain to be fully determined. Future work will aim at providing closed all-$$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$1 expressions for the non-rational ($$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$2) contributions to $$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$3 and $$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$4 channels, as well as extending the full analytic $$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$5-space forms to include all operator mixing, especially for the gluonic sector (Falcioni, 2022, Ruijl et al., 2016).

The robust computational and validation framework now established ensures that N³LO QCD evolution and phenomenology in the singlet sector achieves percent-level theoretical control. Further small-$$\mathbf{P}^{(n)}(x) = \begin{pmatrix} P_{qq}^{(n)}(x) & P_{qg}^{(n)}(x) \ P_{gq}^{(n)}(x) & P_{gg}^{(n)}(x) \end{pmatrix} .$$6 refinements and algorithmic developments will be critical for ultra-high energy collider applications and precision parton fits.