Electroweak pseudo-observables and Z-boson form factors at two-loop accuracy

We present Standard Model predictions for the complete set of phenomenologically relevant electroweak precision pseudo-observables related to the Z-boson: the leptonic and bottom-quark effective weak mixing angles $\sin^2\theta_{\rm eff}^\ell$, $\sin^2\theta_{\rm eff}^b$, the Z-boson partial decay widths $\Gamma_f$, where $f$ indicates any charged lepton, neutrino and quark flavor (except for the top quark), as well as the total Z decay width $\Gamma_Z$, the branching ratios $R_\ell$, $R_c$, $R_b$, and the hadronic cross section $\sigma_{\rm had}^0$. The input parameters are the masses $M_Z$, $M_H$ and $m_t$, and the couplings $\alpha_s$, $\alpha$. The scheme dependence due to the choice of $M_W$ or its alternative $G_\mu$ as a last input parameter is also discussed. Recent substantial technical progress in the calculation of Minkowskian massive higher-order Feynman integrals allows the calculation of the complete electroweak two-loop radiative corrections to all the observables mentioned. QCD contributions are included appropriately. Results are provided in terms of simple and convenient parameterization formulae whose coefficients have been determined from the full numerical multi-loop calculation. The size of the missing electroweak three-loop or QCD higher-order corrections is estimated. We briefly comment on the prospects for their calculation. Finally, direct predictions for the $Z{\bar f}f$ vector and axial-vector form-factors are given, including a discussion of separate order-by-order contributions.


Introduction
In 2018 we celebrated 50 years of the Standard Model of elementary particles. The basics of the model were formulated and experimentally validated in the 1960s/70s. The next decade brought an intensive development of the calculation of quantum field theoretical radiative corrections in that model and in its alternatives. An experimental highlight in this context was the e + e − -collider LEP, which enabled us to check the Standard Model at an accuracy of better than the per-cent level, which corresponds to effects from more than one electroweak and two QCD loop orders. This proved, for the first time in a systematic way, the Standard Model as a quantum field theory. LEP 1 was running, from Summer 1989 to 1995, at and around the Z-boson peak. The expectation for the experimental precision of M Z and Γ Z was 20 MeV in 1986 [1] and reached finally 2 MeV [2]. This precision tag was extremely important because M Z is one of the Standard Model input parameters to the commonly used on-mass-shell renormalization scheme. Indeed, the experimental accuracy of M Z triggered much of the precision loop calculations, including the prediction of the top quark and Higgs masses prior to their discoveries from loop corrections to LEP observables in the Standard Model, see Refs. [3][4][5] (as well as Refs. [6,7] for an overview of the current state of the art). Data from the Z peak and the Z resonance curve (the Z line shape) allow to measure a large variety of observables, such as M Z , Γ Z , cross-sections for different two-fermion final states and their ratios and angular asymmetries, together with radiation of (sufficiently soft) photons, gluons, etc. From the real observables, the so-called electroweak pseudo-observables (EWPOs) are extracted by means of a de-convolution of initial-state radiation and subtraction of backgrounds. The fine details of relating EWPOs to real cross-sections at LEP 1 precision are described in detail in Ref. [8] and references quoted therein.
On occasion of the 50 th anniversary of the Standard Model, many of the related crucial developments were remembered at the conference "SM@50" [9]. This article is devoted to the state-of-the art calculation of the Standard Model (SM) corrections to the Zf f -vertex and their inclusion into the predictions for the various EWPOs. We will mostly focus on recent advances in the calculation of the electroweak two-loop terms. QCD contributions, which are known up to four-loop order, have also been taken into account in the results presented here, but we will refer to the literature for further details.
The Z resonance curve can be described theoretically by writing the S-matrix elements as a Laurent series in the center-of-mass energy squared s (also called Smatrix ansatz). This Laurent series ansatz is worked out up to two loops [10][11][12][13][14][15][16][17][18]. The coefficients of the leading series term contain the Z vertex form factors. Their one-loop corrections were studied in the 1980s; first with massless fermions [19][20][21][22], and slightly later with the full mass dependence of the Standard Model [23][24][25][26].
After several papers on approximate/partial higher-order corrections, the complete two-loop weak corrections were determined in a series of papers from 2004 to 2018 [27][28][29][30][31][32][33][34][35][36][37]. The correct formulation of the interplay of the 2→2 loop corrections with higher order real QED corrections in the S-matrix approach, also called un-folding of the effective 2→2 Born terms from the realistic 2→n observables, is a topic on its own. It was first studied in Refs. [38][39][40][41], but its extension beyond two-loop level will require more work [8,42,43]. The numerically relevant two-loop and partial higher-order corrections were included in the final analysis of LEP 1 data [44]. 1 The theoretical advances described here go beyond the Standard Model theory used for physics at LEP 1 [45,46] but will be needed for the FCC-ee Tera-Z project [6,7,[50][51][52] whose unique experimental precision calls for perturbative predictions at three electroweak loops together with corresponding QCD terms.
In this work, the following pseudo-EWPOs are discussed: The partial widths Γ f for Z-boson decay into ff final states; the total Z width Γ Z ; the branching fractions Γ f /Γ Z ; the total hadronic Z-pole quark-pair production cross-section σ 0 had ; and the effective weak mixing angles sin 2 θ f eff , defined from the ratio of vector-and axial Zboson couplings. The precise definition of all these pseudo-observables will be given below. Pseudo-observables differ from real observables by removing from the former the effects of initial-state and initial-final QED radiation, as well as non-resonant photon-exchange, box and t-channel contributions [6,44].
The Standard Model predictions for Z-pole pseudo-observables can be constructed in terms of the following three theoretical building blocks [29]: where 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 counterparts for the γff vertex. The Σ V 1 V 2 denotes the one-particle irreducible Here I 3 f and Q f are the weak isospin and electric charge (in units of the elementary charge e > 0) of the fermion f , respectively. s W and c W are the sine and cosine of the weak mixing angle, respectively, and the subscript (0) is used to denote tree-level order.
For the theory calculations, these building blocks must be evaluated at the complex Z pole [11][12][13]17] where M Z and Γ Z are the on-shell mass and width of the Z-boson, respectively. Note that M Z and Γ Z differ from the mass M Z and width Γ Z reported in publications of LEP, Tevatron and LHC experiments by a fixed factor [38,53]:

