Photon and Leptons induced processes at the LHC

We study a few basic photon- and lepton-initiated processes at the LHC which can be computed using the recently developed photon and lepton parton densities. First, we consider the production of a massive scalar particle initiated by lepton-antilepton annihilation and photon-photon fusion as representative examples of searches of exotic particles. Then we study lepton-lepton scattering, since this Standard-Model process may be observable at the LHC. We examine these processes at leading and next-to-leading order and, using the POWHEG method, we match our calculations to parton shower programs that implement the required lepton or photon initial-states. We assess the typical size of cross-sections and their uncertainties and discuss the preferred choices for the factorization scale. These processes can also be computed starting directly from the lepto-production hadronic tensor, leading to a result where some collinear-enhanced QED corrections are missing, but all strong corrections are included. Thus, we are in the unique position to perform a comparison of results obtained via the factorization approach to a calculation that does not have strong corrections. This is particularly relevant in the case of lepton-scattering, that is more abundant at lower energies where it is affected by larger strong corrections. We thus compute this process also with the hadronic-tensor method, and compare the results with those obtained with POWHEG. Finally, for some lepton-lepton scattering processes, we compare the size of the signal to the main quark-induced background, which is double Drell-Yan production, and outline a preliminary search strategy to enhance the signal to background ratio.


Introduction
In recent years, precise determinations of the proton Parton-Density Functions (PDFs) of the photon [1,2] (from here on the LUX papers) and leptons [12] (from here on the LUXLep paper) have appeared. These determinations (inspired by Ref. [4]) were based upon the comparison of the computation of a reference photon-or lepton-initiated process (from here on, /γ-initiated process) in the collinear factorization approach, and in terms of the proton electromagnetic hadronic tensor where p denotes the proton momentum, the bra and ket states refer to a single proton state, and we imply an average over its spin. The hadronic tensor can be expressed in terms of the electroproduction structure functions and of the proton electromagnetic form factors. We will denote this method for the computation of /γ-initiated processes in hadronic collisions as the HT (for Hadronic Tensor) method. The calculation in the factorization approach is analogous to the Weizsäcker-Williams approximation in charged particle scattering. It starts with the simple leading order (LO) approach, where one only needs the Born matrix elements for the /γ-initiated process, and is easily extended to next-to-leading order (NLO) and even to higher orders, depending upon the accuracy of the available PDF. It can benefit from available packages for the computation of NLO corrections, and also for implementing NLO calculations interfaced to parton showers (NLO+PS from here on).
So far, most calculations of /γ-initiated processes using the factorization approach have been performed at the LO level, with the noticeable exception of Ref. [5], where the full NLO corrections to the resonant s-channel leptoquark production was presented. On the other hand, the photon and lepton PDF determinations presently available have NLO accuracy. Thus it would be desirable to fully exploit them using NLO calculations. As a related point, while we have considerable experience regarding the relative size of NLO corrections in processes initiated by quarks and gluons, little is known for /γ-initiated processes in hadronic collisions. In the present paper we discuss some of these processes, compute them at LO and NLO, interface the calculations to parton shower generators according to the POWHEG method [6][7][8], and study their scale uncertainties.
We remark that NLO corrections in /γ-initiated processes arise in a rather peculiar way. Rather than adding gluons to the process, as perhaps one would naively expect, corrections that are formally of NLO importance arise when the incoming photons or leptons are resolved as radiation from a parton that has higher abundance in the proton. For example, if we have an incoming photon, we also have a diagram where the photon arises from quark radiation. This diagram leads to a contribution proportional to a quark density times one power of the electromagnetic coupling α. Being a subtracted NLO contribution it does not have any logarithmic enhancement. On the other hand, the corresponding Born process multiplies a photon density, that is of the order of a quark density times one power of the electromagnetic coupling, times a collinear logarithm L = log(Q/Λ), where Λ is a typical hadronic scale. Such logarithm is parametrically of the same size of the inverse of the strong coupling constant α s (Q), so that indeed the NLO correction is suppressed by a power of α s (Q) with respect to the Born term.
An interesting aspect of /γ-initiated processes is the possibility to search for exotic signals, typically resonances, that may preferably or exclusively couple to leptons or photons. We thus consider the production of a scalar via the processes l + l − → φ and γγ → φ, study their NLO corrections and implement corresponding NLO+PS generators using the POWHEG method. In the case of the γγ → φ process, we can interface our POWHEG generator to both Pythia8 [9] and Herwig7 [10]. Incoming leptons are not handled at the moment by Pythia8, while there is a preliminary implementation in Herwig7, so, in the lepton case, we restrict ourselves to consider the latter. We can thus assess the impact of NLO corrections both in the fixed order calculation, and in the full NLO+PS implementation. It is also interesting to compare the full NLO+PS results to LO+PS ones, since often beyond the Standard Model (BSM) studies are carried out at this level.
In the LUXLep paper, it was pointed out that some Standard Model (SM) leptonscattering processes leading to rare final states, like leptons pairs of different flavours or same sign, may actually be detectable at the LHC, and some LO cross sections for these processes were given. Here we examine this problem in more depth, by including NLO corrections and implementing an NLO+PS generator interfaced to Herwig7. The same processes can also be computed directly adopting the HT approach. 1 We are thus in a position to compare the full NLO calculation with the HT result, thereby validating the NLO+PS approach and also the Monte Carlo implementation of lepton initiated processes.
The factorization approach is based upon an extraction of the /γ parton densities from the result of an HT calculation, neglecting terms that are beyond the required accuracy (i.e. typically beyond NLO). Thus, the HT method could be in principle more precise. However it is more cumbersome to implement and, in the case of incoming leptons, it neglects electromagnetic corrections that are included with the factorization method, namely radiation of collinear photons from the leptons entering the scattering process. Often when considering mixed electromagnetic and strong corrections, two powers of α s are considered to be equivalent to a single power of the electromagnetic coupling constant. According to this counting, initial state electromagnetic radiation from incoming leptons should count as αL ∼ α/α s (Q), i.e. as NLO contributions. We are however in a position to assess the size of these effects, by generating a lepton PDF set where initial-state radiation from leptons has been turned off, and thus we can carry out a meaningful comparison of the two methods.
The paper is divided as follows. We start by discussing in Sec. 2 the production of a scalar particle either from two initial state photons, or from a lepton-antilepton pair. We discuss the size of the cross-sections and their scale uncertainties in Sec. 2.1, and present their differential distributions in Sec. 2.2 and 2.3. In Sec. 3 we discuss lepton scattering processes. We restrict ourselves to those final states that do not receive contributions from the overwhelming Drell Yan background at leading order. We discuss the corresponding production cross-sections with their uncertainty bands at LO and NLO in Sec. 3.2. We then discuss in Sec. 3.3 the cuts required to suppress large SM backgrounds and their impact on the signal. We also present a comparison to predictions obtained with the HT method. We give our conclusions in Sec. 4. Appendix A provides details about the calculation in the HT approach, including basic formulae (A.1) and the phase space (A.2). In Appendix B we describe how we include lepton mass effects to cure the final-state collinear divergences arising from photon to lepton splitting contributions in the NLO calculation.
2 γγ → φ and l + l − → φ In order to better understand the impact of the NLO QCD corrections to /γ-initiated processes we focus on a simple case, namely the production of a neutral scalar resonance which couples only either to photons or to leptons. In particular, the objective of this section is to assess the size of the NLO corrections and their dependence upon the scale choice. Furthermore, we compare selected kinematic distributions at LO and NLO matched with parton shower using the POWHEG method. We describe the interaction among photons and the scalar resonance through the effective 1 A result in this framework has been given in Ref. [11]. In the present work we include contributions that were omitted in [11], as will be discussed in the introduction of section 3.

