Leptons in the Proton

As is the case for all light coloured Standard Model particles, also photons and charged leptons appear as constituents in ultrarelativistic hadron beams, and admit a parton density function (PDF). It has been shown recently that the photon PDF can be given in terms of the structure functions and form factors for electron-proton scattering. The same holds for lepton PDFs. In the present work we set up a calculation of the lepton PDFs at next-to-leading order, using the same data input needed in the photon case. A precise knowledge of the lepton densities allows us to study lepton-initiated processes even at a hadron collider, with all possible combinations of same-charge, opposite-charge, same-flavour, different-flavour leptons and leptons-quarks, most of which cannot be realized in any other foreseeable experiment. The lepton densities in the proton are extremely small, so that their contribution to Standard Model processes is generally shadowed by processes initiated by coloured partons. We will show, however, that there are cases where these processes can be relevant, giving rise to rare Standard Model signatures and to new production channels, that can enlarge the discovery potential of New Physics at the LHC and future high energy colliders with hadrons in the initial state.


Introduction
The current LHC research program is on one hand aiming at high precision measurements, to spot any deviation from the Standard Model, and on the other hand at the direct search of particles arising in New Physics scenarios. The vast majority of New Physics searches carried out at the LHC regards processes initiated by coloured partons, and lepton initiated processes are relegated to future colliders involving leptons. On the other hand, the current LHC and its planned high-luminosity (and eventually high-energy) upgrade is the collider that will provide the largest part of new high-energy particle physics data in the next 20 years, and its reach in the lepton-initiated channels should also be exploited. In fact, it is well known that quantum fluctuations can give rise to the presence of leptons inside a proton, although with a much smaller relative abundance with respect to the coloured partons. When leptons initiate a large momentum transfer scattering, the process becomes perturbative, the lepton densities also obey an evolution equation, and a partonic calculation of the process, including higher-order corrections, becomes possible. It is thus natural to explore what is the reach of the LHC as far as lepton-initiated processes are concerned, also considering the fact that at the LHC the three charged leptons contribute democratically, and thus processes (such as, for example, µτ scattering) that are not available at lepton colliders may be accessed.
It is the aim of the present work to derive a precise determination of the lepton densities inside a proton, which can then be used to compute processes involving leptons in the initial state at hadron colliders.
In Refs. [1,2] (we will refer to these references as LUX1 and LUX2 respectively, and as "LUX papers" for both of them) it was first pointed out that the parton distribution function (PDF) of the photon in the proton can be computed with high precision using only information from electron-proton scattering data. In these works it was also pointed out that, for similar reasons, this is also the case for lepton PDFs. In the present paper we undertake the task of computing the lepton PDFs in a framework that is very similar to the one adopted in the LUX papers, and in fact by also heavily using the numerical code developed there.
In order to compute the lepton PDF, we consider here a fictitious Deep-Inelastic-Scattering process involving a lepton in the proton, i.e. the collision of a (fictitious) massless scalar with a proton. The massless scalar interacts with the lepton via the vertex ψ h ψφ + c.c., (1.1) that turns the light lepton ψ into a fictitious heavy lepton ψ h , carrying a mass M much larger than typical hadronic scales. According to the parton model approach, such process can be computed in terms of the light-lepton parton density f (x), using the standard QCD factorization formula. It can also be computed in terms of the electromagnetic current structure functions. In Fig. 1 we show schematically the representations of both computations. By relating the two results, we can obtain the lepton parton density, entering the second computation, in terms of the leptoproduction structure functions, entering the first one. The factorization theorem guarantees that the PDF so obtained is independent of the particular process used in the calculation, as was shown explicitly in the photon case in LUX2. In Ref. [3] the LUX method was applied to the computation of the PDFs of the W and Z bosons. We stress however, that in the case of the W and Z, thanks to their large mass, the electroproduction structure functions are needed only in the perturbative regime, and thus the whole calculation can be carried out in perturbation theory, in a way that closely resembles the computation of the heavy flavour parton density [4]. The case of light leptons (a) Structure function computation (b) Parton level computation Figure 1: Our basic fictitious process, with a scalar of momentum r scattering off a light lepton and turning it into a heavy lepton of mass M , represented by the thick red fermion line. In (a) we show the sum of the two diagrams that relate this process to the ep scattering structure functions. In (b) we show the diagrams that enter the calculation of the same process at next-to-leading order according to the QCD factorization formula. Notice that the second and third diagram in (b) (unlike the (a) diagrams) are computed for an on-shell photon, and the collinear singularity from the photon splitting into leptons that arises there is subtracted.
is instead more similar to the photon one, where the structure functions are needed also in the very low Q 2 region, and thus must be extracted from low Q 2 experimental data.
The paper is organized as follows. In Sec. 2 we present our calculation of the lepton PDFs. We first define our target accuracy, that is based on a careful counting scheme of the strong and electromagnetic coupling constants and of the relevant logarithms. We then proceed to the calculation of the lepton PDF in the limit of zero lepton mass, first in terms of the electroproduction structure function, and then according to the parton model formula at next-to-leading order (NLO). We use the two results to extract a formula for the lepton PDF. We then illustrate how the result changes when the lepton mass is included in the calculation. Finally, we explicitly verify that our lepton PDF satisfies the Altarelli-Parisi evolution equation [5], including QED splitting processes that also involve a term of second order in QED.
In Sec. 3 we explain our procedure to assess the theoretical uncertainty of our final result, which closely follows the one used in LUX2. In Sec. 4 we describe how one can add our lepton PDFs to any full LHAPDF set and we do this in the case of the NNPDF31 nlo as 0118 luxqed set of Ref. [6]. In Sec. 5 we show a number of results that validate our procedure. In Sec. 6 we present a number of phenomenological applications of our lepton PDFs. In particular we consider rare SM signatures of different flavour isolated di-lepton production; the production of leptophilic Z ; the production of doubly charged Higgs; and the production of leptoquarks. For this last case, we show that we can reach unexplored regions of the parameter space using already existing data from the LHC. Finally, we give our conclusions in Sec. 7. In the Appendices A-D we provide further technical details.

