LePDF: Standard Model PDFs for high-energy lepton colliders

The emission of collinear radiation off an elementary lepton can be factorised from the hard scattering process by introducing Parton Distribution Functions of a Lepton (LePDF), which, contrary to protons, can be derived from first principles. In case of multi-TeV lepton colliders, such as the muon colliders currently being proposed, the complete structure of Standard Model interactions must be taken into account. In this work we solve numerically the corresponding DGLAP equations at the double-log order and provide public files with LePDFs for both muons and electrons, including polarisation effects. We discuss several interesting aspects of the resulting PDFs and compare them with the Effective Vector Approximation, showing that the latter fails to describe well the vector bosons PDFs at small momentum fractions, unless it is extended to higher orders.


Conclusions 23
A Inputs and Formalism

Introduction
At present, the Large Hadron Collider at CERN is our only tool for direct exploration of physics at the electroweak scale and above and the high-luminosity phase is planned to last until the early 2040s. It proved to be a formidable machine for both searches of new heavy particles as well as precision studies at the electroweak scale, of which the Higgs discovery and the precision study of its coupling is a prime example. Nevertheless, the next large-scale experiment in high-energy physics is likely to be a lepton collider. The proposed options include circular or linear electron-positron colliders (FCC-ee [1] and CEPC [2] for the former, ILC [3], CLIC [4] and CCC [5] for the latter) as well as µ + µ − colliders [6][7][8][9][10]. Among these, the linear e + e − colliders and the muon colliders could achieve multi-TeV center-of-mass energies.
While leptons are elementary particles in the Standard Model (SM), the process of collinear emission of initial state radiation (ISR), with transverse momentum p T much smaller than the energy of the hard scattering process, can be factorized and a description in terms of parton distribution functions (PDFs) can be introduced, similarly to what is done in case of proton colliders and the parton content of a proton. The case of collinear photon emission from an electron is known since almost a century and at leading order can be described using the effective photon approximation (EPA) [11][12][13][14] f EPA where E is the energy of the initial electron and P γe (x) = 1+(1−x) 2 x is the splitting function, that describes the probability of an electron to emit a photon with a fraction x of its energy and virtuality −p 2 T /(1 − x). When E ≫ m e the large logarithms can be resummed in order to improve the perturbative expansion, and the factorization scale Q is introduced. Invariance of the physics under the factorization scale leads to the DGLAP equations [15][16][17]. For a generic splitting of massless partons A → B + C (as in Fig. 1), choosing p T as factorization scale and working at the leading logarithm (LL) order one has where P C BA (z) are the leading order (LO) splitting functions (listed in App. B), α ABC (Q) the corresponding running coupling, and the term P v B f B describes virtual corrections (see App. C for details). In case of a proton, due to its non-perturbative nature, the initial conditions for this system must be fitted from collider data. For a lepton, instead, the initial condition can be computed perturbatively and the system can be solved from first principles. The initial condition is for initially unpolarised beams 1 , while all other PDFs vanish for Q 2 = m 2 ℓ at this order. Next-to-leading order (NLO) corrections to the initial conditions have also been computed [18] and become relevant when next-to-leading log (NLL) evolution is considered [19,20]. In this work we limit ourselves to LL evolution and LO initial conditions.
In case of multi-TeV lepton colliders one can be interested in factorization scales much higher than the electroweak scale. In this case all SM interactions and fields should be considered [21]. In this aspect, lepton colliders differ qualitatively from hadron colliders. For the latter, QCD interactions are the dominant contributions to the DGLAP evolution in the whole energy range (see however e.g. Refs. [22][23][24][25][26] about the photon and lepton content of the proton). In the case of lepton collider, instead, QED and electroweak (EW) interactions are the leading ones, with QCD playing an important but not dominant role.
The facts that SM gauge group is non abelian, that it is spontaneously broken at the electroweak scale, and that interactions are chiral have several crucial implications for the evaluation of collinear radiation in this regime. The non abelian nature of EW interactions imply a lack of cancellation of infrared (IR) divergencies between virtual corrections and real emission, which generates Sudakov double logarithms [27][28][29]. Electroweak symmetry breaking effects have been shown to provide important contributions and to be the dominant ones in case of longitudinal polarisations of electroweak gauge bosons [30]. Since the SM interactions are chiral, PDFs become polarized above the EW scale [31].
The goal of this work is to numerically solve the DGLAP equations from the initial condition at Q = m ℓ up to multi-TeV scales, taking into account all SM interactions (including the effects mentioned above), and to provide public results with the complete LePDFs. For concreteness, in the following we assume the initial lepton to be a muon, since muon colliders could achieve higher energies, for which our discussion is more relevant. However, all results can be equally applied to e + e − colliders with suitable substitutions thus we provide numerical results for both.
A preliminary study of PDFs for muon colliders by us, mainly focusing on the fermionic degrees of freedom, was included in Ref. [32]. In this work we extend it by including all SM interactions, Sudakov double logs, polarisation, and EW symmetry breaking effects. One similar study has already been performed in the literature, specifically in Refs. [33,34], and we compare our results with the plots presented in those works. We also compare against approximate solutions obtained by solving iteratively, at fixed-order, the DGLAP equations, which provide the analogous of the Effective Vector Boson Approximation (EVA) [35][36][37]. A noteworthy result of this comparison is the realisation that the LO EVA for transverse EW gauge bosons PDFs is insufficient for correctly describing the full result, which is instead well approximated by including contributions up to O(α 2 ). Crucially, Sudakov double logarithms appear at this order due to the gauge boson splitting V → V V and the virtual corrections to the muon PDF.
In Section 2 we present the first part of the evolution below the EW scale, where only QED and QCD interactions are relevant. In Section 3 we discuss the main aspects of the evolution above the EW scale. Our results are collected in Section 4, where several notable features of LePDFs are showed, and a comparison with EVA is presented. We conclude in Section 5, while many details of our computations, the numerical implementation and the formatting of our LePDF files are collected in several Appendices. The numerical results for the LePDFs can be downloaded from GitHub at the link: https://github.com/DavidMarzocca/LePDF.