Lagrangian
where F µν is the electromagnetic field strength tensor, M φ is the mass of the scalar resonance and g γ is a dimensionless effective coupling. For the leptons, we consider a simple effective Yukawa interaction L ( ) Our processes proceed at LO via the partonic sub-processes γγ → φ, in one case, and l ± l ∓ → φ, in the other. According to the discussion given in the introduction, they receive NLO QCD corrections from the partonic processes qγ → qφ (and the corresponding crossed one qq → γφ, which is finite) and l ± γ → l ± φ respectively. In fact, while on one hand the latter processes are down by a power of α, they are enhanced by the ratio f q /f γ ∼ 1/(αL) and f γ /f l ∼ αL/(αL) 2 = 1/(αL) respectively, resulting in an order 1/L ∼ α s correction (for a more detailed discussion see the LUXLep paper).

Total cross sections and scale dependence
In order to understand the relative impact of the NLO corrections and the remaining uncertainties associated to missing higher orders, we show in Figures  γγ → ϕ, pp@LHC13 LO NLO Figure 1: Comparison of LO and NLO cross sections for the photon-induced process in proton-proton collisions at √ s = 13 TeV, normalised to the reference NLO cross section evaluated at µ F = M φ , as a function of the mass of the scalar particle. The three panels correspond to the different choices of the central scale µ 0 = (1, 1/2, 2)M φ . The uncertainty bands are associated to the 3-point variation of the factorization scale in the nominal range 1/2 < µ F /µ 0 < 2.
cross sections for the γ-initiated process in proton-proton collisions at √ s = 13 TeV, as a function of the scalar mass and for different choices of the central factorization scale, normalised to the NLO one computed at a fixed reference scale µ F = M φ . The analogous plot for the -initiated process is reported in Fig. 2. Here and in the following, to obtain our results we have used the PDF set LUXlep-NNPDF31 nlo as 0118 luxqed of Ref. [12], as implemented in LHAPDF [13]. The bands display the associated scale uncertainty cor- 100 1000 σ(µ F )/σ NLO (µ F = M ϕ ) 10 100 1000 10 100 1000 e + e -→ ϕ, pp@LHC13 LO NLO Figure 2: Comparison of LO and NLO cross sections for the lepton-induced process in proton-proton collisions at √ s = 13 TeV, normalised to the reference NLO cross section evaluated at µ F = M φ , as a function of the mass of the scalar particle. The three panels correspond to the different choices of the central scale µ 0 = (1, 1/2, 1/4)M φ . The uncertainty bands are associated to the 3-point variation of the factorization scale in the nominal range 1/2 < µ F /µ 0 < 2.
First, we observe that the LO result is affected by sizeable scale uncertainties which display a qualitatively similar behavior in all cases: they range from about 50% at M φ = 10 GeV to about 5% (20%) at M φ = 2 TeV for the photon (lepton) induced process. The scale uncertainties shrink considerably going from LO to NLO, being of few percent (10%) for the photon (lepton) induced process, with a mild dependence upon M φ . Furthermore, we observe that the different scale-range choices give rather consistent results at NLO. On the other hand, we see that the LO predictions are only indicative of the cross section within a factor of order one, as is typical for LO predictions for collider processes. We find that the choice of the mass of the scalar particle as central scale is appropriate for the photon induced process, leading to K-factors around unity in the whole mass range. The situation is different for the lepton induced process, where this choice leads to large and negative corrections (K-factors ∼ 1/2). This behavior has also been observed in the leptoquark production case [5]. We also observe that in this case for the central scale µ 0 = M φ the LO and NLO bands do not overlap. An optimal choice is provided by a smaller central scale, namely µ 0 = M φ /4. In Tab. 1, we report a collection of LO and NLO cross sections at LHC at 13 TeV for different masses of the scalar resonance, adopting the above optimal choice for the scales, and factoring out the corresponding g γ and y couplings. We observe that photon induced cross sections at a given mass are about a factor 10 4 larger than the corresponding lepton induced ones. This reflects directly the relative abundance of photons and leptons in the proton, since the two LO matrix elements have a constant ratio of order 1.  Table 1: LO and NLO cross sections for a heavy scalar resonance production in /γinitiated processes at the LHC at 13 TeV as a function of the mass M φ . The uncertainties are obtained by performing the 3-point variation around the optimal scale choices µ F = M φ for the γ-induced and µ F = M φ /4 for the -induced process, as discussed in the text.

