Subleading Logarithmic QED Initial State Corrections to $e^+e^- \rightarrow \gamma^*/{Z^{0}}^*$ to $O(\alpha^6 L^5)$

Using the method of massive operator matrix elements, we calculate the subleading QED initial state radiative corrections to the process $e^+e^- \rightarrow \gamma^*/Z^*$ for the first three logarithmic contributions from $O(\alpha^3 L^3), O(\alpha^3 L^2), O(\alpha^3 L)$ to $O(\alpha^5 L^5), O(\alpha^5 L^4), O(\alpha^5 L^3)$ and compare their effects to the leading contribution $O(\alpha^6 L^6)$ and one more subleading term $O(\alpha^6 L^5)$. The calculation is performed in the limit of large center of mass energies squared $s \gg m_e^2$. These terms supplement the known corrections to $O(\alpha^2)$, which were completed recently. Given the high precision at future colliders operating at very large luminosity, these corrections are important for concise theoretical predictions. The present calculation needs the calculation of one more two--loop massive operator matrix element in QED. The radiators are obtained as solutions of the associated Callen--Symanzik equations in the massive case. The radiators can be expressed in terms of harmonic polylogarithms to weight {\sf w = 6} of argument $z$ and $(1-z)$ and in Mellin $N$ space by generalized harmonic sums. Numerical results are presented on the position of the $Z$ peak and corrections to the $Z$ width, $\Gamma_Z$. The corrections calculated result into a final theoretical accuracy for $\delta M_Z$ and $\delta \Gamma_Z$ which is estimated to be of O(30 keV) at an anticipated systematic accuracy at the FCC\_ee of \sim 100 keV. This precision cannot be reached, however, by including only the corrections up to $O(\alpha^3)$.


Introduction
At the planned future e + e − facilities which operate at high energy and at large luminosity, like the ILC, CLIC [1][2][3][4], the FCC ee [5], and also muon colliders [6], tests of the Standard Model are possible at unprecedented accuracy. This concerns a further detailed exploration of the Z peak, beyond what was possible at LEP [7], detailed studies of Higgs boson production using ZH final states, accurate scans of the top-quark threshold region, and various other precision measurements. Due to this the accuracy of the masses and widths of the heaviest particles of the Standard Model will be significantly improved.
One important ingredient to these experimental precision studies are the QED radiative corrections and in particular those due to initial state radiation (ISR). Very recently the O(α 2 ) ISR corrections for the process e + e − → γ * /Z * have been completed in a direct calculation [8][9][10].
Here α denotes the fine structure constant. It has been shown that the result for all channels of a previous calculation [11] needed to be corrected in the non-logarithmic terms at O(α 2 ). Agreement has been found with the results of Ref. [12]. Due to this it has also been proven that one may use the method of massive operator matrix elements (OMEs) for these calculations and that the Drell-Yan process factorizes for massive fermionic states.
In the present paper we take advantage of this method and extend the calculation to the first three logarithmic terms up to the order O(α 5 ) beyond the complete O(α 2 ) corrections and thus reach O(α 5 L 3 ). For comparison we also calculate the leading order contributions of O(α 6 L 6 ) and one more subleading term O(α 6 L 5 ), with L = ln(s/m 2 e ), s denotes the center of mass energy squared of the annihilation process and m e is the electron mass. The universal corrections O((αL) k ) are known in analytic form to order k = 5, cf. [13][14][15][16][17][18][19][20], accounting for the non-singlet and singlet contributions in the unpolarized and polarized case. The method used in Ref. [12] can be extended to higher orders in the coupling constant. For the first subleading term at O(α 3 ) 1 , all ingredients forming the radiators in terms of Mellin transforms are known from the calculation of the anomalous dimensions in QCD [22], the Wilson coefficients of the massless Drell-Yan process [8,23,24] and the massive OMEs in [12]. To obtain the O(α 3 L) and O(α 3 L) correction we need also to calculate the massive OME Γ (1) γe . This also applies to the higher order subleading contributions. This series could be continued straightforwardly to higher and higher order, for the first three terms at each order in the logarithmic expansion in Mellin N space on the expense of longer and longer expressions.
It turns out that it is indeed the case that one has to reach at least corrections of the order O(α 5 L 4 ) for the ISR corrections to satisfy the ambitions goals at the FCC ee of ∼ 100 keV both for the Z mass, M Z , and the Z width, Γ Z , on the theoretical side. The method of massive operator matrix elements [12] makes this calculation possible. In the present approach the constant term of O(α 3 ) and its higher order logarithmic extensions are still missing. They require still higher order massive OMEs and also massless Wilson coefficients in analytic form. Their size is, however, gradually smaller.
The paper is organized as follows. In Section 2 we present the structure of the QED ISR radiative corrections to the Born cross section following from the renormalization group equations (RGEs), for which we present the analytic solutions in terms of anomalous dimensions, massive OMEs and massless Wilson coefficients for all contributions calculated in the present paper. In Section 3 we calculate the massive OME Γ (1) γe . The radiator functions in z space for the O(α 3 ) and O(α 4 ) corrections are presented in Section 4. The radiators can be expressed in terms of harmonic polylogarithms [25], to which the Nielsen integrals and classical polylogarithms form a subset [26,27]. In Section 5, we present numerical results on the corrections in the kinematic region around the Z peak and we determine the corresponding corrections to the pole mass of the Z boson and the Z boson width, Γ Z . Section 6 contains the conclusions. In Mellin N space the radiators are given in terms of harmonic and generalized harmonic sums. If compared to the z space representation they are more compact. Due to the appearance of generalized harmonic sums we derive in Appendix A the singularity structure of the radiators in the complex N -plane and in Appendix B we present all radiators calculated in the present paper in Mellin space. 2 The Initial State Corrections to the e + e − Annihilation Cross Section The initial state QED radiative corrections to e + e − → γ * /Z * can be calculated by applying the method of massive operator matrix elements, cf. [12]. It has been demonstrated by the recent complete calculation in [8][9][10] that the effective method of Ref. [12] is delivering the complete result. Particular non-logarithmic contributions with vanishing massive OME to the Drell-Yan process in the constant term O(α 2 ) could be structurally absorbed in the relations and appear as contributions to the massless Wilson coefficients given in [12], Eqs. (40)(41)(42). It is due to this agreement, that one can now safely apply this method also for subleading logarithms to higher orders in the coupling constant by solving the associated renormalization group equations in the massive case. The method has been used in Ref. [11] before for the logarithmic enhanced contributions to O(α 2 L). 2 The massive effects due to the finite electron mass m e enter here through process-independent massive operator matrix elements. The Wilson coefficients are those of the massless Drell-Yan process [8,23,24]. Note that there are differences between the corrections to the vector and axial-vector coupling.
In z-space the different contributions to the radiators are connected by Mellin convolutions ⊗ which are defined by dz(z N −1 − 1)f (z) (2) for regular functions and +-distributions, respectively. The most recently calculated quantities are the massive OMEs given in [12]. The radiative corrections to the differential scattering cross section is given by in z space. Here we defined L = L + ln(z), and σ (0) denotes the Born cross section, cf. [12], Eq. (8), and Ref. [28]. a 0 = a(µ 2 = m 2 e ) is the normalized fine structure constant with a(µ 2 ) = α(µ 2 )/(4π), which we will widely use in the following. Here it is convenient to refer to a 0 only and account also for all evolution contributions by the functions c k,l .
For the Born cross section for e + e − annihilation, σ 0 , we will consider s-channel e + e − annihilation into a virtual gauge boson (γ, Z) which decays into a fermion pair f f and e = f . This process both describes lower energy γ-exchange and the Z-resonance.
see e.g. [28,29] 3 . Here the final state fermions are considered to be no electrons, to obtain an s-channel Born cross section. In Eqs. (5, 6) the electron mass is neglected kinematically. N C,f is the number of colors of the final state fermion, with N C,f = 1 for colorless fermions, and N C,f = 3 for quarks. The function G(s) = 1 in the case of the pure perturbative calculation. s is the cms energy, Ω is the spherical angle, θ the cms scattering angle, and the effective couplings The reduced Z-propagator is given by where M Z and Γ Z are the mass and the width of the Z boson and m f is the mass of the final state fermion. Q e,f are the electromagnetic charges of the electron (Q e = −1) and the final state fermion, respectively, and the electro-weak couplings v i and a i read v e = 1 sin θ w cos θ w I 3 w,e − 2Q e sin 2 θ w , a e = 1 sin θ w cos θ w I 3 w,e , where θ w is the weak mixing angle, and I 3 w,i = ±1/2 the third component of the weak isospin for up and down particles, respectively.
The inclusive s-channel annihilation scattering cross section σ e + e − is given by with s = sz, σ (0) (s ) the Born cross section and R(z, s/m 2 e ) the distribution-valued radiator [30] describing the initial state radiation of photons and light e + e − pairs. The lower bound s 0 = sz 0 is an invariant mass squared depending on the experiment, cf. [7]. In the later numerical illustrations we will choose s 0 = 4m 2 τ , with m τ the τ lepton mass. The general decomposition of the scattering cross section in Mellin space is given by, cf. [11] The terms in the brackets [...] are Mellin-convoluted. Only massive OMEs of the kind Γ e ± e ± and Γ γe ± contribute because the process considered has electron-positron initial states. The last term in Eq. (16) is only contributing with O(a 4 ). µ denotes the factorization and renormalization scale. As we will see later, it will cancel in the scattering cross section, when performing the expansion consistently to a certain order in a 0 . The massive OMEs, Γ ij , and massless Wilson coefficients,σ ij , obey the following series representations with the logarithms Λ and λ given by The massive OMEs Γ ij and massless Wilson coefficientsσ kl fulfill the following renormalization group equations, cf. [33], and the QED β function has the representation Eq. (23) gives an overview on the orders of the expansion of the radiators in the fine structure constant which are now available, including the results of the present calculation, The expansion coefficients c k,l of Eq. (3) up to the sixth order in a 0 in Mellin space are given by Eqs. (25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40)(41). They are expressed by the anomalous dimensions γ (k) ij [22], the expansion coefficients of the massless Drell-Yan cross section [8,23,24],σ (k) ij , the expansion coefficients of the QED β function and the renormalized massive OMEs Γ (k) ij , where (k + 1) denotes the loop order. In the following we use the notation When working in z space we will use P ij . One obtains c 5,4 = 13 12 where P (0) γγ = −2β 0 . These functions are derived by solving the renormalization group equations for the massive OMEs and the massless Wilson coefficients, (20,21) up to six-loop order. We show also the constant term c 3,0 which contains the term Γ (2) ee , not having been calculated yet. 5 A few remarks are in order. The sub-system cross section for the Drell-Yan process is flavor dependent, cf. Eq. (16), which results in the present case from the to e + e − initial state, the vertex to which the produced neutral current gauge bosons, γ * or Z * , couple. I.e. e.g. in the term describing intermediate photon exchange in Eq. (28) by P γe has to be read as either P and the energy-momentum sum-rule is obeyed. Similar relations hold in higher order and also apply to the corresponding massive OMEs. P (0) ee is flavor conserving, i.e. it does either describe an e − → e − or an e + → e + transition. Furthermore, one has where P (0) γe − (z), like in Quantum Chromodynamics (QCD) setting the color factor C F = 1, cf. [35]. This complication is usually not present in QCD, since there the splitting functions act in the singlet case on fully symmetrized flavor distributions, like the flavor singlet distribution Σ(x, Q 2 ), summing over all quark and antiquark flavors.
The running coupling constant is the solution of the differential equation with β k the expansion coefficients of the β-function, with [ [45][46][47][48] in the single fermion approach, N F = 1, retaining only electrons. Since we are dealing with the first three logarithmic orders from O(α 3 ) onward only terms up to β 2 are contributing. The solution of (45) has been derived in Ref. [49], Eqs. (2,3), in the MS scheme by keeping all terms up to β 2 in closed form. The corresponding perturbative solution for a(µ 2 ) = a(a 0 , L) is then given by One verifies the correctness of (47) by inserting it into Eq. (45). Note that if one does not expand the fine structure constant a(µ 2 ) w.r.t. its reference value a 0 = a(µ 2 = m 2 e ), the µ dependence in the RGE-decomposition given in Ref. [12] is not canceled. However, expanding to the respective order in the coupling constant intended, a scheme-invariant expression is obtained. This is the reason, why we are expanding the logarithmic dependence of the coupling constant and write the scattering cross section in terms of a 0 .
In the above expressions we presented the sub-system cross sectionsσ ij in a genuine way. One should note that these quantities are partly different for vector and axial-vector couplings, cf. [10,23] and so are some of the radiators. The first expressions for c i,j in z space are presented in Section 4. In Appendix B all coefficients c i,j in Mellin N space, which are more compact, are given.
Sumrules do not only hold for the splitting functions (43)(44) but also for the universal unrenormalized massive operator matrix elements, cf. [50,51] to two-loop order which we have verified for N F = 1.

The Operator Matrix Element Γ (1)
γe For the calculation of the operator matrix element Γ (1) γe we follow the notation of [12]. After wave function and mass renormalization we can write the operator matrix element aŝ where Z CT are the counter term contributions due to mass and wave function renormalization, A (i) γe are the unrenormalized operator matrix elements at i-loops andâ is the unrenormalized fine structure constant.
The renormalized operator matrix elements in the MOM-scheme are given by, cf. [51], where the unrenormalized coupling constant can be expressed viâ with and the spherical factor The operator matrix element, after charge and wave function renormalization, can therefore be written aŝ The predicted pole structure serves as a test on the calculation. The renormalized OME in the MS-scheme is given by where we used the relation The Feynman diagrams contributing to Γ (1) γe are shown in Figure 1, with the corresponding symmetrization understood. They can be represented in terms of 18 master integrals by performing the integration-by-parts reduction using the package Litered [52]. Here, as in previous calculations, cf. [53], we resummed the local operators into linear propagators. The master integrals are either calculated using conventional methods like generalized hypergeometric functions, cf. e.g. [54], or can be obtained by solving ordinary differential equations, cf. e.g. [55]. 6 The code Tarcer [56] has been used for checks and to determine initial values. γe . For the notation and the Feynman rules see Ref. [12]. The symbol denotes the counter-term insertion.
The expansion coefficients of the unrenormalized OMEÂ (1) γe are given bŷ see [12]. The coefficient Γ The corresponding result in Mellin N space is the obtained as the N th expansion coefficient in the auxiliary variablex and the z space representation can be obtained by a subsequent inverse Mellin transform. The renormalized OME A (2) γe is given by with the polynomials A (2) γe is expressed by harmonic sums [36,37] and generalized harmonic sums [38,57] has its rightmost pole at N = 1 and is otherwise regular. In particular one may show that applying the algorithms of package HarmonicSums [36][37][38][39][40][41][42][43][44]. In z space the OME is given by with Here we refer to the harmonic polylogarithms [25], and define There are some terms in Eq. (73) which are ∝ 1/z 4 . This is a reflection of terms ∝ 1/(N − 4) in Eq. (59). One derives the small z expansion of A (2) γe (z) which is given by showing that the most singular terms are of O(1/z).

The Radiators in z Space
We express the expansion coefficients (25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37) in terms of harmonic polylogarithms [25]. The corresponding expressions can be obtained using the package HarmonicSums [36][37][38][39][40][41][42][43][44], after reducing the algebraic relations [58]. Because of the occurrence of some denominators 1/(1 ± z) l , l ∈ N, l > 1, one performs the corresponding series expansion and uses summation techniques encoded in the packages Sigma [59,60], EvaluateMultiSums and SumProduction [61]. As usual, one has to separate the radiators into the part ∝ δ(1−z), the contribution due to +-distributions and the regular part, For the inclusive cross section the integral (15) has to account for the fact that the radiator R(z, L) is distribution-valued [30]. For the inverse Mellin transform we use the notion PlusFunctionDefinition → 2 of the package HarmonicSums, We use the following notation The splitting function P (0) ee is thus given by We use the following conventions The expressions up to O(a 2 ) are given in Ref. [10]. In the following we distinguish the expansion coefficients in the vector and axial-vector case. Coefficients c ij without an index v or a apply to both cases. We also display the difference terms The δ-terms are The +-distributions are given by and the regular contributions read c ∆,reg In Figures 2 and 3 we compare the ratios of the three-loop terms with all contributions (25)(26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37) in the kinematic region of √ s ∈ [85, 95] GeV. The O(a) radiative corrections are large and amount to ∼ ±40%, followed by the O(a 2 ) corrections, still varying from +15% to −7%, cf.   Table 1: Shifts in the Z-mass and the width due to the different contributions to the ISR QED radiative corrections for a fixed width of Γ Z = 2.4952 GeV and s-dependent width using M Z = 91.1876 GeV [62] and s 0 = 4m 2 τ , cf. [7]. 8 Finally, we summarize the shifts of the Z peak and the corrections to the Z width, Γ Z , by the different orders of the ISR radiative corrections in Table 1. Here we compare the results for the fixed and the s-dependent width [63]. In Figures 2 and 3 we depicted only the corrections up to O(a 3 ). The other corrections up to O(a 6 L 5 ) are only illustrated w.r.t. its shift of the Z mass and change of the Z width since their behaviour is rather flat, yet they have an impact given the projected experimental resolutions. The difference to the results in Ref. [11] amounts to 4 MeV in Γ Z , [9,10].
The relative shifts in adding the respective order can be positive or negative. at leading order O(α k L k ) the level of 100 keV [64] 9 is only undershoot at O(α 5 L 4 ). Even the O(α 4 L 2 ) corrections are of the order of 50 keV, while at O(α 5 L 3 ) the level of 10 keV is reached. A missing link is still the O(α 3 ) term, which can be estimated to be roughly of O(30 keV), setting the frame of accuracy which is currently reached for the initial state corrections. The QED corrections to 3-loop order are still somewhat larger than the experimental resolution at the the FCC ee [5] making the inclusion of also higher order subleading corrections necessary.
Furthermore, we remark that we have calculated the inclusive ISR corrections only, assuming that the experimental data are extrapolated to the full phase space and only a cut in s is considered. The experimental requirements may be more ambitious, requiring more differential radiative corrections in the future. Due to different cuts, the corrections will turn out to be different and one has to carefully study all the cut dependencies. For a recent summary see [65]. From the size of the corrections it seems that 3-to 4-loop corrections have there to be provided too.
Numerical implementations for harmonic polylogarithms in Fortran needed for the radiators in z space are given in [66] and [67], respectively.

Conclusions
We have calculated the QED initial state corrections to the annihilation process e + e − → γ * /Z * up to the terms O(α 6 L 5 ). They come next to the recently completed O(α 2 ) corrections [8][9][10] and the well-known universal corrections O((αL) k ), which were known to fifth order in explicit form in the non-singlet and singlet case [13][14][15][16][17][18][19][20]. Here we included the first three logarithmic terms for the orders O(α 3 ) to O(α 5 ). The radiators are given by convolutions of splitting functions, the contributions to the Wilson coefficient of the massless Drell-Yan process and massive operator matrix elements. For the present corrections the massive OME Γ (1) γe had to be calculated newly. The other massive OMEs were given in [12] before. The corrections calculated in the present paper can be expressed in terms of harmonic polylogarithms, if one also allows for the argument (1 − x) in case of harmonic polylogarithms containing an index i = −1. In Mellin space the radiators can be represented in terms of harmonic sums and generalized harmonic sums. The present corrections may still miss terms of O(30 keV) for both δM Z and δΓ Z , which can be further improved by calculating the yet missing terms.
It is needless to say that by performing Mellin convolutions, using the quantities calculated in the present paper and the massless quantities available in the literature, one is now in the position to calculate all corrections of O(α k L k ), O(α k L k−1 ) and O(α k L k−2 ) for k ≥ 5 straightforwardly. For the projected experimental accuracies in inclusive measurements at the FCC ee they may not be needed beyond the orders already obtained.

A The Singularities of the Radiators in N space
For the use of complex contour integrals to calculate the radiative corrections the position of the singularities of the radiators in Mellin N space in the complex plane have to be known.
In the case of harmonic sums, except for S 1 (N ), it is known from their representation in terms of factorial series [68] that their singularities are given by the set {−n | n ∈ N}. The sum S 1 (N ) has a representation by the di-gamma function ψ(N ) for which the same holds.
A factorial series is given by The question to be answered is, which functions can be expanded into factorial series. The structure of (106) then provides the singularity structure for x ∈ C. If a function F (x) has the representation and the function ϕ(t) can be expanded into a Taylor series in (1 − t), partial integration will then lead to a factorial series. The radiators are expressed in terms of the following monomials and ξ ∈ {z, 1 − z}, with a i ∈ {0, 1, −1}. To perform the Mellin transform of (108) suitable regularizations have to be chosen. Since the Mellin transform of the complete radiator exists and is unique, these partly arbitrary regularization terms add up to zero. Now one has to assure that H a (ξ) can be expanded into a Taylor series in the variable (1 − z). This will require in case to change the argument from which is a valid operation on H a (ξ) on the expense of introducing the letter The structure in (108) can be generated by applying the following differential operator Furthermore, one has In this way and by the argument mapping (109) one arrives at valid functions ϕ(t) allowing to expand into a factorial series. The above construction has now to be applied to all radiators and one finds the set of singularities to be a subset of {−n | n ∈ N ∪ {−1}}.

B The Radiators in N space
In the following we list all radiators which were calculated in the present paper in Mellin N space for the use in Mellin space programs. The analytic continuation of the respective harmonic sums can be performed as described in Refs. [69,70]. The package HarmonicSums allows to derive the asymptotic representation of these coefficients. Their recursion relations follow from the ones of the harmonic and generalized harmonic sums. As has been shown, there are no singularities right to N = 1, which allows to perform the contour integral to z space with the usual contour in the singlet-case, see e.g. [31], surrounding the singularities of the expression in the complex plane, cf. Appendix A.
The radiators R ij are related to the expansion coefficients from c 3,3 to c 6,5 , also labeling the difference term (85) between the axial-vector and vector contributions, by Alternating sums can be rewritten in terms of Mellin transforms such that the contribution due to ln(2) terms cancel, cf. [71]. We have also checked that the evanescent poles present in the above radiators up to 1/ (N − 4), cancel, leaving the rightmost pole 1/(N − 1). The analytic continuation is performed from the even integers. It can be obtained in the analyticity region for N → C for both the harmonic sums and generalized harmonic sums expressing them for large values of |N | by their asymptotic expansion and by using the recursion relations to map to finite values of N ∈ C, [70].