Physics with Leptons and Photons at the LHC

We present a phenomenological study of photon-initiated (PI) lepton production at the LHC, as implemented in the structure function (SF) approach. We provide detailed predictions for multi-differential lepton pair production, and show that the impact on observables sensitive to the the weak mixing angle, $\sin^2 \theta_W$, and $W$ boson mass, $M_W$, as well as PDFs can be non-negligible, in particular given the high precision being aimed for. The SF calculation can provide percent level precision in the corresponding production cross sections, and is therefore well positioned to account for LHC precision requirements. We consider the pure $\gamma\gamma$ channel, and compare in detail to the NLO collinear calculation. We in addition include initial-state $Z$ as well as mixed $\gamma /Z + q$ contributions, and assess their impact. We also consider photon-initiated lepton-lepton scattering, and again find the SF approach can provide high precision predictions for this process in way that can straightforwardly account for any fiducial cuts imposed. Finally, we provide a publicly available Monte Carlo generator, SFGen, for PI lepton pair production and lepton-lepton scattering within the SF approach, for use by the community.


Introduction
A major aim of the LHC, and the high-luminosity upgrade (HL-LHC) that will follow, is to precisely test the Standard Model (SM) predictions for as wide a range of collider processes as possible. A particularly important element of this involves events with leptons in the final state. To give two key examples, the Drell Yan production of a lepton pair measured multidifferentially allows a measurement of the weak mixing angle, sin 2 θ W , as well as providing constraints on the parton distribution functions (PDFs) of the proton, while precision data on W production via leptonic decays can provide a correspondingly precise measurement of the W mass, M W . Significant amounts of such data have already been collected at the LHC, see e.g. [1][2][3][4][5][6][7][8][9][10]. In more detail, the ATLAS 7 TeV analysis [6] presents a measurement of M W with a total uncertainty of 19 MeV, which is comparable to previous Tevatron determinations. High precision triple-differential measurements of lepton pair production have also recently been reported at 8 TeV by both ATLAS [7] and CMS [8], with the latter reporting a value of sin 2 θ W that is beginning to bear down on LEP precision. As discussed in [11], significant improvements on the precision of these measurements are anticipated and aimed for in the future, in particular from HL-LHC running.
A key ingredient in this programme is the availability of high precision theoretical predictions for the SM production processes. An important and topical element of this is the contribution from photon-initiated (PI) channels to lepton pair production. A great deal of theoretical progress has occurred in the modelling of this process, until recently within the context of including the photon as an additional partonic constituent of the proton, that is via a photon PDF. Following earlier work of [12][13][14][15][16][17], the LUXqed group [18,19] applied the method originally appearing in the context of the equivalent photon approximation (EPA) [20] for determining the photon flux from the elastic and inelastic proton structure functions, providing the first publicly available PDF set within this approach (see also [21][22][23][24] for earlier work). Following from this, photon PDFs have been provided in combination with the global MMHT [25] and NNPDF [26] sets, applying the same approach as LUXqed. In all cases, the experimental input for the corresponding structure function inputs is sufficiently precise that the quoted photon PDF uncertainty is very small, generally at the ∼ 1% level, and thus a high precision photon PDF determination can be claimed.
This however is not the end of the story. While the photon PDF itself is known with high precision, the corresponding PI cross section has at leading order (LO) a sizeable factorization scale variation uncertainty associated with it, that is significantly larger than the small PDF uncertainty. While this situation can be, and should be, improved upon by working to NLO or beyond, as discussed in detail in [27], an alternative approach is to work in the 'structure function' (SF) framework. Here, one calculates the cross section directly in terms of the proton structure functions, in an analagous manner to the calculation of Higgs boson production via vector boson fusion (VBF) [28][29][30][31][32][33]. This contains no reference to the photon PDF, and it was demonstrated in [27] for the specific cases of a toy lepton-hadron scattering processes and off-Z peak lepton pair production at the LHC that this provides predictions which are by construction more precise than calculations within the standard collinear factorization approach. Percent level precision in the predicted cross sections was in particular achieved for the first time, with remaining contributions from e.g. the 'non-factorizable' corrections that break the structure function picture expected to be small, as they are in the VBF Higgs case (see [34,35]).
In this work, we build on this previous study in a number of directions: we investigate in more detail the implications for precision phenomenology at the LHC, calculating the PI contribution to multi-differential dilepton production, which is relevant for determinations of sin 2 θ W as well as PDF constraints, and evaluating the potential impact on W mass determinations, through the tuning that is performed to the Z event sample, from which the PI contribution must be subtracted; we include contributions from Z bosons in the initial state, as well as photons, and the mixed Z/γ + q contributions, within the SF framework, which while generally small can be important in some regions; we calculate the PI cross section for lepton-lepton scattering and consider the interplay this has with the recently developed lepton PDF formalism [36]; finally, we provide the above results in the publicly available SFGen Monte-Carlo (MC) generator, for use by the community.
In more detail, we first investigate the relative impact of PI on the dilepton rapidity distribution at the LHC, which enters at the level of 5−10% of the QCD DY contribution away from the Z peak region. We in addition compare in detail the high precision SF results for the pure γγ PI process, which dominates away from the Z peak, with the calculation in collinear factorisation, at both LO and NLO; the latter order corresponds to the highest currently available precision calculation within the collinear approach. The LO collinear predictions as expected exhibit large scale variation uncertainties, which are significantly larger than the small PDF uncertainty on the photon PDF, and the related experimental uncertainty on the SF calculation, due to the determination of the structure functions. If we take µ = m ll for the renormalization/factorization scale, the agreement between the NLO collinear prediction and the SF is indeed improved with respect to the LO case, and the scale variation uncertainties reduced, although at low to intermediate masses these remain significantly larger than the PDF uncertainty. Moreover, taking a differing scale choice related to the lepton p l ⊥ , the convergence of the NLO collinear result is greatly reduced, and the agreement with the full SF result rather poor. Thus, we find that while in general the inclusion of NLO corrections can improve the precision and accuracy of the collinear result, the overall picture is not completely clear cut; in some cases there remain reasonably large scale variation uncertainties and areas of disagreement with the full SF result, depending on the choice of scale and/or experimental cuts.
We will in addition show that the PI contribution of the off-Z peak ATLAS 8 TeV triple differential lepton pair data [7] is expected to enter at up to the ∼ 6% level in comparison to the QCD DY cross section, with a non-trivial impact on the corresponding cos θ * distribution. Similar results will be expected for other measurements of this process. This can therefore potentially impact on sin 2 θ W determinations as well as PDF constraints. We also consider the impact of PI contributions on dilepton production in the Z peak region, showing that while inclusively these enter at the per mille level, at lowest p ll ⊥ values this is rather larger, at the ∼ 0.5% level. This is a result of the suppression of QCD DY process in this region, and the enhancement of PI production, due in part to the contribution from elastic photon emission; this can only be predicted reliably by using the SF approach directly. Given the level of precision being aimed for in W mass determinations and the fact that the dilepton low p ll ⊥ region is tuned to in order to extract M W from lν events (for which PI production of this sort is absent), this may in principle have an impact on such determinations. That is, the irreducible PI component must be subtracted from the dilepton tuning sample, and the effect of this necessarily depends on the modelling of PI production. We provide a quantitative estimate of this possible impact, and show that while small it is not necessarily negligible with respect to the uncertainties in M W being aimed for at the LHC.
The method for including contributions from Z bosons in the initial state, and the mixed Z/γ + q contributions follows from a straightforward generalisation of the pure γγ case. As mentioned above, the impact from these is generally small, but this is not alway true: for example, at higher dilepton invariant masses the contribution from initial-state Z bosons is not negligible, while when considering the lepton pair p ll ⊥ distribution, Z/γ + q contributions can be important. The treatment of these contributions raises a wider point, namely that the benefit of applying the SF approach directly, while transparent for process where the final state of interest (l + l − , W + W − ...) is directly produced by the γγ initial state, is less clear in the mixed γ + q case. Here, one must deal with the collinear enhancement of the γ → qq splitting, and at this level of precision include QED DGLAP evolution of the quark/antiquark PDFs, which certainly requires one introduce a photon PDF within the LUXqed approach. Nonetheless, one may systematically account for this while still working within the SF approach. We in particular discuss the complication that arises due to double counting with the PI contribution that is already effectively included via the QED DGLAP evolution of the quark PDFs, and provide a simple procedure for removing this. The formal precision of the SF calculation is however then tied to the order at which one does this subtraction.
In a recent study [36] it has been demonstrated how a lepton originating from a γ * → l + l − splitting can be viewed as an initial-state parton in the production process, and hence accounted for via a lepton PDF in the proton, which is determined to high precision in a manner analogous to the LUXqed case for the photon (see [37] for earlier work). As with the photon PDF, one can instead calculate the corresponding cross section in the SF approach directly, including the γ * → l + l − transition via this. We will demonstrate this for the specific example of samesign same-flavour, or different-flavours arbitrary-sign lepton pair production. We find that the SF approach can indeed provide high precision predictions for this process, and the result of applying cuts in order to isolate back-to-back lepton topologies is readily evaluated.
Finally, the SFGen MC provides a publicly available tool for lepton pair and lepton-lepton scattering within the SF approach, including initial-state Z and mixed Z/γ + q contributions. Arbitrary distributions and unweighted events can be generated, and errors due to the experimental uncertainty on the structure function (the equivalent of PDF uncertainties in the photon PDF framework) can be calculated on-the-fly. Unweighted events an be interfaced to Pythia for showering/hadronization of the proton dissociation system, following the approach discussed in [38].
The outline of the paper is as follows. In Section 2.1 we overview the SF approach and discuss the inclusion of initial-state Z contributions. In Section 2.2 we show how mixed Z/γ + q contributions are included. In Section 3.1 we present results for the dilepton invariant mass distributions. In Section 3.2 we present a comparison for these distributions between the SF calculation and the collinear at LO and NLO. In Section 3.3 we consider the impact on the dilepton p ll ⊥ distribution. In Section 3.4 we consider the potential impact of the PI contribution to dilepton production on W boson mass determinations. In Section 3.5 we present results for the PI contribution to triple-differential lepton pair production. In Section 4 we provide predictions for lepton-lepton scattering and compare with the lepton PDF framework. In Section 5 we describe the SFGen MC and its availability. Finally, in Section 6 we conclude.

