# On the Gouy-Chapman-Stern model of the electrical double-layer structure
with a generalized Boltzmann factor

Anis Allagui<EMAIL_ADDRESS>Dept. of Sustainable and Renewable Energy
Engineering, University of Sharjah, Sharjah, UAE  Center for Advanced
Materials Research, Research Institute of Sciences and Engineering, University
of Sharjah, Sharjah, UAE  Dept. of Mechanical and Materials Engineering,
Florida International University, Miami, FL33174, United States  Hachemi
Benaoum  Dept. of Applied Physics and Astronomy, University of Sharjah, PO Box
27272, Sharjah, United Arab Emirates  Oleg Olendski  Dept. of Applied Physics
and Astronomy, University of Sharjah, PO Box 27272, Sharjah, United Arab
Emirates

###### Abstract

The classical treatment of the electrical double-layer (EDL) structure at a
planar metal/electrolyte junction via the Gouy-Chapman-Stern (GCS) approach is
based on the Poisson equation relating the electrostatic potential to the net
mean charge density. The ions concentration in the diffuse layer are assumed
to follow the Boltzmann distribution law, i.e. $\propto\exp(-\tilde{\psi})$
where $\tilde{\psi}$ is the dimensionless electrostatic potential. However,
even in stationary equilibrium in which variables are averaged over a large
number of elementary stochastic events, deviations from the mean-value are
expected. In this study we evaluate the behavior of the EDL by assuming some
small perturbations superposed on top of its Boltzmann distribution of ion
concentrations using the Tsallis nonextensive statistics framework. With this
we assume the ion concentrations to be proportional to
$[1-(1-q)\tilde{\psi}]^{{1}/{(1-q)}}=\exp_{q}({-\tilde{\psi}})$ with $q$ being
a real parameter that characterizes the system’s statistics. We derive
analytical expression and provide computational results for the overall
differential capacitance of the EDL structure, which, depending on the values
of the parameter $q$ can show both the traditional inverse bell-shaped curves
for aqueous solutions and bell curves observed with ionic liquids.

###### keywords:

Double-layer capacitor , Tsallis distribution , Boltzmann distribution ,
Capacitance

††journal: Physica A

##  1 Introduction

The electrical double-layer (EDL) is more than a century-old fundamental topic
of electrochemistry and physical chemistry [1, 2, 3, 4, 5]. It is found in
many applications such as in membrane science and technology [6, 7], water
desalination by capacitive ion removal [8, 9, 5], capacitive mixing [10, 11,
12], energy storage in porous electrodes and fractional-order capacitors [13,
4, 14, 15, 16, 17], colloidal suspensions [18, 19], etc. It refers to the
general phenomenon of charge separation across an interfacial surface, which
can be imagined as a plane that separates two adjacent and distinct phases
(e.g. metal and electrolyte) having excess electronic or ionic charges (see
Fig. 1). At such an electrified interface, the electrode charge is not
necessarily constant, but on the contrary fluctuates in response to the
stochastic thermal motions in the adjacent electrolyte. These fluctuations are
related to a set of configuration variables including the number, positions
and momenta of mobile ions of each species near the electrode, and become
important as the correlation between species increases [20, 21]. The extent of
these correlations is related to the nature of the interfacial liquid, and is
expected to increase with the increase of ions concentration as well as in
nanometer-sized confined geometries at the electrode surface. As a result the
electrical properties of within small scales of the electrode, such as the
differential capacitance, electrostatic potential or accumulated charge
density for instance, are expected to be affected by the statistics of these
microscopic fluctuations [22, 21].

Several computational models such as Monte-Carlo (MC), molecular dynamics (MD)
and density functional theory (DFT) simulations are now very apt to provide
good descriptions of the EDL structure [9, 22, 12, 23]. However, because of
the computational resources needed for the explicit microscopic analysis of
large numbers of interacting solvent and mobile ions, mean-field calculations
remain popular due their simplicity and comparatively good accuracy [5, 13,
4]. The classical Gouy-Chapman-Stern (GCS) approach in which the Poisson-
Boltzmann (PB) equations are used for the treatment of the diffuse part of the
solvent (reviewed in Sect. 1) is widely used because it is verified to be
asymptotically correct in the weak-coupling regime and it provides the
essential of the EDL properties [24, 1, 25, 5, 7, 26, 2]. It assumes the
equilibrium ion concentrations to decay exponentially with the distance from
the interface following a Boltzmann distribution function with respect to the
mean electrostatic potential energy in the EDL. This statistical
representation is considered and widely used as the generic case for thermal
equilibrium of many uncorrelated or weakly-correlated systems. However, when
effects such as ion-ion correlations, ion polarizability, finite size of ions,
electrostriction and dielectric saturation of the solvent cannot be ignored,
the GCS model becomes less reliable in describing the features of the EDL [2,
20, 27]. Furthermore, if the electrode surface is fractal and/or consisting of
small confined geometries, the stochastic nature of state variables associated
with the mobile ions for instance is also expected to deviate from Boltzmann’s
distribution profile [28].

The purpose of this study is to incorporate into the mean-field GCS model the
Tsallis $q$-deformed exponential distribution [29, 30, 31, 32, 33, 34, 28] for
the ion concentrations in the interfacial liquid (Sect. 2). The same analysis
can be applied in principle to other approaches dealing with the EDL
structure. The $q$-exponential function we are considering allows the presence
of random fluctuations of ion concentrations around their mean values, which
in turn provides a sort of a generalized Boltzmann factor. With this approach,
one can still analyze a microscopic system in which correlations are expected
as if they were absent [28]. We provide analytical solutions for the charge
density and differential capacitance from the boundary-value PB problem
parameterized with the value of $q$ that can be attributed to the strength of
elementary noise sources in the EDL structure. A similar approach that came to
our attention during the writing of this paper was reported by Garcia-Morales
et al. [28], in which the focus was mainly on the $q$-dependence of the PB
equation for describing counterion concentration around a charged surface.
However, this study is further concerned with the overall capacitive
performance of the EDL structure by taking into account the series combination
of the Helmholtz capacitance and the diffuse layer capacitance for a symmetric
electrolyte, which was not provided in ref. [28]. We also evaluate the effect
of the values of $q$ (taken between -1 to 1) on the EDL capacitance vs.
voltage profiles which allowed to capture both the traditional inverse bell-
shaped curves for aqueous solutions and bell-like curves measured for ionic
liquids.

##  2 The Gouy-Chapman-Stern Model

The standard model to describe the equilibrium EDL structure at the junction
metal/electrolyte (dilute solutions) is the GCS mean-field model which can
quantitatively explain most experimental results. It consists of subdividing
the EDL into (i) a first sheet of uniform constant electronic charges on the
metal surface, (ii) a charge-free inner compact layer or Helmholtz layer (also
known as Stern layer) of a few angstrom in width and constant charge, and
(iii) an outer, semi-infinite diffuse layer (also known as Gouy-Chapman layer)
that extends into the bulk electrolyte and contains anions and cations of a
certain distribution [24] (see Fig. 1). Ions are assumed to be ideal point
charges in local thermodynamic equilibrium within an isolated planar
interface, and the solvent is assumed to be a dielectric continuum [8]. There
are many subsequent modified and refined versions of the GCS model that apply
for instance corrections for ion crowding or dielectric saturation [35, 36,
37], but here we use the simplest treatment for one electrode only that
captures good enough the essentials. Only physical adsorption in the diffuse
layer is considered to compensate the electronic charges, but no Faradaic
charge-transfer reactions and/or ion adsorption are taken into account.

In the diffuse layer part of the GCS’s model, the electrolyte is treated using
the mean-field PB model. For dilute solutions with negligible correlations,
the ions’ concentrations at a perpendicular distance $x\geqslant 0$ from a
large lateral surface area metal electrode (one-dimensional problem in which
the edge effects are ignored) are assumed to be related to the mean
electrostatic potential $\psi_{x}$ by the continuous Boltzmann distribution
function:

$C_{x_{i}}=C_{x_{i}}^{0}\exp\left(-\frac{z_{i}F\psi_{x}}{RT}\right)$ (1)

where $C_{x_{i}}^{0}$ and $z_{i}$ are the average ion concentration (for
$\psi_{x}=0$) and charge number of species $i$ (cation and anion),
respectively. The term $RT/F=k_{B}T/e$ (with $R$ the ideal gas constant, $F$
the Faraday’s constant, $k_{B}$ the Boltzmann’s constant, $e$ the elemental
electronic charge, and $T$ the absolute temperature) is the thermal voltage.
The molar electrical interaction energy in the Boltzmann factor takes into
account the term $z_{i}F\psi_{x}$ only, and other contributions are ignored in
this study. The approximation of $C_{x_{i}}$ by a Boltzmann distribution is a
standard tool used to evaluate many other thermodynamic properties of systems,
but is valid for some ideal conditions only, i.e. elastic collision between
particles, particles are identical and independent of one another, etc. It
assumes that the ratio of equilibrium ionic concentrations at two locations
$x_{i}$ and $x_{j}$ in the solution, i.e. $C_{x_{i}}/C_{x_{j}}$, at which the
mean electrostatic potentials are $\psi_{{x_{i}}}$ and $\psi_{{x_{j}}}$,
respectively, is related to the work needed to carry an ion of charge $ze$
from $x_{i}$ to $x_{j}$, i.e, $-ze(\psi_{{x_{i}}}-\psi_{{x_{j}}})$, through
the exponential Boltzmann factor given in Eq. 1 [38].

Now the distribution of the mean electrostatic potential $\psi_{x}$ in Eq. 1
is governed by the Poisson equation that connects the net volume charge
density $\rho_{x}$ at the distance $x$ from the interface with the
corresponding potential $\psi_{x}$ as:

$\frac{\mathrm{d}^{2}\psi_{x}}{\mathrm{d}x^{2}}=-\frac{\rho_{x}}{\epsilon}$
(2)

where $\epsilon$ is the uniform dielectric permittivity of the solvent in
units of the vacuum dielectric constant $\epsilon_{0}$. The mean charge
density $\rho_{x}$ can be expressed as the sum
$\rho_{x}=\sum_{i}z_{i}FC_{x_{i}}$ over the ions species, which leads to the
second-order nonlinear differential equation:

$\frac{\mathrm{d}^{2}\psi_{x}}{\mathrm{d}x^{2}}=-\frac{1}{\epsilon}\sum_{i}z_{i}FC_{x_{i}}^{0}\exp\left(-\frac{z_{i}F\psi_{x}}{RT}\right)$
(3)

With the change of variables $p={\mathrm{d}\psi_{x}}/{\mathrm{d}x}$, such that
${\mathrm{d}^{2}\psi_{x}}/{\mathrm{d}x^{2}}=\mathrm{d}p/\mathrm{d}x=p\,\mathrm{d}p/\mathrm{d}\psi_{x}$,
we rewrite Eq. 3 as:

$\displaystyle\frac{1}{2}\left(\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}\right)^{2}$
$\displaystyle=\frac{RT}{\epsilon}\sum_{i}C_{x_{i}}^{0}\left[\exp\left(\frac{-z_{i}F\psi_{x}}{RT}\right)\right]+B$
(4)

The integration constant $B$ is determined using the boundary conditions
$\psi_{x}\,(x\to\infty)=0$ (electroneutrality) and
$\mathrm{d}\psi_{x}/\mathrm{d}x=0$ (electrostatic field must vanish at the
mid-plane) so that:

$\left(\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}\right)^{2}=\frac{2RT}{\epsilon}\sum_{i}C_{x_{i}}^{0}\left[\exp\left(\frac{-z_{i}F\psi_{x}}{RT}\right)-1\right]$
(5)

If we consider the special case of a symmetric electrolyte (i.e. two ionic
species; one positively charged and one negatively charged with $|z_{i}|=z$
and $C_{x_{i}}^{0}=C_{x}^{0},\forall i$), we obtain the analytical solution of
GCS model for a positively charged electrode as:

$\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}=-\sqrt{\frac{8RTC_{x}^{0}}{\epsilon}}\sinh\left(\frac{zF\psi_{x}}{2RT}\right)<0$
(6)

