The photon PDF within the CT18 global analysis

Building upon the most recent CT18 global fit, we present a new calculation of the photon content of the proton based on an application of the LUX formalism. In this work, we explore two principal variations of the LUX ansatz. In one approach, which we designate"CT18lux,"the photon PDF is calculated directly using the LUX formula for all scales, $\mu$. In an alternative realization,"CT18qed,"we instead initialize the photon PDF in terms of the LUX formulation at a lower scale, $\mu\! \sim\! \mu_0$, and evolve to higher scales with a combined QED+QCD kernel at $\mathcal{O}(\alpha),~\mathcal{O}(\alpha\alpha_s)$ and $\mathcal{O}(\alpha^2)$. While we find these two approaches generally agree, especially at intermediate $x$ ($10^{-3}\lesssim x\lesssim0.3$), we discuss some moderate discrepancies that can occur toward the end-point regions at very high or low $x$. We also study effects that follow from variations of the inputs to the LUX calculation originating outside the pure deeply-inelastic scattering (DIS) region, including from elastic form factors and other contributions to the photon PDF. Finally, we investigate the phenomenological implications of these photon PDFs for the LHC, including high-mass Drell-Yan, vector-boson pair, top-quark pair, and Higgs associated with vector-boson production.

With the steady accumulation of copious experimental data at the LHC, we have entered into a high-precision era for hadron-collider physics. In parallel, theoretical computations of higher-order QCD corrections to standard LHC processes have reached to next-to-nextto-leading order (NNLO) in general, and, in some instances, even next-to-NNLO (N 3 LO) accuracy has now become available (see Ref. [1] for an overview). At this level of precision, electroweak (EW) corrections begin to have an observable impact, as α QED ∼ α 2 S . To that end, next-to-leading order (NLO) EW corrections to hard-scattering matrix elements have been computed for many LHC processes of interest, and the automation of the NLO EW corrections has also been achieved in recent years [2,3]. To perform consistent higher-order calculations with EW corrections included in the initial state of parton-scattering processes at the LHC, it is necessary to employ a set of parton distribution functions (PDFs) in which the photon appears as an active, partonic constituent of the proton.
The first such PDF set to include the photon as a parton of the proton was the 2004 release of the MRST group, the MRST2004QED PDFs [4]. The MRST group used a parametrization for the photon PDF based on radiation off of "primordial" up and down quarks, with the photon radiation cut off at low scales governed by constituent-or current-quark masses. Another approach to include the photon PDF is to constrain it in an analogous way to other partons by fitting available Deep Inelastic Scattering (DIS) and Drell-Yan data [5], an approach first developed by the NNPDF Collaboration and released in the NNPDF2.3QED [5] and NNPDF3.0QED [6] PDFs. A similar determination is also performed by the xFitter group [7]. The constraints on the photon PDFs from the data were rather weak, due to the small size of the photon-initiated contributions, so that the NNPDF photon PDF was consistent with zero at the initial scale of √ 2 GeV, with large photon PDF uncertainty at high x.
Contemporaneously, the CT14QED PDFs were constructed by implementing Quantum Electrodynamics (QED)-informed evolution at leading order (LO) along with QCD evolution at NLO within the CTEQ-TEA (CT) global analysis package [8]. The inelastic contribution to the photon PDF was described by a two-parameter ansatz, coming from radiation off the valence quarks, and based on the CT14 NLO PDFs. The inelastic photon PDFs were specified in terms of the inelastic momentum fraction carried by the photon, at the initial scale µ 0 , and they were constrained by comparing with ZEUS data [9] on the production of isolated photons in deeply-inelastic scattering (DIS), ep → eγ + X. The advantage of using this process is that the initial-state photon contributions are at leading order in the perturbation expansion. In contrast, the initial-state photon contribution to Drell-Yan or W and Z production is suppressed by factors of (α/α s ) relative to the leading quark-antiquark production. As discussed by Martin and Ryskin [10], the photon PDF also receives a large elastic contribution in which the proton remains intact, in addition to the inelastic contribution in which the proton breaks into a multihadron final state. 1 In addition to the default CT14QED set, CT therefore also released the CT14QEDinc PDFs, in which the inclusive photon PDF at the scale µ 0 is defined by a sum of an inelastic photon PDF and a corresponding elastic photon distribution obtained from the Equivalent Photon Approximation (EPA) [11]; in this latter formulation, the EPA allows the elastic photon contribution to be determined through the nucleon's elastic form factors, which can themselves be directly measured via elastic scattering experiments. Neither MRST nor NNPDF directly addressed these separated contributions to the photon PDF, although the NNPDF photon can be assumed to be inclusive, containing both inelastic and elastic components, since it was constrained using inclusive Drell-Yan and vector boson data. In all cases, the available data were unable to constrain the photon PDF to a high degree of accuracy.
In order to overcome these deficiencies, more accurate determinations of the photon distribution are necessary. The key to higher accuracy for photon PDFs is again the abovementioned EPA [11], in this case, extended to include inelastic contributions evaluated in terms of inelastic structure functions [12][13][14]. This structure function approach can then be employed to predict photon-initiated production processes of high-mass particles at hadron colliders [15][16][17]. This idea was revived by the LUXqed group [18,19], who demonstrated it in a rigorous framework based upon collinear factorization extending beyond leading or-der via the inclusion of the proper MS matching term; this work additionally released the first public photon PDF based upon this theoretical approach. As the elastic and inelastic proton structure functions have been determined experimentally to high precision, the photon PDFs can be constrained at the level of 1−2% over a wide range of the momentum fraction x. Furthermore, based on the anomalous dimensions in light-cone gauge, 2 the QED splitting kernels in the DGLAP [21][22][23][24] evolution equations have now been calculated up to O(αα S ) [25] and O(α 2 ) [26], whose effects are important to the determination of precise photon PDFs at a large energy scale µ as a function of x. Subsequently, the NNPDF group adopted the LUX formalism and introduced a photon PDF in a global PDF fit, named NNPDF3.1luxQED PDFs [27]. Likewise, the MMHT2015qed PDF set was released a number of years ago. These PDFs were generated in a global fit by adopting the LUX formalism at a low starting scale, µ 0 = 1 GeV, for the photon PDF, and evolving to higher scales using QED-corrected DGLAP evolution equations [28]. These features were inherited in the more recent MSHT20qed [29] analysis, which was based on the latest release of the MSHT20 global fit [30].
In this paper, we follow the original CT14QED strategy of separating the photon PDF into its respective elastic and inelastic components. The photon PDFs are generated based on two different approaches of applying the LUX formalism within the framework of the CT18 NNLO global analysis [31], where the NNLO QCD kernels have been used in the evolution of partons. In the first approach, CT18lux, the photon PDF is calculated directly using the LUX formula at any scale µ. In the second approach, CT18qed, we instead initialize the photon PDF in terms of the LUX formulation at a lower scale, µ ∼ µ 0 , and evolve to higher scales with a combined QED and QCD kernels at O(α), O(αα s ) and O(α 2 ). For convenience, we shall refer to the former as the LUX formalism approach, and the latter as the DGLAP evolution approach. While the PDFs generated by both approaches generally agree, particularly at the intermediate x region (10 −3 x 0.3), they differ in the low x region where the inelastic photon dominates and in the large x region where the contributions from the elastic component of the photon PDF becomes important. Hence, in this work, we have also explored the impact on the photon PDF from various recent updates of the elastic component, which represents the photon contribution from elastic photonproton scattering processes. As we shall discuss in greater detail below, implementation of the LUX formalism, which involves integrations of the proton's unpolarized electromagnetic structure functions, F 2,L , over broad Q 2 , can be sensitive to higher-twist (i.e., twist-4) and other nonperturbative QCD contributions. These effects are unsuppressed at low Q 2 and must be explicitly modeled; this is necessary both for theoretical accuracy as well as uncertainty quantification, for which an estimate of the possible model and parametric dependence is important. We point out that these considerations contrast with the situation in a typical NNLO PDF global analysis like CT18, in which only leading-twist (i.e., twist-2) dynamics are admitted into the relevant calculations. This is achieved by modeling only the twist-2 PDFs at the boundary of QCD evolution, µ = µ 0 , and constraining the resulting parametrization through an admixture of perturbative QCD parton-level cross sections and hadronic data at sufficiently high Q 2 and W 2 (for DIS) to ensure that contributions from sub-leading twist are safely, kinematically suppressed.
Finally, owing to its importance to the LHC phenomenology discussed at the end of this article, in this work we concentrate on the photon PDF of the proton. We note that it is also technically feasible to carry out analogous studies for the photon content of the neutron, as considered in the MMHT2015qed [28] and MSHT20qed [29] analyses. To treat such scenarios in full generality, it is necessary to consider the possibility of explicit chargesymmetry breaking at parton-level, and we defer such issues to later work.
The organization of this paper is as follows. In Sec. II, we present both the LUX formalism and DGLAP evolution approaches to generate photon PDFs. In Sec. III, we discuss the CT18lux photon PDF based on the LUX formalism. Various sources of the photon PDF uncertainty are also discussed. In Sec. IV, we present the result of the DGLAP-driven CT18qed, and compare various photon PDF sets with different choices of the input scale, µ 0 , where the photon PDF is provided by the LUX master formula. The main difference between the CT18qed and CT18lux photon PDFs will also be explored. In Sec. V, we investigate the phenomenological implications of these photon PDFs at the LHC, including high-mass Drell-Yan, vector-boson pair, top-quark pair, and Higgs associated with vectorboson production. In Sec. VI, we discuss our findings and give conclusions. Following the main body of the paper, we defer a number of technical details to a set of appendices. The separation of the photon PDF into elastic and inelastic components is discussed in App. A. In App. B, we detail the physical factorization and MS conversion terms that appear in the LUX formalism.

II. THE LUX FORMALISM VERSUS DGLAP EVOLUTION
As reviewed in Sec. I, the CT14QEDinc photon PDF [8] was comprised of two distinct sub-components. At the initial scale, µ 0 , CT14QEDinc is given by a sum of inelastic (γ inel , i.e., CT14QED) and elastic (γ el ) pieces. The CT14QEDinc photon PDF at any higher energy scale, µ > µ 0 , is obtained by solving the QED-corrected DGLAP evolution equation, The CT14QED photon PDFs were parametrized (by a two-parameter ansatz) and specified in terms of the inelastic momentum fraction carried by the photon at the initial scale µ 0 , which was constrained by comparing with ZEUS data [9] on the production of isolated photons in DIS, ep → eγ + X. The ZEUS experiment imposed a cut requiring at least one reconstructed track, well-separated from the lepton, which removed the elastic component of the photon. Although it is possible that some inelastic photon contribution would also fail this cut, this is expected to be small relative to the experimental uncertainties of the experiment [8]. The elastic component was parametrized by the Equivalent Photon Approximation [11], which involves an integration over the proton electromagnetic form factors. As implied by the Equivalent Photon Approximation, an alternative access to the photon PDF is through its relationship to electromagnetic scattering off the proton [32]. Early investigations on the connection between the photon PDF and the proton structure functions, F 2 (x, Q 2 ) and F L (x, Q 2 ), were presented in Refs. [12][13][14][15]. This idea of obtaining the photon PDF by viewing the ep → e + X scattering process as an electron scattering off the photon parton in the proton was formalized systematically within perturbative QCD/QED by the LUX group in Refs. [18,19]. In such a way, the photon PDF is fully determined by the structure functions, without the need of introducing a non-perturbative parameterization at ] R e s o n a n c e r e g io n ( C L A S / C B ) Low Q 2 continuum region (HERMES GD-11P) High Q 2 continuum region (pQCD CT18NNLO) an input scale, µ 0 . The master formula to determine the LUX photon PDF is 3 which includes all the αL(α s L) n , α(α s L) n , and α 2 L(α s L) n terms, with L = ln Q 2 /m 2 p . In this equation, p γq (z) ≡ [1 + (1 − z) 2 ]/z is the leading order DGLAP splitting kernel, and the O(α 2 , αα s ) term does not contain any logarithmic enhancement for large L. As explained in App. B, the LUX photon consists of two components. The first term inside square brackets corresponds to the physical factorization contribution, γ PF (x, µ 2 ), while the second term involving only F 2 is the MS conversion piece, γ con (x, µ 2 ), after a divergence is cancelled with the corresponding counterterms [19]. To perform the integration in Eq. (2), it is necessary to know the structure functions over the full (x, Q 2 ) plane. Inside the high-Q 2 and high-W 2 region, Q 2 > Q 2 PDF = 9 GeV 2 and W 2 > W 2 high = 4 GeV 2 , over which perturbative QCD is reliably applicable, the structure functions can be calculated directly from quark and gluon PDFs. In the low-W 2 region, W 2 < W 2 low = 3 GeV 2 , the original LUX methodology adopts the structure functions directly from phenomenological fits of CLAS 4 [33] or Christy and Bosted (CB) [34]. In the continuum region, Q 2 < Q 2 PDF and W 2 > W 2 high , the GD11-P fit of HERMES Collaboration [35] based on the ALLM funcion form [36] was adopted. A smooth and continuous transition from the W 2 low to W 2 high (the green band in 3 The MS running coupling α(µ 2 ) is related with the physical coupling α ph (q 2 ) as where q 2 = −Q 2 corresponds to the spacelike region, and Π(q 2 , µ 2 ) is the vacuum polarization. In the large momentum limit, |q 2 | m 2 q or m 2 , where m q (m ) is the masses of light quarks (leptons), We have the freedom to choose the renormalization scale as µ 2 = |q 2 | to make Π(q 2 , µ 2 ) = 0 and, therefore, α(µ 2 ) = α ph (−Q 2 ). 4 The CLAS fit is bounded by a threshold W 2 > (m p + m π ) 2 for nucleon resonance production, illustrated as the lower edge in Fig. 1.
In subsequent analyses, both the NNPDF [27] and MMHT [28] groups released photon PDFs based on the incorporation of the LUX formalism into their respective frameworks, making a number of different technical choices regarding the implementation of the LUX approach. Similar to LUXqed(17) 5 [18,19], NNPDF3.1luxQED [27] initializes the photon PDF with the LUX master formula, Eq. (2), at a high scale, µ 0 = 100 GeV, which falls within the high-Q 2 continuum region shown in Fig. 1. In this approach, the photon PDF is mainly determined by the fitted quark and gluon PDFs with a scale dependence specified by DGLAP evolution. As a consequence, the PDFs evolve bidirectionally in µ 2 .
Representing an alternative general scheme, which we collectively designate the "DGLAP approach" for the purposes of this article, the MMHT2015qed study [28] instead initialized the photon PDF at a low scale, µ 0 = 1 GeV, with an important modification from the default LUX setup. In the LUX formalism embodied by Eq. (2), to determine the photon PDF xγ(x, µ 2 0 ), the upper integration limit µ 2 0 /(1 − z) can become significantly larger than µ 2 0 when z is large. Therefore, MMHT breaks the integration over Q 2 into two parts, as where the dots above represent the expression inside the square brackets appearing in Eq. (2). In the second integration range of this quantity, the F L term is neglected because it is relatively suppressed by one additional order higher in α S compared to F 2 , i.e., F L ∼ O(α S ). Also, given their slow scale dependence, F 2 (Q 2 ) and α(Q 2 ) are approximately stationary, i.e., an observation which permits the sum over the second integration region to be performed analytically: very much like the modified conversion term in Eq. (B6). 5 The LUXqed17 [19] slightly differs from the original LUXqed [18] in the photon PDF calculation and the error estimation. A concise summary of these different methods is listed in Tab. I. In the pure LUX approach 6 , the photon PDF is fully determined by the structure functions which were either extracted from low-energy data, or calculated from the quark or gluon PDFs. In this way, the photon PDF is viewed essentially as an addition to the quark and gluon PDFs. The momentum sum rule will be violated by a small amount 7 , where the singlet PDF appearing above is defined as This violation is remedied in the LUXqed(17), NNPDF and MMHT approaches at the starting scale and maintained at other scales as well through DGLAP evolution. Strictly speaking, however, the momentum sum rule is actually very slightly violated even after being imposed at the starting scale, because the elastic photon satisfies a different evolution [19,28]. We expect this small violation to be only at the 0.01% level, which is fully negligible. Although the photon PDF in the LUX formalism is not obtained through DGLAP evolution, it still runs with the energy scale. It has been demonstrated in Ref. [19] that the LUX expression gives the correct DGLAP evolution kernel up to one order higher than the input coefficient functions. In another words, the one-loop order α s and α of Wilson coefficient functions 8 give the two-loop order αα s and α 2 p γi splitting kernels. As feedback, the splitting diagram (q → qγ) will affect the quark distributions through the QED real correction to the P qq function. Similar diagrams occur for the gluon at higher orders. As a consequence, all parton distributions should run simultaneously, a feature that is not captured by the pure LUX formalism.
In a realistic scattering of the form (p → γ)X → Y as shown in Fig. 2, we can factorize 6 We want to remind the reader that both LUXqed [18] and LUXqed17 [19] initialized the photon PDF in terms of the LUX formula at µ 0 = 10 GeV and evolved the QCD and QED DGLAP equation to obtain the PDFs at other scales. In this sense, they are based on the high-µ 0 DGLAP approach in our language. 7 Typically, this effect is of the order of a few per mille, depending on the quark and gluon PDFs used in calculating the DIS structure functions in the master formula. 8 Pay attention that the leading order of longitudinal structure function F L starts at one loop O(α s ) and the hadronic cross section as Starting at next-to-leading order in the QED coupling, quark-initiated processes begin to participate in the hard scattering as depicted in the middle diagram of Fig. 2. However, part of the inelastic photon originates within the (q → qγ) splitting, which contributes to the first diagram. The overlapping contribution should be subtracted in order to avoid double counting. After this subtraction as well as the cancellation of collinear singularities, the photon and quark (and gluon as well) distributions evolve according to the standard DGLAP equations. A similar situation will also occur when a hard interaction involves the partonic (q → qγ) subprocess in the initial state. In this scenario, inclusion of the QED corrections to the DGLAP evolution (i.e., to P qq in this case) ensures an explicit µ F cancellation in order to get the complete NLO QED precision. In the LUX formalism, this cancellation can only be achieved with proper MS conversion terms, which must be determined order-by-order, even though the difference may not be numerically sizable. In the end, the QED-evolved PDF set is recommended for realizing the corresponding precision in terms of higher-order QED corrections.
In order to accommodate this subtlety, we release two PDF sets in this analysis. In one set, the photon PDF is directly calculated with the LUX master formula, i.e., Eq. (2), at all scales. In this scenario, the quark and gluon PDFs are taken from the CT18 NNLO PDF fit [31] without modification. We designate this PDF set "CT18lux". In this sense, the momentum sum rule of CT18lux PDF will be weakly violated as the photon enters as an additional, small component. In the other set, the photon PDF is obtained through DGLAP evolution, and we call the result "CT18qed". In this set, we initialize the photon at a low starting scale, µ 0 , based on the LUX formalism, similarly to the approach taken by MMHT2015qed [28]. We emphasize that our "CT18lux" differs from the LUXqed(17) PDFs, which were obtained through a high-µ application of Eq. (2), and determined at other scales through DGLAP evolution. As pointed out above, the difference between CT18lux and CT18qed mainly comes from higher-order matchings, which we explore in detail in the next two sections.

III. THE CT18LUX PHOTON PDF
In this section, we first present a photon PDF determined through a more direct implementation of the LUX formalism [18,19] into the framework of CT18 NNLO global analysis [31] and obtain the CT18lux PDF. In the process, we will also illustrate the general use of the LUX master formula and the variety of physics inputs necessary to compute the photon PDF and its uncertainty. These developments will be instructive in Sec. IV, wherein we construct our main recommended photon PDF, CT18qed, which we base on a marriage of the LUX formalism and DGLAP evolution [21][22][23][24].

A. Numerical procedure
The photon PDF in the CT18lux is generated with the LUX master formula, Eq. (2), at all scales above the CT18 starting one µ 0 , according to the numerical prescription described below.
• Instead of taking PDF4LHC15 as input, we use the CT18 NNLO PDFs for the (anti)quark and gluon PDFs needed to calculate the structure functions F 2,L appearing in the LUX master formula, Eq. (2), in the high-Q 2 and high-W 2 region. The quark and gluon PDFs remain unchanged relevant to the default ones fitted in CT18 NNLO. • Similar to the CT18 NNLO PDFs, the CT18lux contains one central set and 58 Hessian error sets, with all generated by applying the LUX master formula to the 1+58 CT18 NNLO PDFs. The error sets quantify a part of the photon PDF uncertainty induced by the quark and gluon partons through perturbative DIS structure functions, F 2,L , in Eq. (2). • Because of this procedure, the proton momentum sum rule is slightly violated. The amount of violation is given by the additional fractional momentum of the proton carried by the photon. In Fig. 15, we show this momentum fraction as a function of the energy scale, µ. It is about 0.22% at µ = 1.3 GeV and grows to about 0.65% at µ = 1 TeV. • In addition to the 58 photon error PDFs encapsulating the quark and gluon PDF uncertainty noted above, a number of other dynamical effects at smaller Q contribute to the photon PDF calculation. We assess uncertainties associated with these in Sec. III B below, and ultimately include them into the final CT18lux photon PDF uncertainty.

B. The CT18lux photon PDF and its uncertainty
Based on the formalism developed in Sec. II, it is clear that xγ(x, µ 2 ) in the master formula of Eq. (2) involves a series of contributions away from the kinematical region dominated by inelastic processes as shown in Fig. 1 -the "DIS" or "high-Q 2 continuum region" in which the most appropriate description is in terms of PDFs. These contributions enter via the direct evaluation of γ(x, µ 2 ) from the phenomenological F 2,L structure functions, and arise from several distinct dynamical sources and scattering processes in the electromagnetic interaction of the photon with the proton. In turn, these contributions generally have a number of potential variations in terms of their calculation and implementation in Eq. (2) which represent a significant source of uncertainty on xγ(x, µ 2 ). We note that, while many of the nonperturbative error sources explored here may have been considered or analyzed in some form in previous work(s), here we systematically revisit them, having considered an expanded set of possible variations and, where relevant, updating inputs into Eq. (2) with more recent parametrizations of structure functions and form factors that have emerged since, e.g., Ref. [18].

Construction of the photon PDF uncertainty
The most straightforward way to determine the full photon PDF uncertainty is to build it sequentially, by first calculating a central PDF through Eq. (2) as detailed above, and then computing eigensets associated with variations among the PDF parameters used to compute F 2,L . This encapsulates the photon PDF uncertainty from the "continuum QCD" or higher-Q 2 region. A realistic assessment of the uncertainties of the extracted photon PDF also depends, however, on the treatment of low-energy (or, "low-Q 2 ") contributions to F 2,L (through and including Q ∼ few GeV), especially at high x; we must therefore account for variations of these low-energy effects within the larger Hessian uncertainty on γ(x, µ 2 ) as well. In particular, we specify the uncertainty on γ(x, µ 2 ) through a number of additional error sets, leading to a collection of PDF sets in a fit with N PDF shape parameters for the parton distributions and n low-Q 2 separate low-Q 2 inputs to F 2,L . For the latter, we symmetrize the uncertainty for each of the low-Q 2 inputs, and rescale the error band to 90% CL, with the convention adopted in the CT PDF family. In such a way, the combined PDF uncertainty can be obtained through a Hessian approach An overall factor 1.645 is applied to get the 68% CL uncertainty, for the comparison with other PDFs throughout this work.

Contributions to the photon PDF and its uncertainty
Here, we describe a number of unique contributions or effects that enter the photon PDF master formula of Eq. (2) and are responsible for the ultimate uncertainty of the photon PDF.
The quark and gluon PDF uncertainty. The structure functions, F 2,L , in the high-Q 2 continuum region are determined through perturbative calculations. The uncertainty of high-Q 2 F 2,L due to the quark and gluon uncertainty will propagate to the inelastic photon component. We show the photon PDF uncertainty purely induced by the high-Q 2 quark and gluon PDFs in Fig. 3. We have presented the same error bands from LUXqed17 and MMHT2015qed as well for comparison. 9 We see that the overall size of the error band agrees very well among these three groups, while MMHT2015qed gives a slightly larger band in the intermediate x region. Instead of directly calculating the photon PDF with the LUX master formula, MMHT2015qed evolves the photon together with the quark and gluon PDFs with the DGLAP equations, which includes the interplay between the photon and other partons, as well. This treatment is similar to that used in CT18qed (1.3GeV), which will be discussed in the next section. Elastic form factors. Another necessary consideration is the interplay of parametric uncertainties in phenomenological fits or models of Sach's electromagnetic form factors of the proton, G E,M , and the ultimate photon PDF uncertainty. Historically, there have been diverse attempts to simulate these form factors in QCD-inspired models [37][38][39][40][41][42], compute them using lattice QCD or other theoretical methods (see, e.g., [43]), or to extract them phenomenologically based on empirical data, analogously to the quark and gluon PDFs themselves. Moreover, a number of simplified prescriptions have also been used for the Q 2 dependence of the electromagnetic form factors; most commonly, this has been a dipole expression based on the observation that the quark model and high-energy QCD prescribe an asymptotic ∼ 1/Q 4 dependence at large enough Q 2 values. At lower Q 2 , the observed form factor is known to deviate significantly from this behavior, and, for the sake of comparisons, we take a dipole ansatz baseline of where the value of the dipole mass parameter is Λ 2 dip = 0.71 GeV 2 and µ p = 2.793 [44] is the magnetic moment of the proton. In this work, we have primarily used the empirical data-driven approach, taking a series of phenomenological fits of the Sach's form factors in order to more fully determine the dependence of our photon PDF on choices for G E,M (Q 2 ).
LUXqed (17) [18,19] has explored variations based on the A1 fit [45] with all the world data involving electron scattering up to the year 2013. We show the A1 fits to the world polarized data and the unpolarized data in the upper pannel of Fig. 4. We notice that the polarized fit include two-photon-exchange (TPE) corrections, while the unpolarized one does not. Corresponding ratios to the commonly-used single-parameter dipole prediction are shown in Fig. 4, middle and lower panels; these ratios highlight the nontrivial Q 2 dependence of the realistic form factors, which can markedly differ from that of the simpler dipole ansatz. In addition, a more recent global fit of the unpolarized world cross-section data from Ye et al. in Ref. [46] similarly incorporated two-photon exchange corrections. We see that G E obtained in this fit agrees better with the A1 fit of polarized data, whereas G M agrees with the A1 unpolarized one. The elastic photon PDF as obtained with these various prescriptions is shown in the panels of Fig. 5. In particular, the corresponding electric and magnetic components of the elastic photon PDF are shown at µ = 100 GeV in the left panel of Fig. 5. For reference, the full calculation of xγ el (x) using the dipole parametrization of Eq. (11) is also shown. We see that the electric contribution dominates the elastic photon in the low-x region, while the magnetic one takes over at large x. This can be understood in terms of the asymptotic behavior of the elastic photon in Eq. (A1). In the small-x region, the Q 2 integration over the electric form factor behaves as ∼ln(1/x) while the analogous integration over the G M -dependent term approaches a constant. On the other hand, at large contribution. Meanwhile, in the right panel of Fig. 5 we plot the ratio at µ = 100 GeV of the elastic photon PDF computed with each form factor parametrization explored in this analysis with respect to that calculated using the A1 polarized fit. We see that the fit of Ye et al. [46] gives a larger uncertainty band than that of A1, due to its broader uncertainty in the low-Q 2 region. At the same time, the central elastic photon associated with the Ye et al. form factor fit agrees better with the A1 unpolarized calculation at large x. The dipole form factor result extends beyond the error bands for most x regions, corresponding to an overestimate of the elastic photon at small x and extremely large x, and an underestimate at moderately large x. Following LUXqed (17) [18,19], we retain the A1 polarized fit as our default choice, while the corresponding uncertainty is estimated by the larger deviation obtained either from the error band of the A1 polarized fit or from the central of A1 unpolarized fit. We also release an extra PDF set based on the Ye et al.'s form factors, but it is not included in the PDF uncertainty estimation.
Higher-twist effects. The structure functions are formally determined as Fourier transforms of hadronic matrix elements of quark-level electroweak current operators, in which the index i on the RHS can be expanded over a complete tensor basis and L i µν represents the set of unique Lorentz structures constructed from the momenta p and q in that basis. At lower energies, various soft quark-parton correlations within the target can be included in the matrix element of Eq. (12), which have the effect of increasing the twist (difference between operator spin and dimension) of the matrix element. For unpolarized matrix elements, the next nonzero term in the twist expansion beyond leading twist appears at twist-4, which is power-suppressed by ∼ 1/Q 2 . Computing or fitting these separately may be an involved task, but one quick alternative would be to use an existing determination such as that from the CJ15 NLO fit [47], which fitted an x-dependent higher-twist (HT) correction, where While there is an obvious model uncertainty associated with this choice of parametric form, there is already a parametric uncertainty that can be used to generate extreme scenarios for the HT correction to F 2 . The absolute higher-twist corrections to F 2 (x, Q), based on the CJ15 NLO fit [47], are shown as solid lines in left plot of Ye et al.
Ye et al.  [45] and Ye et al. [46]. We stress that the dipole calculation is provided only for the sake of reference and is not actively used in subsequent calculations of the photon PDF or related phenomenology, with the exception of the results shown in Fig. 5. Fig. 6. We see clearly that, as expected, the size of the HT correction is maximized at lower values of Q 2 , particularly for large x, but relatively suppressed as the Q 2 scale increases. The sharp dip around x ∼ 0.48 is due to a localized sign change in C HT (x) in the CJ15 NLO fit. Similarly, we expect a higher-twist contribution to F L , which is parameterized with one parameter as Recent studies on the DIS data suggest A HT = 5.5 ± 0.6 [48]. The MMHT group obtained  a smaller value as A HT = 4.3 GeV 2 [49]. Here we take the larger value from Ref. [50], with the impact shown as the dashed lines in Fig. 6, left plot. The corresponding corrections to the inelastic photon PDF arising from these HT corrections are represented by the blue curves in Fig. 8. The overall size of the uncertainty is generally small, a fact which can be understood as originating in the ∼ 1/Q 2 suppression of HT corrections, which here only enter the calculation in the high-Q 2 continuum region. Nevertheless, the HT uncertainty peaks at 1% for x ∼ 0.65, and remains a relevant consideration for precision in the very high-x regime.
Target-mass corrections. At leading twist, additional operators can be inserted into the matrix element associated with the hadronic tensor, W µν of Eq. (12), which, in turn, affect the low-Q behavior of the unfolded structure functions. These additional operators are insertions of covariant derivatives: which increase both the spin and dimension of the matrix element, thus leaving the twist unchanged. These corrections are referred to as "kinematical higher twist" corrections, and can be unfolded in the operator product expansion (OPE) of Georgi-Politzer [51], as done here. The net effect is target-mass dependent corrections (TMCs) that go as ∼ m 2 p /Q 2 and dominate mass-corrected structure functions at high x. There, in principle, can be substantial prescription dependence in the implementation of these corrections, leading to another potential uncertainty. A number of these structure-function level prescriptions have been computed and reviewed in Refs. [52,53]. For example, there can be substantial variation in F L , depending upon the specific prescription. The target mass corrections to F 2,L , defined as are shown in Fig. 6, right plot, calculated according to the standard OPE formalism implemented in APFEL [54]. We remind the reader that the TMCs can be either positive or negative, but we only show the relative absolute deviations here. Similar to the HT case, the several sharp dips are due to sign changes on the log-scale over which we plot. Owing to the kinematical suppression of target-mass effects, which result in a naïve rescaling of computed structure functions in the Georgi-Politzer formalism [51], TMCs introduce a kinematical dependence on the ratio x 2 m 2 p /Q 2 . In particular, at larger x and smaller Q 2 , we expect the most significant impact, especially for quantities that are steeply-falling functions at x → 1; this is true, for instance, of the longitudinal structure functions, F L , which can receive a sizable correction induced by rescaling, x → ξ. The TMCs to the inelastic photon PDF are shown as black lines in Fig. 8. We see the impact can dominate the uncertainty at very high x, e.g., for x > 0.6.
Scattering from nucleon resonances. At lower energies, Q 2 < Q 2 PDF or W 2 < W 2 high , scattering from nucleon resonances (e.g., the ∆, the Roper) dominates the γp cross section and therefore contribute to F 2 and F L (or F 1 equivalently). 10 These are typically described with a combination of Breit-Wigner parametric forms on a smooth background (for instance, informed by Regge Theory). Even more so than the elastic contributions, uncertainties from resonances are likely to be truly data-driven, with little variation from underlying theory or models. The structure functions F 1,2 in the resonance regions at several representative scales are shown in Fig. 7. Here we show the invariant squared mass of the final-state hadronic system, W 2 , which is directly used in the experimental measurements, such as CLAS data [33] or the Christy-Bosted phenomenological fit [34]. The conversion from W 2 to the Bjorken-x can be easily performed as 10 The relation among these three structure functions follows  The "CH" curves in Fig. 7 bridge the CLAS resonance to the HERMES continuum with a smooth transition: where a = 1, 2 (or L) and ρ is The transition points are W 2 low = 3 GeV 2 and W 2 high = 4 GeV 2 , respectively. The original Christy-Bosted fit was released in 2007 [34], which is denoted as "CB07" in Fig. 7. A recent update includes more data from proton and mostly nuclei cross sections [55], denoted as "CB21" in the figure. We show the corresponding inelastic photon of CB21 as the green curves in Fig. 8. We see the variation from our default choice, CLAS, is very mild, and only 1% around x ∼ 0.5. It only deviates significantly when x > 0.85, where the structure functions become unreliable, and so does inelastic photon.
The R L/T . In the resonance and low-Q 2 continuum region, the longitudinal structure function F L is modeled, similarly to LUXqed(17) [18,19], as where R L/T = σ L /σ T . The specific R L/T is provided by LUX group, who adopted the HER-MES convention [35] and the R 1998 fit provided by E143 Collaboration [56]. The uncertainty is assigned conservatively to be ±50%, with the corresponding inelastic photon shown as the red bands in Fig. 8. Around x ∼ 0.45, the uncertainty induced by R L/T can be as large as 2%, which dominates around this region. low-Q 2 HERMES to high-Q 2 pQCD continuum region is set as Q 2 PDF = 9 GeV 2 . LUX varied this scale down to Q 2 PDF = 5 GeV 2 in order to estimate the corresponding uncertainty due to this parametric choice. Different from the LUXqed(17) [18,19] and NNPDF3.1luxQED [27] calculations, the MMHT2015qed photon PDF was initialized at µ 0 = 1 GeV is therefore fully determined by the low-energy SFs, taken from the fits of HERMES [35] and CLAS [33] (or CB [34]). The corresponding uncertainty there is quantified by taking the uncertainty bands of the GDP-11 fit provided by the HERMES Collaboration [35]. At higher scales, µ, γ(x, µ) in MMHT2015qed is entirely determined by DGLAP evolution. For CT18lux, we choose to follow the LUX approach to quantify the matching-scale uncertainty, which is shown as the orange curve appearing in Fig. 8. We find the variation from Q 2 PDF = 9 GeV 2 down to 5 GeV 2 produces a very mild, subpercent effect which is maximally peaked at high x ∼ 0.6.
Missing higher-order (MHO) correction. As explained in App. B, the missing higher-order uncertainty is quantified by varying the separation scale M 2 (z) between the default choice µ 2 /(1 − z) and the alternative one µ 2 . The variation of the inelastic photon PDF from the CT18lux central set is shown as the purple line in Fig. 8. We find that, with the expected cancellation of the divergent log[1/(1 − z)] contribution captured by the modified MS-conversion term in Eq. (B6), we obtain only very mild variations in the resulting photon PDF at nominal x values. The only exception occurs at x > 0.85, where the absolute photon PDF becomes extremely small, as indicated in Fig. 27, such that the calculation becomes unreliable.
To summarize this discussion, in Fig. 9 we present the complete set of contributions to the CT18lux photon-PDF uncertainty, contributed from each of the various sources discussed above. Here, we have summed the elastic and inelastic components and show the total one.

IV. THE PHOTON PDF FROM A DGLAP-DRIVEN APPROACH: CT18QED
In this section, we present the photon PDF based on a DGLAP evolution approach, which consistently includes QED effects. We call the resulting photon PDF sets "CT18qed." After first reviewing our numerical procedure for evaluating PDFs in CT18qed (Sec. IV A) and discussing the CT18qed PDFs themselves and their uncertainties (Sec. IV B), we turn to several issues related to the determination of PDFs in the DGLAP-driven approach, including the enforcement of momentum conservation (Sec. IV C), subtleties comparing the CT18lux and CT18qed calculations at different perturbative orders in the combined QCD and QED expansion, and the impacts of QED evolution within the CT global analysis.

A. Numerical procedure
The CT18qed PDFs are constructed in such a way as to first separate the photon PDF into elastic and inelastic components, in analogy with the CT14QEDinc and CT14QED PDFs [8]. The elastic photon distribution is directly calculated by applying the LUX formula, cf. Eq. (A1), while the inelastic photon distribution, together with the (anti)quark and gluon partons of the proton, are predicted by applying the DGLAP equations to evolve PDFs from an initial scale, µ 0 ∼ [few GeV], to an arbitrary scale at higher energies, µ.
• After separating the photon PDF into elastic and inelastic components as noted above, the elastic photon distribution at any scale, µ, is directly evaluated by applying the appropriate elastic LUX formula, Eq. (A1). At µ 0 = 1.3 GeV, the elastic photon contributes to about xγ el (µ 2 0 ) = 0.15% momentum fraction of the proton. • At µ 0 = 1.3 GeV, we determine the inelastic photon distribution according to the LUX master formula, Eq. (2), excluding the elastic photon component in this stage of the calculation; here, the structure functions F 2,L in the DIS region are calculated directly from the CT18 NNLO PDFs [31]. At our starting scale, the inelastic photon contributes xγ inel (µ 2 0 ) = 0.066% to the total momentum of the proton. We also note that, since CT18 has one central set and 58 error sets, the resulting CT18qed PDFs also contain 1+58 PDF sets corresponding to the central prediction and additional Hessian error sets associated with the underlying uncertainty arising from the quark and gluon PDFs.
• At µ 0 = 1.3 GeV, the shape and normalization of the valence quark and gluon PDF are fixed to those obtained in CT18, except that we require the total momentum fraction carried by the quark sea to be 1 − . This ensures consistency with the momentum sum rule at µ 0 = 1.3 GeV. This procedure is implemented by adjusting the normalizations of sea-quark PDFs in the fashion typically used in CT global fits so as to respect total momentum conservation.
• We then evolve the quark, gluon and inelastic photon PDFs from the starting scale (µ 0 ) to an arbitrary scale µ (> µ 0 ) by applying NNLO QCD plus NLO QED DGLAP evolution equations implemented inside the APFEL package [54]. Consequently, a new set of quark and gluon PDFs, slightly different from those obtained in CT18, are generated with the photon PDF. Crucially, owing to the properties of the combined QCD+QED splitting kernel, the total momentum carried by quark and gluon degreesof-freedom and the inelastic photon will remain invariant, having been fixed to our initial choice of x(g + Σ + γ inel ) (µ 2 ) = 1 − xγ el (µ 2 0 ). Since the momentum fraction carried by the elastic photon will change only very slightly with increasing µ, decreasing by about δ xγ el (µ 2 ) = 0.015% at µ = 10 TeV, the CT18qed PDFs are guaranteed to satisfy the momentum sum rule to greater accuracy than in CT18lux.
• We also compare this scenario outlined above with a different choice for the starting scale of µ 0 = 3 GeV. In this case, the input quark and gluon PDFs at the scale µ 0 = 3 GeV are those predicted by the standard CT18 PDFs (with µ 0 = 1.3 GeV) after having been evolved from 1.3 GeV to 3 GeV. We note that in this case, the charm PDF is non-zero at the starting scale of µ 0 = 3 GeV. For reasons explained in the next subsection, we will from here on take this PDF set with µ 0 = 3 GeV as the default choice of CT18qed, and use the name CT18qed to refer to this set in particular, if no other specifications are given. We will use the name CT18qed1.3GeV, instead, to refer specifically to the PDF set with µ 0 = 1.3 GeV as the starting scale, when a distinction is required. • We stress that the CT18qed PDFs described above are based on a boostrap logic, with an underlying assumption that the effects of QED evolution can be regarded as a minor perturbation about the CT18 global minimum. We validate this assumption in the context of an actual QED refit. Various impacts of the inclusion of QED in a fit, including those upon the χ 2 and fitted PDFs, as well as the redistribution of momentum fractions due to the photonic radiation, are explored. Ultimately, we find these effects to be negligible, such that the CT18qed PDF sets above remain effectively unchanged.
B. The proton's photonic content in CT18qed CT18qed PDFs. As outlined above, it is possible to determine the photon PDF at any perturbative scale through combined QCD+QED evolution from a given initial boundary  [19], NNPDF3.1luxQED [27], and MMHT15qed [28].
condition at the starting scale µ 0 . Various published photon PDF sets were obtained by choosing different µ 0 and their associated boundary conditions. In Fig. 10, we compare the CT18qed (with µ 0 = 1.3 GeV and 3 GeV) and CT18lux photon PDFs with other existing photon PDFs at µ = 100 GeV. It shows that in the intermediate-x region, the CT18lux photon is between LUXqed17 (similar to NNPDF3.1luxQED) and MMHT2015qed, while CT18qed gives a smaller photon distribution than all other sets. In contrast, at extremely large x, the MMHT15qed and CT18qed calculations, both of which are based on a lowscale evolution approach, give relatively smaller photon PDFs than the others. Meanwhile, at small x, CT18qed gives a larger photon PDF than CT18lux. This can be qualitatively understood in terms of the approximation of solution to Eq. (1) as Recall that structure functions at the leading order, i.e., O(α 0 s ), are Hence, in this approximation, the DGLAP solution agrees well with the LUX prediction, when the LO structure functions are used and the MS conversion term is ignored in Eq.
(2). In Fig. 11, left plot, we directly compare the LO with the full NNLO calculation of F 2 employed in the LUX approach. It shows that at high scales µ, the ratio F LO 2 /F NLO 2 is larger than one in the small x region, and becomes much smaller than one in the large x region. This explains why the CT18qed photon PDF is larger than CT18lux at small x values and becomes smaller when x increases. Another noticeable feature of Fig. 10 is that CT18qed photon PDF drops very fast as x approaches to 1. In order to understand the reduction of CT18qed at large x values relative to CT18lux, we need to keep in mind that both the elastic and inelastic photon PDFs drop rapidly when x → 1. Given that the elastic components are the same for CT18lux and CT18qed and essentially scale-invariant (see Fig. 27), the main difference between the two comes from the inelastic contribution, especially from the MSconversion term. In CT18qed with µ 0 = 1.3 GeV, 11 the starting scale of µ 0 falls within the low-Q 2 continuum and resonance regions and therefore receives substantial non-perturbative contributions from the effects reviewed in Sec. III B 2. In the right plot of Fig. 11, we compare a few non-perturbative F 2 (solid lines) with perturbative ones (dashed lines). We see that the non-perturbative F 2 is significantly larger than the perturbative one. 12 Therefore, we expect the corresponding absolute value of low-Q 2 MS-conversion terms to be significantly larger than the one in the high-Q 2 pQCD region. Recalling the negative sign of the MS conversion term, this explains the significant reduction of CT18qed at large x values compared with CT18lux. The same scenario occurs with MMHT2015qed at large x, as shown in Fig. 10. Finally, MMHT2015qed gives a smaller photon PDF in the small x region as compared to CT18qed, which is due to the comparatively lesser value of its charge-weighted singlet distribution, Σ e , which is depicted in Fig. 12. The CT18qed PDF total uncertainty. The next question to ask is "What is the photon PDF total uncertainty in CT18qed?" Before detailing various potential sources of the photon PDF uncertainty, we first compare the photon PDF uncertainties as predicted by various global analysis groups. In Fig. 13, we show the self-normalized uncertainty bands to the corresponding central for each of the photon PDFs examined in this analysis.
The photon PDF uncertainty induced by quark and gluon degrees-of-freedom can be calculated, as described above, in the same way as the PDF uncertainty for the quark and gluon PDFs -by applying the master formula presented in Ref. [57] to the complete set of quark and gluon PDF error eigensets. In fact, the uncertainty induced by the DGLAPevolved quark and gluon PDFs has already been shown in Fig. 3, which roughly gives the same size of error bands as compared with CT18lux and other groups, like LUXqed17 and MMHT2015qed. Similarly to the investigation in Sec. III for CT18lux, we investigate all other low-energy sources of uncertainty, which we summarize and compare side-by-side with CT18lux in Fig. 14. Much as we observed for CT18lux, in the small-x region the CT18qed 11 It is also true for the µ 0 = 3 GeV in the extremely large x region of W 2 > W 2 high in Fig. 1. 12 We remind the reader that the CLAS fit is bounded by a threshold condition for nucleon resonance production, W 2 > (m p + m π ) 2 , which gives a boundary of x at low Q 2 . The F 2 will drop to zero beyond this point, which explains the sharp drop of the solid non-perturbative F 2 in Fig. 11. photon uncertainty is dominated by the uncertainties in the quark and gluon PDFs. In the large-x region, however, most of the low-energy error sources are also significant, inducing uncertainties that are comparable to those seen in CT18lux, with the exception of the M 2 (z) variation ascribed to missing higher-order (MHO) effects. For CT18qed, with µ 0 = 1.3 GeV, the MHO contribution becomes dominant at x 0.2, giving larger error bands than CT18lux in this region. However, due to the much smaller MHO contribution, the total uncertainty band of photon PDF in CT18qed with µ 0 = 3 GeV only becomes large at much larger x values (x 0.6). As discussed in App. B, the MHO uncertainty is estimated by shifting from the separation scale M 2 (z) = µ 2 /(1−z) to µ 2 . The corresponding MS conversion term should be changed as Eq. (B6), which is semi-analytically FIG. 14. Separated contributions to the full uncertainty of the photon PDF as obtained in CT18qed(1.3GeV) (upper) vs. the analogous result in CT18lux (lower), which we again display here for comparison. We note that the contributions discussed here arise from various low-energy sources as enumerated in Sect. III B 2.
integrated over the interval Q 2 ∈ [µ 2 , µ 2 /(1 − z)], using the stationary approximation to replace the Q-dependent F 2 (x/z, Q 2 ) with the fixed-µ 2 contribution F 2 (x/z, µ 2 ). In the pQCD DIS region, F 2 respects the stationary condition with a small logarithmic violation resulted from higher-order corrections. However, this condition is not respected that well in the non-perturbative region, when Q 2 < Q 2 PDF , show in Fig. 11. Hence, the MHO contribution in CT18qed with µ 0 = 1.3 GeV is much larger than that with µ 0 = 3 GeV, particularly in the large x region, as shown in Fig. 14. For this reason, we have chosen to take the PDF set with µ 0 = 3 GeV as the default CT18qed, and the one with µ 0 = 1.3 GeV as an alternative set, which is dubbed as CT18qed1.3GeV for distinction.

C. Momentum conservation and QCD+QED evolution in CT18qed
Photon PDF moments and the momentum sum rule. The proton's partonic constituents are expected to satisfy the momentum sum rule given by Eq. (6). CT18qed, like MMHT2015qed, determines an initial photon distribution at a low scale before consistently evolving all parton flavors under a combined QCD+QED kernel to higher scales in a fashion that closely preserves the total proton momentum. In Sec. IV A we outlined this procedure before showcasing the resulting PDFs in Sec. IV B. Here, we compare the first moments of the photon and other PDFs, and relate this to aspects of the QCD+QED evolution framework. , at a number of scales as obtained in CT18qed and CT18lux (leftmost columns) as well as in several other recent analyses. We remind the reader that µ min here is the lowest energy scale of the corresponding LHAPDF grids, which is same as the DGLAP initialization scales, µ 0 , of CT18lux, CT18qed1.3GeV and MMHT2015qed, but different from the µ 0 of CT18qed, LUXqed(17) and NNPDF3.1luxQED.
In addition, we compare properties of evolution in CT18qed with MMHT2015qed as well as the other fitting groups. As noted before, we enforce the momentum sum rule in CT18qed at µ 0 as making use of the very mild scale dependence of the first moment of the elastic photon PDF, xγ el (µ 2 ). In turn, the DGLAP evolution of the remaining components of Eq. (25) is such that the total inelastic momentum, x(Σ + g + γ inel ) (µ 2 ), remains fixed due to the momentum-conserving properties of the splitting functions: The scale dependence of the separate (in)elastic contributions of the photon momentum are plotted in the plot of Fig. 15. We remind the reader that the combination x(Σ + g + γ inel+el ) (µ 2 ) does not exactly conserve momentum due to the very minor scale variations in γ el (µ 2 ) depicted in the lower set of curves shown in Fig. 15, left plot. As a companion to these plots, in Table II, we display the first moments of the photon PDF as obtained by CT18qed and CT18lux at a number of relevant scales, and compare with the corresponding results published by MMHT2015qed, LUXqed17, and NNPDF3.1luxQED. For completeness, we also show the comparison of the total photon momentum fractions among different PDFs in the right plot of Fig. 15, with the specific numbers at a few typical scales listed in Table II. We see the overall size agrees very well among different photon PDF sets, both for the absolute values as well as for the uncertainties. MMHT2015qed gives 2∼3% larger than other groups, due to its larger elastic photon. (See Fig. 28 in Appendix B.) LUXqed17 gives very much the same as CT18lux, while CT18qed is slightly smaller.
One important feature appears for the inelastic photon that the low-µ 0 DGLAP approach in CT18qed yields a smaller photon at a low scale than LUX one in CT18lux, but gradually exceeds when energy increases up to certain scales.
Here, we note that CT18lux violates the momentum sum rule of Eq. (6) very weakly, as the photon is added on top of the existing quark and gluon distributions without making compensating adjustments to the latter. This small violation can be quantified by the  averaged momentum fraction carried by the photon, for which the scale dependence of the separated inelastic and elastic components are shown in Fig. 15. At the CT18lux starting scale, µ 0 = 1.3 GeV, these separate contributions are xγ el (µ 2 0 ) = 0.15% , xγ inel (µ 2 0 ) = 0.066% (28) which are 0.15% and 0.11%, respectively, when µ 0 = 3 GeV. We note that xγ inel (µ 2 ) receives a contribution from the q → qγ splitting, which grows logarithmically with µ 2 . In contrast, xγ el (µ 2 ) remains nearly constant, with a small decrease that varies as α(µ 2 0 )/α(µ 2 ). In order to examine the size of the impact resulting from enforcing the momentum sum rule, we compare our default treatment, Eq. (25), with an intermediate one, x(g + Σ + γ inel ) (µ 2 0 ) = 1, and a CT18lux-like constraint, x(g + Σ) (µ 2 0 ) = 1, in Fig. 16. Here we show the charge-weighted singlet, Σ e , on the left and the inelastic photon, γ inel , on the right, normalized to the CT18lux PDFs. Due to the reduction of the momentum fraction of the quark sea by the photon, the overall relative decrease in the charge-weighted singlet PDF is roughly 2.5 xγ (µ 2 0 ) ∼ 0.5%, especially in the small-x region, as shown in the left plot of Fig. 16. Reciprocally, the inelastic photon is reduced by roughly the same amount as we explicitly show in the right-hand plot of Fig. 16.
Comparison of CT18lux and CT18qed at a different order of QED evolution. Critical aspects of the scale dependence and agreement among photon PDF calculations follow from the implementation of perturbative QCD+QED evolution in CT18qed vs. the framework in other fits. For example, in comparing different orders of QED evolution in generating CT18qed with the CT18lux inelastic photon in the left plot Fig. 16, we find that the NLO DGLAP kernel gives better agreement than the LO one. It is understood that the LUX formula for the photon PDF includes up to one perturbative order higher than the traditional DGLAP approach [19]. We also explained already that the CT18qed gives a larger inelastic photon at small x but significantly smaller at large x.
The LO and NLO [O(αα s ) and O(α 2 )] QED inelastic photon of CT18qed, together with x Ratio to CT18lux the charge-weighted singlet distribution, Σ e = i e 2 i (q i +q i ), are compared against CT18lux in Fig. 16. We see the redistribution of the proton momentum to the inelastic photon only impacts the inelastic photon and charge-weighted singlet by the corresponding overall factor, which is negligible once compared to the impact of the QED splitting, q → qγ. When turning on the NLO QED evolution, the inelastic photon receives negative corrections, while the remaining Σ e becomes larger as less photon is radiated off quarks. More specifically, the O(αα s ) corrections reduce γ inel by approximately 6% around x ∼ 0.02, and the O(α 2 ) corrections reduces this by another ∼1%. The O(α 2 ) corrections to Σ e are effectively invisible, as they coincide with the O(αα s ) effect in Fig. 16, right plot. This comparison clearly shows that the impact of including higher-order QED effects into the DGLAP evolution equations is much larger than fixing the amount of momentum violation due to the incorporation of the photon PDF.

D. QED evolution in the global fit
In principle, the QED-corrected DGLAP evolution discussed above can play a role in the CT global analysis. In order to examine this possibility, we perform a new global fit, in which the QED corrections to the splitting kernels are included alongside the NNLO QCD contributions, and the QED corrections to DIS structure functions are also applied, similarly to the treatments of NNPDF3.1luxQED [27], MMHT2015qed [28], and MSHT20qed [29]. The resulting χ 2 values for fitting to the CT18 default data sets are presented in the rightmost column of Tabs. III-IV, and the corresponding PDFs in Fig. 17. In order to compare with the CT18 NNLO global fit on the same footing, the initialization scale is chosen to be the same, i.e., we set µ 0 = 1.3 GeV. The results of the unfitted CT18lux (having the same χ 2 as CT18 NNLO) and CT18qed calculations are also shown for comparison. In particular, the χ 2 values for CT18lux and CT18qed are given explicitly in Tabs. III-IV in the fourth and fifth columns from left, respectively. As we demonstrated with the example of the chargeweighted singlet PDF in Fig. 12 and right panel of Fig. 16 already, introducing the photon as a new component within the QED evolution takes away only a very small fractional momentum from the quark PDFs, especially at large x due to the (q → qγ) splitting. As a result, we observe a small increase in χ 2 under the unfitted CT18qed calculation, amounting to ∆χ 2 ≈ 10 units out of ≈ 4293 for the full data set consisting of N pt = 3681   [31]. The CT18lux shares the same χ 2 as the CT18 PDF, as the quark and gluon PDFs remain unchanged.
points. In terms of specific data sets, the theoretical description remains unchanged for most experiments, while the χ 2 values of those data most sensitive to the d/u PDF ratio (e.g., the CMS 8 TeV muon charge asymmetry [86]) or large-x PDF behavior (e.g., LHCb 7 TeV forward W/Z production [84] and BCDMS F p 2 [59]), are shifted by a few units. In the new refit carried out in the presence of QED corrections (corresponding to the 'QED fit' column of Tabs. III-IV), the total χ 2 was reduced marginally, but, in the end, the description of the global data set remains substantially unchanged. In terms of the refitted PDFs of specific parton flavor, we see in Fig. 17 that the reduction to the u-and d-quark PDFs that occurred under CT18qed1.3GeV is mostly restored to match the behavior seen with CT18lux, as preferred by the fitted data. In this case, the reduction to the total partonic momentum carried by the quarks and gluon due to the presence of the photon   PDF derives mainly from the gluon and sea quarks, roughly by an overall factor. This momentum is taken by the photon PDF, which receives a slight enhancement at the level of xγ (µ 2 0 ) ∼ 0.21%, compared with CT18qed, in which the momentum fraction is only taken from the sea quarks manually. In this sense, the role of the momentum sum rule makes the photon momentum fraction in the QCD+QED fit lie in between those obtained by the redistribution of the gluon momentum adopted by LUXqed (17) [18,19] and of the sea quark momentum adopted in our default CT18qed.
The analysis above serves as validation of the boostrap logic leading to CT18qed, in which the QED evolution can be viewed as a minor perturbation giving only a very small variation in the quark and gluon PDFs from the CT18 global minimum. We therefore regard the publicly-released CT18qed set(s) as our primary QED calculation in light of this small change.

V. IMPLICATIONS FOR PHOTON-INITIATED PROCESSES AT THE LHC
Having developed and applied a combined QCD+QED formalism to obtain the photon PDF in the preceding sections, we now examine the sensitivity of observables at the LHC to this photon PDF and explore the phenomenological consequences. Before going to specific SM processes, we first compute the parton-parton luminosities involving the photon. The definition of the parton-parton luminosities can be generically taken as [94] L ij (s, M 2 ) ≡ 1 s where τ = M 2 /s. Usually, the scale is chosen as µ 2 = M 2 . The parton-parton luminosities L γγ , L γΣe , L γg at a 13 TeV pp collider are shown in Fig. 18. First, we examine the CT18lux parton luminosities separately calculated from the elastic and inelastic contributions to the photon PDF compared with the total photon PDF in the upper left of Fig. 18. We see that the elastic photon makes a sizeable contribution in both the low-and high-M limits, due to the relatively large size of the elastic photon at lower scales as well as in the large-x region. We also compare parton-parton luminosities based on the total photon PDF input but using different existing PDF frameworks in Fig. 18. In general, we find CT18lux agrees quite well with LUXqed17 and NNPDF3.1luxQED, whereas CT18qed gives comparatively smaller parton luminosities and MMHT2015qed is somewhat larger. At higher invariant masses, M > 1 TeV, CT18qed and MMHT2015qed give smaller values of L γγ , while NNPDF3.1luxQED is larger, a feature we trace to the difference between the low and high values for the initialization scales in the DGLAP vs. LUX approaches, respectively. In comparison, MMHT2015qed yields a larger L γg while NNPDF3.1luxQED is smaller, as a result of the large-x gluon behaviors. In the following, we examine the impact of the photon PDF upon collider phenomenology as represented by a number of Standard Model processes. For these, we take the production of high-mass Drell-Yan and W + W − pairs, Higgs-associated W + production, and tt pair production as typical examples sensitive to L γγ , L γΣe (or L γΣ ), and L γg .

A. High-mass Drell-Yan production
We start with high-mass Drell-Yan production, which has been extensively measured by both the CMS [95] and ATLAS [96] experiments at the LHC. Representative Feynman diagrams for this process are shown in Fig. 19. In our theoretical predictions, we generally    assume the ATLAS 8 TeV fiducial cuts [96], in order to directly compare with data as illustrated below. The quantity appearing inside paraentheses in Eq. (30) indicates the transverse momentum cut on sub-leading leptons. The NNLO QCD and NLO EW 13 corrections are already examined in the CT18 NNLO global analysis as documented in Ref. [31]. We repeat this same computation here, but now with single and double photon initiated contribution(s)i.e., processes involving ligh-by-quark and light-by-light scattering as shown in Fig. 19-updated according to the various existing photon PDF calculations. The absolute differential cross section for production of Drell-Yan pairs, dσ/dm ¯ , is plotted in the upper left of Fig. 20. In particular, we observe that, with increasing di-lepton invariant masses over the range 116 < m ¯ < 1500 GeV, the absolute cross section monotonically decreases through ∼5 orders-of-magnitude. To better illustrate variations in the cross section of the electroweak corrections, we therefore normalize dσ/dm ¯ to the NNLO QCD calculation as shown in Fig. 20, upper right. In this context, we see that NLO EW corrections are responsibly for an approximate ∼5% negative effect around m ¯ ∼ 1 TeV. The inclusion of double-photon-initiated (DPI) processes, on the other hand, largely counteracts this effect, increasing the absolute cross section by a similar, ∼5% shift within this m ¯ region. Next, we compare the purely DPI cross sections obtained with different assumed photon PDFs in Fig. 20, lower plots. As mentioned before, CT18qed refers to the DGLAP-driven calculation using an initialization scale of µ 0 = 3 GeV. In this case, we see the CT18qed result lies about 2% below the CT18lux prediction, a fact which can be understood in terms of the behavior of the photon PDF itself: the DPI cross section goes approximately as ∼ 2 × γ(x), with the photon PDF being just under ∼1% smaller in CT18qed relative to CT18lux, as shown in Fig. 10. Meanwhile, CT18lux gives almost the same predictions as LUXqed17, while NNPDF3.1luxQED produces a slightly smaller DPI cross section at low m ¯ , but higher at large m ¯ ; this can be attributed to the difference between the DGLAP evolution and the LUX approaches. The MMHT2015qed cross section is 3∼5% larger than that based on CT18lux. We note also that the PDF uncertainty for CT18lux, CT18qed, and MMHT2015qed are roughly the same, lying within the range 1∼2%. In addition to these comparisons based on our CT18 QED results and other recent calculations, we also compare DPI cross sections based on the previous generation of photon PDFs, namely, MRST2004, NNPDF2.3, NNPDF3.0, and CT14qed, which we plot in Fig. 20, lower right, normalized to CT18lux as a reference. The central predictions of CT14qed give quite strong agreement with CT18lux, while the uncertainty is about 20%. The MRST2004qed gives slightly larger uncertainty, which is about 30∼50%. The NNPDF2.3 and NNPDF3.0 calculations give significantly larger predictions at higher invariant mass, and the size of uncertainty bands can be as large as 100∼200% in this latter case.
We remind the reader that a major purpose of this phenomenological study is an exploration of the photon PDF uncertainty. In order to isolate the photon's impact, we focus on double-photon-initiated (DPI) dilepton production and the theoretical calculation is performed with leading-order matrix elements for which the representative diagram is the third graph of Fig. 19. It was pointed out in Ref. [97] that this DPI process suffers from a large scale uncertainty, which is found to be ∼10 − 20%, shown as the inset in the lower-right panel of Fig. 20. Starting at next-to-leading order, we begin to have contributions from the single-photon-initiated process (i.e., qγ scattering, as shown in the last diagram of Fig. 19), which has an overlap with the NLO EW diagrams for Drell-Yan production (e.g., the second diagram of Fig. 19). A systematic treatment to include these dynamics simultaneously at this order requires a proper subtraction to avoid double-counting, which is beyond the scope of this work. Instead of taking a partonic factorization picture, it is also demonstrated in Refs. [98,99] that the photon-initiated processes can be formulated using the structure function approach, which avoids these artificial scale uncertainties and provides a precise definition of the proton's photon content.

B. W H production
At pp colliders with sufficient √ s, W -boson-associated Higgs (i.e., W H) production can proceed through a Drell-Yan-like mechanism mediated by W -boson exchange, as depicted in the diagram shown in Fig. 21, upper left. At one higher EW (or QED) order, we can also have contributions from SPI processes like that appearing in Fig. 21, right diagram. The QCD calculation for W H production can be achieved at NNLO with MCFM [100] or vh@nnlo [101], while the EW corrections considered here can be performed by means of the HAWK package [102] or the general-purpose generator MadGraph aMC@NLO [3]. Here, we consider the total inclusive 14 cross section for W + H production at a 13 TeV pp collider as a demonstration. Following these considerations, we plot the absolute differential cross section, dσ/dM W H , in Fig. 22, upper left, performed here at NNLO in pQCD and with NLO EW effects based on the CT18lux PDF. We observe that the absolute differential cross section drops as drastically as by four magnitudes when the W H invariant mass increases up to 2 TeV. At large invariant mass, the size of NLO EW correction can be as significant as O(1) compared with the pure NNLO QCD prediction.
In order to examine the importance of the photon-initiated contribution, we show the ratio of the SPI cross section to the total one, i.e., NNLO QCD and NLO EW, in Fig. 22, upper right. We see that right above the threshold around M W H ∼200 GeV, the SPI processes only contribute about 1% to the total cross section. In contrast, when the M W H increases up to 2 TeV, the SPI contribution becomes 60% -exceeding even the pure QCD cross section -which highlights the importance of the photon contribution. In addition, we also show the PDF uncertainty, which is about 2%, significantly reduced when compared with the first generation ones shown in Fig. 22, lower plot.
The total cross sections of exclusive or quasi-exclusive W + W − production at the LHC via pp → p ( * ) W + W − p ( * ) → p ( * ) µ ± e ∓ p ( * ) has been measured by both the CMS [103,104] and ATLAS [105] experiments. The experimental cross sections are corrected to the full phase space, with results summarized in Tab. V. The theoretical predictions are evaluated as where pp → pW + W − p represents the elastic scattering process in which the proton remains intact. The Monte-Carlo simulated cross sections are taken as γγ → W + W − → µ ± e ∓ X with photon (γ) distribution function being calculated from the equivalent photon approximation (EPA) [11]. 15 The measured values and theoretical predictions are listed in the second and third column, respectively, of Tab. V. The branching ratio of W + W − pairs decaying into µ ± e ∓ X, including τ leptonic decays was taken as BR=3.23% [107]. 16 The dissociation factor, F , reflects the effect of including also the "quasi-exclusive" or "proton dissociation" production contribution, was extracted from the CMS and ATLAS data of high-mass lepton-pair 14 By "total inclusive," we mean the full phase space without any fiducial cuts; the W and H bosons are assumed to be on-shell final states without subsequent decay. 15 CMS takes CalcHEP and MadGraph, respectively, for the predictions of 7 TeV and 8TeV, while ATLAS adopts Herwig++, resulting in a smaller cross section due to different implementations of EPA. Her-wig++ uses the integration of the dipole form factors [106], while both MadGraph and CalcHEP directly implement the approximated form [11]. 16 In Ref. [108], the branch ratio of W + W − → µ ± e ∓ X was taken as a slightly different value, BR=3.1%.   production γγ → + − . We note that, as found in Ref. [16], underlying event contributions from subsequent rescatterings of the intact proton are expected to be small in the exclusive, non-dissociative channel and are not explicitly calculated here. With the F and BR, we can correct the Monte-Carlo simulated EPA prediction to the exclusive cross section at the W + W − undecayed level, listed as the EPA (fifth) column in Tab. V.
The implications of the CMS (quasi-) exclusive W + W − production data were investigated in Ref. [109], in which the CMS data was found to be in good agreement with a theory prediction including both elastic and single-dissociative contributions predicted by CT14QED and CT14QEDinc PDFs [8]. Here, we compare in Tab. V the predictions of elastic photons in CT18lux (same as CT18qed) and MMHT2015qed to CMS and ATLAS data on the exclusive γγ → W + W − production cross section to the EPA predictions, as shown in the last two columns of the table. We find that MMHT2015qed yields an enlarged (relative to CT18lux) cross section due to its larger elastic photon, as shown in Fig. 28. We also quantify the uncertainty due to potential variations in the elastic photon contribution as discussed  V. The exclusive/quasi-exclusive pp → p ( * ) W + W − p ( * ) → p ( * ) µ ± e ∓ p ( * ) production cross sections σ [fb] measured by CMS [103,104] and ATLAS [105] groups. in Sec. III B, with MMHT2015qed giving a slightly smaller uncertainty. In addition, we also show the invariant mass distribution, dσ/dM W W , for exclusive W + W − production in Fig. 24. At larger invariant masses, the elastic photon of MMHT2015qed results in a similar uncertainty as CT18lux, although the central value is about 10% larger. Recently, the ATLAS group released a new measurement of this process based on 139 fb −1 of 13 TeV LHC data [110], which rejects the background-only hypothesis with a significance of 8.4 standard deviations. Unlike the 7 and 8 TeV measurements, however, this analysis was performed in the fiducial volume. An investigation of the implications of these data is deferred to a future study.

D. Top-quark pair production
Finally, we explore the impact of the photon PDF on tt production. Representative diagrams for tt production proceeding through the leading order QCD channels (e.g., gluon fusion) and single-photon-initiated processes (e.g., photon-gluon fusion) are shown in the left and right panels of Fig. 25, respectively. As before, we consider tt production at a 13 TeV pp machine for specificity, with the inclusion of NNLO QCD and NLO EW (QCD×EW) corrections as calculated in Ref. [111]. In this work, we examine the SPI cross section, whose ratio to the QCD×EW cross section is shown in Fig. 26. We see by its relative size that the SPI process only contributes less than 0.6% of the full QCD×EW cross section. In this sense,   the SPI contribution is safely ignorable in tt production from the perspective of contemporary precision, unlike what we observed for high-mass Drell-Yan and W -boson-associated Higgs production. The SPI uncertainties of CT18lux, CT18qed, and MMHT2015qed are roughly similar and about ∼3% in this range, constituting a significant reduction when compared with the first generation of photon PDFs shown in Fig. 26, right plot. Needless to say that the contribution of photon-photon fusion to the production of top-quark pairs at the 13 TeV LHC is negligible.

VI. CONCLUSION
In the present analysis, we have carried out a first implementation of the recent LUX QED formalism into the larger CT PDF analysis framework. The LUX approach represented a significant advance in the consistent determination of the photon content of the proton with minimal underlying model assumptions. The LUX QED methodology was subsequently interfaced with several QCD global analysis frameworks, each making unique choices regarding implementation strategy and inclusion of physics ingredients in the ultimate calculation. For this reason, as well as the recent development of QED splitting kernels at O(αα s ) and O(α 2 ), it is time to update the CT analysis family with the dedicated QED study presented here. In particular, in this study, we examined possible systematic differences that arise between following each of the two main approaches carried out in recent LUX-based photon PDF calculations. Broadly, these consist either of computing γ(x, µ 2 ) according to the LUX master expression in Eq. (2) at an arbitrary scale often chosen to be µ ∼ 100 GeV for iteration in a global fit, e.g., NNPDF3.1luxQED [27], or of instead evaluating the photon PDF at an initial scale, γ(x, µ 2 0 ), by applying a slight modification of the LUX master formula and relying entirely on QED+QCD evolution to determined the photon PDF at µ > µ 0 , e.g., MMHT2015qed [28]. We adapted the latter of these two approaches to our CT methodology, leading to the CT18qed photon PDF set, which we regard as the primary result of this analysis. For the sake of comparison, we also implement the former approach, in which the photon PDF is everywhere computed according to the unmodified LUX master formula, designating the result CT18lux. We release this alternative calculation alongside CT18qed, and we have compared the two against other recent determinations in the present work.
Together, the CT18qed and CT18lux photon PDFs we have produced involve several novel and unique features, in addition to the fact that we have considered both approaches comprehensively within a single framework. An important aspect of the CT18 photon PDFs is the updated final uncertainty we report, which follows from a critical appraisal of possible error sources following a number of physics updates that have occurred since the original LUX publication. These include a reassessment of the elastic form factor uncertainty, possible higher-twist and target-mass effects, and updates to the description of the proton structure function inside the resonance region as canvassed in Sec. III B 2. While we compared our CT18qed and CT18lux PDFs extensively in Sec. III and IV, a number of qualitative features are worth highlighting. In particular, we obtain a strong general agreement among the calculations most closely aligned with the original LUX approach, as can be seen in Fig. 10, which illustrates the close similarity of CT18lux relative to the LUXqed17 and NNPDF3.1luxQED results. The comparative photon PDF uncertainties among these calculations show some mild differences in Fig. 13, especially at low x, for which our CT18lux calculation is essentially intermediate between the smaller NNPDF3.1luxQED band and the larger low-x uncertainty reported in LUXqed17 [19]. In the CT incarnation of the DGLAP approach, we observe some intriguing differences in the shape and magnitude of our final PDFs, with our final pair of CT18qed PDFs, with initialization scale µ 0 = 1.3 GeV and 3 GeV, somewhat underhanging CT18lux, for example, especially for x > 10 −3 . While we again refer interested readers to Sec. III and IV for an in-depth dissection, we point out that this behavior goes in much the opposite direction compared with that of MMHT2015qed.
In Sec. V, we traced a number of phenomenological consequences of the CT18 photon PDFs, with a particular emphasis on pp collisions at the LHC. We presented TeV-scale parton-parton luminosities, distributions for high-mass Drell-Yan production, as well as W H, W + W − , and tt cross sections. For the parton-parton luminosities shown in Fig. 18, we found generally robust concordance, up to uncertainties, among the various calculations explored in this analysis -albeit with some evidence of deviation in L γγ at M > 10 GeV, especially for MMHT2015qed, and, to a lesser extent, CT18qed. Similarly, we report reasonable agreement among the phenomenological calculations of the invariant-mass distributions computed for high-mass Drell-Yan, W H, W + W − , and tt, production, again finding a frequent excess, especially for MMHT2015qed, as follows from its comparatively larger photon PDF. Of phenomenological importance, the PDF uncertainties in these processes can still be significant and are linked to the many nonperturbative error sources explored in this study. Achieving the elevated precision of the perturbative calculation (to percent-level accuracy) needed for Beyond Standard Model (BSM) searches in, e.g., the tails of invariant-mass dis-tributions will therefore require improvements to these lower-energy inputs to the photon PDF calculation.
As a companion to this article, we publicly release LHAPDF6 [112] grids corresponding to the two main calculations presented above: CT18lux and CT18qed (including CT18qed1.3GeV as well) in the HEPForge repository, https://ct.hepforge.org. Again, we stress that these grids include both our newly-calculated photon PDF within each approach as well as the accompanying (anti-)quark and gluon distributions, with uncertainties quantified according to the Hessian approach combined with low-Q 2 resources as Eq. (10). These uncertainties represent those associated with the underlying quark and gluon degrees-of-freedom, as well as the collection of uncertainty sources reviewed in Sec. III B 2. As pointed out above, we identify our CT18qed (with µ 0 = 3 GeV) calculation as a 'first among equals' result and advocate its primary use in phenomenological calculations like those shown in Sec. V. The DGLAP evolution of photon simultaneously with quark and gluon PDFs provides a consistent description of the photon-initiated processes together with the contribution of quark (and possibly gluon as well) partons in the perturbative expansion. As an alternative, the CT18qed1.3GeV PDFs, in which the photon PDF is initialized at the CT18 starting scale of µ 0 = 1.3 GeV, are more appropriate for describing the photon in the low-energy range 1.3 < µ < 3 GeV but give a larger uncertainty at large x values. The remaining set that we have released simultaneously, CT18lux, uses a method that provides a particularly useful determination of the inclusive photon.
Finally, we remind the reader that a number of aspects of this study, among them, more detailed investigation of parton-level charge-symmetry breaking [113] and simultaneous determinations of the neutron's photon content, leave room for further improvement. Rather than undertaking these in the analysis above, we reserve these issues for future work(s), for which such considerations maybe relevant to achieving still higher electroweak precision in next-generation phenomenology.

ACKNOWLEDGMENTS
We thank Eric Christy for providing updated code to compute the low-W 2 resonanceregion structure functions, as well as for valuable conversations. We also thank our the CTEQ-TEA colleagues for helpful discussions and Sergei Kulagin for useful exchanges.

APPENDIX Appendix A: Separation of elastic and inelastic photon PDF components
As discussed in Sec. I and II, the photon PDF consists of elastic and inelastic components. To illustrate, we show the absolute values and corresponding PDF fractions of these two components for CT18lux at a number of scale choices 17 , namely, µ = 1.3, 10, 10 2 , 10 3 GeV in left and right plot of Fig. 27, respectively. At low energy, e.g., µ 0 = 1.3 GeV, the inelastic photon is mainly determined by the low-Q 2 structure functions F 2,L directly measured in low-energy experiments, such as HERMES [35] in the continuum region and CLAS [33] (or Christy-Bosted [34]) in the resonance region, as illustrated in Fig. 1. As the scale increases, the inelastic photon receives a significant contribution from the quark splitting, q → qγ, with a log(µ 2 ) enhancement. As a result, the absolute value of the inelastic photon PDF increases drastically, and, the corresponding fractional share of the inelastic photon PDF relative to the total comes to dominate at higher µ 2 , as shown in Fig. 27.
Similar to LUXqed (17), the elastic photon in CT18lux and CT18qed is fully determined by the elastic form factor, through in which δxγ el corresponds to the integrand in Eq. (A1) divided by x. The numerical difference between these two prescriptions is given in Fig. 28. Overall, MMHT2015qed gives a larger elastic photon than LUX (and therefore, CT18lux and CT18qed). Furthermore, the ratio increases along with the scale µ. This is mainly a consequence of the fact that MMHT2015qed only includes quark contributions in p γγ , whereas LUX includes both quarks and leptons in the α(µ 2 ) running. We also notice that, at large x and small µ 2 , MMHT2015qed gives a smaller γ el . This is due to the equivalent upper integration limit, µ 2 , in the MMHT2015qed method, which is different from ∞ in Eq. (A1). Especially at large x, the lower limit, x 2 m 2 p /(1 − x), will approach µ 2 , which leaves the integration over [x 2 m 2 p /(1 − x), µ 2 ] significantly smaller than the one over [x 2 m 2 p /(1 − x), ∞]. However, the effect of this relative difference in integration intervals becomes smaller with increasing µ 2 . however, when Q 2 is integrated up to infinity. Following the standard renormalization procedure, we can split the Q 2 integral into two parts, m 2 p x 2 /(1 − z) < Q 2 < µ 2 /(1 − z) and µ 2 /(1 − z) < Q 2 < ∞, which correspond to the physical factorization (PF) and MS conversion (con), respectively, leading to two contributions to the photon PDF, γ(x, µ 2 ) = γ PF (x, µ 2 ) + γ con (x, µ 2 ) . (B2)