Two-loop electroweak threshold corrections in the Standard Model

We study the relationships between the basic parameters of the on-shell renormalization scheme and their counterparts in the $\overline{\mathrm{MS}}$ scheme at full order ${\cal O}(\alpha^2)$ in the Standard Model. These enter as threshold corrections the renormalization group analyses underlying, e.g., the investigation of the vacuum stability. To ensure the gauge invariance of the parameters, in particular of the $\overline{\mathrm{MS}}$ masses, we work in $R_\xi$ gauge and systematically include tadpole contributions. We also consider the gaugeless-limit approximation and compare it with the full two-loop electroweak calculation.


Introduction
The measurement of the Higgs boson mass at the Large Hadron Collider [1] not only fully confirms the validity of the Standard Model (SM) around the electroweak scale, but also opens a possibility for a precise study of the applicability of the SM at energies of the order of the Planck mass. Renormalization group (RG) equations, which determine the dependence on the renormalization scale μ of the running parameters, which are usually defined in the modified minimal-subtraction (MS) scheme of dimensional regularization, play an essential rôle in such analyses. In the SM, the corresponding RG functions, whose knowledge had been limited to the two-loop order for a very long time [2], have recently been evaluated at three loops for all running parameters, including gauge, Yukawa, and scalar self couplings [3].
The other aspect of the problem is the matching between the MS parameters and the physical observables, which gives rise to threshold corrections. 1 These not only include terms of the form ln μ 2 , but also non-logarithmic ones. The relationships between the MS and pole masses of the intermediate bosons were obtained at the two-loop level in Refs. [4,5]. As for the threshold corrections to the top and bottom quark masses and Yukawa couplings, the situation is as follows. The QCD corrections, which are dominant, are available at one [6], two [7,8], three [9], and four [10] loops. The two-loop result in the supersymmetric extension of QCD was obtained in Ref. [11]. The one-loop electroweak corrections, of order O(α), were first considered in Ref. [12]. The two-loop mixed O(αα s ) corrections were provided for the bottom quark in Ref. [13] and for the top quark in Refs. [14][15][16]. Recently, also the two-loop electroweak corrections of order O(α 2 ) have been obtained in the gaugeless-limit approximation [17]. Also the threshold corrections to the self-coupling constant of the scalar field were intensively studied in the literature. The O(α) corrections were evaluated a long time ago in Ref. [18] and the O(αα s ) ones recently in Ref. [16]. As for the O(α 2 ) corrections, the leading term was found in Ref. [19], and an interpolation formula, which also includes subleading contributions, was given in Ref. [20]. These analyses were recently revisited in Ref. [21] by providing precise numerical results.
In this paper, we systematically present the complete two-loop threshold corrections, from the orders O(α), O(αα s ), and O(α 2 ), to all the running parameters of the SM independently obtained by an analytic calculation. This includes the masses of the W , Z, and Higgs bosons (m W , m Z , m H ) and those of the top and bottom quarks (m t , m b ) as well as the gauge couplings (g, g ), the Higgs self-coupling (λ), and the top and bottom Yukawa couplings (y t , y b ). In contrast to Refs. [20][21][22], all our calculations are performed in R ξ gauge keeping the gauge-fixing parameters free, which allows us to explicitly track the ξ dependencies and so to ensure the gauge independence of the threshold corrections and the MS parameters. The tadpole diagrams turn out to play a crucial rôle in this (see Subsection 2.2). This paper is organized as follows. In Section 2, we set up the stage for our calculation of the threshold corrections. In Subsections 2.1-2.4, we discuss the various ingredients entering our analysis. Our results are presented in Section 3 and Appendix A. In Appendix B, we also list the MS renormalization constant of the Higgs boson mass.

Setup
The SM may exist in two different phases: the symmetric phase and the phase with the spontaneously broken symmetry. The phase is determined by the potential of the scalar field φ, where m φ is a mass parameter and λ is the self-coupling constant of the scalar field. While stability requires λ > 0, the term m 2 φ can be either positive (symmetric phase) or negative (broken phase). In the symmetric phase, the SM is naturally parametrized by the following set of parameters: In the broken phase, it is convenient to choose an alternative set of parameters, namely where e is the electromagnetic gauge coupling. At the tree level, the parameters in Eqs. (2) and (3) can be related to each other by for the couplings and for the masses, where we have introduced the vacuum expectation value of the scalar field v ≈ 246 GeV, which characterizes the broken phase. This parameter is not independent, since the parameter sets in Eqs. (2) and (3) each fully determine the theory. 2 In fact, we have in terms of Eq. (2) and in terms of Eq. (3). Eqs. (4)-(7) are subject to radiative corrections. However, we have the freedom to choose 16 out of the 32 parameters in Eqs. (2) and (3) to be "observables" and to consider the remaining ones to be derived parameters of the theory. At energies of the order of v, the parameters in Eq. (3) are intuitively closer to what we would call "observables," especially if we define them in the on-shell renormalization scheme. Then, these include Sommerfeld's fine-structure constant α Th = e 2 Th /(4π) as measured in Thomson scattering and the pole masses M W , M Z , M H , and M f , which are the zeroes of the inverses of the propagators of the respective particles. By contrast, the parameters in Eq. (2) are more suitable for studies of the RG evolution. In this work, we shall discuss the relationships between these parameters and the radiative corrections to Eqs. (4)- (7).
As usual, we consider Eqs. (4)-(7) to be valid through all orders of perturbation theory by definition. For the reason explained above, we define the parameters in Eq. (2), which appear on the right-hand sides of Eqs. (4) and (5), in the MS renormalization scheme, which then carries over to the parameters in Eq. (3), which appear on the left-hand sides of Eqs. (4) and (5) (see Subsections 2.3 and 2.4). The μ dependencies of the parameters in Eq. (2) are governed by the RG equations. These allow us to run the parameters from a few GeV way up to the Planck mass. The relevant β functions are available at two [2] and three loops [3]. The masses in Eq. (3) also obey RG equations. The corresponding RG functions are available at the two-loop order from Refs. [4,5,17] and are partly recovered here by an independent calculation. As explained in Refs. [4,5], we can always relate the RG functions of the unbroken and broken phases, which serves as a welcome check for the correctness of our results.
As for the initial conditions for the RG evolution, the MS parameters are usually expressed in terms of the parameters of the on-shell renormalization scheme, α Th , M W , M Z , M H , and M f . The matching scale where this is done is typically chosen to be of the order of M Z or M t . In this case, one may safely neglect the masses of all fermions, except for the one of the top quark and possibly also the one the bottom quark. We shall retain the full dependence on the bottomquark mass at one loop, but neglect it at two loops. The strong-coupling constant α s (μ) is always defined in the MS renormalization scheme. Finally, v is related to e via Eq. (7). Phenomenologically, it is more convenient to express the MS parameter v(μ) in terms of the Fermi constant G F , which is measured in low-energy processes of the weak interaction, such as muon decay, through the exact relationship where r(μ) is an appropriate variant of Sirlin's r parameter [23] (see Subsection 2.1). Inserting Eq. (8) in Eq. (5) and accommodating the relationships between the MS and pole masses, one may cast the relations of interest in the form In turn, substituting Eq. (8) in the right-hand sides of Eqs. (9)- (12) and equating the outcome with the right-hand sides of Eq. (5), we may write the MS to pole mass relationships as which are valid to all orders of perturbation theory. It is the goal of this work to evaluate r(μ) and δ x (μ) for x = W, Z, H, t, b through order O(α 2 ).

Vacuum expectation value and r(μ)
We start our discussion with the evaluation of the threshold corrections to v(μ), i.e., of r(μ) in Eq. (8). As already mentioned above, this quantity is very similar to r in the on-shell scheme [23], which contains the radiative corrections to the muon lifetime that the SM introduces on top of the electromagnetic corrections in the Fermi model of four-fermion interactions. The difference between r and r(μ) is that the former is defined in the on-shell scheme according to while the latter obeys the analogous relation with the on-shell parameters replaced by their MS counterparts, Consequently, the calculations of r and r(μ) are very similar. In particular, they are both based on the matching of the Fermi model and the SM and exhibit a factorization of the lowenergy scales. The contribution of the light (massless) fermions to r(μ) may be found in Ref. [24]. The MS parameter α(μ) deserves a detailed discussion, which will be presented in Subsection 2.4. The framework for the evaluation of r or r(μ) at any loop order was established in Refs. [25,26]. One may start from any amplitude involving the exchange of the weak charged current in the SM, e.g., A(e + ν e → μ + ν μ ), where e, ν e , μ, and ν μ are the electron, the electron neutrino, the muon, and the muon neutrino, respectively. As demonstrated in Ref. [25], there is a factorization theorem that allows for the convenient separation of the soft scales (see also the discussion in Ref. [27]). Applying this to the evaluation of r(μ), we only need to perform the renormalization procedure in the MS scheme. According to Ref. [25], we may write where Z 2,f is the wave function renormalization constant of the left-handed field component of fermion f in the MS scheme. The subscript hard in Eq. (16) implies that all external fourmomenta and all the light-fermion masses are identically put to zero before the integration over the loop momenta. 3 Obviously, such a procedure generates a lot of infrared divergences in the calculation. These are regularized by the dimensional-regularization parameter ε and cancel on the right-hand side of Eq. (16). This cancellation is nontrivial, and this is exactly the statement of the factorization theorem [25]. Following this procedure, the evaluation of r is completely reduced to vacuum diagrams with one or two loops. In Appendix A, we shall present an analytical expression for r in terms of MS couplings and pole masses.

Role of tadpoles
An important comment is in order here. In a theory with a broken gauge symmetry, it is important to take into account the tadpole diagrams. One-and two-loop tadpole diagrams contributing to the propagator of a particle are shown in Fig. 1. All such tadpole insertions must be made not only in the Feynman diagrams contributing to the counterterms, but also in all the proper Feynman diagrams to ensure the gauge independence of renormalized scattering amplitudes. In particular, the mass counterterms are gauge dependent unless the tadpole contributions are included, as has been known for a long time [28]. This is also true for the threshold corrections, as was first observed for the one-loop electroweak threshold corrections to the Yukawa [12] and Higgs self-couplings [18].
In this work, we perform all calculations at two loops in R ξ gauge with four independent gauge parameters, related to the W and Z bosons, the photon, and the gluon. We verify by explicit calculation that the threshold corrections r in Eq. (8) and δ x in Eqs. (9)-(12) are gauge independent. This serves as strong check for the correctness of our results. Tadpole contributions are singular in the limit M H → 0. The most divergent terms scale as 1/M 2 H at one loop and as 1/M 4 H at two loops. This behavior is fully inherited by r and the MS masses. However, nontrivial cancellations between r and the MS masses take place in the threshold corrections δ W , δ Z , δ t , and δ b , which render the latter finite in the limit M H → 0, i.e., also terms proportional to M 0 H ln M 2 H are canceled. As for the Yukawa couplings, this was noticed at O(α) in Ref. [12] and at O(αα s ) in Ref. [13] for the bottom quark and in Ref. [16] for the top quark. Here, we investigate the behavior for M H → 0 at O(α 2 ) and find that r and the MS masses contain terms proportional to 1/M 4 H , while δ W , δ Z , δ t , and δ b are finite in this limit. The situation is different in the Higgs sector. Already at one loop, the threshold corrections to m H and λ contain terms that diverge as 1/M 2 H for M H → 0 [18]. At two loops, however, the 1/M 4 H terms cancel, so that the leading small-M H behaviors go unchanged. This cancellation ensures that λ stays finite in the limit M H → 0, since λ ∼ M 2 H according to Eq. (11). The tadpole contributions may be quite sizable numerically. In Refs. [15,17], it was noted that one-loop electroweak correction to the MS to pole mass relationship of the top quark roughly compensates the QCD one. In particular, the tadpole contribution contains terms that are enhanced as where N c = 3 is the number of quark colors, at one loop. At two loops, such terms appear in square.

MS masses
The pole mass M of a particle is defined to be the position of the pole of its propagator, i.e., the zero of its inverse propagator. 4 For a scalar boson, such as the Higgs boson, one thus needs to solve the equation where W W,T (p 2 ) is the transverse part of the W -boson self-energy. The Z-boson case is somewhat more complicated because of the γ -Z mixing. Diagonalizing the corresponding propagator matrix, one obtains 4 If the particle is unstable, then the pole position has a complex value. The latter is usually parametrized as p 2 pole = M 2 − i M for bosons and as / p pole = M − i /2 for fermions, where M and are the real pole mass and the total decay width of the particle, respectively.
where ZZ,T (p 2 ), γ Z,T (p 2 ), and γ γ,T (p 2 ) are the respective transverse self-energy functions. Finally, for a fermion f , one has where f (/ p) is the self-energy of f . 5 The solutions of Eqs. (17)-(20) may be found order by order in perturbation theory, by substituting the ansätze (19) and Eq. (20), respectively.
Alternatively, we may perform the mass renormalization in the MS scheme. The relation between the MS mass m(μ) and m 0 has the simple form where the expansion parameter ε = (4 − d)/2 measures the deviation of the space-time dimension d from 4. The mass renormalization constants Z (j ) have double expansions in the weak and strong gauge couplings, g(μ) and g s (μ), respectively. For the purposes of our two-loop analysis, we need to include the following terms As already mentioned in Subsection 2.2, all the relevant mass renormalization constants are gauge independent upon the inclusion of the tadpole contributions. Explicit expressions through order O(α 2 ) may be found for the W and Z bosons in Refs. [4,5], 6 for the top and bottom quarks in Ref. [17], and for the Higgs boson in Appendix B of this paper.
In this paper, we shall evaluate the MS masses m x (μ) for x = W, Z, H, t, b through order O(α 2 ). Besides gauge independence, also the μ dependence, which is dictated by the RG, provides a strong check for the correctness of our results. In fact, the full μ dependence of m x (μ) at two loops may be retrieved from the one-loop result for m x (μ) and its two-loop anomalous dimension. To simplify the discussion in the remainder of this section, we shall only allow for one coupling constant, a(μ). The generalization to the case under consideration in this paper is straightforward.
Adopting standard notations, the generic RG equations for a(μ) and m(μ) read The β and γ functions in Eq. (23) have perturbative expansions of the forms β = a 2 β 1 + a 3 β 2 + · · · and γ = aγ 1 + a 2 γ 2 + · · ·. In the MS scheme, β j are just numbers, while γ j are, in general, 5 Here, fermion mixing is neglected. The renormalization for mixed systems of spin-1/2 fermions was elaborated in Ref. [29]. 6 The expressions for Z (1,2) αα s of the W and Z bosons in Eq. (4.41) of Ref. [5] contain several misprints, which are corrected in Footnote 9 of Ref. [16]. polynomials in squared-mass ratios. Expressing in turn the MS masses in terms of pole ones, we find the expression for γ to assume the form where i and i are independent of μ. On the other hand, the MS to pole mass relationship has the form where The coefficients A 1 , A 2 , and A 2 in front of ln n μ 2 with n = 1, 2, . . . may be derived from the RG. Substituting Eqs. (24) and (25) in Eq. (23) and comparing the coefficients of a m ln n μ 2 , we obtain Consequently, we may independently evaluate the coefficients of ln n μ 2 with n = 1, 2, . . . from one-loop results and the mass anomalous dimension γ . In the MS scheme, the latter is completely determined by the single 1/ε pole in the corresponding renormalization constant in Eq. (21). Taking into account both the weak and strong coupling constants, we have On the other hand, γ may also be related to the anomalous dimension γ G F of the Fermi constant G F and the RG functions in the unbroken phase of the SM [4,5,30].

MS fine-structure constant
In this section, we discuss different definitions of Sommerfeld's fine-structure constant α = e 2 /(4π) and their relationships to the Fermi constant G F .
The relationship between the MS quantity α(μ) and the on-shell quantity α Th defined in the Thompson limit is usually written as [31] where, to first approximation, α(μ) is expressed through the vacuum polarization function of the photon γ (0). The latter quantity is known to be subject to nonperturbative effects due to the hadronic content of the photon. The incorporation of the perturbative O(αα s ) corrections was explained in Ref. [32]. The purely leptonic contributions are known through four loops [33]. A systematic discussion of higher-order corrections may be found in Ref. [34]. Through order O(αα s ), we have 7 where α (5) had (M Z ) = 0.02772(10) [35] is the hadronic contribution to γ (0) and L x = ln(M 2 x /μ 2 ) with M x being pole masses. The four terms inside the curly brackets in Eq. (29) are the bosonic, leptonic, top-quark, and perturbative light-quark contributions, respectively. To the same accuracy as in Eq. (29), we may also write the μ dependence of α(μ) as Using The QED β function is available through five loops including QCD corrections [36].
Eqs. (28) and (29) relate α(M Z ), which belongs to the core of electroweak physics, to the lowenergy-QED parameter α Th , which is insensitive to the weak interaction due to the decoupling of the heavy particles. While this approach allows one to conveniently implement the RG and so to consistently accommodate all the QCD contributions of the orders O(αα n s ), there are obvious drawbacks if one wishes to go beyond this approximation. First, it is less appropriate for the incorporation of electroweak physics, which comes into play by the exchange of W and Z bosons in the loops. Second, the definition of α(μ) from the QED two-point function alone has only restricted meaning. To relate α(μ) to electroweak observables beyond the one-photon-exchange approximation, one has to include more complicated objects.
A very simple consistent definition of α(μ) that avoids these drawbacks is to take Eq. (4) at face value and to install G F as an input parameter via Eqs. (9) and (10), viz.
where δ W (μ) and δ Z (μ) are given in Eqs. (9) and (10), respectively. This definition of α(μ) is exact and gauge independent. Since δ W (μ) and δ Z (μ) themselves depend on α(μ), Eq. (31) provides an implicit definition of α(μ), which may be solved iteratively, e.g., using the Newton-Raphson method. Working through two loops, we thus obtain numerically at μ = M Z where the first number on the right-hand side is the tree-level value with analytic coefficients X ij (μ). Numerically, this gives The difference between the values in Eqs. (32) and (34) may be interpreted as the theoretical uncertainty due to higher-order effects.

Results and discussion
We are now in a position to present our numerical analysis. For this, we adopt the following values of the input parameters from the Particle Data Group [35]:  [37] to the value of α (5) s (M Z ) in Eq. (35), we obtain α (6) s (M t ) = 0.1081.

Threshold corrections to the couplings and masses
We first present the threshold corrections δ x with x = W, Z, H, t , which we parametrize as follows The coefficients Y 1,0 x , Y 1,1 x , and Y 2,0 x may be evaluated numerically using a computer program to be published in a forthcoming paper [38]. In Appendix A, we list all the O(α) and O(αα s ) coefficients, Y 1,0 x and Y 1,1 x , in their full analytic forms and the O(α 2 ) coefficients, Y 2,0 x , in the gaugeless-limit approximation. For the reader's convenience, we present linear interpolation formulae for the two-loop coefficients, Y 1,1 x and Y 2,0 x , at the two most important matching scales, We now present the MS to pole mass relationships, in forms similar to Eq. (36). Specifically, we write for the bosons B = W, Z, H and for the fermion f = t . The coefficients X 1,0 x , X 1,1 x , and X 2,0 x may be evaluated numerically using the computer program to be released [38]. Convenient interpolation formulae for the two-loop coefficients, X 1,1 x and X 2,0 x , at μ = M Z read  Comparing Eqs. (37) and (41), we observe that the coefficients X 1,1 W,Z,t and X 2,0 W,Z,t are much larger than their counterparts Y 1,1 W,Z,t and Y 2,0 W,Z,t , typically by an order of magnitude. This is due to the tadpole contributions that dominate the former coefficients, but largely cancel in the latter, as already discussed in Subsection 2.2. In the case of the Higgs boson, however, the situation is quite different. In fact, X (1,1) H and X (2, Table 1 of Ref. [17]. Here, we improve this analysis by including the full O(α 2 ) contribution and updating the input parameters as specified in Eq. (35). Our new results are presented in Table 1. We observe that the gaugeless-limit approximation successfully predicts the true sign of the O(α 2 ) contribution, but significantly underestimates its magnitude by accounting for only about one half of it. The seemingly poor quality of this approximation may be partly ascribed to the presence of a heavy particle on the external lines of the two-point diagrams involved. On the other hand, the absolute deviation by itself may be considered acceptable.
In Table 2, the analysis of Table 1 is repeated for δ t (M t ). As expected by reason of the tadpole cancellation, the QCD corrections are by far dominant. In fact, the electroweak corrections only make up approximately 2% of δ t (M t ). Furthermore, we observe that the electroweak perturbative expansion exhibits a useful convergence behavior. In contrast to the case of the m t (M t ) − M t shift, the gaugeless-limit approximation works well, overestimating the true O(α 2 ) correction by about 15%.
In Table 3, the corresponding results are presented for λ(M t ). Of course, there are no pure QCD contributions in this case. As already observed in Ref. [16], the O(αα s ) contribution almost doubles the O(α) one. The O(α 2 ) contribution reaches about one quarter of the O(α) one. The gaugeless-limit approximation works here even better than in the case of δ t (M t ), with a deviation of about 10%.
Finally, we study the quantities related to the bottom quark, at μ = M b . For the MS to pole mass relationship, we obtain In Eq. (42), the electroweak corrections are overwhelming and do not exhibit a useful convergence behavior, again because of the uncanceled tadpole contributions. The latter render the electroweak extension of the MS mass definition for the bottom quark quite unfeasible in practice, while they do not affect its pole mass. This situation is unfamiliar from pure QCD, which is tadpole free. Here, the MS mass is frequently preferred to the pole mass because it is only sensitive to short-distance effects, while the latter suffers from a renormalon ambiguity. This is reflected by the relatively large size of the first term on the right-hand side of Eq. (42) and the slow convergence of the O(α n s ) corrections with n = 1, 2, 3 which build it up. We also learn from Eq. (42) that the gaugeless-limit approximation works here remarkably well, within an error of about 3%. On the other hand, the threshold corrections to the Yukawa coupling, are perturbatively stable as for the tadpole contributions. Also here, the gaugeless-limit approximation works at the 3% level.

Conclusions
In this paper, we analytically evaluated, through two loops in the SM, the threshold corrections to the electroweak gauge couplings, the top and bottom Yukawa couplings, the quartic self-coupling, and the vacuum expectation value of the scalar field as well as the MS to pole mass shifts of the top and bottom quarks. Specifically, we included the corrections of orders O(α), O(αα s ), and O(α 2 ). We emphasized the importance of the tadpole contributions to render these corrections gauge independent. Besides comparisons with the literature, UV finiteness and RG invariance served as further checks of the correctness of our results.
The threshold corrections to the gauge and Yukawa couplings are finite in the limit of vanishing Higgs-boson mass M H due to cancellations of leading tadpole contributions, and their perturbative expansions exhibit useful convergence. The threshold correction δ H (μ) to the quartic scalar self-coupling λ(μ) scales as 1/M 2 H for M H → 0, but λ(μ) is finite in this limit. Also in this case, perturbative stability is intact. All these threshold corrections are central ingredients for RG analyses within the SM and, in particular, for the determination of the M H lower bound from the requirement of vacuum stability at the scale of the Planck mass [16,19,20]. By contrast, the MS to pole mass shifts, which do not enter such RG analyses, suffer from sizable tadpole contributions, which are particularly severe in the case of the bottom quark. As a way out of this problem, it was proposed in Refs. [15,17] to define the running fermion masses in terms of the MS Yukawa couplings, as m Y,f (μ) In this appendix, we present analytical results for the threshold corrections r and δ x with x = W, Z, H, t, b defined in Eqs. (8)- (12). The MS to pole mass relationships for these bosons and fermions then follow according to Eq. (13).
Let us first introduce the notations. We work in dimensional regularization with d = 4 − 2ε space-time dimensions and 't Hooft mass μ. N c = 3 is the number of quark colors, C F = (N 2 c − 1)/(2N c ) = 4/3, and n G = 3 is the number of fermion generations. We express our results in terms of the pole masses and use the abbreviations As is well known, any one-loop two-point Feynman integral can be expressed in terms of two types of master integrals, namely where is the integral measure in d-dimensional Minkowski space-time, u j are the squares of the masses of the internal propagators, and q is the external four-momentum. To suppress the appearance of Euler's constant γ E in the Laurent expansions in ε, we pull out the factor e −γ E ε for each loop momentum.
In the following, we shall need the expansions of the Feynman integrals in Eq. (45) in ε through order O(ε). They read where The sign conventions in Eq. (47) have been chosen in accordance with the program library TSIL [39]. Some care must be exercised when a complex-valued function B 0 gets squared or multiplied by some other complex-valued B 0 function. In particular, one needs to distinguish [Re B 0 (. . .)] 2 from Re B 2 0 (. . .). In such cases, we explicitly indicate when the real part should be taken. In particular, we have if t > 0, and if t > w > 0. The former situation occurs, e.g., for δ H and δ Z and the latter for δ t . At two loops, the most general form of the self-energy integral is given by J a 1 a 2 a 3 a 4 a 5 (q 2 ; u 1 , u 2 , u 3 , u 4 , u 5 ) where a j indicate the powers of the respective propagators. The numerator of a Feynman integral may be incorporated in the representation of Eq. (54) by allowing for some of the indices a j to be negative. Alternatively, one may reduce any tensor integral to scalar ones in higher (shifted) space-time dimensions [40]. Then, the recurrence relations of Ref. [41] may be used to shift the indices as well as the space-time dimension to the appropriate values. After the reduction, any two-loop two-point Feynman integral may be expressed as a linear combination of a finite set of master integrals with the coefficients being rational functions of q 2 , u j , and d. However, the choice of the master integrals is not unique. In particular, our set of master integrals differs from the one in Ref. [41]. A drawback of the set in Ref. [41] is that the coefficients in front of master integrals may contain poles in ε. When this happens, then deeper ε expansions of the master integrals are required. 8 Furthermore, the master integrals themselves may posses infrared divergences besides the ultraviolet ones. It is possible, however, to choose the master integrals in such a way that: (a) the coefficients in front of the master integrals always have smooth limits → 0 and (b) the master integrals themselves are infrared finite. Thanks to property (a), we do not need to expand master integrals beyond the order O(ε 0 ).
We keep our notations for the master integrals very close to those of Ref. [39]. Specifically, the different integrals in Eq. (54) are denoted by the bold-faced letters I, S, T, U, V, and M. They always have indices a j = 0, 1, except for T and V, which have a j = 2 one time. These functions generally have ultraviolet poles in ε. We denote the O(ε 0 ) terms of these master integrals by I 0 , S 0 , T 0 , U 0 , V 0 , and M 0 , respectively. 9 The arguments of all these functions are written as in Ref. [39], except that we always include q 2 as the first argument if it is present. As in Ref. [39], we exclude μ 2 from the argument lists.
Since the two-loop functions never appear in products with complex coefficients, we impose the rule that their imaginary parts are to be discarded.
In Ref. [17], the threshold corrections r, δ t , and δ b were calculated at O(α 2 ) in the gaugeless-limit approximation and expressed in terms of the functions H(x) and (x) defined, respectively, in Eqs. (41) and (42) therein. The latter are related to the function I 0 introduced above with specific arguments, namely Similarly, the functions I 1 and I 2 defined in Eq. (18) of Ref. [17] are related to the function B 0 defined in Eq. (50) above, as In the remainder of this appendix, we list our analytic results for the threshold corrections r and δ x with x = W, Z, H, t, b through two electroweak and three QCD loops. The pure QCD corrections, of orders O(α n s ) with n = 1, 2, 3, only arise for δ t and δ b . In the case of δ t , the bottom quark is treated as massless. In the case of δ b , the top quark is decoupled. The O(α) and O(αα s ) corrections are exact, except that the light-fermion masses are neglected. In the O(αα s ) corrections, also M b = 0 is put, except in A 0 (b). Our exact formulae for the O(α 2 ) corrections are too lengthy to be presented here, so that we resort to the gaugeless-limit approximation. The corresponding results for δ w and δ z vanish, which provides a welcome check for our calculation. We are thus left with those for r, δ h , δ t , and δ b . Our master formula reads where δ = r, δ x and the coefficients x 1,0 , x 1,1 , and x 2,0 correspondingly carry the subscripts r and x = W, Z, H, t, b.

A.1. δ QCD
The pure QCD correction in Eq. (58) is given by where n l is the number of massless quarks, n h = 1 is the number of heavy quarks with pole mass M, l μM = ln(μ 2 /M 2 ), and a 4 = Li 4 (1/2) ≈ 0.517479.