On factorisation schemes for the electron parton distribution functions in QED

The electron, positron, and photon Parton Distribution Functions (PDFs) of the unpolarised electron have recently been computed at the next-to-leading logarithmic accuracy in QED, by adopting the $\overline{\rm MS}$ factorisation scheme. We present here analogous results, obtained by working in a different framework that is inspired by the so-called DIS scheme. We derive analytical solutions relevant to the large-$z$ region, where we show that the behaviour of the PDFs depends in a dramatic way on whether running-$\alpha$ effects are included to all orders, as opposed to being truncated to some fixed order. By means of suitable initial and evolution conditions, next-to-leading logarithmic accurate PDFs are obtained whose large-$z$ functional forms are identical to those of their leading logarithmic counterparts.


Introduction
The computation of e + e − cross sections in perturbative QED leads to the emergence of logarithms of the ratio m 2 /E 2 , where m is the electron mass and E a scale of the order of the hardness of the process, such as the collider energy. Owing to the smallness of the electron mass, these logarithms are numerically very large in virtually all situations of phenomenological interest, in particular at the future colliders presently under consideration. This is problematic, since it implies that the perturbative series is not well behaved, with coefficients that grow with the power of α. It is therefore mandatory to re-sum such logarithms in order to obtain meaningful physical predictions. While there are different sources of large logarithms, possibly dependent on the observables one considers, for some classes of them, associated with soft and/or collinear radiation off initial-and final-state particles, general resummation techniques exist, such as YFS [1,2], parton shower [3][4][5][6], and factorisation formulae [7,8] (the latter sometimes improperly referred to as structure function approach).
Factorisation formulae are fully analogous to their QCD counterparts and, as in that case, resum large logarithms of collinear initial-and final-state origin by means of parton distribution functions (PDFs) and fragmentation functions, respectively. At variance with what happens in QCD, these quantities are computable in perturbation theory. As far as the PDFs are concerned, which are the focus of this paper, results of leading-logarithmic (LL) accuracy [9][10][11] have recently been extended to the next-to-leading logarithmic (NLL) accuracy [12,13] 1 . Among other things, this extension has a strong phenomenological motivation, in view of the precision targets relevant to future colliders.
At the NLL and beyond, PDFs depend on the choice of the factorisation scheme. In essence, this is the choice of the finite part in the subtraction of the remaining (after the cancellations stemming from the definition of IRC-safe observables) initial-state collinear singularities, the residue of whose poles in 1/¯ is an Altarelli-Parisi kernel. Such a finite part enters both the short-distance cross sections and the PDFs so that, when both are expanded in α at the same order, it cancels exactly. However, since in practical applications the PDFs must not be expanded (because otherwise they lose the capability of resumming the logarithms), a dependence on the scheme choice is present, although beyond the accuracy of the computation.
While the factorisation-scheme dependence is suppressed by powers of α w.r.t. terms which are under formal control, one may wonder whether numerically this suppression is effective, i.e. whether the coefficients that multiply these α k factors can grow pathologically large. Naively, the MS result of ref. [13] would seem to imply that this is the case. In the z → 1 region, which gives the dominant contribution to the cross section, the electron and photon densities (inside an electron) read as follow 2 : whereas their LL counterparts do not feature any of the log(1 − z) terms that appear on the r.h.s. of eqs. (1.1) and (1.2). However, this is indeed naive. Firstly, the PDFs are unphysical objects: the actual impact of these logarithms must be assessed after the convolution with the cross sections. Secondly, the leading behaviour of the electron PDF is in any case given by the prefactor (i.e. the first z-dependent term on the r.h.s. of eq. (1.1)), that appears at the LL as well and that has an integrable singularity much stronger than logarithmic. However, the question of the impact of the residual scheme dependence on physical observables remains relevant. Since such a dependence is not parametrical, the best way to assess it is that of comparing results obtained in different factorisation schemes; while not completely conclusive, this allows one to see in practice at which level of accuracy the scheme dependence may constitute a phenomenological problem.

Evolution operator
In this section we write the evolution equations for the PDFs [18][19][20][21] in the same way as it has been done in ref. [13], namely by introducing a PDF evolution operator, and by defining it so as it works in a generic factorisation scheme. In the interest of a shorter notation, throughout this paper eq. (x.y) of ref. [13] will be denoted by eq. (II.x.y). The procedure we follow here has significant overlaps with standard works in QCD (e.g. refs. [22,23]); however, a peculiarity of the QED case warrants a compact re-derivation.
The evolution equations for a given particle type 3 are written in Mellin space as in eq. (II. 4.6): Here, we have denoted by Γ N a column vector that collects all of the relevant PDFs:

