Complete electroweak two-loop corrections to Z boson production and decay

This article presents results for the last unknown two-loop contributions to the $Z$-boson partial widths and $Z$-peak cross-section. These are the so-called bosonic electroweak two-loop corrections, where"bosonic"refers to diagrams without closed fermion loops. Together with the corresponding results for the $Z$-pole asymmetries $A_l, A_b$, which have been presented earlier, this completes the theoretical description of $Z$-boson precision observables at full two-loop precision within the Standard Model. The calculation has been achieved through a combination of different methods: (a) numerical integration of Mellin-Barnes representations with contour rotations and contour shifts to improve convergence; (b) sector decomposition with numerical integration over Feynman parameters; (c) dispersion relations for sub-loop insertions. Numerical results are presented in the form of simple parameterization formulae for the total width, $\Gamma_{\rm Z}$, partial decay widths $\Gamma_{e,\mu},\Gamma_{\tau},\Gamma_{\nu},\Gamma_{u},\Gamma_{c},\Gamma_{d,s},\Gamma_{b}$, branching ratios $R_l,R_c,R_b$ and the hadronic peak cross-section, $\sigma_{\rm had}^0$. Theoretical intrinsic uncertainties from missing higher orders are also discussed.


Introduction
The number of Z bosons collected at LEP in the 1990's, 1.7×10 7 , together with SLD data made it possible to determine electroweak pseudo-observables (EWPOs) with high precision: the Z-boson mass M Z , its decay width Γ Z , branching ratios R, forward-backward and left-right asymmetries (or equivalently A f or sin 2 θ f eff ) [1]. At that time, theoretical calculations, which included complete one-loop Standard Model corrections, selected higher order QCD corrections, and partial electroweak two-loop results with intricate QED resummations, were accurate enough to go hand-in-hand with experimental demands [2,3]. However, up to 5 × 10 12 Z-boson decays are planned to be observed at projected future e + e − machines (ILC, FCC-ee, CEPC) running at the Z-boson resonance [4][5][6][7]. These statistics are several orders of magnitude larger than at LEP and would lead to very accurate experimental measurements of EW-POs. Limitations will come from experimental systematics, but they are in many cases estimated to be improved by more than an order of magnitude compared to the LEP experiments [4][5][6][7]. This raises a new situation and theoretical calculations must be much more precise than assumed before [8,9]. The improved precision will provide a platform for deep tests of the quantum structure of nature and unprecedented sensitivity to heavy or super-weakly coupled new physics.
As an important step towards that goal, this article reports on the completion of such calculations at the two-loop level in the Glashow-Weinberg-Salam gauge theory, known as the Standard Model (SM) [10][11][12].
The first non-trivial study of electroweak (EW) loop effects was the calculation of the large quadratic top quark mass contribution to the Z and W propagators at one-loop order [13]. A few years later, the on-shell renormalization scheme as it is used today [14] and the notion of effective weak mixing angles [15] were introduced, and the scheme was used for calculations of the W ± and Z boson masses [16]. The complete one-loop corrections to the Z decay parameters were derived in Refs. [17][18][19][20], and those to the W ± width in Refs. [19,21,22]. Through the years of LEP and SLC studies, the effects of EW corrections became visible in global fits of the SM parameters [1][2][3]23]. Global fits to EW precision measurements allowed to predict the mass of the top quark and the Higgs boson prior to their discoveries at Tevatron in 1995 [24,25] and at the LHC in 2012 [26].
At future e + e − colliders, EWPOs will again play a crucial role. These include the total and partial widths of the Z boson and the Z-boson couplings. The latter can be extracted from measurements of the cross-section and polarization and angular asymmetries of the processes e + e − → (Z) → ff . Here f stands for any SM lepton or quark, except the top quark, whereas the notation (Z) is supposed to indicate that the amplitude is dominated by the s-channel Z-boson resonance, but there is contamination from photon and two-boson backgrounds.
Already for the precision achieved at LEP and SLC, the calculation of loop corrections beyond the one-loop order was necessary to keep theory uncertainties under control. Specifically, these included two-loop O(αα s ) [27][28][29][30][31] and fermionic O(α 2 ) [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46] corrections to the Fermi constant, which can be used to predict the W -boson mass, and to the Z-pole parameters. Here α refers an electroweak loop order, whereas "fermionic" denotes contributions from diagrams with at least one closed fermion loop. In addition, leading three-and four-loop results, enhanced by powers of the top Yukawa coupling y t , were obtained at order O(α t α 2 s ) [47,48], O(α 2 t α s ), O(α 3 t ) [49,50], and O(α t α 3 s ) [51][52][53], where α t = y 2 t /(4π). For the EW two-loop corrections, the calculation of the fermionic contributions was a natural first step, since these are numerically enhanced by the numbers of flavors and colors and by powers of y t . Moreover, the fermionic two-loop diagrams are relatively simpler than the full set. For example, the latter includes non-planar vertex topologies, which are absent in the former. The remaining bosonic two-loop corrections to the Fermi constant and the leptonic effective weak mixing angle, sin 2 θ ℓ eff , have subsequently been presented in Refs. [54][55][56][57][58][59][60], and more recently also for the weak mixing angle in the bb channel [61].
While the numerical effects of the bosonic two-loop corrections are relatively small compared to the current experimental precision from LEP and SLC, their inclusion will become mandatory for future e + e − colliders. Thus the computation of the full two-loop corrections for all Z-pole EWPOs is an important goal. This article completes this goal by presenting the remaining bosonic O(α 2 ) contributions to the Z-boson total and partial widths, and the hadronic Z-peak cross-section within the SM. This has been achieved by using the numerical integration methods discussed in Ref. [61], with some technical improvements.
The paper is organized as follows. After a brief review of the field theoretic definition of the relevant observables in section 2, the technical aspects of the two-loop calculation are described in section 3. The numerical impact of the bosonic EW two-loop corrections is demonstrated in section 4. In particular, results for the total and partial Z widths, several commonly used branching ratios, and the hadronic Z-peak cross-section are given in terms of simple parameterization formulae, which provide an accurate description of the full results within the currently allowed ranges of the input parameters. Finally, the theory uncertainty from missing three-and four-loop contributions is estimated in section 5, before concluding in section 6.

