Ad Lucem: QED Parton Distribution Functions in the MMHT Framework

TL;DR

This paper introduces the MMHT2015qed PDF set, incorporating QED corrections into parton distribution functions, including the photon PDF, and explores their implications for high-energy physics phenomenology.

Contribution

It presents a new PDF set with QED effects, including the photon PDF, and details the methods for coupled QED and QCD evolution at various orders.

Findings

Photon PDF consistent with recent sets

QED corrections impact on other partons analyzed

Photon-initiated processes affect high mass Drell-Yan fits

Abstract

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.

Tables

Table 1. Table 1: The total χ 2 superscript 𝜒 2 \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 χ 2 superscript 𝜒 2 \chi^{2} due to the inclusion of the QED corrections.

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 2. Table 2: The χ 2 superscript 𝜒 2 \chi^{2} breakdown showing χ 2 / N p ​ t ​ s superscript 𝜒 2 subscript 𝑁 𝑝 𝑡 𝑠 \chi^{2}/N_{pts} by data set for NNLO QCD and NNLO QCD + QED PDF fits.

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$ $x{F}_3$ [63]	31 / 42	31 / 42
CHORUS $\nu N$ $x{F}_3$ [42]	19 / 28	19 / 28
CCFR $\nu N\to \mu \mu X$ [64]	77 / 86	78 / 86
NuTeV $\nu N\to \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\overline{p}$ incl. jets [66]	117 / 110	117 / 110
CDF II $p\overline{p}$ incl. jets [67]	60 / 76	60 / 76
CDF II $W$ asm. [68]	16 / 13	15 / 13
DØ II $W\to \nu e$ asym. [69]	31 / 12	30 / 12
DØ II $W\to \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\to e^{+}{e}^{-}$ [76]	22 / 9	22 / 9
LHCb $W$ asymm $p_T>20$ GeV [77]	14 / 10	13 / 10
CMS $Z\to 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\overline{t}}$* [80, 81, 82, 83, 84, 85, 86]	14 / 18	14 / 18
All data	3532 / 3276	3539/3276 (+7)

Table 3. 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.

File number index $\left\{x\right\}$	Corresponding Uncertainty
01-50	The standard PDF uncertainties associated with the $q+\overline{q}$, $q-\overline{q}$ and $g$ distributions for all flavours
51-52	The uncertainty contributions from $A_2^{\prime }$ (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)

Equations

$$\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}$$

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

$$\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}$$

$$\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}$$

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

$$\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}$$

$$\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}$$

$$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})\;,$$

$$Q^2\ge \frac{x^2{m}_p^2}{1-x},$$

$$x\le \frac{-Q^2+Q\sqrt{Q^2+4m_p^2}}{2m_p^2}\equiv x_{\mathrm{cut}},$$

$$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{)},$$

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

$$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}}}\;,$$

$${\int }_0^{1}x(\Sigma (x,Q_0^2)+g(x,Q_0^2)+\gamma (x,Q_0^2))=1,$$

$${\Sigma }_L=u+\bar{u}+d+\bar{d}+s+\bar{s},$$

$$c+\bar{c},b+\bar{b},$$

$$g$$

$$u_V+d_V=u-\bar{u}+d-\bar{d},$$

$$u+\bar{u}-\frac{{\Sigma }_L}{3},-(s+\bar{s})+\frac{{\Sigma }_L}{3},$$

$$\frac{d_V-u_V}{2},(s-\bar{s})-\frac{u_V+d_V}{2},(c-\bar{c})-\frac{u_V+d_V}{2},(b-\bar{b})-\frac{u_V+d_V}{2},$$

$$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)}+...$$

$$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)}+...$$

$$P_{qq}^{(0,1)}$$

$$P_{\gamma q}^{(0,1)}$$

$$\sum _i{e}_i^2=N_C\sum _q^{n_F}{e}_q^2+\sum _l^{n_L}{e}_l^2.$$

$$q_i^{\pm }=q_i\pm \overline{q_i},g,{\gamma }^{(el)},{\gamma }^{(inel)}.$$

$$\frac{\partial q_i^{NS}}{\partial t}=P_{q_i}^{-}\otimes q_i^{NS},$$

$$\frac{\partial q_i^{-}}{\partial t}=P_{q_i}^{-}\otimes q_i^{-}+\sum _{j=1}^{n_F}\Delta P^S\otimes q_j^{-},$$

$$P_{\gamma ,q}^{(0,1)}(z)\to P_{\gamma ,q}^{(0,1)}(z)+\frac{2x^2{m}_p^2}{z{Q}^2}.$$

$$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{)},$$

$${\int }_0^{1}dx{F}_2^{HT}(x,Q^2)=0,$$

Full text

1. 1 Introduction

2. 2 Including QED Effects in the MMHT Framework

3. 2.1 Baseline QCD Fit

4. 2.2 Input Photon Distribution

5. 2.3 Modifications to DGLAP evolution

6. 2.3.1 PDF Basis

7. 2.3.2 Target Mass and Higher Twist Corrections

8. 2.4 Separation of Elastic and Inelastic Components

9. 2.5 Momentum Conservation

10. 3 QED Neutron PDFs

11. 3.1 Modified DGLAP Evolution

12. 3.2 The Photon PDF of the Neutron

13. 4 Results

14. 4.1 Changes to PDFs due to QED corrections

15. 4.2 Results of Global fits with QED corrections

16. 4.2.1 The quality of the global fits

17. 4.2.2 The photon PDF of the proton

18. 4.2.3 QED Corrected Structure Functions

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

20. 4.2.5 Photon-photon luminosity

21. 4.3 Uncertainties on the photon PDF

22. 5 High Mass Drell-Yan

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

24. 5.2 Comparison with ATLAS Drell-Yan data

25. 5.3 Including the ATLAS data in the global fit

26. 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:

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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$,

[TABLE]

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:

[TABLE]

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.

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:

[TABLE]

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

[TABLE]

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)$):

[TABLE]

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:

[TABLE]

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:

[TABLE]

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$):

[TABLE]

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:

[TABLE]

[TABLE]

[TABLE]

and the following non-singlet combinations:

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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:

[TABLE]

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):

[TABLE]

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

[TABLE]

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,

[TABLE]

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.

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.

2.4 Separation of Elastic and Inelastic Components

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):

[TABLE]

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:

[TABLE]

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}$:

[TABLE]

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.

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.

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

[TABLE]

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 [math] 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.

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:

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

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):

[TABLE]

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

[TABLE]

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:

[TABLE]

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

[TABLE]

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:

[TABLE]

[TABLE]

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.

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 [math], 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.

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]:

[TABLE]

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})$):

[TABLE]

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

[TABLE]

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

[TABLE]

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$.

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

[TABLE]

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:

[TABLE]

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

[TABLE]

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.

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

[TABLE]

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:

[TABLE]

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.

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.

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:

[TABLE]

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).

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.

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.

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.

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.

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

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.

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.

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.

4.2.3 QED Corrected Structure Functions

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

[TABLE]

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.

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:

[TABLE]

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).

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.

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.

• •

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:

[TABLE]

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.

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.

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.

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.