Differential distributions for γγ → φ
We now examine differential distributions for the photon-induced process at the level of full NLO+PS generation, that we computed using the POWHEG framework interfaced with Pythia 8.245, and Herwig 7.2 (both these generators implement processes initiated by photons [14]). In our analysis we do not include multi-parton interactions, while we let the events shower and hadronize. We use Pythia8 with default parameters. As for Herwig7, we needed to set it up with customized options 2 in order to obtain a reasonable behavior in the region of small transverse momentum of the scalar particle. In particular, we force the assignment of an intrinsic transverse momentum, whose default value is 2.2 GeV, also to colourless particles, which was not the case with the default Herwig7. Another important aspect, especially relevant when performing the matching at NLO, concerns the recoil scheme (see for example Ref. [15]). Most of the NLO events are qγ → qφ. We remind that Pythia8 adopts by default a global recoil scheme for initial state radiation (ISR), that corresponds to an initial-initial (II) dipole setup and affects all final state particles. On the other hand, Herwig7 implements a strict dipole approach, that allows for both II and initial-final (IF) recoils. We observe that most of the time, Herwig7 treats the initial and final quark as an IF dipole, enforcing a local recoil that does not change the transverse momentum of the scalar. This leads to some distortions in the spectrum of the resonance at small p T . For these reasons, we force a global recoil scheme also in Herwig7.
In the following, we fix the mass of the scalar to M φ = 500 GeV and adopt the optimal scale choice µ 0 = M φ . We show results for the transverse momentum of the resonance and of the leading and second jets. In Figs. 3 and 4 we compare leading order Pythia8 and Herwig7 predictions. Notice that at this order, the selected observables are entirely  generated by the shower. We observe that Pythia8 tends to produce harder spectra with their peak shifted to lower p T -values with respect to Herwig7. Overall, the differences reach about 50% in both the low p T and high p T regions. In the former, this can be attributed to the differences among the shower and hadronization approaches, and in the latter to the different modeling of hard radiation away from the collinear region in the two generators. The same results including NLO corrections are displayed in Figs. 5 and 6. We observe a better agreement in the tails of the p T distributions of the resonance and of the leading jet. This is expected since these variables are now modeled by the real corrections in the fixed-order computation. Differences persist in the tail of the spectrum of the second jet, although they are milder. The shapes at small p T values are qualitatively similar to the pure shower results of Figs. 3 and 4.
PY8/HW7 Figure 5: As in Fig. 3 computed at the NLO+PS level, comparing the prediction obtained with POWHEG+Pythia8 (blue) with POWHEG+Herwig7 (red). In Fig. 7 and 8, we show a comparison between NLO+PS and LO+PS results using Pythia8. We have obtained similar results also using Herwig7. We see that the shower  results provide a good qualitative description of all the considered spectra, with rather flat K-factors for moderate to high p T -values. They tend to give harder distributions, overshooting the NLO results by a factor of almost two for transverse momenta larger than about 5 GeV. Summarizing, we find considerable differences in the transverse momentum distribution of the scalar when comparing Pythia8 and Herwig7 at leading order, the Pythia8 spectrum being considerably harder. This feature is mitigated when NLO corrections are introduced, and in fact very good agreement is found both in the large transverse momentum distribution of the scalar and of the two hardest jets. On the other hand, at very low transverse momenta (i.e. p T 10 GeV) the two shower generators yield a very different modeling of the shape of the distributions. These features are not very relevant when searching for heavy resonances, however they give an indication that a more accurate tuning of the shower generators is needed for the modeling of /γ-initiated processes. We will return to this point after the discussion of the lepton-scattering processes.

Differential distributions for
For this process we can only use Herwig7 since Pythia8 does not currently support initialstate leptons. In Fig. 9 and 10, we show the transverse momentum of the scalar resonance and of the leading lepton, comparing the pure shower results with the NLO+PS ones using Herwig7. In this case, we find a qualitatively good agreement among them. At variance with the photon initiated process, the NLO corrections lead to a harder spectrum in both distributions. Above 20 GeV, we observe a rather flat K-factor of about 1.15. In the : Transverse momentum spectrum of the scalar particle in the -initiated production process at the 13 TeV LHC for M φ = 500 GeV. The right panel shows the distributions at low p T values. We compare the pure Herwig7 parton shower prediction (red) with the one obtained after the matching with the NLO computation (blue).
region of small values of the transverse momentum of the leading lepton, we observe a discontinuity in the NLO+PS prediction related to the minimum p T allowed in the POWHEG generation. Again, this is not a crucial problem in the simulation of the production of high-mass states. Nonetheless, it shows the need for further investigation in the modeling of /γ-initiated processes by the showers. The shower allows us to investigate the associated hadronic activity in the event. In Fig. 11, we compare the transverse momentum of the leading leptons (hardest and next-to-hardest) with that of the leading jet. Notice that at LO+PS all the emissions are provided by the shower, while at NLO+PS the hardest emission (leading lepton) is generated according to the exact matrix element. The second lepton and leading jet, being next-to-hardest emissions, display a much softer spectrum with respect to the leading lepton, and are comparable. The shower results present a slightly harder p t jet spectrum, while the differences are much milder in the NLO+PS computation. The overall picture of the next-to-hardest radiation is consistent with backward evolution starting from a photon and a lepton leg. The probability to have a q → qγ splitting for the first and a γ → ll for the second are in fact parametrically of the same order. e + e -→ ϕ pp@LHC13 l 1 l 2 j 1 Figure 11: Transverse momentum spectrum of the leading lepton, second lepton and leading jet in the -initiated scalar production process at the 13 TeV LHC for M φ = 500 GeV as obtained with Herwig7 (left panel) and after matching at NLO (right panel).

