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 MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} 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-α 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) JHEP07(2021)180 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.

JHEP07(2021)180
In this paper we shall compute the NLL-accurate electron PDFs by adopting a factorisation scheme alternative to the MS employed in ref. [13]; in particular, we shall use a definition whose physical motivation is the same as the one that underpins the so-called DIS scheme in QCD [16]. At variance with ref. [13], where analytical results have been presented for both the z → 1 region (to all orders in α) and the small-and intermediate-z region (up to O(α 3 )), in this paper we shall limit ourselves to studying analytically the z → 1 regime since, as is shown in eqs. (1.1) and (1.2), this is where the dominant effects are expected to appear; a comprehensive phenomenological study, that includes the numerical implementation of the results presented here, will be given elsewhere [17].
This paper is organised as follows. In section 2 we introduce the operator that will be used to solve the PDF evolution equations, and do so in a generic factorisation scheme. The scheme that is actually employed is called ∆, and is defined in section 3. In section 4 we derive the asymptotic z → 1 form of the non-singlet PDF in the ∆ scheme by solving its evolution equation in two different ways, namely by retaining running-α effects to all orders (section 4.1), or by truncating them at the first non-trivial order in α (section 4.2). The comparison between these two results is discussed in section 4.3. In section 5 we apply the same strategy as in section 4.1 to the singlet-photon sector. Finally, in section 6 we present our conclusions. Some Mellin-transform results that are used throughout the paper are collected in appendix A.

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:

JHEP07(2021)180
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 JHEP07(2021)180 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)

JHEP07(2021)180
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 JHEP07(2021)180 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.

JHEP07(2021)180 3 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 section 2, the matrix elements of K whose column index is not equal to the identity of the particle of interest (the electron here) give in any case beyond-NLO effects. By setting them equal to zero as in eq. (3.7) we take into account the trivial role of the positron at this order, and 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]. 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 .

JHEP07(2021)180 4 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 section 3 of ref. [13].

∆ 1 result
According to the definition given at the end of section 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: 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)

JHEP07(2021)180
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:

JHEP07(2021)180
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 eqs. (4.23) and (4.24) after having exchanged the order of the summations over j and p. Thus:

JHEP07(2021)180
with: where, from eqs. (4.23) and (4.24): We have defined: 8 and the S i,p functions are such that: , (4.34) 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. 9 At the NLO, the z → 1 expressions of the electron and the non-singlet PDFs coincide. 10 In addition to those one would obtain by expanding the (1 − z) −1+ξ 1 prefactor. Except for checks on its perturbative behaviour [13], such an expansion must never be carried out.

∆ 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 section 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: (4.44) 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 section 4.1. One first expands the rightmost exponent in eq. (4.48), thus obtaining the analogue of eq. (4.16): From eq. (4.49) we arrive at the analogue of eq. (4.20): where:

JHEP07(2021)180
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 summands on the r.h.s. of eqs. (4.53) and (4.54) have remarkably compact expressions; by construction: 11 The key feature of eqs.
where the functions T i;k (z) are: (4.60)

JHEP07(2021)180
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 .  z) terms, which have been resummed in the procedure carried out above. 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 section 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 13 Likewise, the sign in front of the first terms on the r.h.s. of eq. (2.25) must be changed. However, to the best of our knowledge the analogue of the ∆1 solution has not been considered in QCD.
14 In addition to the space-to time-like exchange, final-state evolution involves matrix kernels that are the transposed of those relevant to the initial state. This does not affect this discussion.

JHEP07(2021)180
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 section 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: the large-N limit by including the leading terms suppressed by powers of 1/N , namely: ≡ K [1,0] S,N + 1 N K [1,1] S,N , 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: 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: eqs. (5.18) and (5.19) we see that the former is identical to the one derived in section 4.1, 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: z) . (5.21) 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),

JHEP07(2021)180
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 we have obtained: 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 JHEP07(2021)180 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)): The j = 5 contribution to eq. (5.46) stems from the second term on the r.h.s. of eq. (5.40). From eq. (5.23) we have, analogously to eq. (II.B.73): 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 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 JHEP07(2021)180 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 JHEP07(2021)180 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.

JHEP07(2021)180
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: g(β)) B q,i g (1) (β), g (2) (β), . . . g (q−i+1) (β) , (A.7) 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: B q (x 1 + y 1 , . . . x q + y q ) =

JHEP07(2021)180
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): .