Photon and Z boson initiated production
The expression for the photon-initiated cross section in proton-proton collisions is given in [27], which itself is based on the analysis of [20]. Here we straightforwardly generalise this to include Z exchange and Z/γ interference: where Y i = γγ, Zγ, ZZ corresponds to pure PI, γZ interference and pure Z-initiated production, respectively. The outgoing hadronic systems have momenta p 1,2 and the photons have momenta q 1,2 , with q 2 1,2 = −Q 2 1,2 . We consider the production of a system of 4-momentum k = q 1 + q 2 = N j=1 k j of N particles, where dΓ = N j=1 d 3 k j /2E j (2π) 3 ≡ N j=1 dΓ j is the standard phase space volume. M µν corresponds to the relevant V V → X(k) production amplitude, with V = γ, Z and the 'Y 1 , Y 2 ' subscript indicating that pure Z, γ or Z/γ interference be included in the appropriate way. ρ is the density matrix, which is given in terms of the well known proton structure functions: where p i is the proton i momentum, for a hadronic system of mass M i and we note that the definition of the photon momentum q i as outgoing from the hadronic vertex is opposite to the usual DIS convention. Here, the integral over M 2 i is understood as being performed simultaneously with the phase space integral over p i , i.e. is not fully factorized from it (the energy E i in particular depends on M i ). The prefactor η is given by For the structure functions, we include the corresponding elastic and inelastic contributions in the γγ case and include an uncertainty due to the experimental inputs on these, as described in [27]. These are evaluated following the procedure discussed in [25], which is closely based on that described in [18,19]. We refer the reader to these references for further details, but in summary we include: an uncertainty on the A1 collaboration [39] fit to the elastic proton form factors, based on adding in quadrature the experimental uncertainty on the polarized extraction and the difference between the unpolarized and polarized; a ±50% variation on the ratio R L/T , relevant to the low Q 2 continuum inelastic region; a variation of W 2 cut , the scale below which we use the CLAS [40] fit to the resonant region, and above which we use the HERMES [41] fit/pQCD calculation (for Q 2 below/above 1 GeV 2 ), between 3-4 GeV 2 ; the symmetrised difference between the default CLAS and Cristy-Bosted [42] fits to the resonant region; the standard PDF uncertainty on the MMHT2015qed nnlo quark and gluon partons in the Q 2 i > 1 GeV 2 continuum region, as calculated via NNLO in QCD predictions for the structure functions in the ZM-VFNS, implemented in APFEL [43]. We do not include any uncertainty due to non-factorizeable corrections: in the case of VBF Higgs (and Higgs pair) production these are found to enter at the percent level or less [34,35], and while no calculation is currently available for PI production, we can expect a roughly similar result here.
For Zγ and ZZ structure functions, bearing in mind that these will only be at all relevant for higher Q 2 ∼ M 2 Z , we can safely ignore any elastic or low Q 2 resonant contributions. We therefore simply set the structure functions to zero for Q 2 i < 1 GeV 2 , and the uncertainty in this case is purely due to the PDF uncertainty on the pQCD prediction. From the above expressions, one naturally expects that at lower Q 2 i the Z/γ contribution will be strongly suppressed by ∼ Q 2 i /M 2 Z , and hence will generally not be significant. However, as we will see at larger values of m ll , where the phase space for larger Q 2 i is increased, and/or in kinematic regions where larger Q 2 i is preferred, such as at larger p ll ⊥ , they can play an important role. For all numerical evaluations we as in [44] write (1) as with m X ⊥ = m 2 X + p 2 X ⊥ , while −q i ⊥ is the transverse momentum imparted to proton beam i, such that p X ⊥ = q 1 ⊥ +q 2 ⊥ . Here X (= l + l − ...) is the produced object,β is as defined in [44], and the integral over M 2 i in the photon density matrices are again understood as being performed simultaneously with the other phase space integrals.
Finally, we have discussed Z-boson initiated contributions above, but in general one can have contributions from initial-state W bosons, via both the W + W − → l + l − and mixed W ± + q subprocesses. However, these will be further suppressed with respect to the Z-initiated channels. In particular, in the W + W − → l + l − case both initial-state particles must be W bosons, and there can of course be no interference with γ-initiated diagrams. In the mixed case, again there is no γ/W interference. This will generally lead to a strong suppression in these channels, and hence we do not consider them here for simplicity. However, these could be readily included within the SF framework.