Lepton scattering
We are interested in lepton scattering processes such that the dilepton signal does not have excessive competition from other more abundant Standard Model (SM) processes. As a minimum requirement, we must focus upon final states not affected by the very large Drell-Yan background. Thus, we consider lepton pairs of different flavours or with the same electric charge. In the following we will refer to these signatures as non-Drell-Yan (NDY) pairs. We have implemented the NDY processes in the POWHEG-BOX framework, in such a way that, besides being able to compute the NLO cross sections we can also generate events to be interfaced with a shower Monte Carlo program. As we will see in the following, the NDY processes may be observable at the LHC as long as the transverse momenta of the leptons are not too large. Thus both higher-order effects and shower details can significantly affect the result. Moreover, lepton-initiated processes have become available in shower programs only very recently, and their implementation is still in a preliminary stage. For these reasons we have also implemented a direct calculation of four lepton production using the HT (Hadronic Tensor) method, in order to better assess the reliability of our approach, and to provide a guide for the tuning of the Monte Carlos. The HT calculation is carried out by first computing the matrix element for the production of two lepton pairs, initiated by two virtual photons. This matrix element is contracted with the hadronic tensors of each incoming proton, and a full phase space integration is performed. This calculation is described in detail in appendix A.
We stress that in our HT calculation we include all contributions to the photon-induced four-leptons production, at order O(α 4 ). This is in contrast with the calculation of Ref. [11], where diagrams of this order with the two incoming photons attached to the same fermion line are instead neglected. These contributions, not included in Ref. [11], enter already in our NLO calculation (see Fig. 13, (c) and (d)) and as discussed later, are not completely negligible. Furthermore, interference terms in the case of identical fermions are never included in Ref. [11].
The HT calculation includes all higher-order QCD corrections to the NDY processes. On the other hand it does not include the effect of initial-state collinear radiation of photons from leptons. These effects are of the order of the electromagnetic coupling constant α multiplied by a collinear logarithm L = log p T /Λ, i.e. of order α/α s , that in our counting (that is α ≈ α 2 s ) should be considered equivalent to an NLO correction. Rather than trying to add these effects to the HT calculation, we have estimated their size by switching them off in the calculation using the lepton PDFs by suitably adapting the Hoppet [16] routines used in the LUXLep paper. We have found that at the transverse momenta that we have considered these effects are below the percent level, and thus we can trust the HT result for a validation of our calculations. 3 We stress however that the HT calculation is in general much more demanding than the POWHEG NLO one. Furthermore, it cannot be straightforwardly interfaced to a shower program, and thus we have not attempted to turn it into a full event generator. 4

The NLO POWHEG calculation
In the left part of Fig. 12, we display one of the Feynman graph for the process of samesign, different-flavours scattering at the Born level. We also consider same-sign identical- Figure 12: Left: the Born level graph for same-sign, different-flavours lepton scattering; right: an example of a real graph associated with radiative corrections to this Born graph. flavours and opposite-sign different-flavours. In the case of identical flavours an exchange graph must also be included. The NLO corrections to these processes arise from the graphs where one incoming lepton originates from an incoming photon, as shown in the right part of Fig. 12. The real graphs for the NLO corrections have several singular regions, that are shown schematically in Fig. 13, where, for the sake of definiteness, we focus upon the Figure 13: Singular regions associated with real graphs in the lepton scattering process.
scattering of leptons with different flavour. The singular regions are (a) A final-state fermion with flavour different from the incoming one becomes collinear to the incoming photon. The underlying Born for this singular region corresponds to a process of different-flavours lepton scattering, i.e. the same process that we want to include at the Born level. The corresponding collinear singularities are removed following the usual factorization procedure.
(b) A final-state fermion of the same flavour of the initial one becomes collinear to it. This singular region yields a final state that is Drell-Yan like. However, this graph does contribute at the NLO level to our signal when far away from the singular region, so it must be included. Its underlying Born configuration is γγ → l + l − .
(c) A final-state fermion with flavour equal to the incoming one becomes collinear to the incoming photon. Also in this case, the underlying Born configuration is Drell-Yan like, but away from the singular region it may contribute at the NLO level to our signal, and must be included.
(d) Two opposite-flavour final-state leptons become collinear to each other. This process can generate a different-flavour pair that is not Drell-Yan like. Yet, it cannot be ascribed to a lepton scattering process, since its underlying Born is the Compton process γl → γl. We must however include it since it is relevant at the NLO level. In order to suppress this contribution we will impose an isolation criterion upon our final leptons.
Configuration (a) is dealt with automatically by the POWHEG-BOX framework.
In order to deal with the remaining singular regions, we introduce also the Born subprocesses γγ → l + l − , l + l − → + − and lγ → lγ. In this way, regions (b), (c), and (d) are correctly detected as singular ones. We then use the Bornzerodamp feature of the POWHEG-BOX [17] to separate those regions as remnants and treat them as pure unsubtracted 5 real contributions. We notice that the divergent contributions of regions (b) and (c) are eventually excluded by the analysis cuts. On the other hand, a realistic cut to exclude region (d) would require complete isolation of the leptons. Since lepton detection cannot go below transverse momenta of the order of a GeV, a complete isolation cannot be achieved. Thus, some logarithmic sensitivity to the lepton mass remains, and we must modify the matrix elements by including mass effects in diagrams that involve the splitting of a photon into a leptonantilepton pair. The details of our procedure are given in appendix B.
In order to deal with the remnants corresponding to regions (b), (c) and (d), that have a divergent or a logarithmically enhanced cross section, we implement their production in POWHEG with a reduced sampling rate by applying a suitable suppression factor, that is divided out from the weight of the generated event. With this procedure, if the analysis cuts are applied, we obtain correct, finite results.