Details of the calculation
We now illustrate our calculation. We first compute our probe process in terms of the electroproduction structure functions. Then we compute the same process in the parton model, and combine the two results to extract the lepton PDF. We finally verify that our PDF satisfy the Altarelli-Parisi equations to the appropriate order.
Before we begin, it is useful to clarify what accuracy we expect from our calculation. To this end, we consider the parton densities for quarks and gluons as being of order one. In fact, perturbatively they may be considered as sum of terms of order (α s L) n , where L = log(µ 2 F /Λ 2 ), µ F is the factorization scale, and Λ is a typical hadronic scale. The strong coupling α s is evaluated at a scale of order µ F . Thus, since α s ≈ 1/L, all these terms are of order 1. NLO corrections to the quark and gluon PDFs have the form α s (α s L) n .
Besides these rules, we should also worry about the possible impact of higher-order electromagnetic effects. These may be relevant if some logarithmic enhancement compensates for the smallness of the electromagnetic coupling constant. We will adopt the same criterium used in the LUX papers, i.e. that α is of the same order as α 2 s . The photon density is of order αL(α s L) n , so, in our counting scheme, it is equivalent to αL. NLO corrections to it are of order α and according to our α ≈ α 2 s rule, we should also include terms of order α 2 L 2 . The lepton PDFs are of order α 2 L 2 , and their NLO corrections of order α 2 L and α 3 L 3 . However, we will see in the following that the terms of order α 3 L 3 are in practice much smaller than the α 2 L terms, and thus the α ≈ α 2 s rule is a quite conservative assumption.
In the present work we aim at NLO accuracy. The computation of the probe process in terms of structure functions involving the graphs of Fig. 1 (a) includes all terms of order α 2 . All possible strong corrections are already included in the electroproduction structure functions. However, in our calculation we also need to include NLO corrections of relative order αL. Corrections of this kind are already present in the electroproduction structure functions (they arise from collinear photon radiation from quarks) and in the self-energy corrections to the photon propagator. We can account for the latter by using the same effective QED coupling used in the LUX papers. The only term of relative order αL that we miss arises from collinear photon radiation from the light lepton. These contributions, however, are easily included using the evolution equations, with a method that will be described in due time. In the parton model calculation (the diagrams in Fig. 1 (b)) the counting of the order goes as follows. The first diagram is of order α 2 L 2 (i.e. the leading order of the lepton PDF). The two remaining diagrams are of order αL (i.e. the leading order of the photon PDF) times α, leading to an NLO contribution of order α 2 L. As one can easily convince oneself, no other NLO corrections arise here, since in the QCD-improved parton model calculations no large logarithms can arise from radiative corrections.
We clarify from the start that throughout this paper we refer to the lepton density as the density of either charge, i.e. not the sum of lepton and antilepton, and our LHAPDF implementation returns the density of each signed lepton. In our approximation the lepton and antilepton densities are equal, and remain equal at higher QCD orders. Differences arise only as subleading electroweak effects that are not considered here.

Calculation in terms of Structure Functions
We begin by considering the scattering process where φ(r) denotes a scalar of momentum r, γ(−q) a photon of momentum −q (Q 2 = −q 2 ), ψ h (k, M ) is the (hypothetical) heavy anti-lepton of mass M and momentum k, and ψ(k, 0) is a massless lepton of momentum k. We define the kinematics of the process in terms of the following variables We introduce the following dimensionless variables: In the parton model language x can be identified with the fraction of momentum of the lepton with respect to the proton; x with the fraction of momentum of the photon with respect to the proton; z with the fraction of momentum of the lepton with respect to the photon that has created it; x bj with the fraction of momentum of the quark with respect to the proton; and z with the fraction of momentum of the photon with respect to the quark that has emitted it. Summing the two diagrams in Fig. 1 (left) we obtain the amplitude for this process from which, upon integration over the two-body phase space, we obtain the leptonic tensor where we have implicitly assumed the sum and averages over the spin of the external particles. The cross section can then be written as 1 where W µν is the standard hadronic tensor, which, for the scattering of a photon of momentum q off a proton of momentum p, has the form (2.8) For ease of notation, here and in the following we will omit the arguments of the structure functions, which will be always evaluated at (x bj , Q 2 ). We also introduce the longitudinal structure function which is of order O(α s ) relative to F 2 . We will write our results using F 2 and F L instead of F 2 and F 1 .
In order to make contact with the results of the LUX papers, using the identity we rewrite Eq. (2.7) as where now the inner integral matches exactly Eq. (1) of Ref. [1], provided M 2 is replaced everywhere by the invariant mass of the heavy-light leptons system, E 2 cm . Hence we can use the same phase space as in Eq. (3.6) of LUX2 and obtain Note that the variable x defined in the LUX papers as x = M 2 /S should be replaced here by x = E 2 cm /S. For the leptonic tensor, writing explicitly the dΦ 2 phase space, we obtain where θ is the angle between k and r in the centre-of-mass (CM) of the scalar-photon system. The leptonic tensor is gauge invariant, and hence it can be written as where L 1 , L 2 are functions of z , Q 2 and M 2 . At this point, we have all the elements to compute L µν W µν in terms of the proton structure functions. The expression obtained is rather lengthy, hence we do not report it here. We note however that the result simplifies considerably if we only retain the terms that are relevant to our approximation. It turns out that the expression for L µν W µν has schematically the form (2.15) where P and R are rational functions of their arguments. We have indicated schematically with F the linear dependence of the result upon the structure functions. Furthermore we have defined . (2.16) The log(M 2 /Q 2 ) arises in the leptonic tensor from the integral in d cos θ. In fact, it is easy to see that the first diagram in Fig. 1 (a), in the limit of small Q 2 has a collinear divergence when the anti-lepton is produced in the forward direction. The P and R coefficients can be separated in the following terms: 1. Terms that behave as m 2 p /Q 2 for small Q 2 .
2. Terms that do not depend upon Q 2 .
3. Terms that vanish at small Q 2 .
The terms of the first item give rise to an integral of the form multiplying the structure functions, and, in the case of P , by L. Neglecting the mild Q 2 dependence in the structure functions and in L, this integral is dominated by small values of Q 2 ≈ m 2 p (since the lower limit of integration in Eq. (2.12) is proportional to m 2 p ). Thus, the contribution proportional to P is of order log M 2 /m 2 p (arising from the L coefficient), while the one proportional to R is of order 1. We only need to keep the former, that in our approximation is NLO.
The terms of the second item give clearly origin to single and double logarithmic enhancement in the contributions proportional to P , and to single logarithm enhancement in the contributions proportional to R.
The terms of the third item lead to integrals (dominated by large values of Q) that are of order one, and thus negligible in our approximation. We also notice that the third logarithm in Eq. (2.16) vanishes for small Q 2 , and thus can be neglected for the same reason. With these simplifications the cross section in Eq. (2.12) becomes where we have introduced the q → qγ and γ →ll splitting functions We stress again that in the above expression the structure functions are evaluated at (x bj , Q 2 ). In Eq. (2.18) we have also dropped terms proportional to F L that are not multiplied by a large logarithm, since F L is down by one power of α s with respect to F 2 , and thus leads to a contribution that is subleading in our counting scheme.

The parton model calculation
We now present the result for the computation of the same cross section using a parton model calculation. The equivalence between the two expressions will allow us to derive a formula for the lepton PDF in terms of the hadronic structure functions. The details of the partonic calculation are reported in Appendix A. The final result is where z is now given as a function of x, z = M 2 /E 2 cm = M 2 /(Sx), and σ B = π.

Extraction of the lepton PDF
In order to extract the lepton PDF we identify the two expressions for σ in Eq. (2.18) and Eq. (2.20).
We obtain where we have replaced dz = z dx/x, and used x = M 2 /S. We recall that F 2 , F L are evaluated at (x bj , Q 2 ), with x bj = x/z. We now recall the expression for the photon PDF, Eq.(6) of Ref. [1]: where again the structure functions are evaluated at (x/z, Q 2 ).
We now observe that if we replace the upper limit in the Q 2 integration in Eq. (2.21) with µ 2 F /(1 − z) we obtain an equivalent expression up to subleading terms, since the difference only involves values of Q 2 near the upper limit, and thus no large logarithms. Proceeding in this way, and substituting Eq. (2.22) for f γ , we obtain where we have used x = M 2 /S and z = x /x. The structure functions F 2 and F L are always evaluated at x bj = x/z and Q 2 .
The above formula has been derived in the massless limit for the physical lepton. As shown in Appendix B, the effect of lepton masses up to our required accuracy are simply accounted for by replacing . (2.24) Including lepton-mass effects in Eq. (2.23) and neglecting higher-order terms proportional to F L , we thus obtain This is our final expression for the lepton PDF. It can be evaluated numerically, similarly to what was done recently for the photon PDF. We observe that, compared to the latter, the lepton PDF requires one extra integration. The structure functions F 2 and F L are the only functions in formula (2.25) that are not known analytically. They depend only upon the two variables x bj and Q 2 . It is thus possible to express formula (2.25) as an integral in x bj and Q 2 , and a third variable, from which the integrand depends analytically. In Appendix C we provide some details regarding this simplification of the integrand. The integration in the third variable is thus simpler to perform, either with numerical methods (e.g. using Gaussian integration) or analytically.