QED and QCD evolution
For factorization scales below the EW scale the relevant degrees of freedom are light quarks and charged leptons, with vectorlike QED and QCD gauge interactions. Neutrinos, while having negligible masses, become relevant only above the EW scale where the W boson can go on-shell. 2 Because the initial condition and the evolution equations are vectorlike, in this regime no polarisation effects are induced, i.e. PDFs will be the same for both fermion chiralities or gauge boson polarisation.
In the DGLAP evolution from the muon mass up to the EW scale one encounters several mass thresholds for each fermion species as well as at the QCD scale Q QCD . At each threshold a matching should be performed. For our purposes, we take all fermions except bottom and top quarks to be massless. 3 The Q QCD scale sets the onset of QCD interactions, which become relevant after the γ → qq splitting. This can be interpreted as the QCD structure of a photon, and can be divided into a perturbative and a non-perturbative component, mainly due to the photon mixing with QCD vector mesons [41][42][43][44][45]. The latter one will be power-suppressed at the large scales we are eventually interested in, so we neglect it. The choice of Q QCD depends on how many resonances are included in the non-perturbative component and a value close to m ρ has been argued to provide a good benchmark [43][44][45]. In practice, we follow the prescription of Ref. [34,43] with Q QCD = 0.7 GeV. 4 Therefore, from m µ to Q QCD we consider only QED interactions, including all charged leptons ad well as the light quarks (u, d, s, c). At Q QCD we match the PDFs and continue the evolution up to the bottom mass scale adding also QCD interactions and setting the 2 A possible effect of neutrinos even below the EW scale is due to neutrinos from muon decay µ − → e −ν eνµ which present IR singularities in the physical region when scattering with the incoming µ + (or viceversa) [38]. In case of muon colliders such singularities are cutoff by the finite width of the muon beam [39,40]. We neglect such effects with our PDF formalism, assuming that it can be described independently of PDFs. 3 In future versions we plan to add also the τ and charm quark mass thresholds. 4 We study the dependence of our results on QQCD by running the evolution also for values QQCD = 0.52(1.0) GeV and interpreting the differences as theory uncertainty on our final results due to nonperturbative QCD dynamics, see dedicated discussion in Sec. 4.6.
initial condition for the gluon PDF as f g (x, Q QCD ) = 0. At Q = m b we perform another matching and continue the evolution up to the EW scale Q EW including also the bottom quark, setting f b (x, m b ) = fb(x, m b ) = 0 as its initial conditions. Given the C and P symmetries of QED and QCD, and the fact that all fermions except for bottom and top quarks are taken massless, below Q EW several PDFs are related: The DGLAP equations, according to Eq. (1.2), are then given by where ℓ = {ℓ sea , µ}, ⊗ indicates a convolution as in Eq. (1.2) and we defined the evolution variable t as such that we start the evolution with the initial condition in Eq. (1.3) at t = 0. The splitting functions are listed in App. B (we regulate the (1 − z) −1 poles using the standard +-distribution, see Eq. (B.1)), while the values of the virtual coefficients P v B can be found in App. C. Finally, the details for the RG evolution of QED and QCD couplings are reported in App. A.1. In Fig. 2 we show the result of our numerical solution of DGLAP equations for a muon, evolved from the muon mass up to the EW scale Q EW = m W (solid lines). As can be seen, at small x the gluon PDF becomes rather important, while quark PDFs have a size similar to the muon itself or sea leptons, as already showed in Refs. [32][33][34]. PDFs of a muon evaluated at the EW scale Q EW = m W . The colored regions corresponds to the variation obtained by modifying Q QCD between 0.5 and 1 GeV (see Sec. 4.6 for details). Here ℓ sea is defined as in Eq. (2.1).

Iterative solution for QED
By solving iteratively the DGLAP equations order by order in (α γ t) = (α γ log Q 2 /m 2 Q ), one gets approximate solutions for the PDFs. The LO contribution for the photon PDF is given by Eq. (1.1) with the crucial substitution E → Q. By including terms up to O(α 2 t 2 ) we get where t is defined in Eq. 2.3 and the I ABBC (x) ≡ 1 x dz z P X AB (z)P Y BC (x/z) integrals are collected in Eq. (B.7). Going up to O(α 2 t 2 ) is required in order to describe the low-x behavior of valence and sea lepton PDFs, that is dominated by the γ → ℓ + ℓ − splitting (i.e. the I f V V f (x) term above). We find good agreement between these expressions and our numerical results.
The gluon PDF starts formally at O(α 2 γ α s t 3 ). However, due to the large value of α s at low energies and the importance of a careful treatment of the matching at the Q QCD , a fixed-order solution would not be a good approximation for the correct result. For the same reason we do not report here also the approximated expression for quark PDFs, even if they formally arise already at O(α 2 γ t 2 ), since they receive large a QCD contribution from gluon splitting into qq. A correct description of the quark and gluon PDFs in a lepton therefore motivates a complete numerical solution of the DGLAP equations [34]. Table 1. Independent degrees of freedom for our DGLAP evolution above the EW scale.