Z-boson decay width, branching ratios and cross-sections
The width of the Z boson, Γ Z , is related to the imaginary part of the Z self-energy. Using the optical theorem, one can derive the following expression for Γ Z [33,35]: Here N f c is the color factor and R V,A are radiator functions that capture final-state QCD and QED corrections, see section 7 in Ref. [56], whereas the remaining electroweak and mixed electroweak-QCD corrections are contained in the form factors F f V,A . Up to two-loop accuracy, the form factors can be written as follows [35]: 3) where Σ Z and Σ Z are shorthand expressions for dΣ Z /ds and d 2 Σ Z /ds 2 , respectively. In addition to the partial widths, certain branching ratios are of phenomenological importance: Here Γ had = f =u,c,d,s,b Γ f . Further, the cross-section for e + e − → hadrons at the Z peak can be expressed in terms of partial widths [35], Here σ 0 had,non−res accounts for non-resonant photon-exchange, box and t-channel contributions. Furthermore, δX occurs from higher-order terms of the Laurent expansion of the full amplitude around the complex pole s 0 . At two-loop order, δX can be written as δX (2) , where subscripts (n) indicate the loop order. In the limit m f M W (f = t), it is given by (2.7) It is important to note that eq. (2.6) assumes that initial-state photon radiation effects have been removed by means of a de-convolution procedure, see e.g. Ref. [57]. Results for the partial and total Z widths, branching ratios and σ 0 had including the full two-loop corrections have first been published in Ref. [37]. They can be expressed in simple parameterization formulae, which are adequate for most phenomenological applications. Here, we present slightly more complicated formulae that cover an extended numerical range of input parameters: Here ∆α is the shift in the running electromagnetic coupling α(q 2 ) from q 2 = 0 to M 2 Z , defined by α(M 2 Z ) = α(0)/(1 − ∆α). It can be divided into a leptonic and a hadronic part, ∆α = ∆α lept + ∆α had . ∆α lept has been computed to three-loop order [59], whereas ∆α had contains non-perturbative hadronic contributions, which are commonly extracted from data [60][61][62]. We neglect the light fermion masses m f , f = t, everywhere besides in ∆α and (at leading power) in the radiator functions R f V/A . The W boson mass M W can be computed from the Fermi constant G µ (see eqs. (6)- (8) in the arXiv version of Ref. [47]) and thus is not listed as an independent input parameter. Both G µ and α, the electromagnetic fine structure constant in the Thomson limit, are known with very small uncertainties, and thus we use their central experimental values [58] without any uncertainty interval. The fitting formulae for the EWPOs have the form The coefficients X 0 and a 1 , . . . a 16 are obtained from fits to a grid of 8750 data points of the full computation. The latter includes • Complete one-loop corrections [23], which have been re-computed for this work, and full two-loop [33,35,37] electroweak corrections; • Corrections of order O(αα s ) to vector-boson self-energies [63][64][65][66][67], which have been re-evaluated for this work; • Higher-loop QCD corrections in the limit of a large top Yukawa coupling y t , of orders O(α t α 2 s ) [74,75], O(α 2 t α s ), O(α 3 t ) [76,77], and O(α t α 3 s ) [78][79][80], where α t ≡ y 2 t /(4π).
Both the Z vertex corrections as well as the prediction of M W from G µ have been computed to this same level of perturbation theory. Numerical values for the coefficients are given in Tab. 1. Some of the numbers for X 0 deviate slightly in the last digit from those in Ref. [37]. This is due to the larger grid of input parameters used here, which can exert a pull on the fit parameters. The differences are well within the accuracy quoted in the last column of Tab. 1.