γ/Z + q initiated production
In addition to pure γ/Z initial states discussed above and shown in Fig. 1 (a), which features a t-channel lepton exchange, for the case of lepton pair production there are also mixed γ/Z + q initiated contributions, as shown in Fig. 1 (b) and (c). These can readily be included within the same approach as above, i.e. with the corresponding initial state γ/Z included via the SF formalism. However, as we will see the inclusion of these type of diagrams is necessarily more subtle, overlapping as it does with QED correction diagrams to the qq initiated process, and hence QED corrections to DGLAP evolution. Therefore, as well will see the benefit of including this type of diagram via the SF approach is less transparent.
Nonetheless, it is useful to explore how we might do so. To include these diagrams, we simply replace Here, at LO we have where A µ is the corresponding γ/Z + q → l + l − + q amplitude, with a collinear initial-state quark/anti-quark carrying proton momentum fraction x B,i . We include both photon and Z initiated production of the final-state lepton pair in the amplitude, while the angle brackets indicate helicity/colour averaging as usual. The above result can be simply read off by matching with the corresponding LO QCD expression for the purely structure function type contribution. We use the G µ EW scheme, with the PDG values [45] for all corresponding inputs. Noting that away from the Z peak region (7) scales as ∼ 1/m 2 ll , one can clearly see that these mixed contributions are suppressed by ∼ Q 2 i /m 2 ll , and hence in the inclusive cross section should only play a minor role. On the other hand, near the Z peak or in kinematic regions, such as at higher p ll ⊥ , where larger values of Q 2 i are preferred and/or in the lower m ll region, the role of these contributions will be larger.
We also note the above expressions do not include any interference between the γγ and γ + q diagrams. These can be effectively included at LO by simply reweighting the cross section according to where the 1, 2 subscript refers to the initiating beam, and the 'LO' indicates that we calculate the corresponding LO quark-level diagrams, such that all initial-state photons are generated via the q → qγ process, i.e. by substituting the corresponding LO QCD expression for the structure functions in (2). We have however checked that the impact of this interference is in general very small, and therefore given that the precise contribution will depend on the details of how we implement such effects (the above prescription is not unique), we do not include them in our results.
To make contact with the usual collinear approximation, we drop Z-initiated contributions in (4) to give We writeβ, defined in [44], asβ where we have decomposed the photon momenta via In the small Q 2 i limit we have where the sum is over the quasi-on-shell photon helicity λ i , and R is the metric tensor that is transverse to the photon and quark momenta, e.g.
when the photon-initiator is from beam 2, and similarly in the alternative case. Here we have taken the small Q 2 i and high energy limit. Then we have as in [27], where ξ i ≥ x i is the photon momentum fraction, and 'PF' indicates this is the physical photon PDF, following the notation of [19], i.e.
A change of variables gives where Γ q = d 3 k q /[2E q (2π) 3 ] is the phase space integral with respect to the outgoing quark and Combining the above we can then indeed write 1 where it is understood that the scale of the coupling in the γq subprocess should be taken as µ in order to match (14). 1 For completeness we also clarify the original discussion on this point in the case of the γγ-initiated cross section, as in Section 4 of [27]. In particular, the relevant change of variables in this case is dq 2 1 ⊥ dq 2 2 ⊥ dx1dx2 = βdQ 2 1 dQ 2 2 dξ1dξ2, which cancels the kinematicβ factor in (9).
Finally, in the PI case we might worry that for the class of diagrams in Fig. 1 (b) there is an IR singularity associated with the region of phase space where the γ → qq splitting becomes collinear. In particular, considering for simplicity the case that the initial-state photon comes from beam 2, in this limit (18) becomes as usual: where P qγ is the usual LO splitting function, k is the momentum of the intermediate quark/antiquark propagator, the sum over q, q is understood to match the qq initial state, and we use the shorthand We can see that the limit on the k integral in (19) extends down to k 2 = 0 and hence we pick up the usual IR singularity that must be regulated and absorbed into the definition of the quark PDF in the usual way; the value of the upper kinematic limit is not important for our purposes, though it will be ∼ m 2 ll . However, if we now consider the corresponding expression in the case of the full structure function calculation (9) then without writing down the rather more complicated general expression, we can simply observe that the relevant integral becomes: which is well defined, and certainly finite. In other words, the IR singularity that is present in the pure collinear calculation is automatically regulated in the SF calculation by the photon virtuality Q 2 2 . Thus naively there is no requirement to absorb a γ → q/q contribution into the quark/antiquark PDF, as there is in the collinear calculation. However, this is not in fact the case, as in the corresponding qq initiated processes, when we consider QED corrections we will pick up IR divergences from q → qγ splittings that must be absorbed. These will generate as usual QED DGLAP evolution of the quark PDFs, and basic consistency requires that ∼ P qγ contributions are also present. As well as being required by this, it is clear that (21), though finite, is logarithmically enhanced, and therefore may be more accurately resummed via DGLAP evolution, at least for larger scales |k 2 | max .
This situation is rather close to the treatment of heavy quark flavours in the initial state, which while technically not requiring a treatment that goes beyond the fixed flavour scheme, often prefer the introduction of heavy quark PDFs in order to provide an accurate description at larger scales Q 2 m 2 q . In the current case, given the smallness of the QED coupling α, the phenomenological importance of this will be smaller, but consistency requires we follow a similar procedure, as discussed above. In analogy to the heavy quark case, we can proceed by including in our calculation the SF result, but subtracting the logarithmic contribution as per (21) in order to avoid double counting with the contribution that is generated from the QED DGLAP evolution of the quark PDFs. In more detail, we note that QED DGLAP evolution of the quark PDFs leads to and similarly for the antiquark, to O(α). Here Q 0 corresponds to the input scale of the PDF set, which we take to be 1 GeV, as in [25]. At O(α) we must then subtract the contribution from our structure function calculation, in order to avoid double counting with the contribution from the DGLAP evolution of the quark PDFs. This µ dependence in this result in particular matches the corresponding µ dependence of the quark PDFs, such that the combined qq and qγ results maintain scale independence at LO in α. Indeed, if one were to estimate the scale variation uncertainty on (23), this would need to be done coherently with the qq channel. Scale variation in the γq contribution alone would overestimate the corresponding uncertainty. Nonetheless, it is the case that this introduces a scale variation uncertainty associated with the γq contribution, which is absent in the pure γγ case. Finally, there is in principle some uncertainty associated with the precise choice of Q 0 , though this is also very small, and the value we take matches naturally that taken in [25]. We note that the generalisation of the above procedure to other processes is straightforward. In particular, we simply include a subtraction by folding the appropriate quark/antiquark initiated cross section with to O(α). If one were to include higher order QCD or EW corrections to the subprocess in the SF calculation, then the order of the above expression should simply be matched to this, provided the order of the DGLAP evolution applied in the non-SF contribution matches this. This subtraction should only be included for suitably inclusive cross section predictions, i.e. in the present case where the qq initial-state can contribute in its LO configuration. In particular, if we consider the dilepton transverse momentum distribution, then no double counting with the purely qq initiated process occurs at the order we consider, and the SF calculation provides the full prediction, with no subtraction required. In the end though, even in the case where technically one can argue as above that it should be included, we will see that the impact of this subtraction on the cross section is relatively minor in most (though not all) phenomenologically relevant regions.
In summary, one can certainly include such mixed diagrams via the SF procedure, but some care is needed to do this consistently. In particular, one needs to take care to suitably subtract any double counting with contributions generated by QED corrections to DGLAP in the pure QCD initial-state contributions. Even at lowest order this is not completely trivial, and certainly if one started to go to higher order in QCD, e.g. including gluon emission from the initial state quark lines in 1 (b), then while the result would still be IR finite, the calculation would become more involved and the subtraction terms would be become more complicated again. Moreover, the formal precision of the SF calculation is clearly then tied to the order at which one does this subtraction. Thus, the benefits of using the SF approach for such a process are less clear, though there is no in principle obstacle to doing this.

Lepton pair production: results
In the following sections we present some selected results for lepton pair production, as produced with the SFGen Monte Carlo generator. We will investigate in particular the impact of quarkinitiated diagrams and initial-state Z contributions by plotting both the pure γγ and the 'total' cross section; the latter includes Z and mixed γ/Z +q contributions (including the corresponding subtraction described in Section 2.2 where relevant), in addition to the γγ. We include the uncertainty due to the experimental uncertainty on the inputs, as discussed in Section 5, for the total cross section, but omit this when showing the γγ result for demonstration purposes.
For comparison we also show the collinear γγ prediction at LO and in places NLO, using the MMHT2015qed nnlo photon PDF [25], with error bands corresponding to the PDF and factorization/renormalization scale variation uncertainty added in quadrature. For the collinear calculation we take µ F = µ R and vary by the usual factor of 2 around the particular choice of µ 0 . In the SF calculation, there is a factorization scale uncertainty in the γ/Z + q contribution, from the quark/antiquark PDF, see (7). For completeness we add this in quadrature to the experimental uncertainty discussed above in the total cross section, but as we will see this is in general a negligible source of uncertainty on the total cross section.
We will in places compare to the LO collinear predictions alone, as used in e.g. [26,46], in order to draw a qualitative comparison. However, the highest publicly available order we can calculate at corresponds to NLO, and this can naturally decrease the corresponding scale variation uncertainties and improve the agreement with the SF result. Nonetheless, as explicitly demonstrated in [27] in the context of lepton-hadron scattering even at NLO the agreement with the SF approach is not perfect. Indeed, with this in mind we will in Section 3.2 present a detailed comparison of the SF result with the highest available precision collinear calculation, namely at NLO. SF, γγ only SF, total Collinear, γγ LO Figure 2: Ratio of the photon-initiated cross sections for lepton pair production to the NLO QCD Drell-Yan cross section at the 13 TeV LHC, as a function of the lepton pair invariant mass, m ll . The leptons are required to lie in the |η l | < 2.5 region. The LO collinear predictions and the structure function result are shown, in the latter case for both the pure γγ initiated and the total. In the former case the uncertainty band due to renormalization/factorization scale variation by a factor of two around the central value µ = m ll , added in quadrature with the much smaller PDF uncertainty, is given. In the latter case, the error band due to the experimental uncertainty on the structure functions and scale variation in the γ/Z + q component is shown added in quadrature (note in some regions this is comparable to the line width of the central value).

Invariant mass distributions
We first consider the dilepton invariant mass distributions, shown in Fig. 2, taking the same mass regions and cuts as in [27]. For the collinear predictions we take µ 0 = m ll . Starting with the low mass region (10 < m ll < 60 GeV) we can see that the initial-state Z contribution is negligible, as at these low masses the average γ/Z virtuality Q 2 i is rather low. In addition, for this inclusive observable we find that the qγ contribution is equally very small.
Turning now to the intermediate mass region (60 < m ll < 120 GeV), then away from the Z peak the picture is rather similar to the low mass case, for the same reasons. Indeed, across the entire region the initial-state Z contribution is again very small. However, in the Z peak region we observe a significant enhancement with respect to the pure γγ that comes from the mixed γq contribution, where the dilepton pair can now be produced by a resonant intermediate Z (see Figs. 1). This provides the dominant contribution, as we might reasonably expect. We also show the result without applying the subtraction (23), to assess its impact, and can see that this results in a non-negligible reduction in the cross section in the Z peak region, while outside of this the effect is tiny. Broadly, it is clear that even the total cross section, while enhanced with respect to the γγ, still constitutes a very small fraction of the DY cross section, representing as it does as genuine NLO EW correction.
We next consider the high mass (120 < m ll < 5000 GeV) region 2 . We can see that as in [27] the pure γγ SF contribution lies somewhat below the collinear prediction, outside the scale variation uncertainty bands of the latter. The total SF cross section is enhanced with respect to the γγ, in particular at higher masses, due to the impact of initial-state Z contributions, which play more of a role as the average γ/Z virtuality Q 2 increases in this region. Clearly though this is not a question of improving the agreement with the collinear γγ calculation, which explicitly does not include such a contribution.
Finally, we can see that the experimental uncertainty on the SF calculation is very small in the low and intermediate mass regions (though shown, it is rather less than the line thickness), while at the highest masses it is a little larger though still small. In more detail, this varies from ∼ 1 − 3% across the mass region, with the scale variation uncertainty in the γ/Z + q contribution being a negligible component of this.
In all cases the LO collinear prediction tends to lie above the SF result, though at lower and intermediate masses the result is in agreement within the large scale variation uncertainties. However, as we have seen at higher mass even then the result lies above the SF prediction. We now investigate the extent to which this situation improves upon the inclusion of NLO corrections to the collinear calculation.

