The two-loop electroweak bosonic corrections to $\sin^2\theta_{\rm eff}^{\rm b}$

The prediction of the effective electroweak mixing angle $\sin^2\theta_{\rm eff}^{\rm b}$ in the Standard Model at two-loop accuracy has now been completed by the first calculation of the bosonic two-loop corrections to the $Z{\bar b}b$ vertex. Numerical predictions are presented in the form of a fitting formula as function of $M_Z, M_W, M_H, m_t$ and $\Delta{\alpha}$, ${\alpha_{\rm s}}$. For central input values, we obtain a relative correction of $\Delta\kappa_{\rm b}^{(\alpha^2,\rm bos)} = -0.9855 \times 10^{-4}$, amounting to about a quarter of the fermionic corrections, and corresponding to $\sin^2\theta_{\rm eff}^{\rm b} = 0.232704$. The integration of the corresponding two-loop vertex Feynman integrals with up to three dimensionless parameters in Minkowskian kinematics has been performed with two approaches: (i) Sector decomposition, implemented in the packages FIESTA 3 and SecDec 3, and (ii) Mellin-Barnes representations, implemented in AMBRE 3/MB and the new package MBnumerics.


Introduction
This paper reports on the calculation of the bosonic O(α 2 ) corrections to sin 2 θ b eff and to the Z → bb decay asymmetry parameter A b . Here bosonic refers to corrections from diagrams without closed fermion loops. This completes the calculation of their O(α 2 ) electroweak corrections: The fermionic two-loop corrections have been given in Ref. [1].
For the leptonic Z decay asymmetry parameter, the calculation of the complete electroweak twoloop corrections was presented in [2,3]. For the other Z-boson parameters -Γ Z , Γ ν Z , Γ q Z , Γ b Z -, and for A b , the fermionic electroweak two-loop corrections have been determined [1,4,5], but the bosonic electroweak two-loop corrections were yet unknown.
We would like to remind the reader that e + e − -annihilation into fermion pairs is described by a gauge invariant, unitary and analytic scattering amplitude [6]: The proper formalism for its perturbative calculation has been derived in [7,8] and its application at two-loop accuracy is described in Ref. [3] and references therein. The amplitude (1) has a Breit-Wigner resonance form with fixed mass M Z and width Γ Z . A Breit-Wigner ansatz with an energydependent width Γ Z as it is used in most experimental analyses leads to a numerically different mass M Z , and the two mass definitions can be translated by [9] The arguments apply to M W as well. While we have used the on-shell masses M in our calculations, the numerical results in section 3 are reported in terms of the commonly used masses M . The residue R in (1) factorizes, in an excellent approximation, into initial-and final state vertex form factors, V Ze + e − µ and V Zbb ν , and Z-propagator corrections, R µν Z . For this reason, the unfolded Z-peak forward-backward asymmetry A bb,0 FB and forward-backward left-right asymmetry A bb,0 FB,LR can be written, also in an excellent approximation, as where P e is the electron polarization and The right part of (4) follows from the definition where Q b = −1/3. Technically, the calculation of A b rests on the calculation of the vertex form factor V Zbb µ , whose vector and axial-vector components can be obtained using the projection operations where D = 4 − 2 is the space-time dimension and p 1,2 are the momenta of the external b-quarks, and k = p 1 + p 2 . As a result, only scalar integrals remain after projection, but they may contain non-trivial combinations of scalar products in the numerator. More specifically, we here calculate the bosonic two-loop contribution to the (complex) ratio The determination of the pseudo-observables in Eq. (3) from true observables requires carefully written interfaces for the unfolding and subtraction of QED, QCD and box contributions, and other contributions not contained in the pseudo-observables; see Refs. [10,11]. In fact, the interfaces implemented in ZFITTER [10,12,13] have proven to be adequate for an analysis of Z-pole pseudoobservables at the O(α 2 ) level [3]. * The experimental values for A b and sin 2 θ b eff from a global fit to the LEP and SLC data are [18,19]: A challenge has been the evaluation of two-loop vertex integrals in the Minkowskian kinematic region. The vertices involve up to three additional mass scales besides s = M 2 Z , and many of them also contain ultraviolet (UV) and infrared (IR) singularities, even though the divergencies cancel in the final result. In general, it is not possible to compute all integrals analytically with available methods and tools, but instead one has to resort to numerical integration strategies. The techniques used in this work are discussed in section 2. Results for the numerical impact of the new corrections are presented in section 3, before the summary in section 4.