2)
3 That is, an electron, a photon, and so forth; in this paper we shall consider explicitly only the case of an electron, but the treatment presented here is general. The reader is encouraged to check sect. 2 of ref. [12] for the definitions of particles and partons used in the present context.
with Γ α i the PDF of parton α i inside the particle of interest. Conventionally, we shall label the partons so that α 1 coincides with the identity of the particle (thus, for an electron particle, α 1 = e − ). The PDFs are regarded as evolved from the values they assume at a given (small) scale µ 0 (initial conditions) by means of an evolution operator (an n × n matrix), thus: Equation (2.1) is then fully equivalent to: More conveniently, in order to take the running of α into account in an easier manner, one uses the variable t of eq. (II.4.10) instead of the scale µ: to re-express eq. (2.4) as follows: In order to simplify the upcoming discussion, we consider the convolution of the PDFs with the short-distance cross sections for only one of the two incoming partons; it should be clear that this entails no loss of generality. Collecting the subtracted partonic cross section again into a column vectorŜ each element of which thus corresponds to a partonic channel, the particle-level cross section in Mellin space reads as follows: The RGE invariance of the cross section is therefore: where b is the power of α that factors out at the Born level, and k is the accuracy of the computation (k = 0 being the LO, k = 1 the NLO, and so forth -in practice, in what follows we shall understand the factor α b , and restrict ourselves to the NLO case k = 1, thus writing O(α 2 ) to denote terms beyond NLO accuracy). So far, we have implicitly assumed to work in the MS factorisation scheme; this is in keeping with what has been done in ref. [13], and applies in particular to the RGEs of the evolution operator. In order to employ a different scheme, we start from the partonic cross sections. In the FKS formalism [24,25] 4 the arbitrariness of the definition of the factorisation scheme is parametrised in terms of a set of generalised functions K ij (z) that enter the so-called (n + 1)-body degenerate cross sections, whose contributions to an NLO partonic cross section can be written as follows in the Mellin space (see e.g. eqs. (3.21)-(3.23) of ref. [26]):σ (2.10) Here, we have included an index "K" to make it explicit the fact that, in general, we work in a factorisation scheme different from MS; we point out that these formulae still encompass the MS case, which corresponds to choosing K ij (z) ≡ 0. By collecting the K ij (z) functions into an n × n matrix: the matrix form of eq. (2.10) reads: As the rightmost side explicitly indicates, the two forms ofŜ (K) N that appear in eq. (2.12) are strictly equivalent at the NLO.
The elementary nature of the asymptotic states in QED renders it possible to compute perturbatively the initial conditions for PDF evolution. The NLO results of ref. [12] for such initial conditions are given in a generic factorisation scheme, and feature the same functions K ij (z) that appear above (for the electron particle, see in particular eqs. (4.121) and (4.189) there). The scheme-dependence of the PDF initial conditions is a beyond-LO effect, and therefore at the LO we can write (by keeping in mind the conventions established in eq. (2.2)): (2.13) It is then easy to see that the NLO results of ref. [12] can be written as follows: where the rightmost side follows from the same considerations as those employed in eq. (2.12). We point out that, owing to eq. (2.13), at the NLO only the first column of the K matrix gives a non-null contribution; in other words, only the K α i α j elements with α j equal to the particle identity (i.e. α 1 ) are relevant. This is in keeping with the elementary nature of that particle, so that at the NLO it is the branchings α 1 → α i + X that are sufficient to control the factorisation scheme. Thus, at this level of precision, the choice of the K α i α j elements with α j = α 1 is essentially arbitrary, since they are associated with NNLO contributions. Note that this is not true for the partonic short distance cross sections, where those functions may contribute at the NLO as well; this is however irrelevant as far as PDF evolution is concerned, which is our goal here. In a generic scheme, we write the analogue of eq. (2.3) with a minimal change of notation, namely: from which the analogue of eq. (2.8) follows: The factorisation-scheme independence of any physical observable is then equivalent to requiring that: In keeping with the fact that it stems from eq. (2.18), eq. (2.20) is understood to hold up to O(α 2 ) terms. In fact, at µ = µ 0 we obtain, from the initial conditions in eqs. (2.3) and (2.16): Although this is within the accuracy at which we are working, it also suggests that the identity of eq. (2.20) can be conveniently modified by observing that, for perturbation theory to work, there must exist a non-zero measure set in the space of the Mellin variable where the following inverse operator exists and is well defined: Since the first two terms in the series on the r.h.s. of eq. (2.22) coincide with the terms within the leftmost brackets in eq. (2.20), it follows that at O(α 2 ) we can also write: which we shall consider as valid to all orders in α (note that the analogue of eq. (2.21) in this case would read I = I). The sought evolution equation for the E (K) N operator can be obtained by deriving both sides of eq. (2.23) w.r.t. the variable t, and by exploiting eq. (2.6). In order to do that, we observe that eq. (2.22) implies: from whence: This result generalises the MS one in the first line of eq. (2.6) and, as it should, coincides with it when K = 0. An O(α 2 )-accurate expansion then leads to: which generalises the second line of eq. (2.6) and coincides with it when K = 0. As consistency demands, eq. (2.26) can also be obtained directly from eq. (2.20), by deriving w.r.t. the variable t the two sides of that equation. The solution of the RGE for the evolution operator that stems from eq. (2.26) is the analogue of what is usually done in QCD (see e.g. refs. [22,27]), whereas to the best of our knowledge the QCD counterpart of eq. (2.25) has never been considered. We shall show in the following that this is justified: while the difference between these two equations is of NNLO in both QCD and in QED, in the latter case it leads to dramatically different behaviours in the z → 1 region; this does not happen in QCD.

