The partonic structure of the electron at the next-to-leading logarithmic accuracy in QED

By working in QED, we obtain the electron, positron, and photon Parton Distribution Functions (PDFs) of the unpolarised electron at the next-to-leading logarithmic accuracy. The PDFs account for all of the universal effects of initial-state collinear origin, and are key ingredients in the calculations of cross sections in the so-called structure-function approach. We present both numerical and analytical results, and show that they agree extremely well with each other. The analytical predictions are defined by means of an additive formula that matches a large-$z$ solution that includes all orders in the QED coupling constant $\alpha$, with a small- and intermediate-$z$ solution that includes terms up to ${\cal O}(\alpha^3)$.


Introduction
The typical cross section relevant to e + e − collisions is in principle entirely computable as a perturbative series in the QED coupling constant α. In practice, however, this is hardly useful, since the coefficients of such a series are very large, thus compensating the suppression due to α -in other words, all terms of the series might be of the same order numerically, which leads to a complete loss of predictive power. The problem stems from the fact that the incoming e ± particles tend to copiously radiate photons (which in turn may convert into e + e − pairs, and so forth) at small angles w.r.t. the beamline. In perturbation theory, any zero-angle emission would induce a divergent cross section, were it not for the screening effect provided by the mass of the emitter and/or the emitted particle. Thus, when integrating over all possible emissions, the cross section will contain logarithms of the ratio m 2 /E 2 , where E is a scale of the order of the hardness of the process, and m is the screening mass (i.e. that of the electron in the case we are interested in). It is these logarithms that, by growing large when m 2 /E 2 1, give the dominant contributions to the perturbative coefficients, and ultimately prevent the series from being well behaved.
Fortunately, such log(m 2 /E 2 ) terms are universal, and because of this they can be taken into account to all orders in α by a process-independent resummation procedure. In the so-called structure-function approach, the physical cross section is then written by means of a factorisation formula such as the following one 1 : dσ e + e − (p e + , p e − , m 2 ) = ij=e ± ,γ dz + dz − Γ i/e + (z + , µ 2 , m 2 ) Γ j/e − (z − , µ 2 , m 2 ) × dσ ij (z + p e + , z − p e − , µ 2 ) . (1.1) The quantities Γ i/e ± are called the Parton Distribution Functions (PDFs) of the electron or the positron, a name that originates from the analogy of eq. (1.1) with its QCD counterpart. PDFs collect and resum all of the log(m 2 /E 2 ) terms; conversely, the short-distance cross sections dσ ij do not contain any such logarithms, and are expected to be well-behaved order by order in perturbation theory. Neither Γ i/e ± nor dσ ij are physical quantities; their definitions always involve some degree of arbitrariness, which is parametrised by the mass scale µ, that is only constrained by the requirement µ ∼ E, and by the chosen factorisation scheme. Fuller details on the usage of the factorisation formula (1.1) in calculations relevant to e + e − colliders and on its physical meaning can be found e.g. in ref. [1]. In particular, we shall adopt the notation of ref. [1], whereby the incoming e ± are called particles (with dσ e + e − thus being a particle-level cross section, defined so as to retain only terms that do not vanish in the m/E → 0 limit), and the objects i and j are called partons (so that dσ ij is a parton-level cross section). This allows one to distinguish easily between an electron that stems from one of the collider beams, and an electron that initiates the hard collisions, and that stems from the PDF Γ e − /e − .
We point out that it is somehow customary in QED to call Γ e − /e − (Γ e + /e + ) the electron (positron) structure function. This is motivated by the fact that, by ignoring the contributions of partons whose species is not the same as that of the incoming particle, and by working at the first order in perturbation theory, a structure function (which is an observable) can be made to coincide with the PDF relevant to the case where particle and parton have the same identity, by means of a suitable definition of such a PDF. This position is not tenable at higher perturbative orders, and when more parton species are allowed in any given particle. Therefore, "structure functions" will not be used in this paper, and we shall only refer to PDFs.
The crucial point that we need to stress here is that in QED e ± PDFs, at variance with hadronic PDFs, are entirely calculable with perturbative techniques. Presently they are known in close analytical forms [2][3][4] which are leading-logarithmic (LL) accurate, and that include all-order in α contributions in the region z ± 1 (which is responsible for the bulk of the cross section), matched with up to O(α 3 ) terms for any values of z ± ; both of these forms exploit leading-order (LO) initial conditions. The goal of the present work is to improve on the results of refs. [2][3][4] by extending them to the next-to-leading logarithm accuracy (NLL) starting from the next-to-leading order (NLO) initial conditions computed in ref. [1]. In keeping with what was done in the literature, we shall present predictions both for all-order PDFs in the z ± 1 region, and for up to O(α 3 ) NLL terms valid for any z ± . By working at the NLL+NLO accuracy, the mixing between the electron/positron and the photon PDFs is taken into proper account, as are running-α effects. Our results are obtained with both analytical and numerical methods, which are compared and used to validate each other. This paper is organised as follows. For those readers who are not interested in the technical procedures which underpin this work but only in their final outcomes, we summarise our final results in sect. 2. The details of the derivations of such results are then given in the remainder of the paper. In sect. 3 we introduce the evolution equations for the PDFs that we are going to solve, and report the associated initial conditions. Sect. 4 briefly describes the evolution-operator formalism. Analytical solutions are computed in sect. 5, for any z ± values in sect. 5.1 (the resulting lengthy expressions are partly collected in appendix A), and for z ± 1 in sect. 5.2 (with additional details reported in appendix B), while a description of the codes employed to obtain numerical results is given in sect. 6. The solutions of sects. 5.1 and 5.2 are combined ("matched") in sect. 7. Our analytical and numerical predictions are extensively compared in sect. 8. Finally, we conclude and give a short outlook in sect. 9. Additional material is collected in the appendices.

Synopsis of results
The e ± PDFs that we shall compute in this paper result from solving the evolution equations of eqs. (3.7)- (3.8), with the initial conditions given in eqs. (3.18)- (3.21). The all-order, large-z solutions and one of the numerical codes we shall employ require the use of an evolution operator, whose RGE is presented in eq. (4.15). The latter can be solved in a closed form in the case of a one-dimensional flavour space; such a closed form is reported either in eq. (4.17) or in eq. (4.21), the two differing by terms of O(α 3 ) and with the latter form suitable to take the fixed-α limit. The two-dimensional flavour space case is discussed in appendix B. The PDF solutions valid for any z and including up to O(α 3 ) contributions are represented in terms of sets of basis functions, which are obtained by solving recursive equations; these are given in eqs. (5.30) and (5.31), and their explicit solutions partly in appendix A, and partly in an ancillary file (see below). Conversely, the all-order, large-z solutions for the PDFs are reported in eq. (5.46) (LL accurate for singlet and non-singlet), eq. (B.88) (LL accurate for photon), eq. (5.63) (NLL accurate with running α for singlet and non-singlet), eq. (B.74) (NLL accurate with running α for photon), and eq. (5.68) (NLL accurate with fixed α for singlet and non-singlet; the corresponding result for the photon is obtained from eq. (B.74), but is not reported explicitly). The all-z and large-z solutions are then matched in an additive manner, as is shown in eq. (7.4). The arXiv submission of the present work will be accompanied by two ancillary files, that will contain the main results for the PDFs as Mathematica formulae, and some analytical results too long to fit in this paper. Furthermore, a numerical code that returns the PDFs will be made public, to be downloaded at: https://github.com/gstagnit/ePDF.

Evolution equations and initial conditions
By working in QED the cases of the electron and of the positron PDFs are identical. Thus, in order to be definite in this paper we shall only consider the PDFs of the electron, which allows us to simplify the notation of ref. [1] in the following way: Γ i (z, µ 2 ) this paper ≡ Γ i/e − (z, µ 2 ) ref. [1] . (3.1) The evolution equations are therefore [5][6][7][8]: Henceforth, we shall omit to write the explicit z and/or µ dependences when no confusion can possibly arise. By working with a single fermion family, eq. (3.2) becomes:

