###### Abstract

We present a detailed benchmarking of different treatments of the QCD
evolution of unpolarized PDFs at approximate ${\rm N}^{3}$LO (a${\rm
N}^{3}$LO) order in the QCD coupling. Namely, the implementations in the
public a${\rm N}^{3}$LO releases of the MSHT and NNPDF global PDF fitters, as
well as that of the theoretical FHMRUVV collaboration are compared. This
follows the same procedure as in previous benchmarking exercises at lower
orders, that is by considering the impact of this evolution on a set of simple
toy PDFs. Excellent agreement between the MSHT and NNPDF results is found at
NNLO, and at a${\rm N}^{3}$LO when the same (FHMRUVV) implementation of the
splitting functions is used. In addition, in the data region only small
differences between these is seen when their original approximations are used
for the splitting functions. The origin of these differences, and the larger
ones observed at lower $x$, is well understood in terms of the differences in
the approximate splitting functions. Good convergence is also observed between
the exact and truncated solution methods for the DGLAP evolution equations.
Overall, this provides confidence in the precision of the existing
implementation of PDF evolution at ${\rm N}^{3}$LO.

DESY-24-088

A Benchmarking of QCD evolution at Approximate N3LO

A. Cooper-Sarkara, T. Cridgeb, F. Giulic, L. A. Harland-Langd, F. Hekhorne,f,
J. Hustong, G. Magnih,i, S. Mochj, R.S. Thorned

a Department of Physics, University of Oxford, Oxford, OX1 3RH, UK
b Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
c CERN, CH-1211 Geneva, Switzerland
d Department of Physics and Astronomy, University College London, London, WC1E
6BT, UK
e University of Jyvaskyla, Department of Physics, P.O. Box 35, FI-40014
University of Jyvaskyla, Finland
f Helsinki Institute of Physics, P.O. Box 64, FI-00014 University of Helsinki,
Finland
g Department of Physics and Astronomy, Michigan State University, East
Lansing, MI 48824, USA
h Department of Physics and Astronomy, Vrije Universiteit, NL-1081 HV
Amsterdam
i Nikhef Theory Group, Science Park 105, 1098 XG Amsterdam, The Netherlands
j II. Institute for Theoretical Physics, Hamburg University, D-22761 Hamburg,
Germany

##  1 Introduction

The scale dependence of parton distribution functions (PDFs) is governed by
QCD evolution equations whose kernels, the splitting functions, are accessible
within QCD perturbation theory. For about 20 years next-to-next-to-leading
order (NNLO) has been the standard for precision predictions in perturbative
QCD, including studies of the proton structure and fits of PDFs. For the
solution of the integro-differential equations describing parton evolution in
QCD, different codes have been developed within the community, making it
necessary to check that there is consistency between the evolution codes used
by the various groups. To this end a set of benchmark tables for parton
evolution up to NNLO accuracy was produced in [1, 2] in order to provide a
means of checking the (numerical) accuracy of any evolution code.

Given the increasing experimental precision and the challenge to push the
accuracy of QCD perturbation theory by one quantum loop beyond the state-of-
the-art, it is now necessary to extend the benchmarking of QCD evolution codes
to the next order, i.e. next-to-NNLO (N3LO) in perturbative QCD. In this
document, we provide a new set of benchmark tables for the evolution of
unpolarized PDFs.111A short summary has been published as part of the
activities of the “Standard Model” working group for the “Physics at TeV
Colliders” workshop (Les Houches, France, 12–30 June, 2023) [3].  As the
required QCD splitting functions are not yet fully known at this order, we
resort to approximations which are valid within a restricted kinematic range.
The outcomes of this study are

  * •

benchmarking of available codes for parton evolution to N3LO.

  * •

discussion and comparisons of different approximations for the QCD splitting
functions at N3LO.

We provide a brief summary of the theoretical status of unpolarized parton
evolution in QCD in Sec. 2, including a list of available tools and codes. We
also discuss the different approximations for the N3LO splitting functions.
The results and tables with benchmark numbers are listed in Sec. 3. We finish
with conclusions and an outlook in Sec. 4.

##  2 Evolution equations and their solutions

###  2.1 The renormalization group equations

The scale evolution of the parton distributions $f_{p}(x,\mu_{\rm
f}^{2})\equiv p(x,\mu_{\rm f}^{2})$, where $p=q_{i},\,\bar{q}_{i}\,,g$ with
$i=1,\,\ldots,n_{f}$, is governed by the $2n_{f}\\!+\\!1$ coupled integro-
differential equations

$\frac{d\,f_{p}(x,\mu_{\rm f}^{2})}{d\ln\mu_{\rm
f}^{2}}\>=\>\sum_{l=0}^{m}\,a_{\rm s}^{l+1}(\mu_{\rm r}^{2})\>{\cal
P}_{pp^{\prime}}^{(l)}\bigg{(}x,\frac{\mu_{\rm f}^{2}}{\mu_{\rm
r}^{2}}\bigg{)}\otimes f_{p^{\prime}}(x,\mu_{\rm f}^{2})\,,$ (1)

where summation over $p^{\prime}$ is understood. The factorization and
renormalization scales are denoted $\mu_{\rm f}$ and $\mu_{\rm r}$, and
$\otimes$ is the standard Mellin convolution in the fractional-momentum
variable $x$. The scale dependence of the strong coupling reads

$\frac{d\,a_{\rm s}}{d\ln\mu_{r}^{2}}\>=\>\beta_{\rm N^{m}LO}(a_{\rm
s})\,=\,-\sum_{l=0}^{m}\,a_{\rm s}^{l+2}\,\beta_{l}\,,$ (2)

where we abbreviate $a_{\rm s}\equiv\alpha_{\rm s}/(4\pi)$. The
$\beta$-function coefficients $\beta_{l}$ are known to five loops in QCD [4,
5, 6, 7], which allows for the solution of Eq. 2 through N4LO for the $a_{\rm
s}$ evolution.

The general splitting functions ${\cal P}^{(l)}$ in Eq. 1 reduce to the
simpler expressions $P^{(l)}(x)$ at equal scales $\mu_{\rm r}=\mu_{\rm f}$. Up
to N3LO the corresponding relations read

$\displaystyle{\cal P}^{(0)}\bigg{(}x,\frac{\mu_{\rm f}^{2}}{\mu_{\rm
r}^{2}}\bigg{)}$ $\displaystyle\\!=\\!$ $\displaystyle P^{(0)}(x)\,,$
$\displaystyle{\cal P}^{(1)}\bigg{(}x,\frac{\mu_{\rm f}^{2}}{\mu_{\rm
r}^{2}}\bigg{)}$ $\displaystyle\\!=\\!$ $\displaystyle
P^{(1)}(x)-\beta_{0}P^{(0)}(x)\,\ln\frac{\mu_{\rm f}^{2}}{\mu_{\rm r}^{2}}\,,$
$\displaystyle{\cal P}^{(2)}\bigg{(}x,\frac{\mu_{\rm f}^{2}}{\mu_{\rm
r}^{2}}\bigg{)}$ $\displaystyle\\!=\\!$ $\displaystyle
P^{(2)}(x)-\bigg{\\{}\beta_{1}P^{(0)}(x)+2\beta_{0}P^{(1)}(x)\bigg{\\}}\ln\frac{\mu_{\rm
f}^{2}}{\mu_{\rm r}^{2}}+\beta_{0}^{2}P^{(0)}(x)\,\ln^{2}\frac{\mu_{\rm
f}^{2}}{\mu_{\rm r}^{2}}\,,$ (3) $\displaystyle{\cal
P}^{(3)}\bigg{(}x,\frac{\mu_{\rm f}^{2}}{\mu_{\rm r}^{2}}\bigg{)}$
$\displaystyle\\!=\\!$ $\displaystyle
P^{(3)}(x)-\bigg{\\{}\beta_{2}P^{(0)}(x)+2\beta_{1}P^{(1)}(x)+3\beta_{0}P^{(2)}(x)\bigg{\\}}\ln\frac{\mu_{\rm
f}^{2}}{\mu_{\rm r}^{2}}$
$\displaystyle+\bigg{\\{}\frac{5}{2}\beta_{0}\beta_{1}P^{(0)}(x)+3\beta_{0}^{2}P^{(1)}(x)\bigg{\\}}\ln^{2}\frac{\mu_{\rm
f}^{2}}{\mu_{\rm r}^{2}}-\beta_{0}^{3}P^{(0)}(x)\ln^{3}\frac{\mu_{\rm
f}^{2}}{\mu_{\rm r}^{2}}\,.$

The splitting functions are currently known completely up to NNLO, i.e.
analytic expression for $P_{pp^{\prime}}^{(2)}(x)$ are available [8, 9]. At
N3LO, this is also the case for non-singlet splitting functions in the planar
limit [10] as well as for some parts of the functions
$P_{pp^{\prime}}^{(3)}(x)$ proportional to powers of $n_{f}$ [11, 12, 13, 14].
In addition, there are a number of Mellin moments available for all splitting
functions at N3LO [15, 16, 17, 18, 19, 20] (at least the first five moments,
often ten or even more for specific colour factors). This information, taken
together with knowledge on the functional form of the splitting functions in
the limits $x\to 1$ and $x\to 0$, to be discussed below, serves as the basis
for approximations to all N3LO splitting functions, allowing for a solution of
Eq. 1 to this order.

###  2.2 Heavy-quark treatment

The renormalization group equations, Eqs. 1 and 2, are viable for a fixed
number of flavors $n_{f}$ and so we need to specify the prescription for
changing $n_{f}$. The transition $n_{f}\to n_{f}+1$ is made at the scales of
the heavy-quark pole masses,

$m_{c}\>=\>\sqrt{2}\mbox{ GeV}^{2}\>,\quad m_{b}\>=\>4.5\mbox{ GeV}^{2}\>.$
(4)

The matching conditions for the strong coupling in the transition $n_{f}\to
n_{f}+1$ derive from the decoupling formulae, given by

$a_{\rm s}^{\,(n_{f}+1)}(k_{\rm r}m_{h}^{2})\>=\>a_{\rm s}^{\,(n_{f})}(k_{\rm
r}m_{h}^{2})+\sum_{n=1}^{m}\,\Big{(}a_{\rm s}^{\,(n_{f})}(k_{\rm
r}m_{h}^{2})\Big{)}^{n+1}\sum_{l=0}^{n}\,c_{n,l}\,\ln k_{\rm r}\>\>.$ (5)

The matching coefficients $c_{n,l}$ are known to four loops [21, 22], so that
five-loop running can be obtained consistently from Eq. 2, taking into account
threshold effects. The benchmark computations presented here use running
$a_{\rm s}^{\,(n_{f})}$ at N3LO with the four-loop $\beta$-function.

The matching conditions at these thresholds $\mu_{\rm f}^{2}=m_{h}^{2}$,
$h=c,\,b,\,t$ for the PDFs involve the operator matrix elements (OMEs) with
heavy-quarks, $A_{pp^{\prime}}(x)$, which are accessible in perturbation
theory in an expansion in $a_{\rm s}$ with coefficients
$A_{pp^{\prime}}^{(l)}$, known [23] to two loops, i.e.
$A_{pp^{\prime}}^{(2)}$. For the benchmark computations, we apply these
conditions up to NNLO(=Nm=2LO) and they read [24] ($\delta_{m2}=1$ at Nm=2LO
below)

$l_{i}^{\,(n_{f}+1)}(x,m_{h}^{2})\>=\>l_{i}^{\,(n_{f})}(x,m_{h}^{2})+\delta_{m2}\>a_{\rm
s}^{2}\>A^{NS,(2)}_{qq,h}(x)\otimes l_{i}^{\,(n_{f})}(x,m_{h}^{2})\,,$ (6)

where $l=q,\,\bar{q}$ and $i=1,\ldots n_{f}$, and

$\displaystyle g^{\,(n_{f}+1)}(x,m_{h}^{2})\>\>\>$ $\displaystyle\\!=\\!\\!$
$\displaystyle g^{\,(n_{f})}(x,m_{h}^{2})+$ (7)
$\displaystyle\delta_{m2}\,a_{\rm
s}^{2}\,\Big{[}A_{gq,h}^{S,(2)}(x)\otimes\Sigma^{\,(n_{f})}(x,m_{h}^{2})+A_{gg,h}^{S,(2)}(x)\otimes
g^{\,(n_{f})}(x,m_{h}^{2})\Big{]}\,,$
$\displaystyle(h+\bar{h})^{\,(n_{f}+1)}(x,m_{h}^{2})\\!$
$\displaystyle\\!=\\!\\!$ $\displaystyle\delta_{m2}\,a_{\rm
s}^{2}\>\Big{[}\tilde{A}_{hq}^{S,(2)}(x)\otimes\Sigma^{\,(n_{f})}(x,m_{h}^{2})+\tilde{A}_{hg}^{S,(2)}(x)\otimes
g^{\,(n_{f})}(x,m_{h}^{2})\Big{]}\,,\quad$

where $h=\bar{h}$ and
$\Sigma^{\,(n_{f})}\equiv\sum_{i=1}^{n_{f}}(q_{i}+\bar{q}_{i})$. The N3LO QCD
corrections to massive OMEs are also available [25, 26, 27]. Those for
$\tilde{A}_{hg}^{S,(3)}$ have recently been completed [28, 29], superseding
previous approximations [30, 31].

The massive OMEs also provide the transition matrix elements in the variable-
flavor number scheme (VFNS), where the transition $n_{f}\to n_{f}+1$ treats
one heavy-quark flavor at the time, i.e. the single-mass case [32]. At higher
orders, both the $c$ and $b$-quark loops appear in massive OMEs [33, 34, 35].
Accounting for those corrections requires the transition $n_{f}\to n_{f}+2$,
i.e. the two-mass case [36]. For the benchmark computations, the VFNS
transition employs the single-mass case at NNLO only.

###  2.3 Available information and tools

As it is not possible to solve Eq. 1 in closed form beyond LO, a numerical
algorithm has to be applied for which several public codes exist: APFEL [37],
apfel++ [38], EKO [39], HOPPET [40], PEGASUS [41], and QCDNUM [42]. The codes
can be split into two different categories according to how they approach the
solution of Eq. 1, by either taking a direct approach in momentum fraction
space (APFEL, apfel++, HOPPET, QCDNUM) or solving in the conjugate Mellin-
space (EKO, PEGASUS).

In the following we give a short review on the different approximation
strategies applied by the MSHT group [43, 44], the NNPDF collaboration [45,
46], and FHMRUVV [10, 17, 18, 19, 20] \- for a detailed discussion we refer to
the respective references.

####  2.3.1 MSHT

The first attempt to produce an approximate set of N3LO PDFs, incorporating
the known information about the N3LO, i.e. ${\cal O}(\alpha_{\rm s}^{4})$,
splitting functions was made by the MSHT group [43]. At the time this study
took place less information was available than is currently the case. The
approximation for the splitting functions was based on the following procedure
(precise details can be found in [43]).

  * •

Information on the Mellin moments of the splitting functions [10, 16] was
used.

  * •

The full information on the non-singlet splitting function presented in [10]
was used. This includes information at high $x$ from [47] as well as
information appropriate in the limit of large flavour and colour number.

  * •

The available information on leading $\ln(1/x)$ resummations from [48, 49, 50,
51, 52, 53, 54, 55] was used.

For each splitting function the terms which are known exactly were input,
together with a set of basis functions chosen to be concentrated over a
variety of $x$ regions. For each splitting function one more function than
there are known Mellin moments was chosen (i.e. a set of 5 basis functions to
accompany the 4 moments known at the time for singlet splitting functions and
a set of 9 basis functions to accompany the 8 moments for the nonsinglet
splitting function). This leaves an unknown parameter, $\rho$ which reflects
the splitting function uncertainty and was chosen to be the coefficient of the
most divergent unknown piece of our set of basis functions at small $x$. The
variation in this parameter was then constrained by two requirements: for
allowed $\rho$ at high $x$, when the other coefficients of basis parameters
are determined by matching the full splitting functions to the known Mellin
moments, then the splitting function corrections at N3LO should not be large
compare to those at lower orders; at sufficiently small $x$, for a fixed value
of $\rho$, the splitting function is required to be contained within the range
of variation predicted from exploring a full range of choice of basis
functions. Once a particular choice of basis functions has been made there is
then a preferred value of $\rho$ and a variation which incurs a $\chi^{2}$
penalty, i.e. there is a prior based on all the available knowledge about the
function and lower orders. This prior was then modified in a fit to data and a
posterior “best-fit” value and uncertainty was obtained for the value of
$\rho$ for each splitting function. The constraints on the splitting functions
were applied in moment space in order to determine the values of $\rho$, but
the splitting function was then expressed in $x$ space and evolution performed
via numerical solution of the renormalization group equations in Eq. 1.

This approach was used in both the first release of the approximate N3LO PDFs
[43] and subsequent work including also QED evolution [44], examining various
aspects of the data used for the high $x$ gluon in the aN3LO PDFs [56], and
the first determination of the strong coupling at aN3LO in a global PDF fit
[57].

####  2.3.2 NNPDF

The default program adopted by the NNPDF collaboration for DGLAP evolution is
the recently published code EKO [39]. It focuses on solving Eq. 1 in terms of
evolution kernel operators (EKOs) $E$ which are independent from the actual
PDF $f_{p}(x,\mu_{\rm f}^{2})$, depending only on the theory setup (such as
perturbative orders or the initial scale $\mu_{\rm f,0}^{2}$):

$f_{p}(\mu_{\rm f}^{2})=E_{pp^{\prime}}(\mu_{\rm f}^{2}\leftarrow\mu_{\rm
f,0}^{2})\otimes f_{p^{\prime}}(\mu_{\rm f,0}^{2})$ (8)

The EKOs are computed in Mellin space and transformed to momentum fraction
space via interpolation techniques. EKO provides several strategies for
solving the respective renormalization group equations which are
perturbatively equivalent, but differ by the number of resummed terms. As EKOs
are independent of the boundary PDF values, they are ideally suited for a PDF
fit and as such EKO is integrated into the pineline framework [58].

The EKO strategy for approximating the N3LO splitting function $P^{(3)}$, can
be summarized as follows:

  * •

The available information from soft gluon resummation in the threshold region
$x\to 1$ [10, 11, 47, 59, 60, 61] is used.

  * •

The available information from BFKL resummation in the high-energy region
$x\to 0$ [48, 49, 50, 51, 52, 53, 62, 63, 64] is used.

  * •

The available information in the limit of large number of flavors [10, 11, 12,
14] is used.

  * •

The available information on the low Mellin moments [10, 16, 17, 18, 19] is
used.

All limits are combined in a unique way and a number of sub-leading terms
corresponding to the number of known Mellin moments are added. Precise details
can be found in Section 2 of [46]. In particular the specific choice of the
form of the sub-leading terms is arbitrary and drawn from a pre-defined list.
This freedom represents the uncertainty associated with the partial knowledge
on splitting functions and is referred to as Incomplete Higher Order
Uncertainty (IHOU). IHOU are estimated independently from other possible
sources of theory errors, such as the incomplete perturbative expansion
(referred to as missing higher order uncertainty (MHOU) [65, 66]).

####  2.3.3 FHMRUVV

Approximate results for the N3LO splitting functions have been constructed
based on the available Mellin moments [15, 16, 17, 18, 19, 20] and all known
results for the large-$x$ and small-$x$ limits. The large-$x$ behavior of the
${\overline{\rm MS}}$-scheme diagonal splitting functions reads

$P_{pp,\,x\to
1\,}^{\,(l-1)}(x)\>\>=\>\;\frac{A_{p}^{(l)}}{(1-x)_{+}}\,+\,B_{p}^{\,(l)}\,\delta(1\\!-\\!x)\,+\,C_{p}^{\,(l)}\,\ln(1-x)\,-\,A_{p}^{(l)}+D_{p}^{\,(l)}\>\>,$
(9)

where the four-loop cusp anomalous dimensions $A_{p}^{(4)}$ [59, 67] as well
as the coefficents $C_{p}^{\,(4)}$ and $D_{p}^{\,(4)}$ [19, 47] are known
analytically, while the virtual anomalous dimensions $B_{p}^{\,(4)}$
multiplying $\delta(1-x)$ have been determined in approximate form [68, 69].
For all N3LO splitting functions further large-$x$ constraints are known [61,
70, 71] and are also used. Information on their small-$x$ behavior is also
available. In case of $P_{gg}^{(3)}$ it is governed by the BFKL logarithms, of
which the next-to-leading logarithmic (NLL) correction has been calculated
[52, 54] and transformed to the ${\overline{\rm MS}}$-scheme some time ago
[72, 73] (for $P_{gq}^{(3)}$, see [55]). Sub-dominant small-$x$ terms for all
splitting functions have also been derived [62].

The procedure is then based on a total of 80 approximations featuring slightly
different functional forms for the splitting functions
$P_{pp^{\prime}}^{\,(3)}(x)$, consistent with their respective known endpoint
behaviour. These approximations define a range and for the two boundaries of
this range selected representatives $P_{pp^{\prime}\,A}^{\,(3)}(x)$ and
$P_{pp^{\prime},\,B}^{\,(3)}(x)$ are presented to provide the error bands for
$n_{f}=3,\,4,\,5$ light flavours. The quality of this procedure is tested
through approximations of the known NNLO splitting functions, demonstrating
good convergence and providing information on the range of $x$ for which
residual uncertainties due to missing higher Mellin moments are small.

The approximate N3LO splitting functions are valid at large-$x$, with the
uncertainties increasing towards smaller values of $x$. The FHMRUVV results
[10, 17, 18, 19] for $P_{qq}^{\,(3)}(x)$ (non-singlet and pure-singlet) and
$P_{qg}^{\,(3)}(x)$, based on ten Mellin moments $N\leq 20$, are reliable
approximations down to $x\simeq 10^{-4}$. The quantities $P_{gq}^{\,(3)}(x)$
and $P_{gg}^{\,(3)}(x)$, based on five Mellin moments $N\leq 10$, are
reasonable approximations down to $x\simeq 10^{-3}$. Their convolution with
PDFs will dampen the uncertainty at small-$x$, so that one gains one
additional order of magnitude on the right hand side of the evolution in Eq. 1
due to fast falling PDFs as $x\to 1$. The N3LO approximations to splitting
functions convoluted with PDFs are, therefore, reliable in the range
$10^{-5}\lesssim x\leq 1$ ($P_{qq}^{\,(3)}(x)$, $P_{qg}^{\,(3)}(x)$) and
$10^{-4}\lesssim x\leq 1$ ($P_{gq}^{\,(3)}(x)$, $P_{gg}^{\,(3)}(x)$),
respectively. The latest results for $P_{gq}^{\,(3)}(x)$ [20], based on ten
Mellin moments ($N\leq 20$), have not been used in this benchmark comparison.

##  3 Benchmark tables

The corresponding lower order benchmark results are available at LO and NLO in
[1] and at NNLO in [2], where LO and NLO polarized evolution has also been
included. A number of minor typos have been reported in [39, 74] and are
listed here for completeness:

  * •

In [1], table headers, the combination $L_{+}$ has to be defined as
$L_{+}=2(\bar{d}+\bar{u})$.

  * •

In [1], header of Tab. 1 $\alpha_{s}(\mu_{\rm r}^{2}=10^{4}\mbox{
GeV}^{2})=0.117574$, as pointed out in [2].

  * •

In [1], Tab. 1, $xs_{+}(x=0.5,\mu_{\rm f}=10^{4}\mbox{ GeV}^{2})=7.3137\cdot
10^{-4}$ for $n_{f}=4$, and $xb_{+}(x=10^{-7},\mu_{\rm f}=10^{4}\mbox{
GeV}^{2})=4.6071\cdot 10^{+1}$ for $n_{f}=3\dots 5$, differ in the last digit.

  * •

In [1], Tab. 4, the value $xu_{v}(x=0.7,\mu_{\rm f}=10^{4}\mbox{
GeV}^{2})=2.0102\cdot 10^{-2}$, differs in the last digit.

  * •