A DIS-like factorisation scheme
The RGEs for the evolution operator, eqs. (2.25) and (2.26), are valid in any factorisation scheme; here, we shall make a definite choice for one such scheme.
The electron PDF initial conditions at the NLO [12], namely: show that it is straightforward to choose the K ij (z) functions so that the O(α) contributions at µ 0 = m vanish, thus making the NLO initial conditions to coincide with their LO counterparts. This condition is reminiscent of the definition of the DIS scheme in QCD [16]; in fact, in view of the stricter similarity between the case of the electron PDFs in QED and that of the b-quark perturbative fragmentation function in QCD [28], we shall call the scheme thus defined as ∆ scheme, for consistency with the notation of ref. [27]. From eqs. (3.2) and (3.3) we thus define: We point out that eqs. (3.5) and (3.6), in addition to fulfilling the requirement above, are also consistent with the momentum-conservation condition for the branchings of an electron (see eq. (5.5) of ref. [12]). The factorisation-scheme matrix K of eq. (2.11) is then defined as follows 5 : As was discussed in sect. 2, the matrix elements of K whose column index is not equal to the identity of the particle of interest (the electron here, first column) give beyond-NLO effects. By setting them equal to zero as in the second column of eq. (3.7) (that corresponds to the photon), we impose conditions that are identical to what one had to enforce if one were to define a DIS-like scheme for the photon particle, per eqs. (5.42) and (5.47) of ref. [12]. As far as the positron is concerned (third column), its matrix elements are motivated by charge-invariance considerations. Although the positron plays an essentially trivial role here, the form of eq. (3.7) leads to simple analytical computations, in particular to the immediate decoupling of the evolution of the singlet and non-singlet components. We shall employ eq. (3.7) when solving the RGE for the evolution operator. For the latter, we shall use both eq. (2.25) and eq. (2.26), and distinguish the solutions thus obtained by referring to the former as ∆ 1 , and to the latter as ∆ 2 .

Non-singlet large-z solutions
As is customary, solving the RGE of the evolution operator is easier if one uses the nonsinglet (NS), singlet (S), and photon PDFs, since the non-singlet decouples and thus one works in a one-dimensional flavour space in that case, which leads to scalar (as opposed to matrix) equations. The reader can find all of the relevant definitions for the evolution in terms of NS, S, and γ in sect. 3 of ref. [13].

∆ 1 result
According to the definition given at the end of sect. 3, the ∆ 1 solution is that obtained by solving eq. (2.25) in the ∆ scheme. For the NS case, the only relevant scheme-changing function is that of eq. (3.5). A direct computation leads to the following Mellin-space result: K with: In eq. (4.1) we have introduced the symbol ∞ = in order to understand that terms which are subleading 6 in N in the N → ∞ limit are discarded. The same symbol will be employed in the case of the inverse Mellin transform, to discard terms subleading when z → 1.
Equation eq. (4.1) is then employed in the solution of eq. (2.25), for which we adopt a notation more in keeping with the scalar nature of the NS case, by means of the following formal replacements: By using eq. (2.5), it is straightforward to see that eq. (2.25) in the NS case has the following solution (note that with eq. (4.1) the inverse operator of eq. (2.22) exists not only in the asymptotic region N → ∞, which would be sufficient, but everywhere in the range of interest, N ∈ [1, ∞). This is convenient for the matching of the analytical large-z solution with its numerical counterpart, relevant to smaller z's): where the first term on the r.h.s. is the MS solution of eq. (II.5.50) with the ξ 1 andξ 1 parameters defined in eq. (II.5.51) and eq. (II.5.53): Here, b 0 and b 1 are the coefficients of the QED β function: and we have defined: With eqs. (2.16) and (3.2), and by bearing in mind that the term in the latter equation that multiplies the log(µ 2 0 /m 2 ) factor is the lowest-order Altarelli-Parisi kernel in the NS sector, we obtain the initial condition in the ∆ scheme in the Mellin space: where we have introduced the shorthand notation: 13) and the sought large-N solution for the NS PDF: (4.14) For consistency with eq. (4.1), we shall use in this section. We are interested in determining the z-space form of the NS PDF in eq. (4.14), which implies computing its inverse Mellin transform. In order to proceed, we re-express its rightmost term as a series in α, thus: With eq. (4.1), we see that the scheme-changing terms in eq. (4.16) are the following polynomials: where, after some algebra 7 : which is expressed in terms of the Pochhammer symbol: By using the algebraic identities above, eq. (4.16) can be recast as follows: where: By using the result of eq. (A.18), we can express the inverse Mellin transforms relevant to eq. (4.20) as follows: In eqs. (4.23) and (4.24) the index p labels the subleading-ness of the corresponding log(1 − z) term: the larger p, the more subleading the associated contribution. The idea is that of computing the inverse Mellin transforms of N −ξ 1 S 1 and N −ξ 1 S 2 by employing 7 It is understood that if the lower entry of a binomial coefficient is larger than the upper one that coefficient is set equal to zero.
eqs. (4.23) and (4.24) after having exchanged the order of the summations over j and p. Thus: with: where, from eqs. (4.23) and (4.24): We have defined 8 : 31) and the S i,p functions are such that: which show again the role of the index p in determining the subleading-ness of a given contribution. In practice, this implies keeping only the first few terms in the sums over p in eqs. (4.25) and (4.26), for which one can compute explicitly the S i,p functions, by carrying out the summations over j that enter their definitions. The results for the lowest p values are:Ŝ In view of a numerical implementation, we also have computed these contributions up to p = 5, which are too long to be reported here; the extension to yet larger p values poses no problem, but it is not particularly useful in view of eqs. (4.32) and (4.33).
In summary, the NLL-accurate large-z expression of the NS PDF in the ∆ scheme obtained by solving eq. (2.25) for the RGE of the evolution operator is: This implies that the strict asymptotic form of eq. (4.40), where all terms that vanish at z → 1 are ignored, reads as follows (by taking eq. (A.22) into account): The comparison of eq. (4.41) with the solution 9 of eq. (1.1) shows that all log(1 − z) terms that appear in the latter 10 (bar for that multiplied by L 0 ; note that usually one sets L 0 = 0) are features of the MS scheme, and disappear if the ∆ scheme is employed. We also point out that eq. (4.41) could have been directly obtained from eq. (4.14), by taking the strict N → ∞ limit of the scheme-dependent contribution to the evolution operator: and by employing eq. (A.18) with q = 0, 1.

∆ 2 result
In this section we shall again obtain the asymptotic form of the NS PDF in the ∆ scheme; at variance with the procedure of sect. 4.1, here this will be done by employing eq. (2.26) rather than eq. (2.25). This leads to the following solution for the evolution operator, which is the analogue of eq. (4.5): where we have used eqs. (4.1) and (4.6), and defined: Equation (4.43) suggests to introduce the following parameters: With the initial condition in the ∆ scheme, eq. (4.12), we therefore arrive at the analogue of eq. (4.14), namely: The computation of the inverse Mellin transform of eq. (4.48) proceeds similarly to what was done in sect. 4.1. One first expands the rightmost exponent in eq. (4.48), thus obtaining the analogue of eq. (4.16): where: In the computation of the inverse Mellin transform of eq. (4.50) one can exploit eqs. (4.23) and (4.24). In this way, we arrive at the analogues of eqs. (4.25) and (4.26), namely: The key feature of eqs. , after some algebra one arrives at the following results: where the functions T i;k (z) are: At this point, one can exploit the explicit results for the d p coefficients, given in terms of their generating functional (eq. (A.21)). This leads to: , (4.63) By replacing these results into eqs. (4.57) and (4.58), respectively, we obtain 12 : . By putting everything back together, we finally arrive at the analogue of eq. (4.40): .
As was the case for the sums over p in eq. (4.40), one can limit oneself to considering only the first few terms in the sum over k in eq. (4.67) to obtain a reliable numerical prediction.