Definition of the observables
The amplitude for e + e − → ff near the Z pole, √ s ≈ M Z can be written in a theoretically well-defined way as a Laurent expansion around the complex pole where M Z and Γ Z are the on-shell mass and width of the Z boson, respectively. According to eq. (1), the approximate line shape of the cross-section near the Z pole is given by σ ∝ It is important to note that this differs from the line shape used in experimental analyses, which is of As a result, the parameters in eq. (1) differ from the experimental mass M Z and width Γ Z from LEP by a fixed factor [62]: Numerically, this leads to The total width, Γ Z , can be extracted from the condition that the Z propagator has a pole at s = s 0 , leading to where Σ Z (s) is the transverse part of the Z self-energy. Using the optical theorem, it can also be written as [45,46] Here the sum runs over all fermion types besides the top quark, f = e, µ, τ, ν e , ν µ , ν τ , u, d, c, s, b, and N f c = 3(1) for quarks (leptons). The radiator functions R V,A capture the effect of final-state QED and QCD corrections. They are known up to O(α 4 s ) and O(α 2 ) for massless external fermions and O(α 3 s ) for the kinematic mass corrections [63][64][65]. For the results shown in this article, the explicit form given in the appendix of Ref. [46] has been used.
The remaining radiative corrections are IR finite and contained in the form factors F f V,A . These include massive EW corrections as well as mixed EW-QCD and EW-QED corrections. The bosonic two-loop contributions, which are of interest for this article, contribute according to [46]: where v f and a f are the effective vector and axial-vector couplings, respectively, which include Zff vertex corrections and Z-γ mixing contributions. Σ ′ Z denotes the derivative of Σ Z , and the loop order is indicated by the subscript (n).
It should be pointed out that v f , a f and Σ Z as defined above Here v Z f and a Z f are the one-particle irreducible Zff vector and axial-vector vertex contributions, respectively, whereas v γ f and a γ f are their counterpart for the γff vertex. Furthermore, Σ V1V2 denotes the one-particle irreducible V 1 -V 2 self-energy.
Another important quantity is the hadronic peak cross section, σ 0 had , which is defined as the total cross section for e + e − → (Z) → hadrons for s = M 2 Z , after removal of s-channel photon exchange and box diagram contributions, as well as after the de-convolution of initial-state and initial-final interference QED effects [1,2]. The impact of the bosonic two-loop vertex corrections on σ 0 had is given by [45,46] σ 0 The form factors F f V,A are understood to include appropriate counterterms such that they are UV finite. Throughout this work, the on-shell renormalization scheme is being used, which defines all particle masses in terms of their (complex) propagator poles and the electromagnetic coupling in terms of the photon-electron vertex in the Thomson limit. A more detailed discussion of the relevant counterterms can be found in Ref. [40].
As a consequence of this renormalization scheme, the EW corrections are organized as a series in the electromagnetic coupling α, rather than the Fermi constant G µ . Instead, G µ will be used to compute M W within the SM, including appropriate two-loop and partial higher-loop corrections. After this step, the remaining input parameters for the prediction of the Z coupling form factors are M Z , M H , m t , G µ , α, α s and ∆α. Here ∆α captures the running of the electromagnetic coupling induced by light fermion loops. It is defined through α(M 2 Z ) = α(0)/(1 − ∆α), where α(q 2 ) is the coupling at scale q 2 . The contribution from leptons to ∆α can be computed perturbatively and is known at the three-loop level [66], ∆α lept (M Z ) = 0.0314976. On the other hand, the quark contribution is nonperturbative at low scales and thus is commonly derived from experimental data. For recent evaluations of ∆α (5) had , see Refs. [67][68][69]. As a reference value, ∆α Additionally, Γ Z and Γ W are needed as inputs to convert M Z and M W to the complex pole scheme, see eq. (2). Furthermore, the radiator functions and m τ to account for kinematic fermion mass effects in the final state, whereas the masses of electron, muon, neutrinos, and u/d/s quarks can be taken as zero to very good approximation. In contrast to all other masses in this work, the MS masses are used for the bottom and charm quarks, since their on-shell counterparts are poorly defined.