Strictly numerical two-loop integration techniques
The complete set of two-loop integrals required for this calculation can be divided into several categories. For the renormalization counterterms one needs two-loop self-energies with Minkowskian external momentum, p 2 = M 2 i + iε, M i = M W , M Z . In addition, there are two-loop vertex integrals with one non-vanishing external momentum squared, s = M 2 Z + iε. Two-loop self-energy integrals and vertex integrals with self-energy subloops have been computed using the dispersion relation techniques described in Refs. [3,20,21]. The remaining two-loop vertex integrals with triangle subloops amount to some 700 integrals, with tensor rank R ≤ 3, Minkowskian external momentum, and up to three dimensionless parameters per integral, from the set The aim is an accuracy of eight significant digits, to be obtained with two completely independent calculations.
A variety of integrals were calculated already for the leptonic Z boson asymmetry parameter A e . Here, up to two dimensionless parameters had to be treated [2,3,22]. In view of the larger number of scales encountered in the Zbb vertex, and also aiming at comparably simple and semi-automatic algorithms with easy re-use, a fully numerical strategy was applied here.
No reduction to a minimal set of master integrals (MIs) was attempted, except for simple cancellations of numerator and denominator terms. There are several reasons; none of them is stringent. One might perform a standard reduction to MIs, which could reduce the number of integrals by about a factor of ten. From the point of view of performance of the project as a whole, this is no important gain in efficiency, because the time of calculating the integrals is not a limiting factor. On the other hand, for cases with many different mass scales, coefficient terms in integral reductions can become very large, which makes this approach cumbersome from a technical point of view. Furthermore, using the program KIRA [23], we observed that the numerical treatment of the coefficient terms becomes difficult for some integrals with propagators of mass M Z at the kinematical point s = M 2 Z + iε. At the same time, we know that a number of MIs will remain to be evaluated numerically. The techniques developed for these can relatively easily be applied to the complete set of (unreduced) integrals, thus obviat-ing the need for integral reductions. Finally, our goal was to create a selfcontained general-purpose numerical package, see Ref. [24] for more details. Reductions and partial analytical solutions are difficult to integrate into this.
As was mentioned above, individual integrals will contain both UV and/or soft and collinear divergencies. We have employed two techniques with an automatic control of these divergencies: sector decomposition (SD) and Mellin Barnes (MB) representations.
It is essential that the numerical methods work sufficiently stable for Minkowskian kinematics. For sector decomposition [25,26] this can be achieved through a complex contour deformation of the Feynman parameter integrals, as implemented in the publicly available packages FIESTA 3 [27] since 2013 and SecDec 3 [28] since 2015. Nevertheless, we observe serious convergence problems for some of our integrals. As a second, independent method we chose the representation of Feynman integrals by Mellin-Barnes integrals [29][30][31]. The MB method has been well developed in recent years and there are useful software packages available at the MBtools webpage in the hepforge archive [32]: MB [33], MBresolve [34], AMBRE 1 [35] and barnesroutines (D. Kosower). Further, one may use PlanarityTest [36], AMBRE 2 [37] and AMBRE 3 [38], as well as MBsums [39], which are available from the AMBRE webpage [40]. For our purposes, we have derived MB representations with AMBRE and used the package MB, aided by MBresolve and barnesroutines, for a derivation of an expansion in terms of = (4 − d)/2. In particular, AMBRE 2 has been employed for planar and AMBRE 3 for non-planar topologies, using PlanarityTest for the automatic identification of the planarity status. For the numerical treatment of massive MB integrals with Minkowskian kinematics, the package MBnumerics is being developed since 2015 [41]. For the final numerical integration, it calls the CUHRE routine of the CUBA library [42,43]. Some general features of both AMBRE 3 and MBnumerics/MB have been described recently [44,45]. For cross-checks, a variety of integrals was also calculated with NICODEMOS [46] and dispersion relation techniques [3]. For one-scale integrals, we could make several comparisons with existing analytical and semi-analytical results [47][48][49][50]. A comprehensive description of our numerical package [41], including a discussion of the numerical derivation of the Z → bb integrals will be given elsewhere [24].

Using sector decomposition
For Euclidean kinematics all the needed integrals can be evaluated straightforwardly with sector decomposition, using the packages FIESTA and SecDec. One obtains Feynman parameter integrals with 4 or 5 dimensions. For Minkowskian kinematics the numerical SD method still works well for most of the integrals, although there were some problematic cases: • For 16 single-scale six-propagator integrals with one massive line and s = M 2 Z , no result at all was obtained with sector decomposition: see Fig. 1   With our implementation of the alternative Mellin-Barnes method, at least 8 significant digits were achieved for all integrals in this list, with exclusion of the last item where we obtain an accuracy of 6 digits.

Using the MBtools suite
The number of dimensions of the Mellin-Barnes integrals increases with the number of mass scales and the complexity of the integral topology. AMBRE 2.1 and AMBRE 3 find the lowest dimensionality of the MB integrals to be solved [38,44]. † The largest number of MB dimensions encountered here is eight: for the constant terms of the non-planar integrals shown in Fig. 1  highly oscillating and slowly vanishing at infinity. There is a variety of methods to improve the convergence of Minkowskian MB integrals. We mention here those which proved to be most efficient, but refer for details to the literature [24,44]: • Integrand mappings. Before applying a standard integration routine like CUHRE, we found a tangent mapping to be efficient, t i → 1/ tan(−πt i ), combined with calculating exp[ i ln(Γ i )] rather than the product Π i Γ i .
• Contour rotations. The transformation z i = x i + it i →z i = x i + (θ i + i)t i may improve the damping of oscillatory terms like [−(s + iε)/M 2 ] f (z i ) at infinity. For multi-dimensional MB integrals, one may try to perform "synchronized" rotations using a universal parameter θ i ≡ θ, in order to avoid crossing of poles by the contour change [51]. However, for singlescale integrals, which depend only on , the contour rotation will not improve the behaviour at infinity.
• Contour shifts. It proved to be extremely efficient to make use of a well-known property of Γ-function: At the negative axis between the pole positions, its value becomes smaller when the function is evaluated at an argument further way from the origin. If a pole gets crossed by an argument shift, one has to add the corresponding residue which by itself is also an integral, but will have a dimension less than the original one. Doing this several times, with several integration variables, the original MB integral gets replaced by several lower-dimensional integrals which may be easier to calculate, plus the original one with shifted integration path. The resulting smaller contribution of the original integral to the net result has the effect that its poor knowledge gets numerically less important. In effect, the procedure consists of a summing over a finite number of residues with a controlled remainder. Shifts of the integration contours were proposed first in Ref. [52].