Comments
The comparison between the ∆ 1 and ∆ 2 solutions given in eqs. (4.40) and (4.67) shows dramatic differences in the z → 1 region. While the former has a fairly smooth behaviour (on top of the integrable-singular prefactor, which is present also at the LL), the latter features a non-integrable singularity (induced by the exponential factor e ∆α log 2 (1−z) ) with alternating signs (induced by the function G d ). It is remarkable that these issues of the ∆ 2 solution stem entirely from the running of α; at fixed α, one has ∆α = 0, which eliminates the problem.
Ultimately, the difference between the evolution operators relevant to the ∆ 1 and ∆ 2 solutions, eqs. (4.5) and (4.47), would be innocuous, were it not for the fact that the schemechanging kernel, which is exponentiated, is driven by a soft-collinear double logarithm, which is the origin of the log 2N term that is dominant in the Mellin space, and that badly spoils the series expansion in α of the evolution equations at N → ∞; only by keeping all terms in such a series one is able to avoid this problem. Lest this fact be misunderstood: we are talking here about the logarithmic behaviour of the PDFs (which is what is relevant in physical applications). The fixed-order expansions of the PDFs (that is not relevant to physics, but useful in QED as a check of the correctness of the computations) give the expected results: namely, that eqs. (4.40) and (4.67) differ by O(α 2 ) terms (i.e. of NNLO), and that their O(αt 0 ) coefficient coincides with the initial condition (in the z → 1 limit) in the ∆ scheme.
We point out that the problem with the ∆ 2 solution cannot be studied numerically in QED in a direct manner. By setting α(µ) = 1/128 and α(µ 0 ) = 1/137, the factor e ∆α log 2 (1−z) is equal to one for z = 1 − 10 −x , with x 34, and no numerical integrator will generate values of z so close to one. However, the integration of e + e − cross sections requires in any case a change of variables that allows one to efficiently probe the z → 1 region (which, owing to the presence of the LL-like integrable-singular prefactor, is relevant to any kind of electron PDF), and this can only be done by knowing the analytical behaviour of the PDF in this region, hence the problem with the ∆ 2 solution.
The QCD case. The evolution operator in QCD obeys formally the same equations as the QED one, with the obvious formal exchange α → α S , and the crucial difference that the sign in front of the b 0 term on the second line on the r.h.s. of eq. (2.26) is negative rather than positive 13 . A simple algebra, that takes into account the fact that the evolution variables t are defined in different ways in QCD and QED, shows that in QCD one arrives at a solution which has the same functional form as that of sect. 4.2, with ∆α → ∆α S . Needless to say, this applies to heavy-quark perturbative fragmentation functions 14 , and not to PDFs, since this is the only case where the initial conditions are the same as their QED counterparts. On top of that, the z → 1 region is known to receive significant contributions from non-perturbative physics, which has no analogue in QED. Finally, and most crucially, while in QED ∆α > 0, in QCD ∆α S < 0, which implies that heavy-quark perturbative fragmentation functions are asymptotically strongly suppressed, rather than enhanced as for PDFs in QED. In conclusion, the QCD case poses no problems, and the only potentially unpleasant fact is that the strong suppression at z → 1 of the NLL solution has no LL counterpart. This has no phenomenological consequences.

Singlet-photon large-z solutions
In this section, we study the large-z behaviour of the singlet and photon NLL PDFs. Given the results of sect. 4, we shall limit ourselves to working in the ∆ 1 framework.
The procedure we shall adopt follows closely that of appendix B of ref. [13]. Because of this, we shall only sketch it, and refer as much as possible to the equations of that paper. While those equations complement the technical information given here, they are not necessary for a general understanding of the strategy used to arrive at the sought solutions. Such a strategy can be summarised as follows: 1. Write the evolution kernels in terms of a double series in α and 1/N . [29,30], which simplifies significantly owing to the kernels being diagonal at this order in N .

Find the leading-N solution in the form of a Magnus expansion
3. Treat the 1/N contributions as perturbations to the leading-N evolution.
4. The solution is finally achieved in a standard way except for the fact that it requires, in order for the inverse Mellin transforms to be computed analytically, a linearisation of the dependence of α on the variable t.
We start by recalling that by keeping only the divergent and constant terms of the Mellin transforms for N → ∞ one fails to obtain a good approximation of the actual z → 1 behaviour of the photon PDF. In MS, this is seen to stem from the fact that the electron-PDF LO initial condition is a Dirac delta that, combined with the relevant off-diagonal term in the Altarelli-Parisi kernel, gives a contribution of the same order as that due to the diagonal kernel element, times the photon-PDF initial condition. This argument is independent of the factorisation scheme, but one needs to verify whether in a different scheme, such as ∆, contributions not present in an MS computation will significantly modify its conclusions. In view of this, the factorisation-scheme matrix that we shall employ will have the following form: which is eq. (3.7) written in the two-dimensional singlet-photon sector, and where the functions K (∆) S,N and K (∆) γ,N are the Mellin transforms of eqs. (3.5) and (3.6), computed in the large-N limit by including the leading terms suppressed by powers of 1/N , namely: Clearly, eq. (5.2) extends eq. (4.1) by adding to it the relevant subleading terms. As was the case for the non-singlet scheme-changing kernel, the matrix of eq. (5.1) is such that the inverse operator of eq. (2.22) exists everywhere. The evolution equation for the evolution operator, eq. (2.25), is re-written by introducing an evolution kernel as was done in eq. (II.B.8), namely: with: and where the Altarelli-Parisi kernels have been written as in eq. (II.B.1) and eq. (II.B.2): By construction, in MS (i.e. with K = 0), M N (t) coincides with that of eq. (II.B.7). The forms in eqs. (5.6) and (5.7) stem from the observation that, as far as the schemechanging terms are concerned, an expansion in α must not be carried out. Conversely, the contributions directly related to the Altarelli-Parisi kernels can be safely expanded.
In keeping with eq. (5.9), we introduce the leading and the first subleading terms of the kernels defined in eqs. (5.6) and (5.7): The pure leading-N evolution will therefore be controlled solely by the kernel: while both eq. (5.12) and the following kernel: will be relevant to computing the subleading corrections to the evolution. With a direct calculation we obtain: with eq. (5.15) coinciding with eq. (II.B.12). Note that eq. (5.15) is non trivial, because it is equivalent to saying that, at this order in N , the left and right multiplications by the scheme-changing matrices in eq. (5.7) do not give any extra contributions w.r.t. the MS result (indeed, this conclusion is not valid starting from the first subleading order in N ). Furthermore, owing to the presence of inverse matrices in the evolution kernels, one may suspect that terms suppressed by higher powers of 1/N in the Altarelli-Parisi and scheme-changing kernels (which have been neglected here) could result, at the level of the evolution kernels, in contributions of the same orders as those kept (i.e., 1/N ); we have explicitly verified that this is not the case. Following ref. [13], the solution for the leading-N evolution operator, which obeys: can be represented by means of the Magnus expansion [29,30] (see eq. (II.B.8)). Since both of the matrices in eqs. (5.14) and (5.15) are diagonal, so is their sum, and therefore only the first Ω term of the Magnus series is non-zero. Furthermore, the exponent of a diagonal matrix is the matrix whose elements are the exponents of the elements of that matrix. This implies that we can write the sought solution as follows: where: i.e. eq. (4.40). As far as the latter is concerned, eq. (5.20) allows one to make use of eq. (II.5.72), that is: Therefore, at this order in N the photon PDFs is: From eqs. (3.3) and (3.6) we thus obtain: which is the analogue of eq. (II.5.73). As is the case for the latter equation, eq. (5.23) tends to a constant at z → 1; but, at variance with what happens in MS, the leading-N result of the photon PDF in the ∆ 1 scheme is identically equal to zero when µ 0 = m. We now turn to considering the first subleading contribution, that stems from the second terms on the r.h.s. of eqs. (5.10) and (5.11), and is treated as a perturbation to the leading-N solution. By proceeding as was done in ref. [13], the evolution operator that includes both the leading and the first subleading contributions is, according to eq. (II.B.30), where, from eq. (II.B.19) and eq. (II.B.21): Here, the linearised version of the leading-N evolution operator has been employed, which stems from eq. (5.17) and reads as follows: By means of an explicit computation one obtains: Note that the leftmost matrix on the r.h.s. of eq. (5.26) is not affected by the linearisation procedure (it is already linear in α(µ)). The matrices M where: From eq. (5.2) we see that the (1, 1) element in eq. (5.34) is suppressed by an extra power of logN w.r.t. to (2, 1) element; thus, it will be discarded henceforth. Given the matrix elements above, one can compute the evolution operator of eq. (5.24). As in the MS case, we are interested here only in the photon PDF (since for the singlet the leading-N solution stemming from eqs. (5.18) and (5.19) is sufficient to give a good description of the actual large-z behaviour), and therefore we shall only consider eq. (II.B.31): The initial conditions are still parametrised as in eq. (II.B.32) and eq. (II.B.33), that is: but in the ∆ scheme we have (eq. (4.12)): γ,j (z) .
We recall that the form of eq. (5.46) is suited to taking the fixed-α limit. In general, the prefactor of that equation is equal to α(µ 0 )/α(µ) (see eq. (II.B.79)), whence eq. (5.47) times that prefactor coincides with eq. (5.23), as it should. As far as the j ≤ 4 contributions to eq. (5.46) are concerned, they all stem from the first term on the r.h.s. of eq. (5.40). They are formally similar to their MS counterparts, eq. (II.B.38), but are significantly different from them in practice. They read as follows: Thus, the largest power of logN in both the numerator and the denominator of eq. (5.48) is equal to 3 (while it is equal to 4 and to 1, respectively, in MS). This is true when α = 0. In the α → 0 limit, i.e. at the LL, the x (j) i and y (j) i coefficients are such that the numerator of eq. (5.48) is N -independent, while the denominator is of O(logN ). Explicitly, one finds that at α → 0 eq. (5.46) coincides with eq. (II.B.88), i.e. with the MS result; this is indeed what must happen, since at the LL any scheme dependence disappears. At the NLL, the behaviour of the numerator and denominator of eq. (5.48) implies that the photon PDF is at most a constant at z → 1. We can easily compute such a constant; terms that vanish as some inverse power of log(1 − z) would require a lengthy calculation (much more so than in the case of MS), and this does not appear to be worth the effort. We obtain what follows: Eq. (5.50) is obtained from eq. (5.49) by undoing, at the purely parametric level, the linearisation introduced to deal with the subleading-N evolution operator; this stems from eq. (5.29). While the undoing of the linearisation is not justified mathematically, it seems to be fairly harmless, and if so eq. (5.50) should be considered as our best prediction for the large-z photon PDF.
The first term on the r.h.s. of eq. (5.50) stems from the leading-N contribution, whereas the second and third terms stem from the subleading contribution. In the strict z → 1 limit, where all terms that vanish need to be ignored, the first and the second terms are equal to non-zero constants of the same order, while the third term vanishes. Thus, by neglecting higher-order terms (in α) originating from the expansion of ξ 1 , we have: This has to be compared with its analogue in MS, eq. (1.2), which is divergent at z → 1, and is so that such a divergence is entirely driven by the subleading-N solution of the evolution equation. While this does not happen in the ∆ scheme, also in this case it is true that subleading-N contributions cannot be neglected for the computation of the photon PDF to be reliable in the z → 1 region.
Finally, note that eq. (5.50) vanishes identically when L 0 = 0. This is a consequence of the fact that, in the ∆ scheme, when L 0 = 0 all initial conditions are equal to zero, except for the O(α 0 ) one of the electron (which is equal to δ(1 − z)). In turn, this implies that the largest power of logN in the numerator of eq. (5.48) is equal to 2, while that in the denominator is equal to 3. Therefore, at L 0 = 0 the photon PDF is in general non-zero for a generic z, but vanishes at least as fast as log −1 (1 − z) at z → 1.