14)
P NS = P e ± e ± − P e ± e ∓ ≡ P V ee − P V eē . (3.15) After solving the evolution equations for the singlet and non-singlet components, one recovers the solutions for the electron and the positron by inverting eq. (3.6): The electron PDFs can be expanded perturbatively. We denote the first two coefficients of such an expansion in the same way as in ref. [1], namely: The evolution equations are supplemented by the initial conditions computed up to O(α) in ref. [1]. These read as follows: where µ 0 m, and m is the electron mass. The rightmost terms on the r.h.s. of eqs. (3.19) and (3.20) are associated with, and fully determined by, the scheme used to subtract the initial-state collinear singularities. In this paper, we work in the MS scheme, which implies: We conclude this section with a general remark on evolution. Throughout this paper, by "evolution" we understand the one governed by RGE's. This implies, in particular, that the contributions of e + e − low-energy data to the running of α are locally neglected (i.e. for scales of the order of the masses of light hadronic resonances). However, nothing prevents one from taking into account such contributions in an inclusive way. Namely, by starting from a precise determination of α = α H that does include low-energy contributions, and that can be associated with a scale µ H (just) larger than the mass of the heaviest hadronic resonance, one can backward-evolve α H = α(µ H ) from µ H down to µ 0 , thus determining the value of α(µ 0 ) that is employed in this work. By doing so, the possible local effects of the resonances on the evolution of the PDFs are still neglected, but this is not important: in the factorisation formulae such as eq. (1.1) where the PDFs are used, the scales are meant to be hard and therefore never assume values comparable to the masses of the light hadronic resonances.

Evolution operator
As far as the evolution in µ is concerned, eqs. (3.7) and (3.8) are identical. We shall thus start dealing with the former one, which has a more involved flavour structure; the results will then be applied to the non-singlet case as well, by simply considering a onedimensional flavour space. We re-write eq. (3.7) by means of a simpler notation, where all of the irrelevant indices are dropped: and Γ is a column vector. We define the Mellin transform of any function f (z) whose domain is [0, 1] as follows: If f (z) is the convolution of two functions g(z) and h(z): then: The evolution kernels are expanded in a perturbative series whose coefficients follow the same conventions as those in eq. (3.17), namely: By denoting by Γ 0,N ≡ Γ N (µ 2 0 ) the PDF initial conditions at the reference scale µ 0 , and by introducing the evolution operator E N (µ 2 , µ 2 0 ) such that: eq. (4.6) becomes: Since eq. (4.8) must be true regardless of the specific choice for Γ 0,N , it is equivalent to: Following ref. [9], it is appropriate to introduce the variable 3 : We use the following definition of the QED β function: 12) and n F the number of active charged fermion families. Equation (4.10) implies that: and thus: With eq. (4.13), eq. (4.9) becomes 4 : Note that, from eq. (4.7), E N (t = 0) = I. If the flavour space is one-dimensional (as for the non-singlet evolution), eq. (4.15) can be solved analytically. Notation-wise, we deal with this case by performing the formal replacements: By exploiting eq. (4.14), one readily obtains: By construction, the O(α 3 ) terms neglected in eq. (4.17) stem from the truncation of the series that gives the evolution kernels in eq. (4.5); conversely, the relationship between α(µ) and α(µ 0 ) is treated exactly thanks to the usage of the variable t. If one wants to expose explicitly the large logarithms that originate from having µ µ 0 , one can use the following series expansions: having defined: (4.20) By employing these results, eq. (4.17) becomes: This result is useful because, at variance with that of eq. (4.17), it allows one to consider the case of a non-running α, which can simply be obtained from eq. (4.21) in the limit b 0 → 0. As a consistency check, it is immediate to verify that, by taking such a limit, one arrives at a form for log E N which could have been directly obtained from eq. (4.9), by working in a one-dimensional flavour space and by freezing α(µ) there.

Analytical solutions for the PDFs
In this section we obtain the NLL-accurate PDFs of the electron in closed analytical forms in two different ways: by solving the evolution equations order by order in perturbation theory (sect. 5.1), and by using the properties of the evolution operator to obtain the asymptotic behaviour in the z → 1 region to all orders in α (sect. 5.2). These two results can then be combined in order to obtain predictions which are numerically well-behaved in the whole of the z range (sect. 7).

Recursive solutions
Following ref. [2], perturbative solutions for the evolution equations can conveniently be obtained by re-writing eq. (4.1) in an integral form: with 5 : 2) and the modified convolution operator defined as follows: which is a valid definition regardless of whether g(x) is a distribution or an ordinary function. Note that F is a column vector, and that eq. (5.1) has a matrix structure, in the flavour space. As was the case for eq. (4.1), this implies that all of the results to be obtained in the following can be applied to the limiting situation of a one-dimensional flavour space as well.
The procedure of ref. [2] is LL-accurate. In order to generalise it to the NLL accuracy we are interested in in this work, it is best to first consider the case of non-running α. With this assumption, the solution of eq. (5.1) can formally be written as follows: From this equation, F can be obtained by representing it by means of a power series: where: and with I LL k and I NLL k two sets of unknown functions 6 . By replacing eq. (5.5) into eq. (5.4), the two sides of the latter equation become two series in η 0 : one then equates the coefficients relevant to the same power of η 0 on the two sides, thereby obtaining equations 5 The use of a Θ function in eq. (5.2) guarantees its validity also when Γ is a distribution, and thus allows one to take into account its possible endpoint contributions. Conversely, while F should also be treated as a distribution, we shall regard it as an ordinary function, because in the large-z region we shall in any case employ the asymptotic solutions whose results, given in sect. 5.2, are more accurate there. 6 More precisely, eq. (5.5) implies that for any k, I LL k and I NLL k are two-dimensional column vectors in the singlet-photon flavour space, whose elements are functions of z, and c-number functions in the non-singlet flavour space. An extra flavour index will be included in the notation when distinguishing the flavour components will be important (see appendix A). that can be solved for I LL k and I NLL k (recursively in k). The r.h.s. of eq. (5.5) is simply an expansion in terms of αL, and thus η 0 is a convenient expansion parameter, irrespective of the logarithmic accuracy one is working at. Indeed, eq. (5.5) can be extended by adding further contributions to its r.h.s., that are suppressed by higher powers of α. Conversely, by keeping only the I LL k contributions, one recovers what was done in ref. [2]. The recursive solutions for I LL k and I NLL k stemming from eq. (5.5) read as follows: with: The quantities in eq. (5.9) must be obtained by direct computation by using the definition in eq. (5.2), with the perturbative expansion of eq. (3.17) and the initial conditions of eqs. (3.18)-(3.21). By doing so, we obtain: The key to the simplicity of the solutions in eqs. (5.7) and (5.8) is the fact that the dependence on µ on the r.h.s. of eq. (5.5) is entirely parametrised by L, which in turn allows one to compute the integral on the r.h.s. of eq. (5.4) in a trivial manner: Unfortunately, things are not so simple when α is running. In this case, as was already done in sect. 4, it is convenient to use the variable t introduced in eq. (4.10). Owing to eq. (4.13), the analogue of eq. (5.4) reads as follows: As a consequence of this, we shall use the representation: 16) rather than that of eq. (5.5). Thus: The r.h.s. of eq. (5.15) therefore features two independent classes of integrals, namely: In order to evaluate eq. (5.19), we make repeated use of eq. (4.14). Then: By direct computation: with: We have thus: where: The r.h.s. of eq. (5.23) can be simplified by means of algebraic manipulations of the summation indices: and: The initial conditions must then be written as follows: By replacing the results of eqs. (5.26), (5.27), and (5.29) into eq. (5.23), and by using the representation of eq. (5.16) for F(z, t), we obtain the sought recursion relations: with: These results generalise those obtained in the case of non-running α, which can be obtained from them. Indeed, in the limit of fixed α, which at the NLL can be achieved by letting b 0 → 0 and b 1 → 0 (with b 1 /b 0 → 0), we have t → η 0 /2, whence eq. (5.16) coincides with eq. (5.5), if one identifies J LL with I LL and J NLL with I NLL . This is justified, since eqs. (5.7) and (5.30) are identical, and the recursive relation of eq. (5.31) coincides with that of eq. (5.8) when α is not running. After solving eqs. (5.7) and (5.8) for I LL k and I NLL k , with the definition in eq. (5.2) one arrives at the following representation of the PDF in the case of fixed α: where Analogously, in the case of running α: We point out that with, for example, α(µ) = 1/128 and α(µ 0 ) = 1/137 we have t 0.1/n F . Furthermore, since: the numerical coefficients in front of the convolution products and of the initial conditions in eq. (5.31) are of order one. Therefore, the series of eq. (5.16) is expected to be poorly convergent only for z → 1 and z → 0, owing to the possible presence of log p (1 − z) and log p z terms in the J LL and J NLL functions. Part of the explicit results for the functions J LL k (with 0 ≤ k ≤ 3) and J NLL k (with 0 ≤ k ≤ 2) are reported in appendix A, and part in one of the ancillary files.