Verifying the Altarelli-Parisi evolution
The lepton PDF given in formula (2.25) must satisfy the Altarelli-Parisi equation [5]. Given that Eq. (2.25) includes accurately terms of order α 2 L 2 and α 2 L, where L = log(µ F /Λ), its logarithmic derivative must contain accurately terms of order α 2 L and α 2 . Taking the derivative of formula (2.25) we obtain where in the first term we have taken a derivative with respect to the explicit µ 2 F dependence of the logarithm, while in the second term we have taken a derivative with respect to the upper limit of the integration. In doing so we have neglected terms of order m 2 /µ 2 F that arise in the argument of the logarithm, and we have replaced and neglected the higher-order terms. In the first term of Eq. (2.26), the z integral corresponds to the LUX expression of the photon parton density (Eq. (2.22)), except that it does not include the term outside the Q 2 integral (this was referred to as the MS correction in the LUX papers). We can thus replace this expression with the photon parton density, adding a term to compensate for the lack of the MS correction. We get where in the second equality we have transferred the subtracted MS correction to the second term. We now rewrite the integral of the second term as where now x = zx bj , and replace F 2 with The expression in the curly bracket is equal to the function p s (ξ) in Eq.(60) of Ref. [7], that enters the P lq splitting function. We stress that, in our case, this term is of order 1/L relative to the first term. In fact, f γ is of order αL, so that the first term is of order α 2 L, while the second term is of order α 2 , since f q is of order 1 in our counting scheme.
Of course there are other terms of order α 2 in the second order QED evolution, but they multiply either f γ or f , and are thus subleading in our counting scheme. We also remind the reader that our expression for the lepton PDF does not include terms of order α 3 L 3 , that, if present, would give rise to the term on the right-hand side of Eq. (2.30). This term is of order α 3 L 2 in our counting scheme (f is of order α 2 L 2 ). Since we assume L −2 to be of order α, this amounts to a NLO correction to the evolution that should be present.

Theoretical error due to missing higher-order effects
In order to estimate the theoretical error due to missing higher-order effects, we parallel the method proposed in LUX2. There, the upper limit in the Q 2 integration yielding the LUX photon parton density was modified using a generic z-dependent form M 2 (z) (see sec. 9.1 in LUX2). This modification was compensated by a corresponding modification of the MS conversion term. As we will discuss in the following section, our determination of the lepton densities is performed together with a determination of the photon density. We thus apply the same method to our lepton-density formula Eq. (2.25), replacing the upper integration limit µ 2 F /(1 − z) with M 2 (z). In the lepton case, this variation is already of one order above our target accuracy, and does not give rise to any modification of the MS conversion term. Thus, following LUX2, we consider two forms for M 2 (z) and take µ M to be a multiple of µ F , to be varied by a factor of two above and below µ F . The corresponding range of results is our estimate for the theoretical error due to missing higher-order effects.

Construction of a PDF set with leptons
We now illustrate our construction of a PDF set including photons and leptons. This set is based upon the NNPDF31 nlo as 0118 luxqed set of Ref. [6], and will be made available as an LHAPDF set under the name LUXlep-NNPDF31 nlo as 0118 luxqed. Here we refer to it simply as the LUXlep set. For brevity, in the following, we will also refer to the NNPDF set upon which it is based as the NNPDF set. The construction of a PDF set with leptons relies upon electroproduction data, both for the elastic and the inelastic case. The same data and fits used in the LUX papers (see Refs. [8][9][10][11][12][13][14][15][16]) are used here.
We constructed a PDF set with leptons starting from the NNPDF set. For better consistency, we generate the photon ourselves using the LUX approach. We implemented our formula for the leptons, Eq. (2.25), by suitably extending the computer code developed for the LUX papers.
We proceeded as follows.
1. We compute the lepton and photon densities at a reference scale µ F (our central value will be µ F = 20 GeV), using formula 2.25, as a function of the electromagnetic structure functions and form factors. The structure functions in the perturbative regime are evaluated using a member m of the NNPDF set.
2. We take the member m at a reference scale µ pdf and evolve it to µ F using Hoppet [17]. The photon evolution is included at order α and αα s , but leptons are not included in any splitting function. This evolution step matches what is done in the NNPDF set, where leptons where simply not included. We do this step ourselves for better stability of the results.
We did not take µ pdf equal to the initial NNPDF evolution scale µ (0) F = 1.65 GeV in order to avoid an excessive sensitivity to the evolution implementation in NNPDF. In fact, we want to use Hoppet for the evolution, and eventual subleading differences with the evolution implemented in NNPDF would manifest themselves especially at low scales, where lower and higher-orders effects become closer in size.
We can also choose µ pdf = µ F , in which case the evolution step mentioned above is avoided. This is in fact our default choice.
3. In the NNPDF set evaluated at the scale µ F , we replace the photon density with the one we computed, and add our computed leptons densities.
4. Using Hoppet, we evolve the set so obtained down to the initial scale µ (0) F . This step of evolution includes the QED splitting functions at order α and αα s , excluding those involving a lepton radiating a photon. We do this because in our calculation of the lepton PDFs, photon radiation from leptons (that is of order α 3 L 3 in our counting scheme) is not included, while, in our approximation, it should be (see the discussion in the introduction of section 2). We implement this radiation by evolving the lepton density from the scale µ (0) F using the full QED evolution. In order to do this we need the lepton densities at the low scale, and we obtain them with the procedure we just outlined. Notice that we cannot compute the lepton densities directly at a low scale, because, at low scales, our calculation is not guaranteed to satisfy the Altarelli-Parisi equation due to the power suppressed effects it includes. We thus perform the computation at a scale that is large enough for power suppressed effects to be negligible, and use the QED evolution (without lepton radiation, since this is not included in our calculation) to produce a partonic lepton density (free from power suppressed effects) at the low scale.
In addition, we also added to Hoppet the splitting function P lq , which is of order α 2 , and included its effects at this step. Notice that the inclusion of the P lq splitting function is mandatory to preserve our accuracy. In fact, as already stated in sec. 2.4, P lq contributes at the NLO level to the evolution of the lepton densities, since it multiplies directly a quark density, that it is of order one in our counting scheme. The leading term in the evolution is of order α times a photon density (that is of order α/α s ) and is thus of order α 2 /α s . Thus the P lq contribution, of order α 2 , is down by a single power of α s with respect to the dominant term.
5. Starting with the set at the low scale obtained in this way, we generate the PDFs at any scales using Hoppet with full α and αα s plus P lq evolution including leptons.
The procedure of item 4 can be avoided by evaluating directly the α 3 L 3 contribution arising from lepton electromagnetic radiation. This calculation is illustrated (and compared with the method of item 4) in Appendix D. In the appendix it is also shown that this effect is small, relative to the NLO effects (of order α 2 L) that are included in our formula (2.25). Thus, the rule α ≈ α 2 s that we have adopted in our counting seems to be quite conservative, and omitting the α 3 L 3 would only lead to a minor error.
The procedure illustrated in items 1 to 5 is applied to all members of the PDF set. Furthermore we also apply it to the central (m = 0) set modified with the addition of the uncertainty variations described in LUX2 in section 10.2, labelled as (EFIT), (EUN), (RES), (R), (M), (PDF), (T) and (HO). We briefly summarize their meaning.
(EFIT) The uncertainty on the elastic contribution induced by the fit of the form factors [8,9], as was done in LUX2.
(EUN) The uncertainty that comes from replacing the fit to the elastic form factors [8] including polarisation data with the fit with only unpolarised data, as in LUX2.
(R) A modification of R (the ratio of the longitudinal to the transverse electroproduction cross section) by ±50% around its central value, as in LUX2.
(M) A modification of the Q 2 PDF scale which governs the transition from fitted data for F 2 and F L to a PDF-based evaluation, as in LUX2.
(PDF) The input PDF uncertainties for Q 2 > Q 2 PDF according to the default prescription for NNPDF, as described in the following.
(T) A twist-4 modification of F L , as in LUX2.
(HO) An estimate of missing higher-order effects, as described in the following.
For the higher-order uncertainties (HO) we cannot use the method described in LUX2, since we did not perform an NNLO calculation of the lepton densities. We thus consider the scale variations described in sec. 3 with the following five choices: , and then take the largest deviation from the central set (in absolute value) as further uncertainty variation. We now must implement our variations in a way that is consistent with the meaning of the NNPDF members, i.e. as replicas. We have used the following approach. We assume that each of our variation, suitably symmetrized, corresponds to a Gaussian error, with variance equal to the variation itself. Thus, for the generic member m of the NNPDF set (excluding m = 0, i.e. the central one), for flavour i and given x and µ F values, we compute the correction where f where N rep is the total number of replicas (100 in the set we are considering). Eq. (4.2) guarantees that the average of all replicas remains the same, i.e. equal to the central member.
We have chosen to compute the lepton and photon PDFs at the initial scale Q = 20 GeV. We do not add an error associated with this choice, since it is much smaller than our estimate of the error due to higher-order perturbative effects. In particular, choosing the much higher value Q = 100 GeV, we get a variation of the electron density at 100 GeV below 0.5% across the whole x range.

