The Lightest Higgs Boson Mass of the MSSM at Three-Loop Accuracy

In this article we provide explicit computations of the quantum corrections to the lightest CP-even Higgs Boson mass in the real version of the Minimal Supersymmetric Standard Model. We compute the contributions coming from the SUSY-QCD sector with a precision of three loops at order $O(\alpha_t\alpha_s^2)$. The renormalization procedure adopted is based on the dimensional reduction scheme in order to preserve supersymmetry to all perturbative orders. The calculation extends the region of validity of previous studies to the whole supersymmetric parameter space. The involved master integrals have been computed by exploiting new proposed approaches based on the dispersion relation techniques for the numerical calculation of three-loop vacuum integrals with arbitrary mass scales.


Introduction
Since the discovery of the Higgs boson particle at the LHC and until the date there are no conclusive observations of new physics beyond the Standard Model (BSM). It seems mandatory the development of new initiatives focused on the search of data that prove the existence of new physical degrees of freedom. This is the main aim of projects at future colliders, like the International Linear Collider (ILC) [1] and Future Circular Collider (FCC) [2]. Even at 250 GeV with a total integrated luminosity of 2 ab −1 there are proposals for the e + e − ILC collider to precisely measure the couplings of the 125 GeV Higgs boson to vector bosons, quarks and leptons [3] with an accuracy of O(1%) compared the 0.2% accuracy based on the SM predicted couplings in terms of the Higgs boson mass. That augmented precision would allow to detect the small deviations of the typical BSM scenarios. In order for the theory to meet experiments a control of perturbation theory beyond two-loops level is required. The possible methods and strategies to match the experimental accuracy expected at the FCC-ee were recently discussed in [4]. For instance the FCC-ee Tera Z for Electroweak Pseudo Observables would require three-loop electroweak calculations complemented with QCD four-loop terms [5]. Quantum corrections must reach this level of accuracy in particular for the Higgs boson mass M h because it affects in a dominant way the value of all the observables involving the interactions of the Higgs. A very precise computation of M h including additional interactions coming from BSM scenarios could give us valuable hints about the existence of new physics. There are good reasons to prefer the Minimal Supersymmetric Standard Model (MSSM) as the BSM scenario which can describe these new interactions. According with the currently accepted analysis of the Higgs boson properties in the SM [6,7,8,9], the assumption that no new physics appears up to the Planck scale (Λ P ) renders the Higgs effective potential meta-stable and forces to accept an unnatural high amount of fine-tuning (10 34 ) for the squared Higgs boson mass at the electro weak scale (Λ EW ) leaving the hierarchy problem without solution. The MSSM can provide for a solution to this unusual fine-tuning by the existence of additional bosonic contributions that cancel the unwanted quadratic divergences. Moreover similar contributions render the effective potential stable, thus the MSSM also cures the vacuum stability problem. It is worth to mention that MSSM also provides a dark matter candidate, a mechanism to explain the neutrino oscillations and a platform to include gravitational interactions.
In most of the benchmark scenarios for MSSM Higgs boson searches, the Higgs boson found at LHC corresponds to the lightest CP-even Higgs boson with a mass M h which is not a free input parameter but it is predicted in the MSSM. The upper bound on its predicted mass at leading order (LO) is given by the Z gauge boson mass, M Z = 91.2 GeV, leading to the exclusion of the MSSM at current collider experiments. However, higher order quantum corrections to M h lead to a considerably large shift on its central value, ∆M h ≈ 40 GeV, making the MSSM Higgs sector compatible with the mass and the detected production rates of the LHC Higgs-like signal over a wide range of the parameter space of the relevant phenomenology scenarios [10]. The state of art of the corrections to the Higgs boson mass in the MSSM is quite advanced and widely studied. At one-loop level the full corrections are found in [11] with real parameters. The detailed results of a Feynman diagrammatic (FD) calculation of the leading two-loop QCD corrections at order O(αα s ) can be found in [12], in particular the O(α b α s ) [13] and O(α t α s ) [14] contributions using the FD approach are known in the limit where the external momentum vanishes and in the MSSM version with complex parameters. In this limit there is an alternative procedure to compute the above corrections, the Effective Potential (EP) approach. A comparison of the corresponding two-loop results in the FD and EP approaches at O(αα s ) can be found in [15,16] and references therein. In contrast to the EP method, the FD approach has the advantage that it can allow for non-vanishing external momentum. An evaluation of the momentum dependence of the two-loop corrections, including all the terms involving the QCD couplings, in the modified dimensional reduction scheme (DR) was presented in [17]. The latest status of the momentum-dependent two-loop corrections was discussed recently in [18,19] using a hybrid on-shell-DR scheme and including corrections of O(p 2 α t α s ) for the real version of the MSSM. A complete two-loop QCD contributions to M h in the MSSM with complex parameters including the full dependence on the external momentum can be found in [20]. Finally, at three-loop level there is a first diagrammatic computation performed by P. Kant and collaborators [21,22,23], where the radiative corrections to M h were computed in the SUSY-QCD sector including non-logarithmic terms of order O(M 2 t α t α 2 s ). They have exploited the methods of asymptotic expansion in order to provide precise approximations in the relevant mass hierarchies. Their results were written in the publicly available code H3m and then implemented into the C++ module Himalaya [24,25] linked to the Mathematica generator FlexibleSUSY [26,27] in a pure DR context. With our work we provide an alternative calculation of the three-loop corrections to the lightest Higgs boson mass in the SUSY-QCD sector of the real version of the MSSM (rMSSM) at order O(α t α 2 s ) [28]. We have followed the Feynman diagrammatic (FD) procedure to obtain a renormalized correction in the DR scheme. Taking in mind that in the MSSM is no clear a priori what are the hierarchies of the masses, we have avoided the application of asymptotic expansions at the integral level and we have obtained a correction in terms of a set of three-loop master integrals whose numerical evaluation is possible for an arbitrary mass hierarchy thanks to the development of new calculation techniques based on the dispersion method [29]. Our results are presented in this paper as follows. Section 2 contains a description of the renormalization procedure to derive the quantum corrections to M h coming from the Higgs sector of the rMSSM. The technical details of our computations are discussed in Section 3. In Section 4 we present a numerical analysis where the effects of the three-loop corrections on the pole mass M h are evaluated at some kinematic limits. We give our conclusions in Section 5.
We have used the short notation s θ = sin(θ) and c θ = cos(θ). At tree level the Higgs sector can therefore be parametrized in terms of the Z gauge boson mass (M Z 0 ), the mass of the CP-odd Higgs boson (M A 0 ) and tanβ, which is the ratio of the two vacuum expectation values. The masses of the five Higgs bosons particles h, H, A and H ± follow as predictions.
The dominant contribution of the three-loop corrections to M h comes from the QCD sector of the MSSM Lagrangian. Having in mind that, we have derived the necessary renormalization conditions by imposing the following procedure. i) We have computed the correction in the EW gaugeless limit at order O(α t α 2 s ). Consequently, all the terms proportional to M Z 0 are disregarded. In this limit the tree-level relation α = β is satisfied. We have also used the non-light fermions limit (NLF) where all the fermion masses are put to zero except the mass of the top quark, M t .
ii) The numerical dominant contributions due to the higher order corrections to the Higgs pole masses can be obtained in the limit of vanishing external momentum, p 2 = 0, where p is the momentum transferred in the external lines of the self-energy corrections. Up to two-loop level the shifts in the mass M h due to the momentum dependence are below 1 GeV in all scenarios studied, the bulk of the two-loop corrections comes from the effective-potential effects which are of the order of 10 GeV [18]. A well-behaved perturbative expansion will make this dependence even weaker at three-loop level. Motivated by these observations we have adopted the approximation of zero external momentum and therefore we have avoided dealing with the Higgs wave function renormalization and also with the renormalization of tanβ. As the wave functions renormalization of the Higgs fields don't make contributions, it is enough to rewrite (5) in terms of the renormalized parameters M A and T j according to iii) We have defined the vevs of the Higgs fields (v 1 , v 2 ) as the minima of the full effective potential. The condition Ω H 0 1,2 Ω = 0 is satisfied order by order in the perturbative expansion of the one-point Green function implying where T (l) j is the l−loop Higgs tadpole contribution. Taking in consideration the above restrictions the expressions for the renormalized three-loop corrections to the neutral CP-even Higgs bosons are reduced to (3) The letter sigma with a hat represents the renormalized Higgs self-energies at zero momenta, the letter sigma without a hat is the unrenormalized Higgs self-energies and the terms with delta are the counter-terms of the CP-even Higgs boson masses whose expressions are: where h t is the top Yukawa coupling, We have imposed an on-shell renormalization condition for the renormalized self-energy of the neutral CP-odd Higgs boson. Finally, the renormalized CP-even Higgs boson masses are obtained by determining the poles of the inverse propagator matrix which is equivalent to solve the determinantal equation The three-loop counter-terms derived from the mass Lagrangian (5) in the above expressions are useful to cancel local UV divergences. However, the unrenormalized self-energies as well as the tadpoles of the neutral CP-even and CP-odd Higgs boson fields can contain also non-local divergences coming from a sub-loop in the three-loop diagrams. It is therefore necessary an additional sub-renormalization procedure to remove these infinities. The procedure consists in the inclusion of one and two-loop vacuum diagrams with counter-term insertions which cancel those divergences arising in a sub-loop. At order O(α t α 2 s ) we need the O(α s )-contribution of the one-loop counterterms coming from the renormalization of the gluino mass, the top quark mass, the squark masses and the stop mixing angles. In addition, we need the two-loop renormalization of the top mass, the stop masses and stop mixing angles at order O(α 2 s ). We have got all the one and two-loop counterterms in the DR scheme, where the UV divergences are minimally subtracted. In order to preserve supersymmetry to all perturbative orders we have used the regularization procedure DRED [32]. In Appendix B we describe the main renormalization conditions to remove the sub-divergences from the three-loop diagrams and we give the expressions of the counter-terms.

Technical Details to the Three-Loop Calculation of M h
In this section we are going to discuss the technical details of our three-loop diagrammatic computation. We restrict the calculations to the SUSY-QCD sector where the Higgs couples just to the top quark or its super-partner. That is, we are interested in the terms of order α 2 s α t in the perturbative expansion of the Higgs mass. The three-loop radiative correction to M h is obtained by evaluating the neutral Higgs boson self-energies ( ) and the tadpoles (T ) for the fields h, H U4(m1, m2, m3, m4) and A according with the equations (10) and (11). All the needed Feynman diagrams and their corresponding amplitudes are generated with the package FeynArts [33]. FeynArts generates the amplitudes without taking the contractions of the Dirac and Color indices. As a result, we obtain a total set of 7738 self-energies and 7180 tadpoles which are moreover not regularized. With the help of the Mathematica package FeynCalc [34] we have written a routine that implements the regularization of Feynman integrals by dimensional reduction. In this scheme all gamma-matrices are of the quasi-four-dimensional space (Q4S) [35,36,37] while loop-momenta remain in the quasi-D-dimensional space. Thus, we have performed the Dirac and Lorentz algebra over the numerators of the integrals using the Q4S algebra. It is worth to mention that the Q4S algebra for diagrams which contain the cubic vertex quark-squark-gluino requires a special treatment. This vertex introduces the γ 5 matrix in the Dirac traces. We have dealt γ 5 as an anti-commuting object, which is allowed in our computation because traces with an odd number of γ 5 always contain less than four gamma matrices [38]. By other side, sum over the Color indices was done with the help of the package SUNSimplify of FeynCalc. After performing the Dirac and the Color algebra, each amplitude can be expressed as a superposition of a set of "scalar integrals" which can contain propagators with negative powers (irreducible numerators). We need in our calculation a total of 3525 of such integrals. We have avoided the use of asymptotic expansions at the integral level in order to obtain a general expression valid for any election of the input masses. For this reason we have exploited the fact that this set of integrals are not independent of each other but related by the integration by parts (IBP) and Lorentz invariance (LI) relations. We have used the IBPs to generate a homogeneous system of linear equations where the scalar integrals are the unknowns.
The system can be reduced to a small set of irreducible integrals, the so called Master Integrals. This is something that cannot be done by hand because there are thousand of equations. Thus, we have used the program Reduze [39], an implementation of the Laporta algorithm, in order to carry out this reduction. We have found a basis of 32 master integrals which correspond to different mass configurations for the five vacuum diagrams depicted in Figure 1, where each topology can contain at most four independent mass scales. The one-loop function A0 and the two-loop function T 3 have very well known analytical expressions for any configuration of masses, while the three-loop functions U 4, U 5 and U 6 have an analytical solution just for the cases where there are one or two independent mass scales but when there are three or four independent scales a numerical evaluation is required. The analysis of convergence and the numerical evaluation of the three-loop master integrals which are unknown analytically were performed with the help of the program TVID developed by A. Freitas [40]. TVID uses the discontinuities coming from the one-loop self-energy (B0) and the one-loop vertex (C0) to produce dispersion relations that are useful in the evaluation of the three-loop vacuum integrals in terms of one and two -dimensional integral representations. In all cases we have been able to get analytical expressions for the divergent part while the finite part could be numerically evaluated with up to 10 digits of accuracy. It is possible to reach this precision because the numerical integrations are at most 2-dimensional and therefore there is a controlled treatment of any singularities (see Appendix A). To sum up, each three-loop amplitude can be expressed as a linear combination of the integrals drawn in Figure 1 with coefficients that are ratios of polynomials in the masses and the spacetime dimension. These coefficients can contain poles of first and second order and therefore the renormalized correction to M h requires also the evaluation of the evanescent terms of the master integrals up to second order, that is to say, the terms at order O(ε) and O(ε 2 ) in the Laurent expansion. It is no possible to evaluate these contributions with TVID because the higher-ε terms of the real and imaginary parts of B0 and C0 and therefore of U 4, U 5 and U 6 are not included in the program. For this reason we have used the code SecDec [41] which admit a numerical evaluation of the evanescent terms. By other side, we also need to generate the amplitudes for the diagrams which are responsible for removing the non-local sub-divergences. We have written a routine in Mathematica that generates all the expressions for the counter-terms listed in Appendix B. The generation of the involved regularized amplitudes was done with the help of the FeynArts and FeynCalc functions. In contrast to the three-loop diagrams, the one and two-loop counter-terms are determined from the evaluation of fermionic and scalar self-energies with the external momentum transferred different from zero, p 2 = 0. The resulting self-energies can be further reduced using the Tarasov method [42], that is implemented in the code TARCER [43], to the basis of one and two -loop master integrals represented by the diagrams of Figure 2. The one-loop counter-terms of the top quark mass, the gluino mass, the sfermion masses and the stop mixing angles can be determined in terms of the Passarino-Veltman functions A0 and B0. The two-loop counter-terms of the stop masses and mixing angles can be expressed as a superposition of the eleven master integrals depicted in Figure 2.
In the DR scheme we just need the divergent parts of the integrals. The one-loop functions A0 and B0 have the analytical expressions found in Appendix C. The functions Dm 1 m 2 J, Dm 1 m 2 m 3 J and D2m 1 m 2 J are finite. D2m 1 J has only a 1/ε divergence with coefficient −1/2m 2 1 , independent of the values of m 2 and m 3 . Dm 1 J has poles of one and second order, its 1/ε 2 divergence is 1/2, mass independent while its 1/ε divergence is 1/2 − ln(m 2 1 /µ 2 r ). The 1/ε 2 divergence coefficient of J is 1/2(m 2 1 + m 2 2 + m 2 3 ) while the 1/ε divergence coefficient is For the cases where one or more masses m 2 j vanish in the expression (17) one should take the zero order term of the Taylor expansion around m 2 j = 0 to get the right expression. The same apply for T 3 but choosing the external momentum p 2 equal to zero. Finally, we need the functions F and V  have a divergence 1/ε 2 with coefficient 1/2 and a 1/ε divergence equal to B0 f in p 2 , m 2 2 , m 2 4 + 1/2, where B0 f in refers to the finite part of the function B0. All these analytical expressions were numerically checked with the code SecDec, where the prefactor Exp[(−2γ E ε)] must be specified in order to get the correct result. In Appendix C the counter-term expressions involved in the renormalization of the non-local ultraviolet divergences are listed. These results can be checked with those of the review [44] and references therein.

Results and Numerical Analysis
Once the local and non-local UV divergences have been subtracted from the Higgs self-energies, we get finite three-loop corrections to the CP-even and CP-odd Higgs boson masses that are useful in the derivation of the renormalized value of M h which is obtained as a solution of the pole equation (16). Our corrections depend upon 26 parameters: the renormalization scale µ r , the SM DR parameters h t , M t , α s and the MSSM parameters µ, tanβ, M A , Mg, θ t ,m f 1,2 and A f , with f = u, d, t, b, c, s. Their values as well as the renormalization group evolution of the SM parameters are determined with the help of the spectrum generator SoftSUSY [55]. We use the package SLAM [56] (Supersymmetry Les Houches Accord with Mathematica) in order to export to Mathematica any needed parameter generated with SoftSUSY. We have performed a numerical comparison between our three-loop predictions and the other higher order corrections currently included in FeynHiggs [45,46] and the three-loop results implemented in H3m [21,22] combined with the lower-order results of FeynHiggs. We discuss our results in three different limits: i) the m f ree h scenario, where µ r , tanβ, M A and A t are left as free input parameters.
ii) The m max h and m mod+ h scenarios analyzed in [10]. In these three different scenarios we don't make any specific assumptions about the soft SUSY-breaking mechanism and we interpret the LHC signal at 125 GeV as the lightest CP-even Higgs boson. We consider values of the SUSY scale in the region M SU SY < 1.2 TeV, where the combined theoretical uncertainty of the fixed-order calculation is lesser than the combined uncertainty of the effective field theory (EFT) calculation [47]. At the critical point M SU SY = 1.2 TeV the fixed-order and EFT combined uncertainties are equal and a hybrid calculation should be used [48,49,50]. Above the critical point the EFT computation  is more accurate and therefore an EFT approach, where an effective SM is used below a supersymmetric scale [51,52,53], should be preferred. In this paper we discuss our numerical results in scenarios where a fixed-order calculation is recommended.
In the m f ree h scenario we fix the parameters Mt L = Mt R =m q 1,2 = M SU SY = 1 TeV, where q denotes any quark different than the top quark. We also set Mg = 1500 GeV, µ = 200 GeV and M A = 1000 GeV. Using this limit we have studied the dependence of the Higgs boson mass M h on the soft-breaking parameter A t (Fig-3(a)) on the input parameter tanβ (Fig-3(b)) and on the renormalization scale µ r (Fig-4(a) and Fig-4(b)). At one and two -loop level we have generated the Higgs mass predictions with the help of the code FeynHiggs. These contributions are represented in the plots with the dot-dashed and dashed curves respectively. This convention has been used in all panels. We have included the full one and two-loop corrections (α s α t , α s α b , α 2 t , α t α b , α 2 b ) in the rMSSM. The one-loop field-renormalization constants and the one-loop tanβ counter-term are set in the DR scheme. We don't assume any approximation for the external momentum value for the one and two-loop corrections, i.e. we set in FeynHiggs a full determination of the propagator matrix poles (15). Besides, the dotted curve represents the three-loop predictions at order O(M 2 t α t α 2 s ) coming from the program H3m while the red solid line represents our three-loop prediction evaluated at the same order. Our three-loop results shown in Figure 3 are quite sizeable, amounting a size between 0.8 to 3.1 GeV compared to the two-loop corrections and −0.373 to 0.418 GeV regarding the three-loop prediction of H3m. The relative size and sign of the corrections depend on our election of the renormalization scheme. Figure 3(b) shows a strong dependence on tanβ for small values close to tanβ = 3, while for large values above tanβ = 10 the variation of M h is marginal and closer to the LHC Higgs mass value. Figure 4 depicts the dependence of M h on the renormalization scale µ r . In Fig-4(a) the dependence is studied for three different values of tanβ, namely tanβ = 3, 10 and 20. The three-loop corrections lead to a more stable dependence of M h with the renormalization scale µ r than the one and two -loop predictions, reducing the scale dependence by a factor between 1.5 and 2.0. This stability increases for lower values of tanβ in the m f ree h scenario. Fig-4(b) shows the RG evolution for three values of soft breaking parameter A t , A t = 0.3 TeV, 1 TeV and 1.7 TeV. The evolution is more  Table 1 for each scenario. The stop mixing angle θ t and the stop massesm t 1 , m t 2 are functions of the parameters specified in Table 1.  [47], shows a combined uncertainty for a maximal stop mixing scenario who varies between 1-4 GeV depending on the SUSY scale. Our results have then very similar numerical values to the ones obtained with H3m with differences which are within the combined theoretical uncertainty.

Conclusions
Following the Feynman diagrammatic procedure we have obtained a finite expression for the renormalized three-loop corrections to the neutral CP-even and CP-odd Higgs boson masses at the fixed-order O(α t α 2 s ) in the rMSSM. We have computed only the dominant contributions coming from the gauge-less and the non-light-fermions limits and the approximation of vanishing external momentum transferred in the calculation of the involved self-energy functions. We have reduced the three-loop corrections, using the integration-by-parts recurrence relations, to a set of five master integrals with at most four independent mass scales whose numerical evaluation is possible thanks to the dispersion relation techniques. The three-loop vacuum integrals are expressed in terms of one and two-dimensional numerical integrals of elementary functions which can be efficiently evaluated for a general mass pattern. In this way we avoided the asymptotic expansions on the amplitudes at the integral level in the realistic mass hierarchies proposed in [22]. Thus, our work represents an independent check of the results implemented in H3m but also provides a numerical evaluation of the three-loop corrections in the whole supersymmetric parameter space without any assumption about the mass hierarchy. We have studied the numerical impact of our three-loop corrections in the value of the pole mass M h for three different benchmark limits: the m f ree h , m max h and m mod+ h scenarios. We have considered scenarios with the SUSY scale M SU SY lower than 1.2 TeV, where a fixed-order calculation have a combined theoretical uncertainty better than the estimated uncertainty from an EFT calculation [47]. We have investigated the dependence of M h on tanβ, the soft-breaking parameter A t , the renormalization scale µ r and the CP-odd mass M A . The three-loop corrections show a good perturbative behaviour, their numerical size is about ten times lesser that the two-loop predictions and the dependence on the renormalization scale is reduced by around a factor two. The contributions yield a shift in the value of the lightest Higgs boson mass of the order of 1 GeV compared with the two-loop results of FeynHiggs and a shift between −0.4 and 0.4 GeV compared with the three-loop predictions of H3m combined with the lower-order contributions of FeynHiggs. Our results are in agreement within the theoretical composited uncertainties with the predictions of M h obtained from H3m if the same renormalization scheme and kinematical limits of the SUSY-QCD mass spectrum are employed. The corrections presented in this paper were computed without taking any assumption about the soft SUSY-breaking mechanism so that a direct comparison with the H3m results in the MSUGRA scenarios requires a reparametrization of our expressions and therefore is beyond the scope of this article. Moreover, there is still missing a comparison with the other three-loop contributions currently found in literature, namely the pure DR three-loop results recently presented in [24] and implemented in FlexibleSUSY+Himalaya, the EFT calculations implemented in HSSUSY [57], the 2-loop + O(α t α 2 s ) 3-loop results as implemented in FeynHiggs via log-resum and also the FeynHiggs full corrections including all logs resummed [46]. We left the numerical analysis of these missing comparisons as the subject of future works.

APPENDICES A Dispersion Method for Three-Loop Vacuum Integrals
The three-loop vacuum integrals are an important block in the determination of the three-loop corrections to the Higgs boson mass as was discussed in Section 3. Besides integrals that factorize into products of one and two-loop contributions, the integrals U4, U5 and U6 defined as where constitute a basis of master integrals for the three-loop corrections to the neutral CP-even and CP-odd Higgs boson self-energies in the approximation of vanishing external momentum transferred. The indices ν j are integers numbers, ε = (4 − D)/2 and D is the number of the space-time dimensions. In this work we have used a method based on dispersion relations [29] in order to get a numerical integration of the finite part of the three-loop vacuum integrals. Using the dispersion relation techniques the four-propagator function U 4 and the five-propagator function U 5 can be represented in terms of one-dimensional integrals. In analogous way the six-propagator function U 6 can be represented by a two-dimensional numerical integral. For the case of the U 4 function the representation is where ∆I db,f in (s, m U 4 div contains the UV divergences of U 4 while ∆B0 and ∆B 0,mj are the discontinuities of the scalar one-loop self-energy function, B0, and its mass derivative, B 0,mj = ∂ ∂m 2 j B0, given by Here λ(x, y, z) is the Källen function defined as λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + yz + zx) and Θ is the Heaviside step function. The finite part of the vacuum integral U 4 is therefore expressed as a numerical integral of a combination of elementary functions, such a logarithms and square roots, which can be efficiently evaluated with numerical methods for general mass patterns without make any assumptions about the mass hierarchy of the SUSY particles. Analogous but more complicated expressions can be derived for the functions U 5 and U 6 and can be consulted in detail in [29,40]. The implementation of these integrals can be found in the file U45.cc of the public code TVID-Version-1.1 (http://www.pitt.edu/∼afreitas/). Internally the program uses the Gauss-Kronrod routine QAG from the Quadpack library [58] to evaluate the dispersion integrals. This routine has been amended to facilitate 30 digit floating point arithmetic from the package doubledouble [59] allowing results with at least ten digits of precision for the U 4 and U 5 integrals and a precision of eight digits for the U 6 vacuum integral.
As we are working at two-loop order in the DR scheme, the left-handed and right-handed components of the squark fields have the same renormalization constants, as a consequence δ (l) Zf 11 = δ (l) Zf 22 and therefore we have used δ (l) Zf 21 = 0 in the derivation of the above two-loop expressions. At one-loop level we have imposed the conditions Dp (1) ii m 2 fi = 0 ; i = 1, 2 , leading to the renormalization constants: At two-loop order the renormalization conditions are derived by imposing that is to say, the finite part of the inverse propagator must be zero. Explicitly the two-loop conditions are: The corresponding two-loop counter-terms in the DR scheme are: δ (2)m2 fj = Dp ; j = 1, 2 , 2 m 2 f1 −m 2 f2 δ (2) θf = Dp δ (2) Z ii = − Dp and