Asymptotic large-z solutions
The electron PDF is equal to δ(1 − z) at the LO (see eq. (3.18)); while the LL evolution of such an initial condition does smooth its behaviour, resulting in a tail that extends down to z = 0 [2][3][4], the PDF remains very peaked towards z = 1, where it has an integrable singularity. This implies that the perturbative expansion of the LL-accurate solution features log(1 − z) terms at each order: if one truncates such a perturbative series, one exposes a non-integrable divergence at z = 1, regardless of the order at which the truncation occurs. The same is true when NLO initial conditions and NLL-accurate evolution are considered, as is explicitly shown by the results in appendix C.
In order to address this issue, the log(1 − z) terms must be resummed. This can conveniently be done by exploiting the evolution-operator formalism presented in sect. 4, whose usage is simplified by the observation that the large-z region corresponds to the large-N region in Mellin space: Thus in this section, when dealing with Mellin transforms and their inverse, we shall often implicitly assume eq. (5.38). A second crucial observation is that the z → 1 asymptotic behaviours of the singlet and non-singlet components are actually identical. This implies that also in the former case one can effectively work as if the evolution operator were a c-number and not a matrix, thereby allowing one to exploit the closed-form solutions of eqs. (4.17) and (4.21). Therefore, we shall start by understanding the non-singlet notation in sects. 5.2.1 and 5.2.2 in order to derive the main results relevant to the asymptotic z → 1 region; we shall return to and comment on the singlet-photon case in sect. 5.2.3 and in appendix B.

LL solution for non-singlet
Given that the LL-accurate result has been available for a long while [6], this case is presented here only to show how the evolution-operator formalism helps find the asymptotic solution in a straightforward manner. At the LL we are entitled to neglect the running 7 of α. Thus, the appropriate form for the evolution operator is obtained by keeping only the O(α) term of eq. (4.21), with α(µ) → α there, supplemented by the LO initial condition: From eqs. (4.7) and (5.39) we obtain: A direct calculation in the large-N region leads to: where all terms suppressed by at least one inverse power of N have been neglected, and we have defined: We point out thatN is a quantity that routinely appears in the computation of Mellin transforms, and which helps retain some universal subleading terms. Therefore: with η 0 defined in eq. (5.6). Equation (5.43), when substituted into eq. (5.40), implies: The inverse Mellin transform can now be evaluted by using the following result, valid for any κ > 0: The comparison of eq. (5.45) with eq. (5.44) allows one to arrive at the final result [6]: This is identical to what is nowadays a rather standard form, except for an exponentiated term of pure-soft origin (stemming from the use of β exp = β, rather than of β exp = η, as defined e.g. in eq. (67) of ref. [10]). Such a term clearly cannot be obtained by means of the collinear resummation carried out here.

MS NLL solution for non-singlet
At the NLL, the PDF initial conditions must be set as given in eqs. (3.18) and (3.19), with K ee = 0 in the latter equation (see eq. (3.22)). By exploiting the property of the Mellin transform of eq. (4.4), we have: with log E N given in eq. (4.17) (where running-α effects are also included). With eq. (5.41) and its NLO analogue: where: we can cast the logarithm of the evolution operator in the same form as in eq. (5.43), namely: having defined: Equation (5.50) implies that we can follow the same steps that have led us to eq. (5.46), and therefore: We must now replace this result into eq. (5.47). In this way, two independent convolution integrals emerge: A tedious but otherwise relatively straightforward procedure leads to the following results: where, inside the square brackets, we have neglected terms that vanish at z → 1. We have introduced the two functions: where: By putting everything back together, we arrive at the final result: which is therefore the NLL-accurate counterpart of eq. (5.46). A couple of observations about eq. (5.63) are in order. Firstly, owing to eqs. (5.51) and (4.19), we have ξ 1 η 0 . With µ 0 and µ of the order of the electron mass and of a few hundred GeV's, respectively, one obtains η 0 ∼ 0.05. Therefore, both the LL and the NLL solutions are still very peaked towards z = 1, where they diverge with an integrable singularity. Furthermore, the z → 1 behavior of eq. (5.63) is worse than that of eq. (5.46) because of the presence of the explicit log p (1 − z) terms in the former equation 8 . Secondly, the small numerical value of ξ 1 just mentioned implies that the following expansions: are rather accurate approximations of the complete results of eqs. (5.60) and (5.61). Equation (5.64), in particular, implies that numerically the log(1 − z) term is much larger than the (formally dominant) log 2 (1 − z) one, even for z values that are extremely close to one. This fact might be significant when performing the integral of the convolution between electron PDFs and short-distance cross sections.
From eq. (5.63) one can also readily obtain a LL-accurate solution, where at variance with that of eq. (5.46) the effects due to the running of α are included. Explicitly: this is again consistent with the findings of ref. [6]. Finally, the running of α can formally be switched off in the NLL-accurate solution. In order to do so, one must repeat the procedure that leads to eq. (5.63); however, rather than using the expression of the evolution operator given in eq. (4.17), one must use that of eq. (4.21). By doing so, one arrives at: where

The singlet and photon cases
The key result relevant to the evolution of the singlet and photon sector in the z → 1 region is the following: that is obtained by means of a direct computation starting from the definitions given in sect. 3. Equation (5.71) implies that the singlet and the photon evolve independently from each other in this limit. Since the kernel evolution is a diagonal matrix, so is the evolution operator, and therefore the solutions for its elements on the diagonal are given by either eq. As far as the photon is concerned, eq. (4.12) and the photon-photon elements in eq. (5.71) imply that the second term on the r.h.s. of eq. (4.17) is equal to zero. Therefore: having used eq. (4.10). The convolution with the initial conditions of eqs. (3.18) and (3.20) is thus trivial, and the final result reads as follows: We can observe the presence of an α(µ) term in the denominator of eq. (5.73) which, typically, will cancel an analogous factor in the short-distance cross sections (given that these, for consistency with the present results, will have to be computed in the MS scheme). Thus, one sees the natural emergence of quantities that are employed in the so-called α(0) scheme (see e.g. ref. [12]) -this is the same mechanism that has been anticipated in ref. [13] in the case of photon fragmentation functions. We stress that this properties stems from the computation of both the PDFs and the cross sections in the same collinear subtraction scheme, and from the solution of the evolution equations for the PDFs. Unfortunately, eq. (5.73) does not give a good description of the true large-z behaviour of the photon PDF. This is because the off-diagonal terms of the evolution kernel imply that such a PDF receives a contribution that primarily stems from the initial conditions of the electron PDF. As we have seen previously, these are much more peaked towards z = 1 than their photon counterparts, so much so that this behaviour compensates the fact that the off-diagonal elements of the evolution kernel are suppressed w.r.t. the diagonal ones, which are the only ones that have been taken into account in eq. (5.71). It then follows that, in order to improve on the solution in eq. (5.73), one needs to solve the evolution equations of the singlet-photon system by including those off-diagonal elements. In turn, this entails a significant increase in complexity, and for this reason we refrain from discussing the relevant procedure here -all of the results are reported in appendix B.

Numerical solutions for the PDFs
The numerical evolution for the PDFs is achieved by first solving the evolution equation for the evolution operator in Mellin space. More specifically, we solve the equation given in the first line of eq. (4.15) without expanding the β function in the denominator. As is discussed in sect. 3, the introduction of the singlet and non-singlet combinations, eq. (3.6), allows for a decoupling of the evolution equations that is well-suited for a numerical implementation. As has been done thus far, in the following we shall implicitly refer to the two-dimensional singlet-photon case, keeping in mind that the non-singlet case is obtained by considering a one-dimensional flavour space.
The numerical solution of eq. (4.15) for the evolution operator E N is obtained by means of a discretised path-ordered product [14]. The evolution range [0, t] is partitioned into n intervals [t i , t i+1 ], with t 0 = 0 and t n = t, and the evolution operator is written as follows 9 : to it can be evaluated by using the trapezoidal approximation: where we have used eq. (4.14). We have found that a number of intervals n = 20 is appropriate for giving stable results for t values relevant to applications to TeV-range colliders. At the perturbative order at which we are working, the sum over k on the r.h.s. of eq. (6.2) is restricted to k ≤ 1. After having obtained the evolution operator, the PDFs at the hard scale µ are computed in the Mellin space according to eq. (4.7). Finally, in order to invert the PDFs from the Mellin to the z space, we employ a numerical algorithm based on the so-called Talbot path. Details on the implementation of this method can be found in ref. [15]. The computer program that implements what has been described thus far was used to obtain all of the numerical results presented in sect. 8 and appendix A.
In the non-singlet case one can devise an alternative procedure. Namely, one can exploit the analytical N -space solution for the evolution operator, given in eq. (4.17), multiply it by the Mellin-transformed initial conditions, and then invert the result thus obtained back to the z space by means of a numerical contour integration. We have implemented this strategy in a computer program 10 fully independent from the one described above, and verified that the two are in perfect agreement.

Matching
The best analytical prediction is obtained by matching the recursive solution of eq. (5.35), that is valid for all z values but in practice can be computed only up to a certain O(α n ) (here, n = 3), with the solutions of eqs. (5.63) (for singlet and non-singlet) and (B.74) (for photon), that retain all orders in α but are sensible only when z 1. In order to distinguish these two classes of solutions, we now denote them as follows 11 : We remind the reader that eq. (5.35) implicitly encompasses the non-singlet, singlet, and photon cases by means of the J LL k and J NLL k functions (see sect. 5.1). We define the matched PDFs with the additive formula 12 :