Validation
We illustrate now how our sources of uncertainties affect the lepton densities. First of all, however, we want to show that the photon PDF, that we compute here using the NLO LUX approach in conjunction with the NNPDF set, has uncertainties that are consistent with what we found in the original LUX papers. The errors are reported in Fig. 2 at the reference scale µ = 100 GeV. The PDF uncertainty is obtained with the usual method x LUXlep, µ = 100 GeV required by PDF sets with N rep replicas. One defines by summing in quadrature the deviation of the photon density of each replica with respect to the central set, dividing the sum by the number of replicas, and then taking the square root. We see that this figure compares well with Fig. 15 of LUX2, the only marked difference being the size of the HO uncertainty, that is much smaller there. On the other hand, this difference is justified if we look at Fig. 13 of the same reference, where both the uncertainty bands of the NLO and NNLO computation of the photon density are illustrated. In Fig. 3 we show the uncertainties for the electron, the muon and the tau at the same reference scale. We see that the uncertainties are in line with those of the photon PDF illustrated in Fig. 2, except for the error due to HO effects, that seems larger for leptons, and slightly larger for larger lepton masses.
We now turn to the validation of the LUXlep set. In order to test the consistency of our evolution machinery, we checked that the QCD partons are consistent within errors in the LUXlep and NNPDF sets. In Fig. 4 we show the comparison for the down quark and the gluon density. The central values and the uncertainty bands perfectly overlap in the central x-region. Only at very small x, one sees a 1% deviation between the central values, well within the uncertainty bands. As expected, the introduction of the lepton densities does not alter significantly the quark and gluon distributions, and gives a negligible contribution to the total proton momentum (we find a momentum sum equal to 1.00065 at µ F = 100 GeV). We thus, at variance with the LUX papers, did not apply any correction to restore the momentum sum rule.
In Fig. 5 we compare the photon densities in the LUXlep set and in NNPDF sets. We find a reasonable agreement, with the two uncertainty bands nicely overlapping within few percents. The LUXlep band is larger since for consistency we recompute the photon PDF with the LUX approach at NLO while in the NNPDF set it was obtained at NNLO. This difference explains also the small deviation of the central values in the central-and high-x region, with NNPDF being bigger.
We conclude this section discussing the impact of the P lq splitting in the evolution of the lepton densities. We compare, for different scale choices, the lepton densities directly computed with the LUX approach, with the ones obtained by using the LUX approach at a fixed reference scale (chosen to be as before µ 2 = 400 GeV 2 ), and then evolving with Hoppet [17] at arbitrary scales. To be consistent with our computation of the lepton densities, we turn off the photon emission from leptons. We consider two evolution options: "without-P lq " which does not include the quark-to-lepton splitting and "with-P lq " which does. In Fig. 6 we report the results for the electron densities at the scales µ 2 = 10000 GeV 2 and µ 2 = 50 GeV 2 . These plots clearly show the relevance of the the inclusion of the P lq splitting, which is crucial to achieve the NLO accuracy. The effects are especially large in the small-x region where the P lq splitting gets logarithmic enhanced contributions (see Eq. (2.30)) and the omission of the P lq splitting would lead to deviations of order 10% for the scales considered. Ratio without-P lq with-P lq Figure 6: Impact of the P lq splitting function in the evolution of the electron density. The reference curve (black line) is the lepton PDF as directly computed at the given scale µ 2 with the LUX approach. The other two curves are obtained by starting with the lepton PDF directly computed at the reference scale with the LUX approach, and then evolved either turning off (red curve) or turning on (blue curve) the P lq splitting. We show results both for the forward evolution at the scale µ 2 = 10000 GeV 2 (left panel) and for the backward evolution, µ 2 = 50 GeV 2 (right panel).
In Ref. [18] a study of lepton PDFs was performed, and a PDF set with lepton was obtained. This study was carried out before the LUX procedure was available. Large variations were found there depending upon the assumptions on the initial conditions for the photon and lepton densities. Nevertheless, at large factorization scales, these studies should capture at least the order of magnitude of the lepton and photon densities, since they are prevalently generated by perturbative radiation. We have compared our set with the set of Ref. [18], and have found that indeed they are compatible in order of magnitude, with differences that range from 10% in the small-x region, up to 50% for large x. These findings are in line with the fact that there is a contribution to the large-x photon PDF coming from the low Q 2 region (see Fig. 18 of LUX2) that amounts to about 50% of the total, and can only be computed with reliable accuracy by exploiting the electron scattering data as we do.

Phenomenology
The precise determination of the leptonic content of the proton allows us to consider the LHC also as either a (broad band beams) high energy lepton-(quark/gluon) or a leptonlepton collider, even including muons and taus in the initial state, which are beyond the current collider accelerator technology. In the next subsections, after a brief illustration of the associated luminosities, we present some physically motivated applications of leptoninitiated processes at the LHC.

Lepton luminosities
We begin by showing in Fig. 7 the luminosities, defined as that we computed for pp collisions at 13 TeV (left), 27 TeV (middle) and 100 TeV (right) using the LUXlep set. The error bands are obtained with the standard method used for PDFs with replicas. In particular we show lepton-gluon (lg, purple, L l − g + L g l − ), leptonup (lu, orange, L l − u + L u l − ) lepton-photon (lγ, green, L l − γ + L γ l − ) and lepton-lepton (ll, blue, L l + l − + L l − l + ) luminosities. As a reference, we also show the photon-photon luminosity (γγ, red, L γγ ). In the upper panels we only show the central values, since the uncertainty band is too small to be appreciated. In the two bottom panels we show the relative uncertainties, obtained with the usual prescription adopted in sets with replicas (see Sec. 4 and 5). The uncertainties are all very similar and below 5% over a large range of M . Only for very high masses (M/ √ s 0.3) the uncertainties exceed 5%. In the case of the lu and lg luminosities it is clear that the uncertainty is dominated by the uncertainty on the QCD partons. If we compare the ll and lu luminosities we note that the former is suppressed by a factor of about 8 · 10 3 with respect to the latter. Similarly the lγ luminosity is suppressed by a factor of about 200-300 with respect to the γγ luminosity.
Next we show in Fig. 8   ones involving electrons plotted in Fig. 7. As expected, the luminosities of heavier leptons are of the same order of magnitude as the ones involving electrons, but they are somewhat suppressed, in particular at lower masses.

NLO corrections
We remark that, in order to exploit the accuracy of our lepton PDFs, NLO calculations of the lepton-initiated processes are needed. In general, these will involve processes where incoming leptons are replaced by incoming photons splitting into a lepton pair, 2 where the associated collinear singularity is subtracted. These subprocesses are suppressed by a single power of L, i.e. the logarithm of the ratio of the process scale to some typical hadronic scale. This is because the photon PDF is larger than the lepton PDF by a factor 1/(Lα), and the NLO correction carries an extra factor of α and no large logarithms (since the collinear divergence has been subtracted). In the following applications, that are given for illustrative purposes, these higher-order corrections are not included.
We can expect that, at the LHC, lepton-initiated processes may become competitive with other production mechanisms for the search of New Physics objects that have a preferential coupling to leptons. Since the New Physics objects we are searching for are expected in general to be very massive, one may also worry about the fact that our lepton PDFs are computed including a photon exchange, but no Z exchange diagrams. A back of the envelope estimate of these effects would lead to an increase of the lepton PDFs of the order of 5% at TeV scales. On the other hand, our definition of the lepton PDFs without Z contributions can consistently be used, as long as Z exchange effects are included as higher-order corrections, that, due to the Z mass, do not present collinear singularities. Thus, the effect of the inclusion of Z exchange should be considered together with the inclusion of NLO effects, that we are neglecting in the following examples.

Lepton-lepton scattering
Signatures that at the LHC have been considered exotic so far, and important to test flavour violating interactions, are two isolated, back-to-back leptons of different flavours with the same or with opposite charge (see e.g. Refs. [19][20][21][22]). Since our parton densities now include lepton PDFs, we are in a position to estimate the Standard Model (SM) contribution to these signatures coming from → scattering mediated by a photon. These SM processes are accompanied by no other significant activity in the event.
We consider here both 13 and 27 TeV collisions and require standard transverse momentum and rapidity cuts on the leptons, Since the processes we are considering are dominated by a photon exchange in the tchannel, we set the factorization scale to the lepton transverse momentum, and estimate the uncertainty on the cross sections by varying the factorization scale by a factor of 2 up and down. In Tab. 1 we give cross sections for the hard scattering among different lepton flavours for both 13 and 27 TeV proton-proton collisions. Rates for charged-conjugated processes are identical. The large errors are due to the fact that the results that we present here are only LO accurate, but the calculation can be easily carried out at NLO order in QCD. One can see that at the end of the High Luminosity program, with an estimated integrated luminosity of 3 ab −1 , about 850 e + µ − , 550 e + τ − and 500 µ + τ − will be produced and pass our basic lepton cuts. At 27 TeV the cross-sections are almost a factor two higher, and the luminosity is a factor of five higher, so that we expect roughly 10 times more events.
The dominant SM background for same sign leptons comes from W + W + (or W − W − ) plus dijet production. For example, the fully inclusive leading order cross section for W + W + , with the W -bosons decaying into a single lepton species is about 443 fb at 13 TeV. However, if one requires that the leptons pass the cuts of Eq. (6.2), vetoes on jets with p t,j > 20 GeV, requires a missing transverse momentum less than p t,miss = 15 GeV and further requires the two leptons to be balanced in transverse momentum and back to back then the background coming from W + W + decays reduces to about 0.1 ab. This means that not a single event is expected to pass the cuts even at the end of the HL-LHC program. These additional cuts have instead no effect on the cross-sections quoted above. Besides same sign W pairs, one should also consider the very abundant background coming from heavy flavour production (c and b flavoured hadrons) that decay leptonically. In order to get rid of these backgrounds, it is crucial to estimate to what extent one can veto events for additional hadronic activity, that is bound to be much smaller for leptoninitiated processes. This requires the availability of a shower Monte Carlo that can handle incoming leptons.

Z searches
As a second application of our lepton PDFs we consider here the production of Z -bosons. Here we make the very generic assumption that we have a flavour diagonal Z that couples only to leptons. Simple models that can account for that have been put forward in the literature [23]. The resonance cross section is given by where J is the spin of the resonance, s 1 and s 2 are the spins of the incoming particles, k = E/2, E, M and Γ are the energy, mass and width of the resonance, and B 1 , B 2 are the branching fractions of the resonance into the initial and final state. In our case (6.5) In the narrow width limit we have The hadronic cross section is then In the following we consider, for reference, a Z that couples vectorially only to muons and taus (and to the corresponding left-handed neutrinos), that is the least constrained in the model of Ref. [23], and consider the µ + µ − final state. In this case, we have B 1 = B 2 = 1/3, and the production proceeds from both µ + µ − and τ + τ − annihilation. The width is given by We now give a very rough estimate of the production rate and significance of the corresponding signal, assuming an irreducible Drell-Yan background. For the muon energy resolution we interpolate the measured points illustrated in Fig. 9 of Ref. [24], and the reconstruction efficiency times the acceptance was obtained by a rough interpolation of Fig. 2 of Ref. [25]. The significance plot for such an object is shown in Fig. 9 (bottom panels), while the upper panels show the number of expected events. The plots were obtained in the following way. For each Z mass we compute the number of signal events falling in a bin centered around M Z , with a size b w = Γ 2 + r 2 M 2 Z , where r is the muon p T resolution. We use formula (6.7) for the cross section, multiplied by a reduction factor 2 arctan(b w /Γ)/π, to account for the loss on the sides of the Lorentzian peak, and by the reconstruction efficiency times the acceptance. We compute the Drell-Yan background using the code of Ref. [26], and integrate the cross section in the given bin, multiplying also by the reconstruction efficiency times the acceptance factor that we used for the signal. The significance is taken as the ratio of the number of signal events to the square root of the number of background events. Comparing our exclusion plot with Fig. 2 of Ref. [27] (for a direct limit see [28]), we see that in the case of current and High Luminosity LHC operation no relevant limit is found that is better than the one coming from neutrino trident production of muon pairs, that (according to Ref. [27,29]) yields M Z /g 540 GeV, and is represented by the region above the red line in the lower plots of Fig. 9. In the case of the 27 TeV LHC, our method can yield exclusions of regions that are still unexplored.
In the present study we have considered the Drell-Yan background as irreducible. This may not be the case, since the hadronic activity accompanying a leptonic collision is much smaller than the one accompanying coloured parton collisions. At the moment, we do not know of any Parton Shower generator that can reliably generate the initial state radiation from leptonic collisions, although there are no serious obstacles to the implementation of such effects [30]. It is however easy to estimate their size. The hardest radiation accompanying an initial state lepton is another lepton, with average transverse momentum given by A hadronic emission will happen only as next-to-hardest radiation, and thus it will be suppressed by two powers of the logarithm. We emphasize also that in our case, roughly 50% of the time the hardest radiation will be a τ , that may decay hadronically. The tau will give rise to a very collimated jet. Thus, in around 50% plus 50 × 0.35% of the cases the hardest accompanying radiation is leptonic, but also in the remaining fraction of hadronic tau decays one may be able to reconstruct the tau by requiring a very narrow jet. So, even if at the moment the effect of a jet or hadronic veto on the efficiency of the signal selection cannot be estimated reliably with a Monte Carlo, we have good reasons to believe that it may yield a considerable reduction of the background.