Calculation of two-loop vertex corrections
For the calculations we followed the strategy developed in Ref. [61], where the two-loop bosonic corrections to the bot-tom quark weak mixing angle, sin 2 θ b eff , were obtained. In fact, the Zbb vertex is the technically most difficult case due to the larger number of mass scales in that problem compared to other flavors. Details are described there and also in [70][71][72]. On the other hand, for the computation of the Z width we are faced not only with ratios v f (2) /a f (2) , but also with sums of powers of v f (2) and a f (2) , see (6) and (7). This leads to the occurrence of extra integrals which cancel out in the ratios v/a.
The complete set of two-loop diagrams required for this calculation have been generated with the computer algebra package FeynArts 3.3 [73]. They can be divided into several categories. The renormalization counterterms require two-loop self-energies with Minkowskian external momenta, In addition, there are two-loop vertex integrals with one non-vanishing external momentum squared, The two-loop self-energy integrals needed for the renormalization procedure and the vertex integrals with selfenergy sub-loops have been computed using the dispersion relation technique described in Refs. [60,74,75]. The remaining bosonic two-loop diagrams amount to about one thousand integrals with a planar or non-planar vertex topology.
Some new classes of integrals compared to the sin 2 θ b eff case are met. They are simpler from a numerical point of view than those solved in Ref. [61]. For instance, there are various oneand two-scale integrals with internal W propagators, which improves the singular threshold behaviour of integrals with only Z propagators. There are altogether about one hundred integrals of this kind with different permutations of propagators, including the tensor integrals. As an example of one of the most difficult cases, the SD method for integrals from Fig. 1 in [61] gives an accuracy of up to four relevant digits. Using the MB method, these diagrams are equivalent to up to 4-dimensional MB integrals, which can be calculated efficiently with eight relevant digits by MBnumerics.
In select cases, like those described above, the MB approach is uniquely powerful. This statement applies to several hundred integrals. In the majority of integrals, though, the SD method is presently more efficient than the MB approach, mainly due to the smaller number of integration variables. For our semiautomatized calculation of massive 2-loop vertices the avail-ability of two complementary numerical methods with a large overlap was crucial.