Collinear predictions: the impact of NLO corrections
In this section we consider the impact of NLO corrections on the collinear calculation for PI production, and compare with the SF calculation. We note that while the fixed-order calculation of the purely QCD contribution to lepton pair production has been performed with very high precision, up to even N 3 LO [47][48][49], for the PI component NLO represents the highest publicly available result. Hence, we will be comparing with the highest precision calculation of PI production in collinear factorization that is currently available.
At NLO in the collinear calculation, as well as the diagrams shown in Figs. 1 (b) and (c), the diagram of the topology (a), but where one of the initial-state photons is produced via a q → qγ splitting, enters. As discussed in Section 2.2, away from the Z peak the (b) and (c) topology diagrams are kinematically suppressed with respect to the t-channel diagram (a), i.e. due to t-channel lepton exchange. Therefore, while these are certainly amenable to being included in a collinear calculation, and indeed can be, for simplicity we choose to focus in this section purely Collinear, LO Collinear, NLO SF Figure 3: Ratio of the LO and NLO collinear predictions to the structure function result, for photon-initiated lepton pair production, as a function of the lepton pair invariant mass, m ll . The leptons are required to lie in the |η l | < 2.5 region. In the collinear case the uncertainty band due to renormalization/factorization scale variation by a factor of two around the central value µ = m ll is shown.
on the topology of diagram (a), along with its corresponding q → qγ NLO contribution. This is of course not strictly the entire NLO contribution, even if away from the Z peak it is by far the dominant one, but for our purposes, i.e. to evaluate the extent to which including NLO corrections stabilises the scale variation of the collinear calculation of diagram (a), we will for clarity only consider this contribution, omitting diagrams (b) and (c) entirely.
In order to achieve this, we make use MadGraph5 aMC@NLO [50,51], suitably modified in order to remove the diagrams of type (b) and (c), as well as those due to final-state photon emission and virtual EW corrections to the γγ → l + l − process, which are not relevant for the current comparison. We in addition correct the results in order to include a running value of α, which is not applied at present in the public version of MadGraph5 aMC@NLO; the collinear photon PDF is extracted under the assumption that one will use a running α and hence this is required to ensure overall consistency, and indeed genuine NLO accuracy.
In Fig. 3 we show the same invariant mass distribution as in the previous section, but now plotting the ratio of the collinear prediction at LO and NLO to the SF result. In the collinear case we show the uncertainty band due to varying the default renormalization/factorization scale µ 0 = m ll by a factor of 2, while in the latter case we only include the pure γγ contribution, i.e. again the t-channel topology of Fig. 1 (a), so that we are comparing completely like-for-like. In the LO case this therefore simply comes from taking the ratio of the corresponding curves in Fig. 2, although we do not include any uncertainty on the SF calculation. This can be read off instead from Fig. 2, though it should be emphasised that this is largely correlated with the corresponding PDF uncertainty in the collinear case, and hence would mostly cancel in the ratios we show here. Moreover, these absolute size of these uncertainties is we recall at the percent Collinear, LO Collinear, NLO SF Figure 4: As in Fig. 3, but for high mass region alone, and with no leptonic loops included in the running of α.
level across most of the mass region. Starting with the low mass region (10 < m ll < 60 GeV), we can see that broadly speaking the inclusion of NLO corrections indeed leads to smaller scale variation uncertainties and to an improved agreement with the SF result. However, the scale variation uncertainty is not negligible: in the lowest mass bin it is over 40% in total (i.e. from the lowest to highest value with respect to the central), decreasing to ∼ 20% in the highest mass bin. This is considerably larger than the PDF uncertainty, and indeed the experimental uncertainty on the SF calculation, which both enter at the level of ∼ 1% in this region. We note that the rather precise level of agreement of the central collinear prediction with the SF is certainly accidental, given the size of the scale variation band; indeed, we will show below that once a cut is imposed on the lepton p ⊥ this precise level of agreement is not maintained. This trend essentially continues for the intermediate mass region (60 < m ll < 120 GeV). That is, we find an improved agreement at NLO, with a reduced scale variation band, but which is nonetheless significantly larger than the PDF uncertainty. The total scale variation band in particular ranges from ∼ 20 − 10% with increasing mass, while the PDF uncertainty is again at the ∼ 1% level. However, it should be noted that as we only include the t-channel topology of Fig. 1 (a), this contribution is of limited phenomenological relevance as we approach the Z peak region.
Turning now to the high mass (120 < m ll < 5000 GeV) region, we observe the same basic trend, however interestingly we find that the NLO collinear result does not agree with the SF within its scale variation band. The same trend was observed in the previous section for the LO result, as is clear here as well. A potential cause of this may lie in the fact that the MMHT2015qed nnlo PDFs we use to evaluate the collinear prediction do not include leptonic contributions in the virtual photon splitting function P γγ . The argument for doing this rests on the fact that lepton PDFs are not included in this set, and hence including such lepton loops in P γγ would lead to a (albeit rather small) violation in the momentum sum rule, see [25] for discussion. However, it is clear from the discussion in [18,27] that this effectively corresponds to evaluating the overall factor of 1/α(µ 2 ) in the definition of the photon PDF, see (14), without leptonic loop contributions to the running of α. There is therefore in principle a mismatch occurring in our results, which as usual do include such leptonic contributions in the running of α. The effect of this is rather small until the highest invariant masses, where the evolution length of α is increased, and the difference is larger, increasing the resultant PI cross section by up to ∼ 10%. To demonstrate this, in Fig. 4 we show the same comparison as before, but now only including quark contributions to the running of α 3 . We can see that the agreement is  indeed improved, although the LO generally still lies outside the SF within scale variation bands, and there are still regions where at NLO the collinear result disagrees within the scale variation bands. In the highest mass bins, this may be due to the increased sensitivity to the z → 1 region of the integral entering the photon PDF, see (14), and in particular in the upper limit of the Q 2 region, which could lead to reduced convergence in the collinear result. In the future, it would be interesting to investigate the comparison using other photon PDF sets [18,26,46] which do including the leptonic contribution to P γγ , and hence should be more consistent with the default running of α. Indeed, as observed in Fig. (24) of [25], the impact of including lepton loops in P γγ is dependent on x as well as scale, and this will not be accounted for in the simple correction applied above. In terms of the size of the scale variation uncertainty, this is at the level of a few percent in this high mass region, and hence is of a similar size to, or even at the highest masses smaller than, the PDF uncertainty. However, as we have seen the difference with respect to the SF baseline is in places larger than that. So far we have taken the dilepton invariant mass, m ll , as the central scale choice. Another possible choice is the lepton transverse momentum, and we therefore investigate the impact of taking this instead. For concreteness we take the scalar average of the lepton transverse momenta, but have checked that for example taking the leading lepton p ⊥ gives a rather similar result, and in fact increases the difference with respect to the SF calculation we observe at NLO. In Fig. 5 Collinear, LO Collinear, NLO SF Figure 6: As in Fig. 3, but for the low and intermediate mass ranges alone, and with a cut on the lepton transverse momenta of p l ⊥ > 5 GeV and 20 GeV imposed, respectively. The leptons are as before required to lie in the |η l | < 2.5 region.The upper (lower) panel corresponds to µ = m ll (p av ⊥ ).
the NLO correction itself rather small with respect to the central LO prediction, the result is systematically lower than the full SF calculation, and well outside these scale variation bands. The reason for this appears to lie in the correlation between the lepton transverse momenta and the transverse momentum of the off-shell initiating photon, due to the q → qγ splitting in the NLO diagram. This can lead to values of the p l ⊥ that are rather higher than the LO kinematic limit p l ⊥ ≤ m ll /2, and this appears to reduce the convergence of the collinear result at lower masses. In the high mass region, the impact of this effect is relatively smaller, and we observe greater stability in the result, and consistency with the µ = m ll case. Thus this scale choice appears to be rather disfavoured. However, it should be emphasised that we can only reach this conclusion after directly performing the calculation in the SF approach, where there is no such ambiguity due to the choice of scale. Without such a comparison, all we could conclude is that the two scale choices lead to rather differing results until fairly high mass, which at NLO largely do not overlap in their scale uncertainty bands.
Finally, we might worry that the above results are all produced with no cut imposed on the lepton p ⊥ , which is certainly not experimentally relevant, and could influence the comparison. We therefore show in Fig. 6 the low and intermediate mass regions again, but now imposing a cut of p l ⊥ > 5 (20) GeV on the lepton transverse momentum in the low (intermediate) cases; at high mass the impact of an equivalent O(10) GeV cut would be less visible, and hence we do not show this here for brevity. The impact of this is relatively moderate, but not negligible. We find in general that the agreement between the NLO collinear and the SF for µ = m ll has deteriorated, with the SF result lying on the edge of the scale uncertainty band or even beyond it. For µ = p av ⊥ the difference with respect to the no p l ⊥ cut case is somewhat smaller. However, in both cases the scale variation band becomes rather larger as one approaches the lower mass region, where the impact of the cuts is more relevant. The agreement between the two scale choices is in addition somewhat improved, although not the agreement with the full SF result.
In summary, we have found that when we take µ = m ll as the central scale choice, the impact of NLO corrections is broadly to improve the matching of the collinear prediction to the full SF result for the t-channel topology of diagram (a) in Fig. 1, with reduced scale variation uncertainties, in comparison to the LO calculation. Nonetheless at low to intermediate masses, m ll 500 GeV, the scale variation is still O(10%) or more, and hence is significantly larger than the PDF uncertainty. At higher masses, m ll 500 GeV, the scale variation uncertainty is rather smaller, at a similar level to the PDF uncertainties, and the matching between the NLO and the SF relatively good, once we account for an apparent mismatch in the treatment of the running of α and the evolution of the MMHT2015qed nnlo PDFs. After correcting for this the agreement is improved, but it is not always within scale variation bands. On the other hand, once we impose realistic cuts on the lepton transverse momenta the agreement in the low to intermediate mass region deteriorates somewhat. Moreover, if we instead take a scale µ = p av ⊥ , given in terms of the lepton p l ⊥ , the agreement between the NLO result and the SF is poor, and well outside the scale variation bands. Therefore, while the stability and the precision of the prediction can be improved by the inclusion of the NLO corrections, this is not guaranteed, and the remaining scale variation uncertainties can remain significantly larger than the percent level PDF uncertainties, as was the case at LO. We re-emphasise here that for this t-channel PI process, which kinematically dominates away from the Z peak, the SF result has no such scale variation uncertainty, or indeed ambiguity due to the choice of dynamical scale, and therefore this issue is bypassed entirely. Indeed, it is only by comparing with the SF result that we can say which scale choice, µ = m ll or µ = p av ⊥ , is preferred. The NLO calculation on the other hand consists in explicitly re-calculating the q → qγ splitting contribution to the PI process, which is already implicitly included in the SF result, but in a manner that introduces a scale ambiguity, and hence source of uncertainty, that is not present in the SF case.

Dilepton p ll ⊥ distribution: revisited
We now consider the impact on the dilepton transverse momentum, p ll ⊥ , distribution. As in [27], we consider the same event selection as the ATLAS 8 TeV measurement [52], which is presented both on and off the Z peak. The leptons are required to have p l ⊥ > 20 GeV and |η l | < 2.4. We now compare to the on-peak region, as well as considering the impact of the γ/Z + q and initial-state Z contributions.
In Fig. 7 we show results in the p ll ⊥ > 30 GeV region, and consider the ratio to NNLO QCD theory, produced with NNLOjet [53]. We do not include EW corrections to the QCD process, but note these are found to reduce the cross section at high p ll ⊥ (see e.g. [54]) and hence will increase the relative contribution from PI production somewhat in this region. In the left figure we show the pure γγ SF predictions, where these contributions are seen to be small, but not necessarily negligible, being at the percent level at lower p ll ⊥ in some mass regions. In the on-peak region the relative contribution is as expected lowest. Now, in the right hand plot we show the total SF cross section, and we can see that this leads to a sizeable increase in the cross section at higher p ll ⊥ , while in the lower p ll ⊥ region the results are relatively stable. This is due both to the γ/Z + q and initial-state Z contributions. The former component plays a significant role due to the ∼ Q 2 i /m 2 ll ∼ (p ll ⊥ ) 2 /m 2 ll scaling discussed in Section 2.2, and is particularly large at lower m ll , for the same reason. The latter component leads to an enhanced at larger p ll ⊥ across all mass regions, due to the larger Q 2 i probed in this region; in the largest p ll ⊥ bin it increases the cross section by a factor of ∼ 2. For the Z peak region, the major contribution is simply the inclusion of the γ/Z + q component, which then allows resonant production of the dilepton system, with the relative contribution being largest at larger p ll ⊥ (though note that the absolute  cross section is of course steeply falling).
In more detail, in Fig. 8 we show the 12 < m ll < 20 GeV mass region of Fig. 7, but now focussing on the relative impact of the γ + q and γγ contributions; we omit the initial-state Z contribution here for clarity. We can see that for the SF results, the trend is as expected, namely at lower p ll ⊥ the γγ contribution dominates, whereas at higher p ll ⊥ the γ + q takes over. The small uncertainty band in the latter case is dominated by the scale variation uncertainty on this mixed contribution due to the initial-state quark/antiquark, with the experimental uncertainty on the pure γγ case being as usual very small, and indeed not visible on the plot if it were included. We in addition however show the γ + q contribution as calculated in the purely collinear LO case, that is via quark/antiquark and photon PDFs; due to the p ll ⊥ > 45 GeV cut that is imposed in this mass region this cross section is free of initial-state collinear singularities. The scale variation uncertainty about µ 0 = m 2 ll + (p ll ⊥ ) 2 is shown, and is significantly larger, due to the introduction of a photon PDF in this case. Interestingly, this lies below the SF prediction, in particular at larger p ll ⊥ . This is as we might expect, given that a non-zero photon q ⊥ is permitted in the SF case, but not in the LO collinear result.
Finally, in Fig. 8 (right) we focus on the low p ll ⊥ region, considering the ratio to NNLO+N 3 LL resummed QCD predictions produced with NNLOjet+RadISH [53]. We consider three mass regions, both and on and off Z peak, and show both the pure γγ (dashed line) and total (solid line) SF predictions. Away from the Z peak, the former component is completely dominant until we approach larger p ll ⊥ values. On the Z peak, as usual the inclusion of the γ + q component allows for resonant dilepton production, and hence a significant enhancement. Of particular interest though is the significant enhancement observed in the pure γγ contribution in the lower 0 < p ll ⊥ < 2 GeV mass region. As discussed in [27], this is explained in part by the Sudakov suppression in the QCD contribution in this region, which is absent in the γγ channel. However, another key factor in this is that the γγ cross section is particularly peaked in this region, due to the significant contribution from elastic photon emission. It is important to point out here that this is an element of the cross section that by construction can never be accounted for within a purely collinear photon PDF framework. By including higher orders in the collinear calculation one can improve the accuracy from the Q 2 continuum region, and in the present case account for the p ll ⊥ dependence of the cross section due to this, which is trivial at LO. However this is not the case for the elastic component, or indeed the low Q 2 resonant and non-resonant components,  Figure 8: Percentage contribution from photon-initiated production to the lepton pair p ⊥ distribution, within the ATLAS [52] off-peak event selection, at 8 TeV. The QCD predictions correspond to NNLO + NNLL QCD theory [53]. The left plot corresponds to the 12 < m ll < 20 GeV mass region, and shows the γγ and γ + q SF contributions, as well as the collinear γ + q. In the collinear case the uncertainty band due to scale variation by a factor of two around the central value µ0 = m ll , is given. The left plots shows the total (solid line) and pure γγ (dashed) contributions for different mass regions. For both plots, in the total SF case, the error band due to the experimental uncertainty on the structure functions and scale variation in the γ + q component is shown.
all of which cannot be modelled differentially in this way, and are the relevant components in the low p ll ⊥ region. Thus the only way to precisely account for this low Q 2 (and hence low p ll ⊥ ) region is via the SF result 4 .
A potentially significant impact of the above discussion is that, while in the inclusive cross section the γγ contribution to the on-peak dilepton cross section is strongly suppressed, this is in part counterbalanced at very low p ll ⊥ by the enhancement in the γγ, and corresponding Sudakov suppression in the QCD DY cross sections. In the lowest mass bin, the γγ contribution is at the ∼ 4 per mille level. Of course this is still rather small, but it is worth recalling that such a level of precision, or even better, is required in the LHC high precision programme. This could in particular be relevant for extractions of the W boson mass, as we will now discuss.