in which physical arguments impose that the appropriate sign should be as
shown. Note that Eq. 6 can be written in the dimensionless form
${\mathrm{d}(\tilde{\psi}_{x}/2)}/{\mathrm{d}\tilde{x}}=-\sinh(\tilde{\psi}_{x}/2)$,
where $\tilde{\psi}_{x}=\psi/(RT/zF)$ and $\tilde{x}=x/\lambda_{D}$ with
$\lambda_{D}=\sqrt{\epsilon RT/2z^{2}F^{2}C_{x}^{0}}$ being the Debye
screening length for the symmetric electrolyte. It can also be linearized in
the limit of low potentials leading to the Debye-Hückel model. With Eq. 6 one
can obtain the total net charge per unit surface $Q_{x}$ within the diffuse
layer as:

$\left(\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}\right)_{x=d}=-\frac{1}{\epsilon}\int\limits_{x=d}^{\infty}\rho_{x}\mathrm{d}x=-\frac{Q_{x}}{\epsilon}$
(7)

Combining Eq. 6 and Eq. 7 gives the total net charge within the diffuse layer
as:

$Q_{x}=\sqrt{{8\epsilon RTC_{x}^{0}}}\sinh\left(\frac{zF\psi_{d}}{2RT}\right)$
(8)

To determine the specific capacitance of the diffuse layer, the net charge
(Eq. 8) is differentiated with respect to the potential, which gives the well-
known nonlinear expression:

$C_{\text{diff}}=\sqrt{\frac{2\epsilon
z^{2}F^{2}C_{x}^{0}}{RT}}\cosh\left(\frac{zF\psi_{d}}{2RT}\right)$ (9)

that depends on both the potential $\psi_{d}$ and the properties of the
electrolyte $C_{x}^{0}$ and $\epsilon$. The expression for $C_{\text{diff}}$
can also be expressed in a dimensionless form as
$\tilde{C}_{\text{diff}}=\lambda_{D}C_{\text{diff}}/\epsilon=\cosh(\tilde{\psi}_{x}/2)$
[4]. The overall EDL capacitance of the electrode is then computed from the
combination in series of the two capacitances arising from the Helmholtz inner
layer (of constant capacitance $C_{\text{H}}=\epsilon/4\pi d$ with $d$ being
its thickness) and the diffuse layer, such that:

$C_{\text{dl}}^{-1}=C_{\text{H}}^{-1}+C_{\text{diff}}^{-1}$ (10)

A plot of $C_{\text{dl}}$ vs. potential ($\psi_{d}>0$) for different values of
the bulk concentrations $C_{x}^{0}$ is shown in Fig. 2. We took the values of
$C_{\text{H}}=\text{28}\,\mu\text{F\,cm}^{-2}$,
$\epsilon=\text{80}\,\epsilon_{0}$, $z=\text{1}$, and $T=\text{298}$ K. At
large voltages and for concentrated electrolytes, the overall capacitance is
determined by that of the Helmholtz layer, independently of the ion
concentrations, whereas as the voltage is decreased and for very dilute
electrolytes, the contribution of the diffuse layer becomes more important.
For $\psi_{d}=0$, we have the reciprocal of the overall capacitance
$C_{\text{dl}}^{-1}=C_{\text{H}}^{-1}+\left({{2\epsilon
z^{2}F^{2}C_{x}^{0}}/{RT}}\right)^{-0.5}$.

##  3 The Generalized Boltzmann Factor

The Boltzmann distribution function, as mentioned above, provides a valid
approximation only under certain assumptions for equilibrium thermodynamic
[39, 40, 25]. While it allows a relatively accurate description of macroscopic
systems in which a very large number of stochastic events take place, when one
goes to smaller scales for instance, such a description breaks down [41]. It
is expected that at small dimensions, fluctuations of some intensive quantity
such as temperature or pressure or chemical potential for example with respect
to that of the surrounding reservoir, and contributions from the surface
energy would play crucial roles in determining the state of the system under
consideration [42]. The system in this case is in stationary quasi-
thermodynamic equilibrium with the reservoir. This situation is closely
related to the case of the EDL structure within small sizes in the nanopores
of porous electrodes. The effects of fluctuations and long-range interactions
are actually more prominent at this level where the number of particles can be
much smaller than Avogadro’s number.

Moy et al. [7] have confirmed that some physical quantities, such as the force
on a test ion, computed using the PB approximation deviates from the
predictions with those of Brownian dynamics simulations (using the Langevin
equation) in which individual ions are treated explicitly. The difference
widens as the distance of the ion from the channel is lower than the Debye’s
length. The same discrepancies on the conductance and concentration profiles
in cylindrical channels and a potassium channel were reported between the
Poisson–Nernst–Planck theory (commonly used for non-equilibrium ion transport
problems) vs. the Brownian dynamics [26]. In another study by Härtel et al.
[9], DFT and MD simulations were applied to the EDL structure and showed very
good agreement for the ions concentration profiles, but clear deviations from
those computed with the PB theory are observed.

Motivated by these observations, the purpose of this study is to evaluate the
behavior of the EDL structure using a more general Boltzmann factor as
proposed by Tsallis [31], and later extended by Beck and Cohen [32, 33, 34]
and others (see ref. [29] and references within). These statistical mechanics
approaches have been introduced as a way to systematically describe complex
statistical systems that behave like the superposition of many Boltzmann
distributions, which makes the single Boltzmann distribution to be nothing but
a special limiting case. Such systems usually involve long-range interactions,
non-Markovian memory effects and anomalous diffusion; examples and
applications can be found in different disciplines of science and engineering
[43, 44]. At the microscopic level, systems relax toward thermodynamic
equilibrium (or more appropriately stationary nonequilibrium [45]) following a
single exponential decay, type Boltzmann, i.e. the probability of finding the
system at some specific state of energy $\epsilon$ is proportional to
$\exp(-\beta\epsilon)$ where $\beta$ is a local inverse temperature. At a
larger scale, the local variable $\beta$ (other variables can be considered)
is not well-defined and experiences fluctuations following a certain
probability distribution function $f(\beta)$, which depends on the underlying
dynamics of the system and is a priori unknown [45]. From Beck and Cohen [32],
the generalized Boltzmann factor can be written as the integral over all
possible fluctuating inverse temperatures $\beta$ as:

$B(E)=\int\limits_{0}^{\infty}f(\beta)\exp({-\beta E})\mathrm{d}\beta$ (11)

For most physical systems a distribution where the random variable $\beta$ is
nonnegative is needed. In particular, if $f(\beta)$ is assumed to follow a
gamma probability distribution function which arises naturally for a
fluctuating environment with a finite number of degrees of freedom, such that:

$f(\beta)=\frac{1}{b\Gamma(c)}\left({\beta}/{b}\right)^{c-1}\exp\left(-{\beta}/{b}\right)$
(12)

where $c$ and $b$ are positive parameters, Eq. 11 leads to the generalized
Boltzmann factor in Tsallis statistics [32, 29, 31, 46, 40, 44]:

$B(E)=(1+bE)^{-c}=\left[1-(1-q)\beta_{0}E\right]^{{1}/{(1-q)}}=\exp_{q}({-\beta_{0}E})$
(13)

with $1/(q-1)=c$, $\beta_{0}=\int_{0}^{\infty}\beta
f(\beta)\mathrm{d}\beta=bc$ is the average of the fluctuating $\beta$, and
$\exp_{q}(y)$ denotes the $q$-exponential function. The parameter $q$ is a
real number that characterizes the system’s statistics, and is defined by the
ratio of standard variation and mean of the distribution $f(\beta)$ ($q=1$ if
there are no fluctuations) [32]. Note that (i) for $q<1$, $\exp_{q}(y)=0$ for
$y<-1/(1-q)$ and $\exp_{q}(y)=\left[1+(1-q)y\right]^{1/(1-q)}$ for
$y\geqslant-1/(1-q)$, (ii) for $q=1$, $\exp_{q}(y)=\exp(y)$ for $\forall y$,
and (iii) for $q>1$, $\exp_{q}(y)=\left[1+(1-q)y\right]^{1/(1-q)}$ for
$y<1/(q-1)$ [44].

##  4 The Extended Gouy-Chapman-Stern Model

By applying generically the generalized Boltzmann factor given by Eq. 13 to
the concentrations of charged ions in the EDL structure (instead of Eq. 1), we
write the $q$-exponential relation [30]:

$C_{x_{i}}^{q}=C_{x_{i}}^{0}\exp_{q}\left(-\frac{z_{i}F\psi_{x}}{RT}\right)$
(14)

Again, this means that the fluctuations of the mean-field value of the ion
concentrations are assumed to follow a gamma distribution with the parameter
$q=1+1/c$ representing the extent of these fluctuations. The adequate
determination of the statistical distribution of such ions concentration is
very difficult to obtain. If $q=1$, the traditional Boltzmann factor given in
Eq. 1 is immediately recovered from Eq. 14. Plots of $C_{x_{i}}^{q}$ as a
function of potential $\psi_{x}$ for different values of $q\leqslant 1$ are
shown in Fig. 3. In a semi-logarithmic scale, it is only the Boltzmann
distribution that would result in a straight line with slope $-zF/RT\ln(10)$.
Otherwise, Decreasing the value of $q$ from unity implies steeper decrease of
the ionic concentration following a power-law profile as shown in the figure.
This means that contrary to the exponential decay in which subsystems are
uncorrelated, values of $q$ different from one implies some form of internal
correlations between subsystems.

The $q$-modified PB model becomes then:

$\mathrm{d}\left(\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}\right)^{2}=-\frac{2F}{\epsilon}\sum_{i}C_{x_{i}}^{0}z_{i}\exp_{q}\left(-\frac{z_{i}F\psi_{x}}{RT}\right)\mathrm{d}\psi_{x}$
(15)

which leads after integration to:

$\displaystyle\left(\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}\right)^{2}$
$\displaystyle=\frac{2RT}{\epsilon(2-q)}\times$
$\displaystyle\sum_{i}C_{x_{i}}^{0}\left[\left(1-(1-q)\frac{z_{i}F\psi_{x}}{RT}\right)\exp_{q}\left(-\frac{z_{i}F\psi_{x}}{RT}\right)-1\right]$
(16)

We verify that for $q=1$, Eq. 16 readily simplifies to the expression given in
Eq. 5. For the case of a symmetric 1:1 electrolyte, we obtain:

$\displaystyle\left(\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}\right)^{2}$
$\displaystyle=\frac{4RTC_{x}^{0}}{\epsilon(2-q)}\times$
$\displaystyle\left[\cosh_{q}\left(\frac{zF\psi_{x}}{RT}\right)+(1-q)\left(\frac{zF\psi_{x}}{RT}\right)\sinh_{q}\left(\frac{zF\psi_{x}}{RT}\right)-1\right]$
(17)

where the $q$-hyperbolic functions $\cosh_{q}$ and $\sinh_{q}$ are defined
respectively as [47]:

$\displaystyle 2\cosh_{q}(y)$ $\displaystyle=$
$\displaystyle{\exp_{q}(y)+\exp_{q}(-y)}$ (18a) $\displaystyle 2\sinh_{q}(y)$
$\displaystyle=$ $\displaystyle{\exp_{q}(y)-\exp_{q}(-y)}$ (18b)

Now using the $q$-identities:

$\displaystyle\cosh_{q}\left(y\right)$ $\displaystyle=$ $\displaystyle
2\sinh_{2q-1}^{2}\left({y}/{2}\right)+\exp_{2q-1}\left[-(1-q){y^{2}}/{2}\right]$
(19a) $\displaystyle\sinh_{q}\left(y\right)$ $\displaystyle=$ $\displaystyle
2\sinh_{2q-1}\left({y}/{2}\right)\cosh_{2q-1}\left({y}/{2}\right)$ (19b)

we obtain the following expression for the electric field:

$\displaystyle\frac{\mathrm{d}\psi_{x}}{\mathrm{d}x}$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{8RTC_{x}^{0}}{\epsilon(2-q)}}\sinh_{2q-1}\left(\frac{zF\psi_{x}}{2RT}\right)h_{q}\left(\frac{zF\psi_{x}}{RT}\right)$
(20)

where the function $h_{q}(y)$ is given by:

$h_{q}(y)=\left[1+(1-q)y\coth_{2q-1}\left(\frac{y}{2}\right)-\frac{1-\exp_{2q-1}\left(-(1-q)\frac{y^{2}}{2}\right)}{2\sinh_{2q-1}^{2}\left(\frac{y}{2}\right)}\right]^{0.5}$
(21)

It is clear that when $q\to 1$, the classical relation given by Eq. 6 for the
electric field is immediately recovered. We can write Eq. 20 in a
dimensionless form as
${\mathrm{d}(\tilde{\psi}_{x}/2)}/{\mathrm{d}\tilde{x}_{q}}=-\sinh_{2q-1}(\tilde{\psi}_{x}/2)h(\tilde{\psi}_{x})$
where $\tilde{x}_{q}$ is the normalized dimension $x$ with respect to a
modified Debye length, $\lambda_{D}^{q}=\lambda_{D}\sqrt{2-q}$, which depends
on the properties of the liquid as well as the extent of the fluctuations
through the parameter $q$. $\lambda_{D}^{q}$ is larger than $\lambda_{D}$ for
$q<1$ which implies that the scale over which the electrolyte screens the
surface charge is larger with the presence of fluctuations. Finally, the
$q$-parameterized net charge $Q_{x}^{q}$ in the diffuse layer using the
generalized Boltzmann factor is found as:

$Q_{x}^{q}=\sqrt{\frac{8\epsilon
RTC_{x}^{0}}{2-q}}\sinh_{2q-1}\left(\frac{zF\psi_{d}}{2RT}\right)h_{q}\left(\frac{zF\psi_{d}}{RT}\right)$
(22)

and its associated diffuse capacitance is found numerically from
$C_{\text{diff}}^{q}=\partial Q_{x}^{q}/\partial\psi_{x}$. The overall EDL
capacitance $C_{\text{dl}}^{q}$ is then computed as per Eq. 10, i.e.:

$(C_{\text{dl}}^{q})^{-1}=C_{\text{H}}^{-1}+(C_{\text{diff}}^{q})^{-1}$ (23)

The total charge $Q_{\text{dl}}^{q}$ is obtained from the voltage integral of
$C_{\text{dl}}^{q}$. Plots of $C_{\text{dl}}^{q}(\psi_{x})$ for different
values of $q$ for a constant bulk concentration of
$C_{x}^{0}=10^{-3}\,\text{mol\,L}^{-1}$, are shown in Fig. 4. The capacitance
at negative voltages exactly duplicate the capacitance at positive voltages
shown in the figure.

##  5 Discussion

From Fig. 4, it is evident that incorporating the generalized Boltzmann factor
for the ions concentrations into the mean-field GCS model affects greatly the
capacitive performance of the EDL structure. It is understood that other
electrical characteristics are also impacted, but the focus of the discussion
here is on the differential EDL capacitance, and thus energy storage
applications. The strength of the fluctuations increases as the value of $q$
deviates further from unity at which the traditional Boltzmann factor is
recovered. In other words, when the parameter $q$ is less than one, it implies
that there are some forms of hidden correlations within the stationary system
in quasi-thermodynamic equilibrium, which means that its constituting
subsystems are not statistically independent anymore as it is the case in
Boltzmann theory [45]. In practice, however, the determination of an
experimental value for $q$ and the particular form of correlations between the
dynamic subsystems is not trivial as mentioned above.

Fig. 4 shows also that the overall capacitance exhibits a local minimum at the
origin for different values of $q>0$, and an asymptotic increase towards the
value of $C_{\text{H}}=\text{28}\,\mu\text{F\,cm}^{-2}$ at higher
polarization. Otherwise, between these two limiting cases, the shapes of the
profiles for $0<q<1$ are in concordance with the shapes obtained with the
traditional GCS model for which $q=1$ (compare with Fig. 2). Now because of
the interdependence of internal subsystems for $q<1$, this results is an extra
contribution to the amount of entropy of the system, i.e. the entropy
$S_{q}(A+B)\geqslant S_{q}(A)+S_{q}(B)$ of two subsystems $A$ and $B$ [40].
This non-additivity of entropy is the basis of a non-extensive Tsallis
statistics. Tsallis proposed the generalized entropy
$S_{q}=-k\sum_{i}p_{i}^{q}\ln_{q}(p_{i})$ [31, 46, 40, 44], where $k$ is a
positive constant, $q\neq 1$, and the quantities $p_{i}=p(E_{i})$ represent
the probabilities for the occurrence of the $i^{\text{th}}$ microstate and
satisfy $\sum_{i}p_{i}=1$. The function $\ln_{q}(y)=(y^{1-q}-1)/(1-q)$ denotes
the $q$-logarithm, inverse of the $q$-exponential (i.e.
$\ln_{q}[\exp_{q}(y)]=\exp_{q}[\ln_{q}(y)]=y$) [40]. In the limit of $q\to 1$,
one recovers the classical form of Boltzmann-Gibbs entropy
$S_{1}=-k_{B}\sum_{i}p_{i}\ln(p_{i})$ where $k=k_{B}$ is the Boltzmann
constant [31]. Such additional amount of entropy, which is related to the
value of $q$, affects negatively the net density of charge in the diffuse
layer of the EDL structure compared to the traditional case when $q=1$. This
in turn lowers the overall EDL capacitance as depicted in Fig. 4.

What is also remarkable from these results is the bell-shaped-like trend of
the total capacitance when $q$ takes values less than zero. A similar profile
is observed for metal/ionic liquid double layer capacitance which exhibit a
local maximum at the point of zero charge [37]. In Fig. 4, we show the example
for $q=-1$ which is typical for all negative values of $q$. The capacitance
shows a peaking value and reverse of curvature close to the origin at
$\psi_{x}=0$, and the same is observed when the capacitance is plotted vs.
charge. It is well known that the non-linear PB theory of points ions is
unable to display a change of curvature observed experimentally in ionic
liquids or highly concentrated electrolytes. This change of curvature can be,
nonetheless, predicted in the so-called primitive model (in which solvent
effects are considered implicitly) by including ion correlations and ionic
excluded volume effects. This has been shown by numerous DFT and integral
equations calculations, as well as Monte Carlo, and Brownian dynamics and
molecular dynamics simulations. The change of curvature of the differential
capacity can be also observed in more detailed atomistic models in which
solvent particles are taken into account explicitly.

Finally, the authors recognize that the main limitation of this theoretical
description is the lack of a proper physical meaning for the parameter $q$ in
terms of a microscopic model. The analysis presented here is unable to provide
an effective description of ion correlations (i.e., the Coulombic attraction
of charged particles of opposite valence and the Coulombic repulsion of
charged particles with the same valence) and ionic excluded volume effects
(i.e., the fact that two ions cannot overlap due to its finite size), but this
is the subject of ongoing investigations. A measurable atomistic model
supporting the physical reality of the parameter $q$ will follow to validate
the use of non-extensive Tsallis statistics in EDL systems. In this sense, the
parameter $q$ can be viewed as mathematically well-fit to capture the change
of curvature of the differential capacity, as it is the case for instance of
fractional integro-differential equations for capturing anomalous time-domain
and frequency-domain features in several physicochemical systems [48, 16, 49,
50].

##  6 Conclusion

In this paper, we analyzed the effect of superposing fluctuations onto the
mean-value Boltzmann distribution of ion concentrations in the GCS model on
the capacitive properties of the EDL structure. This was done via the
$q$-exponential function which allows to embrace a spectrum of empirical
processes connected to the degree of fluctuations in the intensive inverse
temperature parameter. This compact and efficient approach provides an
extended version to the classical GCS model specifically for the description
of the EDL at scales in which correlation between subsystem are possible. The
analytical expression for the EDL capacitance as a function of voltage and ion
concentration, and parameterized with the index $q<1$, predicts lower values
than its GCS counterpart (i.e. for $q=1$) which has been explained by super-
additive entropy feature in Tsallis statics. Additionally, for values of $q$
less than zero, the model reverses the general trend of capacitance at a
metal/electrolyte junction towards that of a bell-shaped-like profiles which
are observed in ionic liquid EDL structure. Further investigations in this
direction with comparison with atomistic models are being developed for future
contributions.

## Acknowledgement

We thank Mathijs Janssen for useful discussions and comments on an earlier
version of the manuscript.

## References

## References

  * Delahay [1965] P. Delahay, Double layer and electrode kinetics, John Wiley & Sons Inc, 1965.
  * Carnie and Torrie [1984] S. L. Carnie, G. M. Torrie, The Statistical Mechanics of the Electrical Double Layer, vol. 58, chap. 2, John Wiley & Sons, Ltd, 141–253, doi:10.1002/9780470142806.ch2, 1984.
  * Parsons [1990] R. Parsons, The electrical double layer: recent experimental and theoretical developments, Chemical Reviews 90 (5) (1990) 813–826.
  * Bazant et al. [2004] M. Z. Bazant, K. Thornton, A. Ajdari, Diffuse-charge dynamics in electrochemical systems, Phys. Rev. E 70 (2) (2004) 021506\.
  * Biesheuvel and Dykstra [2020] P. Biesheuvel, J. Dykstra, Physics of Electrochemical Processes, www.physicsofelectrochemicalprocesses.com, 2020\.
  * McLaughlin [1989] S. McLaughlin, The electrostatic properties of membranes, Annu. Rev. Biophys. Biophys. Chern. 18 (1) (1989) 113–136.
  * Moy et al. [2000] G. Moy, B. Corry, S. Kuyucak, S.-H. Chung, Tests of continuum theories as models of ion channels. I. Poisson- Boltzmann theory versus Brownian dynamics, Biophys. J. 78 (5) (2000) 2349–2363.
  * Biesheuvel and Bazant [2010] P. Biesheuvel, M. Bazant, Nonlinear dynamics of capacitive charging and desalination by porous electrodes, Phys. Rev. E 81 (3) (2010) 031502\.
  * Härtel et al. [2015] A. Härtel, M. Janssen, S. Samin, R. van Roij, Fundamental measure theory for the electric double layer: implications for blue-energy harvesting and water desalination, J. Phys.: Condens. Matter 27 (19) (2015) 194129\.
  * Brogioli [2009] D. Brogioli, Extracting renewable energy from a salinity difference using a capacitor, Phys. Rev. Lett. 103 (5) (2009) 058501\.
  * Rica et al. [2012] R. Rica, R. Ziano, D. Salerno, F. Mantegazza, D. Brogioli, Thermodynamic relation between voltage-concentration dependence and salt adsorption in electrochemical cells, Phys. Rev. Lett. 109 (15) (2012) 156103\.
  * Simoncelli et al. [2018] M. Simoncelli, N. Ganfoud, A. Sene, M. Haefele, B. Daffos, P.-L. Taberna, M. Salanne, P. Simon, B. Rotenberg, Blue energy and desalination with nanoporous carbon electrodes: Capacitance from molecular simulations to continuous models, Phys. Rev. X 8 (2) (2018) 021024\.
  * Biesheuvel et al. [2011] P. Biesheuvel, Y. Fu, M. Z. Bazant, Diffuse charge and Faradaic reactions in porous electrodes, Phys. Rev. E 83 (6) (2011) 061507\.
  * Allagui et al. [2020a] A. Allagui, H. Alnaqbi, A. S. Elwakil, Z. Said, A. Hachicha, C. Wang, M. A. Abdelkareem, Fractional-order electric double-layer capacitors with tunable low-frequency impedance phase angle and energy storage capabilities, Appl. Phys. Lett. 116 (2020a) 013902\.
  * Allagui et al. [2016] A. Allagui, T. J. Freeborn, A. S. Elwakil, B. J. Maundy, Reevaluation of Performance of Electric Double-layer Capacitors from Constant-current Charge/Discharge and Cyclic Voltammetry, Sci. Rep. 6 (38568).
  * Allagui et al. [2018a] A. Allagui, D. Zhang, A. S. Elwakil, Short-Term Memory in Electric Double-Layer Capacitors, Appl. Phys. Lett. 113 (2018a) 253901–5.
  * Allagui et al. [2020b] A. Allagui, D. Zhang, I. Khakpour, A. S. Elwakil, C. Wang, Quantification of memory in fractional-order capacitors, J. Phys. D 53 (02LT03).
  * Hansen and Löwen [2000] J.-P. Hansen, H. Löwen, Effective interactions between electric double layers, Annu. Rev. Phys. Chem. 51 (1) (2000) 209–242.
  * Brown et al. [2016] M. A. Brown, Z. Abbas, A. Kleibert, R. G. Green, A. Goel, S. May, T. M. Squires, Determination of surface potential and electrical double-layer structure at the aqueous electrolyte-nanoparticle interface, Phys. Rev. X 6 (1) (2016) 011007.
  * Outhwaite [1974] C. Outhwaite, Higher-order closures and potential problems in diffuse double layer and strong electrolyte solution theory, Mol. Phys. 27 (3) (1974) 561–575.
  * Uralcan et al. [2016] B. Uralcan, I. A. Aksay, P. G. Debenedetti, D. T. Limmer, Concentration fluctuations and capacitive response in dense ionic solutions, J. Phys. Chem. Lett. 7 (13) (2016) 2333–2338.
  * Limmer et al. [2013] D. T. Limmer, C. Merlet, M. Salanne, D. Chandler, P. A. Madden, R. Van Roij, B. Rotenberg, Charge fluctuations in nanoscale capacitors, Phys. Rev. Lett. 111 (10) (2013) 106102\.
  * Härtel [2017] A. Härtel, Structure of electric double layers in capacitive systems and to what extent (classical) density functional theory describes it, J. Phys.: Condens. Matter 29 (42) (2017) 423002\.
  * Oldham [2008] K. B. Oldham, A Gouy–Chapman–Stern model of the double layer at a (metal)/(ionic liquid) interface, J. Electroanal. Chem. 613 (2) (2008) 131 – 138.
  * Gray and Stiles [2018] C. Gray, P. Stiles, Nonlinear electrostatics: the Poisson–Boltzmann equation, Eur. J. Phys. 39 (5) (2018) 053002.
  * Corry et al. [2000] B. Corry, S. Kuyucak, S.-H. Chung, Tests of continuum theories as models of ion channels. II. Poisson–Nernst–Planck theory versus Brownian dynamics, Biophys. J. 78 (5) (2000) 2364–2381.
  * Moreira and Netz [2000] A. G. Moreira, R. R. Netz, Strong-coupling theory for counter-ion distributions, Europhys. Lett. 52 (6) (2000) 705.
  * García-Morales and Krischer [2011] V. García-Morales, K. Krischer, Superstatistics in nanoscale electrochemical systems, Proc. Natl. Acad. Sci. 108 (49) (2011) 19535–19539.
  * Tsallis [2019] C. Tsallis, Beyond Boltzmann–Gibbs–Shannon in Physics and Elsewhere, Entropy 21 (7) (2019) 696\.
  * Garcia-Morales et al. [2004] V. Garcia-Morales, J. Cervera, J. Pellicer, Coupling theory for counterion distributions based in Tsallis statistics, Physica A 339 (3) (2004) 482 – 490, ISSN 0378-4371.
  * Tsallis [1988] C. Tsallis, Possible generalization of Boltzmann-Gibbs statistics, J. Stat. Phys. 52 (1-2) (1988) 479–487.
  * Beck and Cohen [2003] C. Beck, E. G. Cohen, Superstatistics, Physica A 322 (2003) 267–275.
  * Beck [2004] C. Beck, Superstatistics: theory and applications, Continuum Mech. Thermodyn. 16 (3) (2004) 293–304.
  * Abe et al. [2007] S. Abe, C. Beck, E. G. Cohen, Superstatistics, thermodynamics, and fluctuations, Phys. Rev. E 76 (3) (2007) 031102.
  * Bazant et al. [2009] M. Z. Bazant, M. S. Kilic, B. D. Storey, A. Ajdari, Towards an understanding of induced-charge electrokinetics at large applied voltages in concentrated solutions, Adv. Colloid Interface Sci. 152 (1-2) (2009) 48–88.
  * Biesheuvel and Van Soestbergen [2007] P. Biesheuvel, M. Van Soestbergen, Counterion volume effects in mixed electrical double layers, J. Colloid Interface Sci. 316 (2) (2007) 490–499.
  * Kornyshev [2007] A. A. Kornyshev, Double-layer in ionic liquids: paradigm change?, J. Phys. Chem. B 111 (20) (2007) 5545–5557.
  * Boukamp [1995] B. A. Boukamp, A linear Kronig-Kramers transform test for immittance data validation, J. Electrochem. Soc. 142 (6) (1995) 1885–1894.
  * Hasegawa [2007] H. Hasegawa, Non-extensive thermodynamics of transition-metal nanoclusters, Prog. Mater Sci. 52 (2-3) (2007) 333–351.
  * Tsallis [2004] C. Tsallis, What should a statistical mechanics satisfy to reflect nature?, Physica D 193 (1-4) (2004) 3–34.
  * García-Morales and Krischer [2010] V. García-Morales, K. Krischer, Fluctuation enhanced electrochemical reaction rates at the nanoscale, Proc. Natl. Acad. Sci. 107 (10) (2010) 4528–4532.
  * Bonnet et al. [2012] N. Bonnet, T. Morishita, O. Sugino, M. Otani, First-principles molecular dynamics at a constant electrode potential, Phys. Rev. Lett. 109 (26) (2012) 266101\.
  * Gell-Mann and Tsallis [2004] M. Gell-Mann, C. Tsallis, Nonextensive entropy: interdisciplinary applications, Oxford University Press on Demand, 2004.
  * Abe and Okamoto [2001] S. Abe, Y. Okamoto, Nonextensive statistical mechanics and its applications, vol. 560, Springer Science & Business Media, 2001.
  * Hanel et al. [2011] R. Hanel, S. Thurner, M. Gell-Mann, Generalized entropies and the transformation group of superstatistics, Proc. Natl. Acad. Sci. 108 (16) (2011) 6390–6394.
  * Tsallis [1999] C. Tsallis, Nonextensive statistics: theoretical, experimental and computational evidences and connections, Braz. J. Phys. 29 (1) (1999) 1–35.
  * Borges [1998] E. P. Borges, On a q-generalization of circular and hyperbolic functions, J. Phys. A: Math. Gen. 31 (23) (1998) 5281\.
  * Allagui et al. [2018b] A. Allagui, T. J. Freeborn, A. S. Elwakil, M. E. Fouda, B. J. Maundy, A. G. Radwanh, Z. Said, M. A. Abdelkareem, Review of Fractional-Order Electrical Characterization of Supercapacitors, J. Power Sources 400 (457–467).
  * Fouda et al. [2016] M. E. Fouda, A. S. Elwakil, A. G. Radwan, A. Allagui, Power and Energy Analysis of Fractional-order Electrical Energy Storage Devices, Energy 111 (2016) 785–792.
  * Zhang et al. [2019] D. Zhang, A. Allagui, A. S. Elwakil, A. M. Nassef, H. Rezk, J. Chengi, W. C. Choy, On the Modeling of Dispersive Transient Photocurrent Response of Organic Solar Cells, Org. Electron. 70 (2019) 42–47.

Figure 1: Schematic representation of electric double-layer structure at a
positively-polarized blocking electrode/electrolyte interface. Its thickness
(a measure of its capacity) is a usually in the order of a few nanometers at
most.

Figure 2: Plots of double-layer capacitance (Eq. 10) as a function of the
potential across the double layer for different values of $C_{x}^{0}$ with
$C_{\text{H}}=\text{28}\,\mu\text{F\,cm}^{-2}$,
$\epsilon=\text{80}\,\epsilon_{0}$, $z=\text{1}$, $T=\text{298}$ K

Figure 3: Plots of ion concentration (Eq. 14) as a function of the potential
across the double layer for different values of $q$ and with
$C_{x}^{0}=1\,\text{mol\,l}^{-1}$, $z=\text{1}$, $T=\text{298}$ K. The
straight lines with the slope of ${-zF}/{RT\ln(10)}$ are for $q=1.0$ only,
otherwise the decay-plots follow power-law profiles.

Figure 4: Plots of double-layer capacitance (Eq. 23) as a function of the
potential across the double layer for different values of $q$ (from -1 to 1)
with $C_{x}^{0}=10^{-3}\,\text{mol\,l}^{-1}$,
$C_{\text{H}}=\text{28}\,\mu\text{F\,cm}^{-2}$,
$\epsilon=\text{80}\,\epsilon_{0}$, $z=\text{1}$ and $T=\text{298}$ K


uncertainty on the $\sigma$ prediction,
$\sigma_{\mathrm{aN}^{3}\mathrm{LO}}^{(2)}$ is the NNLO hard cross section
with aN3LO PDFs and $\sigma_{\mathrm{NNLO}}^{(2)}$ is the full NNLO result. A
caveat to this treatment is that the theory uncertainty is sensitive to
unmatched cancellations and should therefore be used with care (and caution),
therefore the NNLO set remains the default in evaluating PDF uncertainties.

##  11 Conclusions

In this paper we have presented the first approximate N3LO global PDF fit.
This follows the MSHT20 framework [3], where the aN3LO PDF set also
incorporates estimates for theoretical uncertainties from missing N3LO
contributions and implicitly some MHOU beyond this. In addition, the framework
presented for obtaining these PDFs provides a means of utilising higher order
information as and when it is available. In contrast, previously, complete
information of the next order was required for theoretical calculations in PDF
fits. This provides a significant advantage moving forward in precision
phenomenology, since as we move to higher orders, this information takes
increasingly longer to calculate. We have analysed the resulting set of PDFs,
denoted MSHT20aN3LO, and made two sets available as described in Section 10.
The aN3LO PDF fits have been performed to the same set of global hard
scattering data and PDF parameterisations included for the MSHT20 NNLO PDF
fits.

The NNLO theoretical framework for MSHT20 PDFs has been extended in Section 2
to include the addition of general N3LO theory parameters into the fit.
Subsequently, we have outlined how these N3LO theory parameters can be
included into the Hessian procedure as controllable nuisance parameters where
they are not yet known. Two methods of handling subsets of the N3LO theory
parameters in the Hessian matrix have then been discussed; i.e. including or
ignoring correlations with aN3LO $K$-factors across distinct processes.

In Section 4 to Section 7 we have presented the N3LO additions to the relevant
splitting functions, transition matrix elements, heavy coefficient functions
and $K$-factors. We present usable and computationally efficient
approximations to N3LO based on known information in the small and large-$x$
regimes and the available Mellin moments (and make these available as
described in Section 10). In all cases the best fit prediction for each N3LO
function is in good agreement with the prior expected behaviour. Also in
Section 7, we find good agreement with recent progress towards N3LO DY and top
production $K$-factors [64, 66]. As more information becomes available
surrounding each of these functions, the framework we present here can be
easily adapted, aiding in the reduction in sources of MHOUs from N3LO. As we
have stressed, we interpret our theoretical uncertainty as being mainly due to
the remaining uncertainty at N3LO, but with some small, but significant
contribution from even higher orders, particularly at small-$x$. Our results
seem consistent with this interpretation. However, in the future we expect the
N3LO description to become more exact. Hence, at some point the remaining N3LO
uncertainty will become comparable to, or smaller than effects beyond N3LO. We
would then have to modify our procedure. However, we expect that once the N3LO
theory becomes very largely known, there will at this point also be more
information known about even higher orders (i.e. N4LO), which could then be
incorporated in a similar manner to maintain an estimate for MHOUs.
Alternatively, in the event that the available information is not suitable to
provide approximations (or indeed to complement these approximations), a
treatment similar in principle, but more sophisticated in practice, to that of
the $K$-factors may be adopted for DIS quantities. On this note, we
acknowledge that the method of constructing aN3LO $K$-factors for
non-(inclusive) DIS processes presented here is a first step towards a more
robust and flexible procedure, which is left for future work.

Combining together all N3LO information, in Section 8 the results of an
approximate N3LO global PDF fit are presented. The new MSHT20 approximate N3LO
PDFs show a significant reduction in $\chi^{2}$ from the MSHT20 NNLO PDF set,
with the leading NNLO tensions between HERA and non-HERA datasets heavily
reduced at aN3LO (most notably with the ATLAS 8 TeV $Z\ p_{T}$ dataset [110]).
With this being said, the aN3LO set does fit selected Jets datasets worse in
an aN3LO global fit than at NNLO, although these are an exception to the
behaviour seen for the other datasets. In performing a fit not including ATLAS
8 TeV $Z\ p_{T}$ data we provide evidence that similar tensions seen at NNLO
(see [3]) remain between this dataset and jet production data at aN3LO.
Further to this, we show that since HERA and ATLAS 8 TeV $Z\ p_{T}$ data are
more in agreement in the form of the high-$x$ gluon at aN3LO, one can observe
that the tension with the jet production data is shared between HERA and ATLAS
8 TeV $Z\ p_{T}$ data. Finally, as discussed, we highlight that in future work
it will be interesting to observe if this increased tension may be alleviated
when considering these jet datasets instead as dijet cross sections.

Investigating the correlations present within an aN3LO PDF fit, a natural
separation between process independent and process dependent parameters can be
observed. With this motivation, a PDF set with decorrelated aN3LO $K$-factor
eigenvectors is constructed. The validity of this is then also verified by
comparison with a second PDF set which includes correlations between all
parameters. Each of these sets exhibits similarly well behaved eigenvectors
and levels of dynamical tolerance.

Considering the form of the individual PDFs, the aN3LO PDFs include a much
harder gluon at small-$x$ due to contributions from the splitting functions as
discussed in Section 8.8. This enhancement then translates into an increase in
the charm and bottom PDFs due to the gluon input into the heavy flavour sector
via the transition matrix elements. At very low-$Q^{2}$ the result of the N3LO
additions is a non-negative charm and gluon PDF at small-$x$. As a consistency
check, the fit dependence on $\alpha_{s}$ and $m_{c}$ has been investigated.
In both of these cases we show a preference for values which suppress the
heavy flavour contributions (slightly lower $\alpha_{s}$ and slightly higher
$m_{c}$ than NNLO). Considering the predicted aN3LO $\alpha_{s}$, we observe a
slightly lower than $1\sigma$ effect when comparing with the NNLO world
average. While an extensive analysis of the aN3LO $\alpha_{s}$ value is left
for further study, since the world average is determined by NNLO results, one
could expect a small systematic effect from moving to N3LO.

Taking this analysis further and using the approximate N3LO PDFs as input to
N3LO cross section calculations, we consider the cases of gluon and vector
boson fusion in Higgs production. We present the first aN3LO calculation for
these cross sections and show how the aN3LO prediction differs from the case
with NNLO PDFs including scale variations, highlighting the importance of
matching orders in calculations. In VBF we provide an example where
cancellation is not realised between orders. However in this case the quark
sector is much more constrained and due to the smaller variation between
orders, there is naturally less scope for cancellation.

In summary, we have presented a set of approximate N3LO PDFs that are able to
more accurately predict physical quantities involving PDFs (given that all
ingredients in these calculations are included at N3LO or aN3LO). In producing
these PDFs, we have provided a more controllable method for estimating
theoretical uncertainties from MHOs in a PDF fit than scale variations. While
some ambiguity remains in this method in how the prior variations are chosen,
we argue that the current knowledge and intuition surrounding each source of
uncertainty can be utilised as and when available. This is therefore much more
in line with what one can expect a theoretical uncertainty to encompass.
Another potential shortcoming is the possibility of fitting to sources of
uncertainty other than higher orders (or higher order corrections elsewhere in
theory calculations included in a PDF fit). Although this is a possibility,
the position of the considered sources of uncertainty in the underlying theory
combined with the prior variations and penalties should act to minimise this
effect. In any case, if a separate source of uncertainty is significantly
affecting the fit, this will present itself as a source of tension with the
N3LO penalties and the $\chi^{2}$ (and PDF uncertainty) will be adapted
accordingly.

In future work it will be interesting to investigate the effects in the
high-$x$ gluon, which is a region of phenomenological importance and where the
interpretation of LHC constraints is not always straightforward. We also note
that there are N3LO results available from di-lepton rapidity in DY processes
[64]. Considering the results in Section 7 which display an agreement with
these recent results, we hope that these approximate N3LO PDFs may be of
interest in this analysis. Similarly for recent results considering top
production [66]. Furthermore, any approximate information from these results
could be included in the N3LO $K$-factor priors, which was not done for this
iteration of the aN3LO PDFs. Finally, in order to continually improve the
description of aN3LO PDFs, the inclusion of more sub-leading sources of MHOUs
could be addressed. With the upcoming wealth of experimental data from future
colliders such as the HL-LHC and the EIC, it will be of interest to gain a
better understanding of the transition matrix elements and also describe
better the charged current and longitudinal structure functions, where
currently theoretical uncertainties are much smaller than the experimental
uncertainties.

## Acknowledgements

J.M. thanks the Science and Technology Facilities Council (STFC) part of U.K.
Research and Innovation for support via Ph.D. funding. T.C. and R.S.T. thank
STFC for support via grant awards ST/P000274/1 and ST/T000856/1. L.H.L. thanks
STFC for support via grant awards ST/L000377/1 and ST/T000864/1. We would like
to thank Xuan Chen, Thomas Gehrmann, Nigel Glover and Alex Huss for providing
details of N3LO calculations. We would like to thank members of the PDF4LHC
working group for numerous discussions on PDFs and theoretical uncertainties.
We would also like to thank Alan Martin for long collaboration on the MSHT
series of PDFs.

##  Appendix A List of N3LO Ingredients

N3LO No. of Moments Small-$x$ Large-$x$ Function Moments (Even only)
$P_{qq}^{\mathrm{NS}}$ 8 $N=2-16$ [21] [21] [21] $P_{qq}^{\mathrm{PS}}$ 4
$N=2-8$ [35, 36] LL [28] N/A $P_{qg}$ 4 $N=2-8$ [35, 36] LL [28] N/A $P_{gq}$
4 $N=2-8$ [35, 36] LL [29, 30, 31] N/A $P_{gg}$ 4 $N=2-8$ [35, 36] LL & NLL
[29, 30, 31, 32, 33] N/A $A_{qq,H}^{\mathrm{NS}}$ 7 $N=2-14$ [50] N/A N/A
$A_{Hq}^{\mathrm{PS}}$ 6 $N=2-12$ [50] [53] [53] $A_{Hg}$ 5 $N=2-10$ [50] LL
[49] N/A $A_{gq,H}$ 7 $N=2-14$ [50] [54] [54] $A_{gg,H}$ 5 $N=2-10$ [50] N/A
N/A

Table A.1: List of all the N3LO ingredients used to construct the approximate
N3LO splitting functions and transition matrix elements. Where only a citation
is provided, extensive knowledge i.e. beyond NLL is used. This table is a non-
exhaustive list of the current knowledge about these functions, however
information beyond that which is provided here is not currently in a usable
format for phenomological studies.

GM-VFNS N3LO Known N3LO Components Function $C_{H,q}$ $C_{H,q}^{(3),\
\mathrm{FF}}\left(Q^{2}\leq m_{h}^{2}\right)$ LL [47, 48, 49],
$C_{H,q}^{\mathrm{VF},\ (3)}$ [57] $C_{H,g}$ $C_{H,g}^{(3),\
\mathrm{FF}}\left(Q^{2}\leq m_{h}^{2}\right)$ LL[47, 48, 49],
$C_{H,q}^{\mathrm{ZM},\ (3)}$ [57] $C_{q,q}^{\mathrm{NS}}$ $C_{q,q,\
\mathrm{NS}}^{\mathrm{ZM},\ (3)}$ [57] $C_{q,q}^{\mathrm{PS}}$ $C_{q,q,\
\mathrm{PS}}^{\mathrm{ZM},\ (3)}$ [57] $C_{q,g}$ $C_{q,g}^{\mathrm{ZM},\ (3)}$
[57]

Table A.2: List of all N3LO ingredients used to construct the approximate N3LO
GM-VFNS coefficient functions. Note that lower order components that
contribute to these functions are also known and are cited in the text. This
table only considers contributing 3-loop functions.

Table’s A.1 and A.2 summarise the available (at the time of writing) and used
information regarding the N3LO splitting functions and coefficient functions
respectively. The formalism presented in Section 2 currently makes use of all
this information and is able to be adapted as and when more information
becomes available.

##  Appendix B $\chi^{2}$ Results without HERA

###  B.1 NNLO

Table B.1 shows the differences in $\chi^{2}$ found when omitting HERA data
from a PDF fit using the MSHT NNLO PDFs. This table is copied here from [3]
for the ease of the reader. We see similarities between these results and the
$\Delta\chi^{2}$’s seen in the case of N3LO PDFs. Specifically the ATLAS $8\
\text{TeV}\ Z\ p_{T}$ displaying a substantial reduction from the global NNLO
fit. This therefore provides evidence that the inclusion of the N3LO
contributions is aiding in reducing tensions between the HERA and non-HERA
datasets.

Dataset $N_{\mathrm{pts}}$ $\chi^{2}$ $\Delta\chi^{2}$ BCDMS $\mu p$ $F_{2}$
[135] 163 174.7 $-5.5$ BCDMS $\mu d$ $F_{2}$ [135] 151 143.9 $-2.1$ NMC $\mu
p$ $F_{2}$ [136] 123 119.6 $-4.5$ NMC $\mu d$ $F_{2}$ [136] 123 96.6 $-16.1$
SLAC $ep$ $F_{2}$ [137, 138] 37 33.0 $+0.9$ SLAC $ed$ $F_{2}$ [137, 138] 38
24.1 $+1.1$ E665 $\mu d$ $F_{2}$ [139] 53 63.5 $+3.9$ E665 $\mu p$ $F_{2}$
[139] 53 68.9 $+4.3$ NuTeV $\nu N$ $F_{2}$ [140] 53 38.0 $-0.3$ NuTeV $\nu N$
$xF_{3}$ [140] 42 27.5 $-3.2$ NMC $\mu n/\mu p$ [141] 148 132.7 $+1.9$ E866 /
NuSea $pp$ DY [80] 184 228.0 $+2.9$ E866 / NuSea $pd/pp$ DY [81] 15 9.1 $-1.3$
CCFR $\nu N\rightarrow\mu\mu X$ [133] 86 66.2 $-1.5$ NuTeV $\nu
N\rightarrow\mu\mu X$ [133] 84 49.0 $-9.5$ CHORUS $\nu N$ $F_{2}$ [146] 42
29.6 $-0.6$ CHORUS $\nu N$ $xF_{3}$ [146] 28 18.2 $-0.3$ CDF II $p\bar{p}$
incl. jets [102] 76 60.9 $+0.5$ DØ II $Z$ rap. [82] 28 16.6 $+0.3$ CDF II $Z$
rap. [83] 28 38.7 $+1.5$ DØ II $W\rightarrow\nu\mu$ asym. [84] 10 17.4 $+0.1$
CDF II $W$ asym. [85] 13 19.0 $+0.0$ DØ II $W\rightarrow\nu e$ asym. [86] 12
30.0 $-3.9$ DØ II $p\bar{p}$ incl. jets [103] 110 119.3 $-0.9$ ATLAS $W^{+},\
W^{-},\ Z$ [87] 30 29.5 $-0.4$

Table B.1: The change in $\chi^{2}$ for a NNLO fit(with negative indicating an
improvement in the fit quality) when the combined HERA data sets including
$F_{L}$ and heavy flavour data are removed, illustrating the tensions of these
data sets with several of the other data sets in the global fit.
$\Delta\chi^{2}$ represents the change from a full global fit at the same
order in $\alpha_{s}$.