DGLAP evolution in the SM
For energies above the EW scale, the splitting processes in the initial states can involve all SM interactions and fields, which must then be included in the DGLAP equations. The chiral nature of EW interactions induces polarisation effects on PDFs, so in this region all gauge bosons and fermions polarisations are treated separately. Splitting functions for SM interactions in the unbroken phase have been computed in Refs. [21,30,31,[46][47][48], with which we agree. Another well known effect in EW PDFs is the interference between photon and transverse Z T , which must be described with a Z/γ mixed PDF, as well as between the Higgs boson and the longitudinal Z L , which induces a h/Z L mixed PDF [21,30,49]. Also, the non-abelian nature of EW interactions induces Sudakov double logarithms, which must be resummed when one is interested at high energies [27][28][29].
Since we neglect all fermion masses except than the top and bottom quarks, several PDFs are related: (3.1) The fact that PDFs for right-handed fermions and their conjugate are different is due to the induced polarisation effects in the gluon, photon, and Z boson PDFs. The independent degrees of freedom are given in Table 1 in the case of PDFs of a muon beam, the final count is 42 independent PDFs. Regarding the mass thresholds, we include each degrees of freedom right at the corresponding mass scale.

Electroweak symmetry breaking effects
EW symmetry breaking effects can be classified as either due to non-vanishing parton masses or as new contributions to the splitting functions that vanish in the limit of unbroken symmetry [30]. The latter are the so-called ultra-collinear splittings, which are particularly relevant for the longitudinal polarisations of EW gauge bosons [30,50], providing the leading contributions for their PDFs even when Q ≫ Q EW . In order to easily keep track of all these effects near the EW scale it is convenient to work in the broken phase and with the mass eigenstates. Then, as the DGLAP equations remain the same also for higher energies, we remain in the broken phase for any Q > Q EW . We compute the splitting functions in the Goldstone equivalence gauge, following Ref. [30]. However, we checked that the results are the same as the ones in [50], in which the theory is formulated in a standard R ξ -gauge, with a 5-dimensional polarisation vector for the longitudinal component of gauge bosons to take into account the Goldstone contribution.
In case any of the particles involved in the splitting A → B +C is massive, the kinematics of the process and the ensuing splitting functions are modified. The relation between the p T and the virtuality of the parton B entering the hard scattering process becomes wherez = 1−z and m 2 /E 2 or p 2 T /E 2 terms can be neglected in the regime where factorization can be applied. In the DGLAP equations, this modified propagator of the virtual parton B effectively corresponds to a rescaling of the splitting functions as [30] Mass effects in matrix elements that are also present in the unbroken (i.e. massless) theory are instead suppressed by powers of m/E, therefore negligible.
Regarding the ultra-collinear splitting functions, the main feature of these contributions is that they do not scale logarithmically with Q 2 but at large factorization scales are suppressed as v 2 /Q 2 . Nevertheless, they provide important contributions that accumulate during the DGLAP evolution in the region Q ∼ Q EW and then remain almost constant at higher scales, see for instance the W L PDF in the right panel of Fig. 3. We report the full set of splitting functions in App. B, while the list of SM DGLAP equations can be found in App. D.
In principle, with massive partons we should care of kinematics bounds: the splittings described by DGLAP equations involve partons A, B and C with energies zE, xE and (z − x)E respectively, all fractions of the energy of the beam E. Since the particle C is emitted on-shell, we need E C ≥ m C , that is z ≥ x + m C E . This means that the lower extreme of integration in Eq. (1.2) should be modified. However, as in Eq. (3.2), m/E terms can be neglected and we can safely start the integration from x.

Electroweak double logarithms
The fact that the initial and final states are EW non-singlets has important implications, even for inclusive processes. The Bloch-Nordsieck theorem [51], that guarantees cancellation of IR divergencies between real emission and virtual corrections in such processes, is violated for non-abelian symmetries, which implies the presence of Sudakov double logarithms. While in the QCD case this effect vanishes upon averaging over color of the initial states, for SU (2) L we do not take such average and the initial state breaks explicitly the symmetry, hence double logs do appear also in inclusive processes [27-29, 49, 52]. In our context they can be seen appearing in the terms of the DGLAP equations containing a 1/(1 − z) pole and can be made explicit by introducing an IR cutoff in the integral, which allows also to resum the Sudakov double logs related to ISR [46-48, 53, 54]. 5 We implement this following Ref. [47], by modifying the boundaries of the integral in Eq. (1.2) as plays the role of an explicit IR cutoff for the 1/(1 − z) poles, and is set equal to 1 (in which case we use the +-distribution to regulate the divergence and perform the numerical evaluation) except for the cases where the soft divergence is not cancelled between real emission and virtual corrections, in which case we set z ABC max (Q) = 1 − Q EW /Q. This modifies the computation of the virtual corrections, that becomes (see App. C) where the ellipses include other finite contributions to the integral, as well as other interactions. In the two parentheses, the first term is due to real emission of a W ± T boson, while the second is the corresponding virtual correction. The fact that the muon and neutrino PDFs are different is an explicit breaking of SU (2) L and implies a non-cancellation of the pole for z → 1 inside P V f f (z), which in turn generates the double log. To see this explicitly, let us isolate on the right-hand side only the terms that are singular in z → 1, fixing z = 1 everywhere else: where ∆f L 2 ≡ f µ L − f νµ . Upon integration in Q 2 , a log 2 Q 2 /Q 2 EW contribution is generated, that tends to deplete the muon PDF and enhance the neutrino one. This example also clearly shows how no such double log is generated for photon or Z T emission from a fermion or W , since the real and virtual contribution would be proportional to the PDF of the same parton. A Sudakov double log is instead expected for splittings such that the splitting function is divergent in the soft limit, z → 1, and the A and B partons are different: otherwise we put z max = 1 and employ the standard +-distribution (see Eq. (B.1)) to regulate the z → 1 divergence, for a more stable numerical evaluation of the DGLAP equations.
In practice, this happens for W ± T emission off any parton (in correspondence to the poles in the P V f f , P V V V , and P V hh splittings) and for Z T boson emission from an initial longitudinal Z L or Higgs (due to the P V hh splitting), since Z emission changes Z L ←→ h. Analogously, for ultra-collinear splittings this takes place for any W L emission and for Z L emission off and initial Higgs or Z L .
This procedure amounts to a double-logarithmic (DL) approximation to the Sudakov factor for initial-state radiation. This could be further improved to LL or NLL resummation by suitably modifying the scale at which the coupling constant α 2 (Q) is evaluated, as discussed in Refs. [31,60]. Since we are interested in energies where α 2 log 2 Q 2 /Q 2 EW ∼ O(1) but α 2 log Q 2 /Q 2 EW ≪ 1, we limit ourselves to the DL approximation for electroweak corrections in the present work.

Top quark as parton
When the energy of the hard process Q is much larger than the top mass, processes with a collinearly emitted top quark can develop a logarithmic enhancement proportional to log Q 2 /m 2 t . It can therefore become useful to resum these logarithms by including the top quark among the other partons [61,62]. The question of whether or not one should include it as a parton depends on the process considered and on the optimal way to rearrange the perturbative series [63,64]. With this in mind, we provide two versions of our PDFs of leptons, one in the 5-flavour-scheme (5FS) and one in a 6FS, where the top quark is added in the DGLAP evolution for scales above m t . While codes that include the top quark in proton PDFs assume it is massless, in our approach we keep a finite top mass in the same spirit in which we keep finite W and Z masses for the weak bosons PDFs. This is justified by the fact that in our case, contrary to proton PDFs, EW interactions and EW symmetry breaking effects are crucial, and the top mass is one of such effects. For a detailed discussion of different schemes for the top mass in the computation of hadron collider observables with a top quark PDF see [64].
The DGLAP equations for the t L ,t L , t R ,t R are reported explicitly in Eqs. (D.11, D.12, D. 19, D.20). We checked numerically that the dominant contributions are those from initial transverse gauge bosons, with electroweak bosons, photon and gluon terms being approximately of similar size. Instead, ultra-collinear contributions are practically negligible. Nevertheless, in our numerical evaluation we keep all terms.

Effective Vector Boson Approximation
The EPA has been generalized to describe EW gauge bosons in high-energy collisions since the '80s [65][66][67][68], in what is now known as the Effective Vector Boson Approximation (EVA) [36,37]. When the hard scattering energy is much larger than the EW gauge boson mass and the p T of the collinearly emitted gauge boson, then the cross section of the process can be factorized into an almost on-shell collinear emission and the subsequent hard scattering [35,[69][70][71].
One can compute the EW gauge bosons PDFs by evaluating the DGLAP equations at fixed order, similarly to what we did in Section 2 for QED, using f and analogously for the Z and Z/γ PDFs The muon mass here serves as an IR cutoff for the logarithm in the transverse case to cure the x → 1 limit, while we neglect it in the other terms. Notably, the W + has no contribution at this order.
In Fig. 3 we show the dependence in p 2 T of the integrands (left), and the resulting PDFs (right), fixing a value x = 0.3 and showing separately the two polarisation of the transverse W − ± . One can see that the integrands are peaked before the EW scale and, while in the case of W T it decreases as ∼ 1/p 2 T inducing the logarithmic grow of the PDF, for the longitudinal W polarisation the contribution to the integral is localised in p 2 T before the EW scale and the PDF tends to a constant at large scales.
In our numerical integration of DGLAP equations, the effects due to EW interactions are introduced only above the Q EW matching scale. Since we employ p T as factorization scale, this effectively corresponds to performing the integration in Eqs. (3.9,3.10) only for To address this issue we match the gauge bosons PDFs at Q EW to the analytically computed one for the same scale and use it as boundary conditions, and then continue the integration numerically to higher scales. The boundary conditions for the PDFs of other heavy states (h, h/Z L , top quark) are instead set to zero at the corresponding mass scales.

Results
Here we discuss several aspects of our LePDFs. The details of our numerical implementation of the DGLAP equations are collected in App. E. Fig. 4 (top panel) collects, as an example, a set of PDFs evaluated at the scale Q = 3 TeV. One first thing to notice is that, as expected, for x ≳ 0.5 the muon PDF dominates, while for smaller x the largest PDF is the photon one. However, the transverse negative W − T PDF is only a factor ∼ 2 smaller and the transverse Z T boson is another factor of 2 smaller than that: they both receive contributions from the emission off an initial muon. Analogously, the muon neutrino ν µ has a large PDF at large x values due to the emission off a µ − L , which has also a Sudakov double-log enhancement. The positive transverse W + T PDF is instead more suppressed because its leading contribution arises from the emission off the muon neutrino and off another gauge boson. The importance of EW gauge bosons PDFs reflects the common lore that a high-energy lepton collider is also a weak boson collider.
In the bottom panel of Fig. 4 we show the PDFs for the longitudinal polarisations of EW gauge bosons and the Higgs, evaluated at scales of 3 TeV (solid lines) and 30 TeV (dashed lines). The PDFs for W − L and Z L are mostly scale independent, since they receive  the dominant contribution from the ultra-collinear splitting off a muon. On the other hand, the ultra-collinear contribution to the W + L PDF comes mostly from the muon neutrino, which has a PDF suppressed with respect to the muon one. Therefore, other contributions from standard splitting functions (e.g. from P h hV and P V hh ) are sizeable and induce a scale dependence. In case of the Higgs boson there is no ultra-collinear contribution from massless fermions, so one does not expect ultra-collinear terms to dominate and indeed its PDF shows a large scale dependence.
The fraction of the momentum carried by each of the partonic components is given by Table 2. Fraction of the momentum carried by each parton at Q = 3, 10, 30 TeV. the n = 2 Mellin transform of the PDF (see App. C for more details) and it can be used to evaluate the relevant role of the various individual components at different scales. To this end, in Table 2 we give three examples for the set of PDFs shown in Fig. 4 at scales 3 TeV, 10 TeV, and 30 TeV. We observe that as the energy of the hard process is increased the percentages of both the left-and right-handed muon components are decreased and those of all other partons increased, which illustrates the importance of electroweak interactions at higher energies.
In Fig. 5 we show some examples of parton luminosities for a 3 and 10 TeV muon colliders where, unless specified, we sum over polarizations. Parton luminosities can be useful for computing cross sections integrated over angular variables. In case of a muon collider they are defined from the convolution of the PDFs of parton i from the muon and parton j from the anti-muon, as follows: where √ s 0 is the collider center of mass energy and √ŝ is the invariant mass of the two-parton system. From Fig. 5 we can notice that, even at small invariant masses, the µμ luminosity is much larger than the W + W − luminosity (also the charged-current µ Lνµ luminosity is sizeable). The impact of this channel in VBF studies at muon colliders should therefore be studied in more details. It is also interesting to point out that the QCD-related luminosities (gluons and quarks) are very small, which is going to strongly suppress QCD-induced backgrounds in electroweak processes.

Polarisation
The chiral structure of SM interactions above the EW scale induces polarisation effects for the PDFs [31]. In Fig. 6 we show polarisation ratios for several PDFs at a scale Q = 3 TeV. The observed behavior can be easily understood as follows. The W − T , W + T , and Z T PDFs receive the dominant contribution from the emission off an initial µ − L,R or ν µ . Since the P f V + f L splitting function goes to zero for z → 1, while P f V − f L tends to a constant, the positive helicity of the EW gauge bosons will be suppressed for x → 1. In case of the photon, the leading contribution comes from the muon splitting and it is vector-like at leading order, so the polarisation effect will be suppressed.
In case of fermions, left-handed chiralities (and their conjugate) receive contributions from W bosons splitting to ψ Lψ ′ L , therefore their PDF is expected to be larger than the right-handed counterparts. The tendency increases with x in the case of the leptons (except for the muon) and down-type quarks as opposed to antileptons and up-type quarks, since the left-handed parts of the former receive contributions from W − T , while the latter from W + T , and the W + T PDF falls faster than the W − T PDF at high x (see Fig. 4) This effect can be of O(1), since the W − PDF is comparable in size to the photon one. The b L PDF is further enhanced compared to b R due to the W − L → b LtR splitting proportional to y t .

Comparison with the Effective Vector Boson Approximation
In Fig. 7 we show our results for EW gauge bosons PDFs, compared with the LO EVA result discussed in Section 3.4 (to reduce the number of lines plotted we show the sum of transverse polarisations). Several things can be noticed.
For Z/γ, the EVA result is two orders of magnitude smaller than what we find with the numerical evolution. This is due to the fact that Q Z µ L + Q Z µ R = − 1 2 + 2s 2 W ≪ 1 is accidentally suppressed (indeed, it becomes zero when evolving the Weinberg angle at a scale of about 3.6TeV). This cancellation takes place because in EVA it is assumed that the initial-state muon is not polarised. However, in the evolution from the EW scale upward, electroweak interactions induce a polarisation of the muon PDF, which becomes up to ∼ 40% at a scale of 3 TeV, as can be seen in Fig. 6 (left panel). Therefore, in the full numerical evolution there is no such tuned cancellation in the Z/γ PDF. This clearly shows that the EVA result is not reliable for this PDF.
For the longitudinal polarisations, the EVA provides a good description of the PDFs, to within ∼ 10% accuracy. In case of the transverse W and Z polarisations, instead, there is a noticeable discrepancy which grows even up to O(50%) at multi-TeV scales for small x values. This is dominantly due to the missing contributions from V → V V splittings, that start at NLO. Such contributions become important due to two effects: the PDF of the initial-state gauge boson at small x is much larger than the muon PDF and they include IR Sudakov corrections that induce a parametric dependence as α 2 log 3 (Q 2 /m 2 W ). We perform two checks to verify this. First, we do a run of our numerical code setting to zero the P V V splitting functions (both in real emission and radiative corrections): the resulting EW gauge bosons PDFs agree well with LO EVA. Second, focussing only on the W − + PDF for simplicity, we compute iteratively the O(α 2 ) contributions to the weak bosons PDFs, adding the real emissions from P V V V and P f V f splittings, and the corresponding virtual  corrections. In practice, we take the DGLAP equations for transverse EW gauge bosons in Section D.3 and use the O(α) results for the gauge bosons PDFs (i.e. the LO EVA of Eqs. (3.9-3.13)) and the muon. For simplicity, at this step we use approximate expressions for the LO EVA by keeping only the log Q 2 m 2 W,Z term, which makes our NLO EVA result reliable only for Q ≫ m W and x ≪ 1. We also neglect contributions from longitudinal modes and ultra-collinear ones. We then perform analytically the convolution with the splitting functions and finally the integral from m W up to the factorization scale Q: It can be noted that Sudakov double logs appear here in the virtual contributions to f , as well as in the P V + V + terms from the neutral gauge bosons.
In Fig. 8 we show the relative deviation of the LO (dashed) and NLO (dot-dashed) EVA results from the complete numerical PDF as function of x (left panel) and as function of the scale (right panel). We observe that the NLO EVA result improves substantially the agreement with the full numerical result, while the LO EVA has large deviations at small x. The missing terms in the LO EVA become more and more important with larger scales, confirming the argument made above.
In Fig. 9, instead, we plot the W − + W + − (left panel) and W − L W + L (right panel) parton luminosities for a 10 TeV MuC. We show a comparison between the LePDF result (solid lines) and the LO EVA expression in the Q ≫ m W approximation (dot-dashed), that is the one implemented in Ref. [37], or with the complete W mass dependence as in Eq. (3.9) (dashed). We see that, at the level of luminosity, the LO EVA with the complete mass dependence provides a good approximation of the resummed LePDF result up to ∼ 15% deviations for the transverse modes. This means that the much larger deviations we observe at the PDF level for small x are diluted when the luminosities are calculated. On the other hand, the massless approximation deviates up to O(1) even at the luminosity level and in particular at partonic center of mass energies of few hundreds of GeV, where the weak-boson fusion process cross sections are the largest.

Muon neutrino PDF
The leading contribution to the muon neutrino PDF arises already at O(α 2 ) from the µ L → W − ν µ splitting, which presents an IR soft divergence that is cutoff by the W mass. As already discussed in Sec. 3.2, the missing counterpart of this divergence in the virtual correction is at the origin of the Sudakov double log. In the same spirit as done for the EVA, we can compute the neutrino PDF by iteratively solving the DGLAP equations up to O(α), using the zeroth-order expression for the µ L PDF: The Heaviside theta is the result of the IR cutoff z max = 1 − m W Q in the integral. Integrating this differential equation from m 2 W up to Q 2 we get: We observe here the single logarithm due to the standard collinear divergence, while the Sudakov log is absent because the initial muon PDF at zeroth order is just a delta function. It will appear, however, in the computation of an inclusive cross section upon integration of the PDF, due to the x → 1 divergence inside P V f f (x), that is not cancelled by a virtual correction at the same order but is instead cut off at the W mass by the theta function.
In Fig. 7 we show a comparison between this O(α) approximation of the muon neutrino PDF (dashed gray line) with the one from LePDF (solid gray), at a scale Q = 3TeV. We obseve a very good agreement at large x values, i.e. where the µ L → W − ν µ splitting dominates. At smaller x values the O(α 2 ) contribution from Z → ν µνµ will instead dominate. Comparing to the results of Refs. [33,34], we observe a different behavior of the ν µ PDF for x → 1: while both our analytic result described above and LePDF show a cutoff (due to the IR m W cutoff) at x ≲ 1, in their result the ν µ PDF increases similarly to the muon PDF up to x = 1.

Mass effects
We showed in Section 3.1 that massive particles modify the virtuality of the particle B as in Eq. (3.2). As already discussed in [30,36], the impact of the latter effect is important since, due to the presence of the masses in the denominators of the DGLAP equations, the PDFs are lowered or enhanced when m B,C ̸ = 0 and m A ̸ = 0 respectively. In Fig. 10 we show the PDFs computed keeping and neglecting the masses in p T , still starting the evolution of a given massive parton at the scale corresponding to its mass.
For instance, in case of EW gauge bosons (both transverse and longitudinal) one can expect that mass effects lower their PDFs, with the biggest effects at small x due to the (1 − x) factor appearing in front of the mass corrections. This can also be verified from the LO EVA in Eqs. (3.9-3.13). In case of the Higgs, its interactions have the form which explains why the Higgs PDF is bigger when masses are neglected.

Top quark PDF
In Fig. 11 we compare the results between the 5FS and the 6FS showing the top PDFs together with the PDFs mostly affected by the inclusion of the top. As already mentioned in Section 3.3 the top PDFs are mainly driven by the collinear emission off the transverse gauge bosons and the differences are more noticeable in high energies (we use the benchmark Q = 3 TeV). Since the PDFs of W − T and W − L are large compared to the Z and W + T , we expect a larger PDF fort L andt R than t L or t R . Also, the same splittings W − T → b LtL and W − L → b LtR will induce a large enhancement of the b L PDF compared to b R in the 6FS. Additionally, we observe that the PDFs of the EW transverse gauge bosons themselves are almost unaffected, since they are dominated by splitting off a muon. The gluon PDF, instead, receives a noticeable further contribution. Shifts of comparable size are also induced in the PDFs of the longitudinal gauge bosons (and even smaller for the Higgs), but in this case the PDFs are decreased due to a mass effect similar to the ones discussed in the previous Section 4.4.

Uncertainties
Here we discuss several sources of uncertainties in our computation. Some are physical, such as the choice of Q QCD or missing higher orders, while others are intrinsic in the numerical implementation of the DGLAP equations and can be improved simply by dedicating more computational time to the task.

QCD matching scale
As already discussed in Section 2 the QCD scale Q QCD is not clearly determined and different choices of this parameter can have a non negligible impact on the PDF, in particular for the colored particles. To study the dependence of our results on Q QCD we repeat the evolution for Q QCD = 0.52 GeV and Q QCD = 1 GeV and we compute the relative differences with respect to the chosen value of 0.7 GeV (4.5) In Fig. 12 we show the results for the colored particles, which are the most affected by the choice of the QCD scale, while for the photon and the leptons the relative differences are smaller than 10 −5 . We report the results at Q = m W after the QED+QCD evolution, since this is the phase in which these effects are stronger.

Discretization
The second source of uncertainty we take into account is the discretization, that is the number of grid points N x . Again we focus for simplicity on the first phase of the evolution, since these effects do not change much with the energy scale 6 . As for Q QCD we repeat the evolution for different values of N x and compute the relative differences of the PDFs obtained. The results at the EW scale m W , obtained varying N x from 300 to 600 and then to 1000, are reported in Fig. 13: it is clear that as we increase N x the relative differences are reduced, as expected since in this way we are approaching the continuum limit. Being the relative difference between N x = 600 and N x = 1000 already of O(10 −2 ), we take the latter as reference value, since a further increase will introduce even smaller corrections.

Integration step
Numerical uncertainties also arise due to the discretization in t, depending on the choice of the integration step dt, as shown in Eq. (E.7). As for the previous cases, we compute the relative differences of the PDFs obtained for two different values of dt, in particular we choose dt = t(m W )/N t , with N t = {100; 200}: this means that we consider N t steps in the first phase of the evolution. We do not report any plot, since we checked that the relative differences are at most of O(10 −3 ), both at Q = m W and at higher scales.

Higher orders
The largest theoretical uncertainties in our results originate from neglecting higher order corrections. In particular, in the DL approximation terms of O(α 2 log(Q/Q EW )) are not consistently resummed [60]. For example, at Q = 3 TeV, these terms already amount to 10%, while at Q = 10 TeV to 14%. We notice that promoting our approximation to the full LL result does not improve the situation [60], since single-log terms of the same size from the NLL expansion are still present. However, performing the NLO matching as prescribed in Refs. [48] can reduce the uncertainties to 4% for Q = 3 TeV. Extending our formalism to NLL order would eventually correspond to < 3% accuracy regardless of the energy scale.

PDFs for electron beams
Our numerical code, with obvious substitutions, can also be used to derive LePDFs for electron beams. While most future projects for e + e − colliders are focussed on EW-scale energies to perform high-precision studies of EW gauge bosons, the Higgs, and top quark, linear collider projects also envisage later stages with TeV-scale center of mass energies. In this case our SM PDFs can provide a useful tool. We therefore provide public PDFs for electron/positron beams alongside those for muons and anti-muons. In Fig. 14  only an example plot for some PDFs of an electron. As for the muon case, in case of EW gauge bosons PDFs we did a comparison with the EVA approximation at LO and NLO, obtaining similar results as shown above for the muon. The main difference in the PDFs of an electron compared to those of a muon is that photon, charged leptons and quarks PDFs are larger. This is due to the longer QED evolution from m e to m µ . A consequence of larger quarks PDF is also a larger gluon one, even if Q QCD > m µ . On the other hand, EW gauge boson PDFs are very similar since, at first order, their PDF is insensitive on physics at scales below the EW one.

Conclusions
Obtaining precise predictions for multi-TeV lepton colliders is a topic of active studies. In fact, while QCD radiation plays a minor role, in comparison to hadron colliders, electroweak corrections in this energy regime can become very large due to single and double logarithmic enhancements. Furthermore, the electroweak sector presents several features that are absent when dealing with QED or QCD radiation. Within this context, the resummation of a subset of large logarithms related to the factorisable emission of initial-state radiation can be viewed as an ingredient for a complete description of collider phenomenology at such machines.
In this paper we solved the set of DGLAP equations for an initial-state lepton, evolving the complete set of PDFs from the infrared up to multi-TeV scales. Our computation is performed with LO splitting functions, keeping into account all relevant mass thresholds, EW symmetry breaking terms, masses of all EW states, and resumming EW Sudakov logs at the double-logarithmic level. The residual uncertainty is dominated by the incomplete single-log EW resummation and can be estimated to be of O(10%), while we show that other systematic uncertainties are fully under control. Improving our results to include single-log resummation is left for future work.
Using this result we discuss several notable features of LePDFs. Polarisation effects, due to the chiral nature of SM interactions, are shown to be of O(1) in Fig. 6, or even larger than that in case of the b quark. This last effect is due to the interaction with the top quark via the large top Yukawa coupling, and becomes much smaller in the 5FS, where the top is not included in the evolution. As already observed in previous studies, we confirm that including EW states masses in the propagators gives sizeable correction to the PDFs, even for large factorisation scales.
Furthermore, we perform a detailed comparison of our results for the EW gauge bosons PDFs with the widely used Effective Vector Approximation. This comparison, presented in Fig. 7, illustrates how the EVA fails to describe correctly the transverse EW gauge bosons PDFs at small x values, with deviations reaching even O(50%) for W ± T and Z T at large scales, or even missing the target by more than one order of magnitude in case of the off-diagonal Z/γ PDF. The cause of the latter large deviation is the well-known accidental cancellation in the vector-like coupling of a lepton to the Z boson, which does not take place in the full result since the muon gains a strong polarisation. The smaller, but still substantial deviations in the PDFs of transverse EW gauge bosons are due to the fact that the EVA, treated at LO, does not include contributions from gauge boson splitting off an initial gauge boson. Such terms, while being formally of higher order, become enhanced due to the large γ, W, Z PDFs and the fact that these splittings come with associated Sudakov double logs. In fact, by extending iteratively the EVA up to O(α 2 ) we obtain a much better agreement with the complete numerical result, as shown in Fig. 8. In light of this, we recommend the use of LePDF to derive precision predictions for SM and BSM processes at high-energy electron or muon colliders (they could also be used for lepton-hadron collisions), when one is interested in being inclusive on radiation emitted at small p T compared to the typical energy of the hard scattering, E ≫ p coll.rad.

T
, and when E ≫ m W , which are the conditions for factorisation to be valid.
Our numerical results for the LePDFs are made public 7 in the LHAPDF6 format (with some modifications due to the necessity of describing independently all helicity states). We provide results for initial-state muons, anti-muons, electrons, and positrons, all in both the 5 and 6 flavour schemes.

A.1 Standard Model inputs
In order to set our notation, we define the covariant derivatives as where t A and T a are SU (3) and SU (2) generators, respectively, while Y is the hypercharge. The electric charge is given by Q = Y + T 3 . In our analysis we neglect the CKM matrix and work with diagonal Yukawas where H c = iσ 2 H * . In particular, in our implementation we keep only y 3 u ≡ y t ̸ = 0. Finally, we define the Higgs potential as To shorten the notation, we define the following quantities where N c = 3. We evaluate the RG evolution for the QCD coupling using the 3-loop result [73] with mass thresholds for the charm, bottom, and top quarks. For the QED and EW couplings, as well as for the top Yukawa, we employ the corresponding 1-loop RGE. The numerical boundary conditions are taken at Q = 200 GeV from [72], we report the relevant ones in Table 3 for convenience.

A.2 Formalism for the DGLAP equations
We consider a process A + X → C + Y mediated by the particle B at tree level, with B carrying a fraction z ∈ [0, 1] of the energy of A, as in Fig. 1. The kinematics, up to quadratic order in the transverse momentum p T , is the following where we remind thatz = 1 − z. In this way the emitted particle C is on-shell (p 2 C = m 2 C ), while for the particle B, neglecting O(p 4 T ) terms, we have Since we are working with high-energy initial beams, we can expand p 2 B neglecting terms of order m/E and p T /E or higher, so that the virtuality of the particle B is given by Invariance of the cross section on the factorization scale Q, that we choose to be the p T of the emitted parton, gives the DGLAP equations [43,74] describes the splitting process. The standard matrix elements, which are also present in the unbroken phase (i.e. in a massless theory), are typically parametrized as where P C BA (z) is the splitting function and α ABC the corresponding coupling. Ultra-collinear matrix elements are instead proportional to v 2 and they can be parametrized with new splitting functions U C BA as The general DGLAP equation for a parton B is then given by where P and U are the splitting functions for massive partons and are obtained from those for massless ones with the redefinition in Eq. (3.3)

B Splitting Functions
Here we list all the splitting functions, including ultra-collinear ones. Some of them have ā z pole and therefore they introduce divergences in the DGLAP equations when z max = 1.
To deal with such divergences, we use the + distribution, defined as (B.1) As already discussed, since the SM is a chiral theory we separate vector polarizations and fermion helicities in the splitting functions.

B.1 Massless splitting functions
We start with the splitting functions of the form in Eq. (A.12). Here f labels a fermion, V a gauge boson and h a scalar. We do not specify the polarization of the particle C, since in the computation we sum over it. (B.5) Since QED and QCD are vectorlike theories, we also define the splitting functions properly summed over the vector polarizations and fermion helicities: Finally, we report here the integrals appearing in Eq. (2.4): (B.7)

B.2 Ultra-collinear splitting functions
The top quark is explicitly written, while for other fermions we write f . s = L, R is the helicity of the fermion, T = ± is a transverse polarization of the gauge bosons. If inside a splitting we write f s f −s it means that the two fermions have opposite helicity (same for the gauge bosons). N f c is 1 for leptons and N c for quarks. (B.10) (B.13) Splitting f → f + V L : (B.14) Splitting t → t + h: Splitting h → t +t: (B.20) Splitting h → V L + V L : Splitting involving the mixed state hZ L : (1 +z) , zz(z − z) . (B.28)

C Radiative corrections
To cancel the IR divergences appearing in the splitting functions in the soft limit z → 1 we need to add virtual corrections. Instead of explicitly computing the corresponding Feynman diagrams, we use that virtual corrections correspond to a process in which the particles B and A are the same and C is absent, so in the splitting formalism they can be treated introducing a new splitting function for each particle of the form P v BB (x, t) = P v B (t)δ(1 − x). Inserting the new term in the DGLAP equations we get The coefficients P v B (t) can be computed using momentum conservation where f (2) i represents the n = 2 Mellin transform of the PDF Deriving in t and using Eqs. (A.10, C.1) we obtain for each particle A From the definitions in Eqs. (A.11, A.12, 3.3, A.13), we get In order to reproduce the non cancellation of IR divergences in SU (2) L interactions, the computation of virtual corrections is modified by changing the boundary of the integral as discussed in Section 3.2 We show explicitly the virtual corrections for the QED+QCD phase, with N ℓ charged leptons, N u up and N d down quarks. All the particles are treated as massless, so momentum conservation equations become just relations between the n = 2 Mellin transform of the splitting functions and the virtual coefficients. Applying Eq. (C.6) respectively to the fermions, the photon and the gluon we get, using Eq. (2.2) where N QED is the effective number of fermions, N q = N u + N d is the number of quarks and δ f,q = 1 for quarks and 0 for leptons. In our specific setup, N ℓ = N d = 3 and N u = 2, so N q = 5 and N QED f = 20 3 . We omit writing the explicit expressions for virtual corrections in the full SM phase, as they are many and too lengthy.

D DGLAP evolution equations in SM
Here we list the full set of DGLAP equations we used above the EW scale. We use the following notation for transverse vector polarization: In the ultra-collinear terms the factor v 2 /(16π 2 Q(t) 2 ) appearing in front of each splitting function is omitted to shorten the notation.

E Numerical implementation
Here we show the details of our numerical implementation of the DGLAP equations. We discretize the x-space from a minimum value x 0 up to 1 in N x bins, x α = {x 0 , x 1 , . . . , x Nx ≡ 1}, with spacing δx α = x α − x α−1 . We choose a spacing that is denser near x = 1 and sparser at small values, in practice we set x α = 10 −6((Nx−α)/Nx) 2.5 for α = 0, 1, . . . , N x to get values from x 0 = 10 −6 to x Nx = 1. This allows us to obtain a set of ODEs, where the integrals are computed using the rectangles method. 8 For the non-divergent terms or when z ABC max ̸ = 1 we obtain where N ABC max is the index of the greatest point of the grid which is smaller than z ABC max . The remaining case is that of the splitting functions with the + distribution, which we can write as With our discretization, the last integral becomes where X β is given by (E.5) With this discretization, starting with N f PDFs we get a set of (N x + 1)N f equations for the variables f Bβ (t) ≡ f B (x β , t), with β = {0, . . . , N x } and B = {1, . . . , x β x α f Aα (t) .

(E.6)
Once we have the equations with the proper initial conditions, we solve them numerically using a fourth order Runge-Kutta with integration step dt. For a set of equations of the form y ′ i (t) = F i (t, y), starting with the solution y i,n at t = t n the step of the algorithm is k 1,i = dtF i (t n , y i,n ) k 2,i = dtF i (t n + dt 2 , y i,n + k 1 2 ) k 3,i = dtF i (t n + dt 2 , y i,n + k 2 2 ) k 4,i = dtF i (t n + dt, y i,n + k 3 ) .

(E.7)
We then reduce the number of variables imposing momentum conservation after every step of the evolution, since performing a numerical integration it will be more and more violated, as done in [34]. With our discretization Eq. (C.2) becomes δx α x α f iα (t) = 1 .

(E.8)
Since the initial conditions on PDFs are given by only f µα will be nonzero for α = N x throughout the evolution. Then we fix 10) reducing by N f the number of variables and solving the remaining equations. The factor L(t) is computed using momentum conservation δx α x α f iα (t) + L(t) , (E.11) that is δx α x α f iα (t) . (E.12) The uncertainties due to the discretizations in x and t spaces are discussed in Section 4.6.

F LHAPDF format
We publish our numerical results in a format inspired by the LHAPDF6 [75] one used for proton PDFs. Some changes are required due to the polarisation of PDFs. The structure of the output is then as follows: • the first three lines just specify the format; • in the fourth and fifth line are reported respectively the grids in x and in Q (the grid in Q is a subset of the grid used in the numerical solution of the DGLAP equations); • the next three lines report the particles' list as in Table 4: name, PDG ID (for the Z/γ and h/Z L interference we join the PDG ID of the two states) and the additional label specifying the helicity (it is understood that for fermions or vectors it will be ±1/2 or ±1, respectively); • in all the remaining lines we report the quantities xf (x, Q): each column corresponds to a particle, following the order of the previous lines. We start at x = x 0 increasing Q at each row and repeating for each x, so that the data have the form reported in Table 5.   Table 4. Names, PDG ID and polarisations of the particles. In particular, the second, third and fourth columns of the tables correspond to the sixth, seventh and eighth lines of the output file.  Table 5. Structure of our PDF data in the LHAPDF6 format.