NLO results and scale dependence
We begin by studying the scale dependence of the NLO cross section for all combinations of different-flavours or equal-sign lepton scattering. For completeness, we include also the contribution of Z exchange and Zγ interference. We notice that, at NLO, scattering processes of leptons of different flavours and equal or opposite charges have different cross sections even if only the electromagnetic interaction is considered. In fact, the amplitude has contributions of order q l q 2 and q q 2 l , with q l and q denoting the electric charges of the two fermions. These two contributions change relative signs if one lepton is replaced with its antilepton. On the other hand, fully charge-conjugated processes yield the same cross sections and distributions even if we include Z exchange, provided we do not look at parity-sensitive observables. Thus we can omit charge-conjugate processes.
In our analysis we always consider proton-proton collisions at 13 TeV. As discussed in the previous section, at the NLO level we must impose an isolation cut on the final-state leptons, in order to avoid the collinear divergence of the region (d).
We thus require that no other lepton with transverse momentum larger than a given cut is inside a cone in the pseudorapidity-azimuthal plane with given aperture ∆r. We have considered the values ∆r = 0.1, 0.3 and 0.6. We have observed a mild dependence on ∆r, thus we present here only the results for ∆r = 0.3. Altogether, we require both leptons to satisfy the cuts p T > p cut , p cut = 5, 10 and 20 GeV , |η| < 2.4 , isolation: ∆r > 0.3 for leptons p T > 0.9 GeV . factor K F = µ F /µ 0 for different choices of the transverse momentum cut. We have chosen the reference scale µ 0 as the sum of the absolute value of the transverse momentum of all final-state leptons divided by two. This choice is collinear-safe, since for both initial-and final-state collinear singularities it reduces to the transverse momentum of the final-state particles in the underlying Born process. We notice that, in general, there is a considerable reduction in the scale dependence when going from the LO to the NLO result. Furthermore, we observe that generally the intersection of the LO and NLO results is slightly above 1 for opposite-sign, different-flavour leptons, slightly below 1 for same-sign, different-flavour, and near 0.5 for same-sign same-flavour leptons. While the scale behaviour of the results for the electrons and muons seems very reasonable, when the τ is involved, we see an anomalously large scale dependence even at the NLO level, the worse case being τ τ scattering. This behaviour is partly understood as a consequence of the fact that for small scales the τ parton density becomes negative, and in fact we see that for small scales the LO cross section becomes very small. Given the differences in the scale behaviour of the different processes, we have chosen the following procedure for selecting a range of scales appropriate to each process. We have scanned the cross sections for 0.5 < K F < 4, with a spacing of 0.1 for 0.5 < K F < 1, 0.2 for 1 < K F < 2, and 0.4 for K F > 2. We have determined the value of K F where the NLO cross section is maximal, and taken half and twice this K F value as extremes for the scale range. We then quote the maximum and minimum value of the cross section in this range, computed at both the LO and NLO level. The results are shown in table 2, where we also report the predictions obtained with the HT method. The range of the LO results are quite large, while it becomes more modest in the NLO case. The NLO result is always included in the LO band. We observe that the cross sections for like-sign different fermions are substantially different than the unlike-sign ones. We remark that the range of values given here for the LO results differs slightly from the one given in the LUXLep paper, because of the different choice of scales.
We notice that, in comparison, the HT result is almost always larger than the upper limit of the NLO band (apart for e + τ − and µ + τ − ), although not by large amounts. We will see in the following, however, that the NLO+PS results are in much better agreement with the HT ones. We remind the reader that fixed-order NLO results in 2→2 parton scattering are typically not well-behaved also in the case of coloured partons, while POWHEG results perform much better [18]. Furthermore, we remark that, in our case, the NLO+PS implementation recovers in an approximate way effects that are included in the HT result, but not included at fixed order NLO.

Study with realistic cuts
NDY pairs can arise in all SM processes where we produce four leptons. The double Drell-Yan (DDY) process qq → l + l − + − represents a particularly serious background to the lepton initiated NDY pair creation. First of all, the process has a double collinear enhancement when the opposite-sign, same-flavour pairs have a small mass. Especially at the low transverse momenta that are required to get a reasonable yield at the LHC, there is a substantial phase space for two of the leptons to be produced below a detection threshold. If, for example, the l − and − leptons have transverse momenta below a GeV, they will escape detection, and the signal will look like a pair of nearly balanced different-flavour leptons, as depicted in Fig. 17. Thus, when setting up our cuts, we will also consider their    Table 2: Cross sections in femtobarns for lepton scattering with acceptance cuts |η| < 2.4 and p T > p cut for both leptons. An isolation cut is also imposed by requiring no leptons in a cone of radius r < 0.3 with p t > 0.9 GeV. For each subprocess, a range in the factorization scale factor K F is chosen by determining the value of K F that maximizes the NLO cross section and then taking half and twice this value as extremes, as detailed in the text. The quoted cross sections correspond to the minimum and maximum in this range.
generator at the parton level, i.e. at the level of the Les Houches events file (LHE from now on), with that obtained with the HT method. The most remarkable difference shows up in the distributions of the invariant mass of the pair and their azimuthal distance, in Fig. 19. From the figure, it is apparent that, in the HT calculation, there is a large contribution  in the region where the invariant mass of the e − µ − pair is below twice the transverse momentum cut. The origin of this contribution is easily identified as the production of a back-to-back lepton-antilepton pair with relatively large transverse momentum, with the lepton subsequently radiating a photon, that in turn splits into a second lepton-antilepton pair, as depicted in Fig. 20. This mechanism is obviously absent in our NLO+PS calculation, since the γγ initial state for our process arises only at the next-to-next to leading order (NNLO) level. The final state has in this case an l − + − low mass system at large transverse momentum, recoiling against a (high transverse momentum) l + lepton, and the l − − system constitutes our signal. This mechanism does not operate at the level of the POWHEG generated Les Houches events, since there no more than three leptons are produced. It is clear that these events are not the ones we are interested in, especially considering that the DDY process qq → l + l − + − can generate the same kind of final state, through the mechanism qq → γ * → l + (l − * → l − (γ * → − + )). By requiring m l > |p t |, a condition that we dub "cut A" in the following, we strongly suppress this effect, as shown in Fig. 21. We see that after cut A is applied, the two calculations are more consistent. We notice, however, that the agreement of the two calculations after cut A is less than perfect, with the POWHEG result overshooting the HT one for small invariant mass and undershooting it for any azimuthal distance away from π. It is interesting to see that, after shower, the agreement becomes really excellent, as shown in Fig. 22. This result is not too surprising, since the shower Monte Carlo, by adding also the backward evolution for the incoming quark and photon, captures the most relevant features already present in the HT calculation. In order to reduce the large background from the DDY process qq → l + l − + − we also impose the following cuts that we have found to be particularly effective. Thus the full set of cuts that we have adopted is summarized below: i. Both signal leptons must have |η| < 2.4, p t > p cut .
ii. We veto events with extra leptons with p t > 0.9 GeV in a cone of radius r < 0.3.
iii. m l > |p v. If there is another lepton passing cuts i and ii with p t > 3 GeV the event is vetoed (cut D).
The cuts (i) and (ii) are always applied. Their combination will be referred to in the tables as "cut T". The e − µ − and e − µ + cross sections obtained by applying various combinations of cuts are reported in tables 3 and 4. Following our findings on the scale dependence in Section 3.2 we choose as factorization scale the sum of the transverse momentum of all final-state leptons divided by two, ignoring the slight preference for smaller (larger) scale of the e − µ − (e − µ + ) processes. First of all, we notice that the NLO+PS and HT results are in good agreement for all combinations of cuts that we considered. We also see that when the simplest cuts are applied, i.e. cut T, all results are in reasonable agreement. On the other hand, as the complexity of the applied cuts increases, the fixed order results departs from the LH-NLO and NLO+PS ones. This is not surprising, and it is a recurring characteristic of fixed order NLO calculations when the events are cut in order to make them more similar to the Born ones. For example, in our case, when we impose the p t balance with cuts B and C, we reduce the size of the real contributions. These are in fact divergent, and are rendered finite by the inclusion of virtual corrections. Asking eventually perfect balance T  TA  TAB  TABD  TAC  TACD    would veto all real events, yielding an infinite negative result. Thus, depending upon the cuts, one can get results that are as small as one pleases, and that can even become negative at fixed order. This effect is quite visible in our table, where, as the cuts become more aggressive, the NLO results decrease and their scale dependence becomes larger, yielding also negative values. We stress that the LH-NLO result by itself is actually incomplete. For example the transverse momentum of the three hardest leptons is strictly zero at the LH-NLO level. We expect that backward evolution in the shower should lift this constraint, and yield a reasonable transverse momentum distribution for this system, and also reduce further the number of balanced events (i.e. those that pass cut B and C). Therefore we expect that after shower the results should be in better agreement with the HT calculation, that has always the incoming leptons resolved into splitting photons, with the photons in turn being produced in association with a jet (that in the elastic case may consist of a single proton), and in fact this is what we observe in our table. We also show in Fig. 23 the transverse momentum of the signal-leptons pair before and after shower compared with T  TA  TAB  TABD  TAC  TACD table 3 for the e − µ + process.
the HT calculation. We see the excess in the first bin of the distribution when the shower is absent, that is spread out at larger transverse momenta when the shower is included, yielding a better agreement with the HT result.
As a last comment, we see from the tables 3 and 4 that the combination of cuts A + B + D and A + C + D are most effective in suppressing the DDY background. The effect of the B/C and D cuts is better understood from the following figures. In Fig. 24, on the left panel, we show the electron transverse momentum distribution computed at the NLO+PS level in comparison with the DDY result when only cut A is imposed. We see that the lepton spectrum is dominated by the DDY process. On the other hand, we see from the right panel that at low transverse momenta of the lepton pair, the lepton scattering process prevails. The same result when applying the A+B cuts is shown in fig. 25. We see that now the signal in the electron p t emerges from the background at prevails over the background as long as the transverse momentum is not too large. This is due to the fact that the DDY process is initiated by one valence quark, while the lepton PDFs are very soft. Some room is left in tightening the cut on the balance, going from B to C, in order to further improve the background rejection, as can be seen in table 3. Before ending this section, we turn our attention to Fig. 23, where there is a comparison of the transverse momentum distribution of the signal pair obtained at the NLO+PS and HT level. It is quite clear that the modeling of the first few bins is less than satisfactory. We stress that in order to obtain this level of agreement we had to lower the minimum radiation scale of POWHEG from 0.8 to 0.1 GeV, and some parameter tuning on the Herwig7 side was also necessary. This, together with similar observations regarding the /γ-initiated scalar production, reminds us that the current status of the shower generators for /γ-initiated processes is perhaps not quite mature, and, as the interest in these processes grows, may require further work.

Conclusions
In this paper we have used the recently developed photon and lepton PDFs to study the impact of NLO corrections to a number of /γ-initiated processes and their matching to a parton shower. While at hadron colliders there is considerable experience in matching higher-order calculations to parton showers for pure QCD processes, little is known about this matching in the case of /γ-induced processes. For this reason, we have first considered two very simple processes, namely the production of a massive scalar particle which couples to incoming photons or leptons. We have found that LO results display a strong scale dependence, which translates into very large uncertainty bands, of the order of O(30%) for γ-and O(50%) for -induced processes. On the other hand, NLO predictions are stable when varying the factorization scale, with missing higher-order uncertainties of the level of 5% for γ-and 10% for -induced production. It is possible to find a choice for the central scale such that the NLO corrections are small, ranging from a few percent to about 20%. For photon induced processes an optimal choice is the mass of the scalar M Φ , while for -initiated processes a scale close to M Φ /4 is preferred. We notice in particular that, in the case of -initiated processes, the scale choice M Φ gives rise to large and negative corrections, with the central LO and NLO results differing by a factor 2 and with no overlap between the two uncertainty bands. Thus limits for new physics based upon searches using LO predictions at a scale M Φ would be overly aggressive, and it would be better to either use a lower scale, or rely on a full NLO calculation.
For photon induced processes we matched the LO and NLO calculations to both Pythia8 and Herwig7 using POWHEG. We find considerable differences between the two predictions at LO, while we obtain a much better agreement after matching to NLO. On the other hand, the two shower Monte Carlo predictions for the very low transverse momentum distribution of the scaler have a very different shape. These features are not relevant for searches of heavy resonances. However, they give an indication that a more accurate tuning of the shower generators is desirable for an accurate modeling of γ-initiated processes. We note that, while in general Monte Carlos are tuned to data, in the present case the tuning could be performed using predictions obtained with the HT method, since they receive no strong corrections.
The state of the art for -initiated processes in shower Monte Carlo programs is even less advanced than for γ-initiated ones, since they are currently implemented only in a preliminary version of Herwig7. When comparing LO and NLO matched predictions for the transverse momenta of the scalar and of the accompanying jets/leptons, obtained with the optimal scale choices (M Φ for γ-and M Φ /4 for -initiated processes), we observe a rather flat ratio as long as the transverse momenta are not too small. However, we find considerable differences at small transverse momenta, where NLO corrections have an opposite effect in γ-and -initiated processes. Concerning the -initiated processes we find, as expected, that Herwig7 generates small associated hadronic activity, with the transverse momentum spectrum of the leading jet being comparable to the one of the sub-leading lepton. This information is important to enhance the signal to background ratio when looking for these processes. Once available also in Pythia8, it would be certainly interesting to compare Herwig7 and Pythia8 matched predictions also for -initiated processes.
We have then considered Standard Model 2 → 2 lepton-scattering processes, focusing our attention on final states which are not Drell-Yan like, and extending to NLO accuracy the LO results presented in the LUXLep paper. While the NLO calculation is per-se straightforward, care must taken in its implementation in POWHEG, since, starting from NLO, the calculation receives contributions from other electromagnetic sub-processes which involve additional initial-and final-state collinear singularities. In particular, because of the latter, at NLO one needs to impose additional isolation cuts for the final-state signal leptons.
For the lepton-scattering processes we have also performed a calculation in the HT approach, which uses the full hadronic tensor and therefore receives no higher-order QCD corrections. Since the HT calculation misses instead collinear emissions of photons from initial-state leptons, we have estimated these effects by deriving a PDF set where such collinear emissions are switched off. We find that their impact is below one percent for the processes at hand. This implies that one can use the HT calculation as a reference prediction. We note however, that the HT calculation is much more involved than an NLO calculation and that its matching to a parton shower is not straightforward.
As in the case of the production of a scalar particle, we find a considerable reduction in the scale uncertainty when going from LO to NLO. In general, setting the factorization scale to the transverse momentum of the outgoing leptons is a reasonable choice, with slight preferences to lower or higher scales according to the process.
As for the comparison of the HT predictions and the NLO and NLO+PS calculations, we have found the following features. When minimal isolation and acceptance cuts are applied, we find discrepancies between the NLO and the HT result due to a production mechanism that is included in the HT calculation but not in the NLO one, i.e. the production of a back-to-back lepton-antilepton pair with relatively large transverse momentum, followed by the collinear radiation of a photon that in turn splits into a second collinear lepton-antilepton pair. This mechanism enters only at NNLO in the collinear factorization approach. It is clear that these events are not the ones of interest, in particular given that the dominant double Drell-Yan (DDY) background process can generate the same kind of final state. Once one imposes fiducial cuts to suppress this region one finds good agreement between the HT and NLO results.
When considering additional cuts, in order to suppress the large DDY background, the NLO results start departing from the HT one. This is a typical feature of fixed-order NLO calculations in collider physics. As one tightens the cuts in order to make the process look more Born-like, large logarithmically-enhanced negative virtual corrections are exposed that spoil the reliability of fixed-order results. However, once the NLO is matched to a parton shower, after a minimal tuning of the latter, we find very good agreement with the HT result even after stringent cuts. This is also due to the fact that the parton shower generates contributions that are present in the HT calculation and missing in the NLO one. Furthermore, there is certainly room for improvements as the implementation of -initiated processes in parton-shower programs is only at its early stage.
As of now, physics with initial photons and leptons have been mostly considered for BSM searches. It would be certainly interesting to experimentally study the leptonscattering processes, since these are the only SM processes that can probe the lepton content of the proton.
A The computation of lepton scattering in the hadronic tensor approach Following the notation of the LUXLep paper, we write the cross section as where P i (i = 1, 2) are the momenta of the incoming protons, S = (P 1 + P 2 ) 2 , q i are the four momenta of the virtual photons in the DIS convention (i.e. directed towards the protons), k l denote collectively the momenta of the outgoing leptons, dΓ is the corresponding phase space, and T µµ ;νν is the squared matrix element for the production of four leptons initiated by the fusion of two off-shell photons. The process is depicted in Fig. 27, where we show two sample Feynman diagrams, involving t-and s-channel production mechanisms.
We computed the phase space by composing nested 2 → 2 scattering processes, loosely following the approach of Ref. [20]. We thus rewrite eq. (A.1) in terms of the three-momenta p i and masses √ s i of the outgoing hadronic systems, as where p 0 i = p 2 i + s i , and the s i must obey the following constraints where the m l are the masses of the outgoing leptons. The hadronic tensor W µν is expressed in terms of the structure functions F 1 and F 2 (comprising here also the contribution of the elastic form factors), that in turn are functions of Q 2 i = −q 2 i and x bj,i . The latter is obtained from the relation The elastic component of F 1 and F 2 carry a factor of We must compute all four combinations of elastic and inelastic contributions. Thus, when considering the elastic case, we have to replace take the p i momentum as having mass m p , and write the elastic F i contributions without the δ(1 − x bj,i ). We computed T µµ ;νν using Madgraph5 aMC@NLO [21] in stand-alone mode. We have generated the process γγ → l + l − + − for different and identical fermions. Madgraph, by default generates the diagram for given photon polarizations. It does not assume, however, that the incoming momenta are on shell. Thus, we replaced the incoming helicities in the generated formulae with simple Lorentz indices. More specifically, we implemented a routine for the vector-fermion-antifermion vertex that, rather than requiring a helicity in input, requires a Lorentz index. This routine was obtained by modifying the VXXXXX(P,VMASS,NHEL,NSV,VC) Helas [22] routine in such a way that, when called with NHEL = 0, 1, 2, 3, it sets all the VC(3 : 6) components to zero except for the NHEL+3 one, that is set to 1.
The tensor T should vanish when contracted with the momenta corresponding with the two incoming photon, and this served as check for our procedure.
The hadronic tensor is written as and where we have defined Q 2 = −q 2 . The fitted/calculated values of F 2 and F L are available in appropriate routines that were developed in the LUX papers. Their value is in part extracted from electro-production data, and in part is calculated in a way that is described there in great detail. The phase space in eq. (A.1) requires great care, since it is logarithmically singular in several regions, the logarithmic singularities being regulated by the proton and the lepton masses. Following Ref. [20], we factorize the phase space as a sequence of nested two-body components. Our main integration variables are thus the s-channel and t-channel invariants of these two-body components, that are chosen in such a way that the logarithmically enhanced regions are adequately sampled.
is decomposed according to the following sequence: P 1 + P 2 → p 1 + (k 1 + k 2 + k 3 + k 4 + p 2 ), (P 1 − p 1 ) + P 2 → (k 1 + k 2 + k 3 + k 4 ) + p 2 , (P 1 − p 1 ) + (P 2 − p 2 ) → k 1 + (k 2 + k 3 + k 4 ), (P 1 − p 1 − k 1 ) + (P 2 − p 2 ) → (k 2 + k 3 ) + k 4 , where the momenta in bracket denote a massive system, with an appropriate bound for its mass. For example, for the first line of the equation, the minimal mass of p 1 is m p , and its maximum is √ S −(m 1 +m 2 +m 3 +m 4 +m p ); the minimal mass of the (k 1 +k 2 +k 3 +k 4 +p 2 ) system if m 1 + m 2 + m 3 + m 4 + m p , and its maximum is √ S − m p . The structure functions have an elastic and inelastic component. In the elastic case, p i has fixed invariant mass equal to m p . Thus, all components, elastic-elastic, inelastic-elastic, elastic-inelastic and inelastic-inelastic must be generated independently. As a further point, the procedure outlined above is adequate for the sampling of the t-channel singularities if lepton 1 is produced mainly along the P 1 direction, and lepton 4 along the P 2 one, a region that we here dub region a. In order to perform an adequate sampling, we define an enhancement factor that singles out region a f a = 1 (P 1 − p 1 − k 1 ) 2 (P 2 − p 2 − k 4 ) 2 , (A. 15) and one that singles out the complementary region b where 1 is produced along the P 2 direction and 4 along P 1 f b = 1 (P 1 − p 1 − k 4 ) 2 (P 2 − p 2 − k 1 ) 2 .
(A. 16) We then multiply the cross section by the factor f a f a + f b , (A.17) thus suppressing region b. The contribution of region b is obtained by generating the phase space with the role of the two flavours exchanged, providing again the similar suppression factor. Several other regions of singularity are present, that however do not contribute significantly to our signal. In order to speed up convergence for the distributions that we are interested in, we suppress these regions with adequate suppression factors. The generated events have this suppression factor removed, so that one obtains the exact cross section, but with less dense sampling in the uninteresting regions. These regions are • The region when the k 1 , k 2 and/or the k 3 , k 4 systems have small masses. These regions do not generate isolated leptons with the appropriate charge. We thus suppress these regions with a factor proportional to these masses.
• The region where the signal leptons k 3 and k 4 are not in a reasonable acceptance range. In particular, they must have a sizeable transverse momentum, and their pseudorapidity should be within acceptance. The suppression factors that we adopt are In this appendix we detail how we include mass effects to cure the collinear divergence arising from squaring the photon splitting diagrams depicted in Fig. 28. First, let us remark that this problem does not show up in a complete massive computation. However, such a computation cannot be performed in the collinear factorization approach. Second, these contributions have as underlying Born the Compton process; thus, in our treatment, they contribute only to the remnant cross section, as discussed in Sec. 3.1. This means that we have to deal only with the real matrix elements, with no need to modify subtraction counterterms.
We proceed as follows. We start by computing the square of the two diagrams in Fig. 28 assuming the green fermion line, associated to lepton l, as massless and the red one, associated to lepton , as massive. Observe that in the limit for a vanishing mass the resulting function, that we dub f c (p; m ), smoothly approaches the massless squared matrix element that exhibits a collinear final-state singularity. Next, we look for all possible ¯ pairs in the final state; there are two in the case of identical leptons and just one otherwise. For each pair, we let the momenta of these two massless leptons acquire mass by retaining their 3-momentum while correcting their energy. The resulting real event has a slightly increased centre of mass energy. In formulae, writing the original set of momenta p i in the partonic centre of mass defined by where √ s is the partonic centre of mass energy, the set of massive momenta is given by In the case of identical leptons, we add the correction factor ∆ c for each of the two pairs found, and divide by the appropriate symmetry factor. For the sake of completeness, we report here the explicit expression of the correction function f c assuming that the collinear leptonic pair has momenta p 3 and p 4 and that p 2 is the incoming photon: