\mathrm{NS}}^{\mathrm{FF},\ (3)}$). This presents an issue as it is no longer
possible to rely on $C_{q,q,\ \mathrm{NS}}^{FF,\ (3)}$ to constrain the
low-$Q^{2}$ limit. Nevertheless, by utilising the experience gained from NNLO,
it is feasible to choose any sensible choice for the low-$Q^{2}$ limit. In
practice, due to the observed weak scaling in $Q^{2}$, the exact form at
low-$Q^{2}$ will not present any noticeable differences.

A naive choice for heuristically placing the $C_{q,q,\
\mathrm{NS}}^{\mathrm{FF},\ (3)}(Q^{2}\rightarrow 0)$ function would be a
constant value i.e. no scaling in $Q^{2}$. We propose to use the intuition
from NNLO and the overall fit quality to give us potentially a more sensible
and viable choice for the GM-VFNS approximation141414The differences in fit
quality for sensible choices are $<0.05\%$ compared to the overall $\chi^{2}$
for the light quark NS coefficient function.. By inserting the high-$Q^{2}$
limit into the NS part of Equation (6.2), the result is a crude approximation
to $C_{q,q,\ \mathrm{NS}}^{\mathrm{FF},\ (3)}(Q^{2}\rightarrow 0)$. Combining
this with Equation (6.25h), we obtain a GM-VFNS parameterisation which is
relatively constant across $Q^{2}$ (similar to the NNLO behaviour) with any
differences arising from the subtraction terms which are known.

Figure 19: The N3LO GM-VFNS function $C_{q,q,\ \mathrm{NS}}^{\mathrm{VF},\
(3)}$ compared with the N3LO ZM-VFNS function $C_{q,q,\
\mathrm{NS}}^{\mathrm{ZM},\ (3)}$ across a variety of $x$ and $Q^{2}$ values.
$C_{q,q,\ \mathrm{NS}}^{\mathrm{VF},\ (3)}$ is parameterised via Equations
(6.25h) and (6.2). Mass threshold is set at the charm quark level
($m_{h}^{2}=m_{c}^{2}=1.4\ \text{GeV}^{2}$).

Fig. 19 shows the result of this approximation for the full $C_{q,q,\
\mathrm{NS}}^{\mathrm{VF},\ (3)}$ function. We notice that the behaviour is
similar to that of NNLO across all $(x,Q^{2})$ and appropriately larger in
magnitude to account for the extra contributions obtained at N3LO compared to
NNLO. By definition, the parameterisation converges well to the ZM-VFNS scheme
with the magnitude at high-$Q^{2}$ (ZM-VFNS regime) remaining similar to that
at low-$Q^{2}$ for each specific value of $x$. This final point gives
assurances that even if this low-$Q^{2}$ guess is not entirely representative
of the actual N3LO function, the effects of including this approximation are
virtually negligible in a PDF fit. Also shown in Fig. 19 is the variation in
the $C_{q,q,\ \mathrm{NS}}^{\mathrm{VF},\ (3)}$ function stemming solely from
the $A_{qq,H}^{\mathrm{NS},\ (3)}$ function.

### $C_{q,q}^{\mathrm{PS}}$

To complete the light-quark GM-VFNS coefficient function picture the pure-
singlet contribution from Equation (3) and Equation (3) is described by,

$C_{q,q}^{\mathrm{FF},\ \mathrm{PS}}=C_{q,q}^{\mathrm{VF},\
\mathrm{PS}}\otimes A_{qq,H}^{\mathrm{PS}}\ +\ C_{q,g}^{\mathrm{VF}}\otimes
A_{gq,H}+\ C_{q,H}^{\mathrm{VF},\ \mathrm{PS}}\otimes A_{Hq}$ (6.26)
$\displaystyle\mathcal{O}(\alpha_{s}^{0}):$ $\displaystyle\hskip
28.45274pt\begin{multlined}C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (0)}\ =\
0\end{multlined}C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (0)}\ =\ 0$ (6.27b)
$\displaystyle\mathcal{O}(\alpha_{s}^{1}):$ $\displaystyle\hskip
28.45274pt\begin{multlined}C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (1)}\ =\
0\end{multlined}C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (1)}\ =\ 0$ (6.27d)
$\displaystyle\mathcal{O}(\alpha_{s}^{2}):$ $\displaystyle\hskip
28.45274pt\begin{multlined}C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (2)}\ =\
C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\ (2)}\end{multlined}C_{q,q}^{\mathrm{FF},\
\mathrm{PS},\ (2)}\ =\ C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\ (2)}$ (6.27f)
$\displaystyle\mathcal{O}(\alpha_{s}^{3}):$ $\displaystyle\hskip
28.45274pt\begin{multlined}C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (3)}\ =\
C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\ (3)}+C_{q,g}^{\mathrm{VF},\ (1)}\otimes
A_{gq,H}^{(2)}.\end{multlined}C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (3)}\ =\
C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\ (3)}+C_{q,g}^{\mathrm{VF},\ (1)}\otimes
A_{gq,H}^{(2)}.$ (6.27h)

As with the non-singlet analysis the heavy flavour contributions to the pure-
singlet appear at higher orders to allow for the possibility of heavy quark
contributions. In the pure-singlet case, the heavy flavour contributions are
pushed one order higher than the non-singlet due to the requirement for an
extra intermediary gluon.

#### N3LO

In the pure-singlet case, the FFNS function is non-existent up until N3LO.
Because of this, we choose to parameterise the pure-singlet with a weak
constraint suppressing the FFNS function $C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\
(3)}$ across all $x$ for very low-$Q^{2}$. The reason for this is that the
coefficient functions acquire more contributions as they exist through higher
orders. If $C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\ (3)}$ is beginning at this
order, then one could expect the low-$Q^{2}$ form to be relatively small
compared to the known ZM-VFNS function [57]. This is somewhat justified by the
low-$Q^{2}$ kinematic restrictions for the singlet distribution which broadly
manifest into a suppression at low-$Q^{2}$. We reiterate here that the
low-$Q^{2}$ form of this function is still essentially around the same
magnitude across all $Q^{2}$. Therefore, as with $C_{q,q}^{\mathrm{FF},\
\mathrm{NS},\ (3)}$, it will be virtually negligible in the overall structure
function.

After constructing the approximation for $C_{q,q}^{\mathrm{FF},\ \mathrm{PS},\
(3)}$ with Equation (6.2), Equation (6.27h) is used to approximate the GM-VFNS
function. The exact form of Equation (6.27h) is chosen based on intuition and
where the best fit quality can be achieved151515The differences in fit quality
for sensible choices are $<0.1\%$ compared to the overall $\chi^{2}$ for the
light quark PS coefficient function..

Figure 20: The N3LO GM-VFNS function $C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\
(3)}$ compared with the N3LO ZM-VFNS function $C_{q,q,\
\mathrm{PS}}^{\mathrm{ZM},\ (3)}$ across a variety of $x$ and $Q^{2}$ values.
$C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\ (3)}$ is parameterised via Equations
(6.27h) and (6.2). Mass threshold is set at the charm quark level
($m_{h}^{2}=m_{c}^{2}=1.4\ \text{GeV}^{2}$).

It can be seen from Fig. 20 that the overall magnitude of $C_{q,q,\
\mathrm{PS}}^{\mathrm{VF},\ (3)}$ decreases substantially towards large-$x$ as
one would expect from a pure-singlet function. Inspecting the predicted values
of $C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\ (3)}$, we can confirm that the non-
singlet function from Fig. 19 begins to dominate at large-$x$. Conversely
towards small-$x$, $C_{q,q,\ \mathrm{PS}}^{\mathrm{VF},\ (3)}$ is much larger
than $C_{q,q,\ \mathrm{NS}}^{\mathrm{VF},\ (3)}$, thereby preserving the
familiar interplay between quark distributions. The suppression of the FFNS
parameterisation towards low-$Q^{2}$ is also seen to give sensible results in
terms of the expected percentage change in magnitude through the range of
$Q^{2}$ values. Specifically we see $<10\%$ difference in magnitude between
low and high-$Q^{2}$. Since scale violating terms become more dominant at
higher orders and we are essentially at leading order in terms of heavy
flavour contributions, a high level of scaling with $Q^{2}$ is not expected at
this order.

### $C_{q,g}$

Finally the gluon-light quark coefficient function is constructed from
Equation (3) and Equation (3) to be,

$C_{q,g}^{\mathrm{FF}}=C_{q,q}^{\mathrm{VF}}\otimes A_{qg,H}\ +\
C_{q,g}^{\mathrm{VF}}\otimes A_{gg,H}+\ C_{q,H}^{\mathrm{VF},\
\mathrm{PS}}\otimes A_{Hg}$ (6.28) $\displaystyle\mathcal{O}(\alpha_{s}^{0}):$
$\displaystyle\hskip 28.45274pt\begin{multlined}C_{q,g}^{\mathrm{FF},\ (0)}\
=\ 0\end{multlined}C_{q,g}^{\mathrm{FF},\ (0)}\ =\ 0$ (6.29b)
$\displaystyle\mathcal{O}(\alpha_{s}^{1}):$ $\displaystyle\hskip
28.45274pt\begin{multlined}C_{q,g}^{\mathrm{FF},\ (1)}\ =\
C_{q,g}^{\mathrm{VF},\ (1)}\end{multlined}C_{q,g}^{\mathrm{FF},\ (1)}\ =\
C_{q,g}^{\mathrm{VF},\ (1)}$ (6.29d)
$\displaystyle\mathcal{O}(\alpha_{s}^{2}):$ $\displaystyle\hskip
28.45274pt\begin{multlined}C_{q,g}^{\mathrm{FF},\ (2)}\ =\
C_{q,g}^{\mathrm{VF},\ (2)}+C_{q,g}^{\mathrm{VF},\ (1)}\otimes
A_{gg,H}^{(1)}+A_{qg,H}^{(2)}\end{multlined}C_{q,g}^{\mathrm{FF},\ (2)}\ =\
C_{q,g}^{\mathrm{VF},\ (2)}+C_{q,g}^{\mathrm{VF},\ (1)}\otimes
A_{gg,H}^{(1)}+A_{qg,H}^{(2)}$ (6.29f)
$\displaystyle\mathcal{O}(\alpha_{s}^{3}):$ $\displaystyle\hskip
28.45274pt\begin{multlined}C_{q,g}^{\mathrm{FF},\ (3)}\ =\
C_{q,g}^{\mathrm{VF},\ (3)}+A_{qg,H}^{(3)}+C_{q,q}^{\mathrm{VF},\ (1)}\otimes
A_{qg,H}^{(2)}\\\ +C_{q,g}^{\mathrm{VF},\ (2)}\otimes
A_{gg,H}^{(1)}+C_{q,g}^{\mathrm{VF},\ (1)}\otimes A_{gg,H}^{(2)}\\\ +C_{q,H,\
\mathrm{PS}}^{\mathrm{VF},\ (2)}\otimes
A_{Hg}^{(1)}.\end{multlined}C_{q,g}^{\mathrm{FF},\ (3)}\ =\
C_{q,g}^{\mathrm{VF},\ (3)}+A_{qg,H}^{(3)}+C_{q,q}^{\mathrm{VF},\ (1)}\otimes
A_{qg,H}^{(2)}\\\ +C_{q,g}^{\mathrm{VF},\ (2)}\otimes
A_{gg,H}^{(1)}+C_{q,g}^{\mathrm{VF},\ (1)}\otimes A_{gg,H}^{(2)}\\\ +C_{q,H,\
\mathrm{PS}}^{\mathrm{VF},\ (2)}\otimes A_{Hg}^{(1)}.$ (6.29j)

For $C_{q,g}$, the FFNS function is non-existent up to NNLO, similar to
$C_{q,q,\ \mathrm{NS}}^{\mathrm{FF},(3)}$. However, the $A_{qg,H}$
contribution at NNLO is sub-leading in $n_{f}$ [39] and is therefore not
considered here.

#### N3LO

Figure 21: The N3LO GM-VFNS function $C_{q,g}^{\mathrm{VF},\ (3)}$ compared
with the N3LO ZM-VFNS function $C_{q,g}^{\mathrm{ZM},\ (3)}$ across a variety
of $x$ and $Q^{2}$ values. $C_{q,g}^{\mathrm{VF},\ (3)}$ is parameterised via
Equations (6.29j) and (6.2). Mass threshold is set at the charm quark level
($m_{h}^{2}=m_{c}^{2}=1.4\ \text{GeV}^{2}$).

At N3LO in Equation (6.29j), no information is available for the
$C_{q,g}^{\mathrm{FF},\ (3)}$ at low-$Q^{2}$. Whereas at high-$Q^{2}$ the ZM-
VFNS function is known [57]. To construct the parameterisation, we apply the
same method described for $C_{q,q,\ \mathrm{PS}}^{\mathrm{FF},\ (3)}$.
Specifically, by applying a suppression to the FFNS parameterisation in the
low-$Q^{2}$ limit. After constructing the parameterisation for
$C_{q,g}^{\mathrm{FF},\ (3)}$ with Equation (6.2), Equation (6.29j) is used to
approximate the GM-VFNS function. Since there is no information in the
low-$Q^{2}$ limit, the parameterisation in Equation (6.2) is chosen roughly
based on how the fit prefers the evolution in $Q^{2}$ to behave.

Fig. 21 illustrates the GM-VFNS function in Equation (6.29j) with Equation
(6.2) as $C_{q,g}^{\mathrm{FF},\ (3)}$ across a range of $x$ and $Q^{2}$.
$C_{q,g}^{\mathrm{VF},\ (3)}$ increases in magnitude when moving to smaller
$x$ and by definition converges onto the ZM-VFNS function. The convergence in
this case is chosen to be less steep than for the light quark convergences due
to some minor tensions in the fit161616The differences in fit quality for
sensible choices of Equation (6.2) are $<0.5\%$ compared to the overall
$\chi^{2}$ for the light quark gluon coefficient function.. The magnitude of
$C_{q,g}^{\mathrm{VF},\ (3)}$ across the entire range of $Q^{2}$ is still
relatively constant, although less flat than the behaviour predicted for
$C_{q,q,\ \mathrm{PS/NS}}^{\mathrm{VF},\ (3)}$. However, considering Equation
(6.29j), some justification for this behaviour can be offered. When comparing
the contributions to the FFNS functions in the NS, PS and gluon cases
(Equations (6.25h), (6.27h) and (6.29j) respectively), the $A_{Hg}$ and
$A_{gg,H}$ contributions involved in $C_{q,g}^{\mathrm{FF}}$ are much larger
than the contributions from $A_{gq,H}$, $A_{Hq}$ and $A^{\mathrm{NS}}_{qq,H}$.
Therefore we can expect a larger difference across $Q^{2}$ for the
$C_{q,g}^{\mathrm{VF},\ (3)}$ function. With this being said, the specific
form at low-$Q^{2}$ is not very important in current PDF fits, only that the
form is continuous and valid.

##  7 N3LO $K$-factors

Thus far the primary concern has been the N3LO additions to the theoretical
form of the DIS cross section. However, to complement these changes it is
necessary to extend other cross section data to the same order. With these
ingredients it is possible to maintain a consistent approximate N3LO treatment
across all datasets. At the time of writing, $K$-factors which provide exact
transformations for each dataset up to NNLO are available171717An exception to
this is the CMS $7\ \text{TeV}\ W+c$ [62] dataset where $K$-factors are
available only up to NLO.. Although there has been progress in N3LO
calculations for various processes including Drell-Yan (DY), top production
and Higgs processes [63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76],
there is still missing information on how these $K$-factors behave above NNLO.
In this section we investigate the effects of the $K$-factors for each dataset
when extended to N3LO. Five process categories are considered separately:
Drell-Yan, Jets, $p_{T}$ jets, $t\bar{t}$ production and Dimuon data. Inside
each of these process categories we assume a perfect positive correlation
between the behaviour of datasets i.e. all Drell-Yan $K$-factor shifts from
NNLO are positively correlated. Clearly this treatment is a simplification,
based on the expectation of a high degree of correlation between datasets
concerned with the same processes. In practice, the uncertainty introduced
from including these $K$-factors is already relatively small compared to the
other sources of MHOUs already discussed, therefore any correction to this is
guaranteed to be small (this will be shown more clearly in Section 8).

###  7.1 Extension to aN3LO

The extension to aN3LO is parameterised with a mixture of the NLO and NNLO
$K$-factors. This allows control of the magnitude and shape of the
transformation from NNLO to aN3LO, using the known shifts from lower orders.

The basic idea is presented as,

$K^{\mathrm{N}^{3}\mathrm{LO}/\mathrm{LO}}=a_{\mathrm{NNLO}}\
K^{\mathrm{NNLO}/\mathrm{LO}}+a_{\mathrm{NLO}}\ K^{\mathrm{NLO}/\mathrm{LO}},$
(7.1)

where
$K^{\mathrm{N}^{3}\mathrm{LO}/\mathrm{LO}},K^{\mathrm{NNLO}/\mathrm{LO}}\
\mathrm{and}\ K^{\mathrm{NLO}/\mathrm{LO}}$ are the relevant $K$-factors with
respect to the LO cross section, and $a_{\mathrm{N(N)LO}}$ are variational
parameters controlling the mixture of NNLO and NLO $K$-factors included in the
N3LO $K$-factor approximation. Hence we have 2 parameters for each of the five
processes included in the fit, and now 20 theory nuisance parameters in total
– 10 controlling aN3LO $K$-factors, 5 controlling aN3LO splitting functions
and 5 controlling heavy flavour aN3LO contributions.

To describe this formalism in terms of physical observables we consider the
cross section,

$\sigma=\sigma_{0}+\sigma_{1}+\sigma_{2}+\dots\equiv\sigma_{\mathrm{NNLO}}+\dots,$
(7.2)

where there is an implicit order of $\alpha_{s}^{p+i}$ absorbed into the
definition of $\sigma_{i}$ beginning at the relevant LO for each process, i.e.
$p=0$ for DY.

$K^{\mathrm{NLO/LO}}$ is then the relative shift from $\sigma_{\mathrm{LO}}$
to $\sigma_{\mathrm{NLO}}$,

$K^{\mathrm{NLO/LO}}=\frac{\sigma_{0}+\sigma_{1}}{\sigma_{0}}=1+\frac{\sigma_{1}}{\sigma_{0}}.$
(7.3)

Similarly for NNLO we have,

$K^{\mathrm{NNLO/LO}}=\frac{\sigma_{0}+\sigma_{1}+\sigma_{2}}{\sigma_{0}}=1+\frac{\sigma_{1}}{\sigma_{0}}+\frac{\sigma_{2}}{\sigma_{0}}.$
(7.4)

Moving to N3LO, we write

$\sigma=\sigma_{0}+\sigma_{1}+\sigma_{2}+\sigma_{3}+\dots\equiv\sigma_{\mathrm{N^{3}LO}}+\dots,$
(7.5)

where $\sigma_{3}=a_{1}\sigma_{1}+a_{2}\sigma_{2}$ is approximated as some
superposition of the two lower orders, with $(a_{1},a_{2})=(0,0)$ reproducing
the NNLO case.

Pushing forward with this approximation and using the definitions for
$\sigma_{1,2}$ in terms of $K$-factors (Equations (7.3) and (7.4)) we have,

$\displaystyle\sigma_{\mathrm{N^{3}LO}}$
$\displaystyle=\sigma_{\mathrm{NNLO}}+a_{1}\sigma_{1}+a_{2}\sigma_{2}$ (7.6)
$\displaystyle=\sigma_{\mathrm{NNLO}}+a_{1}\sigma_{0}(K^{\mathrm{NLO/LO}}-1)+a_{2}\sigma_{0}(K^{\mathrm{NNLO/LO}}-K^{\mathrm{NLO/LO}})$
(7.7)

since,

$\displaystyle\sigma_{1}$
$\displaystyle=\sigma_{0}\left(K^{\mathrm{NLO/LO}}-1\right)$ (7.8)
$\displaystyle\sigma_{2}$
$\displaystyle=\sigma_{0}\left(K^{\mathrm{NNLO/LO}}-\sigma_{1}-\sigma_{0}\right)$
$\displaystyle=\sigma_{0}\left(K^{\mathrm{NNLO/LO}}-K^{\mathrm{NLO/LO}}\right).$
(7.9)

From here one can obtain,

$K^{\mathrm{NNLO/LO}}-K^{\mathrm{NLO/LO}}=\frac{\sigma_{2}}{\sigma_{0}}=\frac{\sigma_{2}+\sigma_{0}}{\sigma_{0}}-1\approx\frac{\sigma_{2}+\sigma_{1}+\sigma_{0}}{\sigma_{1}+\sigma_{0}}-1=K^{\mathrm{NNLO/NLO}}-1,$
(7.10)

assuming $\sigma_{1}\ll\sigma_{0}$, which is in general true for a valid
perturbative expansion. Using (7.10) $\sigma_{\mathrm{N^{3}LO}}$ can be
expressed by,

$\sigma_{\mathrm{N^{3}LO}}\simeq\sigma_{\mathrm{NNLO}}\left(1+a_{1}(K^{\mathrm{NLO/LO}}-1)+a_{2}(K^{\mathrm{NNLO/NLO}}-1)\right),$
(7.11)

where $\sigma_{2}\ll\sigma_{1}\ll\sigma_{0}$.

This defines the proposed approximated N3LO cross section. It is given in
terms of extra contributions from lower order shifts, which are controlled by
variational parameters $a_{1}$ and $a_{2}$. It is also true that the
contributions to N3LO are expected to be suppressed by $\alpha_{s}/\pi$ in the
NNLO case and $(\alpha_{s}/\pi)^{2}$ in the NLO case to account for the
strengths of each contribution. Currently this is taken into account within
the variational parameters $a_{1},a_{2}$. However for the purpose of this
description, it is more appropriate to explicitly redefine
$a_{1},a_{2}=a_{s}^{2}\hat{a}_{1},a_{s}\hat{a}_{2}$ where
$a_{s}=\mathcal{N}\alpha_{s}$ and $\mathcal{N}$ is some normalisation factor.
This then results in,

$K^{\mathrm{N^{3}LO/LO}}=K^{\mathrm{NNLO/LO}}\left(1+\hat{a}_{1}\mathcal{N}^{2}\alpha_{s}^{2}(K^{\mathrm{NLO/LO}}-1)+\hat{a}_{2}\mathcal{N}\alpha_{s}(K^{\mathrm{NNLO/NLO}}-1)\right).$
(7.12)

where the LO cross section $\sigma_{0}$ is cancelled and Equation (7.12) is
written in terms of the $K$-factor shifts only. (7.12) also implicitly
includes the correct order $\mathcal{O}(\alpha_{s}^{3})$ in the
parameterisation through (7.3) and (7.4). We can then choose $\mathcal{N}$ in
order to set the approximate magnitude of our variational parameters
$\hat{a}_{1},\hat{a}_{2}$. Given $\alpha_{s}\sim 0.1$ for the processes
considered, if we neglect $\mathcal{N}$ (i.e. choose $\mathcal{N}\sim 1$),
then our order by order reduction in the magnitude of the $K$-factors would be
$\sim 10\%$ for $\mathcal{O}(1)$ for variational parameters, however from
previous orders we see that typically $K$-factors tend to be $30-40\%$ of the
previous order, therefore we instead choose $\mathcal{N}=3$. This then ensures
the natural scale of variation allowed is also of this order with
$\mathcal{O}(1)$ variational parameters describing the admixture of NLO and
NNLO $K$-factors, with conservative penalties applied accordingly.

Reflecting on this, it is worth noting that these fitted $K$-factors will be
sensitive to all orders, not just N3LO. Considering these $K$-factors as
approximating asymptotic behaviour to all orders in perturbation theory when
assessing the stability of predictions, we can be less concerned with any
somewhat large shifts from NNLO to aN3LO, as we will specifically see in the
case of Fig.’s 25 and 26. Finally, we remind the reader that at higher orders,
new terms with more divergent leading logarithms appear which are missed by
the current theoretical description. Due to this, the all-orders asymptotic
description will still remain approximate up to the inclusion of more
divergent leading logarithms in $(x,Q^{2})$ limits at even higher orders.

###  7.2 Numerical Results

Using this formalism for the aN3LO $K$-factors, we present the global fit
results for each of the five process categories considered.

#### Drell-Yan Processes

Figure 22: $K$-factor expansion up to aN3LO shown for the LHCb 2015 $W,\ Z$
dataset [77, 78]. The $K$-factors shown here are absolute i.e. all with
respect to LO ($K^{\mathrm{N^{m}LO/LO}}\ \forall\ m\in\\{1,2,3\\}$).

For the Drell-Yan processes (all calculated at $\mu_{r,f}=m_{ll}/2$), a
reduction of $\sim 1-2\%$ in the $K$-factor shift is predicted across most of
the corresponding datasets at aN3LO. This is in agreement with recent work
[64]. An example of this reduction is shown in Fig. 22.

Figure 23: $K$-factor expansion up to aN3LO shown for the ATLAS 7 TeV high
precision $W,\ Z$ dataset [79]. The $K$-factors shown here are absolute i.e.
all with respect to LO ($K^{\mathrm{N^{m}LO/LO}}\ \forall\ m\in\\{1,2,3\\}$).

Conversely, Fig. 23 displays an example where the $K$-factor shift has much
less of a contribution at N3LO. This is a feature of the ATLAS datasets
included in the fit due to the impact of chosen $p_{T}$ cuts which reduce the
sensitivity to higher orders.

DY Dataset $\chi^{2}$ $\Delta\chi^{2}$ $\Delta\chi^{2}$ from NNLO from NNLO
(NNLO $K$-factors) E866 / NuSea $pp$ DY [80] 223.3 / 184 $-1.8$ $+2.7$ E866 /
NuSea $pd/pp$ DY [81] 8.4 / 15 $-2.0$ $-1.1$ DØ II $Z$ rap. [82] 17.3 / 28
$+0.9$ $+0.6$ CDF II $Z$ rap. [83] 40.5 / 28 $+3.3$ $+1.3$ DØ II
$W\rightarrow\nu\mu$ asym. [84] 16.5 / 10 $-0.8$ $-1.1$ CDF II $W$ asym. [85]
18.2 / 13 $-0.8$ $-0.8$ DØ II $W\rightarrow\nu e$ asym. [86] 30.7 / 12 $-3.2$
$-3.2$ ATLAS $W^{+},\ W^{-},\ Z$ [87] 30.0 / 30 $+0.1$ $+0.3$ CMS W asym.
$p_{T}>35\ \text{GeV}$ [88] 6.7 / 11 $-1.1$ $-1.1$ CMS W asym. $p_{T}>25,30\
\text{GeV}$ [89] 7.7 / 24 $+0.3$ $+0.3$ LHCb $Z\rightarrow e^{+}e^{-}$ [90]
24.1 / 9 $+1.4$ $+0.7$ LHCb W asym. $p_{T}>20\ \text{GeV}$ [91] 12.6 / 10
$+0.1$ $+0.3$ CMS $Z\rightarrow e^{+}e^{-}$ [92] 17.5 / 35 $-0.5$ $-0.3$ ATLAS
High-mass Drell-Yan [93] 18.1 / 13 $-0.8$ $-0.9$ CMS double diff. Drell-Yan
[94] 129.5 / 132 $-15.1$ $+9.7$ LHCb 2015 $W,Z$ [77, 78] 103.9 / 67 $+4.5$
$-0.1$ LHCb $8\text{TeV}$ $Z\rightarrow ee$ [95] 28.8 / 17 $+2.6$ $+1.7$ CMS
$8\ \text{TeV}\ W$ [96] 11.8 / 22 $-0.9$ $+0.1$ ATLAS $7\ \text{TeV}$ high
prec. $W,Z$ [79] 94.5 / 61 $-22.1$ $-8.3$ DØ $W$ asym. [97] 12.2 / 14 $+0.1$
$+1.3$ ATLAS $8\ \text{TeV}$ High-mass DY [98] 63.0 / 48 $+5.9$ $+1.5$ ATLAS
$8\ \text{TeV}\ W$ [99] 58.0 / 22 $+0.4$ $-0.1$ ATLAS $8\ \text{TeV}$ double
diff. $Z$ [100] 91.6 / 59 $+15.7$ $+6.7$ Total 1065.4 / 864 $-12.8$ $+10.4$

Table 2:  Table showing the relevant DY datasets and how the individual
$\chi^{2}$ changes from NNLO by including the N3LO treatment of $K$-factors,
and theoretical N3LO additions discussed earlier. The result with purely NNLO
$K$-factors included for all data in the fit is also given.

Table 2 demonstrates that in most cases, the new fitted DY aN3LO $K$-factors
are producing a slightly better fit with a moderate cumulative effect. We
remind the reader that we have included a total of 20 extra parameters into
the fit. These extra 20 parameters are fit across all datasets and multiple
processes, whereas the decrease here is for a subset of datasets corresponding
to the DY processes included in a global fit.

Across these datasets, the $K$-factors act to extend the description of these
processes to approximate N3LO. The result of including this procedure is a
better fit in the DY regime while also relaxing tensions with other processes
included in the fit. Comparing the $\Delta\chi^{2}$ results with and without
aN3LO $K$-factors, we can see the extent to which the $K$-factors and all
other N3LO additions are reducing the overall $\chi^{2}$.

In some individual cases, the dataset $\chi^{2}$ becomes somewhat worse
relative to NNLO most notably the ATLAS $8\ \text{TeV}$ double differential
$Z$ [100], whilst in a few others the $\chi^{2}$ improvement upon addition of
the aN3LO splitting functions, transition matrix elements and coefficient
function pieces is seen to deteriorate upon addition of the aN3LO $K$-factors,
e.g. the LHCb 2015 $W,Z$ [77, 78], which exhibits a mild preference for the
N3LO theory with NNLO $K$-factors. The addition of the aN3LO $K$-factors do
nonetheless result in a net reduction in $\chi^{2}$ and for a large number of
cases the aN3LO $K$-factors allow for a slight reduction in the individual
$\chi^{2}$. Of particular note is the $7\ \text{TeV}$ high precision $W,Z$
[79] data, which improves by over 20 points. This may indicate that the
challenge in achieving a statistically good fit to this high precision data
that is observed across all NNLO PDF fits is in part related to the lack of
higher order corrections in the theory. On the other hand, for arguably the
other most precise (and multi-differential) data, namely the ATLAS $8\
\text{TeV}$ double differential $Z$ [100], the fit quality deteriorates. This
may be due to the differing mass binning and/or cuts, but in general it is
difficult to draw firm conclusions here, at least with the current $K$-factor
treatment. The CMS double diff. Drell-Yan [94] also shows a particularly large
reduction when these are added on top of the aN3LO theory, this is a dataset
which shows some tension with the DIS N3LO additions which is then eased by
the addition of the aN3LO $K$-factors.

#### Jet Production Processes

Figure 24: $K$-factor expansion up to aN3LO shown for the CMS 7 TeV jets
dataset ($R=0.7$) [101]. The $K$-factors shown here are absolute i.e. all with
respect to LO ($K^{\mathrm{N^{m}LO/LO}}\ \forall\ m\in\\{1,2,3\\}$).

The jets processes (all calculated for $\mu_{r,f}=p_{T}^{jet}$) show a general
increase in the $K$-factor shifts from NNLO as seen in Fig. 24, which displays
the $K$-factor expansion up to aN3LO for the CMS 7 TeV jets dataset [101]. It
is apparent that there is a mild shift to N3LO from the NNLO $K$-factor. This
behaviour follows what one might expect for a perturbative expansion
considering the forms of the NLO and NNLO functions.

Jets Dataset $\chi^{2}$ $\Delta\chi^{2}$ $\Delta\chi^{2}$ from NNLO from NNLO
(NNLO $K$-factors) CDF II $p\bar{p}$ incl. jets [102] 66.5 / 76 $+6.1$ $+1.7$
DØ II $p\bar{p}$ incl. jets [103] 113.3 / 110 $-6.9$ $-3.0$ ATLAS $7\
\text{TeV}$ jets [104] 215.9 / 140 $-5.6$ $+3.7$ CMS $7\ \text{TeV}$ jets
[101] 186.8 / 158 $+11.0$ $+9.3$ CMS $8\ \text{TeV}$ jets [105] 271.3 / 174
$+10.0$ $+21.5$ CMS $2.76\ \text{TeV}$ jet [106] 109.8 / 81 $+6.9$ $+9.0$
Total 963.6 / 739 $+21.5$ $+42.2$

Table 3: Table showing the relevant jet datasets and how the individual
$\chi^{2}$ changes from NNLO by including the N3LO treatment of $K$-factors.
The result with purely NNLO $K$-factors included for all data in the fit is
also given.

A $\chi^{2}$ summary of the Jets datasets is provided in Table 3. By combining
the N3LO structure function and DGLAP additions (Section’s 4 to 6) with NNLO
$K$-factors, the fit exhibits a substantial increase in the $\chi^{2}$ from
Jets data. Including aN3LO $K$-factors acts to reduce some of this tension
with around half the initial overall $\chi^{2}$ increase still remaining. We
note that in the case of the ATLAS $7\ \text{TeV}$ jets [104], it is well
known that there are issues in achieving a good fit quality across all
rapidity bins (see [107] for a detailed study as well as [108] where the 8 TeV
data are presented and the same issues observed). In [107, 108] the
possibility of decorrelating some of the systematic error sources where the
degree of correlation is less well established, was considered and indeed in
our study we follow such a procedure, as described in [3]. Alternatively,
however, it might be that the issues in fit quality could at least part be due
to deficiencies in theoretical predictions, such as MHOs. To assess this, we
revert to the default ATLAS correlation scenario and repeat the global fit. We
find that the $\chi^{2}$ deteriorates by $+40.7$ points to 256.6, which is
very close to the result found in a pure NNLO fit [3]. In other words, in our
framework the impact of MHOUs does not resolve this issue.

The $\chi^{2}$ results for datasets in Table 3 show evidence for some tensions
with the N3LO form of the high-$x$ gluon. It is also apparent that the CMS
data is in more tension than ATLAS datasets with N3LO structure function and
DGLAP theory. Therefore it will be interesting to see how this behaviour
changes when considering this data as dijets in the global fit [109]. We do
not consider the dijet data here, though this will be addressed in a future
publication.

####  $Z\ p_{T}$ & Vector Boson $+$ Jets Processes

Figure 25: $K$-factor expansion up to aN3LO shown for the ATLAS 8 TeV $Z\
p_{T}$ dataset [110]. The $K$-factors shown here are absolute i.e. all with
respect to LO ($K^{\mathrm{N^{m}LO/LO}}\ \forall\ m\in\\{1,2,3\\}$).

In the case of $Z\ p_{T}$ & vector boson $+$ jet processes (all calculated at
$\mu_{r,f}=\sqrt{p_{T,ll}^{2}+m_{ll}^{2}}$), the $K$-factor shift is almost
completely dominated by the ATLAS 8 TeV $Z\ p_{T}$ dataset [110] (due to the
larger number of data points included in this dataset) shown in Fig. 25. The
gluon is less directly constrained than the quarks in a global fit. Therefore
it can be expected that the significant modifications at small-$x$ will
indirectly affect the high-$x$ gluon, where these processes are most
sensitive. Considering the jet production processes in Table 3, when
performing separate PDF fits not including ATLAS 8 TeV $Z\ p_{T}$ data [110],
we find a reduction of $\Delta\chi^{2}=-7.0$ in CMS $8\ \text{TeV}$ jets data
[105] eliminating most of the tension for this dataset (similar to MSHT20 NNLO
results in Table 16 of [3]). Further to this, when not including HERA and
ATLAS 8 TeV $Z\ p_{T}$ data we find a reduction of $\Delta\chi^{2}=-26.4$ in
CMS $8\ \text{TeV}$ jet data [105] and $\Delta\chi^{2}=-12.7$ in CMS $2.76\
\text{TeV}$ jet data [106].

Although the overall magnitude of the $K$-factor in Fig. 25 may seem large,
this new shift is contained within a $15\%$ increase from NNLO (due to the NLO
and NNLO $K$-factors also being significant). Moreover, not only does the size
of this shift have some dependence on the central scale, but this shift may be
more correctly interpreted as the preferred all-orders cross section rather
than simply the pure ${\rm N}^{3}{\rm LO}$ result.

$p_{T}$ Jets Dataset $\chi^{2}$ $\Delta\chi^{2}$ $\Delta\chi^{2}$ from NNLO
from NNLO (NNLO $K$-factors) CMS $7\ \text{TeV}\ W+c$ [62] 10.8 / 10 $+2.2$
$+0.2$ ATLAS $8\ \text{TeV}\ Z\ p_{T}$ [110] 108.4 / 104 $-80.0$ $-54.5$ ATLAS
$8\ \text{TeV}\ W+\text{jets}$ [111] 18.8 / 30 $+0.7$ $-0.3$ Total 138.0 / 144
$-77.2$ $-54.7$

Table 4: Table showing the relevant $Z\ p_{T}$ & Vector Boson jet datasets and
how the individual $\chi^{2}$ changes from NNLO by including the N3LO
treatment of $K$-factors, and theoretical N3LO additions discussed earlier.
The result with purely NNLO $K$-factors included for all data in the fit is
also given.

The extent of the $\chi^{2}$ reduction in the $Z\ p_{T}$ datasets is shown in
Table 4. Note that around $\sim 68\%$ of the improvement to the ATLAS $8\
\text{TeV}\ Z\ p_{T}$ [110] $\chi^{2}$ is due to the extra N3LO theory
included in the DGLAP and DIS descriptions. It is also known the ATLAS $8\
\text{TeV}\ Z\ p_{T}$ data [110] previously exhibited a significant level of
tension with many datasets (including HERA data) at NNLO [3]. This was
investigated by performing a global PDF fit with and without HERA data and
comparing the individual $\chi^{2}$’s from each dataset. At NNLO it was found
that the ATLAS $8\ \text{TeV}\ Z\ p_{T}$ dataset [110] reduced by
$\Delta\chi^{2}=-39.2$ when fitting to all non-HERA data (see Table B.1). At
aN3LO we observe that the ATLAS $8\ \text{TeV}\ Z\ p_{T}$ dataset [110]
actually increased by $\Delta\chi^{2}=+12.8$ when fitting to all non-HERA data
(see Table B.3). The aN3LO additions therefore eliminate this tension
previously observed at NNLO, suggesting that this issue at NNLO in fitting the
ATLAS $8\text{TeV}Z\ p_{T}$ dataset [110] was a sign of MHOs. This is in
contrast with the result observed for ATLAS $7\ \text{TeV}$ jets [104] where
the issues with fit quality were not alleviated by the inclusion of known
higher order N3LO information and approximations for the remaining missing
pieces.

Finally we remind the reader that the CMS $7\ \text{TeV}\ W+c$ dataset [62]
does not include a $K$-factor at NNLO. To overcome this, we tie the overall
N3LO $K$-factor shift to the NLO value ($K^{\mathrm{NNLO/NLO}}=1$ in Equation
(7.12)), therefore contributing as an overall normalisation effect.

#### Top Quark Processes

Figure 26: $K$-factor expansion up to aN3LO shown for the CMS 8 TeV single
diff. $t\bar{t}$ dataset [112]. The $K$-factors shown here are absolute i.e.
all with respect to LO ($K^{\mathrm{N^{m}LO/LO}}\ \forall\ m\in\\{1,2,3\\}$).

Moving to top quark processes, for the single differential datasets the scale
choice for $\mu_{r,f}$ is $H_{T}/4$ with the exception of data differential in
the average transverse momentum of the top or anti-top,
$p_{T}^{t},p_{T}^{\bar{t}}$, for which $m_{T}/2$ is used. For the double diff.
dataset the scale choice is $H_{T}/4$ and for the inclusive top
$\sigma_{t\bar{t}}$ a scale of $m_{t}$ is chosen. Fig. 26 displays the
$K$-factor shifts up to N3LO for the CMS 8 TeV single diff. $t\bar{t}$ dataset
[112], which shows the greatest reduction in its $\chi^{2}$. A familiar
perturbative pattern can be seen for this process’s $K$-factors, with the
shift at aN3LO increasing by around 3–4% from NNLO. This is in agreement with
a recent $\sim 3.5\%$ predicted increase in the N3LO $t\bar{t}$ production
$K$-factor at $8\ \text{TeV}$ in [66], whereby an approximate N3LO cross
section for $t\bar{t}$ production in proton-proton collisions has been
calculated employing a resummation formalism [113, 114, 115, 116].

Top Dataset $\chi^{2}$ $\Delta\chi^{2}$ $\Delta\chi^{2}$ from NNLO from NNLO
(NNLO $K$-factors) Tevatron, ATLAS, CMS $\sigma_{t\bar{t}}$ [117, 118, 119,
120, 121, 122, 123, 124, 125, 126, 127, 128, 129] 14.3 / 17 $-0.2$ $-1.0$
ATLAS $8\ \text{TeV}$ single diff. $t\bar{t}$ [130] 24.2 / 25 $-1.4$ $-0.8$
ATLAS $8\ \text{TeV}$ single diff. $t\bar{t}$ dilep. [131] 2.7 / 5 $-0.7$
$-0.7$ CMS $8\ \text{TeV}$ double diff. $t\bar{t}$ [132] 23.6 / 15 $+1.0$
$+4.0$ CMS $8\ \text{TeV}$ single diff. $t\bar{t}$ [112] 10.3 / 9 $-2.9$
$-4.0$ Total 75.1 / 71 $-4.2$ $-2.5$

Table 5: Table showing the relevant Top Quark datasets and how the individual
$\chi^{2}$ changes from NNLO by including the N3LO treatment of $K$-factors,
and theoretical N3LO additions discussed earlier. The result with purely NNLO
$K$-factors included for all data in the fit is also given.

The $\chi^{2}$ results in Table 5 display a mildly better fit for top
processes, with most datasets not feeling a large overall effect from the N3LO
additions. Comparing with and without aN3LO $K$-factors, we see a slightly
better fit overall, with most of the reduction in overall $\chi^{2}$ stemming
from CMS $8\ \text{TeV}$ double diff. $t\bar{t}$ data [132].

#### Semi-Inclusive DIS Dimuon Processes

Figure 27: $K$-factor expansion up to aN3LO shown for the NuTeV $\nu
N\rightarrow\mu\mu X$ dataset [133]. The $K$-factors shown here are absolute
i.e. all with respect to LO ($K^{\mathrm{N^{m}LO/LO}}\ \forall\
m\in\\{1,2,3\\}$).

The final set of results to consider in this Section are the aN3LO $K$-factors
associated with semi-inclusive DIS dimuon cross sections (with
$\mu_{r,f}^{2}=Q^{2}$). Although the dimuon cross section is associated with
the DIS process described from our approximate N3LO structure function
picture, it is a semi-inclusive DIS process. Therefore it is sensible to treat
this process as entirely separate from DIS. The NNLO cross-sections used in
this case are a general-mass variable flavour number scheme extension of the
results in [134], as described in more detail in [3]. The $K$-factors shown in
Fig. 27 (for the NuTeV $\nu N\rightarrow\mu\mu X$ data [133]) are somewhat
similar to NNLO. The reason for this is mostly due to these datasets also
including a branching ratio (${\rm BR}(D\to\mu)$) which absorbs any overall
normalisation shifts. This behaviour is not a concern since in practice these
two work in tandem and when combined together it makes no difference where the
normalisation factors are absorbed into.

NLO NNLO aN3LO aN3LO (NNLO $K$-factors) 0.095 0.088 0.082 0.081 ${\rm
BR}(D\to\mu)$

Table 6: Table displaying dimuon branching ratios (BRs) at NLO, NNLO, aN3LO
and aN3LO with NNLO $K$-factors.

Investigating the change in the BR’s with the addition of N3LO contributions
in Table 6, the BR at N3LO decreases substantially from NNLO, with little
difference from the addition of aN3LO $K$-factors. The predicted dimuon BR at
aN3LO is inside the allowed $\pm 1\sigma$ range of $0.092\pm 0.010$. When
performing a fit with the BR fixed at its central value (BR = $0.092$), one is
able to observe the effect of manually forcing the normalisation into the
$K$-factor variation alone. The result of this is a worse global fit quality
$\Delta\chi^{2}=+11.2$, where $+3.9$ units arise from an increased penalty for
the Dimuon $K$-factor description and $+2.3$ units from a slightly worse fit
to the Dimuon datasets listed in Table 7. The rest of the observed increase in
$\chi^{2}_{\mathrm{global}}$ is dominated by a $+4.1$ increase in the ATLAS
$7\ \text{TeV}$ high prec. $W,Z$ [79] data due to a smaller strange quark PDF
(compensating the higher BR in dimuon datasets). Returning to consider the
case of the $K$-factors and BR together, the predicted effect on dimuon
datasets is very similar. However due to the errors accounting for a larger
allowed shift in the BR relative to the $K$-factors, the fit favours moving
the BR by a larger amount to reduce the penalty $\chi^{2}$ contribution from
$K$-factors which explains the results shown in Table 6.

Dimuon Dataset $\chi^{2}$ $\Delta\chi^{2}$ $\Delta\chi^{2}$ from NNLO from
NNLO (NNLO $K$-factors) CCFR $\nu N\rightarrow\mu\mu X$ [133] 68.3 / 86 $+0.6$
$+1.1$ NuTeV $\nu N\rightarrow\mu\mu X$ [133] 56.7 / 84 $-1.8$ $-0.6$ Total
125.0 / 170 $-1.2$ $+0.5$

Table 7: Table showing the relevant Dimuon datasets and how the individual
$\chi^{2}$ changes from NNLO by including the N3LO treatment of $K$-factors,
and theoretical N3LO additions discussed earlier. The result with purely NNLO
$K$-factors included for all data in the fit is also given.

Table 7 further confirms the expectation that the Dimuon datasets are not too
sensitive to N3LO additions. The results with and without a full treatment of
aN3LO $K$-factors are also similar in magnitude. It is therefore clear that
the dimuon BR’s are compensating for any indirect normalisation effects from
the form of the PDFs in the full aN3LO fit, as opposed to the aN3LO
$K$-factors.

##  8 MSHT20 Approximate N3LO Global Analysis

With the inclusion of all N3LO approximations discussed in earlier sections
resulting in 20 extra free parameters from the NNLO MSHT20 fit, we now present
the results for the first approximate N3LO global PDF fit with theoretical
uncertainties from missing N3LO contributions and implicitly some MHOs beyond
this. This includes the results for the best fit for the nuisance parameters
describing the theoretical uncertainty. We remind the reader that these are
parameterised specifically to represent the missing uncertainty at N3LO, which
is currently the dominant source of uncertainty due to missing higher orders.
However, the fit will also be influenced, to a limited extent, by effects at
even higher orders. Later in the Section we discuss this in more detail.

###  8.1 $\chi^{2}$ Breakdown

Dataset $N_{\mathrm{pts}}$ $\chi^{2}$ $\Delta\chi^{2}$ from NNLO BCDMS $\mu p$
$F_{2}$ [135] 163 174.4 $-5.8$ BCDMS $\mu d$ $F_{2}$ [135] 151 144.3 $-1.7$
NMC $\mu p$ $F_{2}$ [136] 123 121.5 $-2.6$ NMC $\mu d$ $F_{2}$ [136] 123 104.2
$-8.4$ SLAC $ep$ $F_{2}$ [137, 138] 37 31.6 $-0.4$ SLAC $ed$ $F_{2}$ [137,
138] 38 22.8 $-0.2$ E665 $\mu d$ $F_{2}$ [139] 53 63.9 $+4.2$ E665 $\mu p$
$F_{2}$ [139] 53 67.5 $+2.9$ NuTeV $\nu N$ $F_{2}$ [140] 53 35.7 $-2.6$ NuTeV
$\nu N$ $xF_{3}$ [140] 42 34.8 $+4.1$ NMC $\mu n/\mu p$ [141] 148 131.6 $+0.8$
E866 / NuSea $pp$ DY [80] 184 223.3 $-1.8$ E866 / NuSea $pd/pp$ DY [81] 15 8.4
$-2.0$ HERA $ep$ $F_{2}^{\text{charm}}$ [142] 79 143.7 $+11.4$
NMC/BCDMS/SLAC/HERA $F_{L}$ [135, 136, 138, 143, 144, 145] 57 45.6 $-22.9$
CCFR $\nu N\rightarrow\mu\mu X$ [133] 86 68.3 $+0.6$ NuTeV $\nu
N\rightarrow\mu\mu X$ [133] 84 56.7 $-1.8$ CHORUS $\nu N$ $F_{2}$ [146] 42
29.2 $-1.0$ CHORUS $\nu N$ $xF_{3}$ [146] 28 18.1 $-0.3$ HERA $e^{+}p$ CC
[147] 39 49.7 $-2.3$ HERA $e^{-}p$ CC [147] 42 64.9 $-5.3$ HERA $e^{+}p$ NC
$820\ \text{GeV}$ [147] 75 84.3 $-5.6$ HERA $e^{-}p$ NC $460\ \text{GeV}$
[147] 209 247.7 $-0.6$ HERA $e^{+}p$ NC $920\ \text{GeV}$ [147] 402 474.0
$-38.7$ HERA $e^{-}p$ NC $575\ \text{GeV}$ [147] 259 248.5 $-14.5$ HERA
$e^{-}p$ NC $920\ \text{GeV}$ [147] 159 243.0 $-1.4$ CDF II $p\bar{p}$ incl.
jets [102] 76 66.5 $+6.1$ DØ II $Z$ rap. [82] 28 17.3 $+0.9$ CDF II $Z$ rap.
[83] 28 40.5 $+3.3$ DØ II $W\rightarrow\nu\mu$ asym. [84] 10 16.5 $-0.8$ CDF
II $W$ asym. [85] 13 18.2 $-0.8$

Table 8: Full breakdown of $\chi^{2}$ results for the aN3LO PDF fit. The
global fit includes the N3LO treatment for transition matrix elements,
coefficient functions, splitting functions and $K$-factor additions with their
variational parameters determined by the fit.

Dataset $N_{\mathrm{pts}}$ $\chi^{2}$ $\Delta\chi^{2}$ from NNLO DØ II
$W\rightarrow\nu e$ asym. [86] 12 30.7 $-3.2$ DØ II $p\bar{p}$ incl. jets
[103] 110 113.3 $-6.9$ ATLAS $W^{+},\ W^{-},\ Z$ [87] 30 30.0 $+0.1$ CMS W
asym. $p_{T}>35\ \text{GeV}$ [88] 11 6.7 $-1.1$ CMS W asym. $p_{T}>25,30\
\text{GeV}$ [89] 24 7.7 $+0.3$ LHCb $Z\rightarrow e^{+}e^{-}$ [90] 9 24.1
$+1.4$ LHCb W asym. $p_{T}>20\ \text{GeV}$ [91] 10 12.6 $+0.1$ CMS
$Z\rightarrow e^{+}e^{-}$ [92] 35 17.5 $-0.5$ ATLAS High-mass Drell-Yan [93]
13 18.1 $-0.8$ Tevatron, ATLAS, CMS $\sigma_{t\bar{t}}$ [117, 118, 119, 120,
121, 122, 123, 124, 125, 126, 127, 128, 129] 17 14.3 $-0.2$ CMS double diff.
Drell-Yan [94] 132 129.5 $-15.1$ LHCb 2015 $W,Z$ [77, 78] 67 103.9 $+4.5$ LHCb
$8\text{TeV}$ $Z\rightarrow ee$ [95] 17 28.8 $+2.6$ CMS $8\ \text{TeV}\ W$
[96] 22 11.8 $-0.9$ ATLAS $7\ \text{TeV}$ jets [104] 140 215.9 $-5.6$ CMS $7\
\text{TeV}\ W+c$ [62] 10 10.8 $+2.2$ ATLAS $7\ \text{TeV}$ high prec. $W,Z$
[79] 61 94.5 $-22.1$ CMS $7\ \text{TeV}$ jets [101] 158 186.8 $+11.0$ DØ $W$
asym. [97] 14 12.2 $+0.1$ ATLAS $8\ \text{TeV}\ Z\ p_{T}$ [110] 104 108.4
$-80.0$ CMS $8\ \text{TeV}$ jets [105] 174 271.3 $+10.0$ ATLAS $8\ \text{TeV}$
sing. diff. $t\bar{t}$ [130] 25 24.2 $-1.4$ ATLAS $8\ \text{TeV}$ sing. diff.
$t\bar{t}$ dilep. [131] 5 2.7 $-0.7$ ATLAS $8\ \text{TeV}$ High-mass DY [98]
48 63.0 $+5.9$ ATLAS $8\ \text{TeV}\ W+\text{jets}$ [111] 30 18.8 $+0.7$ CMS
$8\ \text{TeV}$ double diff. $t\bar{t}$ [132] 15 23.6 $+1.0$ ATLAS $8\
\text{TeV}\ W$ [99] 22 58.0 $+0.4$ CMS $2.76\ \text{TeV}$ jet [106] 81 109.8
$+6.9$ CMS $8\ \text{TeV}$ sing. diff. $t\bar{t}$ [112] 9 10.3 $-2.9$ ATLAS
$8\ \text{TeV}$ double diff. $Z$ [100] 59 91.6 $+15.7$

Table 9: (Continued) Full breakdown of $\chi^{2}$ results for the aN3LO PDF
fit. The global fit includes the N3LO treatment for transition matrix
elements, coefficient functions, splitting functions and $K$-factor additions
with their variational parameters determined by the fit.

Low-$Q^{2}$ Coefficient $c_{q}^{\mathrm{NLL}}$ $=-3.868$ 0.004
$c_{g}^{\mathrm{NLL}}$ $=-5.837$ 0.844 Transition Matrix Elements $a_{Hg}$
$=12214.000$ 0.601 $a_{qq,H}^{\mathrm{NS}}$ $=-64.411$ 0.001 $a_{gg,H}$
$=-1951.600$ 0.857 Splitting Functions $\rho_{qq}^{NS}$ $=0.007$ 0.000
$\rho_{gq}$ $=-1.784$ 0.802 $\rho_{qq}^{PS}$ $=-0.501$ 0.186 $\rho_{gg}$
$=19.245$ 3.419 $\rho_{qg}$ $=-1.754$ 0.015 K-factors
$\mathrm{DY}_{\mathrm{NLO}}$ $=-0.282$ 0.080 $\mathrm{DY}_{\mathrm{NNLO}}$
$=0.079$ 0.006 $\mathrm{Top}_{\mathrm{NLO}}$ $=0.041$ 0.002
$\mathrm{Top}_{\mathrm{NNLO}}$ $=0.651$ 0.424 $\mathrm{Jet}_{\mathrm{NLO}}$
$=-0.300$ 0.090 $\mathrm{Jet}_{\mathrm{NNLO}}$ $=-0.691$ 0.478
$p_{T}\mathrm{Jets}_{\mathrm{NLO}}$ $=0.583$ 0.339
$p_{T}\mathrm{Jets}_{\mathrm{NNLO}}$ $=-0.080$ 0.006
$\mathrm{Dimuon}_{\mathrm{NLO}}$ $=-0.444$ 0.197
$\mathrm{Dimuon}_{\mathrm{NNLO}}$ $=0.922$ 0.850 N3LO Penalty Total 9.201 / 20
Average Penalty 0.460 Total 4957.2 / 4363 $\Delta\chi^{2}$ from NNLO $-154.4$

Table 10: (Continued) Full breakdown of $\chi^{2}$ results for the aN3LO PDF
fit. The global fit includes the N3LO treatment for transition matrix
elements, coefficient functions, splitting functions and $K$-factor additions
with their variational parameters determined by the fit.

Table 8 shows the global $\chi^{2}$ results for an aN3LO best fit, inclusive
of penalties associated with the new theory variational parameters (from
Equation (2.17)). The theory parameters are labelled as: $A_{Hg}(a_{Hg})$,
$A_{gg,H}(a_{gg,H})$, $A_{qq,H}^{\mathrm{NS}}(a_{qq,H}^{\mathrm{NS}})$ for the
transition matrix elements; $P_{qq}^{\mathrm{NS}}(\rho_{qq}^{\mathrm{NS}})$,
$P_{qq}^{\mathrm{PS}}(\rho_{qq}^{\mathrm{PS}})$, $P_{qg}(\rho_{qg})$,
$P_{gq}(\rho_{gq})$ and $P_{gg}(\rho_{gg})$ for the splitting functions; and
$c_{q}^{\mathrm{NLL}}$ and $c_{g}^{\mathrm{NLL}}$ correspond to the NLL
parameters discussed in Section 6.2. These are supplemented by the 10
additional nuisance parameters for the NLO and NNLO $K$-factors for the five
process categories. These 20 additional parameters and their associated
penalties are also shown in Table 8.

The extra N3LO theory and level of freedom introduced has allowed the fit to
achieve a total $\Delta\chi^{2}=-154.4$ compared to MSHT20 NNLO total
$\chi^{2}$ (Table 7 from [3]). Comparing with lower order PDF fits, we find a
smooth convergence in the fit quality which follows what one may expect from
an increase in the accuracy of a perturbative expansion
($\chi^{2}/N_{\mathrm{pts}}$ = LO: 2.57, NLO: 1.33, NNLO: 1.17, N3LO: 1.14).
In part, this is due to the extra freedom in the $K$-factors, which will
almost always act to reduce this $\chi^{2}$ due to the minimisation procedure.
However, even with this freedom, in most cases the N3LO theory (non
$K$-factor) contributions include large divergences from NNLO. With this in
mind, we must conclude that the fit is preferring a description different from
the current NNLO standard.

At NNLO (Table B.1), the tension between HERA and non-HERA datasets accounted
for $\Delta\chi^{2}=-61.6$ reduction in the overall fit quality when the
former was removed, with the majority of this tension between HERA and ATLAS
$8\ \text{TeV}\ Z\ p_{T}$ [110] data. Whereas comparing fit results with and
without HERA data at N3LO, we find $\Delta\chi^{2}=-49.0$. Although the
overall difference is not too substantial we do report a substantial shift in
the leading tensions, where most of the tension with HERA data is now residing
with NMC $F_{2}$ [136] and CMS $8\ \text{TeV}$ jets [105] data. Tensions with
NMC $F_{2}$ [136] data are also seen to some extent at NNLO where we show a
$\Delta\chi^{2}=-20.6$ in a fit omitting HERA data (combining the NMC $F_{2}$
datasets shown in Table B.1). However at N3LO, Table B.3 shows a
$\Delta\chi^{2}=-24.4$ reduction from NMC $F_{2}$ data in a fit omitting HERA
data. Therefore whilst the N3LO additions remove tensions with $Z\ p_{T}$
data, it remains that the HERA data is preferring the high-$x$ quarks to be
lower than favoured by NMC data. This is suggestive of higher twist effects
for NMC data at low-$Q^{2}$ (as we observe a worse fit to low-$Q^{2}$ data).
We also emphasise that when conducting a fit at NNLO with $Z\ p_{T}$ data
removed, an improvement of $\Delta\chi^{2}=-41.3$ is observed in the rest of
the data, whereas at N3LO an improvement of $\Delta\chi^{2}=-65.2$ is observed
in all other datasets without removing $Z\ p_{T}$, therefore these results are
not purely an effect of removing any $Z\ p_{T}$ tension. Considering tensions
with CMS $8\ \text{TeV}$ jets [105] data, as discussed in Section 7, in
general the jets datasets show tensions with the N3LO description (especially
for CMS $8\ \text{TeV}$ jets [105]), therefore it will be interesting to
observe how this picture evolves when considering this data in the form of
dijets.

Since a naturally richer description of the small-$x$ regime is being included
at N3LO, which has a direct effect on the HERA datasets, the reduction of
important tensions from NNLO is even further justification for the inclusion
of the N3LO theory. The extra N3LO additions are allowing the large-$x$
behaviour of the PDFs to be less dominated by data at small-$x$, while also
producing a better fit quality at small-$x$ (i.e. for HERA data). Some of the
above observations are also made in [16, 17] where studies of including
small-$x$ resummation results into a PDF fit have been reported.

Reflecting on the chosen prior distributions for each of the sources of
missing N3LO uncertainty, Table 8 confirms that no especially large penalties
are being incurred in this new description. These results therefore
demonstrate that the fit is succeeding in leveraging contributions (such as
$P_{qq}^{(3)}$ and $P_{qg}^{(3)}$ in the quark evolution part of Equation
(3.3)) to produce a better overall fit.

#### DIS Processes

DIS Dataset $\chi^{2}$ $\Delta\chi^{2}$ $\Delta\chi^{2}$ from NNLO from NNLO
(NNLO $K$-factors) BCDMS $\mu p$ $F_{2}$ [135] 174.4 / 163 $-5.8$ $-5.6$ BCDMS
$\mu d$ $F_{2}$ [135] 144.3 / 151 $-1.7$ $+0.6$ NMC $\mu p$ $F_{2}$ [136]
121.5 / 123 $-2.6$ $-3.8$ NMC $\mu d$ $F_{2}$ [136] 104.2 / 123 $-8.4$ $-11.5$
SLAC $ep$ $F_{2}$ [137, 138] 31.6 / 37 $-0.4$ $-0.0$ SLAC $ed$ $F_{2}$ [137,
138] 22.8 / 38 $-0.2$ $-0.8$ E665 $\mu p$ $F_{2}$ [139] 63.9 / 53 $+4.2$
$+4.5$ E665 $\mu d$ $F_{2}$ [139] 67.5 / 53 $+2.9$ $+2.5$ NuTeV $\nu N$
$F_{2}$ [140] 35.7 / 53 $-2.6$ $-1.3$ NuTeV $\nu N$ $xF_{3}$ [140] 34.8 / 42
$+4.1$ $+2.1$ NMC $\mu n/\mu p$ [141] 131.6 / 148 $+0.8$ $+2.2$ HERA $ep$
$F_{2}^{\text{charm}}$ [142] 143.7 / 79 $+11.4$ $+13.8$ NMC/BCDMS/SLAC/HERA
$F_{L}$ [135, 136, 138, 143, 144, 145] 45.6 / 57 $-22.9$ $-23.2$ CHORUS $\nu
N$ $F_{2}$ [146] 29.2 / 42 $-1.0$ $-0.8$ CHORUS $\nu N$ $xF_{3}$ [146] 18.1 /
28 $-0.3$ $-0.5$ HERA $e^{+}p$ CC [147] 49.7 / 39 $-2.3$ $-1.0$ HERA $e^{-}p$
CC [147] 64.9 / 42 $-5.3$ $-4.9$ HERA $e^{+}p$ NC $820\ \text{GeV}$ [147] 84.3
/ 75 $-5.6$ $-5.1$ HERA $e^{-}p$ NC $460\ \text{GeV}$ [147] 247.7 / 209 $-0.6$
$-0.7$ HERA $e^{+}p$ NC $920\ \text{GeV}$ [147] 474.0 / 402 $-38.7$ $-36.4$
HERA $e^{-}p$ NC $575\ \text{GeV}$ [147] 248.5 / 259 $-14.5$ $-14.1$ HERA
$e^{-}p$ NC $920\ \text{GeV}$ [147] 243.0 / 159 $-1.4$ $-2.1$ Total 2580.9 /
2375 $-90.8$ $-86.2$

Table 11: Table showing the relevant DIS datasets and how the individual
$\chi^{2}$ changes from NNLO by including the N3LO contributions to the
structure function $F_{2}(x,Q^{2})$. The result within purely NNLO $K$-factors
included for all data in the fit is also given.

To complement the discussions in Section 7, we isolate the $\chi^{2}$ results
from DIS data in Table 11. This data is directly affected by the N3LO
structure functions constructed approximately in Section’s 3 to 6. A
substantial decrease in the total $\chi^{2}$ from NNLO is observed across DIS
datasets. Considering the results in Table 11 in the context of Table’s 2 to
7, a better fit quality is observed for all DIS and non-DIS datasets than at
NNLO with the inclusion of N3LO contributions. As the DIS data makes up over
half of the total data included in a global fit, it is the dominant force in
deciding the overall form of the PDFs, especially at small-$x$ (discussed
further in Section 8.4). Table 11 further reinforces the point that the N3LO
description is flexible enough to fit to HERA and non-HERA data, without being
largely constrained by tensions between the small-$x$ (HERA dominated) and
large-$x$ (non-HERA dominated) regions.

###  8.2 Correlation Results

Figure 28: Correlation matrix for all N3LO theory parameters included in the
fit against the subset of the MSHT20 parameters (shown in black) used in
constructing the Hessian eigenvectors. This is shown for the case where the
$K$-factors correlations with the first 42 parameters are included. N3LO
theory parameters associated with the splitting functions are coloured blue,
the parameters affecting the transition matrix elements and coefficient
functions are in red and the $K$-factor parameters are in green.

The correlation matrix shown in Fig. 28 illustrates the correlations between
extra N3LO theory parameters and the subset of the MSHT20 parameters which are
included in the construction of Hessian eigenvectors (see Section 8.3 and [3]
for details). It is apparent that the correlations between $K$-factor
parameters for each process (shown in green) and other PDF and theory
parameters are usually small, with some exceptions e.g. for the
Top${}_{\text{NLO}}$ parameter. Due to this there is an argument that each
process’ $K$-factor parameters could be treated separately from all other
parameters in the Hessian prescription (see Section 4.1) which allows for a
more flexible PDF set that can be decorrelated from a process. By using the
uncorrelated Hessian results for a process NNLO hard cross sections can be
transformed to aN3LO and therefore provide more reliable predictions (more
details in Section 10). This is a fairly intuitive result, since most
correlations are showing a natural separation between the process dependent
and process independent physics in the DIS picture 181818The same pattern can
be seen for $c_{\\{q,g\\}}^{\mathrm{NLL}}$ parameters which are involved in
the DIS hard cross section.. Mathematically, the $K$-factors are directly
associated with the hard cross section, whereas other N3LO theory parameters
($\rho_{ij}$ and $a_{ij}$) are having a direct effect on the PDFs. Fig. 28
therefore begins to motivate the inclusion of the ‘pure’ theory (splitting
functions and transition matrix elements) parameters within the standard MSHT
eigenvector analysis [3], with the decorrelation of the $K$-factor parameters,
as discussed in Section 2.3. We investigate and compare both treatments
(complete correlation and $K$-factor decorrelation) throughout the rest of
this section. We show in Section 8.4 that while the decorrelation of
$K$-factors is not complete, both treatments result in similar uncertainty
bands, therefore confirming that the effect of making the assumption of full
decorrelation is minimal in practice. Note that although the
$c_{i}^{\mathrm{NLL}}$ parameters also show minimal correlation with other
parameters, we include these within the ‘pure’ theory group of parameters
(i.e. correlated with $\rho_{ij}$ and $a_{ij}$) as they are essential
ingredients in the underlying DIS theory.

###  8.3 Eigenvector Results

In the MSHT fitting procedure (described in [3]) the eigenvectors of a Hessian
matrix are found, which encapsulate the sources of uncertainties and
corresponding correlations. Combining these with the central PDFs, forms the
entire PDF set with uncertainties. In this eigenvector analysis a dynamical
rescaling of each eigenvector $e_{i}$ is performed via a tolerance factor $t$
to encapsulate the $68\%$ confidence limit (C.L.).

$a_{i}=a_{i}^{0}\pm te_{i},$ (8.1)

where $a_{i}^{0}$ is the best fit parameter. $t$ is then adjusted to give the
desired tolerance $T$ for the required confidence interval defined as
$T=\sqrt{\Delta\chi_{\mathrm{global}}^{2}}$ (for $68\%$ C.L.). In a quadratic
approximation, for suitably well-behaved eigenvectors, $t=T$ is true. Although
for eigenvectors with larger eigenvalues, it is possible to observe
significant deviations from $t=T$. The standard MSHT fitting procedure
involves allowing all relevant parameters from [3] to vary when finding the
best fit, now including all N3LO theory parameters
($\rho_{ij},a_{ij},c_{i}^{\mathrm{NLL}},K_{\mathrm{NLO/NNLO}}$) discussed in
this work. After accounting for high degrees of correlation between parameters
(described in [42]), the result is a Hessian matrix which in general, depends
on a subset of the parameters that were allowed to vary in a best fit and
provides a set of suitably well-behaved eigenvectors. The standard MSHT NNLO
PDF eigenvectors are based on a set of 32 parameters, reduced from the 52
parameters allowed to vary in the full fit. In the following analysis we are
therefore concerned with a smaller number of parameters, specifically the 32
parameters from the standard MSHT fitting procedure plus an extra 20 N3LO
parameters (shown in Fig. 28).

A standard choice of tolerance $T$ is
$T=\sqrt{\Delta\chi^{2}_{\mathrm{global}}}=1$ for a 68% C.L. limit. However,
this assumes all datasets are consistent with Gaussian errors. In practice,
due to incomplete theory, tensions between datasets and parameterisation
inflexibility, this is known not to be the case in a global PDF fit. To
overcome this, a 68% C.L. region for each dataset is defined. Then for each
eigenvector, the value of $\sqrt{\chi^{2}_{\mathrm{global}}}$ for each chosen
$t$ is recorded (ideally showing a quadratic behaviour). Finally, a value of
$T$ is chosen to ensure that all datasets are described within their 68% CL in
each eigenvector direction. For a fuller mathematical description of the
dynamical tolerance procedure used in MSHT PDF fits, the reader is referred to
[42]. In this section we present a demonstration of how well the resultant
eigenvectors follow the quadratic assumption based on $t=T$, including the
specific choices of dynamical tolerances and which dataset/penalty constrains
this tolerance in each eigenvector direction.

#### PDF + N3LO DIS Theory + N3LO $K$-factor (decorrelated) Parameters

As discussed in Section 8.2, when determining the eigenvectors and therefore
PDF uncertainties, we can choose to either include the correlations between
the 10 $K$-factor parameters added with the other 42 parameters (encompassing
the standard 32 MSHT eigenvector parameters and the 10 new theory parameters
from the splitting functions, transition matrix elements and coefficient
functions) or to decorrelate the 10 $K$-factor parameters.

In this section we address the scenario where we decorrelate the $K$-factors
as

$H_{ij}^{-1}+\sum_{p=1}^{N_{p}}K_{ij,p}^{-1}$ (8.2)

and consider each term individually.

Figure 29:  Correlation matrix of the first 42 (total 52) eigenvectors found
with the N3LO parameters added into the analysis in the case where the
$K$-factors are decorrelated from these first 42 parameters. Parameters
associated with the splitting functions are coloured blue, those affecting the
transition matrix elements and coefficient functions are in red. Figure 30:
Map of the 10 $K$-factor eigenvectors found with the N3LO parameters added
into the analysis in the case where the $K$-factors are decorrelated from
these first 42 parameters. Combined with the 42 eigenvectors shown in Fig. 29,
these form the total 52 eigenvectors in the decorrelated case. Parameters
associated with the $K$-factor parameters are in green.

Fig. 29 shows the map of eigenvectors produced from $H_{ij}$ in Equation
(8.2), where we have included the new N3LO DIS theory parameters (splitting
functions in blue and coefficient functions/transition matrix elements in red)
correlated with the PDF parameters. Eigenvectors 35 and 36 are prime examples
of where the eigenvectors have specifically encompassed the correlation/anti-
correlation between the two NLL FFNS coefficient function parameters
$c_{i}^{\mathrm{NLL}}$ ($i\in\\{q,g\\}$). Whereas the splitting functions
naturally give rise to a much more complicated mixing with other PDF
parameters as these directly affect the evolution of the PDFs. Due to the
direct impact of $\rho_{ij}$’s on the PDFs (via DGLAP evolution), combined
with the large contributions to the evolution shown at N3LO, this result is as
expected.

Another somewhat pleasing aspect is the recovery of a natural separation
between eigenvectors associated with the N3LO coefficient function/transition
matrix elements and our original PDF parameters (incl. N3LO splitting
functions). This separation is reminiscent of our DIS picture, whereby the
splitting functions are much more intertwined with the raw PDFs and the
transition matrix elements have a symbiotic relationship with the coefficient
functions (see GM-VFNS description in Section 3). Due to this, the form of
these eigenvectors has not only some level of physical interpretation
inherited from our underlying theory, but also offers a useful way to access
the different sources of N3LO additions within the PDF set.

In Fig. 30 the eigenvectors resulting from the
$\sum_{p=1}^{N_{p}}K_{ij,p}^{-1}$ terms in Equation (8.2) are shown. These
eigenvectors are constructed in pairs, describing the correlation and anti-
correlation of the two $K$-factor parameters (controlling the NLO and NNLO
contributions to N3LO) for each process $p$ contained within the corresponding
$K_{ij,p}$ correlation matrix.

Table 12 shows further information regarding the $K$-factor parameter limits
from each eigenvector. In most cases the parameter limits are well within the
allowed variation ($-1<a<1$), which is an indication that the data included in
the fit is constraining these parameters rather than the individual penalties
for each parameter191919We remind the reader that the Dimuon datasets also
include a branching ratio factor which is providing some compensation with
these $K$-factor parameters (as discussed in Section 7)..

To assess whether the eigenvectors are violating the quadratic treatment, four
examples displaying this behaviour are shown in Fig. 31, with a full analysis
provided in Appendix C.1. Additionally, Table 13 provides a summary of all
tolerances found within the eigenvector scans.

Matrix | Central Values | Eigenvector | \+ Limit | \- Limit | Scale
---|---|---|---|---|---
$a_{\mathrm{NLO}}$ | $a_{\mathrm{NNLO}}$ | $a_{\mathrm{NLO}}$ | $a_{\mathrm{NNLO}}$ | $a_{\mathrm{NLO}}$ | $a_{\mathrm{NNLO}}$
$K_{ij}^{\mathrm{DY}}$ | -0.282 | 0.079 | 43 | -0.378 | 0.062 | -0.145 | 0.103 | $m/2$
44 | -0.334 | 0.374 | -0.256 | -0.071
$K_{ij}^{\mathrm{Top}}$ | 0.041 | 0.651 | 45 | -0.564 | 0.455 | 0.692 | 0.862 | Section 7.2
46 | 0.026 | 1.210 | 0.070 | -0.456
$K_{ij}^{\mathrm{Jets}}$ | -0.300 | -0.691 | 47 | -0.515 | -0.957 | 0.105 | -0.189 | $p_{T}^{\mathrm{jet}}$
48 | -0.725 | -0.033 | 0.036 | -1.212
$K_{ij}^{p_{T}\ \mathrm{Jets}}$ | 0.583 | -0.080 | 49 | 0.388 | -0.406 | 0.812 | 0.301 | $p_{T}$
50 | 0.480 | 0.624 | 0.680 | -0.742
$K_{ij}^{\mathrm{Dimuon}}$ | -0.444 | 0.922 | 51 | -1.109 | -0.208 | -0.103 | 1.502 | $Q^{2}$
52 | -1.091 | 1.359 | 0.981 | -0.039
Table 12: Limiting values for specific $K$-factor parameters for each of the
processes considered in the decorrelated case. Parameter values are shown in
the positive and negative limits for each eigenvector. The scale choices for
top quark processes are described in Section 7.2 to be $H_{T}/4$ for the
single differential datasets with the exception of data differential in the
average transverse momentum of the top or antitop,
$p_{T}^{t},p_{T}^{\bar{t}}$, for which $m_{T}/2$ is used. For the double diff.
dataset the scale choice is $H_{T}/4$ and for the inclusive top
$\sigma_{t\bar{t}}$ a scale of $m_{t}$ is chosen.

Figure 31: Dynamic tolerance behaviour for 4 selected eigenvectors in the case
of decorrelated $K$-factor parameters. The black dots show the fixed tolerance
relations found for integer values of $t$, whereas the red triangles show the
final chosen dynamical tolerances for each eigenvector direction. For an
exhaustive analysis of all eigenvectors see Fig. 57.

# $t+$ $T+$ Limiting Factor ($+$) $t-$ $T-$ Limiting Factor ($-$) Primary
Parameter 1 3.88 3.86 ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] 2.94 3.13
ATLAS $8\ \text{TeV}$ double diff. $Z$ [100] $a_{S,3}$ 2 5.35 5.56 HERA
$e^{+}p$ NC $920\ \text{GeV}$ [147] 3.10 2.98 NMC $\mu d$ $F_{2}$ [136]
$a_{S,6}$ 3 4.20 4.47 NuTeV $\nu N\rightarrow\mu\mu X$ [133] 2.30 2.17 NMC
$\mu d$ $F_{2}$ [136] $a_{S,2}$ 4 4.87 4.85 ATLAS $8\ \text{TeV}$ double diff.
$Z$ [100] 1.81 1.80 NuTeV $\nu N\rightarrow\mu\mu X$ [133] $a_{s+,5}$ 5 3.30
3.63 ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] 2.95 2.77 NuTeV $\nu
N\rightarrow\mu\mu X$ [133] $a_{s+,3}$ 6 4.89 5.38 ATLAS $7\ \text{TeV}$ high
prec. $W,Z$ [79] 5.47 5.27 NMC $\mu d$ $F_{2}$ [136] $\delta_{g^{\prime}}$ 7
5.30 5.47 DØ II $W\rightarrow\nu e$ asym. [86] 3.44 3.41 DØ $W$ asym. [97]
$\delta_{S}$ 8 4.12 4.05 DØ II $W\rightarrow\nu e$ asym. [86] 2.19 2.64 DØ $W$
asym. [97] $\delta_{S}$ 9 3.88 3.93 BCDMS $\mu p$ $F_{2}$ [135] 6.78 7.11 CMS
W asym. $p_{T}>25,30\ \text{GeV}$ [89] $\delta_{d}$ 10 3.24 3.44 BCDMS $\mu p$
$F_{2}$ [135] 5.64 5.83 DØ II $W\rightarrow\nu\mu$ asym. [84] $a_{u,6}$ 11
2.22 4.41 E866 / NuSea $pp$ DY [80] 5.31 4.79 ATLAS $7\ \text{TeV}$ high prec.
$W,Z$ [79] $a_{s+,2}$ 12 2.62 2.81 DØ $W$ asym. [97] 5.49 5.71 ATLAS $8\
\text{TeV}\ W$ [99] $\delta_{d}$ 13 1.12 3.37 E866 / NuSea $pp$ DY [80] 3.87
3.19 DØ $W$ asym. [97] $a_{\rho,3}$ 14 2.01 2.34 NuTeV $\nu N\rightarrow\mu\mu
X$ [133] 3.54 3.60 NuTeV $\nu N$ $xF_{3}$ [140] $A_{s-}$ 15 3.69 3.75 ATLAS
$8\ \text{TeV}\ Z\ p_{T}$ [110] 3.65 4.05 NMC $\mu d$ $F_{2}$ [136] $a_{g,3}$
16 4.42 5.08 ATLAS $8\ \text{TeV}\ W$ [99] 4.96 4.84 ATLAS $8\ \text{TeV}\ Z\
p_{T}$ [110] $a_{g,3}$ 17 1.23 1.02 $\rho_{qq}^{NS}$ 1.24 1.05
$\rho_{qq}^{NS}$ $\rho_{qq}^{NS}$ 18 3.09 3.50 DØ $W$ asym. [97] 2.93 3.06 DØ
$W$ asym. [97] $a_{d,6}$ 19 4.69 4.57 BCDMS $\mu p$ $F_{2}$ [135] 3.85 5.06
CMS $7\ \text{TeV}$ jets [101] $\eta_{u}$ 20 5.40 5.37 NuTeV $\nu
N\rightarrow\mu\mu X$ [133] 5.04 5.60 HERA $e^{+}p$ NC $920\ \text{GeV}$ [147]
$\delta_{g^{\prime}}$ 21 2.10 2.36 E866 / NuSea $pd/pp$ DY [81] 1.26 1.69 E866
/ NuSea $pd/pp$ DY [81] $a_{\rho,6}$ 22 3.04 3.27 NuTeV $\nu
N\rightarrow\mu\mu X$ [133] 1.90 2.27 DØ $W$ asym. [97] $a_{d,3}$

Table 13: Tolerances resulting from eigenvector scans with decorrelated
$K$-factors for each process. The average tolerance for this set of
eigenvectors is $T=3.34$.

# $t+$ $T+$ Limiting Factor ($+$) $t-$ $T-$ Limiting Factor ($-$) Primary
Parameter 23 5.61 6.20 HERA $e^{+}p$ NC $920\ \text{GeV}$ [147] 5.43 6.13 HERA
$e^{+}p$ NC $920\ \text{GeV}$ [147] $\rho_{qq}^{PS}$ 24 3.53 3.70 E866 / NuSea
$pd/pp$ DY [81] 1.47 1.82 DØ $W$ asym. [97] $a_{d,2}$ 25 1.60 2.03 E866 /
NuSea $pd/pp$ DY [81] 4.37 4.88 CMS W asym. $p_{T}>35\ \text{GeV}$ [88]
$A_{\rho}$ 26 1.00 1.46 E866 / NuSea $pd/pp$ DY [81] 3.02 3.00 DØ $W$ asym.
[97] $a_{\rho,1}$ 27 1.60 2.09 $\rho_{gq}$ 4.16 5.53 ATLAS $7\ \text{TeV}$
high prec. $W,Z$ [79] $\rho_{gq}$ 28 1.51 2.16 $\rho_{gq}$ 3.20 4.07 ATLAS $8\
\text{TeV}$ sing. diff. $t\bar{t}$ dilep. [131] $\eta_{g}$ 29 2.99 3.21 CMS
$8\ \text{TeV}\ W$ [96] 2.04 2.52 NuTeV $\nu N$ $xF_{3}$ [140] $\eta_{S}$ 30
0.97 1.30 DØ $W$ asym. [97] 3.56 4.90 ATLAS $8\ \text{TeV}\ W$ [99]
$\eta_{d}-\eta_{u}$ 31 4.78 6.64 HERA $e^{+}p$ NC $920\ \text{GeV}$ [147] 1.77
1.81 $\rho_{gq}$ $\rho_{qg}$ 32 2.51 7.32 BCDMS $\mu p$ $F_{2}$ [135] 1.80
3.96 DØ $W$ asym. [97] $\eta_{d}-\eta_{u}$ 33 2.71 3.35 CMS $7\ \text{TeV}\
W+c$ [62] 2.92 3.37 NuTeV $\nu N\rightarrow\mu\mu X$ [133] $\eta_{s+}$ 34 3.53
3.91 HERA $ep$ $F_{2}^{\text{charm}}$ [142] 3.87 4.82 CMS $7\ \text{TeV}\ W+c$
[62] $\eta_{s+}$ 35 3.67 3.50 NuTeV $\nu N\rightarrow\mu\mu X$ [133] 2.72 3.62
NuTeV $\nu N\rightarrow\mu\mu X$ [133] $\eta_{s-}$ 36 1.25 1.26
$c_{q}^{\mathrm{NLL}}$ 1.41 1.57 $c_{q}^{\mathrm{NLL}}$ $c_{q}^{\mathrm{NLL}}$
37 0.71 0.77 $\rho_{gg}$ 4.25 5.33 NuTeV $\nu N\rightarrow\mu\mu X$ [133]
$\rho_{gg}$ 38 2.72 2.71 CMS $8\ \text{TeV}\ W$ [96] 3.60 4.21 $\rho_{gq}$
$A_{s+}$ 39 1.56 5.03 ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] 1.60 5.82
ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] $A_{S}$ 40 0.97 1.00
$a_{qq,H}^{\mathrm{NS}}$ 1.04 1.01 $a_{qq,H}^{\mathrm{NS}}$
$a_{qq,H}^{\mathrm{NS}}$ 41 1.56 2.87 HERA $ep$ $F_{2}^{\text{charm}}$ [142]
1.93 2.34 $\rho_{gg}$ $a_{Hg}$ 42 0.63 0.53 $a_{gg,H}$ 2.46 2.72 $a_{gg,H}$
$a_{gg,H}$ 43 1.78 2.30 ATLAS $8\ \text{TeV}$ double diff. $Z$ [100] 2.56 1.94
CMS double diff. Drell-Yan [94] $\mathrm{DY}_{\mathrm{NLO}}$ 44 4.63 4.88
ATLAS $8\ \text{TeV}\ W$ [99] 2.35 2.09 ATLAS $8\ \text{TeV}$ double diff. $Z$
[100] $\mathrm{DY}_{\mathrm{NNLO}}$ 45 5.05 1.46 ATLAS $8\ \text{TeV}$ sing.
diff. $t\bar{t}$ dilep. [131] 5.44 1.46 Tevatron, ATLAS, CMS
$\sigma_{t\bar{t}}$ [117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127,
128, 129] $\mathrm{Top}_{\mathrm{NLO}}$

Table 14: (Continued) Tolerances resulting from eigenvector scans with
decorrelated $K$-factors for each process. The average tolerance for this set
of eigenvectors is $T=3.34$.

# $t+$ $T+$ Limiting Factor ($+$) $t-$ $T-$ Limiting Factor ($-$) Primary
Parameter 46 1.34 1.31 $\mathrm{Top}_{\mathrm{NLO}}$ 2.65 2.68 ATLAS $8\
\text{TeV}$ sing. diff. $t\bar{t}$ [130] $\mathrm{Top}_{\mathrm{NNLO}}$ 47
1.47 1.57 CDF II $p\bar{p}$ incl. jets [102] 2.78 2.96 $a_{gg,H}$
$\mathrm{Jet}_{\mathrm{NLO}}$ 48 3.26 3.38 CMS $2.76\ \text{TeV}$ jet [106]
2.58 3.01 CMS $7\ \text{TeV}$ jets [101] $\mathrm{Jet}_{\mathrm{NNLO}}$ 49
2.20 2.39 ATLAS $8\ \text{TeV}\ Z\ p_{T}$ [110] 2.58 2.56 ATLAS $8\
\text{TeV}\ W+\text{jets}$ [111] $p_{T}\ \mathrm{Jets}_{\mathrm{NLO}}$ 50 2.36
2.43 ATLAS $8\ \text{TeV}\ Z\ p_{T}$ [110] 2.22 2.23 ATLAS $8\ \text{TeV}\
W+\text{jets}$ [111] $p_{T}\ \mathrm{Jets}_{\mathrm{NNLO}}$ 51 1.38 1.55
$\mathrm{Dimuon}_{\mathrm{NLO}}$ 0.71 0.81 $\mathrm{Dimuon}_{\mathrm{NNLO}}$
$\mathrm{Dimuon}_{\mathrm{NNLO}}$ 52 0.85 0.85
$\mathrm{Dimuon}_{\mathrm{NLO}}$ 1.86 1.87 $\mathrm{Dimuon}_{\mathrm{NLO}}$
$\mathrm{Dimuon}_{\mathrm{NLO}}$

Table 15: (Continued) Tolerances resulting from eigenvector scans with
decorrelated $K$-factors for each process. The average tolerance for this set
of eigenvectors is $T=3.34$.

There is relatively consistent agreement between $t$ and $T$ across all
eigenvectors with later eigenvectors (i.e. higher #) generally becoming less
quadratic (a feature which is built into the fit). Eigenvectors 31, 41 and 42
displayed in Fig. 31 are shown in Table 13 to be either dominated or limited
by at least one new N3LO parameter. Conversely, eigenvector 26 is much more
dominated by the original PDF parameters from MSHT20 NNLO. Comparing these
cases, the eigenvectors associated more strongly with the N3LO parameters
exhibit a similar level of agreement (and occasionally better) with the
desired quadratic behaviour as eigenvectors more closely associated with the
original PDF parameters.

The last 5 sets of eigenvectors (i.e. the last 10 where a set contains 2
eigenvectors for a particular process) we see in Table 13 are the decorrelated
$K$-factor eigenvectors, where there are correlated/anti-correlated
eigenvectors for each process. For all $K$-factor cases, Table 13 provides
sensible results with either the dominant datasets or parameter penalties
constraining each eigenvector direction. One interesting feature one can
observe here is a sign of tension between the ATLAS $8\ \text{TeV}\ Z\ p_{T}$
[110] and ATLAS $8\ \text{TeV}\ W+\text{jets}$ [111] datasets where the
limiting factors in Table 13 for eigenvectors 49 and 50 show that these
datasets are preferring a slightly different $K$-factor.

To provide some extra level of comparison between the eigenvectors shown here
and the eigenvectors found in the NNLO case, the average tolerance $T$ for
aN3LO (decorrelated $K$-factors) set is 3.34, compared to the NNLO average $T$
of 3.37.

#### PDF + N3LO DIS Theory + N3LO $K$-factor (correlated) Parameters

In this section we address the scenario,

$H_{ij}^{\prime}=\left(H_{ij}^{-1}+\sum_{p=1}^{N_{p}}K_{ij,p}^{-1}\right)^{-1}.$
(8.3)

Moving to an analysis including aN3LO $K$-factors as correlated parameters
with PDF and other N3LO theory parameters. This provides a comparison to the
case of decorrelated $K$-factors and justification for treating the cross
section behaviour separately to the PDF theory behaviour.

Fig. 32 shows a map of eigenvectors with the extra 10 N3LO $K$-factor
parameters (shown in green) included into the correlations considered. As
expected, the result of including the correlations between PDF parameters and
aN3LO $K$-factors results in a slightly more intertwined set of eigenvectors
(although a high level of decorrelation remains). Specifically, due to the
much higher number of DY datasets included in the global fit, these N3LO
$K$-factor parameters tend to be included across more of a spread of
eigenvectors. On the other hand, the Dimuon $K$-factors are almost entirely
isolated within two eigenvectors, similar to the decorrelated case.

Figure 32: Map of eigenvectors found with the N3LO theory and $K$-factor
parameters added into the analysis. Parameters associated with the splitting
functions are coloured blue, those affecting the transition matrix elements
and coefficient functions are in red and the $K$-factor parameters are in
green.

Figure 33: Dynamic tolerance behaviour for 4 selected eigenvectors in the case
of correlated $K$-factor parameters. The black dots show the fixed tolerance
relations found for integer values of $t$, whereas the red triangles show the
final chosen dynamical tolerances for each eigenvector direction. For an
exhaustive analysis of all eigenvectors see Fig. 64.

# $t+$ $T+$ Limiting Factor ($+$) $t-$ $T-$ Limiting Factor ($-$) Primary
Parameter 1 3.76 3.67 ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] 2.95 3.06
ATLAS $8\ \text{TeV}$ double diff. $Z$ [100] $a_{S,3}$ 2 3.67 3.52 NMC $\mu d$
$F_{2}$ [136] 5.37 5.54 HERA $e^{+}p$ NC $920\ \text{GeV}$ [147] $a_{S,6}$ 3
4.15 4.33 NuTeV $\nu N\rightarrow\mu\mu X$ [133] 2.58 2.39 NMC $\mu d$ $F_{2}$
[136] $a_{S,2}$ 4 4.90 4.82 ATLAS $8\ \text{TeV}$ double diff. $Z$ [100] 2.21
2.16 NuTeV $\nu N\rightarrow\mu\mu X$ [133] $a_{s+,5}$ 5 2.89 2.61 NuTeV $\nu
N\rightarrow\mu\mu X$ [133] 3.34 3.55 CCFR $\nu N\rightarrow\mu\mu X$ [133]
$a_{s+,3}$ 6 5.58 5.23 NMC $\mu d$ $F_{2}$ [136] 5.49 5.84 ATLAS $7\
\text{TeV}$ high prec. $W,Z$ [79] $\delta_{g^{\prime}}$ 7 5.34 5.44 DØ II
$W\rightarrow\nu e$ asym. [86] 3.46 3.37 DØ $W$ asym. [97] $\delta_{S}$ 8 4.11
3.82 DØ II $W\rightarrow\nu e$ asym. [86] 2.18 2.48 DØ $W$ asym. [97]
$\delta_{S}$ 9 6.39 6.59 CMS W asym. $p_{T}>25,30\ \text{GeV}$ [89] 4.20 4.07
BCDMS $\mu p$ $F_{2}$ [135] $\delta_{d}$ 10 3.14 3.11 BCDMS $\mu p$ $F_{2}$
[135] 5.66 5.66 DØ II $W\rightarrow\nu\mu$ asym. [84] $a_{u,6}$ 11 2.07 3.99
E866 / NuSea $pp$ DY [80] 5.88 5.00 ATLAS $7\ \text{TeV}$ high prec. $W,Z$
[79] $a_{g,2}$ 12 2.76 2.76 DØ $W$ asym. [97] 5.60 5.67 CMS $8\ \text{TeV}\ W$
[96] $\delta_{d}$ 13 3.73 3.01 DØ $W$ asym. [97] 1.26 3.96 E866 / NuSea $pp$
DY [80] $a_{\rho,3}$ 14 3.85 3.55 NuTeV $\nu N$ $xF_{3}$ [140] 2.24 2.30 NuTeV
$\nu N\rightarrow\mu\mu X$ [133] $A_{s-}$ 15 3.71 3.51 ATLAS $8\ \text{TeV}\
Z\ p_{T}$ [110] 3.44 3.72 NuTeV $\nu N\rightarrow\mu\mu X$ [133] $a_{g,3}$ 16
1.78 1.18 $\rho_{qq}^{NS}$ 1.78 1.24 $\rho_{qq}^{NS}$ $\rho_{qq}^{NS}$ 17 4.76
4.49 DØ $W$ asym. [97] 4.37 4.81 ATLAS $8\ \text{TeV}\ W$ [99] $a_{g,3}$ 18
3.44 3.81 DØ $W$ asym. [97] 3.18 2.94 DØ $W$ asym. [97] $a_{d,6}$ 19 5.33 5.02
DØ $W$ asym. [97] 3.77 4.42 CMS $7\ \text{TeV}$ jets [101] $\eta_{u}$ 20 6.04
5.82 HERA $e^{+}p$ NC $920\ \text{GeV}$ [147] 5.08 5.63 HERA $e^{+}p$ NC $920\
\text{GeV}$ [147] $\delta_{g^{\prime}}$ 21 2.14 2.02 E866 / NuSea $pd/pp$ DY
[81] 1.29 1.39 E866 / NuSea $pd/pp$ DY [81] $a_{\rho,6}$ 22 1.95 2.03 DØ $W$
asym. [97] 2.88 2.83 NuTeV $\nu N\rightarrow\mu\mu X$ [133] $a_{d,3}$

Table 16: Tolerances resulting from eigenvector scans with correlated
$K$-factors for each process. The average tolerance for this set of
eigenvectors is $T=3.57$.

# $t+$ $T+$ Limiting Factor ($+$) $t-$ $T-$ Limiting Factor ($-$) Primary
Parameter 23 5.33 5.87 HERA $e^{+}p$ NC $920\ \text{GeV}$ [147] 5.66 6.18 HERA
$e^{+}p$ NC $920\ \text{GeV}$ [147] $\rho_{qq}^{PS}$ 24 1.47 1.52 DØ $W$ asym.
[97] 3.51 3.48 E866 / NuSea $pd/pp$ DY [81] $a_{d,2}$ 25 4.34 4.71 CMS W asym.
$p_{T}>35\ \text{GeV}$ [88] 1.64 1.59 E866 / NuSea $pd/pp$ DY [81] $A_{\rho}$
26 1.00 1.24 E866 / NuSea $pd/pp$ DY [81] 3.12 2.84 DØ $W$ asym. [97]
$a_{\rho,1}$ 27 4.86 6.43 CMS double diff. Drell-Yan [94] 1.68 2.01
$\rho_{gq}$ $\rho_{gq}$ 28 2.36 2.66 ATLAS $7\ \text{TeV}$ high prec. $W,Z$
[79] 1.64 1.36 CMS double diff. Drell-Yan [94] $\mathrm{DY}_{\mathrm{NLO}}$ 29
1.95 1.41 $\rho_{qq}^{NS}$ 1.94 1.47 $\rho_{qq}^{NS}$
$\mathrm{Top}_{\mathrm{NLO}}$ 30 1.84 2.22 $\rho_{gq}$ 2.70 2.56 ATLAS $8\
\text{TeV}$ sing. diff. $t\bar{t}$ dilep. [131] $\eta_{g}$ 31 3.53 3.52 ATLAS
$8\ \text{TeV}\ Z\ p_{T}$ [110] 2.60 2.73 NuTeV $\nu N$ $xF_{3}$ [140] $p_{T}\
\mathrm{Jets}_{\mathrm{NLO}}$ 32 2.85 3.15 NuTeV $\nu N$ $xF_{3}$ [140] 3.84
3.57 ATLAS $8\ \text{TeV}\ W+\text{jets}$ [111] $p_{T}\
\mathrm{Jets}_{\mathrm{NLO}}$ 33 5.40 6.16 ATLAS $7\ \text{TeV}$ high prec.
$W,Z$ [79] 1.63 1.68 DØ $W$ asym. [97] $\mathrm{DY}_{\mathrm{NNLO}}$ 34 3.94
4.19 ATLAS $8\ \text{TeV}\ W$ [99] 1.56 1.86 DØ $W$ asym. [97]
$\mathrm{DY}_{\mathrm{NNLO}}$ 35 5.06 6.73 HERA $e^{+}p$ NC $920\ \text{GeV}$
[147] 2.00 1.84 $\rho_{gq}$ $\rho_{qg}$ 36 2.31 2.48 CDF II $p\bar{p}$ incl.
jets [102] 4.73 4.88 ATLAS $7\ \text{TeV}$ jets [104]
$\mathrm{Jet}_{\mathrm{NLO}}$ 37 2.72 7.05 BCDMS $\mu p$ $F_{2}$ [135] 1.67
3.04 DØ $W$ asym. [97] $\eta_{d}-\eta_{u}$ 38 2.63 2.90
$\mathrm{Jet}_{\mathrm{NLO}}$ 3.24 3.35 CMS $2.76\ \text{TeV}$ jet [106]
$\mathrm{Jet}_{\mathrm{NNLO}}$ 39 2.41 2.34 ATLAS $8\ \text{TeV}\
W+\text{jets}$ [111] 2.88 2.99 ATLAS $8\ \text{TeV}\ Z\ p_{T}$ [110] $p_{T}\
\mathrm{Jets}_{\mathrm{NNLO}}$ 40 2.98 3.44 CMS $7\ \text{TeV}\ W+c$ [62] 2.43
2.36 NuTeV $\nu N\rightarrow\mu\mu X$ [133] $\eta_{s+}$ 41 4.59 5.42 HERA $ep$
$F_{2}^{\text{charm}}$ [142] 2.93 3.20 HERA $ep$ $F_{2}^{\text{charm}}$ [142]
$c_{g}^{\mathrm{NLL}}$ 42 2.74 2.83 CMS $8\ \text{TeV}$ double diff.
$t\bar{t}$ [132] 1.36 1.28 $\mathrm{Top}_{\mathrm{NLO}}$
$\mathrm{Top}_{\mathrm{NNLO}}$ 43 0.83 0.71 $\mathrm{Dimuon}_{\mathrm{NLO}}$
1.36 1.46 $p_{T}\ \mathrm{Jets}_{\mathrm{NNLO}}$
$\mathrm{Dimuon}_{\mathrm{NNLO}}$ 44 0.81 0.80
$\mathrm{Dimuon}_{\mathrm{NLO}}$ 2.07 2.05 $p_{T}\
\mathrm{Jets}_{\mathrm{NNLO}}$ $\mathrm{Dimuon}_{\mathrm{NLO}}$ 45 1.45 1.46
$\mathrm{Dimuon}_{\mathrm{NLO}}$ 2.09 2.38 $p_{T}\
\mathrm{Jets}_{\mathrm{NNLO}}$ $\eta_{s-}$

Table 17: (Continued) Tolerances resulting from eigenvector scans with
correlated $K$-factors for each process. The average tolerance for this set of
eigenvectors is $T=3.57$.

# $t+$ $T+$ Limiting Factor ($+$) $t-$ $T-$ Limiting Factor ($-$) Primary
Parameter 46 1.41 1.57 $c_{q}^{\mathrm{NLL}}$ 1.25 1.26 $c_{q}^{\mathrm{NLL}}$
$c_{q}^{\mathrm{NLL}}$ 47 0.71 0.77 $\rho_{gg}$ 4.27 5.21 NuTeV $\nu
N\rightarrow\mu\mu X$ [133] $\rho_{gg}$ 48 3.62 4.07 $\rho_{gq}$ 2.72 2.63 CMS
$8\ \text{TeV}\ W$ [96] $A_{s+}$ 49 1.53 4.81 ATLAS $7\ \text{TeV}$ high prec.
$W,Z$ [79] 1.56 5.54 ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] $A_{S}$ 50
0.97 1.00 $a_{qq,H}^{\mathrm{NS}}$ 1.04 1.00 $a_{qq,H}^{\mathrm{NS}}$
$a_{qq,H}^{\mathrm{NS}}$ 51 1.56 2.69 HERA $ep$ $F_{2}^{\text{charm}}$ [142]
1.93 1.98 $\rho_{gg}$ $a_{Hg}$ 52 0.63 0.54 $a_{gg,H}$ 2.46 2.68 $a_{gg,H}$
$a_{gg,H}$

Table 18: (Continued) Tolerances resulting from eigenvector scans with
correlated $K$-factors for each process. The average tolerance for this set of
eigenvectors is $T=3.57$.

Once again, to investigate deviations from the quadratic behaviour, Fig. 33
illustrates examples of the tolerance behaviours of selected eigenvectors,
with a full analysis provided in Appendix C.2. Further to this, Table 16
displays the tolerances and limiting datasets/parameters for the 52 correlated
eigenvectors. It is difficult to compare and contrast these results with the
decorrelated case, since the eigenvectors are inherently different. However in
both cases, the eigenvectors are similarly well behaved, exhibit relatively
good consistency between $t$ and $T$ and are therefore providing valid
descriptions for a PDF fit.

For most of the 12 eigenvectors with N3LO $K$-factors as primary parameters,
there is expected behaviour, with the eigenvectors constrained either by their
own penalties or by dominant datasets for the associated process. However, due
to the extra correlations considered, there are a small number of eigenvector
directions which are not as trivial to explain (e.g. eigenvector 31). We
therefore recover the lack of correlation between $K$-factor parameters seen
within Fig. 28 in the set of correlated PDF eigenvectors presented here.
Further to this, comparing the $t$ and $T$ values found for eigenvectors
associated with N3LO $K$-factors in Table’s 13 and 16, one can observe clear
similarities between eigenvectors. This suggests that even when correlating
the $K$-factor parameters, the fit succeeds in decorrelating the individual
processes, thereby motivating our original assumption that the correlations
with $K$-factors can be ignored. Another similarity one can observe between
Table 13 and Table 16 is the suggestion of some tension between ATLAS $8\
\text{TeV}\ Z\ p_{T}$ [110] and ATLAS $8\ \text{TeV}\ W+\text{jets}$ [111]
datasets seen in the limiting factors of eigenvector 39 in the correlated
case.

Eigenvectors 27, 29 and 52 displayed in Fig. 33 can be seen from Table 16 to
be associated with the new N3LO theory parameters. Whereas eigenvector 37 is
primarily focused on an original PDF parameter. One can observe a similar
level of quadratic behaviour across all four of these eigenvector tolerances.
Comparing all eigenvectors in the decorrelated/correlated cases, the
behaviours are similarly well behaved. The average tolerance $T$ for the aN3LO
(with correlated $K$-factors) case is 3.57, slightly higher than the NNLO
average of 3.37 and the aN3LO (with decorrelated $K$-factors) average of 3.34.

###  8.4 PDF Results

Figure 34: General forms of NNLO (top) and aN3LO (bottom) PDFs at low (left)
and high (right) $Q^{2}$. Several main features can be compared and contrasted
such as the marked increase in the gluon and charm at small-$x$ (note the
difference in y-axis scale between NNLO (top) and aN3LO (bottom)).

Fig. 34 displays the overall shape of the PDFs including the N3LO additions
compared to the standard NNLO set. We provide this comparison to accompany the
results described in earlier sections. At small-$x$ and low-$Q^{2}$ the gluon
exhibits a marked enhancement due to the large small-$x$ logarithms inserted
at N3LO. The changes induced from specific N3LO contributions are investigated
in Section 8.8.

Figure 35: Low-$Q^{2}$ ratio plots showing the aN3LO 68% confidence intervals
with decorrelated ($H_{ij}+K_{ij}$) and correlated ($H_{ij}^{\prime}$)
$K$-factor parameters, compared to NNLO 68% confidence intervals. Also shown
are the central values at NNLO when fit to all non-HERA datasets which show
similarities with N3LO in the large-$x$ region of selected PDF flavours. All
plots are shown for $Q^{2}=10\ \mathrm{GeV}^{2}$ with the exception of the
bottom quark shown for $Q^{2}=25\ \mathrm{GeV}^{2}$.

Figure 36: (Continued) Low-$Q^{2}$ ratio plots showing the aN3LO 68%
confidence intervals with decorrelated ($H_{ij}+K_{ij}$) and correlated
($H_{ij}^{\prime}$) $K$-factor parameters, compared to NNLO 68% confidence
intervals. Also shown are the central values at NNLO when fit to all non-HERA
datasets which show similarities with N3LO in the large-$x$ region of selected
PDF flavours. All plots are shown for $Q^{2}=10\ \mathrm{GeV}^{2}$ with the
exception of the bottom quark shown for $Q^{2}=25\ \mathrm{GeV}^{2}$.

Figure 37: High-$Q^{2}$ ratio plots showing the aN3LO 68% confidence intervals
with decorrelated ($H_{ij}+K_{ij}$) and correlated ($H_{ij}^{\prime}$)
$K$-factor parameters, compared to NNLO 68% confidence intervals. Also shown
are the central values at NNLO when fit to all non-HERA datasets which show
similarities with N3LO in the large-$x$ region of selected PDF flavours. All
plots are shown for $Q^{2}=10^{4}\ \mathrm{GeV}^{2}$.

Figure 38: (Continued) High-$Q^{2}$ ratio plots showing the aN3LO 68%
confidence intervals with decorrelated ($H_{ij}+K_{ij}$) and correlated
($H_{ij}^{\prime}$) $K$-factor parameters, compared to NNLO 68% confidence
intervals. Also shown are the central values at NNLO when fit to all non-HERA
datasets which show similarities with N3LO in the large-$x$ region of selected
PDF flavours. All plots are shown for $Q^{2}=10^{4}\ \mathrm{GeV}^{2}$.

Shown in Fig.’s 35 and 37 are the ratios for each flavour of aN3LO PDF
compared to the NNLO set with their 68% confidence intervals at low and
high-$Q^{2}$ respectively. The shaded aN3LO regions indicate the PDF
uncertainty produced with the decorrelated
($(H^{-1}_{ij}+\sum_{p=1}^{N_{p}}K^{-1}_{ij,\ p})^{-1}$) aN3LO $K$-factors for
each process. As a comparison to these shaded regions, the bounds of
uncertainty for the fully correlated ($H_{ij}^{\prime}$) N3LO $K$-factor
parameters are also provided (red dashed line).

Considering Fig. 35 we present the aN3LO PDF set at $Q^{2}=10\
\mathrm{GeV}^{2}$ with the bottom quark PDF at $Q^{2}=25\ \mathrm{GeV}^{2}$.
These PDF ratios better display the substantial increase in the gluon at
small-$x$, reminiscent of the gluon PDF presented in [16, 17]. The predicted
harder small-$x$ gluon is then accommodated for by reductions in the PDFs at
large and small-$x$ (particularly the gluon near $x=10^{-2}$) from NNLO.
Another prominent feature is the enhanced charm and bottom quark at N3LO.
Since the heavy flavour quarks are perturbatively calculated in the MSHT
framework, this amplification is a feature of the transition matrix element
$A_{Hg}^{(3)}$ at high-$x$, combined with the increase in the gluon PDF at
small-$x$ (as these two ingredients are convoluted together). Comparing with
Fig. 95 in [3], we observe that the approximate N3LO charm quark now follows a
much closer trend to the CT18 PDF and is therefore even more significantly
different from the NNPDF NNLO fitted charm at large-$x$ than MSHT20 at NNLO.
In the high-$Q^{2}$ setting shown in Fig. 37 we observe similar albeit less
drastic effects to those described above.

Also contained in Fig.’s 35 and 37 are the relative forms of NNLO PDFs when
fit to all non-HERA data (full $\chi^{2}$ results are provided in Appendix B).
Comparing the non-HERA NNLO PDFs with aN3LO PDFs, there are some similarities
in the shapes and magnitudes of a handful of PDFs in the intermediate to
large-$x$ regime, most noticeably the light quarks. At small-$x$ the HERA data
heavily constrains the PDF fit and therefore these similarities rapidly break
down. However, this analysis displays further evidence that including N3LO
contributions, even though approximate, reduces tensions between the HERA and
non-HERA data (when considering the reduction in tension seen in Table B.1).
The aN3LO PDFs are seemingly able to fit to HERA and non-HERA datasets with
superior flexibility than at NNLO.

While in principle the negativity of quarks is possible in the $\overline{MS}$
scheme, it is unlikely to be correct at very high scales and the behaviour can
lead to issues concerning negative cross section predictions [148, 149]. In
the case of the $\overline{d}$, the form of this PDF has a negative central
value above $x\sim 0.5$ with a minimum of $\sim-0.001$ at $x\sim 0.6$. It is
also noted that although the $\overline{d}$ central value becomes negative in
this region, it is still positive within PDF uncertainties. These features are
not uncommon in PDF analyses and are discussed in detail in [8]. The proposed
smoothing of parameterisations employed in [8] ensures the definite positive
nature of PDFs in the high-$x$ region. Comparing the negativity of the
approximate N3LO $\overline{d}$ PDF with that in [8], the $\overline{d}$ PDF
presented here is much less negative and positive within PDF uncertainties.
Due to this and the fact that this effect is only apparent in the
$\overline{d}$, we present these PDFs as they are. We also note that in the
current MSHT20 fit, recent results surrounding the $\overline{d}/\overline{u}$
from the SeaQuest collaboration [150] are not included at the time of writing.
It is therefore only the E866 / NuSea $pd/pp$ DY dataset [81] that is
constraining this ratio, which is not as precise as the more recent results.
However, SeaQuest results suggest a preference for a higher $\overline{d}$ at
large-$x$, therefore including this data may in fact help constrain the
high-$x$ $\overline{d}$ behaviour seen here.

Figure 39: Low-$Q^{2}$ ratio plots showing the aN3LO 68% confidence intervals
with decorrelated and correlated $K$-factor parameters, compared to the aN3LO
central value. Also shown are the central values at aN3LO when fit to all non-
HERA datasets and the central values with all $K$-factors set at NNLO. All
plots are shown for $Q^{2}=10\ \mathrm{GeV}^{2}$ with the exception of the
bottom quark shown for $Q^{2}=25\ \mathrm{GeV}^{2}$.

Figure 40: (Continued) Low-$Q^{2}$ ratio plots showing the aN3LO 68%
confidence intervals with decorrelated and correlated $K$-factor parameters,
compared to the aN3LO central value. Also shown are the central values at
aN3LO when fit to all non-HERA datasets and the central values with all
$K$-factors set at NNLO. All plots are shown for $Q^{2}=10\ \mathrm{GeV}^{2}$
with the exception of the bottom quark shown for $Q^{2}=25\ \mathrm{GeV}^{2}$.

Figure 41:  High-$Q^{2}$ ratio plots showing the aN3LO 68% confidence
intervals with decorrelated and correlated $K$-factor parameters, compared to
the aN3LO central value. Also shown are the central values at aN3LO when fit
to all non-HERA datasets and the central values with all $K$-factors set at
NNLO. All plots are shown for $Q^{2}=10^{4}\ \mathrm{GeV}^{2}$.

Figure 42: (Continued) High-$Q^{2}$ ratio plots showing the aN3LO 68%
confidence intervals with decorrelated and correlated $K$-factor parameters,
compared to the aN3LO central value. Also shown are the central values at
aN3LO when fit to all non-HERA datasets and the central values with all
$K$-factors set at NNLO. All plots are shown for $Q^{2}=10^{4}\
\mathrm{GeV}^{2}$.

Fig.’s 39 and 41 express the aN3LO PDFs with decorrelated (green shaded
region) and correlated (red dashed lines) aN3LO $K$-factors at low and
high-$Q^{2}$ respectively (again with the bottom quark provided at $Q^{2}=25\
\mathrm{GeV}^{2}$ at low-$Q^{2}$) as a ratio to the N3LO central value. For
comparison we also include the level of uncertainty predicted with all N3LO
theory fixed (blue shaded region) i.e. only considering the variation without
N3LO theoretical uncertainty.

Comparing the two different aN3LO sets in Fig.’s 39 and 41, in general there
is good agreement between the total uncertainties considering the cases with
correlated (red dash) and decorrelated (green shaded) aN3LO $K$-factors. The
differences that are apparent between between the two aN3LO cases, are
relatively small across all PDFs, with slightly larger effects only where the
PDF itself tends towards zero i.e. valence quarks at small-$x$.

A larger distinction is observed when comparing the sets with and without
theoretical uncertainty (where N3LO theory is fixed at the best fit value). In
general there is an expected substantial increase in the PDF uncertainties
when taking into account the missing N3LO uncertainty for the gluon (and
therefore the heavy quarks). In particular, the form of the N3LO bottom quark
uncertainty is reminiscent of the $(H+\overline{H})$ prediction from Fig. 12.
One can therefore directly observe the effect of the $A_{Hg}$ theoretical
uncertainty on the bottom quark directly above its mass threshold. In other
areas, the without theoretical uncertainty PDF set exhibits a comparable
uncertainty to aN3LO and is even shown to increase the overall 68% confidence
intervals in certain regions of $(x,Q^{2})$ due to N3LO parameters being fixed
(i.e. $u_{v}$ and $d_{v}$ PDFs in Fig. 39 and Fig. 41). As the fit now resides
in a different $\chi^{2}$ landscape where a best fit has been achieved through
fitting the N3LO theory, fixing the aN3LO theory parameters is likely to have
a substantial effect across all PDFs.

An important point made by Fig.’s 39 and 41 is that that the difference
between the decorrelated and correlated cases is much smaller than the
difference of not including theoretical uncertainties at all (blue shaded
region). This analysis therefore provides evidence to support the original
assumption of being able to decorrelate the cross section (aN3LO $K$-factors)
and PDF theory (including other N3LO theory).

Along with the separate cases of uncertainty illustrated in Fig.’s 39 and 41,
we also display the central values of an aN3LO fit to all non-HERA data and an
aN3LO fit with NNLO $K$-factors. Examining the form of the no HERA aN3LO PDFs
for $x>10^{-2}$, we show some agreement with the standard N3LO central value
across most PDFs (more so at high-$Q^{2}$ than low-$Q^{2}$). Whereas the form
at small-$x$ gives some insight into the importance of HERA data in
constraining PDFs in this region. In slightly better agreement across all $x$
are the aN3LO PDFs with NNLO $K$-factors, which compliment the $\chi^{2}$
results in Section 7 and Section 8.1 arguing that the form (and fit results)
of aN3LO PDFs is mostly determined from the extra PDF $+$ DIS coefficient
function N3LO additions i.e. not aN3LO $K$-factors which prefer a softer
high-$x$ gluon (similar to the N3LO no HERA case – also shown in Fig.’s 39 and
41).

###  8.5 MSHT20aN3LO PDFs at $Q^{2}=2\ \mathrm{GeV}^{2}$

Figure 43:  General forms of NNLO (left) and aN3LO (right) PDFs at $Q^{2}=2\
\mathrm{GeV}^{2}$. Axis are set to the same scale to highlight the main
differences between NNLO and aN3LO. Specifically in the gluon and heavy
flavour sectors.

Fig. 43 compares the MSHT NNLO and aN3LO PDF sets at $Q^{2}=2\
\mathrm{GeV}^{2}$. In this very low-$Q^{2}$ regime, some major differences are
evident between NNLO and aN3LO sets at $Q^{2}=2\ \mathrm{GeV}^{2}$, especially
towards small-$x$. For example, the gluon PDF is predicted to be much harder
across this region, such that it is now positive across all $x$ values
considered here. The effect of this can be immediately seen in the sea and
heavy quarks.

Since the charm quark is directly coupled to the gluon PDF (through a
convolution with $A_{Hg}$), the charm PDF receives a notable enhancement at
small-$x$ and also remains positive across all $x$ values
considered202020Since this is a convolution, it is the higher small-$x$ gluon,
combined with the high-$x$ enhancement of $A_{Hg}$ at N3LO which gives rise to
this increase in the charm PDF.. Another interesting feature is the reduction
in uncertainty of the strange quark at small-$x$. It may seem counter
intuitive to have an uncertainty reduction by adding sources of theoretical
uncertainty, however we should recall that the underlying theory has also been
altered. Although one can expect an uncertainty increase in PDFs across
$(x,Q^{2})$, there are exceptions to this e.g. where tensions are relieved by
introducing the N3LO theory. The shift in the $\chi^{2}$ landscape then has
the potential to result in more precise regions of $(x,Q^{2})$ (in this case
manifesting in an uncertainty reduction for the strange quark towards
small-$x$).

Figure 44: Very low-$Q^{2}$ ratio plots showing the aN3LO 68% confidence
intervals with decorrelated and correlated $K$-factor parameters, compared to
NNLO 68% confidence intervals. All plots are shown for $Q^{2}=2\
\mathrm{GeV}^{2}$.

Figure 45: (Continued) Very low-$Q^{2}$ ratio plots showing the aN3LO 68%
confidence intervals with decorrelated and correlated $K$-factor parameters,
compared to NNLO 68% confidence intervals. All plots are shown for $Q^{2}=2\
\mathrm{GeV}^{2}$.

Fig. 44 displays the ratios of the aN3LO MSHT PDFs to their NNLO counterparts
at $Q^{2}=2\ \mathrm{GeV}^{2}$. Here the specific shifts of each PDF are
displayed more clearly. We note that there are many similar features shown
here to those discussed for Fig.’s 35 and 37. Even in this very low-$Q^{2}$
regime, the uncertainty difference between correlated and decorrelated aN3LO
$K$-factor PDF sets is minimal in all relevant regions of $x$.

###  8.6 Effect of a $x<10^{-3}$ cut at aN3LO

In this section we include results from a global PDF fit with small-$x$
($x<10^{-3}$) data omitted. This analysis is provided to shed some light on
the tensions between regions of $x$ at aN3LO while also providing some context
with regards to the form of the PDFs in different regions of $x$.

Figure 46: Low-$Q^{2}$ PDF ratios showing aN3LO PDFs fitted with and without
small-$x$ ($<10^{-3}$) data included in a global fit. All plots are shown for
$Q^{2}=10\ \mathrm{GeV}^{2}$ with the exception of the bottom quark shown for
$Q^{2}=25\ \mathrm{GeV}^{2}$.

Immediately, in Fig. 46, one can observe that omitting all small-$x$ data
results in a set of less constrained PDFs for $x<10^{-3}$ (most notably in the
gluon sector). However, also in Fig. 46 it can be observed that overall, the
large-$x$ behaviour of these PDFs is very similar across both fits, indicating
that the full fit is able to sufficiently fit both large and small $x$ regions
simultaneously. We provide this analysis as a cross-check to further support
the reliability of our procedure, showing that the small-$x$ behaviour is not
overwhelmingly attempting to fit to any all order, or specifically small-$x$
resummation, result at the expense of the large-$x$ description. We also note
that while there is some definite change in the central values of the PDFs at
small $x$, in most cases this is very well within uncertainties, and at most
at the level of about one standard deviation, particularly for the gluon for
$x$ just below $10^{-3}$. This distinct, but limited shift in the best fit
PDFs suggests that effects beyond N3LO are certainly not insignificant at very
low $x$, they are also not dominating the pull on the fit.

###  8.7 Posterior N3LO Theory Parameters

Following from the previous section, it is also interesting to examine where
the aN3LO theory contributions and their uncertainties reside after a global
PDF fit.

Parameter | Default Fit |  Small-$x$ cut Fit
---|---|---
Central |  $+$ Limit |  $-$ Limit | Central |  $+$ Limit |  $-$ Limit
Low-$Q^{2}$ Coefficient |  |  |  |  |  |
$c_{q}^{\mathrm{NLL}}$ | $-3.868$ | $-1.891$ | $-6.132$ | $-5.822$ | $-4.333$ | $-10.373$
$c_{g}^{\mathrm{NLL}}$ | $-5.837$ | $-4.444$ | $-7.429$ | $-6.995$ | $-3.701$ | $-7.991$
Transition Matrix Elements |  |  |  |  |  |
$a_{qq,H}^{\mathrm{NS}}$ | $-64.411$ | $-38.778$ | $-91.850$ | $-65.103$ | $-40.073$ | $-93.225$
$a_{Hg}$ | $12214.000$ | $12966.856$ | $11279.376$ | $12524.000$ | $13831.976$ | $11286.843$
$a_{gg,H}$ | $-1951.600$ | $-1577.155$ | $-3418.568$ | $-1392.600$ | $-512.817$ | $-2190.354$
Splitting Functions |  |  |  |  |  |
$\rho_{qq}^{NS}$ | $0.007$ | $0.015$ | $-0.002$ | $0.006$ | $0.020$ | $-0.005$
$\rho_{qq}^{PS}$ | $-0.501$ | $-0.254$ | $-0.644$ | $-0.505$ | $-0.285$ | $-0.692$
$\rho_{qg}$ | $-1.754$ | $-1.157$ | $-1.897$ | $-1.309$ | $-0.620$ | $-1.881$
$\rho_{gq}$ | $-1.784$ | $-1.548$ | $-2.212$ | $-1.622$ | $-1.367$ | $-1.877$
$\rho_{gg}$ | $19.245$ | $21.505$ | $9.025$ | $12.997$ | $16.142$ | $6.611$
Table 19: Posterior predicted $\pm 1\sigma$ limits on aN3LO theoretical
nuisance parameters for splitting functions, transition matrix elements and
coefficient functions.

Displayed in Table 19 are the predicted posterior limits on each (non
$K$-factor) aN3LO theory parameter. Here one can directly compare these
variations with the prior variations decided in earlier Sections. Also
provided is a comparison of these posterior limits across a fit with and
without small-$x$ ($x<10^{-3}$) data included. This comparison compliments the
previous section by showing a similar trend in the central values predicted in
both cases (i.e. with an overlap of uncertainties). As with the PDFs, there is
some significant evidence of tensions, but these are not severe, and the
central values of many parameters are extremely stable. Furthermore this is
evidence that the small-$x$ behaviour is influencing, but not likely to be
dominating the behaviour of aN3LO parameters in a manner which is
significantly adversarial to the preference of data at $x>10^{-3}$.

Figure 47: Posterior variations of the aN3LO splitting functions and
transition matrix elements predicted from a full global fit (blue shaded band)
compared to the prior variations in each case (green shaded band).

Fig. 47 displays a comparison of the prior and posterior variations predicted
for the perturbative expansions of the relevant splitting functions and
transition matrix elements discussed in Section’s 4 and 5. We exclude the non-
singlet quantities from this comparison as the variations predicted for these
quantities are very similar to their priors (as can be seen in Table 19) and
have a small overall effect on the PDFs. It is true that once a fit is
performed, the variation of the aN3LO theoretical nuisance parameters becomes
less sensitive to the prior variation, suggesting that the the initial
uncertainty estimate was conservative. Nevertheless in Fig. 47, one can
observe that all posterior variations overlap with their corresponding priors,
in most cases quite considerably. We also note that the most drastic
differences between prior and posterior variations are as expected relating to
the gluon PDF.

Figure 48: Posterior variations of the aN3LO splitting functions and
transition matrix elements predicted from a full global fit (blue shaded band)
compared to a fit with small-$x$ ($x<10^{-3}$) data removed from a fit (red
shaded band).

Fig. 48 contains a comparison between the aN3LO functions posterior variations
with and without small-$x$ ($x<10^{-3}$) data included in a global fit. These
results accompany those presented in Table 19 and further display the
reasonable agreement between the two fits, but there is some degree of tension
occurring mainly in the cases of $P_{gg}$, $A_{Hg}$ and $A_{gg,H}$. In all
cases the predicted variations overlap, with most central values being stable
(i.e. contained well within the uncertainty predictions). However, for
$P_{gg}$, $A_{Hg}$ and $A_{gg,H}$ the fit with the small $x$ cut does result
in posterior functions which are more consistent with the prior functions,
again suggesting that for these functions the posterior values are influenced,
to a significant, but not overwhelming extent by terms beyond N3LO, most
likely those associated with small-$x$ resummation. Hence, as with the PDFs
this provides evidence that the aN3LO predictions are reasonably consistent
across all values of $x$ but are influenced to a limited extent by the
small-$x$ region. This supports our view that while we are explicitly
determining the missing N3LO corrections, which are indeed overall the
dominant part of the missing higher order corrections, the fit is also probing
some even higher order corrections, particularly at small $x$.

###  8.8 N3LO Contributions

In this section all but one N3LO contribution will be switched off, in
particular only splitting functions, or only heavy or light flavour
coefficient functions with their relevant transition matrix elements. In all
cases the aN3LO $K$-factors are left free to allow the fit some freedom in
manipulating the cross sections of other datasets. In practice however, fixing
these $K$-factors at the NNLO values has a minimal effect on the shape of the
PDFs in all cases (as demonstrated in Fig. 39 and 41).

Figure 49: Low-$Q^{2}$ PDF ratios showing aN3LO (with decorrelated
$K$-factors) 68% confidence intervals compared to NNLO 68% confidence
intervals with varying theory contributions. All plots are shown for
$Q^{2}=10\ \mathrm{GeV}^{2}$ with the exception of the bottom quark shown for
$Q^{2}=25\ \mathrm{GeV}^{2}$. The PDFs included are: NNLO (green shaded), All
N3LO contributions (blue shaded), only splitting functions (green dashed),
only heavy flavour coefficient functions and transition matrix elements (dark
grey dash-dot) and only light flavour coefficient functions and transition
matrix elements (red dotted).

Figure 50: (Continued) Low-$Q^{2}$ PDF ratios showing aN3LO (with decorrelated
$K$-factors) 68% confidence intervals compared to NNLO 68% confidence
intervals with varying theory contributions. All plots are shown for
$Q^{2}=10\ \mathrm{GeV}^{2}$ with the exception of the bottom quark shown for
$Q^{2}=25\ \mathrm{GeV}^{2}$. The PDFs included are: NNLO (green shaded), All
N3LO contributions (blue shaded), only splitting functions (green dashed),
only heavy flavour coefficient functions and transition matrix elements (dark
grey dash-dot) and only light flavour coefficient functions and transition
matrix elements (red dotted).

The deconstructed aN3LO PDFs as a ratio to the NNLO MSHT PDFs for various
flavours at $Q^{2}=10\ \mathrm{GeV}^{2}$ (with the bottom quark given at
$Q^{2}=25\ \mathrm{GeV}^{2}$) are shown in Fig. 49. Across the more tightly
constrained light quark PDFs, all contributions lie very close to the aN3LO
$\pm 1\sigma$ uncertainty bands (blue shaded region and solid line). The
additive and compensating nature of these contributions is also clear in a
handful of the ratios from Fig. 49. In other areas the full description is
biased towards a single contribution, for example the charm and bottom quarks
follow the contribution from heavy flavours as one may expect. Conversely, to
some extent the gluon follows the splitting functions much more closely as
these contributions indirectly couple the gluon to the more constraining
data212121An exception to this can be seen around $x\sim 10^{-2}$ where the
contributions act cumulatively. We make this point as this region of $x$ is of
interest for Higgs calculations such as those discussed in Section 9.

###  8.9 $\alpha_{s}$ Variation

Figure 51: Quadratic fit to the total $\chi^{2}$ results from various
$\alpha_{s}(m_{Z})$ starting scales. The minimum of the quadratic fit provides
a rough estimate of $\alpha_{s}(m_{Z})=0.1170$ at aN3LO.

As in the standard MSHT20 NNLO PDF fit, we present the best fit aN3LO PDFs
with $\alpha_{s}(m_{Z})=0.118$, the common value chosen in the PDF4LHC
combination [8]. However, investigating the true minima in
$\alpha_{s}(m_{Z})$, the $\chi^{2}$ profiles in Fig. 51 prefer a value of
around $\alpha_{s}(m_{Z})=0.1170$. This result follows the trend from lower
orders whereby the best fit values are $\alpha_{s}(m_{Z})=0.1174\pm 0.0013$ at
NNLO and $\alpha_{s}(m_{Z})=0.1203\pm 0.0015$ at NLO [151]. Following from
NNLO, the aN3LO $\alpha_{s}(m_{Z})$ prediction is also slightly lower than the
NNLO world average central value at around $\alpha_{s}(m_{Z})=0.1179\pm
0.0010$ [152]. In any case, the preferred aN3LO $\alpha_{s}(m_{Z})$ value
stated here is in agreement with the MSHT20 NNLO result and the world average
within uncertainties. A full analysis is left for a future publication.

###  8.10 Charm Mass Dependence

Figure 52: Quadratic fit to the total $\chi^{2}$ results from various charm
masses ($m_{c}$). The minimum of the quadratic fit provides a rough estimate
of $m_{c}=1.45\ \mathrm{GeV}$ at aN3LO.

Figure 53: Low-$Q^{2}$ PDF ratios showing aN3LO (with decorrelated
$K$-factors) 68% confidence intervals compared to NNLO 68% confidence
intervals with varying fixed values for the charm mass. All plots are shown
for $Q^{2}=10\ \mathrm{GeV}^{2}$ with the exception of the bottom quark shown
for $Q^{2}=25\ \mathrm{GeV}^{2}$. The PDFs included are: $m_{c}=1.40\
\mathrm{GeV}$ (standard MSHT20 choice) (blue solid), $m_{c}=1.30\
\mathrm{GeV}$ (green dashed), $m_{c}=1.45\ \mathrm{GeV}$ (grey dotted dashed)
$m_{c}=1.50\ \mathrm{GeV}$ (red dotted).

Figure 54: (Continued) Low-$Q^{2}$ PDF ratios showing aN3LO (with decorrelated
$K$-factors) 68% confidence intervals compared to NNLO 68% confidence
intervals with varying fixed values for the charm mass. All plots are shown
for $Q^{2}=10\ \mathrm{GeV}^{2}$ with the exception of the bottom quark shown
for $Q^{2}=25\ \mathrm{GeV}^{2}$. The PDFs included are: $m_{c}=1.40\
\mathrm{GeV}$ (standard MSHT20 choice) (blue solid), $m_{c}=1.30\
\mathrm{GeV}$ (green dashed), $m_{c}=1.45\ \mathrm{GeV}$ (grey dotted dashed)
$m_{c}=1.50\ \mathrm{GeV}$ (red dotted).

In a standard MSHT fit [3], aN3LO PDFs are produced with the charm pole mass
$m_{c}=1.40\ \mathrm{GeV}$. Fig. 52 displays the $\chi^{2}$ results when
varying this charm mass. The predicted minimum at NNLO (for MSHT20 PDFs) is in
the range $m_{c}=1.35-1.40\ \mathrm{GeV}$ [151], whereas at aN3LO we show a
minimum in the region of $m_{c}=1.42-1.47\ \mathrm{GeV}$. This aN3LO result
therefore shows a slightly better agreement with the world average
[152]222222There is some ambiguity in this value since the transformation from
$\overline{MS}$ to the pole mass definition is not well-defined (see [151] for
more details). of $m_{c}=1.5\pm 0.2\ \mathrm{GeV}$.

Considering Fig. 53, one is then able to analyse the effect of this slightly
higher charm mass on the form of the PDFs. As one can expect, the charm PDF is
subject to the largest difference and is suppressed by a higher $m_{c}$. The
extra suppression from a higher charm mass allows the fit to suppress the
$c+\bar{c}$ sea contribution. This is then compensated by an increase in the
$\bar{u}$ and $\bar{d}$ distributions which stabilises the overall sea
contribution.

##  9 N3LO Predictions

With the increasing number of hard cross section calculations at N3LO, there
is a growing demand for N3LO accuracy in PDFs. In this section we investigate
the effect of the MSHT approximate N3LO PDFs on Higgs production via gluon
fusion and vector boson fusion (VBF). The hard cross sections for these
processes have been calculated to N3LO accuracy [67, 68, 69, 70, 71, 72, 73,
74, 75, 76, 153, 154]. We present a full N3LO computation for each prediction
with our approximate N3LO PDFs, including theoretical uncertainties. In future
work, the intention will be to expand this analysis to include results for
N3LO DY [63] and approximate N3LO top production [66] cross sections.

Note that in this Section we follow the notation used previously and denote
the aN3LO results with decorrelated $K$-factors as $(H_{ij}+K_{ij})^{-1}$ and
those with correlated $K$-factors with $H_{ij}^{\prime\ -1}$. In all cases,
scale variations are found via the 9-point prescription [11] for results with
NNLO PDFs. Whereas for aN3LO PDFs, although the extra information introduced
is at N3LO, the data (and therefore all relevant theory nuisance parameters)
which are included in the global fit are sensitive to all orders. In
particular, we include theoretical uncertainties into our aN3LO fit which
incorporate MHO effects on the PDFs. Therefore we argue (and in these cases
demonstrate) that the factorisation scale variation is contained within the
PDF uncertainties. Due to this, it is only the renormalisation scale which
requires variation in predictions involving aN3LO PDFs232323This is to
quantify the theoretical MHOU in the hard cross section, whereas the aN3LO
PDFs now come with an estimated MHOU..

###  9.1 Higgs Production – Gluon Fusion: $gg\rightarrow H$

$\sigma$ order PDF order $\sigma+\Delta\sigma_{+}-\Delta\sigma_{-}$ (pb)
$\sigma$ (pb) $+\ \Delta\sigma_{+}-\Delta\sigma_{-}$ (%) PDF uncertainties
N3LO aN3LO (no theory unc.) $45.296+0.723-0.545$ $45.296+1.60\%-1.22\%$ aN3LO
($H_{ij}+K_{ij}$) $45.296+0.832-0.755$ $45.296+1.84\%-1.67\%$ aN3LO
($H_{ij}^{\prime}$) $45.296+0.821-0.761$ $45.296+1.81\%-1.68\%$ NNLO
$47.817+0.558-0.581$ $47.817+1.17\%-1.22\%$ NNLO NNLO $46.206+0.541-0.564$
$46.206+1.17\%-1.22\%$ PDF + Scale uncertainties N3LO aN3LO (no theory unc.)
$45.296+0.723-1.851$ $45.296+1.60\%-4.09\%$ aN3LO ($H_{ij}+K_{ij}$)
$45.296+0.832-1.923$ $45.296+1.84\%-4.25\%$ aN3LO ($H_{ij}^{\prime}$)
$45.296+0.821-1.926$ $45.296+1.81\%-4.25\%$ NNLO $47.817+0.577-2.210$
$47.817+1.21\%-4.62\%$ NNLO NNLO $46.206+4.284-5.414$ $46.206+9.27\%-11.72\%$

Table 20: Higgs production cross section results via gluon fusion (with
$\sqrt{s}=13\ \text{TeV}$) using N3LO and NNLO hard cross sections combined
with NNLO and aN3LO PDFs. All PDFs are at the standard choice
$\alpha_{s}(m_{Z})=0.118$. These results are found with $\mu=m_{H}/2$ unless
stated otherwise, with the values for $\mu=m_{H}$ supplied in Table D.1.

Figure 55: Higgs production cross section results via gluon fusion (with
$\sqrt{s}=13\ \text{TeV}$) at two central scales: $\mu=m_{H}/2$ (left) and
$\mu=m_{H}$ (right). Displayed are the results for aN3LO PDFs with
decorrelated $K$-factors ($(H_{ij}+K_{ij})^{-1}$), correlated $K$-factors
($H_{ij}^{\prime\ -1}=(H_{ij}+K_{ij})^{-1}$) each with a scale variation band
from varying $\mu_{r}$ by a factor of 2. In the NNLO and NLO PDF cases, both
scales $\mu_{f}$ and $\mu_{r}$ are varied by a factor of 2 following the
9-point convention [11].

Table 20 and Fig. 55 (left) show predictions at a central scale of
$\mu=\mu_{f}=\mu_{r}=m_{H}/2$ for the Higgs production cross section via gluon
fusion242424Results are obtained with the code ggHiggs[67, 68, 69, 70, 71, 72,
73, 74, 75, 76, 155]. at the LHC for $\sqrt{s}=13\ \text{TeV}$, where
$m_{H}=125\ \text{GeV}$ is the Higgs mass and no fiducial cuts are applied.
Fig. 55 (right) displays the same analysis for the gluon fusion cross section
with $\mu=\mu_{f}=\mu_{r}=m_{H}$ (numerical results provided in Table D.1).

Considering the $\mu=m_{H}/2$ and $\mu=m_{H}$ central value results displayed
in Table 20 and Fig. 55, it can be observed that aN3LO PDFs predict a lower
central value than NNLO PDFs across all hard cross section orders. One can
also notice an overlap in all cases between predictions from NNLO and aN3LO
PDFs. However for $\mu=m_{H}/2$, whilst the error bands for predictions with
N3LO hard cross section and NNLO and N3LO PDFs overlap, their central values
are outside each other’s respective error bands. Since estimating MHOUs via
scale variations is a somewhat ambiguous procedure (and is therefore estimated
conservatively to reflect this), these results highlight the benefit of being
able to exploit a higher level of control over MHOUs i.e. via nuisance
parameters. By predicting a different central value we include a more accurate
estimation for higher order predictions which may not be contained within
scale variations, especially at unmatched orders in perturbation theory.

Examining the predicted central values further, Fig. 55 suggests that the
increase in the cross section theory at N3LO is compensated by the PDF theory
at N3LO, suggesting a cancellation between terms in the PDF and cross section
theory at N3LO. This point is important to consider when combining unmatched
orders in physical calculations, since we must be open to the possibility that
unmatched cancellations in physical calculations can lead to inaccurate
predictions, as our results suggest here.

Further to this, the change in the gluon PDF is largely driven by the
predicted form of $P_{qg}$ at aN3LO and DIS data. Therefore the relevant
changes in the gluon at aN3LO are most likely due to indirect effects i.e. not
directly related to gluon fusion predictions. Due to this, there is no reason
to believe that the observed level of convergence should happen at aN3LO for
both choices of $\mu$. However, owing to the inclusion of known information at
higher orders, one can be confident that the prediction is more accurate than
NNLO, whichever way it moves.

Comparing PDF uncertainty values calculated using NNLO and aN3LO PDFs, another
prominent feature one can notice in Table 20 is an increase in PDF
uncertainties. We find that the PDF uncertainty without N3LO theory
uncertainties included (i.e. using only the eigenvector description from the
first 32 eigenvectors and with N3LO parameters fixed at the best fit) also
includes a marginal increase in the positive direction compared to NNLO.
Mathematically, the reason for this comes back to the fact that the best fit
is inherently different from the NNLO theory, residing in a completely novel
$\chi^{2}$ landscape. In turn, this means it is not guaranteed that the PDF
uncertainty will remain consistent across the distinct PDF sets252525As we can
see from Section 8.4, the theory uncertainty is also not guaranteed to add to
the total uncertainty (and in fact acts to reduce the uncertainty in some
areas of $(x,Q^{2})$).. In the case of gluon fusion, the leading contribution
to the positive uncertainty direction is an eigenvector primarily dominated by
PDF parameters, while in the negative direction a N3LO splitting function
parameter dominates (eigenvector 9 and 31 in the $(H_{ij}+K_{ij})^{-1}$ N3LO
case – see Table 13). As discussed in Section 8.8, the gluon predominantly
follows the splitting function contributions, therefore it is not surprising
that this eigenvector is having a noticeable effect. Phenomenologically, the
increase in predicted uncertainties from the inclusion of the theoretical
uncertainties is a reflection of the estimated PDF MHOUs in this particular
cross section, and acts to replace factorisation scale variation. As a
consistency check, we find that when performing a 9-point scale variation
procedure with aN3LO PDFs, the values calculated (for both choices of $\mu$)
are within the predicted PDF uncertainties. This is therefore a further
verification of our MHOUs and that the $\mu_{f}$ variation is intrinsic in the
PDF uncertainties.

Finally Fig. 55 also demonstrates the increased stability of predictions when
considering the two different central scales $\mu$ at N3LO. As predicted from
perturbation theory, the scale dependence is reduced and central values become
more in agreement when increasing the order of either the PDFs or hard cross
section. Furthermore, the aN3LO $\sigma$ central predictions for both choices
of $\mu$ are contained within the uncertainty bands of each other. This is
true by definition for the NNLO PDFs since the factorisation scale $\mu_{f}$
variation includes both choices of $\mu$, whereas for aN3LO PDFs this result
is not guaranteed and is therefore intrinsic in the PDF (and renormalisation
scale $\mu_{r}$ variation) uncertainty.

###  9.2 Higgs Production – Vector Boson Fusion: $qq\rightarrow H$

Figure 56: Higgs production cross section results via vector boson fusion
(with $\sqrt{s}=13\ \text{TeV}$) at a central scale set to the vector boson
momentum. Displayed are the results for aN3LO PDFs with decorrelated
$K$-factors ($(H_{ij}+K_{ij})^{-1}$), correlated $K$-factors ($H_{ij}^{\prime\
-1}=(H_{ij}+K_{ij})^{-1}$) each with a scale variation band from varying
$\mu_{r}$ by a factor of 2. In the NNLO and NLO PDF cases, both scales
$\mu_{f}$ and $\mu_{r}$ are varied by a factor of 2 following the 9-point
convention [11].

$\sigma$ order PDF order $\sigma+\Delta\sigma_{+}-\Delta\sigma_{-}$ (pb)
$\sigma$ (pb) $+\ \Delta\sigma_{+}-\Delta\sigma_{-}$ (%) PDF uncertainties
N3LO aN3LO (no theory unc.) $4.1150+0.0638-0.0724$ $4.1150+1.55\%-1.76\%$
aN3LO ($H_{ij}+K_{ij}$) $4.1150+0.0682-0.0755$ $4.1150+1.66\%-1.83\%$ aN3LO
($H_{ij}^{\prime}$) $4.1150+0.0678-0.0742$ $4.1150+1.65\%-1.80\%$ NNLO
$3.9941+0.0558-0.0631$ $3.9941+1.40\%-1.58\%$ NNLO NNLO $3.9974+0.0557-0.0633$
$3.9974+1.39\%-1.58\%$ PDF + Scale uncertainties N3LO aN3LO (no theory unc.)
$4.1150+0.0638-0.0724$ $4.1150+1.55\%-1.76\%$ aN3LO ($H_{ij}+K_{ij}$)
$4.1150+0.0683-0.0755$ $4.1150+1.66\%-1.83\%$ aN3LO ($H_{ij}^{\prime}$)
$4.1150+0.0678-0.0742$ $4.1150+1.65\%-1.80\%$ NNLO $3.9941+0.0560-0.0631$
$3.9941+1.40\%-1.58\%$ NNLO NNLO $3.9974+0.0576-0.0642$ $3.9974+1.44\%-1.61\%$

Table 21: Higgs production cross section results via the vector boson fusion
process (with $\sqrt{s}=13\ \text{TeV}$) using N3LO and NNLO hard cross
sections combined with NNLO and aN3LO PDFs. All PDFs are at the standard
choice $\alpha_{s}(m_{Z})=0.118$. These results are found with $\mu=Q^{2}$
where $Q^{2}$ is the vector boson momentum.

$\sigma$ order PDF order $\sigma+\Delta\sigma_{+}-\Delta\sigma_{-}$ (pb)
$\sigma$ (pb) $+\ \Delta\sigma_{+}-\Delta\sigma_{-}$ (%) N3LO aN3LO $n_{f}=5$
4.1150 + 0.0683 - 0.0755 4.1150 + 1.66% - 1.83% aN3LO $n_{f}=4$ 4.0270 +
0.0685 - 0.0765 4.0270 + 1.70% - 1.90% aN3LO $n_{f}=3$ 2.7248 + 0.0653 -
0.0673 2.7248 + 2.40% - 2.47% NNLO NNLO $n_{f}=5$ 3.9974 + 0.0557 - 0.0633
3.9974 + 1.39% - 1.58% NNLO $n_{f}=4$ 3.9118 + 0.0561 - 0.0634 3.9118 + 1.44%
- 1.62% NNLO $n_{f}=3$ 2.6845 + 0.0539 - 0.0641 2.6845 + 2.01% - 2.39%

Table 22: Higgs production cross section results via the vector boson fusion
process (with $\sqrt{s}=13\ \text{TeV}$) using N3LO and NNLO hard cross
sections combined with NNLO and decorrelated aN3LO PDFs whilst varying the
number of active flavours $n_{f}$. All PDFs are at the standard choice
$\alpha_{s}(m_{Z})=0.118$. These results are found with $\mu=Q^{2}$ where
$Q^{2}$ is the vector boson momentum.

Table 21 and Fig. 56 show the predictions at various orders in $\alpha_{s}$
for Higgs production cross sections via vector boson fusion262626Results are
obtained with the inclusive part of the code proVBFH [153, 154, 156]. at the
LHC for $\sqrt{s}=13\ \text{TeV}$ up to N3LO [153, 154], again no fiducial
cuts are applied in this comparison. The predictions shown are calculated with
$\mu_{f}^{2}=\mu_{r}^{2}=Q^{2}$ as the central scale where $Q^{2}$ is the
vector boson squared momentum.

For this process one can follow the increase in the cross section as higher
order PDFs are used. Contrasting with the case of gluon fusion, Fig. 56
displays little cancellation between the terms added in the aN3LO PDF
description and the N3LO cross section. However, the cross section for VBF
produces around a $\sim 3-4\%$ change order by order and is therefore fairly
constant. Considering this relatively small difference between orders, this
lack of cancellation is not a major concern. Further to this, the vector boson
fusion process is much more reliant on the quark sector which, compared to the
gluon, is relatively constant order by order (see Section 8.4). The reason for
this stems from the more direct data constraints on the shape of quark PDFs.

Comparing the aN3LO VBF cross section (with MHO theoretical uncertainties)
with the NNLO cross section result (with NNLO PDFs) including MHOUs via scale
variations, we see that the scale variation MHOUs are negligible against the
PDF uncertainties at aN3LO. This result is in part due to the fact that the
scale variation for aN3LO is only being included for the renormalisation
scale. However at NNLO, the extra MHOU predicted was still only a small
contribution. Therefore considering these results further, the effects of
higher orders in both cases are expected to be small, which provides some
agreement with the argument that there is little scope for cancellation
between orders for VBF. As for the gluon fusion prediction in Section 9.1, we
confirm that any further factorisation scale variation (i.e. using the 9-point
prescription) is contained within the predicted PDF uncertainties; hence
further motivating our previous argument that factorisation scale variation is
not necessary with aN3LO PDFs.

Another feature of the VBF results is that the level of uncertainty at full
aN3LO is only increased slightly from the calculation involving NNLO PDFs.
Comparing this to the gluon fusion results, where the uncertainty was more
noticeably increased in both directions, it is evident that these approximate
N3LO additions are having a smaller effect on the VBF calculation. Once again,
the origin of this is due to the nature of the process. VBF involves mostly
the quark sector and is therefore much less affected by the extra N3LO theory
we have introduced (due to direct constraints from data). As we have presented
in previous sections, most of the uncertainty in the N3LO theory resides in
the small-$x$ regime which is more directly probed by the gluon sector than in
the quark sector.

Lastly we briefly discuss the $n_{f}$ dependence of the VBF cross section. In
VBF the scaling of contributions follows as $n_{f}^{2}$ due to the presence of
two input quark flavours in the process. In Table 22 we observe that the VBF
cross section receives a large contribution when including the charm quark
($n_{f}=3\rightarrow 4$) due to this scaling. We also show that at aN3LO, this
is where most of the difference in the central value and uncertainty from NNLO
is accounted for. This is a consequence of the predicted enhancement of the
charm PDF at aN3LO, discussed in Section 8.4. Beyond $n_{f}=4$ the bottom
contribution to VBF in the $W^{\pm}$ channel (the dominant channel) is heavily
suppressed, since due to the CKM elements $b$ must transition to $t$ most of
the time. Therefore the VBF cross section only receives a small contribution
moving from $n_{f}=4$ to $n_{f}=5$.

##  10 Availability and Recommended Usage of MSHT20 aN3LO PDFs

We provide the MSHT20 aN3LO PDFs in LHAPDF format [157]:

http://lhapdf.hepforge.org/

as well as on the repository:

http://www.hep.ucl.ac.uk/msht/

The approximate N3LO functions (for $P_{ij}(x)$ and $A_{ij}(x)$) are provided
as lightweight FORTRAN functions or as part of a Python framework in the
repository:

https://github.com/MSHTPDF/N3LO_additions

We present the aN3LO eigenvector sets with and without correlated $K$-factors
as discussed in Section 8, with the default set being provided with
decorrelated $K$-factors272727These grids are updated from a previous version
of the MSHT20 aN3LO PDFs and should be used in favour of any sets downloaded
before the latest upload date of this submission..

MSHT20an3lo_as118
MSHT20an3lo_as118_Kcorr

Both these PDF sets contain a central PDF accompanied by 104 eigenvector
directions (describing 52 eigenvectors) and can be used in exactly the same
way as previous MSHT PDF sets i.e. the MSHT20 NNLO PDFs with 64 eigenvector
directions.

As presented in this work, the aN3LO PDFs include an estimation for missing
N3LO contributions (the leading theoretical uncertainty) and implicitly some
MHOU beyond this within their PDF uncertainties. Due to this, we argue and
motivate in Section 9 that factorisation scale variations are no longer
necessary in calculations involving aN3LO PDFs. However the renormalisation
scale should continue to be varied to provide estimates of MHOUs in the hard
cross section piece of physical calculations.

In the case that the hard cross section for a process is available up to N3LO
the recommendation is to use the aN3LO PDFs, since unmatched ingredients in
cross section calculations can ignore important cancellations (between the
PDFs and hard cross section).

If a process is included within the global fit and the hard cross section is
known only up to NNLO (i.e. those discussed in Section 7), we recommend the
use of the decorrelated version of the aN3LO PDF set. Using these PDFs and the
details provided in Table 12, the hard cross section can be transformed from
NNLO to approximate N3LO. From here the two approximate N3LO ingredients can
be used together to give a full approximate N3LO result.

If a process is not included in the global PDF fit and the hard cross section
is known only up to NNLO, the standard NNLO PDF set remains the default
choice. However, we recommend the use of these aN3LO PDFs as an estimate of
potential MHOUs. In this case the aN3LO PDF set + NNLO hard cross section
prediction should be reflected in any MHOU estimates for the full NNLO
prediction. For example, when the hard cross section is known only up to NNLO
Equation (3.13) from [75] can be adapted to be,

$\delta(\mathrm{PDF}-\mathrm{TH})=\frac{1}{2}\left\lvert\frac{\sigma_{\mathrm{aN}^{3}\mathrm{LO}}^{(2)}-\sigma_{\mathrm{NNLO}}^{(2)}}{\sigma_{\mathrm{aN}^{3}\mathrm{LO}}^{(2)}}\right\rvert$
(10.1)

where $\delta(\mathrm{PDF}-\mathrm{TH})$ is the predicted PDF theory