Impact on M W determination
As is well known, while the W boson mass cannot be determined via its leptonic decay directly due to the invisible neutrino in the final state, one can construct kinematic observables that are sensitive to M W . The two most common such variables are the transverse momentum of the charged decay lepton, p l ⊥ , and the W boson transverse mass where ∆φ is the azimuthal opening angle between the charged lepton and the missing transverse momentum. In particular, both of these variables exhibit Jacobian peak structures that are sensitive to the value of M W , and therefore by fitting the measured distributions with a range of templates for differing mass values, one can extract information about M W . Indeed, this was precisely the technique used by the ATLAS 7 TeV analysis [6], for which the corresponding total uncertainty of 19 MeV is comparable to previous Tevatron determinations, and is starting to bear down on indirect constraints. An important aim for the future LHC and HL-LHC programme is to reduce this uncertainty further, see [11].  To achieve this level of precision in a hadron collider environment, a precise and accurate control over the underlying W production process is essential. The impact of PDF uncertainties is particularly important in this context [55][56][57][58][59], but there are numerous other sources, see [60] for a review. One possible source of contamination that has however not so far been considered in detail is the contribution from PI production. At first glance this might appear completely irrelevant, given lν production can of course not occur directly from the γγ initial state. However, the issue here is that to achieve high precision it is necessary to perform a systematic tuning of certain theoretical inputs to Z/γ → l + l − production, as discussed in detail in Section 6 of [6]. One particularly relevant case in point is the tuning of certain QCD parameters in PYTHIA 8 [61], relating e.g. to the intrinsic partonic k ⊥ , to the measured dilepton p ll ⊥ distribution; these are then used as an input in the W mass measurement. Thus any contribution from γγ-initiated production in the Z/γ case, precisely because it is entirely absent in the W case 5 , may bias the results if it is not accounted for correctly, by being tuned away in the comparison with the dilepton data. In the ATLAS analysis [6], tuning is for example performed to a measurement of the Z boson transverse momentum distribution presented in [62]. Here, the PI contributions is subtracted using the collinear LO predictions with the now outdated MRST2004QED PDF set [12].
As discussed in the previous section the contribution from γγ-initiated production in the Z peak region is small, though at low p ll ⊥ as large as ∼ 4 per mille. The precise quantitive impact of this contribution is non-trivial to evaluate, as it will depend on the details of the tuning performed by any given experiment, through which the effect will propagate. Here, to give some idea of the relative size of effects, we will consider the relative contributions from the γγ channel to the p l ⊥ and Z boson m ⊥ distributions. Note the results are identical with respect to which lepton is considered to correspond to the 'charged' one in defining these observables in this channel. This is shown in Fig. 9, which in particular gives the ratio of the DY prediction including γγ-initiated production to the pure Z/γ DY case. The basic idea in these comparisons is that any γγ contribution that is incorrectly tuned away in the comparison to the Z/γ data could appear as a corresponding shift with respect to the description of the W data. For example,  Figure 10: As in Fig. 9, but with collinear MRST2004QED predictions given. The PDF uncertainty corresponds to the difference between the constituent and current quark models, while scale variation uncertainties are omitted for clarity.
if one did not subtract any PI component from the dilepton data, and one assumed that the underlying lnu data corresponded to e.g. the PDG value of M W , then the template at this value could still be shifted by the amount (and in the direction) shown in the figure. We consider the normalized distributions to be the more relevant observables for such a comparison, although the conclusions are rather similar, though generally larger, if the absolute are instead taken.
For the DY predictions we use MCFM 6.8 [63,64] to provide fixed-order NLO QCD predictions.
Given that the precise shapes of the distributions are sensitive to e.g. parton-shower effects, the following results should therefore be taken as a guide only.
To give an idea of the required precision, we in addition show the impact on the distributions of modifying the Z boson mass by 10−20 MeV away from the PDG value [65]; the corresponding modifications in the W boson case follow these very closely after shifting the distributions to account for the different W mass. Considering first the p l ⊥ distribution in the left plot, we can see that impact of γγ-initiated production on the normalized distribution is at the ∼ 1 per mille level or below around the most relevant M Z /2 region, while at higher values it is as large as ∼ 2 per mille. It should be noted though that the higher p l ⊥ region is likely to have a smaller impact on the corresponding tune, and in addition the underlying DY prediction will be more sensitive to effects beyond fixed-order NLO QCD. In any case, the effect is as we might expect rather small. However, we can see that while the effect of changing the input Z mass by 10 − 20 MeV is generally larger, the γγ contribution is not completely negligible in comparison. A similar conclusion is seen for the m ⊥ distribution in the right figure. Again the impact of the γγ-initiated production is small, but not negligible in comparison to the small impact of shifting the Z mass, and arguably somewhat larger than in the p l ⊥ case. In both plots the shape of the γγ induced modification does not follow that of the mass shifts particularly closely, as we might expect. Now, as discussed above for the ATLAS measurement the PI component is already subtracted from the dilepton data [62] and thus this component will certainly not be fully absorbed into the tune. However, this subtraction is done using LO collinear theory, and with the outdated MRST2004QED set; both of these will lead to imprecision in the corresponding modelling of PI production. To assess the latter effect, we show in Fig. 10 the same plots as before, but now with the LO collinear predictions made using this PDF set, and with µ 0 = m ll . The PDF uncertainties alone, due to the difference between the constituent and current quark models, are shown for clarity; the scale variation uncertainty will as usual be present in addition. We can see that these are non-negligible, and as we would expect significantly larger than the SF uncertainties, though the results overlap within uncertainties in most regions. Now, in future analyses this outdated set will not be used, and hence a potentially more instructive comparison is as in Fig. 9, i.e. with respect to the more precise and newer MMHT2015qed nnlo PDF result. However, here again we can see a clear difference between the SF and collinear LO predictions. We can in particular see that in many regions the SF prediction does not overlap with the LO collinear within scale variation bands of the latter prediction. This is not surprising, given the sensitivity of these distributions to variables such as the dilepton p ll ⊥ , which is only accounted for in a non-trivial way by the SF prediction. As discussed in the previous section, one can of course improve the situation by going to higher orders in the collinear calculation, but the result of Fig. 8 (right) suggests that the contribution from elastic and low scale inelastic production may be particularly relevant at this level of precision, and the full kinematic dependence of these components can never be accounted for by a fixed-order collinear calculation. It is in particular worthy of note that in this comparison the difference between the LO collinear and SF predictions is sometimes comparable to the overall expected impact of PI production on these distributions.
In fact, in ATLAS analyses such as [62] the primary vertex is required to be reconstructed from at least three tracks. This will automatically exclude a significant fraction of PI events, in particular the majority of those due to elastic photon emission from the protons as well as a reasonable fraction of events with inelastic photon emission. A very useful feature of the SF calculation is that one readily isolate the contributions from elastic and inelastic photon emission and hence account for the impact of such a cut; elastic events will to first approximation always fail the vertex requirement, while for inelastic events this will depend on the kinematics of the proton dissociation system and its interface to an all-purpose MC generator for showering and hadronization of the system, all of which is accounted for in SFGen. In more detail however, one must account for the probability of no additional particle production from additional proton-proton interactions, independent of the hard scatter; that is, the so-called survival factor. This effect is not included in SFGen, but these issues are discussed in detail in [38], with the SuperChic MC providing an implementation of the SF calculation of semi-exclusive PI lepton pair production relevant for events in the presence of the vertex requirement imposed by ATLAS.
We end this section by again noting that the above results can only serve as a first look at the possible impact of these effects, in particular given the connection between the tuning to dilepton data and the M W determination is not completely direct, and will as mentioned above depend on the details of the tuning and the precise validation performed, as well as the fact that we only use fixed order NLO QCD predictions here. An alternative, and potentially informative approach would be to consider the effect of modifying the W boson p ⊥ distribution according to Fig. 8 (right) directly on the extracted W mass; such modifications at low p ⊥ have been examined in e.g. [66] in the distinct context of intrinsic parton k ⊥ . In any case, our results certainly verify that γγ-initiated production is potentially relevant at the level of precision being aimed for, and in future experimental analyses the SF approach could provide an invaluable tool to account for this.

Triple-differential Drell Yan: photon-initiated contribution
High precision triple-differential measurements of lepton pair production have been reported at 8 TeV by both ATLAS [7] and CMS [8], with important implications both for PDF constraints and the determination of the weak mixing angle, sin 2 θ W . There are in particular given in terms of the dilepton rapidity, y ll , invariant mass m ll , and cos θ * (and/or the forward-backward asymmetry constructed from this), where θ * is the lepton decay angle in the Collins-Soper frame [67]. We can expect PI production to give an important contribution in certain regions of phase space, and therefore for this to have an impact on any determination of sin 2 θ W as well as PDF constraints from the data. While, as discussed in [8], the former is mostly sensitive to the Z peak region, where the PI contribution is smallest, sensitivity is not completely confined to this, and conversely the largest sensitivity to PDFs comes from the off peak regions. Moreover the approach applied in the CMS analysis is for in-situ constraints on the PDFs to be performed across the considered m ll region in order to maximise the precision of the measurement of sin 2 θ W . Thus even if the relative contribution from PI production is quite small on the Z peak, it could still bias such a procedure if not accounted for correctly.
In this section we therefore show some representative results for the PI contribution to triple differential dilepton production, considering the ATLAS 8 TeV event selection for concreteness. We will in particular present results for the percentage contribution from PI production to the NNLO QCD DY prediction, as calculated using NNLOJET [53]. We leave a comparison to data to future work, as this can most meaningfully be done in the context of a global PDF fit, accounting for all sources of correlated error in the data. In all cases, we only show the pure γγ channel for clarity, but have checked that the contribution from Z bosons and the mixed γ/Z + q channel is small. We only show values for those bins where a non-zero measurement is recorded by ATLAS, in order to avoid regions of phase space where the cross section is very highly suppressed and hence the ratio values may be somewhat misleading. In addition, stable NNLO QCD predictions for some of these bins were not always available.
We first note that in general the forward-backward asymmetry in PI production is of course zero; here σ F and σ B are the cross sections in the forward (cos θ * > 0) and backward (cos θ * < 0) hemispheres, respectively. However, the contribution from this to the overall dilepton sample is not, as it will increase the denominator in the above, while leaving the numerator unaffected. In other words, it will tend to wash out the asymmetry. More generally, it will affect the fully differential cross section in cos θ * in a non-trivial way.  Figure 12: As in Fig. 11, but for 150 < m ll < 200 GeV.  Figure 13: As in Fig. 11, but for 66 < m ll < 80 GeV and the ATLAS [7] central-forward lepton measurement.
We will focus on those mass bins where the PI contribution is largest, namely at the lower and upper ends of the considered range, but will comment on the contribution in other bins at the end of the section. We begin by consider the central lepton (muon and electron combined) measurement, that is with both leptons required to lie in the |η l | < 2.4 region, with p l ⊥ > 20 GeV. Results for the 46 < m ll < 66 GeV and 150 < m ll < 200 GeV mass bins are shown in Figs. 11 and 12, respectively. Again, though the PI contribution is symmetric in the cos θ * > 0 and cos θ * < 0 regions, the DY contribution is not, and hence these ratio plots are not symmetric with respect to this. We can see that the PI contribution ranges from ∼ 1 − 6% in all bins, and is therefore small but certainly not negligible. In more detail, there is a clear tendency for the PI contribution to be larger at low rapidities and forward cos θ * values; in the latter case this is driven by the t-channel nature of the PI process. Thus, as expected this will modify the corresponding cos θ * distribution and asymmetry A F B .
We in addition show the LO collinear result, with corresponding scale variation uncertainties. Interestingly, there is a distinct tendency for the collinear prediction to lie above the SF, as seen in Fig. 2 already. However, here the effect is larger and the LO result lies above the SF outside the scale variation uncertainty bands in most bins. We note that the increased size of this difference results from the event selection being different in comparison to Fig. 2 Figure 14: As in Fig. 13, but for 116 < m ll < 150 GeV.
cut on the p ⊥ of the lepton is placed. Indeed, given the LO collinear prediction corresponds in essence to a cross section integrated over the photon virtualities Q 2 , the less inclusive the observable becomes, the larger we might expect any disagreement with the SF result to be, and this is precisely what we see. In some of the higher y ll bins there is no phase space for production with zero dilepton p ll ⊥ , and hence the collinear prediction vanishes entirely. This is in contrast to the SF result, though these are in rather suppressed regions, where moreover the QCD prediction may be less reliable.
We next consider results for the central-forward electron selection, that is with the central electron required to have p l ⊥ > 25 GeV and |η l | < 2.4, while the forward electron is required to have p l ⊥ > 20 GeV and 2.5 < |η l | < 4.9. Results for 66 < m ll < 80 GeV and 116 < m ll < 150 GeV are shown in Figs. 13 and 14, respectively. Although the precise shapes of the distributions is different, the broad picture is very similar to the central channel. Namely, the PI contribution ranges from ∼ 1 − 6%, with some tendency to decrease with y ll , though this effect is certainly milder here. The LO collinear predictions again have large scale variation uncertainties and lie above the SF predictions, often outside of the uncertainty band. As before, in certain bins the LO collinear result can be zero due to kinematics, while not being the case in the SF calculation.
Finally, we comment on the mass bins not shown here, for which the PI contribution is generally smaller. For the central selection, the off-peak bins of 66 < m ll < 80 GeV and 102 < m ll < 116 GeV have PI contributions of ∼ 1 − 2% and ∼ 1 − 4%, respectively. For the two on-peak bins, 80 < m ll < 91 GeV and 91 < m ll < 102 GeV, the contribution is at the 1 per mille level. For the central-forward selection, the off-peak bin 102 < m ll < 116 GeV has a PI contribution of ∼ 1% or less, while for the on-peak bins this is similar to the central case.

Lepton-lepton scattering and the lepton PDF
In this section we consider an application of the SF formalism to a different class of process to that considered above. Namely, one where a lepton produced in a γ * → l + l − splitting can be viewed as an initial-state parton in the production process. This idea has recently been considered in [36], where it shown that by integrating over the initial-state photon kinematics one can recast the process in terms of a collinear lepton PDF within the proton. As with the original analysis of [18,19] for the photon PDF, they demonstrate how this can be defined consistently such that higher order contributions, i.e. those with an explicit γ → l + l − splitting, where the initial-state photon is treated in the collinear LO approach, can be included in the MS scheme. However, given our previous results in the context of PI production in the SF calculation, we may expect to arrive at a rather precise evaluation of this process by instead using the SF calculation directly.
To demonstrate this, we reconsider the example of lepton-lepton scattering taken in [36]. As discussed there, the observation of an isolated and back-to-back same-sign lepton pair of the same flavour, or a lepton pair of differing flavours and arbitrary signs is an important test of flavour violating interactions in various BSM scenarios, see [68][69][70][71]. A representative diagram for the lepton-lepton scattering subprocess within the SF approach is shown in Fig. 15 (a), while the corresponding LO collinear diagram is given in Fig. 15 (b). The calculation of the SF cross section then requires a relatively straightforward application of (1); we will drop the Z and mixed Z/γ + q contributions, which are both negligible. In a little more detail, some care is needed in performing the phase space integration, as there are two independent logarithmic enhancements in the region of collinear γ → l + l − emission from both initial-state legs, which while regulated by the non-zero photon Q 2 as well as the lepton mass, and so perfectly finite, can destabilise the numerical integration. The solution is simply to express the phase space integration directly in terms of the two leptons we do not require to be observed in the detector, and hence which can/will be emitted favourably in a collinear direction to the initial-state photon. Their angular integration is then alway defined with respect to the photon direction, which for off-shell photons is not in general along the z axis. Once this is done, standard adaptive integration tools and/or a suitable change of variables can stabilise this logarithmic enhancement.
We consider precisely the same observables and event selection as in [36], in order to perform a direct comparison. In particular, the signal leptons are required to have p l ⊥ > 20 GeV and |η l | < 2.4 and results at 13 and 27 TeV are shown. The cross section results are shown in Table 1 for both the SF calculation and the LO collinear case, as given in [36]. Errors due to the experimental uncertainty on the structure function are given in the SF case, and factorization scale variation in the collinear case. We can see that the LO collinear predictions carry significant errors due to this, ranging from ∼ 30 − 100% of the cross section, which are absent in the SF case. The errors on the latter are as expected very small, at the percent level. The SF results lie somewhat below the central LO collinear predictions, but within scale variation uncertainties in all cases. In terms of the scaling with the mass of the leptons in the final state, these are closely matched in the SF and LO calculations, e.g. the ratio of the e + µ − to e + τ − cross section is rather similar in the two cases. As pointed out in [36], a NLO calculation can certainly be performed, which will result in smaller scale variations in the collinear case and, we would expect, a somewhat closer agreement with the SF prediction.     Table 3: Cross sections (in fb) and acceptance corrections (as defined in Table 2 and the text) for same-sign same flavour and opposite-sign different-flavour lepton pairs in the SM, at 13 TeV. Results within the SF approach and for the differing elastic (el.), single dissociation (sd) and double dissociation (dd) cases are given.
Next, we consider the impact of the cuts which as described in [36] can strongly suppress the dominant SM contribution from same-sign W pair plus dijet production. While these leave the LO collinear prediction trivially unchanged, by applying the SF calculation we can get a precise evaluation of their impact, shown in Table 2.
No errors from the experimental uncertainty on the structure functions are given, as these will be tiny. We can see that impact of the cuts is expected to be significant, decreasing the cross section by a factor of ∼ 3 − 10, depending on the process. This is not necessarily surprising, as the pure back-to-back configuration of the LO collinear prediction will always be modified by the non-zero photon Q 2 in the full SF calculation as well the non-zero opening angle in the γ * → l + l − transition. We find that the impact of the cuts is larger for the heavier leptons, as one might expect from the fact that the collinear enhancement in the γ * → l + l − transition is less strong in this case. There is a very mild decrease in the acceptance at the larger √ s = 27 TeV value, perhaps due to the increased phase space in photon Q 2 at the lower photon x value this corresponds to.
To analyse these effects in a little more detail, in Table 3 we show the acceptance (and corresponding cross sections values) for e + µ − and τ + τ + production, but now breaking things down in terms of whether the initial-state photons are emitted elastically or not. We can see that even for purely elastic emission, where the average photon Q 2 is rather low, only 50% of events survive the cuts. As we expect, when one proton dissociates ('sd') this acceptance reduces, and when both protons dissociate ('dd') this reduces further still. The combination of these channels, with their respective acceptances, gives the corresponding values in Table 2.
The above results again highlight a useful feature of the SF calculation, namely that one can readily isolate the contribution from elastic and inelastic photon emission. Indeed, this can provide a useful handle in searches for processes like lepton pair production discussed above. In particular, one can use the fact that additional hadronic activity associated with the event is strongly suppressed for the SM lepton-lepton scattering process, whereas it might not be for a BSM signal, or for other SM contributions to the final state. However, a precise account of this requires that we include the effect of multi-parton interactions, i.e. the so-called soft survival factor, if we require no additional hadronic activity, as well as the fact that the proton dissociation products may be present in the detector for inelastic photon emission. These issues are discussed in detail in [38], where an approach that accounts for these effects in (opposite sign, same flavour) lepton pair production is presented and implemented in the SuperChic 4 MC. Such an approach could straightforwardly be applied to the current scenario, but we leave that for future work.
We note that in the above calculations there is a further possible contribution from γγ → l + 1 l − 1 production, with a subsequent γ → l + 2 l − 2 emission from the lepton type 1 line. For example, for µ + µ + production we would have l 1 = l 2 = µ. Considering this example in more detail, the µ + from the original γγ → l + 1 l − 1 would have to produced at sufficiently high p ⊥ to pass the experimental cuts, and only the µ − emission would receive a collinear enhancement, being dominantly emitted at low p ⊥ . Then, the dominant contribution from the γ → l + 2 l − 2 would be due to collinear emission. To deal with this contribution one would therefore need to impose suitable isolation requirements on the signal leptons. This would remove this collinear enhancement but potentially leave this as a genuine NLO correction. However, once the back-to-back requirements (27) are imposed such a contribution should be very small. In particular, and even absent any lepton isolation conditions, these will only be passed if the lepton l − 2 (= µ − ) is at very low p ⊥ O(GeV) in order that the two signal leptons are back-to-back. This is a rather stringent requirement that should reduce such a contribution significantly below the level one would naively expect from a NLO contribution. We therefore do not include this contribution in the calculation, though it could certainly be included within the SF approach.
Finally, we note that the experimental uncertainties on the structure functions are not the only source of uncertainty in the SF prediction. As in the case of lepton pair production, we should also in general include the uncertainty from non-factorizeable corrections, as well as higher order uncertainties in the pQCD calculation of the structure functions, which as before will enter at the percent level. We in addition should consider the impact of higher order QED (as well as in principle EW/QCD) corrections to the underlying γγ to four lepton process. This raises a possible shortcoming of the SF approach, namely that photon emission from the leptons initiating the lepton-lepton scattering process receives a collinear logarithmic enhancement. This will be resummed via the DGLAP evolution of the lepton PDFs in the collinear calculation, but this would not happen in the SF calculation outlined above.
Such a correction is discussed in detail in [36], Appendix D, where the impact of allowing a single photon emission on the derived lepton PDFs is considered. Even up to rather large scales, µ F ∼ 1 TeV, the effect is small, at the percent level or less apart from at rather high x values. This level of difference, driven by the small size of the QED coupling, is not indicative of a strong degree of logarithmic enhancement and with it the necessity of performing a resummation. Though the impact is larger at high x, here the experimental uncertainty on the structure functions is larger as well, so the phenomenological impact (and indeed relevance) appears to be relatively limited. Moreover, one can see in [36] that there is a ∼ percent level uncertainty in precisely how one assigns such a logarithmically enhanced contribution to the lepton PDF, which is smaller than, but not negligible in comparison to, the overall size of the effect. It is certainly not clear that this is a smaller difference than that between only including one photon emission explicitly, say, in the SF approach, and resumming these contributions, as is done in the DGLAP evolution of the lepton PDFs.
Indeed, one can of course include an additional photon emission explicitly in the SF calculation if the desired level of precision requires it, and in principle one could even resum the effect of such emissions directly in the SF calculation. This would arguably be more cumbersome than achieving the desired resummation in the collinear calculation, i.e. via the DGLAP evolution of the lepton PDFs. However, as the impact of such effect is in any case rather small, and certainly significantly smaller than the factorization scale variations in the LO collinear result, and quite possibly the NLO one, it is almost certainly a question of relatively marginal phenomenological relevance for LHC physics. This is particularly true for the comparisons presented above, where the relevant scale of the hard process is of order p cut l ⊥ = 20 GeV, which is rather low.

SFGen Monte-Carlo: availability
The SFGen MC provides a publicly available tool for lepton pair production and lepton-lepton scattering within the SF approach, including initial-state Z and mixed Z/γ + q contributions. Arbitrary distributions and unweighted events can be generated, and errors due to the experimental uncertainty on the structure function (the equivalent of PDF uncertainties in the photon PDF framework) can be calculated on-the-fly. Unweighted events an be interfaced to Pythia for showering/hadronization of the proton dissociation system, following the approach discussed in [38]. The LO collinear PI predictions can also be calculated for comparison, for arbitrary factorization scales. For the Q 2 > 1 GeV 2 continuum contribution, the code makes use of structure function grids in the LHAPDF [72] format, as calculated using APFEL [43], and first used in [73]. These structure function grids are provided with the MC, and are calculated using the approach discussed in Section 2.1, namely with MMHT2015qed nnlo PDFs [25] in the high Q 2 region. The code and user manual is available at: http://projects.hepforge.org/SFGen

Summary and Outlook
The LHC and its high luminosity upgrade will provide a significant opportunity for stress testing the SM at an unprecedented level of precision. Events with leptons in the final state play a particularly important role in this, enabling for example the determination of the weak mixing angle, sin 2 θ W , and the W boson mass, M W , as well as providing constraints on PDFs. The availability of precise theoretical calculations for SM processes is essential in achieving this, and a key element of this is the contribution from photons in the initial state, that is so-called photon-initiated (PI) production. The study of [27] showed for the first time how the PI production cross section of lepton pairs, away from the Z peak, could be straightforwardly calculated with percent-level precision, by applying the structure function (SF) approach. In this paper, we have explored in more detail the phenomenological implications of this approach for LHC processes with leptons in the final state. We have calculated the PI contribution to multi-differential dilepton production, which is relevant for determinations of sin 2 θ W as well as PDF constraints, and evaluated the potential impact on W mass determinations, through the tuning that is performed to dilepton events. In both cases, a precise account of PI production is found to be highly relevant, in particular given the level of precision being aimed for at the LHC. We have shown how the direct SF calculation automatically provides this. A key element in evaluating this contribution is control over all phase space regions, in particular at low dilepton p ll ⊥ , which cannot be accessed via a purely collinear calculation in terms of a photon PDF, and the SF calculation is therefore uniquely positioned to evaluate.
We have presented a detailed comparison of the SF result to the collinear calculation at the highest available order, namely NLO. We find that broadly an improved agreement with the SF result can be achieved with respect to the LO calculation, with correspondingly smaller scale variation uncertainties. However, these uncertainties (which are absent in the SF result) remain significantly larger than the percent-level PDF uncertainties, and related uncertainty due the experimental determination of the structure functions in the SF result, in many regions. Moreover, the agreement between the NLO prediction and the SF is not always within these scale variation bands, in particular for certain choices of scale. These issues are bypassed for pure γγ PI process, i.e. due to t-channel lepton exchange, which dominates away from the Z peak, in the SF approach.
We have also extended the previous calculation to include contributions from Z bosons in the initial state and mixed Z/γ + q contributions. In addition, the PI production of same-sign same-flavour, or different-flavours arbitrary-sign lepton pairs is calculated. We have found that the SF approach can indeed provide high precision predictions for this process, and the result of applying cuts in order to isolate back-to-back lepton topologies is readily evaluated.
On the other hand, while the benefit of applying the SF approach is clear for process where the final state of interest (l + l − , W + W − ...) is directly produced by the γγ initial state, this is not as clear for the mixed γ + q contributions. Here, one must deal with the collinear enhancement of the γ → qq splitting, and indeed further QCD splittings at higher orders, and at this level of precision include QED DGLAP evolution of the quark/antiquark PDFs. This certainly requires one introduce a photon PDF within the LUXqed approach. Nonetheless, one may systematically account for this while still working within the SF approach. We in particular discuss the complication that arises due to double counting with the PI contribution that is already effectively included via the QED DGLAP evolution of the quark PDFs, and provide a simple procedure for removing this. This however ties the level of precision for the SF calculation to that in the corresponding QED DGLAP evolution. It would interesting in the future to examine the extent to which the two calculations agree or differ for these processes.
Finally, we have developed the SFGen MC, a publicly available tool for lepton pair and lepton-lepton scattering within the SF approach, including initial-state Z and mixed Z/γ + q contributions. Arbitrary distributions and unweighted events can be generated, and errors due to the experimental uncertainty on the structure function (the equivalent of PDF uncertainties in the photon PDF framework) can be calculated on-the-fly.
In summary, the PI mechanism provides a moderate but important contribution to lepton pair production. Given the high precision being aimed for at the LHC, a precise account of it is therefore mandatory. We have shown that the SF approach provides a way to achieve this, and demonstrated how this can be extended to provide high precision predictions for other processes involving leptons, such as lepton-lepton scattering. To enable the application of this approach to calculations of such processes, we have provided a publicly available MC tool for use by the community.