Asymmetries and effective weak mixing angles
The effective weak mixing angle for the Zff vertex is defined, from the theory side, as . (3.1) Here M 2 W and M 2 Z are the real parts of the complex pole of the W and Z propagators, respectively. They are related to the masses commonly reported by experiments at LEP, Tevatron, LHC according to eq (1.6). Moreover, Q f denotes the electric charge of the fermion f .
The effective weak mixing angles can be extracted from a range of asymmetries [8], defined from effective Born two-particle cross-sections, including the left-right asymmetry and the forward-backward asymmetry Here σ e L and σ e R are the cross-sections for e + e − → ff for left-and right-handed polarized electron beams, respectively, whereas σ cos θ>0 and σ cos θ<0 denote the crosssection for f restricted to the forward and backward hemisphere, respectively. Furthermore, A non−res X accounts for the non-resonant photon-exchange, box and t-channel contributions.
The most precisely measured effective weak mixing angles are the leptonic effective weak mixing angle sin 2 θ eff (extracted from A LR ) and the bottom-quark one, [57]. Standard Model predictions for sin 2 θ eff including the full two-loop corrections have been presented originally in Ref. [28,29,31]. We reproduced by an independent calculation the contribution of the bosonic electroweak two-loop corrections using the methods of Ref. [37]. The corrections can be expressed in terms of a weak form factor ∆κ (α 2 ,bos) , where The comparison with Ref. [28] is shown in Tab. 2, which demonstrates that the two calculations agree to an accuracy of O(10 −7 ), which implies an accuracy of better than 10 −7 for sin 2 θ eff . The full two-loop corrections for sin 2 θ b eff have been presented first in Ref. [36].  In the following, we present simple parameterization formulae for sin 2 θ eff and sin 2 θ b eff , which cover the extended range of input parameters of eq. (2.8). The parameterization formula provides a good description of the full result in the parameter region (2.8). Values for the coefficients are obtained by fitting (3.6) to a grid of 8750 data points. Table 3 shows the result of a fit to a calculation that includes all known corrections: • Complete one-and two-loop electroweak corrections, (see Refs. [21,23,27,28,[30][31][32]36] for the original references); • Corrections of order O(αα s ) to vector-boson self-energies [63][64][65][66][67], which we have re-evaluated for this work; • Non-factorizable O(αα s ) Zbb vertex contributions [68][69][70][71][72][73], which do not cancel in the ratio v b /a b ; • Higher-loop corrections in the limit of a large top Yukawa coupling y t , of orders O(α t α 2 s ) [74,75], O(α 2 t α s ), O(α 3 t ) [76,77], and O(α t α 3 s ) [78][79][80] where α t ≡ y 2 t /(4π).
As indicated by the last column in the table, the largest deviation of the fit formulae from the full result is O(few × 10 −6 ), while for most of the parameter region in (2.8) the agreement is better than 10 −6 . The careful reader may realize that the parameterization for sin 2 θ b eff in Table 3 deviates slightly from Eqs. (20,22) in [36]. The difference is due to the larger grid of data points used here. A fit formula is,  Table 4. Goodness of fit for some chosen EWPOs, compared with the envisaged precision measurements for Γ Z and sin 2 θ eff (statistical errors), and sin 2 θ b eff (systematic errors) at the collider projects FCC-ee Tera-Z [83], CEPC [84] and ILC/GigaZ [85]. The values of maximal deviations are taken from Tabs. 1 and 3. The entry "EXP now" gives the present experimental precision, as known since LEP 1 [44].
obviously, not able to reproduce the data points in a grid perfectly. The fitting aims to find the best average agreement between the data points (which are generated with our full numerical calculation) and the fit formula. A larger grid therefore can lead to some shifts of the coefficients. As a consequence, the formula in [36] will probably be more accurate for input values within the ranges in Tab. 1 there. On the other hand, while the formula here may be a little less accurate within these ranges, it covers a much larger range of input values.
It should also be noted that the fit formula for sin 2 θ eff in Ref. [28] does not include the O(α t α 3 s ) corrections from Refs. [78][79][80], but they are included in the formula presented here.
In Tab. 4 it is shown that the technical accuracy of our fit formulae is adequate for the expected experimental precision of several future e + e − colliders, although it will get modified by anticipated future three-loop electroweak corrections.

Vector and axial-vector Z-boson form factors F f
V and F f A The pseudo-observables discussed in the previous sections aim to be closely related to actual observables, such as cross-sections, branching ratios, or asymmetries. On the other hand, for some purposes it is also useful to have numerical results for the underlying vertex corrections themselves [34], For Tab. 6, on the other hand, the Fermi constant G µ is used as an input instead of (4.1b), and M W is computed from where ∆r has been evaluated to the same orders as given in each column of the table.
More details about the calculation of ∆r can be found in Ref. [47]. As before, the dependence of the Standard Model prediction on various input parameters can be expressed in terms of the simple parameterization formula eq. (2.9). Table 7 shows the numerical values for the coefficients obtained by fitting this formula to the currently most precise computation, including the same corrections as in section 2, except for the final-state QED and QCD radiation effects, i.e.
Note that G µ (rather than M W ) has been used as one input in Tab. 7.

Form fact.
Born  Table 5. Contributions of different perturbative orders to the Z vertex form factors. A fixed value of M W has been used as input, instead of G µ . N n f refers to corrections with n closed fermions loops, whereas α 2 bos denotes corrections without closed fermions loops. Furthermore, α t = y t /(4π) where y t is the top Yukawa coupling.

Form fact.
Born  Table 6. Same as Tab. 5, but with M W calculated from G µ .
The form factor results presented here can be easily augmented to include the effects of some new physics model: Here "SM" denotes the SM contributions discussed in the present paper, while "NP" stands for the new physics correction on top of the SM. Since the existing experimental constraints imply that any possible new physics effect is small, it is sufficient to use the tree-level couplings v f (0) and a f (0) in the interference terms and neglect the |v f,NP | 2 and |a f,NP | 2 terms.

Theoretical error estimates for missing higher order corrections
The main theory uncertainty of the results presented in this paper stems from unknown three-and four-loop corrections. The leading missing orders are O(α 3 ), O(α 2 α s ), O(αα 2 s ), and O(αα 3 s ). Partial results for these contributions, in the limit of a large top Yukawa coupling y t , have already been computed [74][75][76][77][78][79][80]. Therefore, when evaluating the impact of theory uncertainties, it is always implied that we refer to these contributions beyond the leading-y t limit.
There are a number of different methods for assessing theory uncertainties from unknown higher orders, none of which is fully reliable. Rather, they should be taken as an order-of-magnitude estimate of the size of these terms. A convenient and widely applicable method is based on the assumption that the first few orders of the perturbation series approximately follow a geometric series [35,37,86]. In this way one obtains as an ansatz where α t = y 2 t /(4π). Since we are only interested in the missing higher orders beyond the leading large-y t limit, the same leading large-y t approximations have been subtracted in the numerators on the right-hand sides.
The contribution of these estimates to the overall theory error evaluation is shown in Tab. 8 for various pseudo-observables, and in Tab. 9 for the Z-boson form factors. Note that the error estimate for sin 2 θ eff is slightly improved compared to Refs. [28,29] due to the inclusion of O(α t α 3 s ) corrections from Refs. [78][79][80]. Nevertheless, we would also like to remind the reader that any estimate of the theory error from missing higher orders is not a precise prediction. Therefore it is generally desirable to ensure that the theory error is sub-dominant in any phenomenological analysis. Comparing the numbers in Tab. 8 to current measurement results [57,58], this is clearly seen to be the case.

Conclusions and outlook
In this study, we have presented some phenomenologically useful applications of the recently completed electroweak two-loop calculation of Z-boson vertex corrections  Table 9. Same as Tab. 8, but for various form factors. [36,37]. The work collects multi-year efforts of several groups for predictions of the EWPOs related to the Z peak up to electroweak full two-loop accuracy, supplemented by leading QCD higher-order terms. We have determined the two-loop electroweak contributions with a net relative numerical accuracy of about four digits. This ensures that these two-loop results will be known with sufficient accuracy even when adding the next perturbative order, as it might be needed for applications at the next generation of e + e − colliders. For practical applications, the results for the EWPOs, as well as for the Zboson vertex form factors, have been presented in terms of simple parameterization formulae, whose coefficients have been fitted to the full numerical computation. It is planned to include these fitting formulae into a new version of the weak library DIZET of ZFITTER [45,46,48,[87][88][89]. 3 The accuracy of the fitting formulae is less than our full numerical two-loop calculation, but more than sufficient for present-day purposes. For the future FCC-ee Tera-Z project, it may be necessary to provide more precise formulae by including more terms with higher powers of the input parameters.
Finally, we would like to make a few comments on the prospects for the calculation of electroweak three-loop corrections, which will be necessary for the level of precision foreseen for FCC-ee and similar e + e − collider proposals. This electroweak third order, by itself, will be needed with only about two digits accuracy [90]. The generation of the amplitudes for O(10 4 −10 5 ) diagrams as well as the evaluation of the Lorentz and Dirac algebra are routine tasks performed by computer algebra programs, and they should be straightforward with increased computing power in the future. Potential specific problems related to the treatment of γ 5 at three-loop level have to be controlled [91].
The most challenging problem will certainly be the stable numerical computation of three-loop Feynman integrals with several different internal mass scales. At two loops, we did not perform any reduction of the Feynman integrals to a smaller number of masters and thus had to calculate about 1000 previously unknown numerical integrals. In the next perturbative order, it may be advantageous to perform such a reduction to masters, given the ever increasing performance of programs like KIRA [92,93], Reduze 2 [94], FIRE [95,95], and LiteRED [96].
For the calculation of the (master) integrals themselves, it is desirable to have a procedure to automatically isolate and treat the ultra-violet and infra-red singularities. Although there is rapid progress in several analytical approaches to complicated loop integrals [6,7], one has to expect that the more complicated ones will have to be done numerically. An additional complication is that for integrals with physical cuts, their stable numerical evaluation becomes more challenging. There are two kinds of software packages available that address these problems, based on either sector decomposition (SD) as realized in the SecDec project [97][98][99][100][101][102] and the FIESTA project [103][104][105][106], or based on Mellin-Barnes (MB) transformations, as implemented in the AMBRE project [107][108][109][110][111][112][113][114][115][116][117][118][119][120]. Sector decomposition is typically advantageous for integrals with many different mass scales, while the MB approach is more efficient for integrals with fewer independent parameters. Both methods certainly have room for crucial improvements. Several other numerical integration methods, as reviewed e.g. in Refs. [6,7], are useful for certain classes of multi-loop integrals, even though they are less general than the SD and MB approaches. Overall, numerical loop integration techniques are well positioned to meet the necessary future precision demands.
An electroweak three-loop result for the Z-peak EWPOs must be accompanied by improved calculations of the corrections needed to translate EWPOs to real observables. These include initial-state and final-state QED corrections and their interference as well as higher-order terms of the Laurent series expansion about the Z resonance pole [6]. The latter will, for example, involve massive two-loop box diagrams. A complete accounting of the required correction terms is still lacking.
To summarize, we have completed the electroweak two-loop predictions for the EWPOs of the Z resonance and collect here an extensive set of new fitting formulae for them. On a longer time scale, the calculation of the next perturbative order for the calculation of the EWPOs will be necessary and, with proper investments, realistically accessible.