Conclusions
In this paper, we have studied analytically the behaviour of the NLL-accurate electron and photon PDFs of the electron in the z → 1 region, by considering a DIS-like factorisation scheme, called ∆. This region gives the dominant contribution to e + e − cross sections, and is therefore important, in view of the precision targets for theoretical calculations to be achieved in the context of the physics programs of future lepton colliders, that the systematics associated with the factorisation-scheme dependence be carefully assessed. This can be done by comparing predictions for any given process obtained with PDFs evolved with different factorisation schemes, such as those derived here and their MS analogues, previously computed at the NLL accuracy in ref. [13].
The main results of this work are given in eqs. (4.40) and (4.67) (for both the singlet and the non-singlet), and eq. (5.50) (for the photon). We have shown that, while the definition of the ∆ scheme aims to obtain NLL-accurate PDFs which are as close as possible to their LL counterparts, this happens only if suitably-defined initial conditions are evolved by including the effects of the running of α to all orders. A truncation of running-α effects to some order in α in the evolution equation causes the electron PDF to develop a nonintegrable singularity at z = 1. This is a result which is peculiar to QED, and its QCD analogue (the fragmentation function of a heavy quark) has no comparable issue.
When the running of α is properly taken into account, the NLL PDFs in the ∆ scheme have the same z → 1 functional forms as their LL analogues. This is not the case for the MS scheme, whose PDFs feature additional log(1 − z) terms w.r.t. those possibly present at the LL. While this is a positive characteristics from a theoretical viewpoint, in view of the fact that PDFs are non-physical quantities its phenomenological implications can be assessed only after computing observable cross sections. Therefore, by assuming uniform convergence: for any positive ρ, then: The expression in eq. (A.6) is well suited to the computation of the derivatives of interest. In fact, one can apply the Faà di Bruno formula for composite derivatives, namely: valid for any two functions f and g. In our case, from eq. (A.6) we have: In eq. (A.7), by B q,i we have denoted the incomplete exponential Bell polynomials. In fact, given eq. (A.8), all of the derivatives of the function f in eq. (A.7) are identical to the function itself, and can therefore be factored out from the summation. In this way, the complete exponential Bell polynomials: emerge naturally in eq. (A.7). By putting all this together, we therefore obtain: where we have introduced the q-dimensional vector Φ q , whose components, according to eq. (A.7), are: With the explicit form of eq. (A.9), this leads to: Ψ q (ξ) = γ E + ψ 0 (ξ), −ψ 1 (ξ), . . . , (−) q−1 ψ q−1 (ξ) , (A.14) where ψ k (ξ) are the polygamma functions: We can now exploit the following sum rule of the complete exponential Bell polynomials: 16) and the identity: This leads to an explicit expression for the inverse Mellin transform of eq. (A.1), namely 15 : The derivation above, and the role that the Bell polynomials play in the Faà di Bruno formula, suggest an alternative form for the coefficients d k . Indeed, one can verify by employing that formula that the quantity: can be used as the generating functional of the d k coefficients, namely: The results of eqs. (A. 19) and (A.21) coincide. We mention the fact that the functions A(ξ) and B(ξ), introduced in ref. [13] and that appear in eq. (1.1), are such that: The form of eq. (A.19) allows one to find relationships among the d k coefficients. Explicitly, by using the following recursion relation for the complete Bell polynomials: we obtain: (A.25) 15 For results of similar kind as that of eq. (A.18), see e.g. refs. [31][32][33].
In fact, eq. (A.25) can be used to compute the d k coefficients with increasing values of k, starting from k = 0 and using the fact that: Furthermore, the expression of the d k (ξ) coefficients in terms of the Bell polynomials allows one to compute the ξ → 0 behaviour of such coefficients, which is directly relevant to the electron PDFs (where ξ = ξ 1 , and ξ 1 0.05 for scales of O(100 GeV)). This is done as follows. Consider: where the rightmost equality follows from the fact that the complete exponential Bell polynomials satisfy the generating-functional relation: (A.28) By construction: We can now use the following representation of the polygamma functions: to write their small-ξ behaviour as follows: The series on the r.h.s. of this equation can easily be summed, leading to: We can use this expression for evaluating the arguments of the Bell polynomials in eq. (A.27); according to eq. (A.14): ∞ (ξ) = γ E δ j1 + (−) j−1 ψ j−1 (ξ) . (A.33) With this, we obtain: Therefore: In the computation of the derivatives of interest, eq. (A.29), it is a simple matter of algebra to see that: This proves that all d k coefficients diverge at most as 1/ξ at small ξ, which is a remarkable fact in view of eq. (A.31), that exhibits a much steeper ξ → 0 behaviour. By using the Faà di Bruno formula in eq. (A.36) we arrive at: Before concluding this section, we remark that eq. (A.18) can be used to arrive at a minor improvement of the result for the MS photon PDF w.r.t. that obtained in ref. [13]. This is due to applying eq. (A.18) in eq. (II.B.39) and eq. (II.B.40) 16 , which leads to the possibility of computing the functions in eq. (II.B.41), eq. (II.B.43), eq. (II.B.44), and eq. (II.B.45) exactly in κ, rather than by means of an O(κ 2 ) expansion as in ref. [13]. We obtain: with the d i coefficients given in eq. (A.19) or eq. (A.21). In practice, the parameter κ actually used in the computation of the MS photon PDF is either very small or equal to zero; thus, the numerical results of eqs. (A.39) and (A.40) do not differ either significantly or at all from those of ref. [13]. Note that eq. (II.B.41) corresponds to limiting the summation on the r.h.s. of eq. (A.39) to i ≤ 2. Conversely, the infinite sum of eq. (A.39) can be formally expressed by means of a Borel integral, thanks to eq. (A.21): .