Dataset $N_{\mathrm{pts}}$ $\chi^{2}$ $\Delta\chi^{2}$ CMS W asym. $p_{T}>35\
\text{GeV}$ [88] 11 6.6 $-1.2$ CMS W asym. $p_{T}>25,30\ \text{GeV}$ [89] 24
7.5 $+0.1$ LHCb $Z\rightarrow e^{+}e^{-}$ [90] 9 24.2 $+1.5$ LHCb W asym.
$p_{T}>20\ \text{GeV}$ [91] 10 12.1 $-0.3$ CMS $Z\rightarrow e^{+}e^{-}$ [92]
35 17.3 $-0.6$ ATLAS High-mass Drell-Yan [93] 13 16.9 $-2.0$ Tevatron, ATLAS,
CMS $\sigma_{t\bar{t}}$ [117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
127, 128, 129] 17 14.2 $-0.4$ CMS double diff. Drell-Yan [94] 132 134.2
$-10.3$ LHCb 2015 $W,Z$ [77, 78] 67 97.4 $-1.9$ LHCb $8\ \text{TeV}$
$Z\rightarrow ee$ [95] 17 24.4 $-1.8$ CMS $8\ \text{TeV}\ W$ [96] 22 13.7
$+0.9$ ATLAS $7\ \text{TeV}$ jets [104] 140 228.0 $+6.5$ CMS $7\ \text{TeV}\
W+c$ [62] 10 9.2 $+0.6$ ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] 61 116.8
$+0.2$ CMS $7\ \text{TeV}$ jets [101] 158 179.5 $+3.8$ DØ $W$ asym. [97] 14
11.3 $-0.8$ ATLAS $8\ \text{TeV}\ Z\ p_{T}$ [110] 104 149.3 $-39.2$ CMS $8\
\text{TeV}$ jets [105] 174 259.5 $-1.8$ ATLAS $8\ \text{TeV}$ sing. diff.
$t\bar{t}$ [130] 25 24.5 $-1.1$ ATLAS $8\ \text{TeV}$ sing. diff. $t\bar{t}$
dilep. [131] 5 2.3 $-1.1$ ATLAS $8\ \text{TeV}$ High-mass DY [98] 48 60.9
$+3.7$ ATLAS $8\ \text{TeV}\ W+\text{jets}$ [111] 30 16.4 $-1.7$ CMS $8\
\text{TeV}$ double diff. $t\bar{t}$ [132] 15 23.3 $+0.8$ ATLAS $8\ \text{TeV}\
W$ [99] 22 54.4 $-3.0$ CMS $2.76\ \text{TeV}$ jet [106] 81 102.9 $+0.0$ CMS
$8\ \text{TeV}$ sing. diff. $t\bar{t}$ [112] 9 10.6 $-2.6$ ATLAS $8\
\text{TeV}$ double diff. $Z$ [100] 59 108.3 $+22.7$ Total 3042 3379.6 $-61.6$

Table B.2: (Continued) The change in $\chi^{2}$ (with negative indicating an
improvement in the fit quality) when the combined HERA data sets including
$F_{L}$ and heavy flavour data are removed, illustrating the tensions of these
data sets with several of the other data sets in the global fit.
$\Delta\chi^{2}$ represents the change from a full global fit at the same
order in $\alpha_{s}$.

###  B.2 aN3LO

Table B.3 shows the differences in $\chi^{2}$ found when omitting HERA data
from a PDF fit using the MSHT aN3LO PDFs. These results show that at aN3LO the
fit no longer experiences large tensions between HERA and ATLAS $8\
\text{TeV}\ Z\ p_{T}$ [110] datasets. The main tensions at N3LO are now
concerning the Jets data with HERA (and most likely some non-HERA datasets).
This result is not unexpected due to the known issues surrounding jets
especially as we move to higher precision [109].

Dataset $N_{\mathrm{pts}}$ $\chi^{2}$ $\Delta\chi^{2}$ BCDMS $\mu p$ $F_{2}$
[135] 163 175.8 $+1.4$ BCDMS $\mu d$ $F_{2}$ [135] 151 144.2 $-0.0$ NMC $\mu
p$ $F_{2}$ [136] 123 113.6 $-7.8$ NMC $\mu d$ $F_{2}$ [136] 123 87.6 $-16.6$
SLAC $ep$ $F_{2}$ [137, 138] 37 30.7 $-0.9$ SLAC $ed$ $F_{2}$ [137, 138] 38
23.2 $+0.4$ E665 $\mu d$ $F_{2}$ [139] 53 65.2 $+1.3$ E665 $\mu p$ $F_{2}$
[139] 53 69.0 $+1.5$ NuTeV $\nu N$ $F_{2}$ [140] 53 35.4 $-0.4$ NuTeV $\nu N$
$xF_{3}$ [140] 42 29.2 $-5.6$ NMC $\mu n/\mu p$ [141] 148 131.1 $-0.5$ E866 /
NuSea $pp$ DY [80] 184 225.4 $+2.1$ E866 / NuSea $pd/pp$ DY [81] 15 8.2 $-0.2$
CCFR $\nu N\rightarrow\mu\mu X$ [133] 86 67.0 $-1.3$ NuTeV $\nu
N\rightarrow\mu\mu X$ [133] 84 47.6 $-9.1$ CHORUS $\nu N$ $F_{2}$ [146] 42
29.0 $-0.2$ CHORUS $\nu N$ $xF_{3}$ [146] 28 18.5 $+0.4$ CDF II $p\bar{p}$
incl. jets [102] 76 65.9 $-0.6$ DØ II $Z$ rap. [82] 28 17.6 $+0.3$ CDF II $Z$
rap. [83] 28 42.0 $+1.5$ DØ II $W\rightarrow\nu\mu$ asym. [84] 10 18.9 $+2.4$
CDF II $W$ asym. [85] 13 19.2 $+1.0$ DØ II $W\rightarrow\nu e$ asym. [86] 12
31.0 $+0.3$ DØ II $p\bar{p}$ incl. jets [103] 110 114.2 $+0.9$ ATLAS $W^{+},\
W^{-},\ Z$ [87] 30 29.4 $-0.6$

Table B.3: The change in $\chi^{2}$ for an N3LO fit (with negative indicating
an improvement in the fit quality) when the combined HERA data sets including
$F_{L}$ and heavy flavour data are removed, illustrating the tensions of these
data sets with several of the other data sets in the global fit.
$\Delta\chi^{2}$ represents the change from a full global fit at the same
order in $\alpha_{s}$.

Dataset $N_{\mathrm{pts}}$ $\chi^{2}$ $\Delta\chi^{2}$ CMS W asym. $p_{T}>35\
\text{GeV}$ [88] 11 7.0 $+0.3$ CMS W asym. $p_{T}>25,30\ \text{GeV}$ [89] 24
7.6 $-0.1$ LHCb $Z\rightarrow e^{+}e^{-}$ [90] 9 22.4 $-1.7$ LHCb W asym.
$p_{T}>20\ \text{GeV}$ [91] 10 12.6 $+0.0$ CMS $Z\rightarrow e^{+}e^{-}$ [92]
35 17.1 $-0.4$ ATLAS High-mass Drell-Yan [93] 13 17.2 $-0.9$ Tevatron, ATLAS,
CMS $\sigma_{t\bar{t}}$ [117, 118, 119, 120, 121, 122, 123, 124, 125, 126,
127, 128, 129] 17 14.4 $+0.1$ CMS double diff. Drell-Yan [94] 132 126.0 $-3.5$
LHCb 2015 $W,Z$ [77, 78] 67 97.6 $-6.3$ LHCb $8\text{TeV}$ $Z\rightarrow ee$
[95] 17 26.6 $-2.2$ CMS $8\ \text{TeV}\ W$ [96] 22 11.9 $+0.1$ ATLAS $7\
\text{TeV}$ jets [104] 140 217.7 $+1.8$ CMS $7\ \text{TeV}\ W+c$ [62] 10 10.8
$+0.0$ ATLAS $7\ \text{TeV}$ high prec. $W,Z$ [79] 61 94.3 $+0.2$ CMS $7\
\text{TeV}$ jets [101] 158 187.8 $+1.0$ DØ $W$ asym. [97] 14 10.1 $-2.1$ ATLAS
$8\ \text{TeV}\ Z\ p_{T}$ [110] 104 121.2 $+12.8$ CMS $8\ \text{TeV}$ jets
[105] 174 259.8 $-11.5$ ATLAS $8\ \text{TeV}$ sing. diff. $t\bar{t}$ [130] 25
24.1 $-0.2$ ATLAS $8\ \text{TeV}$ sing. diff. $t\bar{t}$ dilep. [131] 5 3.0
$+0.3$ ATLAS $8\ \text{TeV}$ High-mass DY [98] 48 65.0 $+2.0$ ATLAS $8\
\text{TeV}\ W+\text{jets}$ [111] 30 18.0 $-0.8$ CMS $8\ \text{TeV}$ double
diff. $t\bar{t}$ [132] 15 22.8 $-0.8$ ATLAS $8\ \text{TeV}\ W$ [99] 22 52.5
$-5.5$ CMS $2.76\ \text{TeV}$ jet [106] 81 103.0 $-6.8$ CMS $8\ \text{TeV}$
sing. diff. $t\bar{t}$ [112] 9 12.3 $+2.0$ ATLAS $8\ \text{TeV}$ double diff.
$Z$ [100] 59 95.5 $+3.9$

Table B.4: (Continued) The change in $\chi^{2}$ for an N3LO fit (with negative
indicating an improvement in the fit quality) when the combined HERA data sets
including $F_{L}$ and heavy flavour data are removed, illustrating the
tensions of these data sets with several of the other data sets in the global
fit. $\Delta\chi^{2}$ represents the change from a full global fit at the same
order in $\alpha_{s}$.

Low-$Q^{2}$ Coefficient $c_{q}^{\mathrm{NLL}}$ $=-3.844$ 0.006
$c_{g}^{\mathrm{NLL}}$ $=-3.875$ 0.004 Transition Matrix Elements $a_{Hg}$
$=17788.000$ 5.607 $a_{qq,H}^{\mathrm{NS}}$ $=-63.950$ 0.000 $a_{gg,H}$
$=-1334.500$ 0.001 Splitting Functions $\rho_{qq}^{NS}$ $=0.007$ 0.000
$\rho_{gq}$ $=-1.647$ 0.001 $\rho_{qq}^{PS}$ $=-0.579$ 0.429 $\rho_{gg}$
$=9.237$ 0.023 $\rho_{qg}$ $=-1.343$ 0.131 K-factors
$\mathrm{DY}_{\mathrm{NLO}}$ $=-0.279$ 0.080 $\mathrm{DY}_{\mathrm{NNLO}}$
$=0.072$ 0.005 $\mathrm{Top}_{\mathrm{NLO}}$ $=-0.204$ 0.042
$\mathrm{Top}_{\mathrm{NNLO}}$ $=0.412$ 0.170 $\mathrm{Jet}_{\mathrm{NLO}}$
$=-0.254$ 0.065 $\mathrm{Jet}_{\mathrm{NNLO}}$ $=-0.861$ 0.741
$p_{T}\mathrm{Jets}_{\mathrm{NLO}}$ $=0.461$ 0.213
$p_{T}\mathrm{Jets}_{\mathrm{NNLO}}$ $=0.016$ 0.000
$\mathrm{Dimuon}_{\mathrm{NLO}}$ $=-0.329$ 0.109
$\mathrm{Dimuon}_{\mathrm{NNLO}}$ $=0.587$ 0.345 Total 3306.8 / 3042
$\Delta\chi^{2}$ from N3LO $-49.0$

Table B.5: (Continued) The change in $\chi^{2}$ for an N3LO fit (with negative
indicating an improvement in the fit quality) when the combined HERA data sets
including $F_{L}$ and heavy flavour data are removed, illustrating the
tensions of these data sets with several of the other data sets in the global
fit. $\Delta\chi^{2}$ represents the change from a full global fit at the same
order in $\alpha_{s}$.

##  Appendix C Dynamic Tolerances

In this section we provide an exhaustive breakdown of the
$\Delta\chi^{2}_{\mathrm{global}}$ behaviour for all eigenvectors found where
N3LO $K$-factor parameters are considered completely decorrelated
($H_{ij}+K_{ij}$) or correlated ($H_{ij}^{\prime}$) with all other parameters.

###  C.1 Case 1: Decorrelated $K$-factor Parameters

Fig. 57 displays the tolerance landscape for each eigenvector found from the
decorrelated ($H_{ij}+K_{ij}$) Hessian described in Section 2. Across all 52
eigenvectors (42 PDF + N3LO DIS theory and 10 N3LO $K$-factor) we show an
overall general agreement with the quadratic assumption similar to that found
at NNLO.

Figure 57: Dynamic tolerances for each eigenvector direction in the case of
complete decorrelation between the theory and PDF parameters, and the
$K$-factor parameters included in the PDF fit.

Figure 58: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete decorrelation between the theory and PDF parameters, and
the $K$-factor parameters included in the PDF fit.

Figure 59: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete decorrelation between the theory and PDF parameters, and
the $K$-factor parameters included in the PDF fit.

Figure 60: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete decorrelation between the theory and PDF parameters, and
the $K$-factor parameters included in the PDF fit.

Figure 61: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete decorrelation between the theory and PDF parameters, and
the $K$-factor parameters included in the PDF fit.

Figure 62: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete decorrelation between the theory and PDF parameters, and
the $K$-factor parameters included in the PDF fit.

Figure 63: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete decorrelation between the theory and PDF parameters, and
the $K$-factor parameters included in the PDF fit.

###  C.2 Case 2: Correlated $K$-factor Parameters

Fig. 64 displays the tolerance landscape for each eigenvector found from the
correlated ($H_{ij}^{\prime}$) Hessian described in Section 2. Across all 52
eigenvectors we show an overall general agreement with the quadratic
assumption similar to that found at NNLO.

Figure 64: Dynamic tolerances for each eigenvector direction in the case of
complete correlation between all theory, PDF and $K$-factor parameters
included in the PDF fit.

Figure 65: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete correlation between all theory, PDF and $K$-factor
parameters included in the PDF fit.

Figure 66: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete correlation between all theory, PDF and $K$-factor
parameters included in the PDF fit.

Figure 67: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete correlation between all theory, PDF and $K$-factor
parameters included in the PDF fit.

Figure 68: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete correlation between all theory, PDF and $K$-factor
parameters included in the PDF fit.

Figure 69: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete correlation between all theory, PDF and $K$-factor
parameters included in the PDF fit.

Figure 70: (Continued) Dynamic tolerances for each eigenvector direction in
the case of complete correlation between all theory, PDF and $K$-factor
parameters included in the PDF fit.

##  Appendix D Higgs Gluon Fusion $\mu=m_{H}$ Results

$\sigma$ order PDF order $\sigma+\Delta\sigma_{+}-\Delta\sigma_{-}$ (pb)
$\sigma$ (pb) $+\ \Delta\sigma_{+}-\Delta\sigma_{-}$ (%) PDF uncertainties
N3LO aN3LO (no theory unc.) $43.803+0.685-0.526$ $42.709+2.81\%-3.14\%$ aN3LO
($H_{ij}+K_{ij}$) $43.803+0.795-0.732$ $42.709+3.30\%-3.17\%$ aN3LO
($H_{ij}^{\prime}$) $43.803+0.785-0.737$ $42.709+3.39\%-3.08\%$ NNLO
$46.243+0.524-0.563$ $46.243+1.13\%-1.22\%$ NNLO NNLO $42.129+0.472-0.510$
$42.129+1.12\%-1.21\%$ PDF + Scale uncertainties N3LO aN3LO (no theory unc.)
$43.803+1.723-2.519$ $42.709+4.68\%-6.44\%$ aN3LO ($H_{ij}+K_{ij}$)
$43.803+1.770-2.570$ $42.709+4.89\%-6.45\%$ aN3LO ($H_{ij}^{\prime}$)
$43.803+1.766-2.571$ $42.709+4.95\%-6.41\%$ NNLO $46.243+1.845-3.078$
$46.243+3.99\%-6.66\%$ NNLO NNLO $42.129+4.989-5.106$ $42.129+11.84\%-12.12\%$

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

Provided in Table D.1 are the results analogous to those in Table 20 but with
the central scale set to $\mu=\mu_{f}=\mu_{r}=m_{H}$. These results show a
higher level of stability for aN3LO PDFs with the chosen central scale. By the
renormalisation group arguments, this scale dependence should disappear at all
orders in perturbation theory. Therefore the results here suggest that the
aN3LO PDFs are following this trend.

## References

  * [1] M. Mojaza, S. J. Brodsky, and X.-G. Wu,  Phys. Rev. Lett. 110, 192001 (2013), 1212.0049.
  * [2] M. Czakon, D. Heymes, and A. Mitov,  JHEP 04, 071 (2017), 1606.03350.
  * [3] S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin, and R. S. Thorne,  Eur. Phys. J. C 81, 341 (2021), 2012.04684.
  * [4] T.-J. Hou et al.,  Phys. Rev. D 103, 014013 (2021), 1912.10053.
  * [5] NNPDF, R. D. Ball et al.,  Eur. Phys. J. C 77, 663 (2017), 1706.00428.
  * [6] NNPDF, R. D. Ball et al.,  Eur. Phys. J. C 82, 428 (2022), 2109.02653.
  * [7] NNPDF, R. D. Ball et al.,  Eur. Phys. J. C 81, 958 (2021), 2109.02671.
  * [8] PDF4LHC Working Group, R. D. Ball et al.,  J. Phys. G 49, 080501 (2022), 2203.05506.
  * [9] S. Alekhin, J. Blümlein, S. Moch, and R. Placakyte,  Phys. Rev. D 96, 014011 (2017), 1701.05838.
  * [10] ATLAS, G. Aad et al.,  Eur. Phys. J. C 82, 438 (2022), 2112.11266.
  * [11] NNPDF, R. Abdul Khalek et al.,  Eur. Phys. J. C 79, 931 (2019), 1906.10698.
  * [12] L. A. Harland-Lang and R. S. Thorne,  Eur. Phys. J. C 79, 225 (2019), 1811.08434.
  * [13] J. Gao, L. Harland-Lang, and J. Rojo,  Phys. Rept. 742, 1 (2018), 1709.04922.
  * [14] R. D. Ball and R. L. Pearson,  Eur. Phys. J. C 81, 830 (2021), 2105.05114.
  * [15] M. Bonvini,  Eur. Phys. J. C 80, 989 (2020), 2006.16293.
  * [16] R. D. Ball et al.,  Eur. Phys. J. C 78, 321 (2018), 1710.05935.
  * [17] xFitter Developers’ Team, H. Abdolmaleki et al.,  Eur. Phys. J. C 78, 621 (2018), 1802.00064.
  * [18] F. J. Tackmann,  Theory Uncertainties from Nuisance Parameters,  SCET 2019 Workshop, 2019.
  * [19] F. J. Tackmann,  Beyond Scale Variations: Perturbative Theory Uncertainties from Nuisance Parameters,  DESY-19-021, 2023.
  * [20] G. Altarelli and G. Parisi,  Nucl. Phys. B 126, 298 (1977).
  * [21] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren, and A. Vogt,  JHEP 10, 041 (2017), 1707.08315.
  * [22] E. G. Floratos, D. A. Ross, and C. T. Sachrajda,  Nucl. Phys. B 152, 493 (1979).
  * [23] A. Gonzalez-Arroyo and C. Lopez,  Nucl. Phys. B 166, 429 (1980).
  * [24] W. Furmanski and R. Petronzio,  Phys. Lett. B 97, 437 (1980).
  * [25] E. G. Floratos, C. Kounnas, and R. Lacaze,  Nucl. Phys. B 192, 417 (1981).
  * [26] S. Moch, J. A. M. Vermaseren, and A. Vogt,  Nucl. Phys. B 688, 101 (2004), hep-ph/0403192.
  * [27] A. Vogt, S. Moch, and J. A. M. Vermaseren,  Nucl. Phys. B 691, 129 (2004), hep-ph/0404111.
  * [28] S. Catani and F. Hautmann,  Nucl. Phys. B 427, 475 (1994), hep-ph/9405388.
  * [29] L. N. Lipatov,  Sov. J. Nucl. Phys. 23, 338 (1976).
  * [30] E. A. Kuraev, L. N. Lipatov, and V. S. Fadin,  Sov. Phys. JETP 45, 199 (1977).
  * [31] I. I. Balitsky and L. N. Lipatov,  Sov. J. Nucl. Phys. 28, 822 (1978).
  * [32] V. S. Fadin and L. N. Lipatov,  Phys. Lett. B 429, 127 (1998), hep-ph/9802290.
  * [33] M. Ciafaloni and G. Camici,  Phys. Lett. B 430, 349 (1998), hep-ph/9803389.
  * [34] M. Bonvini and S. Marzani,  JHEP 06, 145 (2018), 1805.06460.
  * [35] A. Vogt et al.,  PoS LL2018, 050 (2018), 1808.08981.
  * [36] S. Moch, B. Ruijl, T. Ueda, J. A. M. Vermaseren, and A. Vogt,  Phys. Lett. B 825, 136853 (2022), 2111.15561.
  * [37] R. K. Ellis, W. J. Stirling, and B. R. Webber,  QCD and Collider Physics, Vol. 8 (Cambridge University Press, 2011).
  * [38] M. Buza, Y. Matiounine, J. Smith, and W. L. van Neerven,  Eur. Phys. J. C 1, 301 (1998), hep-ph/9612398.
  * [39] M. Buza, Y. Matiounine, J. Smith, and W. L. van Neerven,  Nucl. Phys. B 485, 420 (1997), hep-ph/9608342.
  * [40] M. A. G. Aivazis, J. C. Collins, F. I. Olness, and W.-K. Tung,  Phys. Rev. D 50, 3102 (1994), hep-ph/9312319.
  * [41] R. S. Thorne and R. G. Roberts,  Phys. Rev. D 57, 6871 (1998), hep-ph/9709442.
  * [42] A. D. Martin, W. J. Stirling, R. S. Thorne, and G. Watt,  Eur. Phys. J. C 63, 189 (2009), 0901.0002.
  * [43] R. S. Thorne,  Phys. Rev. D 73, 054019 (2006), hep-ph/0601245.
  * [44] R. S. Thorne,  Phys. Rev. D 86, 074017 (2012), 1201.6180.
  * [45] E. Laenen, S. Riemersma, J. Smith, and W. L. van Neerven,  Nucl. Phys. B 392, 229 (1993).
  * [46] B. W. Harris and J. Smith,  Nucl. Phys. B 452, 109 (1995), hep-ph/9503484.
  * [47] S. Catani, M. Ciafaloni, and F. Hautmann,  Nucl. Phys. B 366, 135 (1991).
  * [48] E. Laenen and S.-O. Moch,  Phys. Rev. D 59, 034027 (1999), hep-ph/9809550.
  * [49] H. Kawamura, N. A. Lo Presti, S. Moch, and A. Vogt,  Nucl. Phys. B 864, 399 (2012), 1205.5727.
  * [50] I. Bierenbaum, J. Blumlein, and S. Klein,  Nucl. Phys. B 820, 417 (2009), 0904.3563.
  * [51] J. Ablinger et al.,  Nucl. Phys. B 886, 733 (2014), 1406.4654.
  * [52] J. Blümlein, P. Marquard, C. Schneider, and K. Schönwald,  Nucl. Phys. B 971, 115542 (2021), 2107.06267.
  * [53] J. Ablinger et al.,  Nucl. Phys. B 890, 48 (2014), 1409.1135.
  * [54] J. Ablinger et al.,  Nuclear Physics B 882, 263 (2014), 1402.0359.
  * [55] J. Blümlein et al.,  PoS QCDEV2017, 031 (2017), 1711.07957.
  * [56] M. A. G. Aivazis, F. I. Olness, and W.-K. Tung,  Phys. Rev. D 50, 3085 (1994), hep-ph/9312318.
  * [57] J. A. M. Vermaseren, A. Vogt, and S. Moch,  Nucl. Phys. B 724, 3 (2005), hep-ph/0504242.
  * [58] W. L. van Neerven and A. Vogt,  Nucl. Phys. B 568, 263 (2000), hep-ph/9907472.
  * [59] W. L. van Neerven and A. Vogt,  Nucl. Phys. B 588, 345 (2000), hep-ph/0006154.
  * [60] J. Ablinger et al.,  Nucl. Phys. B 885, 280 (2014), 1405.4259.
  * [61] J. Ablinger et al.,  Nucl. Part. Phys. Proc. 258-259, 37 (2015), 1409.1435.
  * [62] CMS, S. Chatrchyan et al.,  JHEP 02, 013 (2014), 1310.1138.
  * [63] C. Duhr, F. Dulat, and B. Mistlberger,  Phys. Rev. Lett. 125, 172001 (2020), 2001.07717.
  * [64] X. Chen et al.,  Phys. Rev. Lett. 128, 052001 (2022), 2107.09085.
  * [65] C. Duhr and B. Mistlberger,  JHEP 03, 116 (2022), 2111.10379.
  * [66] N. Kidonakis,  2022 Snowmass Summer Study (2022), 2203.03698.
  * [67] R. D. Ball, M. Bonvini, S. Forte, S. Marzani, and G. Ridolfi,  Nucl. Phys. B 874, 746 (2013), 1303.3590.
  * [68] M. Bonvini, R. D. Ball, S. Forte, S. Marzani, and G. Ridolfi,  J. Phys. G 41, 095002 (2014), 1404.3204.
  * [69] M. Bonvini, S. Marzani, C. Muselli, and L. Rottoli,  JHEP 08, 105 (2016), 1603.08000.
  * [70] T. Ahmed, M. C. Kumar, P. Mathews, N. Rana, and V. Ravindran,  Eur. Phys. J. C 76, 355 (2016), 1510.02235.
  * [71] M. Bonvini and S. Marzani,  Phys. Rev. Lett. 120, 202003 (2018), 1802.07758.
  * [72] M. Bonvini,  Eur. Phys. J. C 78, 834 (2018), 1805.08785.
  * [73] M. Bonvini,  EPJ Web Conf. 60, 12008 (2013), 1306.6633.
  * [74] C. Anastasiou et al.,  JHEP 03, 091 (2015), 1411.3584.
  * [75] C. Anastasiou et al.,  JHEP 05, 058 (2016), 1602.00695.
  * [76] B. Mistlberger,  JHEP 05, 028 (2018), 1802.00833.
  * [77] LHCb, R. Aaij et al.,  JHEP 08, 039 (2015), 1505.07024.
  * [78] LHCb, R. Aaij et al.,  JHEP 01, 155 (2016), 1511.08039.
  * [79] ATLAS, M. Aaboud et al.,  Eur. Phys. J. C 77, 367 (2017), 1612.03016.
  * [80] J. C. Webb,  Measurement of continuum dimuon production in 800-GeV/C proton nucleon collisions,  PhD thesis, New Mexico State U., 2003, hep-ex/0301031.
  * [81] NuSea, R. S. Towell et al.,  Phys. Rev. D64, 052002 (2001), hep-ex/0103030.
  * [82] D0, V. M. Abazov et al.,  Phys. Rev. D76, 012003 (2007), hep-ex/0702025.
  * [83] CDF, T. A. Aaltonen et al.,  Phys. Lett. B692, 232 (2010), 0908.3914.
  * [84] D0, V. M. Abazov et al.,  Phys. Rev. D88, 091102 (2013), 1309.2591.
  * [85] CDF, T. Aaltonen et al.,  Phys. Rev. Lett. 102, 181801 (2009), 0901.2169.
  * [86] D0, V. M. Abazov et al.,  Phys. Rev. Lett. 101, 211801 (2008), 0807.3367.
  * [87] ATLAS, G. Aad et al.,  Phys. Rev. D 85, 072004 (2012), 1109.5141.
  * [88] CMS, S. Chatrchyan et al.,  Phys. Rev. Lett. 109, 111806 (2012), 1206.2598.
  * [89] CMS, S. Chatrchyan et al.,  JHEP 04, 050 (2011), 1103.3470.
  * [90] LHCb, R. Aaij et al.,  JHEP 02, 106 (2013), 1212.4620.
  * [91] LHCb, R. Aaij et al.,  JHEP 06, 058 (2012), 1204.1620.
  * [92] CMS, S. Chatrchyan et al.,  Phys. Rev. D85, 032002 (2012), 1110.4973.
  * [93] ATLAS, G. Aad et al.,  Phys. Lett. B725, 223 (2013), 1305.4192.
  * [94] CMS, S. Chatrchyan et al.,  JHEP 12, 030 (2013), 1310.7291.
  * [95] LHCb, R. Aaij et al.,  JHEP 05, 109 (2015), 1503.00963.
  * [96] CMS, V. Khachatryan et al.,  Eur. Phys. J. C 76, 469 (2016), 1603.01803.
  * [97] D0, V. M. Abazov et al.,  Phys. Rev. Lett. 112, 151803 (2014), 1312.2895,  [Erratum: Phys.Rev.Lett. 114, 049901 (2015)].
  * [98] ATLAS, G. Aad et al.,  JHEP 08, 009 (2016), 1606.01736.
  * [99] ATLAS, G. Aad et al.,  Eur. Phys. J. C 79, 760 (2019), 1904.05631.
  * [100] ATLAS, M. Aaboud et al.,  JHEP 12, 059 (2017), 1710.05167.
  * [101] CMS, S. Chatrchyan et al.,  Phys. Rev. D 90, 072006 (2014), 1406.0324.
  * [102] CDF, A. Abulencia et al.,  Phys. Rev. D75, 092006 (2007), hep-ex/0701051,  [Erratum: Phys. Rev.D75,119901(2007)].
  * [103] D0, V. M. Abazov et al.,  Phys. Rev. D85, 052006 (2012), 1110.3771.
  * [104] ATLAS, G. Aad et al.,  JHEP 02, 153 (2015), 1410.8857,  [Erratum: JHEP 09, 141 (2015)].
  * [105] CMS, V. Khachatryan et al.,  JHEP 03, 156 (2017), 1609.05331.
  * [106] CMS, V. Khachatryan et al.,  Eur. Phys. J. C 76, 265 (2016), 1512.06212.
  * [107] L. A. Harland-Lang, A. D. Martin, and R. S. Thorne,  Eur. Phys. J. C 78, 248 (2018), 1711.05757.
  * [108] ATLAS, M. Aaboud et al.,  JHEP 09, 020 (2017), 1706.03192.
  * [109] R. Abdul Khalek et al.,  Eur. Phys. J. C 80, 797 (2020), 2005.11327.
  * [110] ATLAS, G. Aad et al.,  Eur. Phys. J. C 76, 291 (2016), 1512.02192.
  * [111] ATLAS, M. Aaboud et al.,  JHEP 05, 077 (2018), 1711.03296.
  * [112] CMS, V. Khachatryan et al.,  Eur. Phys. J. C 75, 542 (2015), 1505.04480.
  * [113] N. Kidonakis and G. F. Sterman,  Nucl. Phys. B 505, 321 (1997), hep-ph/9705234.
  * [114] N. Kidonakis,  Phys. Rev. Lett. 102, 232003 (2009), 0903.2561.
  * [115] N. Kidonakis,  Phys. Rev. D 82, 114030 (2010), 1009.4935.
  * [116] N. Kidonakis,  15${}^{\text{th}}$ International Symposium on Radiative Corrections (2021), 2109.14102.
  * [117] CDF, D0, T. A. Aaltonen et al.,  Phys. Rev. D89, 072001 (2014), 1309.7570.
  * [118] ATLAS, G. Aad et al.,  Eur. Phys. J. C71, 1577 (2011), 1012.1792.
  * [119] ATLAS, G. Aad et al.,  Phys. Lett. B707, 459 (2012), 1108.3699.
  * [120] ATLAS, G. Aad et al.,  Phys. Lett. B711, 244 (2012), 1201.1889.
  * [121] ATLAS, G. Aad et al.,  JHEP 05, 059 (2012), 1202.4892.
  * [122] ATLAS, G. Aad et al.,  Phys. Lett. B717, 89 (2012), 1205.2067.
  * [123] ATLAS, G. Aad et al.,  Eur. Phys. J. C73, 2328 (2013), 1211.7205.
  * [124] CMS, S. Chatrchyan et al.,  Phys. Rev. D 85, 112007 (2012), 1203.6810.
  * [125] CMS, S. Chatrchyan et al.,  JHEP 11, 067 (2012), 1208.2671.
  * [126] CMS, S. Chatrchyan et al.,  Phys. Lett. B 720, 83 (2013), 1212.6682.
  * [127] CMS, S. Chatrchyan et al.,  Eur. Phys. J. C 73, 2386 (2013), 1301.5755.
  * [128] CMS, S. Chatrchyan et al.,  JHEP 05, 065 (2013), 1302.0508.
  * [129] CMS, S. Chatrchyan et al.,  JHEP 02, 024 (2014), 1312.7582,  [Erratum: JHEP 02, 102 (2014)].
  * [130] ATLAS, G. Aad et al.,  Eur. Phys. J. C 76, 538 (2016), 1511.04716.
  * [131] ATLAS, M. Aaboud et al.,  Phys. Rev. D 94, 092003 (2016), 1607.07281.
  * [132] CMS, A. M. Sirunyan et al.,  Eur. Phys. J. C 77, 459 (2017), 1703.01630.
  * [133] NuTeV, M. Goncharov et al.,  Phys. Rev. D64, 112006 (2001), hep-ex/0102049.
  * [134] E. L. Berger, J. Gao, C. S. Li, Z. L. Liu, and H. X. Zhu,  Phys. Rev. Lett. 116, 212002 (2016), 1601.05430.
  * [135] BCDMS, A. C. Benvenuti et al.,  Phys. Lett. B223, 485 (1989).
  * [136] NMC, M. Arneodo et al.,  Nucl. Phys. B483, 3 (1997), hep-ph/9610231.
  * [137] L. W. Whitlow, E. M. Riordan, S. Dasu, S. Rock, and A. Bodek,  Phys. Lett. B282, 475 (1992).
  * [138] L. W. Whitlow, S. Rock, A. Bodek, E. M. Riordan, and S. Dasu,  Phys. Lett. B250, 193 (1990).
  * [139] E665, M. R. Adams et al.,  Phys. Rev. D54, 3006 (1996).
  * [140] NuTeV, M. Tzanov et al.,  Phys. Rev. D74, 012008 (2006), hep-ex/0509010.
  * [141] NMC, M. Arneodo et al.,  Nucl. Phys. B487, 3 (1997), hep-ex/9611022.
  * [142] H1, ZEUS, H. Abramowicz et al.,  Eur. Phys. J. C 73, 2311 (2013), 1211.1182.
  * [143] H1, F. D. Aaron et al.,  Phys. Lett. B665, 139 (2008), 0805.2809.
  * [144] H1, F. D. Aaron et al.,  Eur. Phys. J. C71, 1579 (2011), 1012.4355.
  * [145] ZEUS, S. Chekanov et al.,  Phys. Lett. B682, 8 (2009), 0904.1092.
  * [146] CHORUS, G. Onengut et al.,  Phys. Lett. B632, 65 (2006).
  * [147] H1, ZEUS, F. Aaron et al.,  JHEP 01, 109 (2010), 0911.0884.
  * [148] A. Candido, S. Forte, and F. Hekhorn,  JHEP 11, 129 (2020), 2006.07377.
  * [149] J. Collins, T. C. Rogers, and N. Sato,  Phys. Rev. D 105, 076010 (2022), 2111.01170.
  * [150] SeaQuest, J. Dove et al.,  Nature 590, 561 (2021), 2103.04024,  [Erratum: Nature 604, E26 (2022)].
  * [151] T. Cridge, L. A. Harland-Lang, A. D. Martin, and R. S. Thorne,  Eur. Phys. J. C 81, 744 (2021), 2106.10289.
  * [152] Particle Data Group, P. A. Zyla et al.,  PTEP 2020, 083C01 (2020).
  * [153] M. Cacciari, F. A. Dreyer, A. Karlberg, G. P. Salam, and G. Zanderighi,  Phys. Rev. Lett. 115, 082002 (2015), 1506.02660,  [Erratum: Phys.Rev.Lett. 120, 139901 (2018)].
  * [154] F. A. Dreyer and A. Karlberg,  Phys. Rev. Lett. 117, 072001 (2016), 1606.00840.
  * [155] M. Bonvini,  https://www.ge.infn.it/~bonvini/higgs/ , version 4.1 (2020).
  * [156] M. Cacciari, F. Dreyer, A. Karlberg, G. Salam, and G. Zanderighi,  https://provbfh.hepforge.org/ , version 2.0.2 (2018).
  * [157] A. Buckley et al.,  Eur. Phys. J. C 75, 132 (2015), 1412.7420.