Ad Lucem: QED Parton Distribution
Functions in the MMHT Framework

L. A. Harland–Lang1, A. D. Martin2, R. Nathvani3, R. S. Thorne3.

1Rudolf Peierls Centre, Beecroft Building, Parks Road, Oxford, OX1 3PU

2Institute for Particle Physics Phenomenology, University of Durham, Durham,
DH1 3LE

3Department of Physics and Astronomy, University College London, WC1E 6BT, UK

We present the MMHT2015qed PDF set, resulting from the inclusion of QED
corrections to the existing set of MMHT Parton Distribution Functions (PDFs),
and which contain the photon PDF of the proton. Adopting an input distribution
from the LUXqed formulation, we discuss our methods of including QED effects
for the full, coupled DGLAP evolution of all partons with QED at
$\mathcal{O}(\alpha)$, $\mathcal{O}(\alpha\alpha_{S})$,
$\mathcal{O}(\alpha^{2})$. While we find consistency for the photon PDF of the
proton with other recent sets, building on this we also present a set of QED
corrected neutron PDFs and provide the photon PDF separated into its elastic
and inelastic contributions. The effect of QED corrections on the other
partons and the fit quality is investigated, and the sources of uncertainty
for the photon are outlined. Finally we explore the phenomenological
implications of this set, giving the partonic luminosities for both the
elastic and inelastic contributions to the photon and the effect of our photon
PDF on fits to high mass Drell-Yan production, including the photon–initiated
channel.

###### Contents

  1. 1 Introduction
  2. 2 Including QED Effects in the MMHT Framework
    1. 2.1 Baseline QCD Fit
    2. 2.2 Input Photon Distribution
    3. 2.3 Modifications to DGLAP evolution
      1. 2.3.1 PDF Basis
      2. 2.3.2 Target Mass and Higher Twist Corrections
    4. 2.4 Separation of Elastic and Inelastic Components
    5. 2.5 Momentum Conservation
  3. 3 QED Neutron PDFs
    1. 3.1 Modified DGLAP Evolution
    2. 3.2 The Photon PDF of the Neutron
  4. 4 Results
    1. 4.1 Changes to PDFs due to QED corrections
    2. 4.2 Results of Global fits with QED corrections
      1. 4.2.1 The quality of the global fits
      2. 4.2.2 The photon PDF of the proton
      3. 4.2.3 QED Corrected Structure Functions
      4. 4.2.4 Effects of QED on $\alpha_{S}$ determination in the global PDF fit
      5. 4.2.5 Photon-photon luminosity
    3. 4.3 Uncertainties on the photon PDF
  5. 5 High Mass Drell-Yan
    1. 5.1 QED and Photon PDF sensitivity in High Mass Drell-Yan
    2. 5.2 Comparison with ATLAS Drell-Yan data
    3. 5.3 Including the ATLAS data in the global fit
  6. 6 Conclusions

##  1 Introduction

The precision physics program at the Large Hadron Collider (LHC) aims to
observe processes at an unprecedented level of accuracy and experimental
sensitivity. As part of these efforts, the analyses conducted by the LHC
experimental collaborations are increasingly undertaken with theoretical cross
section predictions at next–to–next–to–leading order (NNLO) in QCD, which
includes $\mathcal{O}(\alpha_{S}^{2})$ corrections. At this level of
precision, it is expected that electroweak (EW) corrections, including those
with photon–initiated (PI) processes, will begin to have observable effects as
$\alpha_{QED}\sim\alpha_{S}^{2}$ at the typical scales being probed at the
LHC. These should therefore be incorporated in theoretical predictions. In
particular, electroweak corrected partonic cross sections should be calculated
with corresponding Parton Distribution Functions (PDFs) produced at NLO and
NNLO in QCD and the appropriate order in QED. This is achieved primarily by
modifying the DGLAP [1, 2, 3] factorisation scale evolution of the PDFs to
include QED parton splittings. The most significant effect of this change is
the necessary inclusion of the photon as a constituent parton of the proton.
Subsequently one can also begin to calculate the effect of PI sub-processes as
corrections to the leading QCD cross section for processes such as Drell–Yan
[4], EW boson–boson scattering [5] and Higgs production with an associated EW
boson [6], which are expected to be sensitive to these effects. In a different
context, semi-exclusive [7] and exclusive production of states with EW
couplings are also related to the photon content of the proton, albeit not
directly to the inclusive photon PDF. Here, PI processes play an important
role, see e.g. [8, 9] for recent studies in the context of compressed SUSY
scenarios.

MRST provided the first such publicly available QED set[10], modelling the
photon at the input scale as arising radiatively from the quarks (and their
respective charges) below input, with DGLAP splitting kernels at
$\mathcal{O}(\alpha)$ in QED. Other such sets were subsequently developed that
either adopted similar phenomenological models[11], or sought to constrain the
photon in an analogous way to other partons by fits to Drell-Yan data [12,
13], first developed by the NNPDF Collaboration. These early sets saw
relatively large discrepancies between photon PDFs. Large modelling
uncertainties persisted due to the freedom in the choice of scale above which
photons are produced radiatively, modelled in the MRST set as the difference
between the current and constituent quark masses, while the approach taken by
the CTEQ14QED set [11] was to attempt to fit a parameterisation based on the
total momentum carried by the photon from $ep\rightarrow e\gamma+X$ data. In
the case of NNPDF2.3QED [12], the constraints available directly from data
were rather weak, due to the small size of the PI contributions. This lead to
large photon PDF errors, with a $\mathcal{O}(100\%)$ uncertainty at high $x$.
In all cases the available data was unable to constrain the photon to a high
degree of accuracy.

A final significant drawback of these early sets was that the majority did not
account for the contribution to the photon PDF from elastic scattering, in
which the proton coherently emits electromagnetic radiation without
disintegration, in contrast to photon contributions previously accounted for
from inelastic scattering processes, assumed to arise from quark splittings.
This distinction between the elastic and inelastic photon emission was one
that was seldom systematically treated, if considered at all.

Significant strides have been made in recent years to overcome these
deficiencies. First, more accurate determinations of the photon distribution
at input have been developed by making use of the experimentally well
determined elastic form factors of the proton, as in [14] and further
developed in [7, 15]. More precisely, the photon PDF corresponds to the flux
of emitted photons within the context of the equivalent photon approximation,
and as discussed in some of the early work on this [16], the contributions
from elastic and inelastic emission to the photon PDF are directly related to
the corresponding structure functions ($F_{1,2}^{el}$, $F_{1,2}^{inel}$)
probed in lepton–proton scattering. This idea has been revived in various
works over the previous decades [17, 18, 19, 20], and has most recently been
demonstrated within a rigorous and precise theoretical framework by the LUXqed
group [21, 22], where the first publicly available photon PDF applying this
approach was also provided. As the elastic and inelastic proton structure
functions have been determined experimentally to high precision, this has in
turn allowed for the determination of the elastic and inelastic contributions
to the photon to the level of a few percent. In addition to these
developments, QED DGLAP splitting kernels have now been calculated to
$\mathcal{O}(\alpha\alpha_{S})$ [23] and $\mathcal{O}(\alpha^{2})$ [24], whose
effects, as shown in Section 2.3, are not insignificant to the evolution of
the photon and other partons. In light of this, a greater confidence may be
had regarding the effects of QED modified partons and their impact on cross
section calculations.

In this paper, we outline the efforts undertaken by the MMHT group to develop
a fully consistent set of QED partons, adopting the LUXqed formulation at
input scale $Q_{0}$ for the photon. QED splitting kernels to
$\mathcal{O}(\alpha)$, $\mathcal{O}(\alpha_{S})$ and $\mathcal{O}(\alpha^{2})$
are incorporated into the DGLAP evolution and the effect of this is explored.
Furthermore, we also adopt a model for higher–twist (HT) effects in the quarks
at low $Q^{2}$, as the evolution of the photon PDF is sensitive to these
corrections, due to a lower input scale used in comparison to that of other
PDF sets.

As well as the conventional set of QED altered PDFs, we provide grids for the
photon PDF separated into its elastic and inelastic components, as well as a
consistent set of QED corrected neutron PDFs. Although the phenomenological
implications of a neutron set are limited, their production is necessary for a
consistent fit to deuteron and nuclear fixed target data from neutrino
($\nu$N) DIS scattering experiments used to constrain the PDFs. The QED
corrected neutron PDFs of MRST [10] provided isospin violating partons, with
$u_{(p)}\neq d_{(n)}$, and these were seen to reduce the NuTeV
$\sin^{2}\theta_{W}$ anomaly [25]. The breaking of isospin symmetry may also
have implications for the development of nuclear PDFs, and our current
treatment develops this earlier approach, providing new predictions for the
magnitude of isospin violation.

Finally, we will explore the phenomenological consequences of this set,
demonstrating the effects of QED incorporation on $F_{2}(x,Q^{2})$ as
calculated from PDFs, the partonic luminosities as a function of centre-of-
mass (CoM) energy and the change in fit quality after refitting the partons
with QED. We also explore the consequences of fitting to ATLAS high-mass
Drell-Yan data [26], with both QED effects and PI corrections to the cross
section produced by our set. We find that the effect of a fully coupled QED
DGLAP evolution is non-negligible on the gluon and quark PDFs.

##  2 Including QED Effects in the MMHT Framework

In this section we describe how the MMHT framework has been modified to
incorporate the QED splitting kernels in DGLAP evolution and the form we take
for the input distribution of the photon, and discuss their effect on the
final set of partons and the corresponding PDF uncertainties.

###  2.1 Baseline QCD Fit

Throughout this paper, in order to meaningfully interpret the effects of
including QED effects, we will compare the new partons to a baseline set of
PDFs evolved and fit solely with QCD kernels (at, unless explicitly stated,
NNLO). However, this set differs from the most recent public release of
partons, MMHT2014 [27]. In particular, this more closely corresponds to the
set described in [28], where the HERA Run I + II combined cross section data
[29] have been included in the fit. Furthermore, we now include some
additional data on $t\bar{t}$ production ($\sigma(t\bar{t})$) from the ATLAS
and CMS collaborations. In addition, further small amendments have been made
to the NLO and NNLO QCD kernels in the evolution, as detailed in Section 2.3.
Hence, we refer to this as the MMHT2015 PDF set and the PDFs with the QED
effects included as MMHT2015qed.

###  2.2 Input Photon Distribution

To generate PDFs from QED corrected DGLAP evolution requires an input
distribution for the photon at some starting scale, $Q_{0}$, from which the
PDFs may be evolved to higher scales. In principle, the photon input may be
parameterised in a form similar to other partons, which in MMHT primarily uses
an expansion in a basis of Chebyshev polynomials (as discussed in Section 2.1
of [27] and initially investigated in [30]). Photon input distributions based
on such an approach have been disfavoured by most groups due to the
insufficient constraints provided directly from data when simultaneously
fitting all of the partons. In particular, freely parameterising the photon
(in a suitable expansion basis, analogous to the other partons) in a global
fit is seen to lead to large uncertainties [12].

As discussed above, a significantly more precise approach is to formulate the
photon PDF in terms proton structure functions. This allow a precisely
constrained input PDF to be directly obtained from data for lepton-proton
scattering; i.e. from the experimentally determined values of $F_{2}$ and
$F_{L}$. We are always considering photon exchange in what follows, we will
implicitly be referring to the Neutral Current (NC) structure functions
wherever mentioned (i.e. $F_{2}\equiv F_{2}^{NC}$, $F_{L}\equiv F_{L}^{NC}$).
Moreover, as $Q^{2}_{0}=1\,\textrm{GeV}^{2}\ll M_{Z}^{2}$ we can safely
neglect any contributions from the weak neutral current and related
interference terms. The input expression for the photon PDF used in
MMHT2015qed is derived from that of LUXqed [21] with some modification. At a
given input scale, $\mu^{2}=Q_{0}^{2}$, we take the photon PDF to be:

$\begin{split}x\gamma(x,Q_{0}^{2})=\frac{1}{2\pi\alpha(Q_{0}^{2})}\int_{x}^{1}\frac{dz}{z}\Big{\\{}\int_{\frac{x^{2}m_{p}^{2}}{1-z}}^{\frac{Q_{0}^{2}}{1-z}}\frac{dQ^{2}}{Q^{2}}\alpha^{2}(Q^{2})\bigg{[}\bigg{(}zP_{\gamma,q}(z)+\frac{2x^{2}m_{p}^{2}}{Q^{2}}\bigg{)}F_{2}(x/z,Q^{2})\\\
-z^{2}F_{L}(x/z,Q^{2})\bigg{]}-\alpha^{2}(Q_{0}^{2})z^{2}F_{2}(x/z,Q_{0}^{2})\Big{\\}},\end{split}$
(1)

where $\alpha=\alpha_{\rm QED}$ and $P_{\gamma,q}(z)$ corresponds to the
$\mathcal{O}(\alpha)$ DGLAP splitting kernel given by:

$P_{\gamma,q}(z)=\frac{1+(1-z)^{2}}{z}.$ (2)

Note that the upper limit of the $Q^{2}$ integral introduces a dependency on
terms at scales higher than the input scale. It is more convenient to recast
Eq. (1) such that the photon at input is purely dependent on contributions
from $Q^{2}<Q_{0}^{2}$, with all $Q^{2}>Q_{0}^{2}$ dependence driven by DGLAP
evolution. To achieve this, we separate the $Q^{2}$ range of the integral into
two, with

$\begin{split}x\gamma(x,Q_{0}^{2})=\frac{1}{2\pi\alpha(Q_{0}^{2})}\int_{x}^{1}\frac{dz}{z}\Big{\\{}\int_{\frac{x^{2}m_{p}^{2}}{1-z}}^{Q_{0}^{2}}\frac{dQ^{2}}{Q^{2}}\alpha^{2}(Q^{2})\bigg{[}\bigg{(}zP_{\gamma,q}(z)+\frac{2x^{2}m_{p}^{2}}{Q^{2}}\bigg{)}F_{2}(x/z,Q^{2})\\\
-z^{2}F_{L}(x/z,Q^{2})\bigg{]}+\int_{Q_{0}^{2}}^{\frac{Q_{0}^{2}}{1-z}}\frac{dQ^{2}}{Q^{2}}\alpha^{2}(Q^{2})\bigg{[}\bigg{(}zP_{\gamma,q}(z)+\frac{2x^{2}m_{p}^{2}}{Q^{2}}\bigg{)}F_{2}(x/z,Q^{2})\bigg{]}\\\
-\alpha^{2}(Q_{0}^{2})z^{2}F_{2}(x/z,Q_{0}^{2})\Big{\\}},\end{split}$ (3)

where we have dropped the $F_{L}$ term in the second $Q^{2}$ integrand for
simplicity. This can be justified on the grounds that $F_{L}\ll F_{2}$ and
also by consideration of the fact that $F_{L}\sim\mathcal{O}(\alpha_{S})$ in
the parton model, while the expression given in Eq. (1) is formally only
accurate to $\mathcal{O}(\alpha\alpha_{S},\alpha^{2})$. A more thorough
discussion of this is given in Section 3 of [22].

By taking note of the fact that the scale variation of $F_{2}(Q^{2})$ and
$\alpha(Q^{2})$ may be treated as stationary at the order we are calculating
at ($\partial F_{2}/\partial Q^{2}$, $\partial\alpha/\partial Q^{2}\sim 0$),
we get

$\begin{split}x\gamma(x,Q_{0}^{2})=\frac{1}{2\pi\alpha(Q_{0}^{2})}\int_{x}^{1}\frac{dz}{z}\Big{\\{}\int_{\frac{x^{2}m_{p}^{2}}{1-z}}^{Q_{0}^{2}}\frac{dQ^{2}}{Q^{2}}\alpha^{2}(Q^{2})\bigg{[}\bigg{(}zP_{\gamma,q}(z)+\frac{2x^{2}m_{p}^{2}}{Q^{2}}\bigg{)}F_{2}(x/z,Q^{2})\\\
-z^{2}F_{L}(x/z,Q^{2})\bigg{]}-\alpha^{2}(Q_{0}^{2})\bigg{(}z^{2}+\ln(1-z)zP_{\gamma,q}(z)-\frac{2x^{2}m_{p}^{2}z}{Q_{0}^{2}}\bigg{)}F_{2}(x/z,Q_{0}^{2})\Big{\\}}.\end{split}$
(4)

This is the final expression for the input photon PDF that we will use
throughout this paper, taking $Q_{0}^{2}=1\,{\rm GeV}^{2}$ as the input scale.
We note that this closely resembles Eq.(4.10) of [22], however in our case we
retain the term of order $\mathcal{O}(m_{p}^{2}/Q_{0}^{2})$ as this is more
significant for the lower input scale we consider in comparison to LUXqed,
which uses $Q_{0}^{2}=10\,{\rm GeV}^{2}$. We now elaborate on the composition
of $F_{2,L}$ and how each source contributes to our expression for
$x\gamma(x,Q_{0}^{2})$. As discussed in the previous section, $F_{2,L}$
receive contributions from both elastic and inelastic scattering processes, as
shown in Fig. 1. In other words:

$F_{2,L}=F_{2,L}^{(el)}+F_{2,L}^{(inel)}.$ (5)

As we will discuss below, the elastic and inelastic components of $F_{2,L}$
are obtained from fits to data, largely in the same way as in LUXqed [21, 22].

For $F_{2,L}^{(el)}$ we use the A1 collaboration fit [31] to elastic
scattering data, which is provided in terms of the Sachs electric and magnetic
form factors for the proton:

$\begin{split}F_{2}^{(el)}(x,Q^{2})=\frac{[G_{E}(Q^{2})]^{2}+\tau[G_{M}(Q^{2})]^{2}}{1+\tau}\delta\Big{(}1-x\Big{)},\\\
F_{L}^{(el)}(x,Q^{2})=\frac{[G_{E}(Q^{2})]^{2}}{\tau}\delta\Big{(}1-x\Big{)},\end{split}$
(6)

where $\tau=Q^{2}/(4m_{p}^{2})$. We note that the fits from the A1
collaboration differ from the widely used dipole approximation by about $10\%$
at $x\sim 0.5$; above this the difference increases further but this has
little impact due to the effective kinematic cut at high $x$, discussed below.
However, as discussed in [21], the dipole model’s reasonably good
($\mathcal{O}(5\%)$) correspondence to the data at low $x$ makes it useful in
interpreting the scaling behaviour in this region
($\gamma^{(el)}(x)\sim\alpha\ln(1/x)$).

By substituting Eq. (6) into Eq. (4) we obtain an explicit formula for the
elastic contribution to the photon PDF at a scale $\mu$,

$\begin{split}x\gamma^{(el)}(x,\mu^{2})=\frac{1}{2\pi\alpha(\mu^{2})x}\int_{\frac{x^{2}m_{p}^{2}}{1-x}}^{\mu^{2}}\frac{dQ^{2}}{Q^{2}}\alpha^{2}(Q^{2})\bigg{[}\bigg{(}xP_{\gamma,q}(x)+\frac{2x^{2}m_{p}^{2}}{Q^{2}}\bigg{)}\times\\\
\frac{[G_{E}(Q^{2})]^{2}+\tau[G_{M}(Q^{2})]^{2}}{1+\tau}-x^{2}\frac{[G_{E}(Q^{2})]^{2}}{\tau}\bigg{]},\end{split}$
(7)

which, noting the presence of the $1/\alpha(\mu^{2})$ factor outside the
integral, is equivalent to the order to which we calculate to solving the
coupled DGLAP evolution for $\gamma^{el}$.

Turning to $F_{2,L}^{(inel)}$, this displays two distinct modes of behaviour.
For the continuum $W^{2}\gtrsim 4$ GeV2 region, the $x,Q^{2}$ dependence of
$F_{2,L}$ is seen to be relatively smooth, while in the resonance
$W^{2}\lesssim 3$ GeV2 region, various Breit-Wigner type resonances
contribute, due to the presence of hadronic excited states such as the
$\Delta$ and associated modes. To describe both of these regions, two
different fits are used above and below a threshold of $W_{\rm cut}^{2}=3.5$
GeV2. For the continuum ($W^{2}\geq W_{\rm cut}^{2}$) region, we use the
HERMES GD11-P [32] fit, while for the resonance ($W^{2}<W_{\rm cut}^{2}$)
region we take a fit to data from the CLAS collaboration [33].

The HERMES collaboration [32] provides data for $F_{L}$ by relating it to the
available data for $F_{2}$. In particular, by considering the parameter
$R=\sigma_{L}/\sigma_{T}$, the ratio of the longitudinal and transverse
polarisation cross sections, the two structure functions are related in the
following manner:

$F_{L}(x,Q^{2})=\Big{(}1+\frac{4m_{p}^{2}x^{2}}{Q^{2}}\Big{)}\frac{R(x,Q^{2})}{1+R(x,Q^{2})}F_{2}(x,Q^{2})\;,$
(8)

where the function $R(x,Q^{2})$, following the approach taken by HERMES, is
adapted from the E143 collaboration fit, $R_{1999}$ [34]. Although only
$F_{2}$ data is provided by the CLAS fit, $F_{L}$ is estimated in the
resonance region by using Eq. (8), with the same form of $R(x,Q^{2})$ provided
by HERMES.

Figure 1: Leading order representations of elastic (left) and inelastic
(right) NC lepton–proton scattering processes.

The structure functions themselves exhibit enhanced sensitivity to particular
effects at lower starting scales (1 GeV in the MMHT framework, in comparison
to 10 GeV adopted by LUXqed) such as proton mass corrections
$\mathcal{O}(m_{p}^{2}/Q^{2})$ and higher twist terms. Hence, modifications
are made to account for these during the evolution, as discussed in subsection
2.3.

Finally, we note that the lower bound of the $Q^{2}$ integral in Eq. (4)
introduces a cut on all photon contributions above a certain point in $x$. In
particular, by noting that the integral in $z$ is bounded by $x$, at the
limits of the integral the following inequality is imposed:

$Q^{2}\geq\frac{x^{2}m_{p}^{2}}{1-x},$ (9)

which may be rearranged to express an upper limit on $x$ for
$x\gamma(x,Q^{2})$:

$x\leq\frac{-Q^{2}+Q\sqrt{Q^{2}+4m_{p}^{2}}}{2m_{p}^{2}}\equiv x_{\rm cut},$
(10)

such that all contributions for $x>x_{\rm cut}$ vanish. As the expression at
input, Eq. (4) is valid at all scales, we include this cut at all stages of
the evolution, not just for the input photon. As discussed in Section 2.3,
this leads to a dampening. As $Q\rightarrow\infty$, $x_{cut}\rightarrow 1$ and
this constraint disappears rapidly, e.g. at $Q^{2}=10$ GeV2 we have $x_{\rm
cut}=0.918$. On the other hand, at the starting scale, we have $x_{cut}\simeq
0.62$. This cut has the effect of dampening the effects of other terms
relevant to the photon evolution at high $x$, such as higher twist and target
mass corrections, which are most prevalent at low $Q^{2}$, as well introducing
a source of momentum sum rule violation, as discussed in Section 2.5.

###  2.3 Modifications to DGLAP evolution

####  2.3.1 PDF Basis

In this section we outline the changes made to our evolution procedure to
accommodate the effects of QED. First, we distinguish between the basis (the
linearly independent combinations of partons) as parameterised in the fit and
in the evolution. The partons are parameterised solely at the input scale, the
majority of which, as discussed in subsection 2.1 of [27] and studied in [30],
are based on an expansion in terms of Chebyshev polynomials ($T_{i}^{Ch}(y)$):

$xf(x,Q_{0}^{2})=A(1-x)^{\eta}x^{\delta}\Bigg{(}1+\sum\limits_{i=1}^{n}a_{i}T_{i}^{Ch}(y(x))\Bigg{)},$
(11)

with $y=1-2\sqrt{x}$ and $n=4$. Hence, the free parameters are $A,\eta,\delta$
and $a_{i}$, where some of the $A$ are fixed by sum rules. Distributions of
this form are used for $f=u_{V},d_{V},S,(s+\bar{s)}$, where $S$ denotes the
light-quark sea distribution:

$S=2(\bar{u}+\bar{d})+s+\bar{s},$ (12)

For the differences $\bar{d}-\bar{u}$ and $s-\bar{s}$ a reduced
parameterisation is taken, reflecting the inability of current data to
constrain these distributions to a high degree of precision. The gluon is
provided in a form similar to Eq. (11), but with an additional term:

$xg(x,Q_{0}^{2})=A_{g}(1-x)^{\eta_{g}}x^{\delta_{g}}\Bigg{(}\sum_{n=1}^{2}a_{g,i}T_{i}^{Ch}(y(x))\Bigg{)}+A_{g^{\prime}}(1-x)^{\eta_{g^{\prime}}}x^{\delta_{g^{\prime}}}\;,$
(13)

which is found to significantly improve the quality of the global fit [35],
and essentially provides more freedom for the gluon at low $x$. Here,
$\eta_{g}$, $\eta_{g}^{\prime}$, $A_{g}$ and $A_{g}^{\prime}$ are all
correlated in the fit, and their dependency can be artificially disrupted with
the introduction of QED effects, leading to significant changes in the gluon
if the partons are not refit, see Section 4.1.

In the MMHT framework the QED supplemented evolution of the partons is
unidirectional in $Q^{2}$ from a starting scale of $Q_{0}=1$ GeV, with the
convolution for the partons at each step performed in $x$ space. This is in
contrast to that of NNPDF3.1luxQED case [36], which adopts an iterative
process in its fit that reverses the DGLAP evolution of the partons from
$Q=100$ GeV for a photon produced from the NNPDF partons at high scales (with
an elastic and low $Q^{2}$ contribution whose expression is the same as for
LUXqed) then aims to find a consistent starting scale photon, subject to the
momentum sum rule constraint for the partons, modified to include the photon
($\gamma$):

$\int_{0}^{1}x(\Sigma(x,Q_{0}^{2})+g(x,Q_{0}^{2})+\gamma(x,Q_{0}^{2}))=1,$
(14)

where $\Sigma$ is the total singlet for the quarks. In practice, as we will
see the resulting photon distributions from either approach are in agreement,
differing only on the order of the uncertainties. However, due to certain
higher twist effects and the procedure adopted for the treatment of our
elastic photon distribution, $\gamma^{(el)}$, Eq. (14) is not strictly obeyed
during the evolution, see Section 2.5 for further discussion.

While eqs. 11 to 13 reflect the input distribution parameterisations, a
different and distinct linear combination of the partons are involved in the
evolution procedure itself. Previously in the MMHT framework, the pure QCD
DGLAP evolution of the partons, at all orders, was performed in a basis that
was chosen for computational efficiency. This involved a decoupling of the
partons into a singlet (consisting of the gluon and flavour combinations of
quark and antiquark distributions) and non-singlet distributions which are
evolved separately. Explicitly, the linearly independent combinations of
partons that were evolved consisted of the following singlet (in the space of
quark flavours) combinations:

$\Sigma_{L}=u+\bar{u}+d+\bar{d}+s+\bar{s},$ (15) $c+\bar{c},\hskip
56.9055ptb+\bar{b},$ (16) $g$ (17)

and the following non-singlet combinations:

$u_{V}+d_{V}=u-\bar{u}+d-\bar{d},$ (18) $u+\bar{u}-\frac{\Sigma_{L}}{3},\hskip
56.9055pt-(s+\bar{s})+\frac{\Sigma_{L}}{3},$ (19)
$\frac{d_{V}-u_{V}}{2},\hskip
28.45274pt(s-\bar{s})-\frac{u_{V}+d_{V}}{2},\hskip
28.45274pt(c-\bar{c})-\frac{u_{V}+d_{V}}{2},\hskip
28.45274pt(b-\bar{b})-\frac{u_{V}+d_{V}}{2},$ (20)

where the subscript $L$ in $\Sigma_{L}$ denotes the fact that the singlet
consists only of the light quarks. The charm and bottom singlet distributions
in Eq. (16) are evolved separately since they only become non-zero near the
relevant mass thresholds for production. When considering QCD in isolation,
the $SU(n_{f})$ flavour invariance of the splitting kernels allows such
distributions to be evolved consistently.

Now, the introduction of QED splitting kernels, $P_{ij}^{(QED)}$, in DGLAP
evolution necessarily prohibits such combinations from being used. Writing the
pure QCD splitting kernels via the usual perturbative expansion

$P_{ij}^{(QCD)}=\frac{\alpha_{S}}{2\pi}P_{ij}^{(1)}+\left(\frac{\alpha_{S}}{2\pi}\right)^{2}P_{ij}^{(2)}+\left(\frac{\alpha_{S}}{2\pi}\right)^{3}P_{ij}^{(3)}+...$
(21)

for the QED and mixed QED–QCD case, recent theoretical work [23, 24] enable
the terms

$P_{ij}^{(QED)}=\frac{\alpha}{2\pi}P_{ij}^{(0,1)}+\frac{\alpha\alpha_{S}}{(2\pi)^{2}}P_{ij}^{(1,1)}+\left(\frac{\alpha}{2\pi}\right)^{2}P_{ij}^{(0,2)}+...$
(22)

to be used in the QED supplemented evolution. Here, the first and second
superscript indices denote the order in QCD and QED respectively, and the
second term in this expansion reflects mixed order splitting kernels.

Since the non-abelian nature of QCD does not manifest at leading order in
quark interactions, the majority of the splitting functions in QCD and QED are
simply related at this order:

$\displaystyle P_{qq}^{(0,1)}$
$\displaystyle=\frac{e_{q}^{2}}{C_{F}}P_{qq}^{(1,0)},$ $\displaystyle
P_{q\gamma}^{(0,1)}$ $\displaystyle=\frac{e_{q}^{2}}{T_{F}}P_{qg}^{(1,0)}\;,$
(23) $\displaystyle P_{\gamma q}^{(0,1)}$
$\displaystyle=\frac{e_{q}^{2}}{C_{F}}P_{gq}^{(1,0)},$ $\displaystyle
P_{\gamma\gamma}$
$\displaystyle=-\frac{2}{3}\sum_{i}^{n_{F}}e_{i}^{2}\delta(1-y),$ (24)

The exception is $P_{\gamma\gamma}$, which differs considerably from the
expression for $P_{gg}$, due to the purely gluonic contribution in the latter
case.

A further caveat regarding $P_{\gamma\gamma}$ is that we only include quark
loops, and not those due to leptons. In the latter case, consistency would
require the corresponding introduction of lepton PDFs, which in principle
enter amongst the partons discussed so far, due to splittings of the form
$\gamma\rightarrow l\bar{l}$. More precisely, for $Q^{2}>m_{l}^{2}$, lepton
splittings should also be incorporated into $P_{\gamma\gamma}$, such that the
sum over quarks is modified to include the leptons:

$\sum_{i}e_{i}^{2}=N_{C}\sum_{q}^{n_{F}}e_{q}^{2}+\sum_{l}^{n_{L}}e_{l}^{2}.$
(25)

In our framework, we neglect the latter term which accounts for leptonic
contributions to $P_{\gamma\gamma}$, since the contribution of the photon
itself enters as an $\mathcal{O}(\alpha)$ correction to the PDFs, with the
lepton contributions at $\mathcal{O}(\alpha^{2})$, implying they are extremely
suppressed. This was studied more extensively in [37] where it was found that
the magnitude of the lepton distributions were many orders of magnitude below
those of $x\gamma(x,Q^{2})$, with negligible effects on the PDFs at the scales
considered in this paper.

However, we note that the LUXqed PDF set [22] does include this contribution
in the DGLAP evolution used to develop their $x\gamma(x,Q^{2})$. Since the
right hand side of Eq. 24 is a $\delta(1-x)$ term multiplied by a negative
coefficient, the extra contributions from the lepton splitting terms in DGLAP
are anticipated to slightly reduce the magnitude of a photon whose evolution
accounts for them (as one anticipates from the process $\gamma\rightarrow
l\bar{l}$).

Upon inspection of eqs. 23 to 24, even at leading order it becomes apparent
that the distributions in eqs. 15 to 20 cannot be used since QED couplings no
longer support flavour symmetry, due to the charge separation of up and down
type quarks ($e_{u}\neq e_{d}$). Furthermore, one anticipates based on this
observation the breaking of isospin symmetry when comparing the valence
distributions of the proton and neutron, as discussed in Section 3.1. To
accommodate the requirement of charge sensitivity, the partons are now evolved
in the following basis, which are separable by charge:

$q^{\pm}_{i}=q_{i}\pm\bar{q_{i}},\hskip 14.22636ptg,\hskip
14.22636pt\gamma^{(el)},\hskip 14.22636pt\gamma^{(inel)}.$ (26)

In the following discussion the subscript $i$ denotes any active ($Q>2m_{q}$)
flavour: $i=u,d,s,c,b$ and the +/- superscript denotes the singlet and non-
singlet quark distributions respectively. The gluon and photon components,
$g$, $\gamma^{(el)}$ and $\gamma^{(inel)}$ are then evolved individually in
the flavour space of the partons.

Although the basis given in Eq. (26) is compatible with a joint evolution in
QCD and QED, they require some modification to the form of DGLAP splitting
kernels used. Writing $t=\ln(Q^{2})$ and $(f\otimes
g)(x)=\int_{x}^{1}\frac{dy}{y}f(x/y)g(y)$, the non-singlet distributions
described in Eqs. (18) - (20) may be evolved in the following way in pure QCD:

$\frac{\partial q_{i}^{NS}}{\partial t}=P_{q_{i}}^{-}\otimes q_{i}^{NS}\;,$
(27)

where the expression for $P_{q_{i}}^{-}$ may be found in Eqs. (4.94) to
(4.108) of [38]. The simplicity of this equation arises from the fact that
symmetry allows for evolution of the $p_{i}^{NS}$ distributions to be diagonal
in quark flavour space, such that only the term $P_{q_{i}}^{-}$, which
describes the diagonal elements of the quark-quark and quark-antiquark
splitting functions, is required.

The evolution for the $q_{i}^{-}$ requires an additional component since
although they are also non-singlet functions of the quarks, the non-diagonal
elements in flavour space become necessary to the evolution:

$\frac{\partial q_{i}^{-}}{\partial t}=P_{q_{i}}^{-}\otimes
q_{i}^{-}+\sum_{j=1}^{n_{F}}\Delta P^{S}\otimes q_{j}^{-},$ (28)

where $\Delta P^{S}$ becomes non-zero at NNLO ($\mathcal{O}(\alpha^{3})$) in
QCD and $n_{F}$ is the number of active quarks in the evolution.

Note that this sum over valence-like non-singlet distributions corresponded to
eq. (18) in the original MMHT framework, which neglected the strange, charm
and bottom distributions due to their small relative size. With the release of
the set described in this paper, the contribution from these off-diagonal
splittings for all flavours are now included, which represent minor changes,
$\mathcal{O}(10^{-5})$, in a like-for-like comparison with the original MMHT
partons purely in QCD.

####  2.3.2 Target Mass and Higher Twist Corrections

As previously noted, MMHT2015qed differs in its production of a photon PDF
from other contemporary sets in adopting a straightforward evolution in
$Q^{2}$ space, from a starting scale of $Q_{0}=1$ GeV. However, at low scales
such as these, target mass corrections, which account for the finite mass of
the proton, and higher twist terms have non-negligible contributions to
$F_{2,L}$. Above $Q_{0}$, the $F_{2}^{(inel)}$ contributions to
$\gamma^{(inel)}$, as in eq. (4), are modelled by the parton splittings in
DGLAP, which require some modification to capture the relevant behaviour at
high $x$.

The target mass corrections for the proton are well known, modifying the
$\mathcal{O}(\alpha)$ quark to photon splitting in an identical manner to the
first term in the integrand of Eq. (4):

$P_{\gamma,q}^{(0,1)}(z)\rightarrow
P_{\gamma,q}^{(0,1)}(z)+\frac{2x^{2}m_{p}^{2}}{zQ^{2}}.$ (29)

Further modifications are also required for higher twist terms which lead to
discrepancies between $F_{2}$ as calculated from the partons and experimental
measurements for $F_{2}^{(inel)}$, due to non-perturbative effects at high $x$
and low $Q^{2}$. In a global fit this effect is typically eliminated by
cutting on the low $W^{2}$ region where such corrections are relevant, however
for the determination of the photon PDF which is sensitive to
$F_{2}^{({inel})}$ in the region discussed, we must include this. Therefore a
phenomenological model must be adopted to account for such higher–twist
corrections. We follow the approach of [39], where non-perturbative $\sim
1/Q^{2}$ power corrections to the structure functions are provided, by
characterising the associated infrared divergences in field theory with the
so–called renormalon. In this paper we shall use the term renormalon
synonymously with higher twist corrections of this type. In [39], they provide
at $\mathcal{O}(1/Q^{2})$ a modification to $F_{2}$ that accounts for the
change due to renormalon calculations at high $x$, and this is found to give
an improved description of DIS data [40].

In lieu of $F_{2}^{({inel})}$, during the evolution the contributions to
$\gamma^{({inel})}$ are essentially generated by the quark splittings
($q\rightarrow q\gamma$), where the total quark singlet $\Sigma$ plays the
role of $F_{2}$ in eq. (4). Therefore, to approximate renormalon effects
during the evolution, these modification are instead made to the quarks via

$q(x,Q^{2})\rightarrow
q(x,Q^{2})\Big{(}1+\frac{A^{\prime}_{2}}{Q^{2}}\int^{1}_{x}\frac{dz}{z}C_{2}(z)q(\frac{x}{z},Q^{2})\Big{)},$
(30)

where $A^{\prime}_{2}$ is a parameter not given a priori by the theory and
$C_{2}(z)$ is defined in Eq. 4.1 of [39], and conserves the flavour number
properties of the various $q(x/z,Q^{2})$. As such, higher twist contributions
to $F_{2}$ do not contribute to the Adler sum rule,

$\int_{0}^{1}{d}x\,F^{HT}_{2}(x,Q^{2})=0\;,$ (31)

enforcing that these are well behaved as $x\rightarrow 0$. However no such
restriction applies to $F_{3}$, and renormalon calculations [41] imply that
they become large, necessitating the need for the more stringent cut on
$F_{3}$ data used in the fit (from the CHORUS collaboration [42]) that extend
into this region.

4120414041604180420042204240426000.10.20.30.40.50.6

$\chi^{2}$

$|A^{\prime}_{2}|$Total $\chi^{2}$ for different values of $A^{\prime}_{2}$

Figure 2: The $\chi^{2}$ values obtained in a global fit, with kinematic
constraints on DIS data lowered to $W^{2}=5$ GeV2, with different values of
$A^{\prime}_{2}$ in the renormalon calculations for $F_{2}$ and $F_{3}$ . The
dashed blue line represents a $\Delta\chi^{2}$ = 10 variation from the minimum
to establish an uncertainty band on $A^{\prime}_{2}$.

This is of interest because the parameter $A^{\prime}_{2}$ is not well
determined, and in [39], is fit loosely to structure function data to yield a
value of $A^{\prime}_{2}=-0.2$ GeV2. As discussed above, data sensitive to
renormalon contributions are typically excluded in global fits, to remove any
sensitivity to such non–perturbative effects. In particular, in MMHT kinematic
cuts of $W^{2}>15$ GeV2 (and $W^{2}>20$ GeV2 at LO) are taken, while for those
data sets relating to $\nu(\bar{\nu})N$ experiments to measure $xF_{3}$ a more
stringent cut of $W^{2}>25$ GeV is imposed [43]. However, with the aim of
determining a more precise value of $A^{\prime}_{2}$ we have relaxed these
constraints, lowering the threshold to $W^{2}>5$ GeV2 and modifying $F_{2}$
and $F_{3}$ to include the relevant renormalon contributions as in [39], i.e.
with modifications of the form shown in eq. (30).

In Fig. 2 the fit quality for different values of $A^{\prime}_{2}$ is shown.
We find that $A^{\prime}_{2}=-(0.3\pm 0.1)$ GeV2, with uncertainties
determined from a generous $\Delta\chi^{2}=\pm 10$ variation in the fit (to
one significant figure). This is motivated by the dynamical tolerance scheme
used in our framework, as outlined in Section 6 of [43], where it was found
that in order to provide reasonable uncertainties when fitting to many
disparate data sets in tension with one another, one typically requires
tolerances $T=\sqrt{\Delta\chi^{2}_{global}}\sim 3$ rather than the $T=1$ one
would obtain from a standard ‘parameter-fitting’ criterion. We note this
choice also corresponds to the fixed tolerance uncertainty schemes adopted by
early CTEQ sets [44]. The uncertainty on this is then propagated as an
independent source of uncertainty for the photon, as discussed in Section 4.3.
This represents a slightly larger renormalon contribution than predicted from
[39], though the data are unable to provide significant constraints in either
case.

As seen in Fig. 3, the target mass corrections lead to a $\sim 3\%$ increase
in the photon at high $x$, while the renormalon contributions, which provide
an increasingly positive contribution to $F_{2}$ at high $x$, correspondingly
enhance the photon at moderate to high $x$. Note that the turn around in both
figures at $x\simeq 0.5$ occurs due to the previously mentioned effective
kinematic cut on all photon contributions at high $x$ and low $Q^{2}$. This
cut itself is also a function of the proton mass $m_{p}$, though for our
purposes we consider the kinematic cut imposed due to the target mass (i.e.
the cut in $x$) as independent from the term introduced in the evolution and
it is seen that the two have opposite effects on the high $x$ photon, with the
kinematic cut ultimately dominating and the effect of the corrections to the
splitting function and the renormalon contribution being suppressed as $x\to
1$.

Since the target mass and renormalon contributions are both $\sim 1/Q^{2}$
corrections, their relative importance at higher scales is seen to decrease
slightly, as shown by a comparison of the red ($Q^{2}=100$ GeV2) and green
($Q^{2}=10^{4}$ GeV2) curves. We note that both the proton mass term and the
modification to the quarks in eq. (30) introduce small, independent sources of
momentum violation in the evolution, as discussed in Section 2.5.

11.0051.011.0151.021.0251.031.03500.10.20.30.40.50.60.70.8

$\frac{\gamma}{\gamma_{\frac{x^{2}}{Q^{2}}}}(x,Q^{2})$, $Q^{2}=100,10^{4}$
GeV2

$x$Ratio of the Photon PDF with and without $O(\frac{x^{2}}{Q^{2}})$
corrections$Q^{2}=100$GeV2$Q^{2}=10^{4}$GeV2

11.011.021.031.041.051.061.0700.10.20.30.40.50.60.70.8

$\frac{\gamma}{\gamma_{HT}}(x,Q^{2})$, $Q^{2}=100,10^{4}$ GeV2

$x$Ratio of the Photon PDF with and without HT
corrections$Q^{2}=100$GeV2$Q^{2}=10^{4}$GeV2

Figure 3: (Left) Ratio of the photon PDF with ($\gamma(x,Q^{2})$) and without
($\gamma_{\frac{x^{2}}{Q^{2}}}(x,Q^{2})$) target mass corrections and (right)
Higher Twist (renormalon) corrections.

###  2.4 Separation of Elastic and Inelastic Components

00.010.020.030.040.050.060.071e-050.00010.0010.010.11

$x\gamma(x,Q^{2})$, $Q=1$ GeV

$x$Comparisons of the Elastic and Inelastic Photon PDFElasticInelasticTotal
Photon

00.020.040.060.080.10.120.140.160.181e-050.00010.0010.010.11

$x\gamma(x,Q^{2})$, $Q=10$ GeV

$x$Comparisons of the Elastic and Inelastic Photon PDFElasticInelasticTotal
Photon

00.050.10.150.20.250.30.350.40.450.51e-050.00010.0010.010.11

$x\gamma(x,Q^{2})$, $Q=100$ GeV

$x$Comparisons of the Elastic and Inelastic Photon PDFElasticInelasticTotal
Photon

Figure 4: The Elastic and Inelastic Photon components at different values of
$Q$. Top left: $Q=1$ GeV, top right: $Q=10$ GeV, bottom: $Q=100$ GeV.

As noted in Section 2.2, the photon PDF actually comprises of two component
distributions,
$\gamma(x,Q^{2})=\gamma^{(el)}(x,Q^{2})+\gamma^{(inel)}(x,Q^{2})$, which
represent photon contributions from elastic and inelastic proton scattering
events, respectively. Separating $\gamma^{(el)}$ and $\gamma^{(inel)}$ from
one another while consistently performing the evolution for all the partons
required certain changes to be made from the standard procedure for performing
DGLAP, due to the fact that the generation of $\gamma^{(el)}$ in the evolution
is independent of parton splittings, as detailed below.

For $\gamma^{(inel)}$, the evolution is analogous to that of the other
partons. The contributions from the HERMES (continuum) and CLAS (resonance)
data for $F_{2}^{(inel)}$ are present only at input, above which DGLAP
evolution is performed. We emphasise that all photon contributions that arise
from the splitting of other partons (the quarks, antiquarks and both photon
components themselves, but also the gluon at $\mathcal{O}(\alpha\alpha_{S})$)
in DGLAP are absorbed into the definition of $\gamma^{(inel)}$ (using the
notation of the previous section):

$\frac{d\gamma^{(inel)}}{dt}=\sum_{j}^{n_{F}}P_{\gamma q_{j}}\otimes
q_{j}+\sum_{j}^{n_{F}}P_{\gamma\bar{q}_{j}}\otimes\bar{q}_{j}+P_{\gamma
g}\otimes g+P_{\gamma\gamma}\otimes\gamma^{(inel)}.$ (32)

This reflects the fact that scattering processes that are sensitive to the
partons are themselves inelastic and that therefore any photon contributions
that arise from their evolution in DGLAP are necessarily inelastic
contributions.

While $\gamma^{(el)}$ is included at input and passed to the other partons
during evolution, its own evolution requires consideration of the
contributions it receives above $Q_{0}$ from $F_{2}^{(el)}$, since our
expression for $\gamma^{(el)}$ given in Section 2.2, Eq. (7), holds generally
above the input scale. Incorporating this and splittings of the form
$\gamma\rightarrow q\bar{q}$ and $\gamma\rightarrow q\bar{q}g$ at
$\mathcal{O}(\alpha\alpha_{S})$, the evolution for $\gamma^{(el)}$ is given
as:

$\frac{d\gamma^{(el)}}{dt}=P_{\gamma\gamma}\otimes\gamma^{(el)}+\delta
x\gamma^{(el)}.$ (33)

The expression for $\delta x\gamma^{(el)}$ is given by taking the derivative
of the expression for the elastic photon, eq. (7), w.r.t $Q^{2}$:

$\begin{split}\delta
x\gamma^{(el)}(x,Q^{2})=\frac{\alpha(Q^{2})}{2\pi}\frac{1}{x}\bigg{[}\bigg{(}xP_{\gamma,q}(x)+\frac{2x^{2}m_{p}^{2}}{Q^{2}}\bigg{)}\frac{[G_{E}(Q^{2})]^{2}+\tau[G_{M}(Q^{2})]^{2}}{1+\tau}\\\
-x^{2}\frac{[G_{E}(Q^{2})]^{2}}{\tau}\bigg{]}.\end{split}$ (34)

As discussed in the next section, including the term introduced in Eq. (34) as
an external contribution (not generated from parton splittings but added into
the evolution from $F_{2}^{(el)}$ data) introduces a small amount of momentum
violation, as do subsequent splittings of the form $\gamma^{(el)}\rightarrow
q\bar{q}$.

Although the provisions outlined above are needed for the evolutions of
$\gamma^{(el)}$ and $\gamma^{(inel)}$, i.e. those contributions from splitting
functions of the form $P_{\gamma\\{q,\bar{q},g,\gamma\\}}$, the treatment for
the rest of the partons remains broadly unchanged. Since the quark, antiquark
and gluon contributions from $P_{\\{q,\bar{q},g,\gamma\\}\gamma}$ splittings
do not distinguish between $\gamma^{(el)}$ and $\gamma^{(inel)}$, the entire
photon contribution,
$\gamma(x,Q^{2})=\gamma^{(el)}(x,Q^{2})+\gamma^{(inel)}(x,Q^{2})$, is passed
to the relevant splitting kernels during evolution.

As $\gamma^{(el)}$ and $\gamma^{(inel)}$ distinguish between the photon in two
distinct categories of scattering processes, there is a phenomenological
interest in comparing the two. At input, the elastic contribution dominates
over that of the inelastic, as $F_{2}^{(el)}>F_{2}^{(inel)}$ in the region
$Q\lesssim 1$ GeV. However, evolution quickly enhances the contributions of
$\gamma^{(inel)}$, particularly at low $x$, predominantly due to quark
splittings, as shown in Figs. 4 and 5. As discussed above, the only
contributions $\gamma^{(el)}$ receives during the evolution are those from Eq.
(34). Since $G_{E,M}(Q^{2})$ are known to diminish with increasing $Q^{2}$ and
$1/\tau\sim 1/Q^{2}$, an inspection of the form of eq. (34) reveals that it
will be of diminishing importance in a significant range of $x$. In fact,
investigating the effects of leaving out this term in eq. (33) entirely
yielded a $\gamma^{(el)}$ with differences of just $\mathcal{O}(10^{-3})$ from
the form with the contributions included. However, the elastic distribution’s
contribution at input, and above, is proportionally large at high $x$, even at
high $Q^{2}$ (Fig. 5), and due to the kinematic cut all contributions to the
photon at the highest $x$ are from $Q^{2}>Q_{0}^{2}$. Indeed, in this region
the elastic contribution even above $Q_{0}^{2}$ dominates the photon, as can
be seen in Fig. 6, which shows the effect on the photon of turning this
contribution off. One slight caveat, however, is that as $\lim_{x\to
1}\gamma^{(el)},\gamma^{(inel)}\rightarrow 0$, and ultimately uncertainties
become large in this region (see Section 4.3), making it difficult to make
very strong predictive statements about either distribution in this region.

00.10.20.30.40.50.60.70.80.911e-050.00010.0010.010.1

$x\gamma(x,Q^{2})$, $Q=1$ GeV

$x$Relative contributions of the Elastic and Inelastic Photon
PDFElasticInelastic

0.10.20.30.40.50.60.70.80.91e-050.00010.0010.010.1

$x\gamma(x,Q^{2})$, $Q=10$ GeV

$x$Relative contributions of the Elastic and Inelastic Photon
PDFElasticInelastic

0.10.20.30.40.50.60.70.80.91e-050.00010.0010.010.1

$x\gamma(x,Q^{2})$, $Q=100$ GeV

$x$Relative contributions of the Elastic and Inelastic Photon
PDFElasticInelastic

Figure 5: The relative contributions of the Elastic and Inelastic Photon
components at different values of $Q$. Top left: $Q=1$ GeV, top right: $Q=10$
GeV, bottom: $Q=100$ GeV.

0.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.6

$\gamma_{Q^{2}<Q_{0}^{2}}/\gamma(x,Q^{2})$

$x$Ratio of the Photon PDF with and without $\delta\gamma^{(el)}$ for
$Q^{2}>Q_{0}^{2}$$Q^{2}=100$GeV2$Q^{2}=10^{4}$GeV2

00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.6

$\gamma^{(el)}_{Q^{2}<Q_{0}^{2}}/\gamma^{(el)}(x,Q^{2})$

$x$Ratio of the Elastic Photon PDF with and without $\delta\gamma^{(el)}$ for
$Q^{2}>Q_{0}^{2}$$Q^{2}=100$GeV2$Q^{2}=10^{4}$GeV2

Figure 6: A plot showing the ratio of the photon distribution with the elastic
contribution for $Q^{2}>Q_{0}^{2}$ removed and the total distribution. The
right plot shows the same ratio for only the elastic contribution to the
photon distribution.

###  2.5 Momentum Conservation

The inclusion of the photon PDF requires that the photon be included in the
momentum sum rule (14), naturally leading to a redistribution of momentum in
the other partons in order to obey eq. (14) at input. However, due to the
procedure adopted for the inclusion of $\gamma^{(el)}$, outlined in the
previous section, as well as higher twist terms, this equation is not strictly
obeyed during the evolution. This reflects the discrepancy between effects of
non-perturbative corrections, such as that of target masses, and the parton
model. In this section we outline the consequences of such changes.

First we discuss the effect of a kinematic cut on the photon, as introduced by
the lower limit of the integral in $Q^{2}$ in the expression for
$x\gamma(x,Q^{2})$, which as discussed in Section 2.2 has the effect of
introducing an effective cut on the photon PDF at high $x$ during the
evolution. In essence, this removal of photon contributions at high $x$ is a
target mass correction (since the cut has a dependence on $m_{p}^{2}$), which
is not required to obey the momentum conservation of the partons ordinarily
found in DGLAP evolution and therefore introduces small amounts of violation
(in the form of a reduction of total momentum carried by the partons) into eq.
(14) of $\mathcal{O}(10^{-3}\,\%)$. This is seen in Fig. 7 (left), where we
display the ratio of the total momentum of the partons with and without this
cut applied.

0.9999840.9999860.9999880.999990.9999920.9999940.9999960.9999981

1

10

100

$10^{3}$

$10^{4}$

$10^{5}$

$10^{6}$

$10^{7}$

$10^{8}$

$10^{9}$

Total momentum ratio ($Q^{2}$)

$Q^{2}$Ratio of total partonic momentum with and without cuts

0.99550.9960.99650.9970.99750.9980.99850.9990.99951

$1$

$10$

$100$

$10^{3}$

$10^{4}$

$10^{5}$

$10^{6}$

$10^{7}$

$10^{8}$

$10^{9}$

Photon momentum ratio ($Q^{2}$)

$Q^{2}$Ratio of $x\gamma$ momentum with and without cuts

Figure 7: (Left) A plot showing the ratio between the total momentum carried
by all the partons, at a given point in $Q^{2}$ of the evolution with the
kinematic cut $Q^{2}\geq x^{2}m_{p}^{2}/(1-x)$ applied to $x\gamma(x,Q^{2})$
and without. The right plot shows an identical plot but focusing solely on the
proportional difference in momenta caused by the photon, where the effects on
the evolution are seen to peak at $Q^{2}\sim 10$ GeV2.

In particular, Fig. 7 (right) indicates that the reduction to the total
momentum carried by the photon is, as anticipated, most strongly affected by
the kinematic cut at low scales until $Q^{2}\sim 10$ GeV2 (with total changes
of less than 1%). Since the overall momentum carried by the photon is small,
$\sim 2-3\times 10^{-3}$, at low scales where momentum violating effects are
most prevalent, this leads to the minuscule amount of change observed in the
total momentum of the partons.

We now discuss other effects during the evolution which contribute to
violation of the momentum sum rule. The momentum sum rule is constrained to be
obeyed by all the partons at the input scale, and both the inelastic and
elastic photons are considered when imposing the momentum sum rule for the
parameterisation of the quarks, as in eq. (14). However, above the input scale
the contribution to $\gamma^{(el)}$ that comes from the second term in (33),
that is due to elastic photon emission, will lead to some momentum sum rule
violation, as this contribution does not originate from standard DGLAP
evolution, and is not balanced by a corresponding loss of quark and antiquark
momentum, i.e., any $\gamma$ contribution from the quarks during evolution,
e.g. $q\rightarrow q+\gamma$ is absorbed into the definition of
$\gamma^{(inel)}$. ($\gamma^{(el)}$ is not entirely decoupled from the
evolution of the quarks, since $\gamma^{(el)}\rightarrow q\bar{q}$ splitting
are still permissible.) In practice, this effect is negligible, with momentum
violating effects of $O(10^{-4})$ observed in the sum rule during evolution,
and in fact stabilises at higher $Q^{2}$ where the elastic contribution is
less significant.

Similarly, the proton mass term given in eq. (29) naturally breaks the form of
momentum conservation usually obeyed between splitting functions of this type,
implied by the equation

$\int_{0}^{1}x\Big{[}P_{q,q}^{(0,1)}(x)+P_{\gamma,q}^{(0,1)}(x)\Big{]}=0.$
(35)

In essence, the proton mass term invalidates this relationship, though in
rapidly diminishing amounts as $1/Q^{2}\rightarrow 0$, leading to changes of
$\mathcal{O}(10^{-5})$ in the total momentum carried by the partons.

Likewise, other higher twist terms included in the evolution for the purposes
of QED lead to small amounts of momentum violation. Since the quark
distributions, $q_{i}(x,Q^{2})$, passed to both $P_{q,q}^{(0,1)}(x)$ and
$P_{\gamma,q}^{(0,1)}(x)$ differ due to the inclusion of renormalon
corrections for the latter but not the former, this aspect of the evolution
also invalidates momentum violation to a small degree, also shown in Fig. 8,
creating a small amount of violation of $\mathcal{O}(2\times 10^{-5})$.

Overall, even in conjunction, the combined magnitude of momentum sum rule
violation is less than $10^{-4}$. In practice, this total effect is less than
the momentum violation coming from the ‘leakage’ of the partons that occurs
due to the fact that the integration range during DGLAP does not strictly
begin at $0$ for the convolutions of $xf(x/z,Q^{2})$ with the splitting
functions, which are instead defined in the MMHT framework only to a finite
level of precision (defined at a lower bound of $x\sim 10^{-12}$). Therefore,
we do not consider any of the effects described above as serious invalidations
of the parton model, even with the full spectrum of effects due to QED
included.

02$\times 10^{-5}$4$\times 10^{-5}$6$\times 10^{-5}$8$\times
10^{-5}$0.00010.00012

1

10

100

$10^{3}$

$10^{4}$

$10^{5}$

$10^{6}$

$10^{7}$

$10^{8}$

$10^{9}$

Momentum change

$Q^{2}$Induced change in total momentum carried by the partonsAll
includedProton Mass in splitting$\delta\gamma^{(el)}$Renormalon/HT

Figure 8: The absolute change induced in the total momentum carried by the
partons due individually to the target mass corrections,
$\mathcal{O}(m_{p}^{2}/Q^{2})$, the inclusion of $\delta\gamma^{(el)}$ in eq.
(33) (elastic contributions above input) and Higher Twist/Renormalon
contributions to $\gamma^{(inel)}$.

##  3 QED Neutron PDFs

Neutron PDFs are necessary for interpreting the results of deuterium
scattering experiments that are still widely used in PDF fits to constrain the
flavour decomposition of the proton. The most widely adopted approach is to
assume isospin symmetry between hadrons in the valence distributions:

$u_{V,(p)}(x,Q^{2})=u_{(p)}(x,Q^{2})-\bar{u}_{(p)}(x,Q^{2})=d_{V,(n)}(x,Q^{2})=d_{(n)}(x,Q^{2})-\bar{d}_{(n)}(x,Q^{2}),$
(36)
$d_{V,(p)}(x,Q^{2})=d_{(p)}(x,Q^{2})-\bar{d}_{(p)}(x,Q^{2})=u_{V,(n)}(x,Q^{2})=u_{(n)}(x,Q^{2})-\bar{u}_{(n)}(x,Q^{2}),$
(37)

where the subscripts {(p),(n)} denote the proton and neutron respectively. In
practice, this is seen to produce a good agreement with the observed data and
is well motivated by the $SU(n_{f})$ flavour symmetry of QCD as well as the
fact that the evolution treats both quark flavours as essentially massless
($m_{u}^{2}/Q^{2},m_{d}^{2}/Q^{2}\sim 0$ for $Q>1$ GeV). However, as discussed
in Section 2.3, QED splitting kernels such as those in eqs. (23), (24) no
longer uphold this symmetry and are expected to generate $\mathcal{O}(\alpha)$
violations in the above relations. Therefore, to relate the distributions of
the proton to those of the neutron in a manner consistent with QED evolution,
one needs to carefully account for the effects of the relevant quark charges
$e_{u},e_{d}$ in the evolution and to allow for small amounts of isospin
violation to be introduced.

###  3.1 Modified DGLAP Evolution

QED corrections automatically result in isospin violating effects such that at
given $x$ and $Q^{2}$ values, the valence distributions can no longer be
related to one another by eqs. (36), (37). However, any modification to these
relations must still preserve the flavour quantum numbers of the proton and
neutron via the usual sum rules

$\begin{split}\int_{0}^{1}\mathrm{d}x\,u_{V,(p)}(x)=\int_{0}^{1}\mathrm{d}x\,d_{V,(n)}(x)=2\\\
\int_{0}^{1}\mathrm{d}x\,d_{V,(p)}(x)=\int_{0}^{1}\mathrm{d}x\,u_{V,(n)}(x)=1.\end{split}$
(38)

Above input, one can in principle keep track of all contributions to the
quarks that arise from QED splittings. In the case of the valence
distributions, the evolution is governed by eq. (27), where the splitting
kernels are separated into QED and QCD contributions via

$P_{i,j}=P_{i,j}^{(QCD)}+P_{i,j}^{(QED)}\;.$ (39)

Therefore, one can distinguish between two contributions to the valence
distributions in the proton (which we refer to as $q_{V}$ in the following
discussion):

$q_{V}(x,Q^{2})=q^{(QCD)}_{V}(x,Q^{2})+q^{(QED)}_{V}(x,Q^{2}),$ (40)

where $q^{(QED)}_{V}$ is defined as:

$q^{(QED)}_{V}(x,Q^{2})=\int_{Q_{0}^{2}}^{Q^{2}}\frac{d\mu^{2}}{\mu^{2}}\frac{\alpha(\mu^{2})}{2\pi}\Big{(}P_{q_{i}}^{-(QED)}\otimes
q_{V}\big{(}\frac{x}{z},\mu^{2}\big{)}\Big{)},$ (41)

and the integrand contains all QED splitting contributions to the valence
distributions. Note that an implicit overall factor of the quark electric
charge, $e_{q_{i}}^{2}$, is contained in $P_{q_{i}}^{-(QED)}$.

To parameterise the isospin violating components between the proton and the
neutron, we define:

$\begin{split}\Delta
d_{V,(n)}(x,Q^{2})=d_{V,(n)}(x,Q^{2})-u_{V,(p)}(x,Q^{2})\\\ \Delta
u_{V,(n)}(x,Q^{2})=u_{V,(n)}(x,Q^{2})-d_{V,(p)}(x,Q^{2}),\end{split}$ (42)

where naïve pointwise isospin conservation would lead both of these
expressions to evaluate to $0$. For isospin violation generated by QED
splittings, we assume that

$\Delta d_{V,(n)}(x,Q^{2})\propto u_{V,(p)}^{(QED)}(x,Q^{2}),\hskip
28.45274pt\Delta u_{V,(n)}(x,Q^{2})\propto d_{V,(p)}^{(QED)}(x,Q^{2}).$ (43)

In particular, we assume that provided that the momentum and number
conservation rules (eqs. (14), (38)) are obeyed by the constant of
proportionality, the only further step needed in relating the valence
distributions of the proton to that of the neutron is the charge re-weighting
of the relevant valence distributions, $q_{V,(p)}^{(QED)}$, to correct for
charge proportional terms in the evolution. Then, we may rewrite eq. (42) in
the form of the following equations:

$\Delta
d_{V,(n)}(x,Q_{0}^{2})=\epsilon\Big{(}1-\frac{e_{d}^{2}}{e_{u}^{2}}\Big{)}u_{V,(p)}^{(QED)}(x,Q_{0}^{2}),$
(44) $\Delta
u_{V,(n)}(x,Q_{0}^{2})=\epsilon\Big{(}1-\frac{e_{u}^{2}}{e_{d}^{2}}\Big{)}d_{V,(p)}^{(QED)}(x,Q_{0}^{2}).$
(45)

where $\epsilon$ is fixed to conserve momentum at input.

In order to satisfy momentum conservation, eq. (14), at input for the neutron,
one needs the neutron photon distribution at input. This defines the constant
of proportionality, $\epsilon$, by:

$\epsilon=\frac{\int_{0}^{1}dxx(\gamma_{(p)}(x)-\gamma_{(n)}(x))}{\int_{0}^{1}dxx(\frac{3}{4}u_{V,(p)}^{(QED)}(x)-3d_{V,(p)}^{(QED)}(x))}$
(46)

where all the distributions are evaluated at $Q_{0}^{2}=1$ GeV2. This follows
a procedure similar to that adopted in [10].

This expression implicitly depends on the assumption that the remaining
partons are then related to one another in the standard manner, assuming that
the antiquark (or sea) distributions are still well approximated by

$(\bar{u})_{(n)}(x,Q_{0}^{2})=(\bar{d})_{(p)}(x,Q_{0}^{2}),\hskip
28.45274pt(\bar{d})_{(n)}(x,Q_{0}^{2})=(\bar{u})_{(p)}(x,Q_{0}^{2})$ (47)

with all other quark flavours and the gluon being related identically between
hadrons.

Using eqs. (44) and (45) the $u$ and $d$ singlet distributions are then
related to one another between hadrons by:

$(d+\bar{d})_{(n)}(x,Q^{2})=(u+\bar{u})_{(p)}(x,Q^{2})+\Delta
d_{V,(n)}(x,Q^{2})$ (48)
$(u+\bar{u})_{(n)}(x,Q^{2})=(d+\bar{d})_{(p)}(x,Q^{2})+\Delta
u_{V,(n)}(x,Q^{2}),$ (49)

where $\Delta\\{d,u\\}_{V,(n)}$ are as defined above. Though of less apparent
interest, these relations pertain to the discussion in Section 3.2, where the
neutron photon PDF is considered as primary sensitive to distributions of the
type $q+\bar{q}$ during the evolution. In anticipation of this, we note that
$\Delta\\{d,u\\}_{V,(n)}$ lead to differences between the isospin related $u$
and $d$ singlet distributions between hadrons of only $\mathcal{O}(1\%)$,
since the $\Delta q_{V}$ terms are proportional to the contributions to the
valence quarks that arise solely from QED evolution, which are
$\mathcal{O}(\alpha)$ suppressed. In practice, relating these distributions to
one another by isospin symmetry still remains a good approximation. This will
underpin our development of a photon PDF of the neutron in the next section.

0.990.99511.0051.011.0151.021.0251.031.0351e-050.00010.0010.010.1

$\frac{d_{v,(n)}}{u_{v,(p)}}(x,Q^{2})$, $Q^{2}=1$ GeV2

$x$Comparisons of valence distributions in the Proton and Neutron

0.820.840.860.880.90.920.940.960.9811.021e-050.00010.0010.010.1

$\frac{u_{v,(n)}}{d_{v,(p)}}(x,Q^{2})$, $Q^{2}=1$ GeV2

$x$Comparisons of valence distributions in the Proton and Neutron

Figure 9: The ratio of valence quarks, related to one another by isospin, of
the neutron to that of the proton at the input scale $Q_{0}^{2}=1$ GeV2. On
the left is $u_{V,(n)}/d_{V,(p)}$, and on the right is $d_{V,(n)}/u_{V,(p)}$,
both as functions of $x$.

For the valence distributions, in practice, the magnitude of isospin violation
is seen to be a few percent, becoming significant especially at low and high
$x$, where all distributions tend towards $0$, as shown in Fig. 9. Of note is
the fact that the discrepancy between the predicted ratio of valence quarks
and the naïve isospin assumption remains at the $\sim 1\%$ level, even for the
peak of the valence distributions (at $x\sim\frac{1}{3}$,
$x\sim\frac{2}{3})$). This effect is seen to increase during the evolution,
with differences of $\sim 5\%$ at $Q=100$ GeV2.

Finally, although the primarily interest in this paper for the development of
QED corrected neutron PDFs is to provide a manner of relating the PDFs to
deuterium scattering experiments used to constrain the partons, we also wish
to highlight the potential relevance of this set in the determination of
nuclear PDFs. In particular, the assumption made in modern determinations of
nuclear PDFs (such as those of EPPS [45] and nCTEQ [46]) is to fit to data
with the assumption that the $u$ and $d$ quark type distributions in the
neutron and proton are related to one another by isospin symmetry. With the
development of this set, we propose that this assumption need not be applied
strictly and that with the introduction of QED effects, the small amounts of
isospin violation shown in Fig. 9 may be of relevance when the determination
of nuclear PDFs reach the $\mathcal{O}(5\%)$ level. While current
determinations do not reach this level of precision, a QED corrected
relationship between proton and neutron PDFs may provide better fits to the
available data, and is of interest given that recent work has begun to adopt
quark flavour dependence in fits [47].

###  3.2 The Photon PDF of the Neutron

As in the case of the proton, there is also a corresponding photon PDF of the
neutron, $\gamma_{(n)}(x,Q^{2})$, which should in general be included. At
input, the expression for this is adapted from that of the proton, eq. (4),
with the proton mass replaced by that of the neutron and the relevant form
factors substituted or approximated in the manner discussed below.

As in the case of the proton, the input distribution is due to both inelastic
and elastic photon emission.

1e-081e-071e-061e-050.00010.0010.010.11e-050.00010.0010.010.1

$x\gamma(x,Q^{2})$, $Q=1.0$ GeV

$x$Comparisons of the Elastic and Inelastic Neutron Photon
PDFElasticInelasticTotal Photon

Figure 10: The Elastic and Inelastic Photon components at $Q_{0}^{2}=1$ GeV2.

The neutron elastic distribution $\gamma_{(n)}^{(el)}$ at input is given in
terms of the Sachs form factors of the neutron, $G_{E,(n)},G_{M,(n)}$. We
adopt the phenomenological Galster parameterisation [48]:

$G_{E,(n)}=\frac{A\tau}{1+B\tau}G_{D}(Q^{2}),$ (50)

where $\tau=Q^{2}/4m_{n}^{2}$ and $G_{D}(Q^{2})$ is the dipole form factor for
hadrons (in the form commonly used to approximate $G_{E,(p)}$ when multiplied
by the proton’s magnetic moment, $G_{E,(p)}=\mu_{p}G_{D}(Q^{2})$):

$G_{D}(Q^{2})=\frac{1}{(1+\frac{Q^{2}}{\Lambda^{2}})^{2}},$ (51)

with $\Lambda^{2}=0.71$ GeV2. Values for $A$ and $B$ are then taken from a fit
to deuterium and ${}^{3}He$ scattering experiments provided by [49], for which

$A=1.70\pm 0.04,\quad B=3.30\pm 0.32\;.$ (52)

For $G_{M,(n)}$ meanwhile, a simple dipole approximation of the form:

$G_{M,(n)}=\mu_{n}G_{D}(Q^{2})\;,$ (53)

is used. These are found [50] to give a reasonably good fit to data provided
from deuterium scattering experiments.

Due to the net neutral charge of the neutron, both form factors are
significantly smaller in magnitude than those of the proton, and one therefore
expects the relevant elastic contribution to $\gamma_{(n)}$ at input to be
significantly smaller. In fact, as seen in Figs. 10 and 11 it is found to
scarcely contribute at all, comprising $\mathcal{O}(1\%)$ of the total photon
over a large range of $x$, becoming significant only at $x\sim 0.5$, where the
magnitude of the PDF itself is of vanishing importance. Therefore, given the
uncertainties associated with both models adopted for both $G_{E,(p)}$ and
$G_{M,(p)}$, $\gamma^{(el)}_{(n)}$ may reasonably be omitted for
phenomenological purposes. This is even more true for contributions above
input, since further elastic contributions are attenuated as
$1/Q^{2}\rightarrow 0$ such that $\gamma_{(n)}^{(el)}/\gamma_{(n)}\rightarrow
0$.

0.0010.010.111e-050.00010.0010.010.1

$\gamma_{(n)}^{(el,inel)}/\gamma_{(n)}^{Total}(x,Q^{2})$, $Q=1$ GeV

$x$Propotions of the Elastic and Inelastic Neutron Photon PDF
contributionsElasticInelastic

Figure 11: The relative proportions of the Elastic and Inelastic Photon
component contributions to the total Neutron Photon PDF at $Q_{0}^{2}=1$ GeV2

For the proton $F_{2}^{(inel)}$ can be divided into resonance and continuum
contributions. In the resonance region, the fit provided by CLAS for
$F_{2}^{(inel)}$ is also given for the neutron, and so can be
straightforwardly applied here. For the continuum region however, the HERMES
fit for $F_{2}$ is provided solely for the proton. Therefore, for
$F_{2,(n)}^{(inel)}$ we instead relate this approximately to the proton case.
In particular, we re-weight the proton continuum contribution according to

$F_{2,(n)}^{(inel)}=r_{F_{2}}\times F_{2,(p)}^{(inel)}\;,$ (54)

where $r_{F_{2}}$ is the ratio of the charge weighted singlet of partons
$\Sigma$, at input, for the neutron to that of the proton:

$r_{F_{2}}=\frac{4(d+\bar{d})+(u+\bar{u})+(s+\bar{s})}{4(u+\bar{u})+(d+\bar{d})+(s+\bar{s})}.$
(55)

It should be noted that in the expression above, all the distributions refer
to those of the proton, where we have used the assumption (which as discussed
in the previous section holds to a high degree of accuracy) that

$(u+\bar{u})_{(n)}=(d+\bar{d})_{(p)},\quad(d+\bar{d})_{(n)}=(u+\bar{u})_{(p)}.$
(56)

Note that any attempt to improve the accuracy of the expression in eq. (55) by
using eqs. (48), (49) would not be feasible in the current framework since
those equations depend on the parameter $\epsilon$ from the previous section,
which in turn is determined from $\gamma_{(n)}^{(inel)}$ itself.

By approximating the ratio of structure functions between the hadrons by their
respective quark singlets, the form of $F_{2,(n)}^{(inel)}$ substituted in eq.
(4) for $\gamma_{(n)}^{(inel)}$ in the continuum region at input is simply
given as $F_{2,(n)}^{(inel)}=r_{F_{2}}\times F_{2,(p)}^{(inel)}$. In Fig. 12,
one sees a broad correspondence between $r_{F_{2}}$ and the ratio
$\gamma^{(inel)}_{(n)}/\gamma^{(inel)}_{(p)}(x,Q_{0}^{2})$, particularly as
$x\rightarrow 0$, where the continuum region is dominant. Some discrepancy
between the two plots exists due to the presence of the resonance region
contribution, which as stated above is reformulated based on available neutron
data, rather than being re-scaled by $r_{F_{2}}$. The general similarly
however persists because at low $x$, the behaviour of the light quark singlets
are dominated by the sea quarks, and $u\simeq\bar{u}\simeq d\simeq\bar{d}$,
such that the effect of swapping flavours via isospin leaves the PDFs roughly
invariant in this region.

0.50.550.60.650.70.750.80.850.90.9511e-050.00010.0010.010.1

$r_{F_{2}}(x,Q^{2})$, $Q^{2}=1$ GeV2

$x$$r_{F_{2}}$

0.550.60.650.70.750.80.850.90.951e-050.00010.0010.010.1

$\frac{\gamma^{(inel)}_{(n)}}{\gamma^{(inel)}_{(p)}}(x,Q^{2})$

$x$Ratio of inelastic neutron and proton photon PDFs, $Q^{2}=1$ GeV2

Figure 12: (Left) the ratio of the charge weighted light quark singlets
between the neutron and proton. (Right) the ratio of $x\gamma^{(inel)}$
between the neutron and proton, for comparison.

Above input, $\gamma_{(n)}$ is approximated from the evolution of
$\gamma_{(p)}$ in a manner analogous to that of the quarks as described in
Section 3.1. We also distinguish between the flavour of the quark whose
splitting leads to the evolution of the photon. One can label the
contributions to $\gamma$ from the originating quark or antiquark flavour to
obtain $\gamma_{q}$, given by the following expression

$\gamma(x,\mu^{2})_{q}=\int_{Q_{0}^{2}}^{\mu^{2}}\frac{\alpha(Q^{2})}{2\pi}\frac{dQ^{2}}{Q^{2}}\int_{x}^{1}\frac{dz}{z}\Big{(}P_{\gamma,q}(z)q^{+}(\frac{x}{z},Q^{2})\Big{)}.$
(57)

where the $+$ superscript once again denotes singlet type distributions of the
form $q+\bar{q}$.

Assuming isospin symmetry, which as shown in the previous section holds to a
good approximation, one can make assumptions based on the predicted splittings
in the neutron evolution to re-weight the contributions of each $\gamma_{q}$
of the proton, based on the scheme laid out in eq. (47) to obtain:

$\begin{split}\gamma(x,Q^{2})_{(n)}^{(inel)}=\frac{e_{d}^{2}}{e_{u}^{2}}\gamma_{u,(p)}(x,Q^{2})+\frac{e_{u}^{2}}{e_{d}^{2}}\gamma_{d,(p)}(x,Q^{2})+\gamma_{\\{s,c,b,g\\},(p)}(x,Q^{2})\;,\end{split}$
(58)

where the final term accounts for all other flavours, whose contributions are
assumed to be identical for the neutron and the proton.

At the level of approximation adopted, the expression given above is expected
to be accurate to $\mathcal{O}(\alpha)$, with errors of
$\mathcal{O}(\alpha\alpha_{S}+\alpha^{2})$. Anticipating results from the next
section, it is seen that these higher orders induce changes in the resultant
photon of $\sim 3\%$ at high $x$, while the uncertainties on the CLAS fit and
the PDFs themselves each introduce a $\sim 1\%$ uncertainty on the photon PDF
at low and high $x$ respectively. Therefore, one can conservatively estimate
the uncertainty of the photon PDF of the neutron to be $\mathcal{O}(5\%)$ at
high $x$ and $\mathcal{O}(2-3\%)$ at low $x$ where the PDF and higher order
uncertainties dominate.

00.010.020.030.040.050.060.071e-050.00010.0010.010.1

$x\gamma_{(n)}(x,Q^{2})$, $Q^{2}=1.0$ GeV2

$x$ProtonNeutron

00.050.10.150.20.250.30.350.40.450.51e-050.00010.0010.010.1

$x\gamma(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$ProtonNeutron

Figure 13: A comparison of the total Neutron and Photon PDFs at $Q=1$ GeV
(left) and $Q=100$ GeV (right).

As seen in Fig. 13, at the input scale the photon PDF in the neutron is a
factor of $\sim 2$ smaller than in the proton case, while for $Q=100^{2}$ GeV2
the PDFs are comparable in size. This is as expected, since the ratio of
charges used to re-weight the proton contributions are $\mathcal{O}(1)$, and
as $\gamma^{(el)}_{(p)}$ becomes less significant in the evolution, as seen in
Fig. 4, the inelastic contribution dominates.

0.40.50.60.70.80.911e-050.00010.0010.010.1$x$Ratios at $Q^{2}=100$ GeV2Charge-
weighted Singlet Ratio $\gamma^{(inel)}_{(n)}/\gamma^{(inel)}_{(p)}$

Figure 14: The upper (red) curve is the ratio of
$\sum_{i}e_{q_{i}}^{2}(q_{i}+\bar{q}_{i})$, that is, the charge weighted sum
of quark singlets, in the neutron to that of the proton. The lower (green)
curve is the ratio of inelastic photon PDFs between the neutron and photon. At
low $x$, these are both seen to tend towards unity as the flavour invariance
of the sea (which obeys isospin symmetry maximally) dominates.

This is seen in Fig. 14, which shows the ratio of the charged-weighted quark
singlets ($\Sigma_{C}$) between the proton and neutron, and the ratio of
$\gamma^{(inel)}_{(n)}/\gamma^{(inel)}_{(p)}(x,Q^{2})$ at the same scale. As
shown above for the input, the isospin invariance demonstrated at low $x$ in
the sea quarks means that the valence properties of the hadrons are less
relevant at higher scales, leading to a photon PDF of the neutron that is
comparable to that of the proton.

##  4 Results

We now discuss the effect of adding QED corrections to the global PDF
analysis. First, in subsection 4.1, we present the changes to the PDFs due to
including the QED corrections into the input and the evolution and we show the
proton PDFs obtained from the new global analysis. In subsection 4.2.1 we
discuss the global fit quality and in 4.2.2 we present the photon PDF and
compare with other contemporary analyses. Also, in subsection 4.2.3 we show
the QED corrected structure functions. In subsection 4.2.4 We briefly outline
the impact of QED corrections on the best-fit value of $\alpha_{S}$.
Subsection 4.2.5 then finishes with a presentation of the photon-photon
luminosities in $pp$ collisions. In subsection 4.3 we quantify the
uncertainties in our determination of the photon PDF of the proton.

###  4.1 Changes to PDFs due to QED corrections

Here we show the changes in the parton distributions that are produced as a
result of the changes given in the preceding sections. We include the
$\mathcal{O}(\alpha)$, $\mathcal{O}(\alpha\alpha_{S})$ and
$\mathcal{O}(\alpha^{2})$ QED corrections, unless otherwise stated, and
compare against the baseline PDF set without QED effects described in Section
2.1.

In Fig. 15 (left) we present the percentage change for the $u,d$ and $s$
distributions as well as the gluon when QED kernels are included, against a
default of pure QCD kernels at NNLO. The effect of QED evolution on the
quarks, prior to refitting, is relatively modest, as expected due to the
$\mathcal{O}(\alpha/\alpha_{S})$ relevance of the QED splitting kernels in
comparison to those of QCD. Although the change appears to grow at low $x$,
this is in fact an artefact of the gluon PDF parameterisation, the expression
for which is reproduced here for convenience:

$xg(x,Q_{0}^{2})=A_{g}(1-x)^{\eta_{g}}x^{\delta_{g}}\Bigg{(}\sum_{n=1}^{2}a_{g,i}T_{i}^{Ch}(y(x))\Bigg{)}+A_{g^{\prime}}(1-x)^{\eta_{g^{\prime}}}x^{\delta_{g^{\prime}}}.$
(59)

Here, the two competing contributions to this expression dominate the form of
the input distribution (and therefore subsequent effects in the evolution of
the sea) at low $x$. In particular, there is a strong correlation between the
coefficients of the first and second terms, with the former term tending to
increase the gluon at low $x$ during a fit, while the latter tends to decrease
it. A delicate balance and cancellation between these effects is seen to
provide the best fit quality. However, unlike the other parameters in this
expression, $A_{g}$ is determined solely from the requirement that the
momentum sum rule (14) be satisfied. If all other parameter values are taken
from a fit using purely QCD kernels, the extra momentum provided by
$x\gamma(x,Q_{0}^{2})$ at input is compensated by a reduction of $A_{g}$,
which diminishes the gluon contribution at low $x$. Such an effect disrupts
the delicate cancellation between the terms described above. This is seen to
reduce the overall gluon momentum during the evolution, as well as that of the
quark singlet distributions, as the latter at low $x$ are primarily driven by
DGLAP emission from the gluon. Therefore, a reduction in $g$ is expected and
observed to have a knock-on effect in the same region, as shown in Fig. 15
(left).

-5-4-3-2-101231e-050.00010.0010.010.1

$\delta f$(%)

xUDSGLUON

-10-8-6-4-202461e-050.00010.0010.010.1

$\delta f$(%)

xUDSGLUON

Figure 15: The percentage change in the $u,d,s,g$ partons at $Q=100$ GeV due
to QED evolution with (right) and without (left) refitting to data.

In Fig. 15 (right) show the effect of on the quarks of refitting, described in
more detail in Section 4.1. We can see that the exaggerated effects of the
evolution at low $x$ are compensated by the other parameters of the gluon, as
discussed above. On the other hand, the behaviour of the partons at high $x$,
which shows a small reduction in the singlet distributions are a genuine
effect due to the inclusion of the QED contribution to $P_{qq}$. In
particular, this reduction is primarily a natural consequence of the
$q\rightarrow q+\gamma$ emission, which at high $x$ has the effect of reducing
the quark singlet momenta, with corresponding increases in $x\gamma(x,Q^{2})$.

We note that although the $s$ distribution experiences a larger magnitude of
change due to QED than that of the other partons, this effect is a consequence
of the $s+\bar{s}$ distribution being less well constrained by the data, and
therefore more sensitive to the effects of refitting, rather than having an
enhanced sensitivity to the effects of QED.

In Figs. 16 \- 18 we show the ratio of the PDFs with and without QED effects,
including the corresponding PDF uncertainties. We can see that upon refitting
the singlet ($q+\bar{q}$) and gluon PDFs all lie within the PDF uncertainties
of the pure QCD fit, with the central values and uncertainties remaining only
modestly affected, with $\mathcal{O}(2\%)$ reduction for the $s+\bar{s}$
distribution, (with a slight increase in the reduction at high $x$, due to the
effect of QED splittings mentioned above). The up valence quark, $u_{V}$ and
to a lesser extent the down valence quark $d_{V}$, are most sensitive to QED
effects, with a $\mathcal{O}(2-5\%)$ change at low $x$ in their central
values, though this is relatively marginal given the large uncertainties
($\sim 20\%$) in the valence quark PDFs in this region.

0.940.960.9811.021.041.061e-050.00010.0010.010.1

$(u+\bar{u})_{QED}/(u+\bar{u})_{QCD}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$($u+\bar{u}$) QCD($u+\bar{u}$) QED refit

0.940.960.9811.021.041.061e-050.00010.0010.010.1

$(d+\bar{d})_{QED}/(d+\bar{d})_{QCD}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$($d+\bar{d}$) QCD($d+\bar{d}$) QED refit

Figure 16: The ratio of the $(u+\bar{u})$, $(d+\bar{d})$ distributions (with
uncertainties) fit with and without the effects of QED in the evolution (both
at NNLO in QCD) at $Q=100$ GeV.

0.750.80.850.90.9511.051.11.151.21e-050.00010.0010.010.1

$(s+\bar{s})_{QED}/(s+\bar{s})_{QCD}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$Ratios at $Q^{2}=10^{4}$ GeV2($s+\bar{s}$) QCD($s+\bar{s}$)S QED refit

0.940.960.9811.021.041.061e-050.00010.0010.010.1

$g_{QED}/g_{QCD}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$Gluon QCDGluon QED refit

Figure 17: The ratio of the $(s+\bar{s})$, $g$ distributions (with
uncertainties) fit with and without the effects of QED in the evolution (both
at NNLO in QCD) at $Q=100$ GeV.

0.80.911.11.21e-050.00010.0010.010.1

$(s+\bar{s})_{QED}/(s+\bar{s})_{QCD}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$($u-\bar{u}$) QCD($u-\bar{u}$) QED refit

0.80.911.11.21e-050.00010.0010.010.1

$(s+\bar{s})_{QED}/(s+\bar{s})_{QCD}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$($d-\bar{d}$) QCD($d-\bar{d}$) QED refit

Figure 18: The ratio of the $(u-\bar{u})$, $(d-\bar{d})$ distributions (with
uncertainties) fit with and without the effects of QED in the evolution (both
at NNLO in QCD) at $Q=100$ GeV.

In Fig. 19 we see the details of the momentum carried by each of the partons
as a function of $Q^{2}$ for both the proton and neutron. At input the
fractional momentum carried by the photon in the proton is $0.00196$, and this
increases to about $0.007$ at very high $Q^{2}$. In the neutron the input
figure is much smaller, i.e. $0.0003$, but the rate of increase at higher
$Q^{2}$ is comparable to the proton, though a little lower due to the dominant
radiation at high $x$ being from down quarks rather than up quarks.

0.00010.0010.010.11110100$10^{3}$$10^{4}$$10^{5}$$10^{6}$$10^{7}$$10^{8}$$10^{9}$

Fraction of total momentum

$Q^{2}$Gluon$u+\bar{u}$Mom.$d+\bar{d}$Mom.$s+\bar{s}$Mom.$c+\bar{c}$Mom.$b+\bar{b}$Mom.Photon
Momentum

0.00010.0010.010.11110100$10^{3}$$10^{4}$$10^{5}$$10^{6}$$10^{7}$$10^{8}$$10^{9}$

Fraction of total momentum

$Q^{2}$Gluon$u+\bar{u}$Mom.$d+\bar{d}$Mom.$s+\bar{s}$Mom.$c+\bar{c}$Mom.$b+\bar{b}$Mom.Photon
Momentum

Figure 19: The percentage of the momentum carried by the partons, including
the photon PDF, in the proton (left) and neutron (right) as a function of
$Q^{2}$ when QED evolution is included.

In Fig. 20 we show the effect of the higher and mixed order corrections to the
evolution on $x\gamma(x,Q^{2})$. We can see that the
$\mathcal{O}(\alpha\alpha_{S})$ and $\mathcal{O}(\alpha^{2})$ kernels are seen
to reduce the photon distribution by $\sim 1-3\%$, particularly at high $x$.
The effect induced by the $\mathcal{O}(\alpha^{2})$ kernels is of
$\mathcal{O}(0.5-1\%)$, and further changes associated with the exclusion of
yet higher orders in perturbation theory are expected to be even smaller.
Since other sources of uncertainty, discussed in Section 4.3 are somewhat
larger, it is not thought that such scale uncertainties will be significant
for the photon at the level of accuracy being discussed in this paper.

0.9650.970.9750.980.9850.990.99511.0051e-050.00010.0010.010.1

$\frac{\gamma_{\alpha+\alpha\alpha_{S}+(\alpha^{2})}}{\gamma_{\alpha}}(x,Q^{2})$,
$Q^{2}=10^{4}$ GeV22

$x$$\alpha+\alpha\alpha_{S}$$\alpha+\alpha\alpha_{S}+\alpha^{2}$

Figure 20: Ratio of the Photon PDF with and without
$\mathcal{O}(\alpha\alpha_{S})$,$\mathcal{O}(\alpha^{2})$ corrections, at
$Q^{2}=10^{4}$ GeV2.

Similarly, as shown Fig. 21, the QCD order of the DGLAP evolution is found to
have a modest effect on the resultant photon PDF produced. The photon
experiences a slight reduction for intermediate values of $x$,
$\mathcal{O}(1-2\%)$, with a slight increase at high and low $x$. This is
largely due to differences in the underlying quark singlet, which as
previously noted, have a strong role in influencing the form of
$x\gamma(x,Q^{2})$ at higher scales.

0.9750.980.9850.990.99511.0051.011.0151e-050.00010.0010.010.1

$\gamma_{NNLO}/\gamma_{NLO}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$Ratio of the Photon PDF at NLO and NNLO in QCD

Figure 21: Ratio of the Photon PDF at NLO and NNLO in QCD during DGLAP
evolution, at $Q^{2}=10^{4}$ GeV2.

###  4.2 Results of Global fits with QED corrections

The fitting procedure is broadly similar to that of MMHT14 , with the
exception of changes to the structure function fits described below, and the
inclusion of some new data as described in Section 2.1. In Section 5 we detail
the results of an alternative fit that also includes high mass Drell-Yan data,
which is seen to have some sensitivity to the inclusion of QED effects.

####  4.2.1 The quality of the global fits

Change in $\chi^{2}$ due to QED evolution compared to MMHT14+HERA I+II
---
NLO before fit |  NLO after fit |  NNLO before fit |  NNLO after fit
4180 (+41) |  4151(+12) |  3574 (+42) |  3539 (+7)
Table 1: The total $\chi^{2}$ for partons with the effects of QED, both prior
to and after refitting the parton parameters, at NLO and NNLO. Before the fit,
the parameters derived from the QCD fits described in Section 2.1 are used.
The NLO fit contains 3609 data points, while the NNLO contains a total of 3276
(since the later omit some jet data). The numbers in brackets show the change
in $\chi^{2}$ due to the inclusion of the QED corrections.

In Table 1, we provide the change in the total $\chi^{2}$ in the fit to all
data after the inclusion of all QED effects at NLO and NNLO, before and after
refitting the partons. The full breakdown between individual datasets is given
in Table 2. We can see that the change in fit quality due to both the changes
in the evolution and $\mathcal{O}(\alpha)$ corrections to the structure
functions detailed below are modest. While the effects purely driven by the
evolution naturally lead to an increase in $\chi^{2}$ (where the pure QCD fit
parameters are used), this is somewhat reduced after refitting the partons to
data with the full QED effects included. However, some increase with respect
to the original QCD fit is still observed after refitting. From Table 2 we can
see that this is primarily due to tension with the BCDMS $F_{2}$ and ZEUS CC
data, with the former responsible for a $\sim+6$ increase in the total
$\chi^{2}$ and the latter $\sim+2$. This is somewhat compensated for by a
$\sim-2$ reduction from a slightly improved fit to $F_{2}$ and $F_{3}$ data
from the NuTeV experiment, which see a mild improvement to the fit.

Data set |  $\chi^{2}/N_{pts}$ NNLO before fit (QCD) |  $\chi^{2}/N_{pts}$ NNLO after fit(QCD+QED)
---|---|---
BCDMS $\mu p$ $F_{2}$ [51] |  178 / 163 |  182 / 163 (+4)
BCDMS $\mu d$ $F_{2}$ [52] |  142 / 151 |  144 / 151 (+2)
NMC $\mu p$ $F_{2}$ [53] |  124 / 123 |  125 / 123
NMC $\mu d$ $F_{2}$ [53] |  108 / 123 |  108 / 123
NMC $\mu n/\mu p$ $F_{2}$ [54] |  128 / 148 |  127 / 148
E665 $\mu p$ $F_{2}$ [55] |  65 / 53 |  65 / 53
E665 $\mu d$ $F_{2}$ [55] |  61 / 53 |  61 / 53
SLAC $ep$ $F_{2}$ [56][57] |  31 / 37 |  31 / 37
SLAC $ed$ $F_{2}$ [56][57] |  26 / 38 |  25 / 38
NMC/BCDMS/SLAC/HERA $F_{L}$ [53, 51, 57, 58, 59, 60] |  66 / 57 |  66 / 57
E866/NuSea $pp$ DY [61] |  224 / 184 |  223 / 184
E866/NuSea $pd/pp$ DY [62] |  11 / 15 |  11 / 15
NuTeV $\nu N$ $F_{2}$ [63] |  37 / 53 |  36 / 53 (-1)
CHORUS $\nu N$ $F_{2}$ [42] |  29 / 42 |  29 / 42
NuTeV $\nu N$ $xF_{3}$ [63] |  31 / 42 |  31 / 42
CHORUS $\nu N$ $xF_{3}$ [42] |  19 / 28 |  19 / 28
CCFR $\nu N\rightarrow\mu\mu X$ [64] |  77 / 86 |  78 / 86
NuTeV $\nu N\rightarrow\mu\mu X$ [64] |  42 / 40 |  41 / 40
HERA I+II CC $e^{+}p$ [29] |  52 / 39 |  52 / 39
HERA I+II CC $e^{-}p$ [29] |  63 / 42 |  65 / 42 (+2)
HERA I+II NC $e^{+}p$ 920 GeV [29] |  510 / 402 |  510 / 402
HERA I+II NC $e^{-}p$ 920 GeV [29] |  239 / 159 |  240 / 159 (+1)
HERA I+II NC $e^{+}p$ 820 GeV [29] |  88 / 75 |  88 / 75
HERA I+II NC $e^{-}p$ 575 GeV [29] |  261 / 259 |  262 / 259
HERA I+II NC $e^{-}p$ 460 GeV [29] |  246 / 209 |  246 / 209
HERA $ep$ $F_{2}^{charm}$ [65] |  80 / 52 |  80 / 52
DØ II $p\bar{p}$ incl. jets [66] |  117 / 110 |  117 / 110
CDF II $p\bar{p}$ incl. jets [67] |  60 / 76 |  60 / 76
CDF II $W$ asm. [68] |  16 / 13 |  15 / 13
DØ II $W\rightarrow\nu e$ asym. [69] |  31 / 12 |  30 / 12
DØ II $W\rightarrow\nu\mu$ asym. [70] |  16 / 10 |  16 / 10
DØ II $Z$ rap. [71] |  17 / 28 |  17 / 28
CDF $Z$ rap. [72] |  40 / 28 |  40 / 28
ATLAS $W^{+}$, $W^{-}$, $Z$ [73] |  41 / 30 |  41 / 30
CMS $W$ asymm $p_{T}>35$ GeV [74] |  7 / 11 |  7 / 11
CMS asymm $p_{T}>$ 25 GeV, 30 GeV [75] |  8 / 24 |  8 / 24
LHCb $Z\rightarrow e^{+}e^{-}$ [76] |  22 / 9 |  22 / 9
LHCb $W$ asymm $p_{T}>20$ GeV [77] |  14 / 10 |  13 / 10
CMS $Z\rightarrow e^{+}e^{-}$ [78] |  23 / 35 |  22 / 35
ATLAS high-mass Drell-Yan [4] |  17 / 13 |  18 / 13
CMS double diff. Drell-Yan [79] |  152 / 132 |  152 / 132
Tevatron, ATLAS, CMS $\sigma_{t\bar{t}}$* [80, 81, 82, 83, 84, 85, 86] |  14 / 18 |  14 / 18
All data |  3532 / 3276 |  3539/3276 (+7)
Table 2: The $\chi^{2}$ breakdown showing $\chi^{2}/N_{pts}$ by data set for
NNLO QCD and NNLO QCD + QED PDF fits.

####  4.2.2 The photon PDF of the proton

In Fig. 22 we compare our photon PDF with those of LUXqed and NNPDF3.1luxQED,
and the agreement is found to be quite good, i.e. they are within $\sim 2\%$
over a broad range of $x$, diverging somewhat at high $x$ where uncertainties
are seen to be large, and are close to the LUXqed photon. In particular the
MMHT photon displays a very slight tendency to be somewhat larger in the
intermediate range of $x$ and predicts a somewhat smaller photon at lowest
$x$. This is to some extent explained by the fact that the charge-weighted
singlet ($\sum_{i}e_{q_{i}}^{2}(q+\bar{q})$) differs between MMHT2015qed and
those of PDF4LHC15$\\_$nnlo$\\_$100 [87] (which are the underlying partons
used for LUXqed in the higher $Q^{2}$ representation of $F_{2,L}$) and
NNPDF3.1 [36] as shown in Fig. 23. Note that the MMHT2015qed baseline PDFs are
a percent or more larger than those of MMHT2014 at low $x$ and high-$Q^{2}$ as
a consequence of fitting to the updated HERA data. This $\sim 2-4\%$ reduction
of the charge weighted singlet between the sets as compared with ours in the
range $10^{-4}<x<10^{-1}$ then leads to a reduction in the relevant photon PDF
ratios, as the evolution of $x\gamma(x,Q^{2})$ is sensitive to this
combination of partons. We also note that the largest discrepancy in the
photon PDFs between MMHT2015qed and LUXqed at $x\sim 0.5-0.6$ is also a common
feature of the charge-weighted quark difference in the relevant $x>0.5$
region.

0.90.9511.051.11.151e-050.00010.0010.010.1

$\frac{\gamma_{LUXqed}}{\gamma_{MMHTqed}}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$MMHTQEDLUXQED

0.90.9511.051.11.151e-050.00010.0010.010.1

$\frac{\gamma_{NNPDF31qed}}{\gamma_{MMHTqed}}(x,Q^{2})$, $Q^{2}=10^{4}$ GeV2

$x$MMHTQEDNNPDF

Figure 22: The ratio of Photon PDFs between the LUXqed and NNPDF3.1luxQED sets
with that of MMHT, at $Q^{2}=10^{4}$ GeV2.

Another reason why we anticipate that the $x\gamma(x,Q^{2})$ as outlined in
this work may be somewhat greater in value, in an intermediate range in $x$,
compared to that of LUXqed is due to the exclusion of lepton splitting
contributions in our DGLAP evolution, which are included in the evolution used
to develop the LUXqed set. In Section 2.3 we explicitly neglected the sum over
lepton charges in Eq. 25. In general, since $\gamma\rightarrow l\bar{l}$
splittings should reduce the photon distribution (nearly uniformly since it
occurs as a coefficient to $\delta(1-x)$ in $P_{\gamma\gamma}$), one expects
that excluding this term should lead to a somewhat increased photon. To
estimate the effect of including this term, in Fig. 24 we draw a comparison to
$x\gamma(x,Q^{2})$ evolved with $\mathcal{O}(\alpha)$ lepton splittings
included in evolution and as anticipated find that this term does lead to a
$\mathcal{O}(1-2\%)$ reduction, which becomes more pronounced at higher
$Q^{2}$. Along with the ratio of the charged singlet used in the evolution,
neglecting lepton splittings111Note that excluding this term is still a
reasonable approximation given that a fully consistent treatment with a
coupled DGLAP evolution would require the development of lepton PDF
distributions which as discussed in [37] are found to have a negligible impact
on the evolution of the PDFs on the whole. leads to an independent source of
enhancement for our $x\gamma(x,Q^{2})$, further accounting for the difference
seen in Fig. 22.

Common to all the sets are errors of $\mathcal{O}(1\%)$, displaying the
remarkable improvements in accuracy seen in photon PDFs developed on the
strategy outlined in this paper and that of [15] and [21], in comparison to
that of older sets. A full breakdown of the contributing sources of error are
explored in Section 4.3.

0.80.850.90.9511.051.11.151.21.251e-050.00010.0010.010.1$x$Ratios at
$Q^{2}=10^{4}$ GeV2Singlet Ratio LUXSinglet Ratio NNPDFMMHTQED

Figure 23: The ratios of the charged singlet,
$\sum_{i}e_{q_{i}}^{2}(q+\bar{q})$, between the LUXqed (which in turn adopts
the quark and antiquark PDFs of PDF4LHC15) and NNPDF3.1LUXqed, against that of
our set.

0.980.9820.9840.9860.9880.990.9920.9940.9960.99811e-050.00010.0010.010.1

$\frac{\gamma_{l\bar{l}}}{\gamma}$, $Q^{2}=100,10^{4}$ GeV2

$x$Ratio of the Photon PDF with and without lepton contributions in
$P_{\gamma\gamma}$$Q^{2}=100$GeV2$Q^{2}=10^{4}$GeV2

Figure 24: The ratios of the photon with and without the sum over lepton
charge $\sum e_{l}^{2}$ contribution (the right hand side of eq. 25) in the
${\cal O}(\alpha)$ contribution to $P_{\gamma\gamma}$.

####  4.2.3 QED Corrected Structure Functions

The PDFs are related to the measured structure functions by the standard
formulae

$\begin{split}F_{i}(x,Q^{2})=x\sum_{q,\bar{q}}e_{q}^{2}\int_{x}^{1}\frac{dz}{z}q(z,Q^{2})\Big{\\{}\delta\left(1-\frac{x}{z}\right)+\frac{\alpha_{S}}{2\pi}C_{i,q}\left(\frac{x}{z}\right)+...\Big{\\}}+\\\
x\sum_{q,\bar{q}}e_{q}^{2}\int_{x}^{1}\frac{dz}{z}g(z,Q^{2})\Big{\\{}\frac{\alpha_{S}}{2\pi}C_{i,g}\left(\frac{x}{z}\right)+...\Big{\\}}.\end{split}$
(60)

where $i$ ($=2,3,L...$) labels the structure function and $C_{q,g}$ are the
corresponding coefficient functions. The introduction of a photon PDF and of
QED corrections to the DGLAP splitting kernels requires that we also modify
the expression for these to include $\mathcal{O}(\alpha)$ corrections, in
particular introducing terms of the form
$C_{\gamma}^{(\alpha)}\otimes\gamma(z,Q^{2})$, for both Neutral Current (NC)
and Charged Current (CC) processes.

0.980.9850.990.99511.0051.010.0010.010.1

$F_{NNLOQCD+QED}(x,Q^{2})/F_{NNLOQCD}(x,Q^{2})$

$x$Structure Function (F) ratios with and without
QED$F_{2}$CC$F_{2}$NC$F_{3}$CC

0.980.9850.990.99511.0051.010.0010.010.1

$F_{NNLOQCD+QED}(x,Q^{2})/F_{NNLOQCD}(x,Q^{2})$

$x$Structure Function (F) ratios with and without QED
(refitted)$F_{2}$CC$F_{2}$NC$F_{3}$CC

Figure 25: The ratio of the Charged and Neutral Current $F_{2}$ and Charged
Current $xF3$ for the proton, with and without the effects of QED, both at
$Q^{2}=10^{4}$ GeV2. (Left) the effects of naïve inclusion of QED splitting
kernels without the refitting of the partons (in which the artificial
reduction in the low $x$ gluon and hence the sea quarks has an enhanced
effect, as discussed in the text). (Right) The ratio of Structure Functions
after refitting the partons, with modest effects observed in $F_{2}$ CC and
NC.

In Fig. 25 we show the effect of these changes with and without refitting.
Again, the sensitivity introduced by the gluon parameterisation is seen to
have an effect at low $x$, reducing $F_{2,3}$ somewhat, while after fitting,
the CC structure functions $F_{2,3}$ are moderately decreased at low $x$. In
the NC case however, $F_{2}$ is generally reduced by $\mathcal{O}(0.5\%)$, as
anticipated by the fact that the introduction of QED in the evolution is seen
in general to diminish the quark singlet content, see Fig. 15.

####  4.2.4 Effects of QED on $\alpha_{S}$ determination in the global PDF fit

In addition to the fit described above, we have also performed a simultaneous
fit to the strong coupling , $\alpha_{S}(M_{Z})$. The value typically used
during the evolution and the comparison to data is taken as a fixed value
$\alpha_{S}(M_{Z})=0.118$, which reflects a combination of both the best fit
value exclusively from our fit to data, and the independent inclusion of the
world average of $\alpha_{S}(M_{Z})=0.1181\pm 0.0011$ [88], as discussed in
Section 5.1 of [27].

In principle, one might expect that the value of $\alpha_{S}(M_{Z})$ found
after refitting with the effects of QED included will be somewhat less than
that in a pure QCD fit. This is because at leading order, the effect on the
$q+\bar{q}$ distributions during the evolution, particularly at high $x$, is
due to gluon emission, $q\rightarrow qg$, which leads to a slight reduction of
the singlet. In a pure QCD fit, the parameters that provide the best fit are a
combination of both the input distribution and a value of $\alpha_{S}(M_{Z})$
which drives gluon emission at a rate (determined by $P_{qq}^{(QCD)}$) in the
evolution such that the PDFs at higher scales are best fit to the data.

At LO in QED however, the electromagnetic coupling $\alpha$ plays virtually
the same role in the evolution of the singlet distributions, diminishing the
high $x$ content due to photon emissions $q\rightarrow q\gamma$. Therefore at
LO, one can consider the inclusion of QED as an enhancement to $P_{qq}$ with
an increased effecting coupling:

$\alpha_{S}\rightarrow\alpha^{\prime}=\Big{(}\alpha_{S}+\frac{e_{q}^{2}\alpha}{C_{F}}\Big{)}.$
(61)

In a fit that includes the coupling constants as free parameters, one expects
that $\alpha^{\prime}$, rather than $\alpha_{S}$ would tend towards a value
that best models the loss of the singlet during evolution to emission (whether
to a photon or gluon). Since $\alpha_{S}$ is the only free parameter in the
fit (where we adopt the world best measurement value for $\alpha$ [88]), one
naturally expects the best fit value for $\alpha_{S}$ to be reduced to
accommodate the modification in eq. (61). Naïvely, one may expect the
magnitude of this reduction to compensate for the magnitude of the
modification term $e_{q}^{2}\alpha/C_{F}\sim 10^{-3}$. Though small, this is
similar to the global fit uncertainty on $\alpha_{S}$, and the effects of QED
may therefore be significant in its determination.

This was also investigated in the development of the original MRST QED set
[10], where it was found that despite the above considerations, between the
pure QCD and QCD+QED fit, $\alpha_{S}(M_{Z})$ remained essentially unchanged.
The reason found for this was that the fit (especially the NMC and HERA data)
preferred a larger value for the gluon at low $x$, which is sensitive to
$\alpha_{S}(M_{Z})$ since $dF_{2}/d\ln Q^{2}\propto\alpha_{S}P_{qg}\otimes
g(Q^{2})$. However, the momentum carried by the photon detracts from that
carried by the small-$x$ gluon and as a result, the change to the gluon at
small $x$ has a tendency to require a larger value of $\alpha_{S}(M_{Z})$ than
would otherwise be obtained. This pulls in a direction opposite to the
reduction of $\alpha_{S}(M_{Z})$ as described above, and reduces the magnitude
by which one might anticipate a change after refitting with the effects of
QED.

With the updated QED parton framework, we find that $\alpha_{S}(M_{Z})$
experiences a reduction from 0.1181 in the pure NNLO QCD case to 0.1180 in the
fit with QED, while at NLO the result is unchanged within the numerical
precision of the fit. Although at NNLO this does represent a small reduction,
in neither case is allowing $\alpha_{S}$ to be free seen to improve the total
fits by any significant degree, with $\Delta\chi^{2}<1$. However, in future
global fits, the inclusion of QED effects in the partons may come to be
significant as the accuracy of such measurements are improved.

####  4.2.5 Photon-photon luminosity

A sense of the relevance of the photon PDF to particle production at colliders
such as the LHC may be determined from an inspection of the $\gamma\gamma$
luminosity expected at these energies (14 TeV), shown in Fig. 26. As seen in
Fig. 22, our photon and that of other sets based on the LUXqed formulation
show good agreement, and therefore our predicted $\gamma\gamma$ luminosity,
${\rm d}L_{\gamma\gamma}/{\rm d}\ln M^{2}$, bears a strong resemblance to
others in the literature (see e.g. Fig. 19 in [22]). Also shown in Fig. 27 is
the expected luminosity for a High-Energy LHC proposal with (CoM) energy
$\sqrt{s}=$27 TeV, and a Future Circular Collider with $\sqrt{s}=$100 TeV,
where the total $\gamma\gamma$ luminosity is comparable to that of
$\Sigma_{i}(q_{i}\bar{q}_{i}+\bar{q}_{i}q_{i})$ at present LHC CoM energies
(14 TeV).

Furthermore, as our photon PDF is separable by its elastic and inelastic
components, we are able to distinguish between
$\gamma^{(inel)}\gamma^{(inel)}$ and $\gamma^{(el)}\gamma^{(el)}$
contributions to the overall luminosity. The latter is of particular interest
in the context of photon-initiated central exclusive production (CEP – see
e.g. [89, 90]). In this process the protons collide peripherally, exchanging
only photons while remaining in tact, such that they can be detected and their
kinematic properties reconstructed in dedicated proton tagging detectors
installed in association with ATLAS [91] (AFP) and CMS [92] (CT–PPS).

1$\times 10^{-7}$1$\times 10^{-6}$1$\times
10^{-5}$0.00010.0010.010.11101001000

Partonic Luminosity

M
[GeV]$\gamma\gamma$$\gamma^{(inel)}\gamma^{(inel)}$$\gamma^{(el)}\gamma^{(el)}$

Figure 26: $\gamma\gamma$ partonic luminosities as a function of invariant
mass at Centre-of-Mass energies of 14 TeV. Note that for the elastic
$\gamma^{(el)}\gamma^{(el)}$, the multi-particle interaction (MPI) effects are
not included, and their inclusion would reduce (lower) the blue curve to some
degree (discussed in text).

1$\times 10^{-7}$1$\times 10^{-6}$1$\times
10^{-5}$0.00010.0010.010.111101001000

Partonic Luminosity

M [GeV]$\gamma\gamma$(14 TeV)$\gamma\gamma$(27 TeV)$\gamma\gamma$(100 TeV)

Figure 27: $\gamma\gamma$ partonic luminosities as a function of invariant
mass at Centre-of-Mass energies of 14, 27 and 100 TeV.

The cross section for this CEP process can be calculated within the so–called
Equivalent Photon Approximation [16], in which the photon flux associated with
the colliding beam of charged particles may be expressed in terms of the
elastic structure functions $F_{2,L}^{(el)}$, in a manner similar to that
considered in this paper. The $\gamma^{(el)}\gamma^{(el)}$ luminosity,
represented in Fig. 26, corresponds to precisely the luminosity that could be
delivered in this approach.

However, this interpretation must be qualified with an important caveat, which
is that for an exclusive production process, where both protons remain intact
after scattering, one needs to multiply the final result obtained from the
naïve use of $\gamma^{(el)}$ as an incoming parton by a ‘soft survival’
factor, corresponding to the probability of no additional particle production
due to multi-particle interactions (MPI) [93]. Furthermore, the luminosities
shown in Fig. 26 can not directly be applied to the calculation of cross
sections for more exclusive final states, such as when explicit cuts are
placed on the presence of additional tracks within the central portion of the
detector, but require suitable modification as in [7].

We conclude this section with a discussion of the effect of higher and mixed
orders of QED, $\mathcal{O}(\alpha\alpha_{S})$ and $\mathcal{O}(\alpha^{2})$
during the evolution and the significance of their impact on the total
luminosity at present CoM energies at the LHC. As previously observed in Fig.
20, the inclusion of these higher order splitting functions in the evolution
of $x\gamma(x,Q^{2})$ have a tendency to reduce its magnitude, particularly at
the higher range in $x$. In Fig. 28, the proportional effects of such changes
in $dL_{\gamma\gamma}/d\ln M^{2}$ are shown. We see that, above the
electroweak and near TeV scales, the importance of these higher orders become
significant, inducing a $\mathcal{O}(5\%)$ reduction in the total
$\gamma\gamma$ luminosity.

0.930.940.950.960.970.980.991110100100010000M [GeV]Ratio of total
$\gamma\gamma$ luminosity with and without
$\mathcal{O}(\alpha\alpha_{S}(+\alpha^{2}))$$\mathcal{O}(\alpha\alpha_{S})$included$\mathcal{O}(\alpha\alpha_{S}+\alpha^{2})$included

Figure 28: Ratio of $\gamma\gamma$ partonic luminosity as calculated from
$x\gamma(x,Q^{2})$ with and without $\mathcal{O}(\alpha\alpha_{S})$ and
$\mathcal{O}(\alpha^{2})$ DGLAP splitting kernels during evolution as a
function of invariant mass, at a proton-proton Centre-of-Mass energy of 13
TeV.

###  4.3 Uncertainties on the photon PDF

Our treatment of the contributions to the photon PDF uncertainty are in some
cases identical to LUXqed. However, as discussed in Section 2.3, due to the
lower starting scale adopted in our evolution procedure, we also include
higher twist corrections in the form of a renormalon model, for which the
undetermined coefficient $A^{\prime}_{2}$ in eq. (30) is fit to the data,
introducing an independent source of uncertainty.

For completeness, a full description of the uncertainty contributions is given
below. The size of the different sources of uncertainty as a function of $x$
and for different scales $Q^{2}$ is shown in Figs. 29-31.

00.0020.0040.0060.0080.010.0120.0140.0160.0180.021e-050.00010.0010.010.1

$\delta x\gamma/x\gamma$

$x$Photon PDF Uncertainty contributions, $Q_{0}^{2}$ = 1 GeV2Total
UncertaintyElasticR$W^{2}$ResonanceContinuumRenormalonPDFs

Figure 29: Photon PDF Uncertainty contributions (added in quadrature to give
the total uncertainty), $Q_{0}^{2}=1$ GeV2. Note that the upper $x$ range has
been restricted in this plot due to the effect of the kinematic cut given in
(10)

00.010.020.030.040.050.060.070.081e-050.00010.0010.010.1

$\delta x\gamma/x\gamma$

$x$Photon PDF Uncertainty contributions, $Q^{2}$ = 10 GeV2Total
UncertaintyElasticR$W^{2}$ResonanceContinuumRenormalonPDFs

Figure 30: Photon PDF Uncertainty contributions (added in quadrature to give
the total uncertainty), $Q^{2}=10$ GeV2.

00.010.020.030.040.050.060.071e-050.00010.0010.010.1

$\delta x\gamma/x\gamma$

$x$Photon PDF Uncertainty contributions, $Q^{2}$ = 104 GeV2Total
UncertaintyElasticR$W^{2}$ResonanceContinuumRenormalonPDFs

Figure 31: Photon PDF Uncertainty contributions (added in quadrature to give
the total uncertainty), $Q^{2}=10^{4}$ GeV2.

  * •

Elastic: The uncertainty contributions from the A1 fit for $F_{2}^{(el)}$ are
twofold. In particular, the fits provided by A1 are given in the unpolarized
and polarised forms, where the latter accounts for potential two photon
exchange (TPE) processes between the lepton probe and the proton in DIS
experiments. Following the approach of LUXqed, we use the latter for our
estimate precisely because it provides constraints on TPE. As well as the
intrinsic uncertainty provided by the A1 collaboration for this fit
$\delta(F_{2}^{(el)})_{a}$, similarly to LUX, we adopt the symmetrised
difference between the polarised and unpolarized fit as an independent source
of error, $\delta(F_{2}^{(el)})_{b}$. The total uncertainty on $F_{2}^{(el)}$
is then simply given by the sum of these two contributions in quadrature.

  * •

R: The contributions from $F_{L}$ are modelled in precisely the same manner as
that of LUXqed, using the parameterisation of the form:

$F_{L}(x,Q^{2})=F_{2}(x,Q^{2})\Big{(}1+\frac{4m_{p}^{2}x^{2}}{Q^{2}}\Big{)}\frac{R_{L/T}(x,Q^{2})}{1+R_{L/T}(x,Q^{2})},$
(62)

where $R_{L/T}=\sigma_{L}(x,Q^{2})/\sigma_{T}(x,Q^{2})$ represents the ratio
between the absorption cross sections for longitudinal and transversely
polarised photons. Our expression for this ratio is provided by the LUXqed
group, who, following the procedure used by the HERMES collaboration [33], in
turn adapt the expression from the $R_{1998}$ fit[34] provided by the E143
Collaboration for use in low $Q^{2}$ regions and assign it a conservative $\pm
50\%$ uncertainty, which we also adopt.

  * •

W2: As discussed in Section 2.2, two distinct fits for $F_{2}^{(inel)}$ are
used above (HERMES [33]) and below (CLAS[32] and Cristy-Bosted[94]) a
threshold of $W_{cut}^{2}=3.5$ GeV2. Since $W_{cut}^{2}$ is defined somewhat
arbitrarily and theoretically induces some small amount of discontinuity in
the contributions to $\gamma^{(inel)}$, we treat the cut value as an
independent source of uncertainty, varying it the region $3<W_{cut}^{2}<4$
GeV2. Even with this relatively conservative approach, the uncertainty on
$W_{cut}^{2}$ is seen to be vastly dominated by other sources.

  * •

Resonance: The uncertainty of $F_{2}^{(inel)}$ in the resonance region is
taken as the symmetrised difference between the CLAS fit, which is used as the
standard for our input, and that of the Cristy-Bosted, similar to the
procedure used by LUXqed.

  * •

Continuum: The uncertainty of $F_{2}^{(inel)}$ in the continuum region is
adapted directly from the uncertainty bands of the GDP-11 fit provided by the
HERMES collaboration. This is a different type of uncertainty estimate from
this source as that adopted by LUXqed, who vary the scale at which $F_{2}$
goes from being described by the GDP-11 fit to calculated in terms of the
PDFs. However, each estimation of uncertainty is very small.

  * •

Renormalon: For the fitting and uncertainty of the coefficient
$A^{\prime}_{2}$ in eq. (30), we implemented the original renormalon model of
[39] into the calculation of the structure functions themselves, $F_{2,3}$, as
used in the fit . $A^{\prime}_{2}$ was then varied to induce a
$\Delta\chi^{2}=\pm 10$ change in the overall fit quality of the partons (as
seen in Fig. 2 in Section 2.3), creating a generous uncertainty band of
$-0.4<A^{\prime}_{2}<-0.2$, with a best fit value of -0.3. We note that our
global fit to the data favours a renormalon contribution $\sim 50\%$ greater
than the value used in the original model by Dasgupta and Webber [39]. At high
$x$, this is seen to be a comparable source of uncertainty with that of
$\delta(F_{2}^{(el)})$. Unlike all other terms discussed so far, the
uncertainty in $A^{\prime}_{2}$ enters during the evolution, rather than at
input.

  * •

PDFs: Above the input scale $Q_{0}^{2}=1$ GeV2 the $\gamma^{(inel)}$
contributions are modelled solely from the splittings of other partons during
the DGLAP evolution. Hence, the intrinsic uncertainty on the other PDFs
propagate into the form of the photon PDF as it evolves. This reflects the
standard 50 eigenvector uncertainties associated with the fit of the free
parameters in the MMHT parameterisation (see eqs. (11) and (13)), which
generate the uncertainty bands for all flavours of parton ($q,\bar{q},g$),
naturally generating uncertainties in the photon during splittings of the form
$q\rightarrow q\gamma$ and $g\rightarrow q\bar{q}\gamma$. At low $x$, as is
the case of LUXqed, this dominates as the primary source of uncertainty.

Since our $\gamma^{(inel)}$ is evolved from a common starting scale, we are
alleviated of the consideration of matching scales between the photon and
other partons (though this is seen to be negligible even when necessary, as
shown for (M) in Fig. 15 of [22]). Furthermore, in comparison to that of
LUXqed, our set neglects certain contributions to the photon uncertainty. In
particular, rather than the Twist-4 uncertainties considered by LUXqed for
$F_{L}$ (which an inspection of (T) in Fig. 15 of [22] reveals to be
overwhelmingly dominated by other sources), our treatment of the Higher-Twist
(HT) corrections to the structure function in the form of the renormalon lead
to a more significant uncertainty at high $x$, consistent with our choice of a
lower starting scale for the evolution.

Indeed, since our starting scale is at $Q_{0}^{2}=1$ GeV, as shown in Fig. 29,
the uncertainties at input have a markedly different form to the kind that
arises during the evolution. Naturally, effects that pertain to the evolution,
(the PDF eigenvector uncertainties and the renormalon) are absent at this
scale, and the dominating effects are seen to be the uncertainty on the
resonance contribution to $F_{2}^{(inel)}$, the uncertainties on the Sachs
form factors provided by the GD-11 fit ($\delta F_{2}^{(el)})$ and the
uncertainty on $R_{L/T}$. As the evolution occurs however, the PDFs
overwhelmingly dominate as the source of uncertainty at low $x$, and in
conjunction with the uncertainty on the renormalon parameter $A^{\prime}_{2}$,
become significant contributions along with those of the Sachs form factors at
higher $x$.

It is noted that we do not account for the uncertainty that arises from the
Higher Order (HO) terms missing from the QCD components of the evolution, as
estimated in LUXqed. Although we have given an indication of the magnitude of
the change in order from QCD (from NLO to NNLO) in the evolution in Fig. 21 of
the previous section (which broadly corresponds to the (HO) band in fig. 15 of
[22]), we do not treat this difference as an independent source of
uncertainty, since PDFs have typically been provided at both NLO and NNLO in
QCD, each with independently derived uncertainty bands. Despite not being
included as a default, recent work [95] has begun to explore the possibility
of incorporating such uncertainties into the PDF fitting framework of MMHT in
a standard manner.

Overall, we note the similarity between the form of our uncertainty with
others, being less than 2% for $10^{-5}<x<0.5$, demonstrating a drastic
improvement with early photon PDF sets such as MRST2004QED [10] and NNPDF2.3
[12].

We provide the photon PDF along with the quark, antiquark and gluon PDFs in
grids which also contain all information about the uncertainties. PDF sets are
typically provided as grids in the LHAPDF6 format, with each grid representing
either the central value of the PDFs, or the PDFs at a given $\pm$ eigenvector
direction in the independent parameter space PDFs. As noted above, as well as
the uncertainties that are routinely given in such sets associated with the
non-photon PDF parameters, the set that is produced as a result of the work
described here now contains uncertainties associated with the photon
parameters at input and the $A^{\prime}_{2}$ parameter for the renormalon in
the evolution. The grids will be discussed in more detail in the Appendix.

##  5 High Mass Drell-Yan

###  5.1 QED and Photon PDF sensitivity in High Mass Drell-Yan

In order to explore the phenomenological implications of our photon PDF set,
we calculate the effects on the double differential cross section for lepton
pair (Drell-Yan) production at the LHC. This process is of particular
interest, since the effects of QED, especially in the partons, is expected to
be of non-negligible significance, particularly due the inclusion of
$x\gamma(x,Q^{2})$ as a contribution to the cross section. Below, we will
consider the impact of both including QED effects in the evolution of the PDFs
as well as the addition of photon-initiated (PI) contributions, as shown in
Fig. 32, where the photon PDF enters as a direct input for the colliding
partons.

Figure 32: Leading order Drell-Yan production (left), with diagrams (centre,
right) indicating $\mathcal{O}(\alpha)$ photon-initiated (PI) contributions to
the total cross section.

###  5.2 Comparison with ATLAS Drell-Yan data

In order to gauge the magnitude (and phenomenological significance) of these
effects we compare to data provided by the ATLAS collaboration [26] for high
mass (116 GeV $<m_{ll}<$ 1500 GeV) Drell-Yan lepton pair production. The focus
on production at high mass is chosen in order to reduce the effects of the $Z$
production peak, $Q\sim M_{Z}=91$ GeV, since the relative contribution of the
PI processes are greater in the regions dominated by the $\gamma$ channel.
Therefore, the effects of PI contributions are anticipated to be more readily
observable at low, $m_{ll}\ll M_{Z}^{2}$, or high, $m_{ll}\gg M_{Z}^{2}$,
lepton pair invariant masses.

ATLAS provides double differential cross section measurements in 5 bins of
lepton pair invariant mass, $m_{ll}$ and 12 or 6 pseudo-rapidity $\eta$ bins,
depending on the mass region. Fig. 33 shows as a comparison for a range of
cases: (a) a standard QCD fit partons at NNLO as outlined in Section 2.1, (b)
with QED modified partons to provide cross section calculations at NNLO in QCD
and (c) with QED modified partons and additional contributions to the cross
section from $\mathcal{O}(\alpha)$ photon initiated processes as shown in Fig.
32.

To calculate cross sections, we use grids provided by the xFitter
collaboration [13], at NLO in QCD (generated with MadGraph5$\\_$aMC@NLO [96],
aMCfast [97] and FEWZ [98]), and including PI processes at LO in QED. NNLO QCD
corrections are included via $K$-factors. Such grids were developed and used
in [13] with the aim of determining $x\gamma(x,Q^{2})$ from the same ATLAS
data. These are then interfaced with a modified version of APPLgrid that we
have adapted to include $\gamma\gamma$ processes for the final calculation.

In the following analysis it is emphasised that the contributions of PI
processes implemented in the comparison to data will be most sensitive to
$x\gamma^{(inel)}(x,Q^{2})$, due to the prevalence of this contribution in
comparison to $x\gamma^{(el)}(x,Q^{2})$ at higher scales (as was seen in the
lower part of Fig. 5 in Section 2.4).

First, it is observed that the addition of QED in the process of DGLAP leads
to a tendency to decrease the dominantly $q\bar{q}$ contribution to the cross
section, increasingly so at higher rapidity. This is expected, as from Fig. 15
one observes that the quarks experience a reduction at high $x$ of $\sim 1\%$
due to $q\rightarrow q+\gamma$ type splittings. Second, the inclusion of PI
contributions to the cross section is seen, as expected, to lead to an
increase in the cross section relative to the QED corrected partons across all
bins, as the inclusion of $x\gamma(x,Q^{2})$ opens up a new channel for lepton
pair production, unaccounted for in pure QCD calculations. Since the magnitude
of the photon PDF is seen to become larger at low $x$, particularly at high
scales ($Q^{2}=10^{4}\sim 10^{8}$ GeV2) and
$\eta\simeq\frac{1}{2}\ln{(x_{1}/x_{2})}$ where $1$ and $2$ denote the
incoming photons, the predominance of the photon at low $x$ manifests as an
enhanced cross section contribution in the lower and intermediate $\eta$ bins,
an effect seen to hold across all mass bins. At high rapidities the smallness
of the large-$x$ photon makes this photon contribution smaller than the
decrease due to the quark suppression noted above.

0.920.940.960.9811.021.041.061.081.100.511.522.5Centre of Rapidity Bin(ATLAS 8
TeV 2016), 116 GeV $<M_{ll}<$ 150 GeVQCD DGLAP PartonsQCD+QED DGLAP
PartonsQCD+QED DGLAP Partons with PIData
0.850.90.9511.051.11.1500.511.522.5Centre of Rapidity Bin(ATLAS 8 TeV 2016),
150 GeV $<M_{ll}<$ 200 GeVQCD DGLAP PartonsQCD+QED DGLAP PartonsQCD+QED DGLAP
Partons with PIData

0.90.9511.051.11.151.200.511.522.5Centre of Rapidity Bin(ATLAS 8 TeV 2016),
200 GeV $<M_{ll}<$ 300 GeVQCD DGLAP PartonsQCD+QED DGLAP PartonsQCD+QED DGLAP
Partons with PIData 0.90.9511.051.11.151.21.251.31.351.400.511.522.5Centre of
Rapidity BinData (ATLAS 8 TeV 2016), 300 GeV $<M_{ll}<$ 500 GeVQCD DGLAP
PartonsQCD+QED DGLAP PartonsQCD+QED DGLAP Partons with PIData

0.60.811.21.41.600.511.522.5Centre of Rapidity Bin(ATLAS 8 TeV 2016), 500 GeV
$<M_{ll}<$ 1500 GeVQCD DGLAP PartonsQCD+QED DGLAP PartonsQCD+QED DGLAP Partons
with PIData

Figure 33: The theory prediction/Data Ratio for ATLAS 8 TeV Drell-Yan data as
a function of rapidity in different mass bins. Shown are the predictions using
QCD only, QED included in PDF evolution but photon-initiated processes not
included, and full QCD plus QED including photon-initiated processes.

At high $\eta$, however, the change due to QED effects in the evolution is
seen to be of comparable in magnitude to that of PI contributions. In
particular, we wish to highlight that for precision calculations of
electroweak effects, one requires that all the partons be consistently treated
(i.e. to contain all QED splittings for the quarks and gluons in an
interdependent and coupled fashion) with QED in the evolution, as well as
including the photon for a consistent treatment. This is especially noteworthy
since the general trend of the partons after refitting with QED has an
opposing effect on the cross section compared to that of PI contributions (due
to a reduction of the total quark singlet), and as such, neglecting them can
in principle lead to an over-estimation of the cross section where PI
contributions are simply added on top of the standard QCD result, without the
compensating effect in the other partons.

In fact, at high $x,\eta$, where PI contributions are less relatively
important as $x\gamma(x,Q^{2})$ rapidly diminishes, the effect of refitting
the partons with QED is such that even the inclusion of PI contributions after
accounting for QED in the evolution leads to a cross section less than that of
the standard NNLO QCD prediction. In other words, the reduction of the total
quark singlet content has a greater impact than the additional cross section
contributions that are available from PI processes.

###  5.3 Including the ATLAS data in the global fit

In the aforementioned analysis, the cross section calculations are performed
using a set of PDFs which has not included the Drell-Yan data from ATLAS
itself in the global fit for the determination of parton parameters. In the
remainder of this section, we discuss the effects of including these data in
the fit itself and the subsequent effect on the recalculation of the cross
section. In Fig. 34 we present the ratio of the cross section calculation from
the QED corrected partons, including the contributions of PI processes, both
before and after refitting to the data with these effects. We can see that
there is no substantial improvement in data description after refitting.

Of note however, is the fact that the PDF contributions to the uncertainties
of the predicted cross sections (the sole contribution to the uncertainty
bands in Fig. 34) are incrementally reduced when refitting with the effects of
QED included. This is best observed in the bins for high $\eta$, especially in
the lower mass bins. In particular, we note that this incremental reduction is
seen when refitting with the effects of QED in the evolution and with the
inclusion of PI effects, but not when refit with purely with NNLO QCD parton
evolution (and with no PI contributions). This indicates a weak preference to
the effects of QED in the partons themselves and more accurate data may yet
provide a better indication of how sensitive the comparison to the theory is
with and without the effects outlined in this paper.

0.920.940.960.9811.021.041.061.081.100.511.522.5(ATLAS 8 TeV 2016), 116 GeV
$<M_{ll}<$ 150 GeVBefore RefitAfter RefitData
0.850.90.9511.051.11.1500.511.522.5(ATLAS 8 TeV 2016), 150 GeV $<M_{ll}<$ 200
GeVBefore RefitAfter RefitData

0.90.9511.051.11.151.200.511.522.5(ATLAS 8 TeV 2016), 200 GeV $<M_{ll}<$ 300
GeVBefore RefitAfter RefitData 0.911.11.21.31.41.500.511.522.5(ATLAS 8 TeV
2016), 300 GeV $<M_{ll}<$ 500 GeVBefore RefitAfter RefitData

0.60.811.21.41.61.8200.511.522.5(ATLAS 8 TeV 2016), 500 GeV $<M_{ll}<$ 1500
GeVBefore RefitAfter RefitData

Figure 34: The theory prediction/Data Ratio for ATLAS 8 TeV Drell-Yan data as
a function of rapidity in different mass bins. Shown are the predictions using
a fit in which the Drell-Yan data are not included (before refit) and once the
Drell-Yan data are included (after refit).

##  6 Conclusions

In this paper, we have presented the updated MMHT partons, modified to include
the effects of QED in their evolution. Our resultant photon PDF,
$x\gamma(x,Q^{2})$, based on a similar methodology for the input to that of
LUXqed is seen to closely resemble others in the literature, despite several
modifications made to take into account our lower starting scale for the
evolution and the fact that we use our own PDFs.

We have also outlined the procedure developed to provide an approximate QED
corrected DGLAP evolution for the PDFs of the neutron, leading to a neutron
photon PDF and isospin violating valence quark PDFs, which may hold
significance for the future development of neutron PDFs. The photon PDF of the
neutron is seen to be of a similar magnitude to that of the proton at higher
$Q^{2}$. We provide the PDFs in grids which contain the central sets and
uncertainties. PDF sets are provided as grids in the LHAPDF6 format. Details
are contained in the Appendix.

Finally, although the fit quality remains broadly unchanged after refitting
with these effects, we have observed that for the process of high-mass Drell-
Yan production, the effects of both photon initiated processes, as well as
changes in the quark and antiquark PDFs due to the effects of evolution, may
become significant with the advent of precision measurements in this kinematic
region and that the effects of QED in the evolution may be as significant as
that of the photon, highlighting a need for a fully consistent set of QED
corrected partons.

## Acknowledgements

We would like to thank Valerio Bertone and Francesco Giuli. LHL thanks the
Science and Technology Facilities Council (STFC) for support via grant awards
ST/P004547/1. RST thanks the Science and Technology Facilities Council (STFC)
for support via grant award ST/P000274/1. RN thanks the Elizabeth Spreadbury
Fund.

## Appendix – PDF Grids

As noted earlier, the set of QED corrected partons, MMHT2015qed, developed in
this paper will be released in the LHAPDF6 format for public use. The exact
nature in which the grids are provided is clarified here along with the
numbering of the grids and their associated uncertainties.

The LHAPDF6 format requires that in each file for a given grid, each column,
which represents a given PDF distribution, be labelled with an associated
number from the Monte Carlo Particle Numbering Scheme as described in [88],
where every flavour of particle is associated with an integer. This represents
an obstacle for the photon distributions as represented in this paper, since
only one such number is allocated for the $\gamma$, $22$, while we wish to
distinguish between the total, the elastic and the inelastic components.

To provide users with the ability to call upon $\gamma$, $\gamma^{(el)}$,
$\gamma^{(inel)}$, as needed, we provide three separate PDF sets for each use
case. Each set contains the full 62 eigenvector uncertainties as well as the
central values described in Section 4.3,. The ‘MMHT2015qed_nnlo_total’ set
provides the full $\gamma=\gamma^{(el)}+\gamma^{(inel)}$ distribution in the
column reserved for the photon (22). The ‘MMHT2015qed_nnlo_inelastic’ set
provides the $\gamma^{(inel)}(x,Q^{2})$ distribution while the
‘MMHT2015qed_nnlo_elastic’ set provides the $\gamma^{(el)}(x,Q^{2})$
distribution. Users should therefore distinguish by name the appropriate
LHAPDF6 variables in code for each distinct photon component as needed,
calling each from the sets as labelled above.

Each grid is a file labelled as, at NLO,
‘mmht2015qed_nlo_{type}_00$\\{x\\}$.dat’ or, at NNLO,
‘mmht2015qed_nnlo_{type}_00$\\{x\\}$.dat’, where {type} is a label denoting
which photon contribution is included in the set and $\\{x\\}$ represents
numbers in the range $\\{01,02,...,62\\}$. The particular uncertainties (as
described above) associated with the numbers denoting each set are detailed in
the Table 3.

File number index $\\{x\\}$ |  Corresponding Uncertainty
---|---
01-50 |  The standard PDF uncertainties associated with the $q+\bar{q}$, $q-\bar{q}$ and $g$ distributions for all flavours
51-52 |  The uncertainty contributions from $A^{\prime}_{2}$ (51: -0.4, 52:-0.2)
53-54 |  The uncertainty contributions from the Continuum contributions (53: Upper band, 54: Lower band)
55-56 |  The uncertainty contributions from the Resonance contributions (53: Upper band, 54: Lower band)
57-58 |  The uncertainty contributions from $W_{cut}$ (57: 3 GeV2, 58: 4 GeV2)
59-60 |  The uncertainty contributions from R (59: +50%, 60: -50%)
61-62 |  The uncertainty contributions from the Elastic contributions (53: Upper band, 54: Lower band)
Table 3: A table denoting how the numbering of the grid files (produced in the
LHAPDF6 format) corresponds to the uncertainties listed in the text.

## References

  * [1] G. Altarelli and G. Parisi,  Nucl. Phys. B126, 298 (1977).
  * [2] Y. L. Dokshitzer,  Sov. Phys. JETP 46, 641 (1977),  [Zh. Eksp. Teor. Fiz.73,1216(1977)].
  * [3] V. N. Gribov and L. N. Lipatov,  Sov. J. Nucl. Phys. 15, 438 (1972),  [Yad. Fiz.15,781(1972)].
  * [4] ATLAS, G. Aad et al.,  Phys. Lett. B725, 223 (2013), 1305.4192.
  * [5] A. Bierweiler, T. Kasprzik, J. H. Kühn, and S. Uccirati,  JHEP 11, 093 (2012), 1208.3147.
  * [6] A. Denner, S. Dittmaier, S. Kallweit, and A. Muck,  JHEP 03, 075 (2012), 1112.5142.
  * [7] L. A. Harland-Lang, V. A. Khoze, and M. G. Ryskin,  Eur. Phys. J. C76, 255 (2016), 1601.03772.
  * [8] L. A. Harland-Lang, V. A. Khoze, M. G. Ryskin, and M. Tasevsky,  (2018), 1812.04886.
  * [9] L. Beresford and J. Liu,  (2018), 1811.06465.
  * [10] A. D. Martin, R. G. Roberts, W. J. Stirling, and R. S. Thorne,  Eur. Phys. J. C39, 155 (2005), hep-ph/0411040.
  * [11] C. Schmidt, J. Pumplin, D. Stump, and C. P. Yuan,  Phys. Rev. D93, 114015 (2016), 1509.02905.
  * [12] NNPDF, R. D. Ball et al.,  Nucl. Phys. B877, 290 (2013), 1308.0598.
  * [13] xFitter Developers’ Team, F. Giuli et al.,  Eur. Phys. J. C77, 400 (2017), 1701.08553.
  * [14] A. D. Martin and M. G. Ryskin,  Eur. Phys. J. C74, 3040 (2014), 1406.2118.