4)
10 This builds upon the code originally written by the authors of ref. [16]. 11 We shall henceforth consider the case of NLL solutions with running α, which constitutes our most accurate scenario. However, the procedure is unchanged in the case of NLL solutions with fixed α, or in the case of LL solutions. 12 Additive matching has been considered in refs. [2,17,18]; refs. [3,4] adopt a multiplicative one.
where G(z) is a largely arbitrary function that must obey the following condition lim z→1 G(z) = 1 , (7.5) and that might optionally be used to power-suppress at small z the difference in round brackets in eq. (7.4) -we shall give more details on this point later. The quantity Γ subt (z) (that we call "subtraction term") is responsible for removing the double counting, i.e. the contributions which are present both in the recursive and in the asymptotic solutions. We shall eventually construct it explicitly, but we anticipate the obvious fact that it must feature the dominant z → 1 contributions to the PDF (which, in turn, are present in both the recursive and the asymptotic solutions, as is discussed in appendix C). Before proceeding we stress that, although general, the arguments that follow from eq. (7.4) are best understood if the PDFs are strongly peaked at z → 1, which is indeed what happens for the singlet and non-singlet components, but not for the photon (at least to a certain extent). Thus, we shall first understand the two former cases, and return to the latter one only towards the end of this section.
We want the matched PDF to coincide with the asymptotic or the recursive solution for those z values appropriate for either of the latter two quantities. This is equivalent to requiring: Given eq. (7.5), eq. (7.6) is satisfied when: Conversely, there are two ways in which the behaviour in eq. (7.7) can be achieved.
a) G(z) can be expanded in series around z = 0, and is such that: in addition to satisfying eq. (7.5).
b) One has: for small and intermediate z values. When eq. (7.10) holds, one can set: The option of item a) stems from the observation that both Γ asy (z) and Γ subt (z) are only sensible when the log p (1 − z) terms they contain are large. When this is not the case, i.e. at small-and intermediate-z values, one can suppress them (in fact, one must, if eq. (7.7) is to be fulfilled) by means of power-suppressed terms, here parametrised by G(z), without affecting the formal accuracy of the matched PDF. However, this has the potential drawback of introducing in Γ mtc (z) a dependence on the arbitrary quantity G(z), which must remain small in order not to lose predictive power. This issue is avoided if the option of item b) is viable. This has the drawback that it relies on the condition in eq. (7.10), that might be problematic since it constrains Γ subt (z) in a z region which is not the natural domain of such a function.
Although there is significant freedom in the construction of the subtraction term, the recursive and asymptotic solutions provide us with two obvious candidates. Namely, we can set either In other words: with eq. (7.12) we use all of the contributions to the recursive solution which are non-vanishing when z → 1, while with eq. (7.13) we employ the O(α 3 ) expansion of the asymptotic solution. Therefore, as it follows from the discussion in appendix C, Γ A subt (z) essentially contains a subset of the terms present in Γ R subt (z). More precisely: By construction (see eq. (A.7)), we have and therefore eq. (7.8) automatically holds when the subtraction term is defined by means of the recursive solution. Conversely, However, in spite of this, eq. (7.8) holds also in this case, since: The conclusion is that eq. (7.8) is satisfied with both choices of the subtraction term. The difference between adopting Γ R subt (z) or Γ A subt (z) is that with the former function the matched PDF will converge towards the asymptotic solution at z values relatively smaller than those relevant to the latter function. This can be seen in fig. 1, where the ratio of the l.h.s. over the r.h.s. of eq. (7.8) (without the absolute values) is plotted as a function of − log 10 (1 − z) for the two choices of the subtraction term considered here, and for three different hard scales µ. Note that the scale on the y axis of the plots in fig. 1 is in units of 10  to those shown here. Apart from being in keeping with the expectations emerging from eqs. (7.16)-(7.18), fig. 1 shows that, even in the case of eq. (7.12), the matched PDF will attain its asymptotic form only for values of z which are extremely close to one; in other words, non-logarithmic contributions are almost always important. This being the case, by choosing Γ A subt (z) as a subtraction term rather than Γ R subt (z) (which, as was anticipated, "delays" the onset of the asymptotic regime in the matched PDF) one renders the transition between the asymptotic and recursive solutions less abrupt; this turns out to be beneficial in order to reproduce the results of the numerical evolution.
As far as the small-and intermediate-z region is concerned, we observe that: Equation (7.19) implies that it is unlikely that, if the subtraction term is defined by means of the recursive solution, one can avoid the use of the G(z) function. Conversely, eq. (7.20) implies that the definition by means of the asymptotic solution has a better chance of satisfying eq. (7.10), thus bypassing the need to introduce G(z). Note that the difference in eq. (7.20) is of O(α 4 ) as a direct consequence of the fact that we have computed Γ A subt (z) to O(α 3 ) (see eq. (C.6)). In fig. 2 we plot the ratio of the l.h.s. over the r.h.s. of eq. (7.10) (without the absolute values), by using the same layout as in fig. 1. In order to be definite, we have considered again the non-singlet component in fig. 2, and have verified that the singlet one gives results which are extremely similar to those of the non-singlet. Figure 2 confirms our expectation based on eqs. (7.19) and (7.20).
We now turn to discussing the case of the photon PDF, which is quite different from that of the singlet and the non-singlet. The starting point is the same as for the latter PDFs, namely the definition of the subtraction term with either eq. (7.12) or (7.13), since  those formulae are the general parametrisations of the perturbative expansion of the recursive or the asymptotic solutions, respectively, whose actual values are determined by the parameters (given in appendices A.1, A.2, C.1, and C.2) specific to the particle which is being considered. Indeed, in the case of the photon, the analogues of eqs. (7.14) and (7.15) are: Actually, because of eqs. (C.25) and (C.26), one can make a stronger statement, namely: We stress that eq. (7.23) is not a property inherent to the photon PDF, but a consequence of having been able to keep the relevant subleading terms in the computation of its large-z form as carried out in appendix B. In order to be definite, for consistency with the case of the singlet/non-singlet we shall label the subtraction term with "A" here. The analogues of the right-hand side panels of fig. 1 and of fig. 2 are presented in the right and left panels of fig. 3, respectively. We start from the right-hand side one, in order to assess the validity of eq. (7.8). Unfortunately, it turns out that at large z's the NLL photon PDF becomes negative in a certain range, and it crosses twice the zero. For this reason, we cannot consider the ratio of the two sides of eq. (7.8) as was done in fig. 1, but only plot separately Γ rec (z) − Γ subt (z) and Γ asy (z); these two quantities are displayed in fig. 3 by adopting identical patterns (each associated with a different hard scale µ), with the curves relevant to Γ asy (z) overlaid by full circles. Furthermore, in order for the latter curves to fit into the layout, they have been multiplied by a constant factor equal to 10 −3 . The plot clearly shows how eq. (7.8) is safely fulfilled 13 . Photon at NLL, µ 0 = m e Γ rec -Γ A subt , µ = 0.01 GeV Γ rec -Γ A subt , µ = 1 GeV Γ rec -Γ A subt , µ = 100 GeV Γ asy /1000, µ = 0.01 GeV Γ asy /1000, µ = 1 GeV Γ asy /1000, µ = 100 GeV We now consider the left panel of fig. 3, in order to assess the validity of eq. (7.10); given that for all of the z values employed in the plot the photon PDFs is positive, we can compute the ratio of the two sides of eq. (7.10) (without the absolute values) as was done previously in fig. 2. It is immediate to see that the conclusions are the opposite of those valid in the cases of the singlet and non-singlet -namely, in a very large range in z the subtraction term and the asymptotic solution do not agree with each other 14 . Thus, in the case of the photon the use of a damping function G(z) is unavoidable. In order to define it, it is useful to introduce the function: by means of which we set: This is a smooth function that obeys eqs. (7.5) and (7.9), and whereẑ 0 ,ẑ 1 , and p are free parameters. The physical meaning of the parametersẑ 0 andẑ 1 is that, for z such thatẑ(z) <ẑ 0 (ẑ(z) >ẑ 1 ), the matched PDF coincides with the recursive (the asymptotic) solution. As a matter of fact, eqs. (7.24) and (7.25) stem from the observation that it iŝ z(z), and not z, the natural variable to carry out the matching, and this is because the large-z behaviour of the PDFs is achieved when logarithmic terms grow much larger than non-logarithmic ones.
In principle, the parametersẑ 0 ,ẑ 1 , and p are unconstrained. In order to choose them sensibly, we plot in fig. 4 the asymptotic and recursive solutions as solid black and red curves, respectively (both are multiplied by a factor of 10 −2 , for reasons that will soon  become clear). For the matching to work reasonably well, the transition between the recursive and the asymptotic solutions must occur in a region where these two predictions are as close as possible to each other (which we interpret as the signal that both give a reasonable description of the "true" photon PDF). From fig. 4, we gather that such a region is 2 ẑ 6; this suggests to setẑ 0 = 2 andẑ 1 = 6. However, it is clear that there is (and there must be) a certain flexibility in these choices. The agreement between the asymptotic and recursive solutions quickly worsens in the range z ∈ (1, 2), which implies thatẑ 0 = 1 must be considered as an extreme choice. On the other hand, forẑ > 6 the asymptotic and recursive solutions do tend to stay relatively close to each other, to the extent that evenẑ 1 = 10 appears to be an acceptable choice. As far as p is concerned, the larger this parameter the more abrupt is the transition between the two regimes. We have therefore opted to employ p = 2, which essentially corresponds to the slowest transition compatible with the derivatives of G(z) being continuous. In order to assess the impact of the choices ofẑ 0 andẑ 1 on the matched PDF, we have computed the latter for several values of these parameters. In fig. 4 we display as dashed curves the differences between any of the matched predictions (relevant to (ẑ 0 ,ẑ 1 ) = (1, 5), (3,5), (1,7), (3,7), and (1, 10)) minus the one obtained with (ẑ 0 ,ẑ 1 ) = (2, 6). For comparison, we also show the differences between the asymptotic and recursive solutions minus the (ẑ 0 ,ẑ 1 ) = (2, 6) matched PDF as black and red dot-dashed curves, respectively. We see that the differences between any two pairs of matched predictions are roughly in the range (−2, 3) · 10 −4 , i.e. at least a factor 25 smaller than the recursive and the asymptotic solutions. While this statement progressively loses validity when moving towardsẑ = 8 (where the asymptotic solution, which is the appropriate one in this region, crosses zero), it also becomes less relevant, since indeed all quantities of interest tend to become negligible in absolute value. Having said that, it is important to bear in mind that the dependence on the matching-function parameters is a genuine uncertainty that affects the matched predictions; plots such as those in fig. 4 help assess its size, and should be re-produced whenever new conditions become relevant (specifically, for hard-scale values significantly different w.r.t. those considered in fig. 4). We finally point out that we have repeated the exercise by using p = 3 and p = 4; the overall uncertainties are similar to those obtained with p = 2.
In summary, our best analytical results are obtained with the matching formula of eq. (7.4). In the case of the singlet and the non-singlet, we employ eq. (7.13) for the definition of the subtraction term, and a constant G(z) function as in eq. (7.11). In the case of the photon, the definition of the subtraction term is still given by eq. (7.13) -however, there are more limited possibilities here, owing to eq. (7.23). The matched photon PDF does need a non-trivial matching function: we adopt that of eq. (7.25), withẑ 0 = 2,ẑ 1 = 6, and p = 2.