The package
MBnumerics is yet under development. Currently, it can treat Minkowskian MB integrals with up to four dimensions and Euclidean ones with up to five dimensions with a good precision, and with more dimensions at reduced precision. The package is currently limited to integrals with few different mass scales. There is the opportunity of internal cross-checks by integral reductions with the package KIRA [23], followed by a numerical evaluation with MBnumerics.m. This procedure was used for few integrals where a second calculation was difficult, and an accuracy of at least 6 digits was reached in these cases for the second calculations. Further improvement is possible, but the result is more than sufficient for the purposes of the present calculation.
We will complete this section with a few numerical examples. Using 24 hours on the same computer, we obtained with SecDec: I 1d,SD = 1.541 + 0.2487 i + 1 (0.123615 − 1.06103 i) (10) With the MB method, we solved all the 100 integrals which depend on only one parameter, The one-scale integrals have up to four MB dimensions and were the testing ground during the development of MBnumerics [41]. One may calculate all these integrals using the results of [47][48][49][50]. Unfortunately, the authors did not provide a ready-to-use implementation for numerical evaluation. So it was more efficient for us to apply our numerical packages and to perform additional checks for some selected cases. We like to mention here two examples.
Integral (a) of Fig. 1 is a finite planar five-propagator Feynman integral. We derived a representation with several 2-dimensional MB integrals, requiring about 300 seconds for an accuracy of 14 digits with AMBRE 2/MB/MBnumerics, and got The integral is analytically known from [49] as a combination of generalized harmonic polylogarithms: The difference in sign compared to (11) is due to different metrics.
The non-planar finite integral (c) of Fig. 1 with two massive lines may be written as a fourdimensional MB integral and was solved with a series of contour shifts. The needed computer time to determine 11 digits amounts to few minutes: The integral is known from [50]. It is one of three master integrals, calculated by solving a system of differential equations (DEQ) numerically. At s/M 2 Z = 1 + iε it is: To reach an accuracy better than the 11 digits shown above with MBnumerics would require some effort, but is feasible.
The 1/ 2 poles have been verified to cancel analytically and numerically for g b V (M 2 Z )/g b A (M 2 Z ), with more than 12 digits precision. The cancellation of the 1/ poles has been checked numerically with 8 digits precision. For the finite part, we obtain a net precision of better than 7 digits, which is more than sufficient for practical purposes.
As evident from the table, the new bosonic two-loop result is about a factor of four smaller, but of similar order of magnitude, as the fermionic electroweak two-loop corrections [1]. § For varying input parameters, the new result is best expressed in terms of a simple fitting formula, This parameterization reproduces the full calculation with average and maximal deviations of 5×10 −8 and 1.2 × 10 −7 , respectively, for the input parameter ranges indicated in Tab. 1.
Combining this result with the already known corrections (see above) the currently most precise prediction for sin 2 θ b eff is obtained. Additionally, one free parameter can be eliminated by using the Standard Model prediction of M W from the Fermi constant G µ . The W -boson mass has been calculated previously including the same perturbative higher order contributions as listed above [73]. To a very good approximation, this result can be written as Here ∆α is the shift of the electromagnetic fine structure constant due to light fermion loops between the scales q 2 = 0 and M 2 Z . The best-fit numerical values for the coefficients are given by With these values, the formula (20) approximates the full result with average and maximal deviations of 2 × 10 −7 and 1.3 × 10 −6 , respectively, within the ranges in Tab. 1.

Summary
The determination of the electroweak two-loop corrections to A b and sin 2 θ b eff has been completed. We have shown by explicit calculation that the numerical result for their bosonic corrections is expectedly small compared to the presently available experimental accuracy. However, the anticipated measurements at a future accelerator of the ILC/FCC-ee/CEPC generation aim for an accuracy comparable to electroweak two-loop effects [74][75][76].
Applications related to Drell-Yan processes at the LHC are also of the single-particle resonance type and may be envisaged with the technique developed here. The numerical packages FIESTA, SecDec and the MBtools suite with the new packages AMBRE 3 and MBnumerics will be sufficient, in combination, for calculating the whole class of massive two-loop self-energy and vertex integrals in the Standard Model and beyond.