Doubly charged Higgs production
From Fig. 9, we see that it is the large Drell-Yan background, rather than the lack of signal events, that limits the reach for the detection of a Z . This suggests that we should turn to signals that are essentially background free. Opposite-sign leptons of different flavours suffer for the presence of a large W + W − background, while same-sign leptons may be considered to be essentially background free. The production and decay of a doubly charged Higgs H ±± via lepton-lepton collisions may give a relevant signal for large enough values of the coupling. Such an exotic particle usually arises in extensions of the Standard Model which aim to accommodate neutrino masses. For example, in the context of a type-II see-saw mechanism [31][32][33][34], a (non-renormalizable) dimension-5 operator added to the Standard Model Lagrangian that can give rise to a Majorana mass term for the neutrino, can be effectively generated by renormalizable interactions with a triplet of scalar particles with SU(2) L × U(1) Y quantum numbers (3,2). The triplet comprises a doubly charged H ±± , a single charged H ± and a neutral component H 0 .
The doubly charged state can couple both to leptons and to W bosons. Therefore, its main production mechanisms are the pair production via an s-channel intermediate Z boson or photon, and the associated production with a singly charged Higgs H ± . Various experimental searches have been carried out by both the ATLAS [20] and CMS [35] collaborations focusing on multi-lepton final states. Typically, these searches assume a scenario in which the doubly charged Higgs decays predominantly into leptons, assuming that the coupling to the W bosons is negligible. The background is usually small since events with two prompt and well isolated leptons with the same electric charge are produced very rarely by Standard Model processes. A bump search is performed for a narrow resonance in the same-sign lepton pair invariant mass, which allows to put constraints on the lower value of the mass of the doubly charged Higgs. The current limits exclude a doubly charged Higgs with mass M H ±± 800 GeV at 95% CL. These searches are insensitive to the coupling between the H ±± and the leptons, and only the leptonic branching ratios matter, provided the coupling is large enough for the decay to occur inside the fiducial volume of the detector.
We consider here the direct resonant production of a single H ±± from lepton-lepton annihilation, whose rate is proportional to the square of the y l 1 l 2 Yukawa coupling. This is complementary to the searches mentioned earlier. We observe that this search strategy for the signal is analogous to the previous study on the Z , so that we can effectively apply the same procedure. At variance with that case, the Standard Model background is drastically reduced, and we assume here that the process is essentially background free. The signal signature is indeed given by a pair of same-sign leptons with no missing energy and very limited activity in the event. We consider, for concreteness, a H ±± that couples only to muons. In fig. 10 we plot the number of detected events in the [M (H ±± ), y µµ ]-plane. We don't consider the mass region below 0.8 TeV, that has already been ruled out by ATLAS and CMS [20,35]. Projections to higher luminosities of these analyses have been given in Ref. [36], and are reported in the figure. From the plot, we see that at the present LHC with 300 fb −1 a few events will be available for masses above the projected exclusion limit of Ref. [36], that corresponds to 1.168 TeV at 139 fb −1 . The same is true at the High Luminosity LHC, where the projected limit is 2.276 TeV, and at the High Energy LHC, where the projected limit is 4.56 TeV. We thus conclude that, for sufficiently large coupling, the s channel production of a doubly charged Higgs may have a mass reach comparable to analyses relying upon pair production. an appealing explanation for tensions in flavour physics. For a review of the various aspects of the LQ physics we refer to the Refs. [37][38][39][40].
As an illustrative model, we consider a chiral R 2 LQ of charge 5/3 (following the nomenclature in Ref. [41]), which couples to the conjugate of the left charged leptons to right u-type quarks. A summary of LQ searches based on both single and pair LQ production can be found in [42]. According to this analysis, based on 36 fb −1 data at 13 TeV, the point m LQ = 2 TeV, y eu = 0.3 in the mass-Yukawa coupling parameter space is still allowed for a LQ which couples only to the quark up and the positron 3 . For the case of a LQ which couples only to the quark up and the muon and for the same value of the LQ mass, slightly larger couplings remain unconstrained, as for example y µu = 0.5. In the following, we consider these two points as our benchmark scenarios. As for the width, we assume it is dominated by the 2-body decay and it is given by neglecting all fermion masses. Having at our disposal a precise determination of the lepton densities, we can investigate the sensitivity reach of a completely unexplored mechanism: the production of a single leptoquark in the s−channel via the process + q → + q. Here, we consider the two subprocesses e + + u → e + + u and µ + + u → µ + + u separately. As for the background, we assume as main source the associated production of a jet and a W boson decaying leptonically. We require that the lepton and the jet are both central, |η| < 2.5, and with transverse-momentum larger than 500 GeV. Furthermore, since the signal is a lepton+jet system balanced in the transverse plane, for the background estimate we veto missing transverse-momentum associated to the neutrino larger than 50 GeV.
In Figs. 11 and 12 we plot the invariant mass of the lepton+jet system and the charged 3 It is not excluded even if one considers the more stringent bounds coming from the recast of the dσ/dp T,e [fb/GeV] pp @ LHC13 LQ WJ Figure 11: Invariant mass of the lepton plus jet system (left) and electron transversemomentum (right) distributions. The signal (red line) is due to a LQ of mass M LQ = 2 TeV and Yukawa coupling y eu = 0.3 in the hypothesis of minimal coupling of the LQ between first-generation quarks and leptons. W (lν) + 1 jet background is shown in blue. The selection cuts are p T,l , p T,j > 500 GeV, |η l |, |η j | < 2.5 and p T,miss < 50 GeV.
lepton transverse-momentum distributions produced in pp collisions at 13 TeV for the positron-up and the antimuon-up processes respectively. A good sensitivity to the LQ is reached with clear peaks in both distributions. For a more quantitative and immediate comparison to Ref. [42], in tab. 2 we report the cross sections and the number of expected events in the lepton+jet mass window 1950 GeV < m j < 2050 GeV corresponding to an integrated luminosity of 36 fb −1 . Our finding is that the signal-to-background ratio S/ √ B is very good, about 6.2 and 16 for the electron-up and the muon-up processes respectively. From this preliminary analysis, the single s-channel LQ mechanism outlined above seems able to reach regions of the parameter space that cannot be accessed with current methods. experimental results on the measurement of the weak charge in atomic systems (see Appendix B in Ref. [42]).

Conclusions
In this work we have carried out a computation of the lepton densities in the proton, up to NLO accuracy. The computation relies upon the good quality of the available data on the electroproduction structure functions and the proton form factors, and is carried out in full analogy with what was done for the photon density in Refs. [1,2]. As in the photon case, the inclusion of NLO corrections is mandatory in order to reach an accuracy at the few percent level. It is in principle possible to increase the precision of our lepton PDFs up to the NNLO-QCD level, as was done in Refs. [2] for the photon.
We have used our result to extend the pdf set NNPDF31 nlo as 0118 luxqed of Ref. [6] with the inclusion of leptons. The set so obtained will be soon made available in the LHAPDF [43] format, under the name LUXlep-NNPDF31 nlo as 0118 luxqed.
With respect to the corresponding calculation in the photon case, some novel requirements have emerged for the lepton PDFs that were not present there. These are better understood if we remind that we have adopted the following rules for the determination of the perturbative order of the calculation: quark densities are of order 1 (actually of order (α s (µ)L) n , with L = log(µ/Λ) and Λ a typical hadronic scale for all n); photon densities are of order αL relative to the quark densities; and lepton densities are of order αL relative to the photon densities, and thus of order α 2 L 2 relative to the quark densities. NLO accuracy requires that terms of order α should be retained for the photons, and terms of order α 2 L should be retained for the leptons, thus adopting the criterion L ≈ α −1 s (µ). According to this counting, leading order electromagnetic plus mixed QED-QCD (of order αα s ) splitting functions are all what is needed to maintain NLO accuracy of the photon density, while in the case of lepton densities a term of order α 2 is also needed. In fact on the right hand side of the evolution equation for the lepton densities we expect terms of the form ∂f l log µ 2 = p lγ ⊗ f γ + p lq ⊗ f q , (7.1) where in the first term we have contributions of order α 2 L and α 2 from the leading and NLO contributions to f γ , and the second term, of order α 2 , needs the inclusion of the α 2 splitting function p lq . This splitting function is not normally included in the electromagnetic evolution of parton densities, but is needed in our case, and we implemented it in Hoppet in order to complete our work. Another aspect is the treatment of higher-order electromagnetic terms, giving rise to contributions of order α 3 L 3 to the lepton PDFs. If we adopt the commonly used rule that α ≈ α 2 s , and thus αL 2 ≈ 1, these terms should be counted as NLO in our calculation. In the framework of the photon PDFs, similar consideration lead to the inclusion of terms of order α 2 L 2 , that in that case arise only from the electromagnetic charge renormalization. In the lepton PDF case, the α 3 L 3 terms arise from collinear radiation of photons from the final lepton. We have devised different methods to include these terms in our calculation.
The inclusion of the lepton parton densities in the proton adds new production mechanisms for Standard Model and New Physics processes at the LHC. Lepton-lepton electromagnetic scattering can give rise to final states with different flavour and/or same sign leptons. We found that these phenomena may be observable at the LHC. Besides being useful to test the underlying theory, these processes may also be used to assess the structure of the underlying event in lepton-initiated processes. In view of the large contribution of the elastic component to the lepton PDF, there must be a substantial fraction of events with large rapidity gaps, containing only the matching leptons of opposite charge arising from the photon splitting process. Even if we let aside the possible presence of rapidity gaps, we know that the hadronic activity in lepton-initiated processes must be greatly reduced. At present, to our knowledge, there are no shower Monte Carlo that can simulate lepton-initiated processes, although we may assume that they will become available in the near future. It is likely that when these tools will become available, it will be possible to optimize methods for rejecting processes initiated by coloured partons (with respect to those initiated by leptons) based upon the accompanying event activity, thereby reducing potential backgrounds to searches targeting lepton-initiated New Physics processes.
In this work we have also considered few applications of the lepton PDFs to basic Standard Model scattering processes and to some selected New Physics processes. We have found that lepton scattering at low transverse momentum (above 20 GeV) is likely to be observable at the LHC, with rates increasing from a handful of events with the current LHC settings, up to several hundreds for the High Luminosity LHC, and few thousands for a High-Energy LHC.
As example of searches of New Physics, we have considered the case of leptonic production of a hadrophobic Z that couples only to muons and taus. By considering the µ + µ − final states, we have found non-negligible production rates, and a relevant significance over the Drell-Yan background in a large region of the mass-coupling plane, although, for the case at hand, neutrino scattering data already excludes a large fraction of this region.
Production of doubly charged resonances coupled to leptons can yield a signal that is essentially background free, consisting of same-sign leptons. We have considered the production of a doubly charged Higgs that couples only to muons. This object can also be pair-produced, or produced in association with a singly charged partner at the LHC, and in this case the signal does not depend upon its coupling to leptons, but only upon its branching ratios. With the current LHC and its High-Luminosity and high-energy upgrades we find that, for sufficiently large couplings, the leptonic production can reach the limits that can be set using pair production.
Incoming leptons colliding with quarks can give rise to s-channel production of leptoquarks. In this case, the rates are higher with respect to the lepton collision processes, since they benefit from the valence density of the quarks. We have considered two benchmark points for the production of a scalar leptoquark coupled either to up quarks and positrons or to up quarks and anti-muons, with couplings h ue = 0.3 for the first case and h uµ = 0.5 for the second case, which are currently not excluded. In both cases we observe a relevant signal over the W + j background, both in the invariant mass of the lepton-jet final state and in the lepton spectrum. This production mechanism is thus very promising, and further investigation in this direction are undergoing [44].
The applications that we have reported here only serve as examples of what could be achieved in the framework of lepton-initiated processes. We remark again that this framework is quite new, and in order to develop it further, extensions of Monte Carlo generators that can handle lepton densities in the proton are needed. Likewise, NLO calculations of lepton-initiated processes should be performed, and implemented in NLO+PS frameworks. As stressed several times in this work, these NLO corrections are of the same order as the typical NLO corrections in hadronic collisions, i.e. from 1 to few 10%, and thus they are necessary in order to achieve a reasonable accuracy. Interfacing them to parton showers is also necessary in order to understand to what extent vetoing over hadronic activity affects lepton-initiated processes, and thus can be used to limit potential backgrounds. In the studies that we have performed, we have found that lepton-initiated processes in proton collisions have the potential to increase the reach of New Physics searches at the LHC, thus justifying the theoretical effort needed for their study.