In [1], Tab. 3, the values $xL_{-}(x=10^{-5},\mu_{\rm f}=10^{4}\mbox{
GeV}^{2})=1.0121\cdot 10^{-4}$ and $xL_{-}(x=10^{-1},\mu_{\rm f}=10^{4}\mbox{
GeV}^{2})=9.8435\cdot 10^{-3}$, contain wrong exponents.

  * •

In [2], Tab. 15, the values $xd_{v}(x=10^{-7},\mu_{\rm f}=10^{4}\mbox{
GeV}^{2})=1.0699\cdot 10^{-4}$ and $xg(x=10^{-7},\mu_{\rm f}=10^{4}\mbox{
GeV}^{2})=9.9694\cdot 10^{2}$, contain wrong exponents.

###  3.1 Initial conditions

The initial conditions for the reference results are taken from [1, 2] with
the evolution starting at

$\mu_{\rm f,0}^{2}\>=\>2\mbox{ GeV}^{2}\,,$ (10)

and the parametrizations of input distributions are chosen as

$\displaystyle xu_{v}(x,\mu_{\rm f,0}^{2})$
$\displaystyle=5.107200\>x^{0.8}\>(1-x)^{3}\,,$ $\displaystyle
xd_{v}(x,\mu_{\rm f,0}^{2})$ $\displaystyle=3.064320\>x^{0.8}\>(1-x)^{4}\,,$
$\displaystyle xg\,(x,\mu_{\rm f,0}^{2})$
$\displaystyle=1.700000\,x^{-0.1}(1-x)^{5}\,,$ (11) $\displaystyle
x\bar{d}\,(x,\mu_{\rm f,0}^{2})$
$\displaystyle=0.1939875\,x^{-0.1}(1-x)^{6}\,,$ $\displaystyle
x\bar{u}\,(x,\mu_{\rm f,0}^{2})$ $\displaystyle=(1-x)\>x\bar{d}\,(x,\mu_{\rm
f,0}^{2})\,,$ $\displaystyle xs\,(x,\mu_{\rm f,0}^{2})$
$\displaystyle=x\bar{s}\,(x,\mu_{\rm
f,0}^{2})\>=\>0.2\,x(\bar{u}+\bar{d}\,)(x,\mu_{\rm f,0}^{2})\,,$

where the valence distributions are defined as $q_{i,v}\equiv
q_{i}-\bar{q}_{i}$ and the value for the running coupling at the input scale
is fixed as

$\alpha_{\rm s}(\mu_{\rm r}^{2}\\!=\\!2\mbox{ GeV}^{2})\>=\>0.35\>\>.$ (12)

These initial conditions are used regardless of the perturbative order of the
evolution and of the ratio of the renormalization and factorization scales,
which take the values

$\mu_{\rm r}^{2}\>=\>k_{\rm r}\mu_{\rm f}^{2}\>\>,\quad\quad k_{\rm
r}\>=\>0.5\,,\>\>1\,,\>\>2\,,$ (13)

except for LO, where this ratio is fixed to unity.

We assume two different setups for the heavy quark treatment discussed in
Section 2.2:

  1. 1.

We assume a fixed number $n_{f}=4$ of quarks participating in DGLAP evolution
and we refer to this as the Fixed Flavor Number scheme (FFNS). We assume that
all distributions in Section 3.1 are given directly in this scheme and that
the charm distribution vanishes at the initial scale $\mu_{\rm f,0}^{2}$.
Similarly, we assume that the boundary condition of the strong coupling Eq. 12
is given in this scheme.

  2. 2.

We assume a dynamic range $n_{f}=3\ldots 5$ of participating quark flavors for
our Variable Flavor Number scheme (VFNS), where each quark (here charm and
bottom) gets activated when the evolution scale matches the quark mass. We
assume Sections 3.1 and 12 apply for $n_{f}=3$, such that an immediate
matching procedure is triggered at the beginning of the evolution.

###  3.2 Results

The benchmark values are reported in the same basis adopted in [2], and we
define:

$q_{v}=q-\bar{q}\,,\quad L_{-}=(\bar{d}-\bar{u})\,,\quad
L_{+}=2(\bar{d}+\bar{u})\,,\quad q_{+}=q+\bar{q}\,.$ (14)

In order to simplify the comparison, only NNLO matching conditions are
adopted, although the N3LO coefficients are now available in literature and
implemented in the evolution tools. The evolved PDFs are then evaluated at a
final scale $\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$, which is chosen as
representative for LHC processes. As a preliminary condition to the N3LO
benchmark, we have tested that the different evolution codes are able to
reproduce the NNLO benchmark tables of [2] with an accuracy of $0.01\%$ or
below. We adopt the same notation as in the previous study and write
$x^{a}=x\times 10^{a}$.