Numerical and analytical predictions
In this section we present our predictions for the PDFs, by computing them both with the numerical code described in the first part of sect. 6, and by evaluating the analytical formulae; these are always the matched ones. We compare these two classes of predictions, mutually validating them in the process. Unless explicitly indicated, all results are NLLaccurate with running α, and all have been obtained by setting µ 0 = m.
We begin by plotting in fig. 5 the electron, photon, and positron PDFs, computed at µ = 100 GeV with the numerical code. The left panel shows these quantities in the full  z ∈ (0, 1) range, while the right panel is a zoom to the large-z region, where we consider z ∈ (1, 15) (see eq. (7.24) for the definition ofẑ). Owing to the much faster growth of the electron PDFs in this region w.r.t. that of the other two partons, we have multiplied this PDF by a factor equal to (1 − z), in order for all of the three curves to fit into the same layout. Figure 5 renders it manifest that the production of heavy (relative to the collider c.m. energy) objects is dominated by the partonic lepton whose charge is the same as that of the particle lepton it stems from 15 (in eq. (1.1), one has the implicit constraint z + z − ≥ M 2 /S, with M and √ S the mass of the object produced and the collider c.m. energy, respectively). Note that from the right panel of fig. 5, given that the solid-red and dashed-blue curves are roughly of the same order, and that the former includes the (1 − z) factor, one can immediately see that the photon PDF is smaller than the electron PDF by a number of orders of magnitude equal to the value ofẑ. Conversely, by producing lighter objects and/or by increasing the collider energy, the contribution(s) of the incoming photon(s) become(s) important.
In view of the smallness of the positron PDF as is documented in fig. 5, it is more convenient to present our findings in terms of the singlet and the non-singlet PDFs rather than by means of the electron and positron ones. This is what we shall do in the remainder of this section. In order to establish the level of agreement between our numerical and analytical predictions, we plot in figures 6, 7, and 8 the ratios of the latter over the former, minus one, in the cases of the non-singlet, singlet, and photon, respectively. In each plot, there  are three curves, corresponding to three different choices of hard scales: µ = 0.01 GeV (dot-dashed green curves), µ = 1 GeV (dashed blue curves), and µ = 100 GeV (solid red curves). An overarching observation is that, in all of the cases bar for the photon at large z's (an exception to which we shall return later), the µ = 100 GeV results are those which display the largest analytical-numerical disagreements.
However, even in this worst-case scenario, the level of agreement is excellent, being typically of the order of 10 −5 -10 −4 (relative); the largest disagreements are to be found at small z's in the case of the singlet (because of the presence of un-resummed log z terms 16 ). In keeping with the previous remark relevant to the hard-scale dependence, the case of the photon at z 1 15 The reader must bear in mind that all our results are obtained by assuming an electron particle. In the case of a positron particle, the roles of the electron and positron partons are simply reversed. 16 Techniques to resum such logarithms exist, see e.g. refs. [19,20].   constitutes an exception: from the right panel of fig. 8 we see that the analytical and numerical predictions agree at the level of 10 −3 (10 −2 ) at µ = 100 GeV (µ = 0.01 GeV) for 2 ẑ 6 ; furthermore, the behaviour atẑ > 6 might seem to suggest that the z → 1 limits of the analytical and numerical computations are different. We shall show in the following (see fig. 12) that this is in fact not the case. For the time being, the crucial thing to bear in mind is that, in the z region we are discussing, the photon PDF is very small in absolute value and, more importantly, smaller than the electron PDF by several orders of magnitude. Thus, even a relatively large discrepancy of 0.1-1% between the numerical and analytical photon PDFs will be quite irrelevant. The general conclusion, which applies to all partons, is therefore that the analytical formulae appear to be perfectly adequate, and can be employed in calculations of cross sections for phenomenological purposes.
We now turn to assessing the effects on the PDFs of the NLL corrections, by comparing the NLL results with their LL counterparts. While this will fully account for the predictions obtained here for the first time, it is important to bear in mind that the PDFs are unphysical quantities, and that beyond LL cancellations do occur (in particular, in MS) between them     fig. 9, for the photon. and the short-distance cross sections. Thus, an increase or decrease by a factor X of an NLL PDF w.r.t. an LL one will most definitely not translate into a corresponding increase or decrease of the NLO physical cross section w.r.t. its LO counterpart.
In the main frames of figs. 9, 10, and 11, we plot the ratios of the NLL PDFs over the LL ones, both computed with the numerical code. As was done previously, all figures feature three curves, that correspond to different choices of hard scales; the same scale values, and the same graphical patterns, are used here as in figs. 6-8. All the figures have an inset, which displays the double ratio (minus one): The agreement between the numerical and analytical predictions is again extremely good, especially at large z's; once more, the photon in this region is the (relative) exception to that general rule, on which we shall comment later. The agreement becomes marginally worse with increasing µ, but this effect is less evident w.r.t. that in the case of the absolute predictions of figs. 6-8. Interestingly, the size of the NLL effects decreases with the hard scale. This is particularly easy to understand in the large-z region in the case of the singlet (or non-singlet), since it can be directly read from eq. (5.63). As was already remarked there, this behaviour is driven by eq. (5.64), which implies that: a) the coefficient of the log(1 − z) term is much larger than that of the log 2 (1 − z) term up to extremely large values of z; b) such a coefficient, being proportional to 1/α(µ), decreases with µ. These two effects can clearly be seen in the main frames of the right panels of figs. 9 and 10, where the various lines are almost straight ones, but relatively less so at larger values of the hard scales. Keeping in mind the general observation made above on the unphysical nature of the PDFs, we point out that the natural applications of the quantities computed here involve scales that are large. We now go back to commenting on the large-z behaviour of the photon PDF. We have already remarked in fig. 8 that such a PDF in this region is close to zero in absolute value, and orders of magnitude smaller than its electron counterpart. On top of that, for the specific issue of the NLL vs LL results, the comparison between eqs. (B.74) and (B.88) (or between their over-simplified forms of eqs. (B.87) and (B.91)) shows that, at variance with the case of the electron (eqs. (5.63) and (5.66)), the NLL asymptotic photon PDF does not factorise the functional form relevant to its LL version. Hence, larger differences in the matching region have to be expected between the NLL and LL photon PDF, which are larger than those for the electron.
At the right end of the z range in fig. 11 we see again the kind of pattern as in the same region of fig. 8, which might cast doubts on the agreement between the z → 1 limits of the analytical and numerical predictions. In order to address this concern, in fig. 12 we plot the photon PDF in a much more extended z range w.r.t. what was done so far. The blue and red solid curves are the differences between the analytical and numerical results computed at the LL and NLL, respectively; the dashed curves of the same colours are the corresponding PDFs, multiplied by an overall constant factor equal to 10 −3 ; finally, the blue dot-dashed curve is the rescaled ratio of the analytical over the numerical LL results, minus one, which can be sensibly computed owing to the fact that the LL PDF does not vanish for values of z = 1. Apart from the similarity between the LL and NLL differences, the key message of fig. 12 is that at z → 1 the analytical and numerical predictions do tend to the same value, but in a much slower way w.r.t. the case of the singlet/non-singlet. In other words, the onset of the true asymptotic regime occurs at much larger z values for the photon than for the singlet or non-singlet. This needs not be surprising, owing to the mechanism that governs the asymptotic photon behaviour, as is documented in appendix B. An improvement of the analytical large-z PDF computed here would require keeping all terms suppressed by powers of N −2 in Mellin space, an extremely involved computation which is not justified in view of the smallness of the photon PDF in this region.

Conclusions
In this paper we have computed the electron, positron, and photon PDFs relevant to an incoming unpolarised electron; the case of an incoming positron is trivially obtained by charge conjugation. Our predictions include up to next-to-leading logarithmic (NLL) terms, and are obtained by evolving the initial conditions that have recently been calculated [1] at the next-to-leading order (NLO). We thus improve upon the long-standing leadinglogarithmic (LL) PDFs of refs. [2][3][4]; this is crucial to help achieve the high-precision results needed at future e + e − colliders. All of the calculations are perturbative in the QED coupling constant α, which by default we take as running.
The PDFs have been obtained by means of both numerical and analytical methods, which are shown to agree extremely well with each other (typically, at the 10 −4 level or better for those z values where the relevant PDFs are not vanishingly small). As far as the analytical results are concerned, they stem from an additive matching formula (eq. (7.4)), which combines a prediction that is accurate to all orders in α but only at z → 1, with a prediction that is accurate up to O(α 3 ) in the whole z range; these are referred to as the asymptotic and the recursive solutions, respectively. The latter are thus called because they stem from recursive equations (derived here for the first time at the NLL accuracy), whereby the O(α k ) contributions are obtained from the O(α p ) ones (with p < k), the NLO initial conditions, and the Altarelli-Parisi kernels. Although we have limited ourselves to considering k ≤ 3 in this work, nothing prevents one from employing the recursive equations to achieve an even higher precision if need be.
The electron LL PDF is extremely peaked towards z → 1, where it features an integrable singularity. The NLL result has the same qualitative behaviour, and in fact the z → 1 singularity is even more pronounced than at the LL because of the presence of additional log(1 − z) terms. While the photon LL PDF vanishes at z → 1, its NLL counterpart grows logarithmically there, thus exhibiting the same enhanced growth at higher orders as the electron. However, one must bear in mind that PDFs are unphysical quantities: in particular, beyond LL they become dependent on the adopted collinear subtraction scheme. In this paper we have worked in MS, and some of the logarithms mentioned above stem directly from this scheme choice; as such, they will cancel against their counterparts in the short-distance cross sections, so as to have scheme-independent predictions for physical observables.
The analytical knowledge of the PDFs is important to better understand the details of QED collinear dynamics. However, in view of the rapidity of the growth of the electron PDF at z → 1, such a knowledge is also crucial in the context of numerical computations, because it gives one the possibility to adapt in a very tailored manner the integration procedure, which would otherwise be hardly converging.
We finally point out that the PDFs of the photon (understood as an incoming particle), and/or those of any polarised particle, can be dealt with similar techniques as those we have employed here.

A Results for the recursive solutions
In this appendix we report the results for the J LL k and J NLL k basis functions that appear in eq. (5.35), which we have computed for 0 ≤ k ≤ 3 and 0 ≤ k ≤ 2, respectively, i.e. up to O(α 3 ). We write the actual recursive solution that constitutes one of the main results of this paper, and which we have used in our numerical studies, as follows: We also remind the reader that from eq. (A.1) one can obtain the solution in the case of non-running α, by replacing J LL k with I LL k and J NLL k with I NLL k , where: It is convenient, also in view of the matching with the large-z solution, to present the results for the basis functions by writing them as follows: By definition,Ĵ LL k andĴ NLL k collect all of the terms of J LL k and J NLL k , respectively, that vanish at z = 1: lim It then follows thatJ LL k andJ NLL k include all contributions that are either divergent (which then feature all the log p (1 − z) terms) or equal to a non-null constant at z = 1. Because of this, it is useful to introduce the following auxiliary functions: and write:J In addition to this, one must take into account that, at O(α 0 ): The contribution to Γ(z) that does not vanish at z → 1 is then written as follows: coefficients for the non-singlet, singlet, and photon PDFs will be presented in appendix A.1 (LL results), and in appendix A.2 (NLL results). The expressions of the functionsĴ (N)LL (z) are lengthy (with some of them receiving contributions that we have not computed analytically, as detailed below), and not relevant to the matching; for these reasons, they are only reported in an ancillary file that will accompany the submission of this paper to the arXiv. We remind the reader that the recursive solutions are obtained by following the procedure outlined in sect. 5.1. Namely, one first computes the J LL k and J NLL k functions, by employing eqs. (5.30) and (5.31). These equations must be applied recursively, by working one's own way up in k from the known k = 0 results (given in eqs. (5.10)-(5.13)). The expressions for the Altarelli-Parisi kernels are taken from ref. [21]. Finally, the J LL k and J NLL    included in the recursive solutions, we plot in fig. 13 the ratio of the result of eq. (A.1) over the numerical predictions minus one; eq. (A.1) is computed by setting: The ratios are displayed as green dot-dashed lines (k LL max = 1), blue dashed lines (k LL max = 2), and red solid lines (k LL max = 3). In order for the results to fit into the layout of the figures, the green and blue curves are multiplied by a constant factor equal to 10 −2 and 10 −1 , respectively. In keeping with what has been discussed in sect. 8, we see that our most accurate recursive predictions (k LL max = 3) agree with the numerical results at the level of a few 10 −4 at the worst. Note that since here we are dealing only with the recursive solutions we have limited ourselves to plotting the PDFs in the range z ∈ (0, 0.9) -at the upper end of the range, the absence of the contribution from the asymptotic solution starts to be felt. The new information stemming from fig. 13 is that, if we had only computed either the first term or the first two terms in the sums of eq. (A.1), the O(10 −4 ) agreement remarked above would actually have been roughly equal to, but generally worse than, 10 −2 and 10 −3 , respectively. The figure also shows that, for any given accuracy of the recursive solution, the agreement with the numerical prediction marginally worsens towards z → 0 in the case of the singlet, owing to the presence of log z terms which are not resummed.
In the course of the recursive procedure, we have found that some integrals relevant to J NLL 2 (i.e. the function associated with the O(αt 2 ) term in the representation of the PDFs) are not easily computable analytically. We have therefore opted to limit ourselves to obtaining their z → 1 leading terms analytically, while evaluating all of the remaining terms numerically, so that the latter contribute only toĴ NLL 2 (we point out that an analogous strategy has already been adopted in ref. [2]). More precisely, let us consider the generic modified-convolution integral of eq. (5.3). We distinguish two possibilities: either g(x) is a plus distribution, or it is an ordinary function. Notation-wise, these two cases are written as follows: In the case of eq. (A.17), we have: where we have defined the endpoint and bulk contributions, respectively, as follows: These equations can also be used in the simpler case of eq. (A.18): one simply sets the endpoint contribution equal to zero, and computes eq. (A.21) by removing the subtraction term h(z) and with the formal replacementĝ → g there.
The endpoint contribution of eq. (A.20) is always computed analytically, and its results are included inJ NLL 2 (z) and/orĴ NLL k (z), according to the behaviour of h(z) at z → 1. As far as eq. (A.21) is concerned, for the sake of the forthcoming discussion let us re-write it more compactly as follows: If the integral in eq. (A.22) were strongly convergent, then we might obtain its contribution to the PDF (see eq. (5.2)) by means of a derivation under the integral sign, namely: Unfortunately, the strong convergence of the integral is not guaranteed, given that F (z) in general is logarithmically divergent at z → 1. However, the contributions that are non vanishing at z = 1 are also easy to compute analytically; such computation can be carried out directly at the differential level of eq. (A.23), and stems from expanding the integrand on the r.h.s. of that equation in a series of z around 1. The latter must include all terms that result in either a logarithmically-divergent or a constant non-null term at z → 1, which typically implies up to (1 − z) 0 contributions. In this way we arrive at the following identity: with: ∂F (z) ∂z asy = Conversely, the quantity in square brackets in eq. (A.24), where the rightmost term is regarded as a regularising factor, is computed numerically 17 , and eventually included in We list here the pairsĝ (or g) and h that we handle in the way we have just described: We stress again that each of these pairs will contribute to both eq. (A.26) and (A.27). We denote generically either of these contributions as follows: These will enterJ NLL 2 (z) andĴ NLL 2 (z) as linear combinations with identical coefficients (owing to eq. (A.24)), which however do depend on the flavour structure. Explicitly: non − singlet : (A.56) As was anticipated, the results on the r.h.s. of eqs. (A.54)-(A.56) do not contribute to any of the b NLL 2,i coefficients, while they enter the coefficients c NLL 2,1 and c NLL 2,0 (singlet and non-singlet), and c NLL 2,3 , c NLL 2,1 , and c NLL 2,0 (photon).

A.1 LL coefficients
In this appendix we report the results for the coefficients that enter eq. (A.10); all of the coefficients that do not appear below are understood to be equal to zero.

A.2 NLL coefficients
In this appendix we report the results for the coefficients that enter eq. (A.11). Note that these do already include the r.h.s. of eqs. (A.54)-(A.56); all of the coefficients that do not appear below are understood to be equal to zero. We employ the following shorthand notation: B Asymptotic large-z solution for photon PDF beyond leading N In this appendix we consider the problem that has been anticipated in sect. 5.2.3, namely the improvement of the asymptotic behaviour of the photon PDF given in eq. (5.73) stemming from the inclusion of the off-diagonal elements in the evolution kernels. In order to do this, we start from writing the O(α) expressions of the Altarelli-Parisi kernels as follows (see eq. (4.5)): having introduced, at each order in α, the leading-and subleading-N contributions. They read as follows: Note that, by considering only eqs. (B.3) and (B.5), one recovers eq. (5.71). According to eq. (4.15), the Altarelli-Parisi kernels define the evolution kernel as follows: whence one can write the evolution equation and its formal solution as follows: The solution in eq. (B.8) is based on the so-called Magnus expansion [22] (see also ref. [23]), which is constructed solely in terms of the evolution kernel: with Ω k,N (t) featuring k instances of M N , all appearing in commutators. Thus, in the case of a one-dimensional flavour space or of a diagonal evolution kernels, eq. (B.8) is identical to the solutions given in eq. (4.17) and in sect. 5.2.3. As far as the singlet-photon sector is concerned, we can indeed recover the solutions we have found previously in terms of the quantity introduced in this appendix. We define the leading-N evolution kernel: and denote by E   where: ΣΣ,N = exp −ξ 1 logN +ξ 1 , (B.14) .  12)). This is useful when one considers the limit of non-running α of the formulae presented here. An expression equivalent to eq. (B.15), as well as the QED "limit" of both, is given in eq. (B.16).
We stress that the case of non-running α is problematic, as it might lead to inconsistencies. By switching off the running, one effectively neglects bubble-diagram contributions which are exactly the same as those that lead to the γγ entries in eqs. (B.3) and (B.5). In this paper we ignore such potential inconsistencies, but then we need to carefully distinguish the γγ contributions to the Altarelli-Parisi kernels (which we always parametrise by means of n F ) from those to the QED β function (which we parametrise by means of the β-function coefficients b i ). We shall return to this point with one explicit example later in this appendix (see eqs. (B.75) and (B.76)).
In order to improve on the leading-N results, we shall introduce the subleading-N contributions to the evolution kernel, and treat them as a perturbation to the solution of eq. (B.13). This entails writing: having defined: N there. We then observe that Ω k,N ∝ 1/N k , and thus for consistency with eq. (B.2) we are allowed to discard all contributions with k ≥ 2. Therefore:  N . Apart from rendering the t 1 integral in eq. (B.20) non trivial, this will also induce functional forms in the N -space whose analytical inverse Mellin transforms will be extremely hard to compute. We shall therefore resort to simplifying the expression of E where: In equation (B.23) we have introduced the quantities ξ 1,0 andξ 1,0 which we have defined as follows: with ξ 1 andξ 1 given in eqs. (5.51) and (5.53). By means of an explicit computations from the latter two equations we obtain: As far as eq. (B.24) is concerned, its expression stems from that of eq. (B.15); in particular: from whence: In summary, the evolution operator we shall use is the following: Having established that the asymptotic solutions presented in sect. 5.2 are perfectly adequate for the case of the singlet, we shall now focus on the implications of eq. (B.30) on the photon PDF. We obtain: where: 1 } four sets of N -independent quantities, whose specific forms are unimportant here. For any given j, the five terms in the numerators on the r.h.s. of eq. (B.37) can be re-expressed algebraically (i.e. without any approximations) in terms of the corresponding denominators. In this way, one arrives at the following forms (note that E (0,L) γγ,N is independent of N ): Here, we have introduced the inverse Mellin transforms relevant to eq. (B.37) which are linearly independent from each other, namely 18 : Explicit computations give: 18 Here and elsewhere, some quantities are denoted to depend on parameters which do not actually enter their functional forms. This allows one to write eq. (B.38) in a compact, and formally correct, way.
where, consistently with eqs. (B.39) and (B.40), in eqs. (B.41)-(B.45) some terms that vanish at z → 1 have not been included. This is of course arbitrary to some extent, and the logic we have followed is that of keeping those terms which, when expanded in series, either contribute to the same monomials t n and αt n as the recursive solutions considered in this paper, or have the same power of κ as the former ones. On top of this, one has the special case of eq. (B.41) which has the structure of a series in M k−1 When z → 1, these terms are progressively more suppressed with increasing k. Unfortunately, this hierarchy is not valid at intermediate z's; in fact, for the values of M 1 and M 2 relevant to our computation there is a singularity at z 0.65 which is dominated by increasingly large values of k. This is what prevents the asymptotic solution of the photon PDF from being well-behaved in all of the z range, at variance with its electron counterpart. This has significant implications for the matching, which are discussed in sect. 7. A solution to this problem would be that of resumming the series on the r.h.s. of eq. (B.41); we have computed its first seven coefficients, but have not been able to identify the corresponding generating function. Numerically, the use of all of these seven contributions instead of the three reported in eq. (B.41) does not change the behaviour at large z's, and does not improve that at intermediate z's.
The R i functions that appear in eq. (B.38) are:  , and to eq. (B.33), it is immediate to see that this contribution, up to terms vanishing in the z → 1 limit, is identical to that of eq. (5.73), bar for an α(µ 0 )/α(µ) prefactor that here needs to be written according to eq. (B.24). Thus, by introducing the quantity: we can write the sought large-z expression of the photon PDFs in a compact form: with Γ γ,j (z) given in eq. (B.38) for j ≤ 4 and in eq. (B.73) for j = 5. The results presented above allow one to obtain their counterparts in the case of nonrunning α, by means of the following formal replacements (see eq. (4.19)): with η 0 defined in eq. (5.6). This procedure is consistent with its analogue relevant to the recursive solutions (see sect. 5.1 and appendix A). We can also see that, by using the replacements above in the expression for E That being said, we point out again that the limit of non-running α must be interpreted with some care (see the comments that follow eq. (B.16)).
When not considering the case of non-running α, one can re-expressed the exponential prefactors in eq. (B.74) and in eqs. (B.51)-(B.70), and their combinations, in simpler ways, namely: Equation (B.74) is the asymptotic solution that emerges from solving the evolution equation by keeping the dominant off-diagonal terms in the Altarelli-Parisi kernels. As we shall discuss in appendix C, it shares with its singlet and non-singlet counterparts the nice property that its perturbative expansion lead to the same coefficients as those of the recursive solutions (for certain classes of basis functions in the z space). However, its functional form is rather involved, but it is fortunately possible to simplify it, by keeping only the truly dominant terms in the z → 1 limit at each order in α. In order to do so one starts by observing that, in such a limit, one has: There is a certain similarity between eq. (B.87) and eq. (5.63) which is worth stressing. In particular, the dominant term at z → 1 in both equations (proportional to log(1 − z) 3 and log(1 − z) 2 , respectively) is suppressed w.r.t. the subdominant one (log(1 − z) in both cases) by a factor proportional to α (owing to eq. (5.64) for eq. (5.63)). This implies that numerically the onset of the behaviour driven by the most divergent terms occurs only at z values which are exceedingly large, and in fact hardly relevant to any phenomenological applications -we have commented further on this point in sect. 8. Equation (B.74) simplifies considerably when one retains only the LL terms. A direct calculation leads to the following result: with ξ 0 andξ 0 defined in eq. (5.67), and: We point out that, consistently with the results of appendix A.1, the LL photon PDF is of O(t) (i.e. it does vanish with α → 0): the two terms on the r.h.s. of eq. (B.88) cancel each other at t = 0. From eq. (B.41), we also see that the LL-accurate photon PDF of eq. (B.88) vanishes in the z → 1 limit: By comparing eqs. (B.87) and (B.91) we observe that the photon PDF has a behaviour analogous to that of the electron PDF, namely that its NLL form grows faster than its LL counterpart at z → 1; to a good extent, this is an artifact of the MS scheme. We finally point out that eq. (B.88) can be directly obtained from solving the evolution equation of eq. (B.8), by using there: Since the kernel of eq. (B.92) is independent of t, eq. (B.8) can be simply solved by diagonalisation. After that, one multiplies the results by the LO initial conditions, and performs the inverse Mellin transform. The fact that by doing so one recovers eq. (B.88) is a rather powerful check on the procedure adopted in this appendix.