Numerical results
In this section, numerical results for bosonic two-loop corrections are compared to and combined with all other known corrections to the Zff vertices. These are  (6) and (7); • Non-factorizable O(αα s ) vertex contributions [95][96][97][98][99][100], which cannot be written as a product of EW form factors F V,A and final-state radiator functions R V,A , but instead are added separately to the formula in (5).
These are applied to a range of EWPOs: The partial Z widths, Γ f ≡ Γ(Z → ff ), as well as total width, Γ Z , various branching ratios, and the hadronic peak cross-section σ had . The full electroweak two-loop corrections for the leptonic and bottomquark asymmetries have been published previously [58][59][60][61] and are not repeated here. Nevertheless, as a cross-check we reproduced the result for the leptonic asymmetry and found agreement with Refs. [58,59] within intrinsic numerical uncertainties. Moreover, with the methods described here we can produce results for the bosonic two-loop corrections to sin 2 θ ℓ eff with four robust digits of precision, which exceeds the accuracy obtained with asymptotic expansions as in Ref. [58]. As discussed above, the gauge-boson mass renormalization has been performed in accordance with the complex-pole scheme in eq. (1). However, for the sake of comparison with the wider literature, the numerical results below are presented after translating to the scheme with an s-dependent width. In other words, results are shown for un-barred quantities, such as Γ Z in eq.

Partial widths
Let us begin by presenting results for a fixed value of M W as input, instead of calculating M W from G µ . This more clearly illustrates the impact of the newly completed O(α bos ) 2 corrections. Table 2 shows the contributions from different loop orders to the SM prediction of various partial Z widths. As is evident from the table, the two-loop EW corrections are significant and larger than the current experimental uncertainty (2.3 GeV for Γ Z [1]). The newly calculated bosonic corrections O(α 2 bos ) are smaller but still noteworthy. They amount to half of all known leading three-loop QCD corrections O(α t α 2 s , α t α 3 s , α 2 t α s , α 3 t ), even though the latter are enhanced by powers of α s , α t and N f . Table 3 shows the SM predictions obtained if one uses G µ as an input to compute M W , based on the results of [39,40,[54][55][56][57]. Each line of the table adds an additional order of perturbation theory to the previous line, using the same order for the Zff vertex corrections and the calculation of the W mass 1 .
The O(α 2 bos ) correction to Γ Z , corresponding to the difference between the last two rows in Table 3, amounts to 0.34 MeV, which is more than three times larger than its previous estimation [46]. An updated discussion on how this knowledge changes the intrinsic error estimations will be given in section 5.

Ratios
The experimental results from LEP and SLC are typically not presented in terms of partial widths for the different final   Table 3: Results for Γ Z and σ 0 had , with M W calculated from Gµ using the same order of perturbation theory as indicated in each line. In all cases, the complete radiator functions R V,A are included.
states. Instead, this information is captured in the form of various branching ratios. The most relevant ones are (12) where Γ ℓ = 1 3 (Γ e + Γ µ + Γ τ ), and Γ had is the partial width into hadronic final states, which at the parton level is equivalent to q Γ q (q = u, d, c, s, b).
In addition, the hadronic peak cross-section (11) is, to a good approximation, defined as a ratio of partial widths and the total Z width.
Numerical results for σ 0 had and the ratios in (12) are given in Tab. 3 and Tab. 4, respectively, again broken down to different orders of radiative corrections. These quantities are less sensitive to higher loop effects than Γ Z , since there is a partial cancellation between the corrections in the numerators and denominators of the ratios. Thus the influence of the new bosonic corrections on all branching ratios R ℓ , R c , R b and on σ 0 had is about 0.02% or less, which is far below the current experimental errors: R ℓ = 20.767 ± 0.025, R c = 0.1721 ± 0.0030, R b = 0.21629 ± 0.00066, and σ 0 had = 41.541 ± 0.037 nb [1]. However, these are at the level of sensitivity of proposed measurements of R b at future e + e − colliders [4][5][6][7]

Parameterization formulae
While the tables above only contain numbers for a single benchmark point, the results for a range of input values can be conveniently expressed in terms of simple parameterization formulae. The coefficients of these formulae have been fitted to the full calculation results on a grid that spans the currently allowed experimental ranges for each input parameter. Here the full calculation includes all higher-order corrections listed at the beginning of section 4 for the partial widths, branching ratios and the peak cross-sections, and with M W calculated from G µ to the same precision 2 . For all EWPOs reported here, the same form of parameterization formula is utilized: As before, M H , M Z , m t and ∆α are defined in the on-shell scheme, using the s-dependent width scheme for M Z (to match the published experimental values), while α s is defined in the MS scheme. The dependence on m b , m c and m τ is negligible within the allowed ranges for these quantities.
The fit values of the coefficients for the different EWPOs are given in Tab. 5. With these parameters, the formulae provide very good approximations to the full results within the ranges M H = 125.1 ± 5.0 GeV, m t = 173.2 ± 4.0 GeV, α s = 0.1184 ± 0.0050, ∆α = 0.0590 ± 0.0005 and M Z = 91.1876 ± 0.0042 GeV, with maximal deviations as quoted in the last column of Tab. 5. As can be seen from the latter, the accuracies of the fit formulae are sufficient for the forseeable future.
Following Refs. [46,102], the size of these terms may be estimated by assuming that the perturbation series approximately is a geometric series. In this way one obtains where the known leading large-m t approximations have been subtracted in the numerators. For the example of the total Z width, these expressions lead to An additional source of theoretical uncertainty stems from the unknown O(α 5 s ) final-state QCD corrections and three-loop mixed QED/QCD final-state corrections of order O(αα 2 s ) and O(α 2 α s ). In [46] they were found to be sub-dominant, and the estimates can be taken over from there without change. Combining these findings with eqs. (15) in quadrature, the total theory error adds up to δΓ Z ≈ 0.4 MeV. Compared to the previous theory error estimate δΓ Z ≈ 0.5 MeV [46] one observes a slight decrease due to the knowledge of the bosonic corrections calculated in this work.
In addition to the elimination of an uncertainty associated with the previous unknown O(α 2 bos ) corrections, the values in the first and second rows of (15) also shifted since the full O(α 2 ) corrections used in (14) were not available before. These shifts conspire to result in a reduction of the uncertainty estimate for these two error sources. The corresponding error estimates for the partial widths are shown in Table 6. For the ratios (R ℓ , R c and R b ), the theory uncertainty has been estimated from the partial widths using simple Gaussian error propagation.
The theory uncertainty for the hadronic peak cross-section is dominated by a non-factorizable contribution stemming from the imaginary part of the Z-boson self-energy [46]. This nonfactorizable term does not receive any bosonic two-loop corrections, so that its error estimate can be taken from Ref. [46] without change: Adding these in quadrature leads to the overall uncertainty estimate of δσ 0 had ≈ 6 pb.

Summary
In this work the bosonic two-loop electroweak corrections, O(α 2 bos ), to Z boson production and decay parameters are presented for the first time. Updated theory error estimations are given, which are slightly reduced due to the newly available full two-loop corrections. We expect that the numerical integration methods used here can be extended to compute the full three-loop corrections to Z-pole EWPOs. For a more detailed discussion of future projections, see Ref. [8,9]. However, this is very demanding and needs more effort and resources. Further, at this level of complexity independent cross-checks by different groups, using independent calculations and approaches, are welcome.
It should be noted that the O(α 2 bos ) correction for the total Z decay width appears to be relatively large compared to previous estimates based on the knowledge of the lower order result O(α bos ). A similar observation concerns the bosonic two-loop corrections to A b . This means that all estimations at this level of accuracy should be taken with a grain of salt. Therefore, explicit calculations are important even for contributions that were previously estimated to be subdominant.
At this point we should mention that we did not consider the theoretical efforts needed to unfold the large QED corrections from the measured real cross sections in the Z peak region and to extract the EWPOs studied here in detail. For LEP, this was based on tools such as the ZFITTER package [103][104][105] and was discussed carefully e. g. in Refs. [1,2,106]. The correct unfolding framework for extracting 2 → 2 observables at accuracies amounting to about 1/20 of the LEP era certainly has to rely on the correct treatment of Laurent series for the Z line shape as is discussed e. g. in [107][108][109][110].
The 1-loop corrections to the Z boson parameters were determined in the 1980s [17]. Today, 33 years later, while the present study finalizes the determination of the electroweak twoloop corrections to the Z-boson parameters, we are already faced with the need of more precision in the future. eff and indirect determination of the Higgs boson mass, Phys. Rev. Lett. 93 (2004)