The $O(\alpha^2)$ Initial State QED Corrections to $e^+e^- \rightarrow \gamma^*/Z_0^*$

We calculate the complete $O(\alpha^2)$ initial state radiation corrections to $e^+ e^-$ annihilation into a neutral vector boson in a direct analytic computation without any approximation. The corrections are represented in terms of iterated incomplete (elliptic) integrals over alphabets of square--root valued letters. Performing the limit $s \gg m_e^2$, we find discrepancies with the earlier results of Ref.~\cite{Berends:1987ab} and confirm results obtained in Ref.~\cite{Blumlein:2011mi} where the effective method of massive operator matrix elements has been used, which works for all but the power corrections in $m^2_e/s$. In this way, we also confirm the validity of the factorization of massive partons in the Drell--Yan process. We add non--logarithmic terms at $O(\alpha^2)$ which have not been considered in previous calculations. The final results in the limit $s \gg m_e^2$ can be given in terms of Nielsen integrals.

1 Introduction e + e − colliders operating at high energy and at large luminosity measure the fundamental parameters of the Standard Model with high precision and perform crucial tests on the structure of the Standard Model.In the past the experiments at LEP obtained very precise results on the parameters of the Z-boson [3].The future large scale facilities like the ILC, CLIC [4][5][6][7], the FCC ee [8], and muon colliders [9] are planned to operate at even higher energies and luminosities.There one can perform in addition also precise scans of the t t-threshold measuring the properties of the top quark in detail and produce the Higgs boson under very clean conditions in ZH-final states, which will finally allow to understand more properties of the Higgs boson in great detail.
One important condition to perform these highly precise measurements is the exact knowledge of the QED radiative corrections for the process e + e − → γ * /Z * , which has to be known to twoloop order in the fine structure constant α, adding further logarithmic contributions in higher orders.A first calculation of the O(α 2 ) initial state radiative corrections to this process has been performed in Ref. [1].In this reference various approximations have been made in the integrands of the Feynman diagrams, to simplify the integration process.In 2011 it has been noticed, however, in a second calculation based on massive operator matrix elements (OMEs) [2] that the results deviated in all channels for the constant term at O(α 2 ), while the O(α) result and the logarithmic terms at O(α 2 ) agreed.In the latter calculation it was assumed that the Drell-Yan process with massive external lines factorizes.At that time, the new results did not yield a thorough counter argument against the results in [1], since one might have argued that there is no factorization in the massive Drell-Yan process.
There is actually only one way to decide which of the results is correct.One has to perform the complete calculation of the scattering cross section without any approximation or assumption analytically.In the final result one will of course expand in the ratio m 2 e /s ≈ 3 • 10 −11 , where m e denotes the electron mass and s the cms energy squared to obtain a compact result.
In Ref. [1] some processes which only contribute to the O(α 2 ) term and have no logarithmic contributions were not considered.A first calculation, however, in the massless case, was performed in [10] and later in [11] for the Drell-Yan process in Quantum Chromodynamics (QCD).Furthermore, differences appearing in the calculation of the contributing vector and axial-vector terms were not considered in [1].
In the present paper we perform a thorough analytic calculation of all contributing terms.We confirm the results given in [2] before and we add the pure O(α 2 ) terms, which cannot be derived using the method of massive OMEs [12].The calculation in Ref. [2] has been performed for vector couplings.Here we add also the axial-vector contributions, whenever they are not suppressed by power corrections of O(m 2 e /s).Our final results are expanded in m 2 e /s and we maintain all terms up to O((m 2 e /s) 0 ).The final radiators can be expressed by Nielsen integrals [13] S p,n (x) = (−1) n+p−1 (n − 1)!p!
The radiator functions have the general structure The respective differential cross sections are then given by dσ e + e − ds = 1 s σ e + e − (s )R z, α, s m 2 , ( with σ e + e − (s ) the scattering cross section without the initial state radiation (ISR) corrections and where s is the invariant mass of the produced (off-shell) γ/Z boson.Here and in the following the mass m denotes the electron mass m e , if not stated otherwise.
The paper is organized as follows.In Section 2 we present the Born cross section for the process.The O(α) corrections are given in Section 3. General aspects of the integration at two-loop order are discussed in Section 4. In Section 5 we present the results for the different processes contributing to the two-photon corrections.The non-singlet process of the e + e − pair radiation process is discussed in Section 6.The contribution due to the radiation of heavier final states in the non-singlet process is calculated in Section 7, followed by those due to the puresinglet process, Section 8, and the interference term between the non-singlet and the pure singlet terms, Section 9, in the case of e + e − emission.In Section 10, we give the results for processes that have no logarithmic contributions at O(α 2 ).The axial-vector contributions are discussed for the processes they contribute to.In all other cases the radiators are the same as in the vector case.Finally, we discuss the soft-photon exponentiation contributions beyond the radiative corrections to O(α 2 ) in Section 11, and Section 12 contains the conclusions.Numerical results on the Z peak, for the ZH production process and tt production have already been presented in Ref. [15].There and in [16] we also line out the numerical differences to [1].In Appendix A we present details on phase-space integrals which have been performed in the present paper.

The Process
We consider the process of e + e − annihilation into a virtual photon γ * or virtual Z * 0 boson above a mass threshold of s ≥ 4m 2 µ or larger, with m µ the muon mass and s the cms energy squared of the annihilation process.Also the production of other fermionic final states can be considered such as τ + τ − , massless q q and the corresponding heavy quark pairs.The phase space limit on s is s ≥ 4m 2 f .We will usually assume s ≥ 4m 2 µ or a more conservative cut.The differential Born cross section is given by where σ (0) (s ) denotes the integrated cross section of one of the above processes.It corresponds to the annihilation diagram in Figure 1.For s-channel e + e − annihilation into a virtual gauge boson (γ * , Z * ), which decays into a fermion pair f f , the scattering cross section reads The Born cross section for the process e + e − → Z * /γ * .
see e.g.[17,18]. 1 Here the final state fermions are considered not to be electrons, to obtain an s-channel Born cross section.In Eqs.(8, 9) the electron mass is neglected kinematically.α denotes the fine structure constant, 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) is set to 1 in the case of the pure perturbative calculation.s is the cms energy, Ω 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 with 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 are given by 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.
For the radiative corrections studied below, we will consider the integrated cross section (9) in the energy region of the Z-peak.
In the following we will use the fine structure constant with the normalization The scattering cross section including the contributions due to initial state radiation can be expressed as follows where R(z; a, L) is the distribution-valued [19] radiation function, with The different radiators calculated in the present paper sum to the following distribution (z, L) 3 The One-Loop Corrections At one-loop order only photonic corrections contribute.We will work in D = 4 dimensions.This allows to treat axial-vector couplings without an additional finite renormalization.
The O(α) e + e − annihilation graphs into a photon and a virtual gauge boson.
The O(α) virtual corrections to e + e − annihilation into virtual gauge boson.
The photon radiation and virtual diagrams are given in Figures 2 and 3, respectively.We parameterize the different contributions by where the indices S, V, and H stand for soft, virtual and hard.Let ∆ denote the energy cut-off for soft-photon bremsstrahlung.Then Here λ denotes a soft-photon mass, with λ m e .This parameter can be introduced in an Abelian gauge theory calculation and it cancels in the final result.
One obtains and is the Riemann ζ-function at integer argument.The virtual corrections are obtained from the one-loop massive Dirac form factor in the limit with [34][35][36]] The Pauli form factor is power suppressed in the limit s m 2 e .Finally, the hard contribution is given by Here there is also a lower bound on z from the minimal value of s ≥ 4m 2 f .One obtains for (22) the well-known result cf. [1], [2], Eqs.(40,96), where denotes the first order electron-electron splitting function.It can be promoted to a +-distribution where the +-prescription is defined by Then one can drop the θ-function and all contributions proportional to ln(ε).The above relations were derived for the vector case.The same corrections are obtained in the axial-vector case in the limit s m 2 e since the difference at one-loop order is suppressed by power corrections in m 2 e /s.
4 Analytic Integration of the O(α2 ) Corrections For the O(α 2 ) terms, non-trivial phase space integrals are occurring given by fourfold integrals.Details on their calculation are presented in Appendix A. They are given by two angular integrals and two further integrals over invariants.In course of these integrations one obtains square-root valued arguments in logarithms and polylogarithms, which are nested in part and have to be rationalized or transformed to single roots to perform the next integration. 2The principal way to obtain the corresponding square root-valued iterated integrals, containing real parameters in the letters, has been described already in Refs.[39,40].Square root-valued iterated integrals based on rational parameters have been considered earlier in [30].They occur as Mellin inversions of finite binomial and inverse binomial sums.In total, up to weight w = 3 iterated integrals emerge.We aim on an analytic iterated integral representation over an alphabet also containing square root-valued letters and will keep the complete dependence on Using special variables, it is also possible to expand in ρ 1 prior to the last integration is carried out.The integration has been performed using Mathematica.After this expansion one obtains a large number of logarithms and classical polylogarithms Li 2 (g i (z)), Li 3 (g i (z)) with involved, partly complex arguments.They have to be mapped to logarithms and polylogarithms of the convenient arguments z and 1 − z.For this we use associated differential equations.Some of the polylogarithms also depend on G-functions [30,[41][42][43], containing square root-valued letters.
Here the G-functions are defined by over an alphabet A .The different letters h({c i }, x) are not yet independent w.r.t. the associated differential field.Adding all contributions, the G-functions cancel.
We finally compare the exact analytic result, not expanded in ρ, with the expanded result including the O(ρ 0 ) terms, numerically.The expansion in ρ can also be performed starting with the complete result.This requires the introduction of suitable regularizations.We have performed the last step in the non-singlet case and obtained the same result as expanding below the last integral, using appropriate variables.In all cases the numerical comparison shows a relative agreement of O(10 −7 ) at s = M 2 Z , which is the expected result in this approximation.It shows that the formulae expanded in the light fermion mass can be used for experimental analyses.
The complete analytic results will be given in terms of iterative integrals over a certain alphabet of letters A, which are mostly square-root valued and contain real parameters.Some of them were occurring already in earlier investigations [30].They are labeled by v i .Other letters are new and are named by d i .The iterative integrals are given by The different letters of the alphabet, f k (t; z, ρ) ≡ f k are : We call the iterated integrals containing both Kummer-type letters [29] and the above letters, Kummer-elliptic integrals, since some of them integrate to elliptic structures, although they do so as indefinite integrals, 3 which are still iterative if compared to complete elliptic integrals and their extensions, cf.[46].We label the processes using the same scheme as in Ref. [1], i.e. process I: the photon emission case; process II: the non-singlet case for fermion-pair production; process III: the pure singlet process; process IV: the interference term between the non-singlet case for e + e − pair emission and the pure singlet case.Furthermore, we denote contributions not covered in [1] but belonging to the O(a 2 ) QED corrections as process B, in accordance with Ref. [10].
A few remarks on the size of the present calculation are in order.i) the size of the amplitudes amounts to 10 Gb (process I), 25 kb (process II), 56 kb (process III) and 124 kb (process IV).The calculation of process I required several months of code design and 30h of computation time.The reduction to the basis of iterative integrals took 1 day (process II), 1 month (process III), and 2 months (process IV).The integration time for processes II-IV amounted to minutes, 2h and 5h.No essential resources were necessary to perform the calculation for process B. The size of the project mainly resulted from the fact that only in the last step an expansion in the parameter ρ has been performed.The required computer power was not available at the time when Ref. [1] was worked out.This also applies to several computer-algebraic and mathematical calculation techniques we were able to apply, which became available only recently.Now we turn to the calculation of the individual sub-processes at two-loop order.

The Photonic Two-Loop Corrections
The photonic two-loop corrections consist out of the following six contributions with 1. T S 2 2 : both emitted photons are soft, Figure 4 2. T V 2 2 : both photons are virtual, Figure 5 3.
: one photon is soft, one is virtual, Figure 6 4. T S 1 H 1

2
: one photon is soft, one is hard, Figure 4 5.
: one photon is virtual, one is hard, Figure 6 6.T H 2 2 : both emitted photons are hard, Figure 4. We will calculate the different contributions in this order.
Due to the energy cut-off ε on the photon energy all contributions are functions.However, one can introduce +-distributions and drop the dependence on the photon cut-off.We still use the abbreviation besides of the δ(1 − z)-distribution.However, using the ε cut-off the +-distributions can be understood as simple functions.The double soft photon correction is obtained calculating the graphs given in Figure 4 in the soft photon approximation.Since the soft corrections are factorizing, one obtains cf. [1,47,48].Here the last term stems from the integral The diagrams for the double virtual corrections are shown in Figure 5, and is given by cf. (26) with cf. [35,36,49,50].Again, the Pauli Form Factor does not contribute in the limit s m 2 .The soft corrections given in [49] were corrected in [50].The virtual-soft corrections, Figure 6, are given by The virtual-hard corrections, Figure 6, read The O(a 2 ) virtual corrections to e + e − annihilation graphs into one photons and a virtual gauge boson.The external self-energy corrections are not shown.
The two hard photon corrections, Figure 4, yield All corrections but the virtual-hard corrections agree with the results in [1].
The complete photonic corrections are given by We now turn to the fermion-pair emission contributions in different channels.
6 O(α 2 ) Non-Singlet Corrections due to e + e − Emission The non-singlet contributions can be given by where R e + e − ,NS,S

2
and R e + e − ,NS,H 2 denote corrections due to soft, virtual and hard fermion-pair radiation respectively.The first two contributions were correctly given in [1] and read In the following we will calculate the hard contributions.They can be expressed by iterative integrals H a (u) ≡ H a up to weight w = 2 over the alphabet given in Section 4 and u = 4ρ/(1 − √ z) 2 .The contributing diagrams are shown in Figure 7.
The correction to the scattering cross section for electron pair emission is given by with and Figure 7: The e + e − annihilation graphs into a fermion pair and a gauge boson (process A).
One performs the transformation in order to introduce dimensionless quantities.This yields with Integrating this expression exactly, one obtains To obtain the expansion in the case of the emission of an e + e − pair one cannot simply set m e → 0 in (88) as has been done in Refs.[1,51]; see, however, Section 7. One first rewrites the integral (88) by Here f ρ (t, ρ, z) denotes the integrand f (t, ρ, z), expanded in ρ, including the ρ 0 term.For the first term in Eq. ( 92) the variable transformation is performed leading to t ∈ [0, 1].After that, the integrand can be expanded in ρ.A further variable transformation is necessary to rationalize the root One may choose The integral can now be performed.This leads to the correct result in the limit m 2 e s.Of course we can also expand Eq. (90) in ρ.We checked that both methods agree.Including the term ρ 0 one obtains This result differs from the one presented in Refs.[1,51] exactly by the term given in [16], Eq. ( 8) and agrees with the result obtained in Ref. [2] based on massive operator matrix elements.The reason for this disagreement lays in the neglection of some of the electron mass terms before all integrals have been performed.The full radiator therefore reads 7 Heavier Fermionic Final States in the Non-Singlet Process If in e + e − annihilation a heavier fermion pair f f with m f m e is radiated via a virtual photon from the initial state electrons one may use, cf.[51] and [1], Erratum, Here we consider the case of heavier charged lepton pairs and heavy quark (c, b, t) pairs.One obtains This result agrees with those of Refs.[1,51].One may derive it also using the method of massive operator matrix elements.Here, the external lines in Figure 7, [2], have to be taken massless, since m e m, and m denotes the mass of the internal fermion line.Γ (1),II ee in (41), [2], reads then [52] A MS,II ee,µ = and Eq. ( 75) in [2] has to be replaced by Eq. (4.16) in [52] for QED.Here the relation where N F denotes the number of massless fermions, which is N F = 1 here.The term T µ + µ − II corresponding to the one in (95) for µ + µ − pair radiation reads with L µ = ln(s/m 2 µ ).

The pure singlet corrections
The diagrams contributing to the pure singlet case are shown in Figure 8.Here one has to distinguish between the vector and axial-vector case since different corrections are obtained.The radiator in the vector case is given by and consists of the direct terms and the interference terms Figure 8: The e + e − annihilation graphs into a fermion pair and a gauge boson (process A).
In the axial-vector case one has to replace T III,v,int by T III,a,int in Eq. ( 103) The result (103) differs form that in [1] by the term in [16], Eq. ( 8).Note that the interference term between the diagrams in the upper line and the ones in the lower line of Figure 8 appeared in [1] with the wrong sign, see [16], Eq. ( 9) in the vector case.This term is scheme independent and it has been calculated in [53] correctly.The axial-vector contribution is newly given.This contribution is not contained in Ref. [1] but agrees with that of Ref. [10].
9 The interference term between non-singlet and pure singlet corrections The radiator for this process is given by It is the same in the vector and axial-vector case, since the amplitude squared has one closed fermion line only.The result differs from that in [1] by the term in [16], Eq. ( 9).
10 Further terms with no logarithmic enhancement at O(α 2 ) For pure vector couplings there are as well fermion-pair production contributions corresponding to the terms |B| 2 , BC and BD in the case of electrons and for |B| 2 for heavier radiated fermions, cf.[10].They have no logarithmic enhancement and were not considered in [1].The B-diagrams are shown in Figure 9 Figure 9: The e + e − annihilation graphs into a fermion pair and a gauge boson (process B).
The corresponding radiator is given by in the vector and axial-vector cases, with and The contributions of the diagrams AB vanish due to Furry's theorem in the vector case.We performed the calculation in D = 4 dimensions keeping the fermion masses, which were set to zero at the end of the calculation.We agree with the results of Ref. [10].The massive operator matrix element vanishes for these processes and therefore the massive and the massless result have to agree according to the factorization theorem postulated in Ref. [12].Due to this, these contributions have not been included in Ref. [2].
11 Contributions due to Soft Photon Exponentiation beyond O(a 2 ) The resummation of the soft corrections has been considered early [1].Here the idea is to resum first the leading distribution-valued contributions and The inclusive cross section reads then cf. [47].The soft-resummed contributions for O(a 3 ) and higher are One may extend the soft photon exponentiation by including as well the soft production of e + e − pairs according to the non-singlet process described in Section 5 in the region z → 1. this modifies the term cf. [54].The Mellin inversion leading to (117) has been calculated in Ref. [55], see also [47].These are both leading order resummations.One may as well resum the logarithms ln k (z)/z in the small z region, cf.e.g.[56] leading to associated Bessel functions, as has been known in QCD before, see e.g.[57].
The leading logarithmic orders O((aL) k ), which are process independent, can be treated rather straightforwardly to rather high orders, accounting both for the non-singlet and singlet contributions, cf.[56,[58][59][60][61][62][63][64][65].These corrections include the resummations mentioned and do even account for more contributions by resumming as well all collinear contributions according to the QED evolution equations.Note, however, that the sub-leading contributions always require the inclusion of mass effects due to massive OMEs, cf.[2].

Conclusions
The O(a 2 ) initial state radiative corrections to the process e + e − → γ * /Z * have been computed in a direct calculation without neglecting the electron mass against the s-channel energy of the process.The expansion in the ratio m 2 e /s has been only performed in a very late stage of the calculation by controlling the result based on precision numerics in mathematica with the complete result.The corrections can be grouped into four main processes, I-IV, as already done in Ref. [1], with the addition of non-logarithmic terms and terms due to soft-exponentiation for the contributions beyond O(a 2 ).Furthermore, one has to account for differing axial-vector contributions in some of the channels.For the processes I-IV we find differing results for the non-logarithmic terms of O(a 2 ) given in [1], while we agree in the logarithmic contributions and those of O(a).On the other hand, we agree with the results of Ref. [2].In the case of process II we agree with Refs.[1,51] if the initial state fermion radiation concerns µ + µ − or heavier lepton or quark pairs.We also agree for the pure-singlet interference terms with a result in Ref. [53].Furthermore, we agree with non-logarithmic corrections derived first for the Drell-Yan process [10].
The present calculation proofs, here for QED, that the massive Drell-Yan process factorizes, and we revise an earlier doubt in Ref. [2].The present rather voluminous calculation has been the only way to establish this.Fortunately, mathematical methods are now available to perform the corresponding integrals analytically and allow to represent them as iterated integrals of square root-valued alphabets, carrying real parameters.It is this representation which finally allows the controlled limit m 2 e /s → 0 for the power corrections.A part of the integrals are incomplete elliptic integrals and generalizations thereof, which does not lead to a further sophistication, since the corresponding integrals are still iterative.Numerical illustrations of the present results have been given in Ref. [15].already.
The numerical accuracy to which both the Z boson mass and width are planned to be measured at the FCC ee is rather high.It amounts to ∼ 100 keV systematic uncertainty, with a much higher statistical precision.It is clear from Ref. [15] that the O(a 2 ) ISR corrections will not yet be sufficient to cope with this accuracy.Therefore, even higher order corrections have to be calculated to sub-leading levels, cf.[66].

A The phase space integrals
In the following we describe the parameterization of the phase space integrals both for the case of fermion pair and photon pair radiation, followed by explicit expressions obtained after the angular integrations.Here the setup is similar to that used in [1,53].

A.1 Fermion Pair Radiation
For massive fermion pair radiation we only encounter 2 → 3 scattering with the kinematics We also introduce the invariants which satisfy the identity The phase space integral is given by In deriving these relations, the identities were used.The integration variables are transformed according to and the symmetry of the angular integration allows to transform The phase space boundaries are given by where the explicit expressions for s − 3 and s + 3 are given by We can also change the order of integration in which case we obtain with the explicit expressions We can use the following parameterization of the vectors: with the abbreviation c(x) = cos(x) and s(x) = sin(x).The missing components of the vectors are given by The direction of the 3-vector component of k 2 is achieved by rotating k 1 with angle χ 0 around the x-axis and then with angle φ 0 around k 1 .It is convenient to transform to the dimensionless variables in the explicit calculations.Since all involved particles are massive, the phase space integrals are convergent and do not need any kind of regularization.

A.2 Photon Radiation
The 2 → 3 scattering can be very similarly parameterized as before.However, the replacements k − → k 1 and k + → k 2 with have to be made.Therefore, the limit m → 0 has to be taken in the expressions given in the previous section.We will give the explicit expressions for completeness.The phase space integral reads The angle between the two photons is given by cos(χ 0 ) = 1 − 2ss (s − s 3 )(s − s 4 ) .
The phase space boundaries simplify to s ≤ s 3 ≤ s.
They are symmetric in s 3 and s 4 .
It is also possible to only radiate one additional photon.In this case the phase space for 2 → 2 scattering is needed.Using the kinematics p − + p − = q + k (146) with k 2 = 0, it is given by dPS 2 = d 4 q d 4 kδ(s − s )δ(k 2 )δ (4)

A.3 The Angular Integrals
For the photon emission graphs we find the following denominators For some denominator structures we have to use partial fractioning.Some cases are trivial, like 1 The more involved ones read For some combinations of denominators we have to interchange the parameterizations of k − and k + in order to arrive at angular integrals of the form (150).
If either l or k are negative we can use the relations given in Eqs.(??,??) for D = 4 to arrive at the angular integrals.If both indices are negative we were not able to find a closed form in D dimensions.For D = 4 we find with X = (aA − bB) 2 − (a 2 − b 2 )(A 2 − B 2 − C 2 ).Note that we agree with the results given in [10,67,68].

Figure 4 :
Figure 4: The e + e − annihilation graphs into two photons and a virtual gauge boson.

Figure 5 :
Figure 5: The double virtual corrections of O(a 2 ) to e + e − annihilation into a virtual gauge boson.The external self-energy corrections are not shown.