C Expansion of large-z solutions
In view of the matching between the asymptotic large-z solutions and the recursive solutions, it is useful to consider the expansion of the former ones in a series of α; this will also allow us to perform some consistency checks on them. We can formally write the result of such an expansion for the NLL-accurate, running-α solutions of eqs.
As the notation with an overline suggests, we only take into account contributions that do no vanish at z → 1. We point out that we consider the expansion up to O(α 3 ), i.e. we use the values in eq. (A.2), for the sole reason of consistency with what has been done for the recursive solutions in appendix A. The flavour structure of eq. (C.1) is the same as that in eqs. (7.2) and (7.3), and can therefore be accounted for by K LL k and K NLL k , precisely as is the case of the J LL k and J NLL k functions for the recursive solutions; in practice we shall omit flavour indices here, in order to simplify the notation, since no confusion is possible. In fact, one must bear in mind that the large-z solutions of the singlet and non-singlet PDFs coincide, and that the one of the photon has a functional behaviour significantly different from the former two. Therefore we shall first deal with the singlet non-singlet cases together, and with that of the photon afterwards.
We have determined the coefficients A LL k and B LL k,i for k ≤ 3, and A NLL k and B NLL k,i for k ≤ 2, by means of a direct computation. The results for k = 0 are particularly interesting since, in view of eq. (C.1), they must be related to the initial conditions of eqs. (3.18) and (3.19). We have obtained: where L 0 has been defined in eq. (A.77). With the result of eq. (C.9), K LL 0 (z) is indeed identical to eq. (3.18). However, by replacing the results of eqs. (C.10)-(C.12) into eq. (C.5), K NLL 0 (z) turns out not to coincide with Γ [1] e − (z) of eq. (3.19). This is hardly surprising: when working in the large-z region, one is entitled to set z = 1 in all of the polynomial terms that appear in the numerators. Therefore, while K NLL 0 (z) should not necessarily be equal to Γ [1] e − (z), it must be equal to the z → 1 asymptotic form of the latter -if that were not the case, the large-z solution would not be compatible with the initial conditions from which it supposedly originates. In order to obtain the asymptotic expression of the initial condition, one cannot set z = 1 in all of the numerators of the latter right away, since Γ [1] e − (z) is not an ordinary function, but a distribution. Before doing so, one must first pull out the 1 + z 2 factors from the plus distributions in eq. (3.19). This can be done by exploiting the following identities: (C.14) After having done this, one can finally let 1 + z 2 → 2 in the numerators. It is a matter of simple algebra to show that this procedure leads to the expected result: In summary, we have thus proven that the solution of eq. (5.63) embeds the initial conditions of eqs. (3.18) and (3.19).
We conclude this appendix by reporting the results for the coefficients with k > 0. We have obtained what follows: 20) and: with b LL S, k,i and b LL NS, k,i given in appendix A.1, and b NLL S, k,i and b NLL NS, k,i in appendix A.2. We point out that eqs. (C.21) and (C.22) hold for all values of k and i we have considered here. This is remarkable, because it tells one that with the expressions obtained in this paper all of the log p (1 − z)/(1 − z) terms in the PDF are the same regardless of whether one obtains them from the recursive solution, or by expanding the asymptotic solution. In general, one expects the logarithms from the latter to coincide with those of the former only for the larger values of p at any given k. The result obtained here ultimately stems from keeping some formally subleading contributions in the procedure of sect. 5.2.2; in particular, it is important that the numerators in eqs. (5.56) and (5.57) be 1 + x 2 rather than 2 (which would be a perfectly fine choice in the asymptotic region) 20 .

C.2 Photon
In the case of the photon one needs to employ eq. (B.74). We start by observing that the Taylor series in t and α of such a quantity leads order by order to integrable singularities; as expected, there is therefore no endpoint contribution, and the expansion of the large-z solution can be expressed in terms of ordinary functions. Before turning to the explicit form of the latter, we point out that the t 0 term in the expansion of eq. (B.74) is equal 20 It turns out that the use of 1 + x 2 is essential in the determination of the endpoint contributions in the plus distributions of eqs. (5.56) and (5.57), which in turn induce (some of) the z-independent terms in eqs. (5.58) and (5.59). Conversely, away from the endpoints the replacement of 1 + x 2 with 2 leads to power-suppressed terms at z → 1.
to Γ γ,5 (z), since the contributions of the Γ γ,j (z) with j ≤ 4 terms mutually cancel (that of j = 1 (j = 3) against that of j = 2 (j = 4)). One thus recovers the initial conditions of eqs. (3.18) and (3.20), which is a first consistency check on eq. (B.74). We now write the expansion of the large-z photon PDF in the same way as was done in eq. (C.6), but with the K k functions defined as follows: having introduced the q i (z) functions in eq. (A.9). It is a matter of algebra to arrive at the final results: with c LL γ, k,i given in appendix A.1, and c NLL γ, k,i in appendix A.2. As was the case for their singlet and non-singlet counterparts (eqs. (C.21) and (C.22)), eqs. (C.25) and (C.26) have the property of holding for all of the k and i values considered here. Thus, the same remarks done previously are valid here as well (with the obvious exception that they apply to the log p (1 − z) terms rather than to the log p (1 − z)/(1 − z) ones relevant to the singlet and non-singlet cases).

D Alternative z-space derivation of asymptotic large-z solutions
In this appendix we show how some of the asymptotic results of sect. 5.2 can be obtained directly in configuration space, that is without resorting to Mellin-space techniques, and thus providing one with a cross-check on the results of the latter. We have considered this alternative procedure starting from a couple of simplifying assumptions 21 : namely, we only deal with the non-singlet case, and we neglect the running of α. We point out that this method has already been used to obtain the LL solution of eq. (5.46) -see e.g. ref. [24]. Here, we extend it to the NLL accuracy.
In essence, the procedure works as follows. One makes an ansatz for the z-space functional form of Γ(z, µ 2 ), where the µ 2 dependence is parametrised by unknown functions. The PDF evolution equations, simplified in the z → 1 limit, are then turned into differential equations for such unknown functions, where the independent variable is µ 2 . By solving these equations, one is left with arbitrary integration constants, whose values are finally determined by matching the solutions to the known PDF initial conditions.
In order to proceed, we start by observing that the assumption of non-running α implies that the dependence on µ 2 can be entirely parametrised by means of the quantity 21 We do not make any claims as to whether this z-space approach remains viable if either of these assumptions is relaxed. η 0 , introduced in eq. (5.6); thus, we shall use the latter as our independent variable. At the LL, this implies that the evolution equation of eq. (3.8) reads as follows: one finds a system of differential equations: = e(η 0 ) c 0 − d 0 ψ 0 (η 1 + a 0 + 1) + ψ 0 (η 1 + a 0 + 1) 2 − ψ 1 (η 1 + a 0 + 1) , (D. 24) where η 1 andη 1 have been defined in eqs. (5.69) and (5.70), respectively. The arbitrary integration constants a 0 , . . . e 0 can be found by matching with the initial condition. We observe that at µ = µ 0 the α → 0 NLL result for the PDF must coincide with the LL one; this implies that eq. (D.9) must still hold true. Because of this, one can expand eq. (D.12) by using the techniques employed in appendix C (see in particular eq. (C.2)), to obtain at O(α) the same functional form as in eq. (3.19), which leads to the following results: