The Photon Content of the Proton

The photon PDF of the proton is needed for precision comparisons of LHC cross sections with theoretical predictions. In a recent paper, we showed how the photon PDF could be determined in terms of the electromagnetic proton structure functions $F_2$ and $F_L$ measured in electron-proton scattering experiments, and gave an explicit formula for the PDF including all terms up to next-to-leading order. In this paper we give details of the derivation. We obtain the photon PDF using the factorisation theorem and applying it to suitable BSM hard scattering processes. We also obtain the same PDF in a process-independent manner using the usual definition of PDFs in terms of light-cone Fourier transforms of products of operators. We show how our method gives an exact representation for the photon PDF in terms of $F_2$ and $F_L$, valid to all orders in QED and QCD, and including all non-perturbative corrections. This representation is then used to give an explicit formula for the photon PDF to one order higher than our previous result. We also generalise our results to obtain formul\ae\ for the polarised photon PDF, as well as the photon TMDPDF. Using our formula, we derive the $P_{\gamma i}$ subset of DGLAP splitting functions to order $\alpha \alpha_s$ and $\alpha^2$, which agree with known results. We give a detailed explanation of the approach that we follow to determine a photon PDF and its uncertainty within the above framework.


Contents
modelling. The main idea in our derivation is that a photon initiated process can be computed two different ways: in terms of a photon PDF, using the factorisation theorem, or in terms of the deep inelastic scattering (DIS) hadronic tensor. Equating the two expressions leads to our result. 2 All information needed to extract f γ is thus contained in ep scattering data. This point of view is implicit also in Refs. [24,25].
In this paper, we give additional details on the derivation of the photon PDF formula. We report the explicit calculations performed for the two different hard probes that we consider: heavy lepton production by a flavour violating magnetic moment interaction, and scalar production via γγ collisions. We also give an independent derivation of the formula by defining the photon PDF as the light-cone Fourier transform of the two-point function of electromagnetic field-strength tensors. Our derivation gives a representation for the photon PDF which is exact, i.e. valid to all orders in the strong and electromagnetic interactions. We use this exact representation to obtain expressions for the photon PDF including all terms of order αα s and α 2 , one order higher in α s (µ) or α(µ) than our original result [22].
In the present work we have slightly revised the photon PDF calculation and error estimate compared with Ref. [22], also taking into account the new results presented here (see Sec. 9.2). For this reason, we refer to the photon PDF computed using the method in this paper as LUXqed17, and that computed using the previous procedure as LUXqed. The ratio of the two is shown in Fig. 14. As one can see the ensuing changes are very minor, and confirm the error estimate given in Ref. [22].
Our formula for the photon PDF can be differentiated w.r.t. µ to give the P γi subset of the DGLAP splitting functions. The result in Ref. [22] was used to obtain the two-loop order αα s splitting kernels, which agree with a recent computation [26]. Here, we present this calculation in detail, including also the computation of the splitting kernels of order α 2 , which agree with Ref. [27].
We also give some natural extensions of our results. We obtain the photon transverse momentum dependent PDF (TMDPDF), and by considering spin-dependent scattering, we obtain the polarised photon PDF f ∆γ (x, µ 2 ) in terms of the polarised structure functions g 1,2 (x, Q 2 ).
The outline of the paper is as follows. Section 2 introduces our notation and reviews some known results on QED corrections to the hadronic tensors, which are important for obtaining our all-orders formula. Sec. 3 gives the derivation of the photon PDF using heavy-lepton production via a magnetic moment interaction. Section 4 derives the order αα s and α 2 splitting functions from our photon PDF formula. The extension of our results to polarised PDFs is given in Sec. 5. Section 6 gives the derivation of the photon PDF for the polarised and unpolarised cases in terms of PDF operators. The result in this section is exact, and is also used to obtain the photon TMDPDF. Section 7 uses the exact result in Sec. 6 to obtain the PDF to order α, αα s and α 2 s relative to the leading order term. In Sec. 8 we discuss the experimental data inputs used for our numerical evaluation of the photon PDF. In Sec. 9 we give our procedure for estimating errors from missing higher-order corrections. In Sec. 10 we give all details about the practical implementation of the photon PDF formula and how the photon is matched to other partons. Some results on the photon PDF are given in Sec. 11. We present our conclusions in Sec. 12. Useful technical results are given in the Appendices. Appendix A gives technical details on QED corrections, App. B discusses the kinematics for our exact representation, App. C derives the photon PDF using scalar production in γγ collisions, App. D presents the parton model NLO calculation for the parton level process ql → qL, App. E and F discuss the low-energy behaviour for F 2/L and R, App. G summarises the DIS coefficient and splitting functions and App. H discusses issues in PDF4LHC15 at low scales.

Definitions and notation
To aid the reader, it is useful to introduce some elements of our notation. We will use dimensional regularisation in D = 4 − 2 dimensions. According to the MS prescription, the dimensional regularisation scale µ is replaced by µ → (Sµ), with . (2.1) Subtracting only 1/ poles then gives the renormalised result in the MS scheme. The notation i ∈ {q, l, g, γ}, means a sum over quarks and antiquarks, leptons and antileptons, photons, and gluons. Similarly for i ∈ {q}, i ∈ {q, g}, i ∈ {q, l, γ}, etc. Colour multiplicities, when needed, will be denoted by n i and written explicitly, n i = 1 for leptons and n i = 3 for quarks.
The perturbative expansion of an object O, such as a partonic cross section σ or a splitting kernel P will be written as One of the results of this paper is a representation for the photon PDF which is exact, including all radiative and non-perturbative QED and QCD corrections. This requires a careful treatment of electromagnetic radiative corrections. 3 A brief summary of known results on electromagnetic radiative corrections and the renormalisation of the electromagnetic current is given in App. A. We will use the results in the appendix to motivate the definitions in this section.
The proton hadronic tensor, which enters the cross section formula for deep-inelastic scattering, is defined as W µν (p, q, s) = 1 4π d 4 z e iq·z p, s| [j µ (z), j ν (0)] |p, s , (2.4) where p is the incoming proton momentum, q is the photon momentum (transferred to the proton) and s denotes the proton spin. W µν is the discontinuity of the time-ordered product T µν (p, q, s) = i d 4 z e iq·z p, s|T {j µ (z), j ν (0)} |p, s . (2.5) The discontinuity for q 0 > 0 gives the proton hadronic tensor, corresponding to the j µ (z)j ν (0) term in Eq. (2.4). The discontinuity for q 0 < 0 is related by crossing to the anti-proton hadronic tensor and corresponds to the second term of the commutator, which will not be needed here.
The conventional decomposition of W µν into structure functions, after requiring Lorentz invariance, time reversal, parity, and current conservation is ig 2 (p · q) 2 µνλσ q λ (p · q s σ − s · q p σ ) . (2.6) The proton spin four-vector s is normalised so that p · s = 0, s · s = −m 2 p , where m p is the proton mass. The structure functions F 1 , F 2 , g 1 and g 2 are functions of and Q 2 . We also introduce the longitudinal structure function F L (x bj , Q 2 ) ≡ 1 + 4x 2 bj m 2 p Q 2 F 2 (x bj , Q 2 ) − 2x bj F 1 (x bj , Q 2 ), (2.8) and write our results using F 2 and F L instead of F 2 and F 1 .
The usual textbook analysis of deep inelastic scattering is at lowest order in the electromagnetic coupling. In this paper, we are also interested in higher-order electromagnetic corrections, so we should be careful about the definition of W µν . We define the hadronic tensor in Eq. (2.4) to be given by one-photon-irreducible graphs in the current channels, and including all electromagnetic corrections to the hadronic matrix element. With this definition, W µν is µ-independent (see App. A), so the structure functions in the decomposition Eq. (2.6) are also µ-independent, and depend only on x bj and Q 2 . 4 The structure is an O(α) contribution to the quark coefficient, but is one-photon-reducible, and with our prescription does not contribute to the coefficient functions C 2/L . functions and W µν defined in terms of one-photon irreducible graphs will be denoted by the label 1γI in App. A. In the main text of this paper, we drop this label for simplicity since we are always referring to the one-photon irreducible hadronic tensor.
The structure functions F 2/L can be computed using the QCD improved parton model formula if Q 2 is large, F 2/L (x, Q 2 ) = i∈{q,l,g,γ} zy,x C 2/L,i (z, Q 2 , µ) f i (y, µ 2 ) , (2.9) where we have introduced the notation z 1 ···zn,x = 1 0 dz 1 · · · dz n δ(z 1 · · · z n − x) , (2.10) as the convolution of short distance coefficient functions C 2/L,i (x, Q 2 , µ) and PDFs f i (x, µ 2 ), and the notation ∈ {q, l, g, γ} is defined at the beginning of Sec. 2. A few sample graphs for C 2/L,i (x, Q 2 , µ) are shown in Fig. 1 (a-d). The graph in Fig. 1(e) with our prescription does not contribute to the coefficient functions C 2/L,i , as it is one-photon reducible.

Photon distribution via a photon-probing process
In this section we will assume that we have at our disposal a BSM (beyond standard model) photon probe that couples to SM particles only through photon exchange. We neglect weak p p l(k) l(k) q q X L(k′) Figure 2. The process l + p → L + X to leading order in 1/Λ, but including QED corrections to all orders. The blob in the photon propagators represent the full vacuum polarisation correction. The lower hadronic blob also includes QED corrections.
interactions, and thus do not worry about making this interaction SU (2) × U(1) invariant. The BSM photon probe combined with factorisation will allow us to determine the photon PDF, which must be independent of the probe we have chosen. We verify this by performing the same computation with a different probe in App. C. The BSM process we consider here involves two neutral spin-1/2 leptons, an incoming lepton l and an outgoing lepton L with masses zero and M m p , and momenta k and k = k − q, respectively, with a transition (i.e. flavour violating) magnetic moment interaction L = c e(µ) Λ L σ µν l F µν + h.c. . where σ µν = i 2 [γ µ , γ ν ]. We formally assume the limit Λ → ∞, i.e. we work to lowest order in 1/Λ. The QED Ward identity Eq. (A.3) implies that c is independent of µ to all orders in α. With the interaction Eq. (3.1), the interaction vertex is (not including the overall i) c e(µ) Λ u(k − q) / , / q u(k) .