A Partonic calculation
We need to compute the process shown in Fig. 1 (right), where now q = xp and q 2 = 0. We have where y is the cosine of the scattering angle in the photon-scalar CM frame. The squared amplitude can be written as where the overall 1/2 factor is the spin average for the photon, and two minus signs, one for the fermion loop and one for the photon spin projection, compensate each other. The coupling constant for the scalar has been omitted, since it cancels when dividing by the Born cross section. The real phase space is dy, (A.1) and the (partonic) flux factor is 1 4q · r = 1 2E 2 cm . (A. 2) The Born cross section is where the first factor is the spin average, the second factor is the flux factor, the third factor is the result of the trace and the 2πδ(E 2 cm − M 2 ) is the single particle phase space. The real cross section σ r is computed in four dimensions by standard means, multiplying the square amplitude, the flux factor and the phase space. It has the form σ r σ B = dy A(y, z) + 1 1 + y B(z) , According to the FKS prescription [45], the parton model result is where the second term is the collinear remnant (see Eq.(2.102) of ref [46]). We get The full parton model formula at NLO is thus where z = M 2 /E 2 cm = M 2 /(Sx).

B The effect of the lepton mass
There are cases when the effect of the lepton mass on the lepton PDF cannot be neglected. For the muon, for instance, it turns out that the smallest values of Q 2 that contribute to the lepton parton density is of the same order of the muon mass squared, and for the tau there is a contribution from a range of Q 2 values below the tau mass squared. For the electron, one would be inclined to believe that the mass should not matter. However, we recall that, since the proton has an overall electric charge, at very high energy it carries an accompanying electromagnetic field that can be described as a superposition of virtual photons, that in turn can materialize into pairs of nearly massless leptons. This implies that for very small x the lepton PDF for a truly massless lepton should diverge, due to the elastic contribution.
In the case M 2 (z) = µ 2 F we do not have any restriction on z from the first theta function, and the z integrals are given by In this case, if the condition (C.10) holds, we have z (2) min = ξ, and the first integral vanishes. Thus, in this case, the minimum value of Q 2 is and the integration is written as The second term in Eq. (C.15) is present if the Q 2 integration is interpreted as an oriented integral, and absent otherwise. The z integration of Eq. (2.25) written in terms of the variables illustrated here was easily performed using MAXIMA [47]. We do not report here the lengthy result. In order to evaluate it with sufficient accuracy it must be implemented in quadruple precision in the fortran code, in order to avoid sizable rounding errors.

D The O(α 3 ) term
Terms of order α 3 arise from graphs where one further photon emission is allowed from the lepton, as the one illustrated in Fig. 13. We need to compute the leading double-logarithmic Figure 13: A diagram contributing at order α 3 with up to three powers of log(M 2 /Λ 2 ) to the probe process. term in the leptonic tensor for this process. We have The double logarithmic region requires that |q 2 | M 2 , and from analyticity considerations we infer that L 2 must be proportional to q 2 for small q 2 (up to logarithms), since L µν is the expectation value of a product of currents, and thus cannot have poles in q 2 . Likewise, we must have that in L 1 q µ q ν q 2 + L 2 q · r q µ q ν q 4 , (D.2) the q 2 singularity must cancel. Thus, for small q 2 Therefore, in our limit In order to evaluate L 1 , we notice that the cross section for an almost on-shell photon to inclusively produce the heavy fermion is given in terms of the leptonic tensor as where the first factor is for the spin average, the second is the flux factor, and the third is the contraction of the leptonic tensor with the spin projection. We now write σ using the factorization formula twice where σ B (r, qzy) is the Born cross section σ B (r, qy 1 y 2 ) = πδ((r − qy 1 y 2 ) 2 − M 2 ) = π (−2r · q) δ(y 1 y 2 − z ), (D.7) and we have used z ≈ −M 2 /(−2r · q) in the low Q 2 limit. Thus − L µ µ = 4π α 2π 2 1 2 log 2 M 2 Q 2 P γ (y 1 )dy 1 P ll (y 2 )dy 2 δ(y 1 y 2 − z ). (D.8) We define the convolution of splitting functions as P (2) γ (z ) ≡ dy 1 dy 2 P γ (y 1 )P (y 2 )δ(y 1 y 2 − z ) = 1 z dy 2 y 2 P γ z y 2 P (y 2 ) = 1 z dy 2 y 2 P γ z y 2 P (y 2 ) − In fact, with the same procedure used in this section we could have computed the leading logarithmic term in formula (2.23), the only difference being that we would have had a single logarithmic integral, and thus a single log instead of half a log squared, and a single splitting function to find a lepton in the photon, instead of the convolution of two splitting functions that we computed here. Eq. (D.15) can also be written, performing a change of variables, and using the parton model formula for F 2 , in the form This form is very suggestive, since it is the convolution of three leading order splitting functions combined with the (ordered) integration of three intermediate scales. In fact, this equation can be obtained by iterating the first three rungs of the integral form of the leading-order Altarelli-Parisi equation, under the assumption that f (1) vanishes when µ F ≈ m p . This is certainly the case, since f (1) is of order α 3 without logarithmic enhancement in this limit. The calculation described in this appendix could be used to compute directly the leptons PDFs at any (perturbative) scale, without making use of step 4 of section 4. It is therefore interesting to compare what we obtain with this method with respect to our default one. The comparison is shown in Fig. 14 Figure 14: Comparison of the different methods to include the effects of order α 3 , see text for more details. of our main formula for the lepton densities, Eq. (2.25), with (black) or without (red) the terms of order α 3 given in this appendix. We show results at three different values of the factorization scale, and normalize them to our default procedure, described in item 4 in Sec. 4, corresponding to the central value of the LUXlep set. We have implemented the α 3 contribution using three different variations (denoted as "types" in the figure), that should all agree up to higher-order corrections. In all the three types we have used the integration limits of formula (2.12). In type 1 we use formula (2.12) as is, in type 2 we replace one of the two powers of logarithms in Eq. (2.12) with 18) and for type 3 we make the same replacement for both powers of the logarithm. This form of the logarithm is the one that we have in formula (2.25). Of course, without making a more detailed calculation we cannot know the true form of the argument in the logarithm.
We use this artifact just to gauge the sensitivity to subleading effects. First of all, we notice that the α 3 correction, irrespective of the method used to compute it, does not seem to be as important as the NLO correction included in formula (2.25), that is suppressed by the absence of a large logarithm (which in our counting is equivalent to a power of α s ). While the latter is of order 10%, the correction that we computed here is not much larger than 1%, if one excludes the very large x region. All the four methods that we have considered to include α 3 effects (i.e. type 1, 2 and 3 plus our default method) seem to reduce the density evaluated with formula (2.25) at large x. This is consistent with the expectation that the electromagnetic radiation from leptons should be particularly important at large values of x, where it tends to soften the lepton distributions. In all cases, our default method seems to be more effective in doing so than the calculation performed in this appendix. This is perhaps due to the fact that mixed QED-QCD splitting kernels are included with our default method, but not in types 1-3. We also notice that subleading terms, estimated as the difference between types 1, 2 and 3, are not negligible if compared to the full magnitude of the effect (which is the difference with respect to the result without α 3 effects) being several tens of percent of the overall α 3 effect. The same can be said about the difference of the default result with respect type 1, 2 and 3. There we observe a deviation that is also not negligible in the small x region.
The variation between the default method and type 1, 2 and 3 could be added to our calculation as a further source of uncertainty arising from higher-order effects. Since the uncertainty that we find here is considerably smaller than the one obtained with the method described as item (HO) of Section 4 we do not include it, also reassured by the fact that even if we did include it, the final error estimate would not change substantially.