Table 1:  Results for the FFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the FHMRUVV splitting
functions approximation and the MSHT20aN3LO code. FHMRUVV in MSHT20aN3LO
$\,n_{f}=4$, $\,\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{v}$ | $xs_{+}$ | $xc_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 1.1073-4 | 5.3078-5 | 7.8440-6 | 1.4895+2 | -2.7409-5 | 7.3405+1 | 7.2681+1 | 1.0858+3
1.0-6 | 5.7851-4 | 3.1908-4 | 3.0969-5 | 7.1890+1 | -4.2596-5 | 3.5094+1 | 3.4520+1 | 5.0367+2
1.0-5 | 2.9147-3 | 1.6974-3 | 1.2301-4 | 3.3841+1 | 7.7558-6 | 1.6245+1 | 1.5791+1 | 2.2274+2
1.0-4 | 1.3657-2 | 8.0056-3 | 4.7181-4 | 1.5175+1 | 2.3197-4 | 7.0536+0 | 6.6950+0 | 9.0183+1
1.0-3 | 5.9557-2 | 3.4379-2 | 1.7184-3 | 6.2980+0 | 4.0044-4 | 2.7362+0 | 2.4588+0 | 3.1320+1
1.0-2 | 2.3118-1 | 1.2953-1 | 5.7449-3 | 2.2665+0 | -2.3112-4 | 8.5192-1 | 6.6316-1 | 8.1503+0
1.0-1 | 5.5109-1 | 2.7129-1 | 1.0031-2 | 3.9000-1 | -3.2796-4 | 1.1381-1 | 5.9754-2 | 9.0578-1
3.0-1 | 3.5031-1 | 1.3010-1 | 2.9715-3 | 3.5600-2 | -3.8277-5 | 9.0840-3 | 3.3480-3 | 8.4418-2
5.0-1 | 1.2107-1 | 3.1503-2 | 3.6016-4 | 2.4562-3 | -3.2055-6 | 5.8459-4 | 1.7823-4 | 8.1560-3
7.0-1 | 2.0068-2 | 3.0873-3 | 1.0076-5 | 6.3848-5 | -1.1067-7 | 1.3107-5 | 3.9940-6 | 3.9250-4
9.0-1 | 3.5111-4 | 1.7834-5 | -5.9504-8 | 1.7816-7 | -1.3937-10 | 6.9228-9 | 8.6110-9 | 1.2313-6
Table 2:  Results for the FFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the FHMRUVV splitting
functions approximation and the NNPDF code.  FHMRUVV in NNPDF $\,n_{f}=4$,
$\,\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{v}$ | $xs_{+}$ | $xc_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 9.8369-5 | 4.5170-5 | 7.5011-6 | 1.4885+2 | -2.9106-5 | 7.3353+1 | 7.2639+1 | 1.0846+3
1.0-6 | 5.6405-4 | 3.0895-4 | 3.0730-5 | 7.1894+1 | -4.6738-5 | 3.5096+1 | 3.4528+1 | 5.0358+2
1.0-5 | 2.8946-3 | 1.6810-3 | 1.2302-4 | 3.3858+1 | -3.5753-6 | 1.6253+1 | 1.5803+1 | 2.2282+2
1.0-4 | 1.3633-2 | 7.9832-3 | 4.7274-4 | 1.5182+1 | 2.1123-4 | 7.0571+0 | 6.7005+0 | 9.0222+1
1.0-3 | 5.9567-2 | 3.4382-2 | 1.7231-3 | 6.3007+0 | 3.9314-4 | 2.7376+0 | 2.4611+0 | 3.1336+1
1.0-2 | 2.3130-1 | 1.2962-1 | 5.7645-3 | 2.2668+0 | -1.9644-4 | 8.5223-1 | 6.6372-1 | 8.1534+0
1.0-1 | 5.5131-1 | 2.7140-1 | 1.0085-2 | 3.8981-1 | -3.1812-4 | 1.1389-1 | 5.9850-2 | 9.0623-1
3.0-1 | 3.5044-1 | 1.3015-1 | 3.0144-3 | 3.5450-2 | -3.8409-5 | 9.1019-3 | 3.3625-3 | 8.4444-2
5.0-1 | 1.2112-1 | 3.1518-2 | 3.7781-4 | 2.4085-3 | -3.3051-6 | 5.9065-4 | 1.8276-4 | 8.1586-3
7.0-1 | 2.0077-2 | 3.0900-3 | 1.3437-5 | 5.6818-5 | -1.1797-7 | 1.4269-5 | 4.9009-6 | 3.9235-4
9.0-1 | 3.4849-4 | 1.8135-5 | 1.3838-8 | 1.0982-7 | -2.1740-10 | 4.2627-8 | 3.4628-8 | 1.3275-6
$\mu_{\rm r}^{2}=2\ \mu_{\rm f}^{2}$
1.0-7 | 1.1400-4 | 6.3301-5 | 6.5066-6 | 1.4118+2 | -7.8056-6 | 6.9517+1 | 6.8801+1 | 1.0533+3
1.0-6 | 5.9327-4 | 3.4617-4 | 2.7505-5 | 6.9885+1 | 6.1035-7 | 3.4090+1 | 3.3522+1 | 4.9859+2
1.0-5 | 2.9066-3 | 1.7186-3 | 1.1268-4 | 3.3325+1 | 6.6447-5 | 1.5986+1 | 1.5535+1 | 2.2228+2
1.0-4 | 1.3478-2 | 7.9084-3 | 4.4131-4 | 1.5014+1 | 2.3633-4 | 6.9723+0 | 6.6151+0 | 9.0028+1
1.0-3 | 5.8973-2 | 3.4010-2 | 1.6441-3 | 6.2494+0 | 2.7889-4 | 2.7109+0 | 2.4339+0 | 3.1216+1
1.0-2 | 2.3048-1 | 1.2918-1 | 5.6473-3 | 2.2597+0 | -2.8956-4 | 8.4756-1 | 6.5859-1 | 8.1296+0
1.0-1 | 5.5246-1 | 2.7213-1 | 1.0081-2 | 3.9105-1 | -2.8317-4 | 1.1407-1 | 5.9775-2 | 9.0774-1
3.0-1 | 3.5229-1 | 1.3092-1 | 3.0315-3 | 3.5676-2 | -3.2738-5 | 9.1494-3 | 3.3664-3 | 8.4779-2
5.0-1 | 1.2205-1 | 3.1785-2 | 3.8118-4 | 2.4251-3 | -2.8167-6 | 5.9294-4 | 1.8094-4 | 8.2004-3
7.0-1 | 2.0291-2 | 3.1256-3 | 1.3604-5 | 5.6613-5 | -1.0201-7 | 1.3991-5 | 4.5001-6 | 3.9492-4
9.0-1 | 3.5422-4 | 1.8440-5 | 1.4099-8 | 9.3900-8 | -1.9148-10 | 3.4475-8 | 2.6320-8 | 1.3407-6
$\mu_{\rm r}^{2}=0.5\ \mu_{\rm f}^{2}$
1.0-7 | 4.4192-5 | -1.1184-5 | 9.0493-6 | 1.6657+2 | -8.7990-5 | 8.2215+1 | 8.1501+1 | 1.1828+3
1.0-6 | 4.3474-4 | 1.7517-4 | 3.5176-5 | 7.5683+1 | -1.8556-4 | 3.6992+1 | 3.6425+1 | 5.2150+2
1.0-5 | 2.6807-3 | 1.4586-3 | 1.3589-4 | 3.4392+1 | -2.3604-4 | 1.6521+1 | 1.6071+1 | 2.2470+2
1.0-4 | 1.3512-2 | 7.8398-3 | 5.1075-4 | 1.5206+1 | 3.9416-5 | 7.0699+0 | 6.7139+0 | 9.0326+1
1.0-3 | 5.9946-2 | 3.4688-2 | 1.8231-3 | 6.2896+0 | 6.1057-4 | 2.7333+0 | 2.4572+0 | 3.1414+1
1.0-2 | 2.3195-1 | 1.3013-1 | 5.9287-3 | 2.2606+0 | 1.9683-4 | 8.5037-1 | 6.6210-1 | 8.1731+0
1.0-1 | 5.5088-1 | 2.7116-1 | 1.0132-2 | 3.8936-1 | -3.0333-4 | 1.1393-1 | 5.9982-2 | 9.0498-1
3.0-1 | 3.4987-1 | 1.2993-1 | 3.0140-3 | 3.5489-2 | -4.1682-5 | 9.1433-3 | 3.4151-3 | 8.4372-2
5.0-1 | 1.2088-1 | 3.1453-2 | 3.7730-4 | 2.4190-3 | -3.6076-6 | 5.9739-4 | 1.9043-4 | 8.1626-3
7.0-1 | 2.0030-2 | 3.0824-3 | 1.3408-5 | 5.8137-5 | -1.2487-7 | 1.4967-5 | 5.6241-6 | 3.9298-4
9.0-1 | 3.4739-4 | 1.8076-5 | 1.3770-8 | 1.3969-7 | -2.2060-10 | 5.7582-8 | 4.9611-8 | 1.3299-6
Table 3:  Results for the FFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the MSHT20aN3LO prior
splitting functions approximation. MSHT20aN3LO prior $\,n_{f}=4$, $\,\mu_{\rm
f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{v}$ | $xs_{+}$ | $xc_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 1.2189-4 | 6.7096-5 | 2.4727-6 | 5.2287+1 | -9.9201-6 | 2.5070+1 | 2.4349+1 | 4.0508+2
1.0-6 | 5.8593-4 | 3.3362-4 | 1.2769-5 | 3.4289+1 | -1.9932-5 | 1.6292+1 | 1.5721+1 | 2.6377+2
1.0-5 | 2.8965-3 | 1.6914-3 | 6.8967-5 | 2.1730+1 | 1.3958-5 | 1.0188+1 | 9.7359+0 | 1.5546+2
1.0-4 | 1.3619-2 | 7.9736-3 | 3.4675-4 | 1.2464+1 | 2.0204-4 | 5.6972+0 | 5.3396+0 | 7.7766+1
1.0-3 | 5.9587-2 | 3.4383-2 | 1.5513-3 | 6.1230+0 | 3.7152-4 | 2.6481+0 | 2.3710+0 | 3.0816+1
1.0-2 | 2.3129-1 | 1.2959-1 | 5.7459-3 | 2.3265+0 | -2.3233-4 | 8.8203-1 | 6.9326-1 | 8.4107+0
1.0-1 | 5.5130-1 | 2.7140-1 | 1.0091-2 | 3.8948-1 | -3.2809-4 | 1.1372-1 | 5.9658-2 | 9.0604-1
3.0-1 | 3.5044-1 | 1.3015-1 | 3.0182-3 | 3.5400-2 | -3.8291-5 | 9.0869-3 | 3.3463-3 | 8.4342-2
5.0-1 | 1.2111-1 | 3.1516-2 | 3.7949-4 | 2.3927-3 | -3.2067-6 | 5.8587-4 | 1.7784-4 | 8.1760-3
7.0-1 | 2.0075-2 | 3.0885-3 | 1.3847-5 | 5.3271-5 | -1.1071-7 | 1.2962-5 | 3.5378-6 | 3.8989-4
9.0-1 | 3.5125-4 | 1.7841-5 | 1.6416-8 | 8.6161-9 | -1.3944-10 | 6.3097-9 | 2.0583-9 | 1.1858-6
Table 4:  Results for the FFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the MSHT20aN3LO
posterior splitting functions approximation. MSHT20aN3LO posterior
$\,n_{f}=4$, $\,\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{v}$ | $xs_{+}$ | $xc_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 1.2189-4 | 6.7096-5 | 2.4727-6 | 8.4227+1 | -9.9201-6 | 4.1040+1 | 4.0322+1 | 7.1749+2
1.0-6 | 5.8593-4 | 3.3362-4 | 1.2769-5 | 4.4674+1 | -1.9932-5 | 2.1485+1 | 2.0915+1 | 3.6931+2
1.0-5 | 2.8965-3 | 1.6914-3 | 6.8967-5 | 2.3876+1 | 1.3958-5 | 1.1261+1 | 1.0810+1 | 1.8224+2
1.0-4 | 1.3619-2 | 7.9736-3 | 3.4675-4 | 1.2426+1 | 2.0204-4 | 5.6784+0 | 5.3212+0 | 8.1561+1
1.0-3 | 5.9587-2 | 3.4383-2 | 1.5513-3 | 5.9728+0 | 3.7152-4 | 2.5730+0 | 2.2960+0 | 3.0671+1
1.0-2 | 2.3129-1 | 1.2959-1 | 5.7459-3 | 2.3273+0 | -2.3233-4 | 8.8244-1 | 6.9366-1 | 8.3365+0
1.0-1 | 5.5130-1 | 2.7140-1 | 1.0091-2 | 3.9002-1 | -3.2809-4 | 1.1399-1 | 5.9926-2 | 9.0673-1
3.0-1 | 3.5044-1 | 1.3015-1 | 3.0182-3 | 3.5367-2 | -3.8291-5 | 9.0707-3 | 3.3300-3 | 8.4315-2
5.0-1 | 1.2111-1 | 3.1516-2 | 3.7949-4 | 2.3999-3 | -3.2067-6 | 5.8951-4 | 1.8148-4 | 8.1793-3
7.0-1 | 2.0075-2 | 3.0885-3 | 1.3847-5 | 5.2544-5 | -1.1071-7 | 1.2598-5 | 3.1752-6 | 3.8961-4
9.0-1 | 3.5125-4 | 1.7841-5 | 1.6416-8 | -1.9342-9 | -1.3944-10 | 1.0345-9 | -3.2071-9 | 1.1788-6
Table 5:  Results for the FFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the NNPDF splitting
functions approximation.  NNPDF $\,n_{f}=4$, $\,\mu_{\rm f}^{2}=10^{4}\mbox{
GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{v}$ | $xs_{+}$ | $xc_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 9.1711-5 | 4.1520-5 | 6.3455-6 | 1.3633+2 | -2.9134-5 | 6.7094+1 | 6.6379+1 | 9.9074+2
1.0-6 | 5.4825-4 | 3.0047-4 | 2.8366-5 | 6.8193+1 | -4.6755-5 | 3.3245+1 | 3.2678+1 | 4.7660+2
1.0-5 | 2.8653-3 | 1.6659-3 | 1.1996-4 | 3.3061+1 | -3.5228-6 | 1.5855+1 | 1.5404+1 | 2.1712+2
1.0-4 | 1.3605-2 | 7.9705-3 | 4.7399-4 | 1.5100+1 | 2.1129-4 | 7.0162+0 | 6.6596+0 | 8.9608+1
1.0-3 | 5.9601-2 | 3.4403-2 | 1.7396-3 | 6.3082+0 | 3.9309-4 | 2.7415+0 | 2.4650+0 | 3.1379+1
1.0-2 | 2.3135-1 | 1.2964-1 | 5.7805-3 | 2.2682+0 | -1.9645-4 | 8.5296-1 | 6.6444-1 | 8.1650+0
1.0-1 | 5.5130-1 | 2.7140-1 | 1.0084-2 | 3.8980-1 | -3.1812-4 | 1.1388-1 | 5.9847-2 | 9.0618-1
3.0-1 | 3.5044-1 | 1.3015-1 | 3.0144-3 | 3.5451-2 | -3.8409-5 | 9.1020-3 | 3.3626-3 | 8.4445-2
5.0-1 | 1.2112-1 | 3.1518-2 | 3.7782-4 | 2.4085-3 | -3.3051-6 | 5.9064-4 | 1.8275-4 | 8.1587-3
7.0-1 | 2.0077-2 | 3.0900-3 | 1.3438-5 | 5.6818-5 | -1.1797-7 | 1.4269-5 | 4.9012-6 | 3.9230-4
9.0-1 | 3.4849-4 | 1.8135-5 | 1.3877-8 | 1.0977-7 | -2.1739-10 | 4.2631-8 | 3.4632-8 | 1.3315-6
$\mu_{\rm r}^{2}=2\ \mu_{\rm f}^{2}$
1.0-7 | 1.1006-4 | 6.1143-5 | 5.8224-6 | 1.3396+2 | -7.8173-6 | 6.5905+1 | 6.5189+1 | 9.9727+2
1.0-6 | 5.8391-4 | 3.4115-4 | 2.6103-5 | 6.7757+1 | 6.1119-7 | 3.3026+1 | 3.2458+1 | 4.8250+2
1.0-5 | 2.8891-3 | 1.7096-3 | 1.1086-4 | 3.2868+1 | 6.6482-5 | 1.5757+1 | 1.5306+1 | 2.1888+2
1.0-4 | 1.3461-2 | 7.9008-3 | 4.4203-4 | 1.4967+1 | 2.3636-4 | 6.9487+0 | 6.5916+0 | 8.9657+1
1.0-3 | 5.8993-2 | 3.4022-2 | 1.6540-3 | 6.2537+0 | 2.7886-4 | 2.7131+0 | 2.4361+0 | 3.1241+1
1.0-2 | 2.3051-1 | 1.2919-1 | 5.6569-3 | 2.2605+0 | -2.8957-4 | 8.4798-1 | 6.5901-1 | 8.1366+0
1.0-1 | 5.5246-1 | 2.7213-1 | 1.0080-2 | 3.9105-1 | -2.8317-4 | 1.1407-1 | 5.9773-2 | 9.0772-1
3.0-1 | 3.5229-1 | 1.3092-1 | 3.0315-3 | 3.5676-2 | -3.2738-5 | 9.1495-3 | 3.3665-3 | 8.4780-2
5.0-1 | 1.2205-1 | 3.1785-2 | 3.8118-4 | 2.4251-3 | -2.8167-6 | 5.9293-4 | 1.8093-4 | 8.2005-3
7.0-1 | 2.0291-2 | 3.1256-3 | 1.3604-5 | 5.6612-5 | -1.0201-7 | 1.3991-5 | 4.5003-6 | 3.9489-4
9.0-1 | 3.5422-4 | 1.8440-5 | 1.4122-8 | 9.3870-8 | -1.9148-10 | 3.4478-8 | 2.6323-8 | 1.3434-6
$\mu_{\rm r}^{2}=0.5\ \mu_{\rm f}^{2}$
1.0-7 | 3.1361-5 | -1.8232-5 | 6.8212-6 | 1.4192+2 | -8.8075-5 | 6.9888+1 | 6.9175+1 | 1.0040+3
1.0-6 | 4.0428-4 | 1.5878-4 | 3.0618-5 | 6.8336+1 | -1.8570-4 | 3.3318+1 | 3.2752+1 | 4.6992+2
1.0-5 | 2.6243-3 | 1.4295-3 | 1.3000-4 | 3.2801+1 | -2.3602-4 | 1.5726+1 | 1.5276+1 | 2.1381+2
1.0-4 | 1.3458-2 | 7.8155-3 | 5.1314-4 | 1.5042+1 | 3.9611-5 | 6.9878+0 | 6.6319+0 | 8.9170+1
1.0-3 | 6.0011-2 | 3.4728-2 | 1.8548-3 | 6.3047+0 | 6.1053-4 | 2.7409+0 | 2.4648+0 | 3.1498+1
1.0-2 | 2.3206-1 | 1.3018-1 | 5.9593-3 | 2.2634+0 | 1.9681-4 | 8.5182-1 | 6.6355-1 | 8.1950+0
1.0-1 | 5.5087-1 | 2.7116-1 | 1.0130-2 | 3.8935-1 | -3.0333-4 | 1.1393-1 | 5.9975-2 | 9.0489-1
3.0-1 | 3.4987-1 | 1.2993-1 | 3.0140-3 | 3.5489-2 | -4.1682-5 | 9.1436-3 | 3.4154-3 | 8.4375-2
5.0-1 | 1.2088-1 | 3.1453-2 | 3.7731-4 | 2.4189-3 | -3.6076-6 | 5.9738-4 | 1.9041-4 | 8.1628-3
7.0-1 | 2.0030-2 | 3.0824-3 | 1.3409-5 | 5.8135-5 | -1.2487-7 | 1.4967-5 | 5.6246-6 | 3.9290-4
9.0-1 | 3.4739-4 | 1.8076-5 | 1.3844-8 | 1.3960-7 | -2.2058-10 | 5.7593-8 | 4.9622-8 | 1.3364-6
Table 6:  Results for the VFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the FHMRUVV splitting
functions approximation and the MSHT20aN3LO code.  FHMRUVV in MSHT20aN3LO,
$\,n_{f}=3\ldots 5\,$, $\,\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{+}$ | $xc_{+}$ | $xb_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 1.1639-4 | 5.5630-5 | 8.3396-6 | 1.5915+2 | 7.8505+1 | 7.6906+1 | 6.4245+1 | 1.0973+3
1.0-6 | 6.0389-4 | 3.3257-4 | 3.2692-5 | 7.6153+1 | 3.7227+1 | 3.6062+1 | 2.9625+1 | 5.0544+2
1.0-5 | 3.0183-3 | 1.7563-3 | 1.2872-4 | 3.5490+1 | 1.7071+1 | 1.6248+1 | 1.3056+1 | 2.2175+2
1.0-4 | 1.4024-2 | 8.2158-3 | 4.8879-4 | 1.5724+1 | 7.3290+0 | 6.7755+0 | 5.2766+0 | 8.8931+1
1.0-3 | 6.0636-2 | 3.4976-2 | 1.7596-3 | 6.4308+0 | 2.8038+0 | 2.4567+0 | 1.8204+0 | 3.0517+1
1.0-2 | 2.3325-1 | 1.3055-1 | 5.8026-3 | 2.2752+0 | 8.5774-1 | 6.6785-1 | 4.5172-1 | 7.8160+0
1.0-1 | 5.4876-1 | 2.6967-1 | 9.9483-3 | 3.8453-1 | 1.1219-1 | 6.4546-2 | 3.7236-2 | 8.5051-1
3.0-1 | 3.4472-1 | 1.2772-1 | 2.9024-3 | 3.4706-2 | 8.8422-3 | 4.0495-3 | 2.1095-3 | 7.8671-2
5.0-1 | 1.1798-1 | 3.0615-2 | 3.4653-4 | 2.3789-3 | 5.6438-4 | 2.4370-4 | 1.2144-4 | 7.6252-3
7.0-1 | 1.9335-2 | 2.9654-3 | 9.2473-6 | 6.2100-5 | 1.2548-5 | 5.9902-6 | 2.9754-6 | 3.7046-4
9.0-1 | 3.3148-4 | 1.6782-5 | -6.3382-8 | 1.8337-7 | 6.2892-9 | 8.4763-9 | 3.5330-9 | 1.1753-6
Table 7:  Results for the VFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the FHMRUVV splitting
functions approximation and the NNPDF code.  FHMRUVV in NNPDF,
$\,n_{f}=3\ldots 5\,$, $\,\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{+}$ | $xc_{+}$ | $xb_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 1.0347-4 | 4.7296-5 | 7.9869-6 | 1.5846+2 | 7.8162+1 | 7.6575+1 | 6.3898+1 | 1.0941+3
1.0-6 | 5.8850-4 | 3.2158-4 | 3.2458-5 | 7.5933+1 | 3.7117+1 | 3.5959+1 | 2.9511+1 | 5.0475+2
1.0-5 | 2.9963-3 | 1.7380-3 | 1.2878-4 | 3.5427+1 | 1.7039+1 | 1.6221+1 | 1.3021+1 | 2.2169+2
1.0-4 | 1.3998-2 | 8.1908-3 | 4.8992-4 | 1.5705+1 | 7.3201+0 | 6.7688+0 | 5.2666+0 | 8.8942+1
1.0-3 | 6.0651-2 | 3.4981-2 | 1.7651-3 | 6.4269+0 | 2.8020+0 | 2.4559+0 | 1.8183+0 | 3.0527+1
1.0-2 | 2.3339-1 | 1.3065-1 | 5.8248-3 | 2.2744+0 | 8.5749-1 | 6.6796-1 | 4.5145-1 | 7.8173+0
1.0-1 | 5.4900-1 | 2.6978-1 | 1.0007-2 | 3.8425-1 | 1.1225-1 | 6.4658-2 | 3.7274-2 | 8.5068-1
3.0-1 | 3.4485-1 | 1.2776-1 | 2.9491-3 | 3.4522-2 | 8.8517-3 | 4.0601-3 | 2.1150-3 | 7.8667-2
5.0-1 | 1.1803-1 | 3.0629-2 | 3.6563-4 | 2.3176-3 | 5.6611-4 | 2.4371-4 | 1.2203-4 | 7.6246-3
7.0-1 | 1.9344-2 | 2.9683-3 | 1.2839-5 | 5.2471-5 | 1.2749-5 | 5.9801-6 | 3.0847-6 | 3.7020-4
9.0-1 | 3.2955-4 | 1.7232-5 | 1.0350-8 | 5.5700-8 | 1.9893-8 | 1.4418-8 | 7.7074-9 | 1.2669-6
$\mu_{\rm r}^{2}=2\ \mu_{\rm f}^{2}$
1.0-7 | 1.1909-4 | 6.5970-5 | 6.8741-6 | 1.4822+2 | 7.3036+1 | 7.1713+1 | 6.0917+1 | 1.0484+3
1.0-6 | 6.1538-4 | 3.5855-4 | 2.8837-5 | 7.2880+1 | 3.5589+1 | 3.4609+1 | 2.8797+1 | 4.9380+2
1.0-5 | 2.9935-3 | 1.7685-3 | 1.1716-4 | 3.4488+1 | 1.6568+1 | 1.5861+1 | 1.2847+1 | 2.1891+2
1.0-4 | 1.3780-2 | 8.0806-3 | 4.5472-4 | 1.5398+1 | 7.1650+0 | 6.6722+0 | 5.2108+0 | 8.8043+1
1.0-3 | 5.9855-2 | 3.4495-2 | 1.6768-3 | 6.3396+0 | 2.7570+0 | 2.4315+0 | 1.7991+0 | 3.0254+1
1.0-2 | 2.3216-1 | 1.3001-1 | 5.6938-3 | 2.2639+0 | 8.5089-1 | 6.6127-1 | 4.4740-1 | 7.7845+0
1.0-1 | 5.5058-1 | 2.7082-1 | 1.0018-2 | 3.8658-1 | 1.1275-1 | 6.3045-2 | 3.7010-2 | 8.5589-1
3.0-1 | 3.4779-1 | 1.2900-1 | 2.9789-3 | 3.4935-2 | 8.9514-3 | 3.8484-3 | 2.0908-3 | 7.9589-2
5.0-1 | 1.1956-1 | 3.1068-2 | 3.7136-4 | 2.3531-3 | 5.7380-4 | 2.2346-4 | 1.1959-4 | 7.7331-3
7.0-1 | 1.9700-2 | 3.0274-3 | 1.3120-5 | 5.3288-5 | 1.2861-5 | 5.2550-6 | 2.9616-6 | 3.7634-4
9.0-1 | 3.3899-4 | 1.7735-5 | 1.0779-8 | 5.3917-8 | 1.8648-8 | 1.2598-8 | 6.8401-9 | 1.2945-6
$\mu_{\rm r}^{2}=0.5\ \mu_{\rm f}^{2}$
1.0-7 | 4.5549-5 | -1.3718-5 | 9.7980-6 | 1.8191+2 | 8.9889+1 | 8.7771+1 | 7.1481+1 | 1.2220+3
1.0-6 | 4.5429-4 | 1.7961-4 | 3.7731-5 | 8.1754+1 | 4.0029+1 | 3.8512+1 | 3.1146+1 | 5.3370+2
1.0-5 | 2.7921-3 | 1.5132-3 | 1.4423-4 | 3.6683+1 | 1.7669+1 | 1.6627+1 | 1.3311+1 | 2.2754+2
1.0-4 | 1.3969-2 | 8.0954-3 | 5.3559-4 | 1.5967+1 | 7.4522+0 | 6.7828+0 | 5.3273+0 | 9.0292+1
1.0-3 | 6.1375-2 | 3.5488-2 | 1.8845-3 | 6.4753+0 | 2.8278+0 | 2.4397+0 | 1.8381+0 | 3.0881+1
1.0-2 | 2.3477-1 | 1.3154-1 | 6.0188-3 | 2.2728+0 | 8.5849-1 | 6.6835-1 | 4.5573-1 | 7.8549+0
1.0-1 | 5.4786-1 | 2.6903-1 | 1.0032-2 | 3.8184-1 | 1.1169-1 | 6.7883-2 | 3.7592-2 | 8.4316-1
3.0-1 | 3.4245-1 | 1.2676-1 | 2.9274-3 | 3.4233-2 | 8.8003-3 | 4.5516-3 | 2.1528-3 | 7.7581-2
5.0-1 | 1.1677-1 | 3.0270-2 | 3.6111-4 | 2.2949-3 | 5.6319-4 | 2.8913-4 | 1.2584-4 | 7.5139-3
7.0-1 | 1.9055-2 | 2.9204-3 | 1.2615-5 | 5.1944-5 | 1.2734-5 | 7.4168-6 | 3.2604-6 | 3.6473-4
9.0-1 | 3.2199-4 | 1.6832-5 | 9.9709-9 | 5.9105-8 | 2.1920-8 | 1.7299-8 | 8.9014-9 | 1.2443-6
Table 8:  Results for the VFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the MSHT20aN3LO prior
splitting functions approximation.  MSHT20aN3LO prior, $\,n_{f}=3\ldots 5\,$,
$\,\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{+}$ | $xc_{+}$ | $xb_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 1.2800-4 | 7.0312-5 | 2.6923-6 | 5.4816+1 | 2.6337+1 | 2.4743+1 | 2.1448+1 | 3.9567+2
1.0-6 | 6.1135-4 | 3.4762-4 | 1.3705-5 | 3.5854+1 | 1.7076+1 | 1.5914+1 | 1.4042+1 | 2.5956+2
1.0-5 | 2.9987-3 | 1.7496-3 | 7.2860-5 | 2.2642+1 | 1.0646+1 | 9.8250+0 | 8.5356+0 | 1.5346+2
1.0-4 | 1.3984-2 | 8.1825-3 | 3.6107-4 | 1.2894+1 | 5.9135+0 | 5.3608+0 | 4.4296+0 | 7.6535+1
1.0-3 | 6.0671-2 | 3.4982-2 | 1.5927-3 | 6.2577+0 | 2.7167+0 | 2.3699+0 | 1.7961+0 | 3.0052+1
1.0-2 | 2.3338-1 | 1.3062-1 | 5.8105-3 | 2.3375+0 | 8.8903-1 | 6.9912-1 | 4.6967-1 | 8.0754+0
1.0-1 | 5.4900-1 | 2.6978-1 | 1.0014-2 | 3.8396-1 | 1.1210-1 | 6.4442-2 | 3.7165-2 | 8.5071-1
3.0-1 | 3.4487-1 | 1.2777-1 | 2.9532-3 | 3.4489-2 | 8.8453-3 | 4.0478-3 | 2.1095-3 | 7.8598-2
5.0-1 | 1.1803-1 | 3.0628-2 | 3.6736-4 | 2.3104-3 | 5.6565-4 | 2.4326-4 | 1.2147-4 | 7.6445-3
7.0-1 | 1.9343-2 | 2.9667-3 | 1.3259-5 | 5.0879-5 | 1.2404-5 | 5.5290-6 | 2.8363-6 | 3.6793-4
9.0-1 | 3.3162-4 | 1.6789-5 | 1.5757-8 | 6.7806-9 | 5.6859-9 | 1.8831-9 | 1.7371-9 | 1.1316-6
Table 9:  Results for the VFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the MSHT20aN3LO
posterior splitting functions approximation.  MSHT20aN3LO posterior,
$\,n_{f}=3\ldots 5\,$, $\,\mu_{\rm f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{+}$ | $xc_{+}$ | $xb_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 1.2800-4 | 7.0312-5 | 2.6923-6 | 8.9533+1 | 4.3695+1 | 4.2104+1 | 3.7572+1 | 7.1965+2
1.0-6 | 6.1135-4 | 3.4762-4 | 1.3705-5 | 4.6946+1 | 2.2622+1 | 2.1462+1 | 1.9346+1 | 3.6755+2
1.0-5 | 2.9987-3 | 1.7496-3 | 7.2860-5 | 2.4866+1 | 1.1758+1 | 1.0938+1 | 9.7679+0 | 1.8033+2
1.0-4 | 1.3984-2 | 8.1825-3 | 3.6107-4 | 1.2832+1 | 5.8825+0 | 5.3303+0 | 4.5420+0 | 8.0208+1
1.0-3 | 6.0671-2 | 3.4982-2 | 1.5927-3 | 6.0992+0 | 2.6374+0 | 2.2908+0 | 1.7687+0 | 2.9886+1
1.0-2 | 2.3338-1 | 1.3062-1 | 5.8105-3 | 2.3390+0 | 8.8977-1 | 6.9985-1 | 4.6773-1 | 8.0041+0
1.0-1 | 5.4900-1 | 2.6978-1 | 1.0014-2 | 3.8449-1 | 1.1236-1 | 6.4711-2 | 3.7240-2 | 8.5137-1
3.0-1 | 3.4487-1 | 1.2777-1 | 2.9532-3 | 3.4457-2 | 8.8291-3 | 4.0316-3 | 2.1054-3 | 7.8573-2
5.0-1 | 1.1803-1 | 3.0628-2 | 3.6736-4 | 2.3177-3 | 5.6933-4 | 2.4693-4 | 1.2236-4 | 7.6478-3
7.0-1 | 1.9343-2 | 2.9667-3 | 1.3259-5 | 5.0141-5 | 1.2035-5 | 5.1605-6 | 2.7464-6 | 3.6765-4
9.0-1 | 3.3162-4 | 1.6789-5 | 1.5757-8 | -3.6339-9 | 4.7866-10 | -3.3152-9 | 4.9166-10 | 1.1249-6
Table 10:  Results for the VFNS aN3LO evolution for the initial conditions and
the input parton distributions given in Sec. 3.1, with the NNPDF splitting
functions approximation.  NNPDF, $\,n_{f}=3\ldots 5\,$, $\,\mu_{\rm
f}^{2}=10^{4}\mbox{ GeV}^{2}$
---
$x$ | $xu_{v}$ | $xd_{v}$ | $xL_{-}$ | $xL_{+}$ | $xs_{+}$ | $xc_{+}$ | $xb_{+}$ | $xg$
$\mu_{\rm r}^{2}=\ \mu_{\rm f}^{2}$
1.0-7 | 9.6344-5 | 4.3393-5 | 6.7727-6 | 1.4525+2 | 7.1554+1 | 6.9967+1 | 5.8515+1 | 1.0001+3
1.0-6 | 5.7167-4 | 3.1255-4 | 3.0007-5 | 7.2074+1 | 3.5187+1 | 3.4029+1 | 2.7983+1 | 4.7809+2
1.0-5 | 2.9653-3 | 1.7220-3 | 1.2572-4 | 3.4609+1 | 1.6630+1 | 1.5812+1 | 1.2709+1 | 2.1616+2
1.0-4 | 1.3969-2 | 8.1776-3 | 4.9163-4 | 1.5624+1 | 7.2793+0 | 6.7280+0 | 5.2371+0 | 8.8369+1
1.0-3 | 6.0687-2 | 3.5003-2 | 1.7829-3 | 6.4348+0 | 2.8060+0 | 2.4600+0 | 1.8214+0 | 3.0571+1
1.0-2 | 2.3345-1 | 1.3068-1 | 5.8414-3 | 2.2757+0 | 8.5822-1 | 6.6869-1 | 4.5199-1 | 7.8282+0
1.0-1 | 5.4900-1 | 2.6978-1 | 1.0007-2 | 3.8424-1 | 1.1224-1 | 6.4654-2 | 3.7272-2 | 8.5064-1
3.0-1 | 3.4485-1 | 1.2776-1 | 2.9491-3 | 3.4522-2 | 8.8518-3 | 4.0602-3 | 2.1151-3 | 7.8668-2
5.0-1 | 1.1803-1 | 3.0629-2 | 3.6564-4 | 2.3176-3 | 5.6610-4 | 2.4370-4 | 1.2202-4 | 7.6247-3
7.0-1 | 1.9344-2 | 2.9683-3 | 1.2840-5 | 5.2470-5 | 1.2749-5 | 5.9804-6 | 3.0846-6 | 3.7016-4
9.0-1 | 3.2955-4 | 1.7232-5 | 1.0383-8 | 5.5649-8 | 1.9895-8 | 1.4420-8 | 7.7216-9 | 1.2702-6
$\mu_{\rm r}^{2}=2\ \mu_{\rm f}^{2}$
1.0-7 | 1.1488-4 | 6.3669-5 | 6.1597-6 | 1.4072+2 | 6.9288+1 | 6.7964+1 | 5.7807+1 | 9.9296+2
1.0-6 | 6.0542-4 | 3.5321-4 | 2.7394-5 | 7.0695+1 | 3.4496+1 | 3.3517+1 | 2.7917+1 | 4.7808+2
1.0-5 | 2.9750-3 | 1.7590-3 | 1.1536-4 | 3.4025+1 | 1.6337+1 | 1.5629+1 | 1.2667+1 | 2.1563+2
1.0-4 | 1.3763-2 | 8.0727-3 | 4.5575-4 | 1.5351+1 | 7.1417+0 | 6.6490+0 | 5.1938+0 | 8.7698+1
1.0-3 | 5.9877-2 | 3.4508-2 | 1.6874-3 | 6.3441+0 | 2.7592+0 | 2.4338+0 | 1.8009+0 | 3.0279+1
1.0-2 | 2.3219-1 | 1.3002-1 | 5.7037-3 | 2.2647+0 | 8.5130-1 | 6.6169-1 | 4.4771-1 | 7.7911+0
1.0-1 | 5.5058-1 | 2.7082-1 | 1.0017-2 | 3.8658-1 | 1.1275-1 | 6.3043-2 | 3.7008-2 | 8.5586-1
3.0-1 | 3.4779-1 | 1.2900-1 | 2.9789-3 | 3.4936-2 | 8.9514-3 | 3.8485-3 | 2.0909-3 | 7.9590-2
5.0-1 | 1.1956-1 | 3.1068-2 | 3.7137-4 | 2.3530-3 | 5.7380-4 | 2.2346-4 | 1.1958-4 | 7.7332-3
7.0-1 | 1.9700-2 | 3.0274-3 | 1.3120-5 | 5.3287-5 | 1.2861-5 | 5.2552-6 | 2.9616-6 | 3.7631-4
9.0-1 | 3.3899-4 | 1.7735-5 | 1.0800-8 | 5.3886-8 | 1.8649-8 | 1.2599-8 | 6.8475-9 | 1.2968-6
$\mu_{\rm r}^{2}=0.5\ \mu_{\rm f}^{2}$
1.0-7 | 3.1561-5 | -2.1401-5 | 7.4030-6 | 1.5487+2 | 7.6366+1 | 7.4248+1 | 6.0692+1 | 1.0369+3
1.0-6 | 4.2125-4 | 1.6183-4 | 3.2892-5 | 7.3820+1 | 3.6062+1 | 3.4545+1 | 2.8065+1 | 4.8122+2
1.0-5 | 2.7315-3 | 1.4819-3 | 1.3816-4 | 3.4999+1 | 1.6827+1 | 1.5785+1 | 1.2679+1 | 2.1672+2
1.0-4 | 1.3911-2 | 8.0700-3 | 5.3883-4 | 1.5800+1 | 7.3686+0 | 6.6992+0 | 5.2677+0 | 8.9197+1
1.0-3 | 6.1446-2 | 3.5531-2 | 1.9191-3 | 6.4917+0 | 2.8361+0 | 2.4480+0 | 1.8444+0 | 3.0971+1
1.0-2 | 2.3489-1 | 1.3159-1 | 6.0512-3 | 2.2756+0 | 8.5998-1 | 6.6984-1 | 4.5683-1 | 7.8758+0
1.0-1 | 5.4786-1 | 2.6903-1 | 1.0031-2 | 3.8183-1 | 1.1168-1 | 6.7876-2 | 3.7587-2 | 8.4307-1
3.0-1 | 3.4245-1 | 1.2676-1 | 2.9275-3 | 3.4234-2 | 8.8006-3 | 4.5518-3 | 2.1530-3 | 7.7584-2
5.0-1 | 1.1677-1 | 3.0270-2 | 3.6112-4 | 2.2948-3 | 5.6318-4 | 2.8912-4 | 1.2583-4 | 7.5141-3
7.0-1 | 1.9055-2 | 2.9204-3 | 1.2616-5 | 5.1942-5 | 1.2735-5 | 7.4174-6 | 3.2603-6 | 3.6466-4
9.0-1 | 3.2199-4 | 1.6832-5 | 1.0033-8 | 5.9017-8 | 2.1925-8 | 1.7304-8 | 8.9339-9 | 1.2496-6

In Tables 2, 1, 3, 4 and 5 we report the central values for the FFNS aN3LO
evolution for the 3 different splitting function approximations from FHMRUVV,
MSHT and NNPDF respectively. The first part of the table refers to the central
scale $\mu_{\rm f}=\mu_{\rm r}$ while the second and the third parts refer to
the scale varied results. For MSHT we report only the values at the central
scale $\mu_{\rm f}=\mu_{\rm r}$. Taking the same splitting function
approximations from FHMRUVV as input, Tables 1 and 2 show a comparision of the
two aN3LO evolution codes used by MSHT20 and NNPDF in the FFNS scheme at the
central scale $\mu_{\rm f}=\mu_{\rm r}$. In the range $10^{-5}\lesssim
x\lesssim 0.5$ there is generally very good agreement between the two
evolution codes, with slightly larger deviations observed in regions of $x$
where PDFs are becoming very small, i.e., in the limits $x\to 0$ (vanishing
valence PDFs) and $x\to 1$. This very good agreement among the evolution codes
provides the basis for quantifying the aN3LO FFNS evolution, now with the
splitting function approximations from MSHT20 as input (Table 3 with MSHT20
prior and Table 4 with MSHT20 posterior) using the MSHT20 code, and for the
NNPDF approximations (Table 5) with the NNPDF code, cf. Section 2.3 for
details. Due to the different methodology in perparing those approximations,
larger deviations are observed as a result of the aN3LO evolution, especially
for smaller $x$ values, e.g., for the gluon PDF at $x\lesssim 10^{-3}$, when
sufficient constraints on the splitting fuction $P_{gg}^{\,(3)}(x)$ from
Mellin moments are lacking. In essence, these deviations demonstrate the
current uncertainties inherent in the knowledge of the evolution kernels.
However, these deviations would reduce with the use of more up to date input
information for each of the groups determination of the splitting functions
(particularly for the older MSHT determination).

The corresponding values for the VFNS evolution are listed in Tables 7, 6, 8,
9 and 10. Tables 7 and 6 show again very good agreement in the same $x$ range
as before between the two aN3LO evolution codes (MSHT20 and NNPDF) with the
FHMRUVV splitting function approximations, now employing the VFNS scheme,
while Tables 8 and 9 (MSHT20 prior and posterior with the MSHT20 code) and
Table 5 (NNPDF approximations with NNPDF code) in the VFNS scheme display the
same similarities and differences as observed in the FFNS scheme.

For simplicity we only give the central values in Tables 2, 1, 3, 4, 5, 7, 6,
8, 9 and 10, while the complete set containing all evolved PDFs with all the
splitting functions variations are available under the following URL:

https://www.hep.ucl.ac.uk/pdf4lhc/aN3LObenchmarking.shtml

Figure 1:  Absolute difference of the aN3LO evolution with respect to the
NNLO, for the initial conditions and the input parton distributions given in
Sec. 3.1 with the FFNS settings. We display results for FHMRUVV (blue), MSHT
prior (orange), MSHT posterior (green) and NNPDF (red) approximations. The
displayed FHMRUVV result is averaged between the values of the MSHT20 and
NNPDF evolution code.  Figure 2: Same as Fig. 1, but now displaying the
relative difference with respect to the NNLO evolution. Figure 3: Same as Fig.
1, but now for VFNS evolution settings. Figure 4: Same as Fig. 2, but now for
VFNS evolution settings.

The benchmark numbers are illustrated in Figs. 1, 2, 3 and 4, where we show
the absolute and relative difference with respect to the NNLO evolution for
the different aN3LO evolutions, for each PDF flavor combination, reported in
the previous tables. We display the FFNS and VFNS settings respectively with
the central scale $\mu_{\rm r}=\mu_{\rm f}$. The figures include an estimate
of the uncertainties due to the approximated splitting functions, which is
obtained by varying a single splitting function at a time during the
evolution, and taking the uncertainty, for each point in $x$, as $1/2$ of the
spread of the final PDF. The displayed FHMRUVV result is averaged between the
values obtained with the MSHT20 and NNPDF evolution codes. Since MSHT and
NNPDF parametrizations do not include uncertainties for the non-singlet
combinations, errors for the combinations $d_{v},u_{v}$, $s_{V}$ and $L_{-}$
have been set to 0, though the difference from different approximations (which
is largely at small $x$) can be regarded as a measure of the uncertainty. This
is generally very small on an absolute scale due to the smallness of the
corresponding non-singlet flavor combinations. Figs. 1 and 2 for the evolution
in the FFNS scheme show that all the N3LO approximations provide consistent
results in the data region of $x\gtrsim 10^{-3}$, with at most few percent
differences relative to NNLO and often less for most PDF flavor combinations,
whilst absolute differences are smaller still. There are some remaining
differences in the small-$x$ region where the uncertainties of the individual
approximations also grow due to unknown small-$x$ logarithms and so the
theoretical uncertainties are larger. Thus in the kinematic range of the LHC,
the benchmark numbers for the aN3LO evolution demonstrate good perturbative
convergence and a significant reduction of the residual theoretical
uncertainty. Hence, they should be sufficient for most collider-physics
applications. Similar findings are observed in Figs. 3 and 4 for the evolution
in the VFNS scheme.

As a final study, in order to check the stability of the DGLAP kernels at
different QCD orders, we investigate the effect of different DGLAP solution
methods [41], called exact and truncated, which differ by the inclusion of
higher order terms. The former is utilised by MSHT in their $x$-space
implementation and the latter by NNPDF in their Mellin space implementation.
In the NNPDF code either method of solution is implemented and we utilise this
to compare the impact of this difference. From Fig. 5 we observe that, as
higher orders in the splitting functions are included the difference between
the different solution methods is smaller, indicating a good perturbative
convergence. This is visible both for non-singlet like distributions ($L^{-}$)
and singlet-like ($L^{+}$ and $g$).

Figure 5: Relative difference in percent between the exact and truncated
solutiong of the DGLAP evolution at different pertubative orders. We adopt the
same evolution and boundary condition as in the benchmark exercise and display
the result for the non-singlet quark combinations $L^{-}$, $L^{+}$ and the
gluon $g$.

##  4 Conclusions and Outlook

The study has shown excellent agreement between MSHT and NNPDF evolution when
exactly the same versions of the N3LO splitting functions are used. In
addition, it is observed that the differences for $x>10^{-3}$ are very small
when the distinct original estimates used in [43] and [46] are used. Moreover,
the origin of the small differences in this region, and larger ones for
smaller $x$ is well understood in terms of differences in the approximate
splitting functions. Hence, the study has been a success in this regard.
Furthermore, we can see in Fig. 5 that there is also good convergence in the
difference between the exact DGLAP evolution and the truncated solution method
as the order of the equation used increases, and specifically that at N3LO
this is also small fractions of a percent except at very small $x$ or in very
small non-singlet quark combinations. Thus, overall, we can have high
confidence in our precision in PDF evolution at N3LO. However, there are
obvious further avenues for exploration.

As mentioned in Section 2.2, the ${\cal O}(\alpha_{S}^{3})$ versions of the
OMEs for transition across heavy flavour transition points are now complete
[27, 29, 30, 35, 75, 76, 77, 78, 79, 80], though the final parts only became
available while our study was already well underway. Also, not all parts are
currently in an easily usable form. Hence we only included the transition
matrix elements up to ${\cal O}(\alpha_{S}^{2})$, i.e. NNLO in this study. In
the future we can consider the effects of improving upon this.

There are also effects of additional information [14, 17, 18, 19], such as
moments, computed since the original MSHT estimates of the splitting
functions. We have studied the effect of the benchmark evolution of MSHT
including instead the FHMRUVV determination of the splitting functions, but it
would be of at least as much interest to examine the effect of this update in
a full PDF fit. Indeed, MSHT have preliminary results [81] which show changes
mainly in the gluon distribution, but remaining within the uncertainty band of
the original MSHT determination. These bring the change in the MSHT aN3LO
gluon compared to that at NNLO somewhat closer to that observed by NNPDF,
which also observes a dip in the gluon around $x\sim 0.01$ but of reduced
magnitude. To be explicit the change in the dip at $x\sim 0.01$ is reduced by
about 1.5$\%$. Finally very recently additional moments for $P_{gq}$ were
determined in [20], these are therefore included in neither of the PDF groups
approximations and so it would be interesting to see the effect on the PDFs,
and corresponding improvements in agreement. As such, the level of variation
in N3LO evolution illustrated in Figs. 1, 2, 3 and 4 is very much an upper
estimate, which will be reduced when all up to date information which can be
used to estimate N3LO splitting functions now available is used by all groups.
Uncertainties from other sources, e.g. limitations in the precision of cross
sections to use together with N3LO evolution in PDF extractions, is a possible
study for the future. In addition, there are also methodological differences
in the implementation of the associated theoretical uncertainties.

However, given that the results in this report already show a good level of
convergence in the most important regions for N3LO PDF evolution
phenomenological studies, we can conclude that at least as far as an
understanding of this evolution is concerned, we are reaching a point where
relevant uncertainties are very small and are estimated reliably.

## Acknowledgements

T. C. acknowledges that this project has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (Grant agreement No. 101002090 COLORFREE). L. H.-L. and
R.S. T. thank STFC for support via grant awards ST/T000856/1 and ST/X000516/1.
F. H. is supported by the Academy of Finland project 358090 and is funded as a
part of the Center of Excellence in Quark Matter of the Academy of Finland,
project 346326. S. M. acknowledges the ERC Advanced Grant 101095857 Conformal-
EIC.

## References

  * [1] W. Giele et al., _The QCD / SM working group: Summary report_ , hep-ph/0204316.
  * [2] M. Dittmar et al., _Working Group I: Parton distributions: Summary report for the HERA LHC Workshop Proceedings_ , hep-ph/0511119.
  * [3] J. Andersen et al., _Les Houches 2023: Physics at TeV Colliders: Standard Model Working Group Report_ , in _Physics of the TeV Scale and Beyond the Standard Model: Intensifying the Quest for New Physics_ , 6, 2024, 2406.00708.
  * [4] P. A. Baikov, K. G. Chetyrkin and J. H. Kühn, _Five-Loop Running of the QCD coupling constant_ , _Phys. Rev. Lett._ 118 (2017) 082002 [1606.08659].
  * [5] F. Herzog, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt, _The five-loop beta function of Yang-Mills theory with fermions_ , _JHEP_ 02 (2017) 090 [1701.01404].
  * [6] K. G. Chetyrkin, G. Falcioni, F. Herzog and J. A. M. Vermaseren, _Five-loop renormalisation of QCD in covariant gauges_ , _JHEP_ 10 (2017) 179 [1709.08541].
  * [7] T. Luthe, A. Maier, P. Marquard and Y. Schroder, _The five-loop Beta function for a general gauge group and anomalous dimensions beyond Feynman gauge_ , _JHEP_ 10 (2017) 166 [1709.07718].
  * [8] S. Moch, J. A. M. Vermaseren and A. Vogt, _The Three loop splitting functions in QCD: The Nonsinglet case_ , _Nucl. Phys. B_ 688 (2004) 101 [hep-ph/0403192].
  * [9] A. Vogt, S. Moch and J. A. M. Vermaseren, _The Three-loop splitting functions in QCD: The Singlet case_ , _Nucl. Phys. B_ 691 (2004) 129 [hep-ph/0404111].
  * [10] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt, _Four-Loop Non-Singlet Splitting Functions in the Planar Limit and Beyond_ , _JHEP_ 10 (2017) 041 [1707.08315].
  * [11] J. Davies, A. Vogt, B. Ruijl, T. Ueda and J. A. M. Vermaseren, _Large-nf contributions to the four-loop splitting functions in QCD_ , _Nucl. Phys. B_ 915 (2017) 335 [1610.07477].
  * [12] T. Gehrmann, A. von Manteuffel, V. Sotnikov and T.-Z. Yang, _Complete ${N}_{f}^{2}$ contributions to four-loop pure-singlet splitting functions_, _JHEP_ 01 (2024) 029 [2308.07958].
  * [13] T. Gehrmann, A. von Manteuffel, V. Sotnikov and T.-Z. Yang, _The $N_{f}C_{F}^{3}$ contribution to the non-singlet splitting function at four-loop order_, _Phys. Lett. B_ 849 (2024) 138427 [2310.12240].
  * [14] G. Falcioni, F. Herzog, S. Moch, J. Vermaseren and A. Vogt, _The double fermionic contribution to the four-loop quark-to-gluon splitting function_ , _Phys. Lett. B_ 848 (2024) 138351 [2310.01245].
  * [15] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt, _On quartic colour factors in splitting functions and the gluon cusp anomalous dimension_ , _Phys. Lett. B_ 782 (2018) 627 [1805.09638].
  * [16] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren and A. Vogt, _Low moments of the four-loop splitting functions in QCD_ , _Phys. Lett. B_ 825 (2022) 136853 [2111.15561].
  * [17] G. Falcioni, F. Herzog, S. Moch and A. Vogt, _Four-loop splitting functions in QCD – The quark-quark case_ , _Phys. Lett. B_ 842 (2023) 137944 [2302.07593].
  * [18] G. Falcioni, F. Herzog, S. Moch and A. Vogt, _Four-loop splitting functions in QCD – The gluon-to-quark case_ , _Phys. Lett. B_ 846 (2023) 138215 [2307.04158].
  * [19] S. Moch, B. Ruijl, T. Ueda, J. Vermaseren and A. Vogt, _Additional moments and x-space approximations of four-loop splitting functions in QCD_ , 2310.05744.
  * [20] G. Falcioni, F. Herzog, S. Moch, A. Pelloni and A. Vogt, _Four-loop splitting functions in QCD – The quark-to-gluon case_ , 2404.09701.
  * [21] Y. Schröder and M. Steinhauser, _Four-loop decoupling relations for the strong coupling_ , _JHEP_ 01 (2006) 051 [hep-ph/0512058].
  * [22] K. G. Chetyrkin, J. H. Kühn and C. Sturm, _QCD decoupling at four loops_ , _Nucl. Phys. B_ 744 (2006) 121 [hep-ph/0512060].
  * [23] M. Buza, Y. Matiounine, J. Smith, R. Migneron and W. L. van Neerven, _Heavy quark coefficient functions at asymptotic values $Q^{2}\gg m^{2}$_, _Nucl. Phys. B_ 472 (1996) 611 [hep-ph/9601302].
  * [24] M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven, _Charm electroproduction viewed in the variable flavor number scheme versus fixed order perturbation theory_ , _Eur. Phys. J. C_ 1 (1998) 301 [hep-ph/9612398].
  * [25] J. Ablinger, J. Blumlein, S. Klein, C. Schneider and F. Wissbrock, _The $O(\alpha_{s}^{3})$ Massive Operator Matrix Elements of $O(n_{f})$ for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_, _Nucl. Phys. B_ 844 (2011) 26 [1008.3347].
  * [26] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round et al., _The transition matrix element $a_{gq}(n)$ of the variable flavor number scheme at $\cal{O}(\alpha_{s}^{3})$_, _Nuclear Physics B_ 882 (2014) 263–288.
  * [27] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider, _The 3-loop pure singlet heavy flavor contributions to the structure function $F_{2}(x,Q^{2})$ and the anomalous dimension_, _Nucl. Phys. B_ 890 (2014) 48 [1409.1135].
  * [28] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider et al., _The first–order factorizable contributions to the three–loop massive operator matrix elements $A_{Qg}^{(3)}$ and $\Delta A_{Qg}^{(3)}$_, _Nucl. Phys. B_ 999 (2024) 116427 [2311.00644].
  * [29] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, C. Schneider et al., _The non-first-order-factorizable contributions to the three-loop single-mass operator matrix elements $A_{Qg}^{(3)}$ and $\Delta A_{Qg}^{(3)}$_, 2403.00513.
  * [30] H. Kawamura, N. A. Lo Presti, S. Moch and A. Vogt, _On the next-to-next-to-leading order QCD corrections to heavy-quark production in deep-inelastic scattering_ , _Nucl. Phys. B_ 864 (2012) 399 [1205.5727].
  * [31] S. Alekhin, J. Blümlein, S. Moch and R. Placakyte, _Parton distribution functions, $\alpha_{s}$, and heavy-quark masses for LHC Run II_, _Phys. Rev. D_ 96 (2017) 014011 [1701.05838].
  * [32] M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven, _Charm electroproduction viewed in the variable flavor number scheme versus fixed order perturbation theory_ , _Eur. Phys. J. C_ 1 (1998) 301 [hep-ph/9612398].
  * [33] J. Ablinger, J. Blümlein, A. De Freitas, C. Schneider and K. Schönwald, _The two-mass contribution to the three-loop pure singlet operator matrix element_ , _Nucl. Phys. B_ 927 (2018) 339 [1711.06717].
  * [34] J. Ablinger, J. Blümlein, A. De Freitas, A. Goedicke, C. Schneider and K. Schönwald, _The Two-mass Contribution to the Three-Loop Gluonic Operator Matrix Element $A_{gg,Q}^{(3)}$_, _Nucl. Phys. B_ 932 (2018) 129 [1804.02226].
  * [35] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Goedicke, A. von Manteuffel et al., _The unpolarized and polarized single-mass three-loop heavy flavor operator matrix elements A gg,Q and $\Delta$Agg,Q_, _JHEP_ 12 (2022) 134 [2211.05462].
  * [36] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, C. Schneider and F. Wißbrock, _Three Loop Massive Operator Matrix Elements and Asymptotic Wilson Coefficients with Two Different Masses_ , _Nucl. Phys. B_ 921 (2017) 585 [1705.07030].
  * [37] APFEL collaboration, _APFEL: A PDF Evolution Library with QED corrections_ , _Comput. Phys. Commun._ 185 (2014) 1647 [1310.1394].
  * [38] V. Bertone, _APFEL++: A new PDF evolution library in C++_ , _PoS_ DIS2017 (2018) 201 [1708.00911].
  * [39] A. Candido, F. Hekhorn and G. Magni, _EKO: evolution kernel operators_ , _Eur. Phys. J. C_ 82 (2022) 976 [2202.02338].
  * [40] G. P. Salam and J. Rojo, _A Higher Order Perturbative Parton Evolution Toolkit (HOPPET)_ , _Comput. Phys. Commun._ 180 (2009) 120 [0804.3755].
  * [41] A. Vogt, _Efficient evolution of unpolarized and polarized parton distributions with QCD-PEGASUS_ , _Comput. Phys. Commun._ 170 (2005) 65 [hep-ph/0408244].
  * [42] M. Botje, _QCDNUM: Fast QCD Evolution and Convolution_ , _Comput. Phys. Commun._ 182 (2011) 490 [1005.1481].
  * [43] J. McGowan, T. Cridge, L. A. Harland-Lang and R. S. Thorne, _Approximate N 3LO parton distribution functions with theoretical uncertainties: MSHT20aN3LO PDFs_, _Eur. Phys. J. C_ 83 (2023) 185 [2207.04739].
  * [44] T. Cridge, L. A. Harland-Lang and R. S. Thorne, _Combining QED and Approximate N 3LO QCD Corrections in a Global PDF Fit: MSHT20qed_an3lo PDFs_, 2312.07665.
  * [45] F. Hekhorn and G. Magni, _DGLAP evolution of parton distributions at approximate N 3LO_, 2306.15294.
  * [46] NNPDF collaboration, _The Path to N 3LO Parton Distributions_, 2402.18635.
  * [47] Y. L. Dokshitzer, G. Marchesini and G. P. Salam, _Revisiting parton evolution and the large-x limit_ , _Phys. Lett. B_ 634 (2006) 504 [hep-ph/0511302].
  * [48] V. S. Fadin, E. A. Kuraev and L. N. Lipatov, _On the Pomeranchuk Singularity in Asymptotically Free Theories_ , _Phys. Lett. B_ 60 (1975) 50.
  * [49] E. A. Kuraev, L. N. Lipatov and V. S. Fadin, _Multi - Reggeon Processes in the Yang-Mills Theory_ , _Sov. Phys. JETP_ 44 (1976) 443.
  * [50] L. N. Lipatov, _Reggeization of the Vector Meson and the Vacuum Singularity in Nonabelian Gauge Theories_ , _Sov. J. Nucl. Phys._ 23 (1976) 338.
  * [51] E. A. Kuraev, L. N. Lipatov and V. S. Fadin, _The Pomeranchuk Singularity in Nonabelian Gauge Theories_ , _Sov. Phys. JETP_ 45 (1977) 199.
  * [52] V. S. Fadin and L. N. Lipatov, _BFKL pomeron in the next-to-leading approximation_ , _Phys. Lett. B_ 429 (1998) 127 [hep-ph/9802290].
  * [53] T. Jaroszewicz, _Gluonic Regge Singularities and Anomalous Dimensions in QCD_ , _Phys. Lett. B_ 116 (1982) 291.
  * [54] M. Ciafaloni and G. Camici, _Energy scale(s) and next-to-leading BFKL equation_ , _Phys. Lett. B_ 430 (1998) 349 [hep-ph/9803389].
  * [55] S. Catani and F. Hautmann, _High-energy factorization and small x deep inelastic scattering beyond leading order_ , _Nucl. Phys. B_ 427 (1994) 475 [hep-ph/9405388].
  * [56] T. Cridge, L. A. Harland-Lang and R. S. Thorne, _The impact of LHC jet and Z $p_{T}$ data at up to approximate N3LO order in the MSHT global PDF fit_, _Eur. Phys. J. C_ 84 (2024) 446 [2312.12505].
  * [57] T. Cridge, L. A. Harland-Lang and R. S. Thorne, _A first determination of the strong coupling $\alpha_{S}$ at approximate N3LO order in a global PDF fit_, 2404.02964.
  * [58] A. Barontini, A. Candido, J. M. Cruz-Martinez, F. Hekhorn and C. Schwan, _Pineline: Industrialization of High-Energy Theory Predictions_ , 2302.12124.
  * [59] J. M. Henn, G. P. Korchemsky and B. Mistlberger, _The full four-loop cusp anomalous dimension in $\mathcal{N}=4$ super Yang-Mills and QCD_, _JHEP_ 04 (2020) 018 [1911.10174].
  * [60] C. Duhr, B. Mistlberger and G. Vita, _Soft integrals and soft anomalous dimensions at N 3LO and beyond_, _JHEP_ 09 (2022) 155 [2205.04493].
  * [61] A. A. Almasy, G. Soar and A. Vogt, _Generalized double-logarithmic large-x resummation in inclusive deep-inelastic scattering_ , _JHEP_ 03 (2011) 030 [1012.3352].
  * [62] J. Davies, C. H. Kom, S. Moch and A. Vogt, _Resummation of small-x double logarithms in QCD: inclusive deep-inelastic scattering_ , _JHEP_ 08 (2022) 135 [2202.10362].
  * [63] R. D. Ball and S. Forte, _The Small x behavior of Altarelli-Parisi splitting functions_ , _Phys. Lett. B_ 465 (1999) 271 [hep-ph/9906222].
  * [64] M. Bonvini and S. Marzani, _Four-loop splitting functions at small $x$_, _JHEP_ 06 (2018) 145 [1805.06460].
  * [65] NNPDF collaboration, _Parton Distributions with Theory Uncertainties: General Formalism and First Phenomenological Studies_ , _Eur. Phys. J. C_ 79 (2019) 931 [1906.10698].
  * [66] NNPDF collaboration, _Determination of the theory uncertainties from missing higher orders on NNLO parton distributions with percent accuracy_ , 2401.10319.
  * [67] A. von Manteuffel, E. Panzer and R. M. Schabinger, _Cusp and collinear anomalous dimensions in four-loop QCD from form factors_ , _Phys. Rev. Lett._ 124 (2020) 162001 [2002.04617].
  * [68] G. Das, S.-O. Moch and A. Vogt, _Soft corrections to inclusive deep-inelastic scattering at four loops and beyond_ , _JHEP_ 03 (2020) 116 [1912.12920].
  * [69] G. Das, S. Moch and A. Vogt, _Approximate four-loop QCD corrections to the Higgs-boson production cross section_ , _Phys. Lett. B_ 807 (2020) 135546 [2004.00563].
  * [70] G. Soar, S. Moch, J. A. M. Vermaseren and A. Vogt, _On Higgs-exchange DIS, physical evolution kernels and fourth-order splitting functions at large x_ , _Nucl. Phys. B_ 832 (2010) 152 [0912.0369].
  * [71] A. Vogt, _Leading logarithmic large-x resummation of off-diagonal splitting functions and coefficient functions_ , _Phys. Lett. B_ 691 (2010) 77 [1005.1606].
  * [72] M. Ciafaloni and D. Colferai, _Dimensional regularisation and factorisation schemes in the BFKL equation at subleading level_ , _JHEP_ 09 (2005) 069 [hep-ph/0507106].
  * [73] M. Ciafaloni, D. Colferai, G. P. Salam and A. M. Stasto, _Minimal subtraction vs. physical factorisation schemes in small-x QCD_ , _Phys. Lett. B_ 635 (2006) 320 [hep-ph/0601200].
  * [74] M. Diehl, R. Nagar and F. J. Tackmann, _ChiliPDF: Chebyshev interpolation for parton distributions_ , _Eur. Phys. J. C_ 82 (2022) 257 [2112.09703].
  * [75] I. Bierenbaum, J. Blumlein and S. Klein, _Mellin Moments of the O(alpha**3(s)) Heavy Flavor Contributions to unpolarized Deep-Inelastic Scattering at Q**2 $>$$>$ m**2 and Anomalous Dimensions_, _Nucl. Phys. B_ 820 (2009) 417 [0904.3563].
  * [76] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel et al., _The 3-Loop Non-Singlet Heavy Flavor Contributions and Anomalous Dimensions for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_, _Nucl. Phys. B_ 886 (2014) 733 [1406.4654].
  * [77] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald, _The three-loop unpolarized and polarized non-singlet anomalous dimensions from off shell operator matrix elements_ , _Nucl. Phys. B_ 971 (2021) 115542 [2107.06267].
  * [78] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round et al., _The $O(\alpha_{s}^{3}T_{F}^{2})$ Contributions to the Gluonic Operator Matrix Element_, _Nucl. Phys. B_ 885 (2014) 280 [1405.4259].
  * [79] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round et al., _3-loop Massive $O(T_{F}^{2})$ Contributions to the DIS Operator Matrix Element $A_{gg}$_, _Nucl. Part. Phys. Proc._ 258-259 (2015) 37 [1409.1435].
  * [80] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round et al., _The transition matrix element A gq(N) of the variable flavor number scheme at O($\alpha$s3)_, _Nuclear Physics B_ 882 (2014) 263 [1402.0359].
  * [81] R. S. Thorne, T. Cridge and L. A. Harland-Lang, _Updates of MSHT PDFs (to appear)_ , 24xx.xxxxx.