(3.2)
We define the spin-averaged leptonic tensor for the l(k) → L(k ) transition as L µν (q, k) = 1 2 Evaluating the traces gives which satisfies the conditions q µ L µν = q ν L µν = 0. We now show how to express a photon initiated process in terms of hadron structure functions. The cross section for the process l + p → L + X, depicted in Fig. 2, including all QED corrections, is This expression is exact. The vertices on the lepton line in the figure give the leptonic tensor L µν of Eq. (3.4), and on the lower line give the hadronic tensor 4πW µν , defined in Eq. (2.4). There are vacuum polarisation corrections, represented as dashed blobs, to each photon propagator connecting the lepton and hadron sides of the graph, which give a factor of 1/[1 − Π(q 2 , µ)] 2 . In addition, there can be arbitrary QED corrections included in the hadron part of the diagram (except for the one-particle reducible ones) which are included in our definition of W µν . QED corrections on the leptonic line are suppressed by powers of Λ, and disappear in our limit, as do corrections associated with the exchange of more than one photon. Eq. (2.3) has been used to convert the powers of e 4 (µ) together with the factor 1/[1 − Π(q 2 , µ)] 2 into e 4 ph (q 2 ). The δ and θ functions for the second line in Eq. (3.5) ensure that L is on-shell, and has positive energy. The θ functions on the third line of the same equation ensure that X has positive energy, and m 2 X ≥ m 2 p , since the proton is the lightest baryon. The phase space is given by . (3.6) Here z and x are defined as where x bj = Q 2 /(2p · q) is the usual Bjorken variable. The limits on the Q 2 integration are given by Expanding the limits in m 2 p /M 2 gives Details of the phase space and kinematics computation are given in App. B.
A straightforward calculation, combining the definition of the hadronic and leptonic tensors Eqs. (2.6) and (3.4), yields that, together with Eq. (3.5) and the phase space Eq. (3.6), gives The final cross section depends on x, which is fixed by the external kinematic variables M (the mass of L), and 2p · k = s − m 2 p , where √ s is the centre-of-mass energy. We stress again that Eq. (3.11) is exact as long as Λ suppressed terms are neglected, α ph is exact, and F 2/L are defined from the full one-particle irreducible hadronic tensor including QED corrections.
We now turn to the calculation of the cross section using the QCD improved parton model formula σ lp (p) = j∈{q,l,g,γ} dy σ lj (yp) f j (y, µ 2 ) , (3.13) where f j are the PDFs of the proton. The PDFs are in the MS scheme if the partonic hard scattering cross sections σ are computed in the MS scheme. Power corrections of the form m 2 p /M 2 or m 2 p /(2p · k) are neglected. Expanding this formula in the coupling constants yields The diagrams corresponding to σ (0,0) lγ and σ (0,1) lj are shown in Fig. 3. The first diagram is the lowest order term, involving directly the proton photon-density. The second diagram is the next-to-leading correction of relative order α, involving the proton charged-fermions densities. The contribution corresponding to the rightmost graph yields zero in the MS scheme since we consider massless fermions, since at this order in collinear factorisation the photon is implicitly on shell and so the vacuum polarisation is evaluated at zero virtuality. Thus the index j in the sum is restricted in practice to the charged fermions only. . Leading and next-to-leading graphs for the process l + γ → L in the QCD improved parton model. At this point a comment is in order. We can systematically compute the cross section assuming that α and α s are of the same size, and that the parton densities themselves are formally all of the same order. We dub this counting of the order "democratic", and adopt it here in what follows, since it is more transparent. In the democratic order-counting, the index i appearing in Eq. (3.14) should also run over leptons. Furthermore, neglected terms are of second order in both α and α s , i.e. of order α 2 and αα s (the α 2 s term being absent), relative to the Born term.
For phenomenological applications, however, we will take into account the fact that α is smaller than α s , using as a guideline the relation α ≈ α 2 s . We dub this counting "phenomenological". According to it, the photon density of the proton is of order αL with respect to a quark density, L being a log of µ 2 over some typical hadronic scale. We can assume L ≈ 1/α s . In this framework the contributions corresponding to the first and second diagram in Fig. 3.14 are respectively of order α 2 L, α 2 , while the last graph is formally of order α 3 L ≈ α 2 α s (but is zero in the MS scheme). The next-to-leading correction is of relative order 1/L ∼ α s , rather than of order α (as in the democratic counting), with respect to the Born term. In the middle diagram of Fig. 3 light leptons can be excluded, since their PDF is of order L 2 α 2 , and their contribution is of order α 4 L 2 . 5 The cross section for the process σ(l + q → L + q), illustrated in the middle graph of Fig. 3, is easily computed with standard methods. Details of the calculation are given in App. D. We get where σ 0 is given in Eq. (3.12),ŝ = ys, z = M 2 /ŝ = x/y and This yields We now define a "physical" photon PDF The reason for the upper limit on the integration will become clear later. The terminology "physical" is associated with the fact that when we consider the transverse-momentum dependent photon PDF in Sec. 6.3, the latter's integral up to k 2 ⊥ ≤ µ 2 will coincide with Eq. (3.21). As we will see below, Eq. (3.21) includes the αL contribution to the photon PDF, but not the total α piece.
Combining Eq. (3.21) with Eq. (3.20), we obtain Since the remaining Q 2 integral is now dominated by large Q 2 values, we can evaluate it with the same approximations used earlier, and get Now consider the parton model formula for the same cross section, Eq. (3.18). Using the lowest order expression in α s (µ), α(µ) The l.h.s. is the desired result, the photon PDF in the MS scheme. The r.h.s. expresses it as an integral over lepton-proton scattering structure functions. The first term is given in Eq. (3.21), and we refer to the second term as the "MS-conversion term". The upper integration limit in Eq. (3.21) was chosen to be µ 2 /(1 − z) so that the logarithms cancelled between Eq. (3.23) and Eq. (3.18). Writing the first term explicitly gives Formula Eq. (3.26) is independent of the particular probe process that we have chosen. In App. C we verify this by carrying out the same derivation using as a probe the production of a scalar particle via photon-photon fusion and using PDF operators in Sec.6. A more direct method of obtaining Eq. (3.26) is illustrated in Sec. 6, and will be used in Sec. 7 to extend its accuracy to higher orders. Examining Eq. (3.26) according to our "phenomenological" counting, the dominant terms under the Q 2 integration are of order αL, and the MS-conversion term is of order α. For the contribution of terms under the Q 2 integration, we should be careful to include terms of relative order αL in both α ph and F 2 . These, together with the L coming from the logarithmic integration, yield corrections of order α 2 L 2 ≈ αα 2 s L 2 ≈ α. Thus α ph (−Q 2 ) can be replaced with α(Q 2 ), since they differ by terms of relative higher order in α without powers of L enhancement. We plan to extract F 2/L at low and moderate Q 2 directly from data, while for high Q 2 we will compute it using available PDF parametrisations. Thus the discussion of its accuracy is more delicate, and we postpone it to Sec. 8.

Connection with splitting functions
From our formula for f γ , Eq. (3.26), we can derive formulae for the splitting functions of the photon. To this end we will adopt the "democratic" counting, taking Eq. (3.26) at full leading order accuracy in α and α s , and treating the two couplings as if they were of the same order. Neglected terms are therefore of order α 2 and αα s , while α 2 s is obviously absent.
In order to simplify our notation we will use the following conventions. When we write the PDFs or the MS couplings without a scale argument, we imply that the scale is µ 2 . Furthermore we will adopt the normalisation convention and the convolution is defined in Eq. (2.10). The DIS coefficient functions have the expansion where We also define The factor of 1/2 inside the parenthesis is due to the fact that we sum over all charged fermions and antifermions. In the following we will also use the notation The RG equation for the QED coupling is We have where n i is 3 for quarks and 1 for leptons. Furthermore we define We now begin by re-writing Eq. (3.26) with a suitable change in the upper limit of integration and dropping power-suppressed terms, in order to ease the derivation of the evolution equation: To find the evolution equation, we compute where we have neglected terms suppressed by powers of m 2 p /µ 2 and, for ease of notation, we have omitted the µ-dependence in the coupling and parton densities. In the first line we must use for α ph (−µ 2 ) and F 2/L expressions that are accurate at order α s and α: where the coefficient functions are given in App. G. In the second line, the derivatives of F 2 and α are only needed to leading order. We have dF 2 (y) d log µ 2 = y vw,y i∈{q} Inserting these expressions we get We can now read off from Eq. (4.16) the splitting function coefficients, 20) L,q (y) , (4.21) . The order α kernels P (0,1) agree with the known expressions given in Ref. [29]. Notice that P (0,1) γγ has the correct sign, because our photon PDF Eq. (3.26) is proportional to 1/α(µ), not α(µ). The order α s and α coefficient functions for F 2/L are summarised in App. G.
Evaluating the convolutions for Eqs. (4.21) and (4.23), we get full agreement with the recent calculation in Ref. [26], Eqs. (35,36) and (30). Similarly, evaluating the integrals of the second order QED splitting functions, Eqs. (4.20) and (4.22), we find full agreement with formulae (3.9) and (3.21) of Ref. [27]. It is remarkable that our expression for the photon PDF gives the DGLAP evolution kernel to one higher order in the couplings than was present in the input coefficient functions. Specifically we obtained the (two-loop) order αα s and α 2 P γi splitting kernels using coefficient functions at (one-loop) order α s and α.

Spin dependent case
There is an extensive experimental and theoretical effort to understand the polarised gluon distribution in the proton, and ∆g, the gluon contribution to the proton spin [30][31][32][33]. ∆γ is the photon analogue of ∆g. Since photons and gluons both couple to quarks via a γ µ interaction, ∆γ could shed some light on ∆g.
The results of Sec. 3 can be readily generalised to the spin-dependent case, to obtain the polarised photon PDF of a proton of helicity H, the difference between the probabilities to find photons with spin parallel and anti-parallel to the proton spin in a longitudinally polarised proton target.
The derivation here of the polarised photon PDF follows that for the unpolarised PDF, using the polarisation asymmetry in the cross sections instead of the spin-averaged cross sections. We start with the same interaction Lagrangian as before, Eq. (5.1). The difference between the cross sections for right-handed and left-handed leptons to scatter off a longitudinally polarised proton with helicity H to order ασ 0 is, using Eq. (2.6) for W µν , To obtain the photon polarised parton density f ∆γ , we use the factorisation formula for the scattering of a lepton l with helicity h off a target p with helicity H. The sum on i s is over all partons and their helicities. We use the notation f ∆q = f q R /p R − f q L /p R for the polarised parton distributions. Note that f q R /p R = f q L /p L and σ l R p R = σ l L p L , etc. by parity invariance. Since f ∆γ is of order α, we need the photon hard-scattering cross section to lowest order, where σ 0 is given in Eq. (3.12). Observe that the total (spin averaged) cross section is given by 4) and this implies, together with Eq. (5.3), that σ l L γ R = 0. This result is easily understood. We can look at the cross section in the centre-of-mass frame where we collide a right moving left-handed lepton onto a left moving right-handed photon, forming a heavy lepton at rest. By angular momentum conservation, the heavy lepton must have J z = 3/2 along the collision axis, which is not possible for a spin-1/2 particle. This simple argument confirms the correctness of the sign in Eq. (5.3).
The polarisation asymmetry to order ασ 0 is where the 1/ term is an infrared divergence, and we have used the 't Hooft-Veltman scheme for γ 5 . Combining with Eqs. (5.2) and (5.3) and performing the MS subtraction we get We now have all the ingredients necessary to determine the polarised photon PDF. We follow the derivation for the unpolarised case given in Sec. 3. Define, as before, the polarised PDF in a "physical factorisation" scheme by dropping the 1/M 2 suppressed terms in Eq. (5.1). As in Sec. 3, the integral in Eq. (5.1) can be written as the integral for the terms included in Eq. (5.7), divided into the range Q 2 min → µ 2 /(1 − z) and µ 2 /(1 − z) → Q 2 max , plus the integral over the remaining terms. The remaining integrals only get contributions from Q 2 of order M 2 , since µ ∼ M , and so contain no large logarithms. Since Q 2 is large over the integration region, we can replace g 1 (x/z, Q 2 ) by g 1 (x/z, µ 2 ) and α 2 ph (−Q 2 ) by α(µ 2 ) up to corrections of order α s (µ 2 ) and α(µ 2 ). The integral is now trivial, and gives The g 1 structure function is where C ∆q = δ(1 − z) to lowest order in α s . Comparing Eqs. (5.6, 5.8) we get for the polarised photon PDF. Note that the log terms cancel. Taking the first moment of our result, Eq. (5.10) gives ∆γ, the photon contribution to the proton spin. As in the unpolarised case, we have checked that the process γγ → S gives the same result. Furthermore, we can easily derive from these expressions a number of polarised splitting functions involving the photon. The expressions are similar to those in Eqs. (4.17-4.23), but we don't give them here.
6 Alternative derivation using PDF operators

Collinear photon PDF
We have derived the photon PDF by using factorisation applied to two different processes, l + p → L + X and γ + p → S + X in Secs. 3, 5 and App. C, and shown that we get an MS PDF that is process independent. PDFs can also be defined as the matrix element of a PDF operator, which is a light-cone Fourier transform of a two-point function [34], and does not depend on using a probe process. This has the advantage that it yields much simpler calculations. In this section, we derive the photon PDF using PDF operators. We will exploit the simplicity of this method when going to higher orders in the next section.
The operator PDF definition is manifestly process independent, and directly gives the MS PDF. The photon PDF operator can be obtained by analogy with the gluon PDF  operator in Ref. [34] by replacing the gluon field-strength tensor G µν by the photon fieldstrength F µν , and dropping the Wilson line, since F µν is gauge invariant: We have introduced a light-like vector n, and we have defined p + = p·n. We use the notation F nλ ≡ n σ F σλ , etc. in this section. The subscript c is a reminder that only the connected matrix element contributes. The operators are at the same x + and x ⊥ coordinates, and Fourier transformed along the x − direction. The µ dependence arises because the matrix element has divergences, and is renormalised in the MS scheme (so that the PDF is in the MS scheme). All other results in this section are also renormalised in the MS scheme and, as in Eq. (6.1), we do not explicitly show the counterterms. The polarised photon PDF can be obtained in the same way from the polarised gluon PDF [35,36], F αβ = 1 2 αβλσ F λσ , with 0123 = +1, and we conventionally assume that the proton has positive helicity. The -tensor is defined to live in the four physical dimensions, according to the 't Hooft-Veltman scheme.
The PDFs can be computed from Eqs. (6.1, 6.2). The leading order graphs are shown in Fig. 4. The key observation is that the lower part of the diagram, the interaction of the photon with the proton target, is precisely the definition of the hadronic tensor W µν (p, q), because the photon interacts with the proton through the electromagnetic current j µ that appears in Eq. (2.6). 6 Note that Eq. (6.1) involves the ordinary product of operators, not 6 It is this factorisation property which allows us to compute the photon PDF. There are expressions analogous to Eqs. (6.1,6.2) with Fµν replaced by Gµν which give the gluon PDF. Since the gluon is a colour adjoint, there is also a colour Wilson line between the two field strength tensors in the operator product. The gluon analogue of graphs Fig. 4 cannot be written in terms of experimentally measured proton structure functions, since the gluon does not couple to the proton via gauge invariant colour singlet operators. In addition, there are graphs with gluon exchange between the Wilson line and the proton. the time-ordered product. The matrix element in Eq. (6.1) can be evaluated by comparing with the expression for W µν in Eq. (2.4), which corresponds to taking the discontinuity of the time-ordered product, T µν in Eq. (2.5), i.e. evaluating the cut graphs in Fig. 4. The diagrams give where S is defined in Eq. (2.1), and Π D (q 2 , µ 2 ) is the MS subtracted polarisation function computed in D dimensions (i.e. the D = 4 limit is not yet taken).
The two δ-functions arise from the two operator orderings in Eqs. (6.1,6.2), the next factor in each equation is the Feynman rule for the field-strength tensors, and the last factor is from the photon propagators and the proton matrix element, which, by definition is the one-photon irreducible hadron tensor W (D) µν . The two contributions are from the direct and crossed graphs. The formula Eq. (6.3) is exact, and has no QCD or QED corrections, provided that all QED corrections on the hadronic side are included in the definition of the structure functions. We stress that also the hadronic tensor appearing in Eq. (6.3) is MS subtracted, but is evaluated in D dimensions. The limit D → 4 can only be taken after the overall MS counterterm is added to Eq. (6.3).
The vacuum polarisation factors can be written in terms of α ph using Eq. (2.3). The two W (D) µν terms can be combined by letting q → −q. The hadronic tensor W (D) µν is non-zero for 0 ≤ −q 2 /(2p · q) ≤ 1, so only the q + < 0 piece contributes, and we find µν n µ n ν and nµqν = αµβν n α q β . The physical coupling constant e ph,D , in analogy to Eq. (2.3), is defined as The D dimensional tensor vanishes if its indices are outside 4 dimensions. The contrac- where x bj = Q 2 /(2p · q) and the structure functions (F 2/1/L,D and g 1/2,D ) are evaluated at x bj , Q 2 in D dimensions with standard MS subtraction (i.e., as in the case of the physical electromagnetic coupling, the limit D = 4 is not yet taken). The second expression depends on Q 2 4 involving only the component of Q in D = 4 dimensions. We have where Q 2 −2 is the fractional dimension part of Q. Let q + = −xp + and z = x/x bj then, for a right-handed proton moving along the z axis, and Eq. (6.6) becomes where the structure functions are evaluated at x/z, Q 2 . The q integral can be performed by switching to light-cone coordinates, We also have which is used to replace the integration variable q − by z. The q + integral is evaluated using the δ-function, leaving the integrals over z and Q 2 . Q 2 −2 is in the ⊥ direction, and all ⊥ directions are equivalent, so we can make the replacement in evaluating the integrals. The kinematically allowed region in the z, Q 2 plane is The first inequality follows because W µν vanishes for w ≥ 1, and the second from Q 2 ⊥ ≥ 0. The resulting integrals are with α ph,D (q 2 ) = e 2 ph,D (q 2 )/(4π). This expression for the PDF is exact, and includes QED radiative corrections if the structure functions are the one-photon-irreducible ones.
We split the The first part is finite, and gives (by definition of f PF γ,∆γ ) which is the result obtained previously for what was referred to as the "physical factorisation" PDF, Eqs. (3.21) and (5.7). We automatically get the correct integrand for the PDF, without dropping any terms, as had to be done for the BSM probes in Sec. 3. The "MS -conversion" term is the remaining integral over µ 2 /(1 − z) → ∞ plus the counterterms. Since Q 2 µ 2 , we can set m 2 p → 0, up to neglected power corrections, and use Eq.
Eq. (6.17) is still exact. The integral is over large values of Q, so the q dependence of α ph,D , F 2/L,D and g 1/2,D can be computed in perturbation theory. The integral does not involve any large logarithms, since it only involves the scale µ. 7 Introducing the dimensionless variable s, In fact, the integral can instead generate a further 1/ pole, but not logarithms.
the integrals become The factor of (1 − z) D/2−2 has cancelled, which is why µ 2 /(1 − z) was chosen as the intermediate Q 2 value to split the integrals.
For the remainder of this section, we restrict to lowest order in α and α s . At this order, the structure functions do not depend on Q 2 , and can be evaluated at µ 2 without incurring any large logarithms, and α ph,D (q 2 ) ≈ α(µ 2 )(µS) 2 from Eq. (6.5). Since F L is order α s , it can be dropped. Evaluating the s integral, and expanding in gives The 1/ term is absorbed by the MS counterterm, and one obtains in the MS scheme The PDFs in the MS scheme are thus The results agree with our previous results using the lγ → L and γγ → S processes.
The derivation in this section shows the PDFs are process independent, since it does not use any physical process. It also shows how to systematically extend the result to higher orders. The PF term Eq. (6.16) is exact, and the only corrections are to the MS-conversion terms Eq. (6.19), which does not contain any large logarithms. The result Eq. (6.22) can be extended to higher orders using the perturbation expansion for α ph,D and for the structure functions in D dimensions. The result to one higher order is given in Sec. 7. Unlike the method of the earlier sections, we do not have to compute any hard scattering cross sections to higher order when using PDF operators.

Derivation using an abstract probe
Using the operator definition of the photon PDF avoids the complication of having to compute a specific probe process. On the other hand, the freedom of choosing an appropriate probe process can be exploited to achieve a similar simplicity. In particular we may consider the following photon probe tensor to be contracted with the hadronic tensor in Eq. (3.5). This probe does not refer to any specific BSM process, and indeed it is difficult to imagine a physical process that would give rise to this tensor. However, this does not matter as long as the tensor has the appropriate transversality properties, which implies that the "cross section" can be computed both in terms of the proton electromagnetic structure functions and in terms of the parton model formula.
The calculation in terms of structure functions yields The 1/(2p 0 ) is from the incident proton flux, since the proton state is normalised to 2p 0 . The parton model calculation yields instead In the σ PM term, the partonic cross section was evaluated by setting q = −py in Eq. (6.23), contracting it with a −g µν , dividing by the number of photon polarisations in D dimensions, and dividing by the incident parton flux 2yp 0 . The σ (HO) PM term must be computed in D dimensions, for incoming massless particles. This leads to the presence of collinear singularities, that are removed according to the usual factorisation procedure. In order to write the HO term in the form of Eq. (6.25) we have exploited the fact that all higher-order corrections must begin with a dressed photon attached to the L µν tensor on one side, and to the W µν tensor on the other side. Again we must assume that the renormalisation and factorisation procedure has been carried out for the W µν tensor and the dressed photon, but their expressions must remain in D dimension until we perform the MS subtraction of the divergence arising from the q integration. Our final result is We now notice that the σ SF term corresponds exactly to f PF γ , while the σ (HO) PM term is directly related to f con γ . The first correspondence follows from the fact that the θ(µ 2 −|q 2 ⊥ |) condition is equivalent to θ(µ 2 − Q 2 (1 − z)). In fact, using Eq. (6.11), we find 27) up to corrections of order m 2 p /µ 2 . The second correspondence is slightly more subtle, since in f con On the other hand, the unrestricted integration evaluated in the parton model limit (i.e. neglecting all masses) is zero in perturbation theory, since it leads to scaleless integrals which vanish in dimensional regularisation. Thus demonstrating the equivalence of the two procedures.

The photon TMDPDF
The analysis of the previous section can readily be generalised to obtain the photon transverse momentum dependent PDF (TMDPDF) f γ (x, k ⊥ , µ 2 ). The TMDPDF (see Ref. [37] for a review on TMDPDFs) is given by Eq. (6.1), with the two field-strength tensors separated by x ⊥ and Fourier transforming in x ⊥ , In momentum space, it is given by Eq. (6.3) with an insertion of (2π The derivation proceeds as in the previous section. The q integral is written as and the q ⊥ and q + integrals are done using the δ-functions, leaving only the q − integral.
In terms of the variable z = x/w one has so that the q − integral can be replaced by one over z. Defining this leads to Eq. (6.33) is an exact expression for the TMDPDF if one-photon-irreducible structure functions are used. The TMDPDF is connected with our physical factorisation component of the collinear PDF through the following simple relation: The Photon PDF to higher order In the previous sections, we obtained the photon PDFs f γ and f ∆γ up to corrections of order α(µ)α s (µ) and α 2 (µ), as given in Eqs. (3.26) and (5.10). In Sec. 6, we wrote an exact formula for the photon PDFs, Eq. (6.14), and we expressed it as the sum of two terms: f PF γ,∆γ in Eq. (6.16) and f con γ,∆γ in Eq. (6.19). The physical PDF f PF γ,∆γ is a finite integral that can be evaluated explicitly using measured DIS structure functions, while f con γ,∆γ is an infinite integral, that needs to be renormalised in the MS scheme to obtain the MS PDF. The renormalised f MS con γ,∆γ was computed in perturbation theory in Eq. (6.21) to order α(µ) to reproduce the earlier results in Eq. (3.26) and Eq. (5.10). In principle, one can extend the computation of f MS con γ,∆γ to arbitrarily high orders in perturbation theory. In this section, we do this to order α(µ)α s (µ) and α 2 (µ), and adopt our "democratic counting" for the order in perturbation theory. For the sake of simplicity, we will limit ourselves to the unpolarised case.
We start with the expression Eq. (6.19a) for f con γ (x, µ 2 ), which we can write in compact form as In the sum, I is over 2, L, while a is over all partons. The functions p I and r I are We stress that the expression [p I + r I ] is exact, i.e. it has no higher order corrections in . The coefficient functions F are defined in such a way that consistently with the notation of Ref. [38]. They are perturbatively calculable, since µ is a hard scale. We now focus upon the O(αα s ) term. We only need a ∈ {q, g}. As usual we define series expansions for F 2/L,a and f con γ , At lowest order the only non-vanishing coefficient function is For the computation of the O(α s α) corrections we need the coefficient functions to order α s , including terms of order , The "c.t." suffix on the square bracket is there to remind us that the enclosed expression is an MS counterterm. The C I,a and a I,a coefficients (taken from Ref. [38]) are given in App. G. We can substitute the coefficient functions in Eq. (7.1) and evaluate the integral. Note that so that the ds integration yields a 1/(2 ) for all terms of Eq. (7.8) with the exception of the terms in square brackets (the counterterms), where it yields a 1/ . The leading term is order α(µ)/(2π): Explicitly, the finite part is and gives the result obtained earlier in Eq. (6.21a) when F 2 there is replaced by its leading order approximation. The infinite pieces are cancelled by the counterterms.
The infinite terms of Eq. (7.13) are cancelled by the MS counterterms for the photon PDF. The finite part of Eq. (7.13) gives instead the MS-conversion factor at order αα s . In explicit form, we get: Differentiating the finite parts of the photon PDF at order αα s gives the αα 2 s splitting functions, following the procedure in Sec. 4. The expression is lengthy, involving triple convolutions, and is not given explicitly here.
The α 2 corrections are given by the same equations, with the QCD corrections to the splitting functions and to the coefficient functions replaced by the QED ones. The only new feature in the QED correction is the vacuum polarisation contribution from α ph , which gives an additional contribution The final result for the f (0,2) γ , reinstating the sum over all flavours, and remembering that also leptons contribute in the democratic counting, is for a ∈ {q, l} p qγ i∈{q,l} n i e 4 i for a = γ.
In summary, the final result for the photon PDF is 19b) The f PF γ contribution was given in Eq. (6.16a), while the f terms are given in Eqs. (7.12, 7.14, 7.17). For s ≥ 2, the r.h.s. contains the photon PDF f γ , and Eq. (7.19) can be solved iteratively for f γ . The QCD correction can be included by simply adding the f (r,s) γ term to the final result. Including QED corrections is more difficult. The formalism presented here uses the one-photon-irreducible structure functions. To include QED corrections in a systematic way without double-counting requires experimental knowledge of the elastic form factors and DIS structure functions in the one-particle irreducible definition. This requires removing two-photon-exchange (TPE) contributions to the scattering cross sections. However, electromagnetic corrections only on the hadronic side should not be removed. This makes the analysis simpler than the usual one presented in the literature, in which QED corrections on the hadronic side are also removed (see, e.g. [39]). Furthermore, α ph (Q 2 ) should also be evaluated including corrections of order α 2 (in our NLO result it was enough to include corrections of order α 2 L). This requires the knowledge of the hadronic vacuum polarisation, that, on the other hand, is extracted from e + e − data (see the review on electroweak physics in Ref. [40]).
In analogy with our analysis in Sec. 4, the f coefficients, together with the two-loop QCD and QED coefficient functions could be used to obtain the P (r,s) γi splitting functions for r + s = 3 and s ≥ 1. We note that the complementary P (2,1) qγ and P (2,1) gγ splitting functions have been given in Ref. [41].

Inputs for the unpolarised photon distribution
To evaluate the photon parton density we require information on the F 2 and F L structure functions over the full x, Q 2 kinematic range. Our evaluation will be up to the accuracy outlined in section 3, i.e. using L ∼ ln µ 2 /m 2 p , we include terms αL(α s L) n at lowest order and α(α s L) n , α 2 L 2 (α s L) n corrections at higher order. As we proceed we will highlight issues that arise if one wishes to go to higher accuracy.
The F 2 and F L structure functions are most commonly determined from electronproton scattering data. A first comment concerns the treatment of electromagnetic corrections to the structure functions. We use the prescription described in Sec. 2, whereby all electromagnetic corrections to the interaction of a proton with a single photon are included, with the exception of the photon vacuum polarisation contributions. The treatments of electromagnetic corrections to data differ in their details according to the kinematic region considered. However two features emerge: the experimental analyses always correct for electromagnetic radiation from the incoming electron and they always take out the photon vacuum polarisation contributions. These are the two elements that are required for consistency with our accuracy. A more detailed discussion of QED corrections to DIS measurements is given in section 9 of Ref. [42].
We separate the data inputs according to the kinematic region and the corresponding final state in ep scattering. The main kinematic variables for the separation will be Q 2 and is the squared invariant mass of the outgoing system associated with the hadronic side of the collision.

Elastic contribution
In our definition, the elastic contribution corresponds to the region of W < m p + m π 0 .
In particular it includes configurations where one or more photons are radiated from the proton. 8 Experimental data on elastic scattering is usually corrected for radiation from the 8 For the determination of the structure functions we find it useful to think of a process in which there can be at most one exchanged photon between the probe and the proton, as in the process of Sec. 3. In actual electron-proton scattering experiments there can be two or more exchanged photons, either real of virtual. These corrections are beyond our accuracy and cannot be classified in terms of the usual electromagnetic structure functions, since they correspond to a more complex tensor structure.
proton, since the measurements are performed with the goal of extracting the electric and magnetic Sachs form factors of the proton, G E and G M respectively. The correspondence between the form factors and the structure functions (see e.g. Eqs. (19) and (20) of Ref. [43]) is given by where τ = Q 2 /(4m 2 p ). The approximation of F 2 and F L as δ-functions neglects precisely the photon radiation from the proton. However, most of the radiation is soft and so cancels in inclusive quantities and all that changes when going beyond the δ-function approximation is a relative O (α) correction (free of any logarithms) to the photon distribution, which is beyond our accuracy. Collinear logarithms are absent since the proton form factors fall off rapidly with Q 2 .
Substituting Eq. (8.2) into Eq. (3.26), one obtains for the elastic component of the photon distribution. Note that we have set the upper limit of the integration to infinity, rather than make it µ 2 -dependent as in Eq. (3.26). For large µ 2 values, because of the 1/Q 4 scaling of the form factors (see below), the resulting difference is a higher-twist effect. Using an infinite upper limit ensures the absence of higher-twist contamination. For the same reason, the last term in Eq. A widely used approximation for the G E,M form factors is the dipole form, where m 2 dip = 0.71 GeV 2 and µ p 2.793 is the anomalous magnetic moment of the proton. For Q 2 = 0 this form yields the exact results G E (0) = 1 and G M (0) = µ p , while elsewhere it is an approximation.
The dipole form is of interest in part because it provides insight into the asymptotic behaviours of the elastic component of the photon distribution. At small x, to leading order in α, one finds   For accurate results, the dipole approximation, Eq. (8.4), is not sufficient. The most recent extensive experimental study of the form factors comes from the A1 collaboration [39]. The A1 data itself is limited to Q 2 1 GeV 2 , however the work includes fits to global data up to Q 2 ∼ 10 GeV 2 . The electric and magnetic form factor fits are shown in Fig. 6, normalised to the dipole form. The A1 paper includes two classes of fits, one for just unpolarised data, and one that also includes polarised data. 9 Both fits show clear deviations from the dipole form. The fit with polarised data is the recommended default because it provides additional constraints on two-photon-exchange (TPE) contributions and we take it as our default. Given the delicacy of the treatment of the TPE contribution, we use the central value of the unpolarised fit as an error estimate that comes in addition (in quadrature) to the quoted uncertainty on the global polarised fit. Note that the fits extend only up to Q 2 = 10 GeV 2 and beyond this point we use the dipole shape, normalised  Figure 7. Elastic contribution to f γ (x, Q 2 ) using various fits for the form factors, normalised to the result obtained with the A1 world fit including polarised data. The ratio for the the A1 world fit freezes above x = 0.9 because the A1 fits extends only up Q 2 = 10 GeV 2 and beyond that scale we simply extrapolate the results for G E/M (Q 2 ) using the standard-dipole shape, normalised to G E/M (10 GeV 2 ).
to the fitted G E/M (Q 2 ) at Q 2 = 10 GeV 2 . We treat the fit uncertainties on the elastic and magnetic components as 100% correlated, which is the most conservative assumption, because they both enter with the same sign in Eq. (8.3). The impact of the A1 fits on the elastic contribution, relative to the dipole form, is shown in Fig. 7. The deviation from the dipole form is manifest, with an impact of up to 10% for x 0.5. In this x-range, the difference between the results using fits with and without polarised data is smaller than the intrinsic uncertainty bands for the individual fits. Note that since the data constrains the form factors for Q 2 10 GeV 2 , the elastic contribution is effectively known only for x 0.9, cf. the form of the lower limit in Eq. (8.3).
We have also computed the elastic contribution using a preliminary unpublished fit to the elastic form factors by Lee, Arrington and Hill. 10 There is some change in the elastic contribution to f γ , which is negligible for small x, and becomes comparable to our overall total error (for the elastic plus inelastic components) for x 0.5.

Low-Q 2 region
In addition to the elastic component, we need the inelastic part. This corresponds to the region of W > m p + m π 0 . We split the inelastic part into several sub-regions. At low Q 2 , the structure functions cannot be computed from parton distribution functions and we must rely on direct measurements and theoretical model-independent constraints.
Results for low-Q 2 scattering tend to be given either as a differential scattering cross section or in terms of some subset of F 2 , F L , R L/T , σ T and σ L , only two of which are independent. They are related through the following equations (using the same normalisation convention adopted by the HERMES collaboration in Ref. [45]), As discussed in App. E, F 2 and F L vanish respectively as Q 2 and Q 4 in the limit of small Q 2 and fixed W 2 . In this limit, as one can see from Eq. (8.1), x goes to zero in proportion to Q 2 . The quantity σ T (x, Q 2 ) becomes a function W 2 only and it follows from Eq. (8.7) that σ L (x, Q 2 ) must vanish as Q 2 , leading to where σ T (W 2 ) is the photoproduction cross section. As a result, in Eq. (3.26), Q 2 scales much below the proton mass will not give a substantial contribution. There is wealth of data covering the low Q 2 region, including also photoproduction data through the constraint Eq. (8.8). Rather than using these data directly, we will rely on existing fits of those data. Fits generally focus upon either the resonance region, W 2 3 GeV 2 or the continuum region, W 2 4 GeV 2 . Fig. 8 (left) shows data from the CLAS experiment [46], compared to two global fits to resonance region data. The figure (which includes only a small subset of the available data) illustrates the coverage in Q 2 and the quality of the available data. The data is shown as a function of W 2 in order to clearly show the resonance peaks starting with the ∆ resonance and beyond. The CLAS fit is intended for use only for Q 2 > 0.5 GeV 2 , while the Christy-Bosted [47] fit is intended for use down to Q 2 = 0 and explicitly includes photoproduction data. Comparing the   Figure 8. Left: illustration of a subset of the CLAS data [46] in the resonance region, compared to two fits, one from the CLAS paper and the other from Christy and Bosted (CB) [47]. The CLAS dataset covers the range 0.225 ≤ Q 2 /GeV 2 ≤ 4.725 in steps of 0.05, with a quality comparable to that shown in the plot across the whole Q 2 range. The errors on the data correspond to the sum in quadrature of statistical and systematic components. The CLAS data is only a small part of the data that is available in the resonance region and used for the fits (see also Fig. 6 of Ref. [46]). Right: illustration of the GD11-P fit from the HERMES collaboration [45] and corresponding data in the continuum region (plot reproduced with kind permission of the HERMES collaboration).
CLAS and Christy-Bosted fits at Q 2 values below the quoted validity range of the CLAS fit shows, however, that they are relatively similar. For the continuum region, the HERMES collaboration has provided a fit, GD11-P [45], using a wide range of data and the ALLM [48] functional form. Fig. 8 (right), taken from Ref. [45], illustrates the good quality of the fit. Careful inspection of the figure reveals that at each Q 2 value the fit consists of three lines, whose separation represents the uncertainty.
Electron-proton scattering cross sections give information on both F 2 and F L , with the former often dominating the cross section. As a result the knowledge of R L/T in Eq. (8.7a) is often much poorer than the knowledge of F 2 . Some of the F 2 fits (CLAS, GD11-P) rely on earlier independent fits for R L/T , notably the R 1998 [49] fit and that of Ref. [43]. Instead the Christy-Bosted article carried out an independent fit for R L/T . Fig. 9 (left) shows the R 1998 and Christy-Bosted fits for R L/T , compared to a subsequent extraction of R L/T by the E94-110 collaboration [50]. The moderate degree of agreement between data and fits  Figure 9. Left: the R 1998 and Christy-Bosted fits for R L/T as a function of W 2 , compared to data from the E94-110 experiment at JLAB [50]. The data corresponds to version 2 of the arXiv submission of the E94-110 article, which postdates the Christy-Bosted fit. Also shown is the uncertainty band that we adopt for the R 1998 fit. Right: the Q 2 dependence for the R 1998 fit at various W 2 values, including our extension to low Q 2 , Eq. (8.9).
will motivate us to assign a generous ±50% uncertainty to R L/T in the low-Q 2 region. As a default we will adopt the R 1998 fit for Q 2 > Q 2 b with Q 2 b = 0.34 GeV 2 . The fit is, however, not intended for use at low values of Q 2 , and so for Q 2 < Q 2 b we will use the form that is continuous and reasonably smooth at Q 2 b and vanishes as Q 2 for Q 2 → 0. This is shown in Fig. 9 (right).
No single low-Q 2 fit for F 2 is simultaneously adequate in the resonance and continuum regions. Accordingly we will combine F 2 resonance (F res 2 ) and continuum (F cont 2 ) fits using two transition scales W 2 lo = 3 GeV 2 and W 2 hi = 4 GeV 2 where W 2 is evaluated as in Eq. (8.1) and ρ(W 2 ) is This ensures a continuous and smooth transition between the low-and high-W 2 regions. We will consider two combinations of low and high-W 2 fits for F 2 : our main one will be CLAS with GD11-P, using R 1998 for both (including the modification Eq. (8.9)). As a  Figure 10. The sum of σ T +σ L , cf. Eq. (8.7b), shown as a function of W 2 (lower axis) and x (upper axis) for various Q 2 values. Our default prescription corresponds to the curve labelled "GD11-P + CLAS", which transitions smoothly between the GD11-P and CLAS fits in the region of W 2 between 3 and 4 GeV 2 .
cross-check of uncertainties we will have a subsidiary combination, composed of Christy-Bosted and GD11-P for F 2 , and Christy-Bosted R L/T at low W 2 and R 1998 at high-W 2 , with the two R L/T fits combined in a manner identical to Eq. (8.10).
Our main combination is shown in Fig. 10, compared to the individual CLAS, Christy-Bosted and GD11-P results. The need for a combination of more than one fit is evident from the fact that the GD11-P fit misses the resonance structures below 4 GeV 2 , while the resonance fits miss the rise at large-W 2 . For all Q 2 values shown, the transition region between 3 and 4 GeV 2 is reasonably covered. Note also the reduced significance of the resonance peaks at high-Q 2 , though the Christy-Bosted fit appears to show artefacts in this regard near W 2 = 3 GeV 2 . For this reason, when we use the Christy-Bosted fit in the region of Q 2 > Q 2 0,CB = 8 GeV 2 we will modify its large-Q 2 behaviour as follows Figure 11. Value of the structure functions F 2 (left) and F L (right) divided by Q 2 as a function of x bj and Q 2 , using our CLAS+GD11-P prescription. They are divided by Q 2 to reflect the structure of the integrand for f γ (x, µ 2 ), Eq. (3.26). The latter receives contributions only from the region from x bj > x. For a full understanding of the contribution of different regions of x bj and Q 2 to the integral, one must also take into account the z = x/x bj -dependent factors that multiply the structure functions and the exact limits in x bj and Q 2 .
with Q 2 1,CB = 30 GeV 2 . The Christy-Bosted curves shown in Fig. 10 do not include this modification.
To close our discussion of the low-Q 2 region, we show in Fig. 11 F 2 (left) and F L (right) divided by Q 2 , as a function of x bj and Q 2 . This gives an indication of the size of the contribution of different x bj and Q 2 values in Eq. (3.26). The main features to note are the importance of the resonance region and the relevance of Q 2 values down to zero over nearly the whole x bj range. One also sees that F L brings a substantially smaller contribution than F 2 , and tends to vanish for small Q 2 , as dictated by our use of Eq. (8.9). The figure focuses on the region Q 2 < 0.6 GeV 2 . For larger values of Q 2 , the usual scaling behaviour of F 2 sets in, with, to first approximation, a uniform contribution in ln Q 2 .

High-Q 2 continuum
For sufficiently large Q 2 and W 2 one can calculate F 2 and F L from parton distribution functions (PDFs) using the known perturbative expansion of the DIS coefficient functions. This is more reliable than using a fit to available data (e.g. GD11-P also includes some high-Q 2 data), because the extension to arbitrarily large Q 2 is provided by DGLAP evolution rather than the extrapolation provided by some a priori arbitrary parametrisation. Furthermore, in recent years there has been extensive progress in the extraction of PDFs from DIS and collider data, including detailed and well-tested uncertainty estimates.
Our prescription for evaluating F 2 and F L is as follows. Our choice of PDF and associated uncertainties will be PDF4LHC15 nnlo 100 [51]. This is based on a combination of the CT14nnlo [52], MMHT2014nnlo [53] and NNPDF30 nnlo as 0118 [54] global PDF fits. We will use NNLO coefficient functions [38,[55][56][57], implemented in HOPPET [58][59][60]. All quark flavours will be treated as massless and we will correspondingly use a zero-mass Figure 12. Values of the structure functions F 2 and F L as a function of x bj and Q 2 , using a PDF4LHC15 nnlo 100-based NNLO ZM-VFNS prescription for Q 2 > 9 GeV and W 2 > 4 GeV 2 , and the CLAS+GD11-P combination elsewhere.
variable flavour-number scheme (ZM-VFNS). The heavy flavour contribution to F 2 and F L is of order α s and by taking the massless approximation for a quark of mass m q we mistreat that order α s contribution in a region of Q 2 ∼ m 2 q . Examining Eq. (3.26), one sees that this will translate to an inaccuracy in f γ of order αα s , i.e. beyond the order that we aim to reproduce in our analysis with data.
Determinations of structure functions from PDFs are potentially subject to corrections from higher-twist effects. Recent studies of DIS data [61,62] suggest that higher-twist effects could be substantial for F L , at least at small values of x bj . In line with the observations of this study, we will account for such a possibility when evaluating our uncertainties, using a multiplicative correction taking the larger of the two normalisations in Refs. [61,62]. Our default domain for using a PDF-based evaluation of the structure functions will be Q 2 > Q 2 PDF = 9 GeV 2 and W 2 > 4 GeV 2 . In the rest of the kinematic plane we will use the resonance and low Q 2 continuum fits. The breakup of the x bj −Q 2 plane is summarised in Fig. 12, analogous to Fig. 11, but showing a larger range of Q 2 . The Q 2 scale is now logarithmic and we no longer divide F 2 and F L by Q 2 . The colour-coding once again provides a visualisation of the density of the integrand in Eq. (3.26). At large x bj one sees the gradual reduction with increasing Q 2 of the structure functions. A consequence of this is that the resonance part of the low-Q 2 region plays an especially important role in the determination of the large-x photon distribution. At moderate x bj , the structure functions are largely independent of x bj , i.e. they display Bjorken scaling, while at small x bj the structure functions increase rapidly with Q 2 , a consequence mainly of double logarithms of x bj and Q 2 in the scaling violations.
Careful inspection of Fig. 12 reveals that F 2 and F L are not perfectly continuous at the transition scale of Q 2 trans = 9 GeV 2 , i.e. the results of the GD11-P fit and PDF-based structure function evaluations do not quite match up. This is most visible if one inspects the band of yellow colour in the plot for F 2 , around x bj = 2 · 10 −5 . To probe the impact of this discontinuity in our f γ determinations we will introduce, as one of our uncertainty sources, a variation of the Q 2 threshold for switching between the low-Q 2 GD11-P fit and the high-Q 2 PDF-based determination. The alternative Q 2 threshold that we will use is Q 2 trans = 5 GeV 2 .

Uncertainty from missing higher orders
The integral in our photon PDF was broken up into two pieces, an all-Q 2 part which gave the "physical" PDF, and a high-Q 2 which gave the MS-conversion piece, cf. Eqs. (6.16a, 7.19). In this section, we examine the uncertainty in the final photon PDF result depend on the precise way in which the PDF integral is split up, i.e. on the upper limit of the Q 2 integral in Eq. (6.16a).

Dependence on physical factorisation separation scale
To probe the impact of the choice of separation between the physical PDF and the MSconversion, we use the following alternative definition of the physical and MS-conversion terms where [M ] denotes the functional dependence on M (z). We easily find that In the following we will consider two alternative choices for M (z): and The variation of the answer as a function of µ M will provide us with an estimate of uncertainties from missing higher-order (QCD) contributions. This is similar in spirit, but different in its details, from standard scale variation. M (z) is a large scale, where we can use QCD expressions for the structure functions in terms of parton distributions. We have where, up to NLO QCD accuracy for the coefficient functions,     (x, µ 2 ) term, Eqs. (7.12) and (9.7), in its series expansion, Eq. (7.6). The NNLO result additionally includes the f (1,1) γ (x, µ 2 ) term, Eqs. (7.14) and (9.8). Higher-order QED corrections are not included. The high-Q 2 contribution to the f PF term is always evaluated including a NNLO PDF (PDF4LHC15 nnlo 100 [51]), NNLO splitting kernels [63,64] and massless NNLO QCD coefficient functions [38,[55][56][57] as implemented [59,60] in HOPPET [58].

MS-conversion terms:
The upper plot shows the substantial reduction in the scale uncertainty in going from LO to NLO with both scale choices. There is an especially large uncertainty in the LO photon PDF for small x. In this region the photon PDF is dominated by the high-Q 2 integration due to the DGLAP evolution of the quark distribution. Changing the scale by a factor of two around our central value of µ = 100 GeV causes a 20% shift in ln(µ/Q trans ). This helps explain the size of the uncertainty in the upper figures. Recall that Q trans = 3 GeV is the scale above which we use PDFs and coefficient functions to evaluate the structure functions, as explained at the end of Sec. 8.3. Previous estimates of the photon PDF were only accurate to LO, and so had an intrinsic uncertainty comparable to that of our LO band. 11 In the lower left-hand plot, the curve labelled FSF-NLO ("Full Structure Function") corresponds to Eq. (3.26), with F 2 computed at NNLO in both the physical-factorisation and the MS-conversion terms. The PRL variant of the FSF-NLO curve corresponds to the result of Ref. [22], whose implementation had a formally subleading bug in the scale of α that multiplied the MS-conversion term. The corresponding curves in the lower right-hand plot are analogous but with M 2 (z) = µ 2 as the upper limit of the physical-factorisation component.
For our final central result, to be shown in the next section, we will take the NLO curve with M 2 (z) = µ 2 /(1 − z), i.e. using f con γ (x, µ 2 ) that includes just the f (0,1) γ (x, µ 2 ) term, Eq. (7.12). As an estimate of the uncertainty from missing higher-order (MHO) contributions, we will take the largest deviation of any of the NNLO M 2 (z) = µ 2 M /(1 − z) scale choices and symmetrise it with respect to the central NLO result.
This differs from the prescription adopted in Ref. [22], where we used the FSF-NLO evaluation with M 2 (z) = µ 2 /(1 − z). There the MHO uncertainty was estimated from the difference between FSF-NLO evaluations with M 2 (z) = µ 2 /(1 − z) and M 2 (z) = µ 2 , corresponding to the difference between the orange (dot-dashed) curves in the left and righthand plots of Fig. 13. The difference between M 2 (z) = constant and M 2 (z) ∝ 1/(1 − z) is a log(1 − z) contribution, which gets large as z → 1. These "endpoint" logarithms arise because a new scale Q 2 (1 − z), the invariant mass of the final state hadronic system, enters the problem. While formally, this spurious logarithm is cancelled by a corresponding one in the MS-conversion term, there is a left-over higher-order log(1 − z) with the choice M 2 (z) = constant. For our default photon-PDF prediction, we prefer not to add a spurious endpoint logarithm to f con γ via our choice of M 2 (z). Figure 14 compares the old and new versions of our photon PDF.
Note that we have chosen not to use a full NNLO result. While this is in principle straightforward to achieve, several further practical elements would be needed beyond those we have implemented so far. The one that is potentially most problematic relates to the fact that one must take into account QED corrections to the structure functions. This would call for a detailed assessment of the way in which QED corrections have been removed from published data. 11 One should keep in mind that many current predictions for processes involving an incoming photon are only available at LO for the photon-induced process. They will therefore have a substantial relative scale uncertainty. This uncertainty should mostly be cancelled in calculations that include NLO QED corrections to the photon-induced processes, e.g. Refs. [16,65].  Figure 14. Ratio of the photon PDF in this paper (LUXqed17), shown as a blue solid line, to our previous result [22] (LUXqed) at µ = 100 GeV. The solid orange area is the original error band, and the dashed blue lines indicate the new error band.

Evaluation of the photon distribution and prescription for other partons
We evaluate the photon distribution using Eqs. (6.16a, 7.12, 7.19), keeping only the f (0,1) γ term in Eq. (7.19b). The elastic component is given by Eq. (8.3). A number of choices need to be made in the evaluation, both for the central value of the photon distribution and its uncertainties. For practical usage one also requires a consistent set of parton distribution functions for the other flavours. In this section we describe the approach we take and considerations for future evaluations. We recall that, as in Sec. 8, our accuracy aim for the photon distribution will be to control terms up to α(α s L) n and α 2 L 2 (α s L) n .

Evaluation scale and evolution
A first choice relates to the question of the µ 2 value(s) where the photon distribution is evaluated and the issue of higher-twist corrections. When µ 2 is not very large, the upper limit of Eq. (6.16a) is in a region where both the inelastic and elastic parts of structure functions may themselves contain higher-twist corrections.
For the elastic part, in Eq. (8.3) we have deliberately set the upper limit in the Q 2 integration to infinity. This is possible because the contribution to the µ 2 dependence of the elastic part comes from the running of the QED coupling in front of Eq. (8.3), which is not affected by the choice of upper limit in Q 2 . Leaving the upper limit as µ 2 /(1 − z) would instead have resulted in higher twist contributions associated with the 1/Q 4 scaling of the form factors.
In the inelastic part, there is no easy way, at low µ 2 scales, of separating leading and higher twist components. This is because at the corresponding upper limit of the Q 2 integration, the inelastic structure functions mix both leading and higher-twist effects.
Accordingly we choose to evaluate Eq. (7.19) at a scale µ 2 eval m 2 p and then determine the photon distribution at other scales through DGLAP evolution. We take µ eval = 100 GeV so as to ensure that potential residual higher-twist effects, of order m 2 p /µ 2 eval , are much smaller than the precision we will be seeking, which will be roughly at the 1% level.
For the DGLAP evolution we will include the O (α) and O (αα s ) splitting kernels [26], as well as the QCD kernels up to O α 3 s [63,64]. Given the precision that we use in the evaluation of the photon distribution, and our assumption α 2 s ∼ α, this forms a consistent set of evolution terms.

Uncertainties
The final uncertainty on our PDF distribution is taken to be the sum in quadrature of many individual uncertainty sources, because they are uncorrelated. The individual contributions have already been discussed in Sec. 8, where we explained the input data to our analysis. The various uncertainties with labels as in Fig. 15 are: (EFIT) The uncertainty on the elastic contribution that comes from the uncertainty on the A1 world polarised form factor fits, as shown in Fig. 7. This band is asymmetric and we symmetrise it using the largest deviation to obtain a (more conservative) symmetric band.
(EUN) The uncertainty that comes from replacing the A1 world polarised fit (which includes a two-photon-exchange correction) with just the world unpolarised data (which does not). This provides a one-sided uncertainty, which we again symmetrise.
(RES) We replace the CLAS resonance-region fit with the Christy-Bosted fit (modified as in Eq. (8.12a). This replacement gives a one-sided uncertainty, which we once again symmetrise.
(R) A modification of R 1998 L/T by ±50% around its central value, as shown in Fig. 9. Recall that this R choice is only used in the regions where we take F 2 from one of the GD11-P, CLAS or Christy-Bosted fits. Of the ±50% variations of R, the one with the larger impact on the photon distribution is identified and the resulting uncertainty symmetrised.
(M) A modification of the Q 2 PDF scale which governs the transition from the GD11-P structure function fit to a PDF-based evaluation. The default choice of Q 2 PDF = 9 GeV 2 is reduced to 5 GeV 2 and since this is a one-sided uncertainty, the resulting effect is symmetrised.
(PDF) The input PDF uncertainties for Q 2 > Q 2 PDF according to the default prescription for the PDF (PDF4LHC15 nnlo 100).
(T) A twist-4 modification of F L as in Eq. (8.13). This is a one-sided modification that is then symmetrised. Note that we do not include the quoted uncertainty from the GD11-P fit. When studying that uncertainty we found on one hand that its impact was negligible compared to the other uncertainties, and on the other hand that the resulting uncertainty band did not always overlap with F 2 as calculated from PDFs in regions at small x and high Q 2 that lack direct F 2 data. This observation motivated our choice to restrict the use of the GD11-P fit to Q 2 values with sufficient data coverage and to vary the Q 2 trans transition scale (and R) for the uncertainty estimate. We also did not include any uncertainty associated with the prescription for matching other parton flavours.
The impact of the different sources of uncertainty is shown in Fig. 15, and our final uncertainty, shown by the black line, is given by adding the contributions in quadrature. 12 The overall uncertainty is less than 2% for 10 −4 < x < 0.1. For small values of x, the uncertainty is dominated by the uncertainties in the parton distributions of quarks (and gluons), which enter the high-Q 2 part of the photon PDF integral. As x → 1, the uncertainties are dominated by the low-Q 2 part of the photon PDF integral from elastic form factors, the resonance contribution, and σ L . This is a reflection of the fact that non-perturbative effects (such as higher twist corrections) grow like Λ 2 QCD /[Q 2 (1−x)] as x → 1, and that, for x close to 1, quark PDFs fall off rapidly as Q 2 increases, so the low-Q 2 region contributes significantly to f γ as x → 1.

Matching to other partons
The introduction of a photon component of the proton necessarily implies modifications of other partons relative to a set that has been determined without QED corrections. QED corrections to the quarks start at O (αL(α s L) n ), i.e. the same order as the leading photon contribution. These terms are generated by the order α QED contribution in DGLAP evolution. If one wishes to know the full set of partons to the same accuracy as the photon, i.e. O (α(α s L) n ) relative corrections, then it is effectively necessarily to have a reliable estimate of the QED corrections to the initial conditions for the DGLAP evolution. This would require that one repeat a global PDF fit with QED corrections, e.g. the order α corrections to the DIS coefficient functions, and with information about the photon PDF as an input (e.g. because it affects the momentum sum rule at order α at the starting scale). Such a fit is beyond the scope of our article, though below we outline a procedure for how it might be performed.

Our procedure
The prescription that we adopt is as follows. At a scale µ 2 match , we assume that the quarks are unchanged relative to a global fit without QED contributions. At this scale we rescale the gluon by a factor as follows: where ω i (µ 2 ) is the momentum fraction carried by parton flavour i at scale µ 2 , As we shall see below, Fig. 20, the photon momentum fraction, ω γ , is a fraction of a percent. The above procedure ensures that the momentum sum rule, including the photon contribution, is satisfied. The reason for absorbing the momentum into an adjustment of the gluon is that the gluon is the parton least directly constrained from DIS data. The choice of µ 2 match is somewhat arbitrary. Ideally, it should be close to the Q 2 scales in DIS that provide the greatest constraint on the quark PDFs, given that our procedure leaves the quark distribution unchanged at scale µ 2 match . Since we use a PDF fit that was determined without QED corrections, when we convolute with the QCD coefficient functions (without QED corrections), we should reproduce the true experimental F 2/L at the scales where the DIS data is most constraining. There is however a practical problem that we should consider. The PDF4LHC15 nnlo 100 set, which we take as our base PDF set, is designed for use at large µ 2 values. At low µ 2 we encountered two issues, discussed in detail in App. H: one is that it is a merger of sets with different underlying heavy-quark thresholds, and as a result does not strictly satisfy DGLAP evolution. The second relates to the way the PDF4LHC15 nnlo 100 underlying PDF sets are encoded in LHAPDF [20] files. These issues prevent us from starting a DGLAP evolution from scales below about 6 GeV. Putting together the various considerations, we opted for the choice µ match = 10 GeV. The   Figure  16.
Ratio of the gluon and u-quark PDF in this paper (LUXqed17 plus PDF4LHC15 nnlo 100) to the corresponding PDF4LHC15 nnlo 100 [51] distributions at µ = 100 GeV, shown as the solid blue line. The solid orange region is the original PDF4LHC15 nnlo 100 error band, and the dashed blue lines represent our error band. In the case of the up-quark distribution, we also show in green the impact of using µ match = 6 GeV instead of our default of 10 GeV. resulting set, LUXqed17 plus PDF4LHC15 nnlo 100, is valid only for scales µ ≥ 10 GeV (below this scale, LHAPDF interpolates the LUXqed17 partons to zero).
Since we include QED effects in the DGLAP evolution, at scales other than µ match all partons acquire QED-induced modifications, of order (αL) n (α s L) m for the quarks, with n ≥ 1, m ≥ 0. The changes in the gluon and up-quark PDFs at scale µ = 100 GeV are illustrated in Fig. 16. In the case of the up-quark PDF, we also show the effect of reducing µ match from our default of 10 GeV to 6 GeV, demonstrating that the impact is minimal compared to the overall uncertainty on the up-quark distribution.
As discussed in Sec. 10.1, it is not advisable to directly use Eqs. (6.16a, 7.12, 7.19) to evaluate the photon PDF at scales as low as µ match , because of the presence of higher twist effects that are difficult to control. On the other hand, if we evaluate the photon PDF at a scale µ eval µ match , the O (αL) modifications to the quark distributions from DGLAP evolution must be taken into account in order to correctly treat contributions of order (αL) 2 (α s L) n (n ≥ 0) in the photon distribution. We do so technically as follows, keeping in mind that our base PDF4LHC15 nnlo 100 PDF is based on fits that have neither QED evolution, nor QED corrections in coefficient functions. First we evaluate the photon distribution at scale µ eval µ match , using a high-Q 2 F 2 calculated without QED effects. We then evolve the photon distribution down to the scale µ match , using a special variant of DGLAP evolution in which all QED contributions to P qi and P gi are set to zero, for all parton flavours i. With this procedure, the quark and gluon densities remain identical to those of the original PDF set. These unchanged distributions are then used in the evaluation of P γq terms for the photon evolution. We perform the matching at scale µ match , as described above, and finally evolve back up in scale using DGLAP evolution including the full set of QED contributions. We stress that when evolving back up to the µ eval scale, all partons will acquire corrections of relative order (αL) n , with n ≥ 1. In particular the photon  Figure 17. The impact on the photon distribution from leaving out QED effects in the quark evolution (blue curve), compared to the LUXqed17 plus PDF4LHC15 nnlo 100 uncertainty. The irregularities in the curve at large x values are an artefact associated with the use of different underlying x grids between the results with and without QED effects in the quark evolution.
PDF will acquire relative correction of order αL, that are required by our aimed accuracy. As illustrated in Fig. 17, without the QED effects in the quark evolution (which tend to reduce the quark distribution), the photon distribution comes out slightly higher. However, the effect is minimal compared to the overall LUXqed17 plus PDF4LHC15 nnlo 100 uncertainty. Note that an O (1) change in the choice of µ match corresponds to a NNLO, O α 2 L(α s L) n , effect on the photon distribution. The smaller effect in the photon density in Fig. 17 relative to that in the quark densities in Fig. 16 (right) is because the higher-order QED modification of the photon is driven by the average of the QED modification of the quarks across all scales from µ match up to µ. However the quarks are modified mostly at scales close to µ, i.e. the average modification is suppressed.
A final comment concerns lepton distributions. For consistency with the running of the QED coupling, photon splitting to leptons should be included, which generates non-zero lepton distributions. These are of order (αL) 2 . Their feedback on the photon distribution is an (αL) 3 effect and so beyond our accuracy. 13 In practice we set the lepton distributions to zero at scale µ match and include the full set of P lγ , P ll and P γl splitting functions in the evolution. A more complete approach would determine the lepton distributions directly from F 2 and F L , in a manner analogous to that used for the photon distribution in this paper. Note that while we include leptons in our DGLAP evolution, their distributions are not included in the LHAPDF files that we make available.

Iterative procedure
The above procedure allows us to obtain NLO QED accuracy for the photon distribution and LO QED accuracy for the other partons. Recall that LO corresponds to αL and NLO to an extra suppression of αL or α s , with α s L considered to be of order 1. For compactness, below we will use δ to refer to any quantity of order of αL or α s . The photon distribution and QED effects in the quarks both start at order δ.
To go beyond the NLO QED accuracy for the photon, and to achieve better than LO QED accuracy for the other partons, requires a more sophisticated procedure than that discussed in the previous section.
One conceivable difficulty in constructing a general approach is that the photon evaluation depends on one's knowledge of the quark distributions, while the extraction of the quark distributions from data is itself affected by the photon distribution.
This apparent circularity can be circumvented using an iterative procedure that exploits the smallness of the QED coupling, or equivalently of δ. First one determines a photon distribution, f (0) γ as in Sec. 10.3.1, based on QCD parton distributions obtained from a PDF fit without QED effects, f no-QED q . As discussed in the previous subsection, the resulting photon PDF includes terms up to the order (αL) 2 (α s L) n ∼ δ 2 and α(α s L) n ∼ δ 2 , i.e. of relative order δ with respect to the leading term. One then carries out a global fit with QED effects in the evolution and coefficient functions, 14 and a photon distribution f (0) γ that is not changed during the fit. This gives a zeroth iteration of the QCD parton distributions with QED effects, f (0) q . Their accuracy depends on the hard processes being used in the QCD fit: • For Drell-Yan production, i.e. qq → + − and γγ → + − , the relative contribution from the photon distribution starts at order (f γ /f q ) 2 ∼ (αL) 2 ∼ δ 2 (the γγ and qq Drell-Yan production channels differ only by the sizes of the incoming parton distributions). Since the photon distribution is known up to relative order δ, the QED corrections to the quark distribution are known up to order δ 3 . In general if the photon distribution is known up to relative order δ n , the quark distribution can be extracted to order δ n+2 .
• For DIS structure functions, the relative contribution from the photon distribution starts at order αf γ /f q ∼ α 2 L ∼ δ 3 . Since the photon distribution is known up to relative order δ, the QED corrections to the quark distribution are known up to order δ 4 . In general if the photon distribution is known up to relative order δ n , the quark distribution can be extracted to order δ n+3 .
It is the Drell-Yan process, an important input in all modern global PDF fits, that will limit the accuracy of the iteration, so we concentrate on that process. We now use the QCD partons f (0) q to determine an improved approximation to the photon distribution, f (1) γ . Since the QCD partons are known to QED accuracy δ 3 , this new estimate of the 14 We assume, for the principle of the demonstration, that the evolution and coefficient functions are known to all perturbative orders. We also assume that lepton distributions are treated correctly to the relevant order. photon will be accurate to order αL × δ 3 ∼ δ 4 . One can then insert f (1) γ into a renewed determination of the QCD partons to obtain f (1) q , which will be accurate to order δ 5 . In general, at iteration i, f (i) γ will be accurate up to and including order δ 2+2i , while f (i) q will be accurate up to and including order δ 3+2i . The convergence of this procedure is therefore fast, and the limiting consideration will be the availability of coefficient and splitting functions of the required order. 15,16 11 Numerical results In this paper, we have made some small changes to the procedure used to evaluate the photon PDF compared to our previous paper, as discussed in Sec. 9.2. Figure 14 in that section showed the difference in the PDF using the two methods. The change is very small, and well within our errors.
The photon PDF computed using the method of this paper is given in the left-hand panel of Fig. 18 for µ = 100 GeV. It is rescaled by a factor 10 3 x 0.4 /(1−x) 4.5 to facilitate the simultaneous study of different x regions. The plot includes a breakup into the different contributions discussed in Sec. 8. There is a sizeable elastic contribution, with an important magnetic component at large values of x. The resonance and continuum regions are also quantitatively relevant. The white line represents contributions arising from the Q 2 < 1 region of all the structure functions, including the full elastic contribution, and this serves to illustrate the importance of a proper inclusion of the low Q 2 region, given the accuracy we aim for. The PDF contribution, in the physical factorisation scheme, is from the bottom of the grey region to the blue dashed curve. The MS-conversion term, Eqs. (7.6, 7.12), is negative and corresponds to the difference between the blue dashed curve (PF result) and the top edge of the grey region (final full result in the MS scheme).
The right-hand plot of Fig. 18 illustrates how the components evolve when increasing the factorisation scale µ. The main change is associated with the log Q 2 growth of the "PDF" contribution and is most important at small x values, a consequence of the fact that the quarks distributions themselves increase rapidly with Q 2 at small x. The elastic, resonance and (low-Q 2 ) continuum contributions to the photon PDF all depend slowly on µ via the overall 1/α(µ 2 ) factor in Eq. (6.16a). These components, though formally NLO, remain a significant fraction of the overall photon PDF, even at this large value of µ.
The right panel in Fig. 18 also shows the impact of scale variation on the contributions at µ = 500 GeV. The blue dashed curves are for µ M = µ/2 and µ M = 2µ in Eqs. (9.1,9.4a) for the PF photon. The total MS photon PDF (the top edge of the grey band) uses our 15 As discussed in Sec. 7, starting from order α 2 = δ 4 , the photon distribution is itself a direct input to its own determination; an iterative type approach also addresses this case. 16 The above discussion assumes that PDF fits for the QCD partons are carried out without an imposed momentum sum rule and are constrained exclusively by data. However, it is quite common for the momentum sum to be imposed. In this case the momentum carried by the photon acts as an effective input to the fit, since it determines the momentum allowed for the QCD partons. For the f (0) q determination, the highest known QED term in the momentum sum is δ 2 . Thus, insofar as the overall accuracy of the fit is limited by the input(s) with worst accuracy, after one iteration the absolute accuracy for both f   central choice of µ M = µ in Eq. (9.4a). The impact of a change of µ M on the MS photon PDF would be barely visible in the plot, because the substantial scale dependence of the PF result is largely cancelled by that of the MS-conversion term, as was noted earlier in the discussion of Fig. 13. Previous photon PDFs were at best accurate to leading order, and hence had much larger scale uncertainties than LUXqed. Figure 19 shows the γγ luminosity compared with the gg and total qq luminosities for √ s = 13 TeV and 100 TeV, where the luminosity dL ij /d ln m 2 for partons i and j in pp collisions is defined by For the i q i q i luminosity, we have included a factor of two in the sum, since either quarks or antiquarks can come from each beam. The γγ luminosity is about three orders of magnitude smaller than the gg and qq luminosities over a wide range of masses. The impact of the γγ luminosity is however enhanced in processes with leptons and electroweak gauge bosons, such as pp → W W , where the γγ process has a t-channel exchange diagram which is enhanced in some kinematic regions. Figure 20 shows the photon momentum fraction of the proton as a function of µ. The momentum fraction is ∼ 0.43% at µ = 100 GeV, and increases by about 0.1% for each factor of 10 increase in µ. Eventually, neglecting electroweak corrections, the photon momentum fraction saturates at a value that is independent of α and α s , however this occurs at trans-Planckian scales.

Conclusions
In this paper, we have provided a detailed derivation and explanation of the results in Ref. [22]. In addition, we have discussed several extensions of the results. We have computed the polarised photon PDFs and the photon TMDPDF using the same method, and given the corresponding formulae.
We have given an alternative derivation of the photon PDF directly from the operator definition of the photon parton density, and we have used this definition to compute the MS conversion term in the photon PDF to one higher order in α s and α than our previous result. Our formalism allows for the computation of the MS conversion term to yet higher orders using MS DIS coefficient function expansions in D dimensions. On the phenomenological side, we have used the higher-order α s correction to give a more detailed analysis of the theoretical error.
We have highlighted what is needed as experimental input for the proton structure functions in order to make phenomenological use of our new higher-order calculations, in particular the QED contributions. Specifically, our photon PDF formula requires experimental data on DIS structure functions without removing the QED corrections on the hadronic side. Two-photon exchange contributions, which couple the hadronic and leptonic sides, must still be removed. We also stress that, rather than the measurement of the proton form factor for the elastic contribution, it would be useful to have a measurement of the structure functions also below the hadronic inelastic threshold W 2 = (m p + m π ) 2 , in order to account for final states consisting of a proton accompanied by an arbitrary number of photons and electron-positron pairs.
The evolution of the photon PDF can also be used to obtain higher order P γx DGLAP splitting functions. We verified that the resulting αα s and α 2 P γx unpolarised splitting functions agree with recent computations [26,27]. From our results it is also possible to derive the αα s polarised splitting function, as well as the unpolarised splitting functions to one higher order. We have not, however, given explicit results for them.
We have given a detailed explanation of the data inputs for our photon PDF, and the resulting uncertainties. The detailed treatment of the higher-order contributions used to evaluate the PDF and estimate uncertainties is slightly different from the previous version [22], as explained in the text, cf. also Fig. 14. To distinguish this new set from our previous determination, we call it LUXqed17 plus PDF4LHC15 nnlo 100 (or LUXqed17 for brevity), and it is valid for scales µ ≥ 10 GeV.
The few percent precision that we have achieved for the photon PDF is more than adequate for the computation of photon induced corrections to Standard Model processes at present. It is possible and natural to adopt this method in the context of global PDF fits. 17 Since we have shown that a proton collider can be viewed as a broad beam photon collider, and that the broad-beam distribution can be computed with high precision, the question remains whether such precision can be fully exploited in experimental measurements. The first thing that comes to mind is the possibility to search for totally hadrophobic BSM particles, whose cross section could be computed with very high precision using the results presented here. It is tempting to accompany such searches or measurements with a veto on the event hadronic activity, in view of the large elastic component that we found. We must point out, however, that even for the elastic component one cannot guarantee the absence of hadronic activity. While this is certainly the case for lepton-proton collisions, one should remember that in the proton-proton case the colliding protons are likely to pass close, or even to cross each other, sometimes giving rise to the production of secondary hadrons [66]. The precision of the photon PDF input would then be spoiled by a less well known survival probability factor for the protons.
The discussion in this paper has been for the photon PDF of the proton, but the method can be applied to derive a photon PDF formula for any hadron. Eqs. (6.16a, 7.19) for the unpolarised PDF hold without change for any hadron. The corresponding polarised formulae hold only for hadrons of spin 1/2. Evaluating the PDF from the formula requires, of course, experimental data on the hadronic form factors and structure functions. For the neutron, the high Q 2 data is available in terms of quark and gluon PDFs, but low-Q 2 data is not as accurate as for the proton. It is also important to keep in mind that at low Q 2 , a nuclear target such as a deuteron is not simply the sum of a neutron and proton. Finally, we remark that the methods developed here can also be extended to derive the lepton distributions in the proton.
Data files corresponding to the figures in this paper are available through Zenodo [67]. We will provide LHAPDF [20] files for the LUXqed17 plus PDF4LHC15 nnlo 100 and LUXqed17 plus PDF4LHC15 nnlo 30 sets through LHAPDF.

A QED corrections
In this Appendix, we summarise some known results on the renormalisation of the electromagnetic current [68][69][70] which will be useful for the subsequent discussion. The formulae are given for the case of QED for simplicity, but they are trivially generalised to the case of QCD with electromagnetic interactions.
The QED Lagrangian, including the gauge-fixing term is (A.1) Figure 21. Penguin graph which renormalises the electromagnetic current.
Bare quantities are labelled with a subscript 0. The renormalisation constants are defined through where S 2 = e γ E /4π, and the QED Ward identity implies that The gauge-fixing parameter ξ is not renormalised, by BRST invariance. The QED β function is given in terms of the anomalous dimension γ A of the photon field, where e(µ) is the MS coupling, and The Noether current is The usual (incorrect) textbook statement is that the electromagnetic current is a conserved current, and not renormalised, so that j µ N has finite matrix elements. This is false, because of the well-known penguin graph Fig. 21. As shown in Ref. [68][69][70], the renormalised current in the MS scheme is which satisfies the renormalisation group equation where the Gupta-Bleuler condition ∂ · A = 0 has been used for physical matrix elements. There is a unique local electromagnetic current Figure 22. (a) Graphs that are one-particle-irreducible in the current channel for insertion of a current vertex in a Green function or matrix element; the standard non-renormalisation argument applies only to these. (b) These graphs, one-particle-reducible in the current channel, also contribute to matrix elements of the current and to its renormalisation. The two subgraphs that are cross hatched are irreducible in the photon line, while the other subgraph gives the full propagator corrections to the photon propagator. (c) Counterterm to (b). The filled square corresponds to an operator proportional to ∂ ν F νµ .
which satisfies the usual QED Ward identities, and is not renormalised, Here Π(q 2 , µ) is the vacuum polarisation, 18 defined so that the renormalised photon propagator is In physical matrix elements, the equation of motion for the photon field gives iΓ (2) µν (q 2 , µ) = P µν q 2 1 − Π(q 2 , µ) − and Eq. (A.4,A.14) imply that e ph (q 2 ) is independent of µ. e ph (q 2 = 0) is the usual electron charge used in classical physics. e ph (q 2 ) depends on the exact vacuum polarisation function, and is non-perturbatively related to e(µ) in the presence of strong interactions, unless q 2 is at large Euclidean values, where perturbation theory holds.
We can express an arbitrary Green function or matrix element of j µ in terms of G µ 1γI , which is the 1-photon-irreducible part in the current channel. The total Green function is from Graphs Fig. 22(a), (b) and (c), respectively, giving 16) In physical matrix elements, the last term vanishes by current conservation, so Since j µ is µ-independent, so is G µ j and hence also G µ 1γI . The 1γI graph is the same whether one uses the current j µ N or j µ MS or j µ , since the three operators only differ in their 1γI part.
The one-photon-irreducible part is given by matrix elements of the non-local current We can now relate physical matrix elements of the current j µ defined in Eq. (A.9) and the MS current to one-photon-irreducible matrix elements, using Eq. (2.3). X|j µ |Y is independent of µ, as is X|j µ |Y 1γI , but not X|j µ MS |Y . Similar expressions hold for matrix elements with multiple insertions of the electromagnetic current.
In scattering processes with external (i.e. on-shell) photons, the LSZ reduction formula says that the S-matrix is given by where the wavefunction factor is so that Xγ|S|Y = e ph (0) X|j µ |Y 1γI , (A. 22) and external photons couple with strength e ph (0). Internal photon lines still have the coupling e(µ) (Sµ) . It is convenient in our discussion of the photon PDF to use one-photon-irreducible matrix elements, even though it is the matrix element of a non-local current Eq. (A.18). One can convert expressions to those of the local current Eq. (A.9) using Eq. (A.19).
The proton hadronic tensor is defined in Eq. (2.4). We will define three versions of this equation, W µν which uses the current j µ in Eq. (A.9), W µν MS , which uses the MS current Eq. (A.7), and W µν 1γI , which uses the non-local current Eq. (A.18), or equivalently, is given diagrammatically in terms of one-photon irreducible graphs. The relation between these is from Eq. (A. 19). All three hadronic tensors are renormalised using the MS scheme. The difference is which current is used in the hadronic matrix element in Eq. (2.4). One has to be careful about the definition of W µν and the structure functions in the presence of electromagnetic corrections. W µν defined using the µ-independent current Eq. (A.9), and including all electromagnetic corrections to the matrix element is µ-independent, so the structure functions in the decomposition Eq. (2.6) are also µindependent, and depend only on x bj and Q 2 . Note that if one had instead used the MS current, the structure functions would depend on µ as well as x bj and Q 2 .
The one-photon-irreducible F 2 structure function is denoted F 1γI 2 , and is related to F 2 by Similarly, the MS structure function is and depends on µ. Similar results hold for the other structure functions. The structure function F 2 (and similarly for the others) can be computed using the operator product expansion (OPE) as the convolution of short distance coefficient functions C 2,i (x, Q 2 , µ) and PDFs f i (x, µ). Graphically, it is easiest to compute C 1γI 2,i (x, Q 2 , µ), and then obtain C 2,i (x, Q 2 , µ) from Eq. (A.24).

B Kinematics
In this section, we summarise the kinematics for the scattering process k + p → k + X, where p is the incoming proton momentum with p 2 = m 2 p , k is an incoming massless particle momentum with k 2 = 0, and k is an outgoing particle momentum with (k ) 2 = M 2 . In deep inelastic scattering (DIS), M = 0, but we will need to consider cases where M = 0 as well. The momentum transfer q is defined as q = k − k , so that the final hadronic state has invariant mass m 2 X = (p + q) 2 ≥ m 2 p . We adopt light-cone coordinates, and introduce two null vectors n and n which define the collision axis, with n along the k direction, such that and n 2 =n 2 = 0 and n ·n = 2. Then p + = n · p, p − = n · p, etc. We adopt the convention and use the standard definitions We also introduce the variable The d 4 q integration measure in light-cone coordinates is Eq. (B.5) uses q + , q − and q 2 ⊥ as the independent variables. In a fixed Lorentz frame, with p ± given, q ± and q 2 ⊥ can be used to determine Q 2 , x bj and χ using Eqs. (B.3) and Eq. (B.4). We can express q ± , q 2 ⊥ as a function of χ, x bj and Q 2 by inverting these equations, We thus have a one-to-one mapping of the variables q 2 ⊥ , q + , q − into Q 2 , x bj and χ. However, the inverse mapping yields a valid kinematic point if and only if q 2 ⊥ is non-negative. Therefore we must add the constraint The relevant Jacobian is easily computed to be Q 2 /x 2 bj , so that the phase space in terms of the new variables becomes In computing the total scattering cross section, the integral over q is restricted by θ and δ functions (see Eq. (3.5)) which ensure that the final state has positive energy and the correct invariant mass. The energy theta functions are easily evaluated if we choose n andn such that p + = p − = m p , i.e. we work in the proton rest frame.
Since 2E ≥ 2E ≥ E + p ≥ 0, the factor in square brackets is positive, and We have 14) The invariant mass inequality m 2 X = (p + q) 2 ≥ m 2 p gives  Figure 23. Graph for the process γ + p → S + X at lowest order in α.
γ S q q Figure 24. Lowest order graph for γ + q → S + q.
The cross section coefficient has again been called σ 0 , so the formulae can be easily compared with the l → L case. The γ + p → S + X cross section, from Fig. 23 including the QED vacuum polarisation bubbles and corrections to the hadronic tensor is From Eq. (C.5) and Eq. (C.6), we obtain the same photon PDF as before, Eq. (3.26). Even though the non-logarithmic terms in Eqs. (C.5) and Eq. (C.6) are different from those in Eqs. (3.18) and (3.23), the difference of the two remains the same, and leads to the same MS-conversion term. This shows explicitly that, as expected, our derivation leads to a process-independent result for the MS photon PDF. The structure of radiative corrections for γ + p → S + X is more complicated than for l + p → L + X, so it is more difficult to extend the γ + p → S + X result to higher orders. Nevertheless, it must continue to give the same result for the photon PDF, even at higher orders. An alternate derivation using PDF operators that does not rely on any BSM process is given in Sec. 6. We have also checked the derivation of the polarised photon PDF using the γγ → S probe.
E Low Q 2 behaviour of F 2 and F L W µν is the discontinuity of a forward amplitude in W 2 = (p + q) 2 , and should be analytic in Q 2 and W 2 for W 2 away from the thresholds at (m p + nm π ) 2 , and for Q 2 < (2m π ) 2 . In particular, it should be analytic as Q 2 → 0 at fixed W 2 away from thresholds. This implies that the coefficients of its independent tensor structures should be analytic. Looking at the tensor structure q µ p ν in Eq. (2.6) we immediately see that F 2 must vanish as Q 2 .
Considering instead the tensor structure q µ q ν one finds that F 1 − F 2 /(2x) should vanish at least as Q 2 for small Q 2 . At small Q 2 and fixed W 2 , x behaves as Q 2 , so that vanishes at least as Q 4 as Q 2 → 0. where the proportionality constant Φ depends on the convention used for the incident photon flux for an off-shell (unphysical) photon with momentum Q 2 . In terms of structure functions, Pick a frame where p = (m p , 0, 0, 0) and the photon has momentum q = (q 0 , 0, 0, q 3 ) with longitudinal polarisation L = 1 Q (q 3 , 0, 0, q 0 ) .
The order coefficient functions are [38] a (1,0) 2,i∈{q} = e 2 i C F 2,q has a coefficient 9 + 3ζ(2)/4, rather than 9 + 3ζ(2)/2 given in Ref. [38]. With this modification, the Adler sum rule ∆i∈{q,g} coefficients are easily obtained from the above using the replacement rules Eqs. (G.6).  Figure 25. Left: the up-quark distribution at x = 0.1 as a function of the factorisation scale µ in the PDF4LHC15 nnlo 100 set and in its three underlying input PDF sets, CT14nnlo, MMHT2014nnlo68cl and NNPDF30 nnlo as 0118. One sees anomalous behaviour for µ < 1.295 GeV in PDF4LHC15 nnlo 100 and CT14nnlo, associated with LHAPDF's extrapolation of the CT14nnlo below its range of validity. Right: similarly, but for the b-quark distribution (multiplied by 10 3 ). It illustrates the different locations of the b-quark thresholds in the various sets. It also shows that only two of the sets display the expected discontinuous threshold at the b-quark mass.

H Issues at low and moderate µ 2 in PDF4LHC15
Our base PDF set, PDF4LHC15 nnlo 100 [51], is a combination of three underlying sets, CT14nnlo [52], MMHT2014nnlo68cl [53] and NNPDF30 nnlo as 0118 [54]. This combination is intended for LHC applications, which mostly involve high µ 2 values (e.g. (10 GeV) 2 upwards). In our work here, we need access to the PDF at low µ 2 values in order to then evolve it upwards with supplemental QED corrections. 19 In doing so, we have encountered some issues, which we document here.
A first point is that the set is quoted, within LHAPDF, as being valid from µ = 1 GeV. However if one uses it at this scale, one encounters unexpected behaviour such as a momentum sum of 0.94 rather than 1 and inconsistency with other µ values if one evolves up with an independent DGLAP code. The origin of the problem turns out to be trivial, namely that one of the input PDFs, CT14nnlo, is valid only from µ = µ 0 = 1.295 GeV (while the other two are valid from 1 GeV). When CT14nnlo is used below its starting scale, LHAPDF appears to extrapolate it such that f i (x, µ 2 ) ∼ µ 2 /µ 2 0 f i (x, µ 2 0 ) as µ 2 → 0. This is illustrated in Fig. 25 (left), which shows the up quark PDF versus µ. One sees the CT14 curve dropping rapidly for µ 2 < µ 2 0 , whereas the other input PDFs vary much more 19 With many standard methods for numerical DGLAP evolution, upwards evolution is stable, while downwards evolution is often less so. In particular parton distributions from LHAPDF tend to have small irregularities (typically below a part in 10 3 ) that make this especially problematic. One could imagine developing methods to make the downwards evolution more stable, however we have not investigated them.
slowly down to µ 2 = 1 GeV 2 . Since the PDF4LHC15 combination effectively averages the central values of the three sets, it inherits part of the extrapolated CT14 behaviour. This problem is common to all light flavours. It would be trivial to fix, by increasing the lower µ 2 limit of the PDF set to the largest of the lower limits of the underlying sets. A second issue concerns flavour thresholds. The three underlying sets do not share the same thresholds. For example CT14nnlo and MMHT2014nnlo68cl have their b-quark threshold at 4.75 GeV, while NNPDF30 nnlo as 0118 has its threshold at 4.18 GeV, as can be seen in Fig. 25 (right). Consequently, no single evolution that starts below the (highest) b threshold can reproduce the results of the combined set across all µ 2 values. 20 This suggests that the lowest scale from which one may start the evolution is just above the highest of the different b thresholds. However yet another issue arises. In Fig. 25 one sees that the MMHT and NNPDF b-distributions have a discontinuity at the respective b masses, associated with a second order threshold term in the evolution [75,76]. The CT14 curve is instead continuous there. For µ > 6 GeV it clearly approaches the MMHT2014 result, which suggests that the underlying evolution contains the correct mass threshold. It seems likely, therefore, that the issue is related to the way that the CT14nnlo set is included in LHAPDF. 21 This issue seems to be reflected also in the PDF4LHC15 set.
Overall, therefore, we can only use PDF4LHC15 nnlo 100 as a starting point for our evolution from a scale of about µ = 6 GeV upwards.