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

TL;DR

This paper presents a three-loop quantum correction calculation for the lightest Higgs boson mass in the MSSM, extending previous results to the entire parameter space with advanced numerical techniques.

Contribution

It provides the first three-loop order calculation of the Higgs mass in the MSSM using dimensional reduction and dispersion relation techniques for arbitrary mass scales.

Findings

Extended the validity of Higgs mass calculations to the full MSSM parameter space.

Developed new numerical methods for three-loop vacuum integrals with arbitrary masses.

Achieved higher precision in Higgs mass predictions within supersymmetric models.

Abstract

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.

Figures

Tables

Table 1. Table 1: Input parameters for the m h m ​ a ​ x superscript subscript 𝑚 ℎ 𝑚 𝑎 𝑥 m_{h}^{max} and m h m ​ o ​ d + superscript subscript 𝑚 ℎ limit-from 𝑚 𝑜 𝑑 m_{h}^{mod+} scenarios.

	$M_t$	$M_{SUSY}$	$X_t$	$M_{\tilde{g}}$	$\mu$
$m_h^{max}$	173.2 GeV	1000 GeV	2$M_{SUSY}$	1500 GeV	200 GeV
$m_h^{mod+}$	173.2 GeV	1000 GeV	+1.5$M_{SUSY}$	1500 GeV	200 GeV

Equations

$$\displaystyle H_{1}=\left(\begin{array}[]{c}H_{1}^{0}+\frac{v_{1}}{\sqrt{2}}\\ H_{1}^{-}\end{array}\right)$$

$${\mathcal{L}}_{\mathcal{H}}={\partial }_{\sigma }{H}_1^{†}{\partial }^{\sigma }{H}_1+{\partial }_{\sigma }{H}_2^{†}{\partial }^{\sigma }{H}_2-V(H_1,H_2).$$

$$V(H_1,H_2)=\left({\left|\mu \right|}^2+m_{H_1}^2\right){\left|H_1\right|}^2+\left({\left|\mu \right|}^2+m_{H_2}^2\right){\left|H_2\right|}^2+b\left({ϵ}_{ab}{H}_1^a{H}_2^b+h.c.\right)$$

$$\text{leavevmode}\text{leavevmode}+\text{leavevmode}\text{leavevmode}\frac{1}{2}{g}^2{\left|H_1^{†}{H}_2\right|}^2+\frac{1}{8}\left(g^2+g^{\prime 2}\right){\left({\left|H_2\right|}^2-{\left|H_1\right|}^2\right)}^2.$$

$$\displaystyle\left(\begin{array}[]{c}\chi_{1}^{0}\\ \chi_{2}^{0}\end{array}\right)=\frac{1}{\sqrt{2}}D^{\dagger}(\beta)\left(\begin{array}[]{c}\pi^{0}\\ A\end{array}\right),$$

$$\displaystyle\left(\begin{array}[]{cc}H&h\end{array}\right)D(\alpha)\left(\begin{array}[]{cc}M_{Z_{0}}^{2}c_{\beta}^{2}+M_{A_{0}}^{2}s_{\beta}^{2}+\frac{T_{01}}{\sqrt{2}\>v_{1}}&-\left(M_{A_{0}}^{2}+M_{Z_{0}}^{2}\right)s_{\beta}c_{\beta}\\ -\left(M_{A_{0}}^{2}+M_{Z_{0}}^{2}\right)s_{\beta}c_{\beta}&M_{Z_{0}}^{2}s_{\beta}^{2}+M_{A_{0}}^{2}c_{\beta}^{2}+\frac{T_{02}}{\sqrt{2}\>v_{2}}\end{array}\right)D^{\dagger}(\alpha)\left(\begin{array}[]{c}H\\ h\end{array}\right)$$

$$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ +\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \left(\begin{array}[]{cc}\pi^{0}&A\end{array}\right)\left(\begin{array}[]{cc}m_{\pi_{0}}^{2}&0\\ 0&M_{A_{0}}^{2}+\frac{T_{01}}{\sqrt{2}\>v_{1}}c_{\beta}^{2}+\frac{T_{02}}{\sqrt{2}\>v_{2}}s_{\beta}^{2}\end{array}\right)\left(\begin{array}[]{c}\pi^{0}\\ A\end{array}\right),$$

$$M_{A_0}^2=b\left(cot\beta +tan\beta \right),$$

$$tan(2\alpha )=tan(2\beta )\frac{M_{A_0}^2+M_{Z_0}^2}{M_{A_0}^2-M_{Z_0}^2}$$

$$X_0\to X+\delta X$$

$$T_{1,2}^{tree}=0,$$

$${\hat{\sum }}_{{\psi }_i{\psi }_j}^{(3)}={\sum }_{{\psi }_i{\psi }_j}^{(3)}\text{leavevmode}-\text{leavevmode}{\delta }^{(3)}{M}_{{\psi }_i{\psi }_j}^2$$

$${\delta }^{(3)}{M}_{hh}^2\text{leavevmode}=\text{leavevmode}{\delta }^{(3)}{M}_{AA}^2-\left(\frac{h_t}{2M_t}{s}_{\beta }\right)\left(\frac{T_h^{(3)}}{\sqrt{2}}{s}_{4\beta }+\frac{T_H^{(3)}}{\sqrt{2}}{s}_{2\beta }^2\right),$$

$${\delta }^{(3)}{M}_{HH}^2\text{leavevmode}=\text{leavevmode}\left(\frac{h_t}{2M_t}{s}_{\beta }\right)\left(\frac{T_h^{(3)}}{\sqrt{2}}{s}_{4\beta }-\frac{T_H^{(3)}}{\sqrt{2}}{c}_{2\beta }^2\right),$$

$${\delta }^{(3)}{M}_{hH}^2\text{leavevmode}=\text{leavevmode}\left(\frac{h_t}{2M_t}{s}_{\beta }\right)\left(\frac{T_h^{(3)}}{\sqrt{2}}{c}_{4\beta }+\frac{T_H^{(3)}}{2\sqrt{2}}{s}_{4\beta }\right),$$

$$\displaystyle\left(\begin{array}[]{c}T_{H}\\ T_{h}\end{array}\right)=D\left(\beta\right)\left(\begin{array}[]{c}T_{1}\\ T_{2}\end{array}\right)$$

$$\displaystyle\left(\Delta_{H}\right)^{-1}=-i\left(\begin{array}[]{cc}p^{2}-M_{H}^{2}+\left.\widehat{\sum}\right._{HH}&\left.\widehat{\sum}\right._{hH}\\ \left.\widehat{\sum}\right._{hH}&p^{2}-M_{h}^{2}+\left.\widehat{\sum}\right._{hh}\end{array}\right),$$

$$\left[p^2-M_H^2+{\hat{\sum }}_{HH}\right]\left[p^2-M_h^2+{\hat{\sum }}_{hh}\right]-{\left[{\hat{\sum }}_{hH}\right]}^2=0.$$

$$\frac{3}{2}\left(m_1^2+m_2^2+m_3^2\right)-\frac{p^2}{4}-\sum _{j=1}^3{m}_j^{2}ln\left(m_j^2/{\mu }_r^2\right).$$

$$\displaystyle U4\left(m_{1},m_{2},m_{3},m_{4}\right)\equiv{\rm INT}\left[2,1,1,1,0,0\right],\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\displaystyle U5\left(m_{1},m_{2},m_{3},m_{4},m_{5}\right)\equiv{\rm INT}\left[1,1,1,1,1,0\right],\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$U6\left(m_1,m_2,m_3,m_4,m_5,m_6\right)\equiv \mathrm{INT}\left[1,1,1,1,1,1\right],$$

$$\mathrm{INT}\left[{\nu }_1,{\nu }_2,{\nu }_3,{\nu }_4,{\nu }_5,{\nu }_6\right]=i\frac{e^{3{\gamma }_E\epsilon }}{{\pi }^{3D/2}}\int \prod _{j=1}^3{d}^D{q}_j\frac{1}{{\left[q_1^2-m_1^2\right]}^{{\nu }_1}{\left[{\left(q_1-q_2\right)}^2-m_2^2\right]}^{{\nu }_2}}\times$$

$$\frac{1}{{\left[{\left(q_2-q_3\right)}^2-m_3^2\right]}^{{\nu }_3}{\left[q_3^2-m_4^2\right]}^{{\nu }_4}{\left[q_2^2-m_5^2\right]}^{{\nu }_5}{\left[{\left(q_1-q_3\right)}^2-m_6^2\right]}^{{\nu }_6}},$$

$$U4\left(m_1,m_2,m_3,m_4\right)=-\frac{e^{{\gamma }_E\epsilon }}{i{\pi }^{D/2}}\int d^D{q}_3{\int }_0^{\infty }ds\frac{\Delta I_{db}\left(s\right)}{q_3^2-s+i\epsilon }$$

$$\text{leavevmode}\text{leavevmode}=\text{leavevmode}U{4}_{div}-{\int }_0^{\infty }dsA{0}_{fin}\left(s\right)\Delta I_{db,fin}\left(s\right),$$

$$\Delta I_{db,fin}(s,m_1^2,m_2^2,m_3^2,m_4^2)=\Delta B_{0,m_1}\left(s,m_1^2,m_2^2\right)\mathrm{Re}\left[B0\left(s,m_3^2,m_4^2\right)-B0\left(s,0,0\right)\right]$$

$$\text{leavevmode}-\text{leavevmode}\Delta B_{0,m_1}\left(s,m_1^2,0\right)\mathrm{Re}\left[B0\left(s,m_3^2,0\right)+B0\left(s,m_4^2,0\right)-2B0\left(s,0,0\right)\right]$$

$$\text{leavevmode}+\text{leavevmode}\mathrm{Re}\left[B_{0,m_1}\left(s,m_1^2,m_2^2\right)\right]\left(\Delta B0\left(s,m_3^2,m_4^2\right)-\Delta B0\left(s,0,0\right)\right)$$

$$\text{leavevmode}-\text{leavevmode}\mathrm{Re}\left[B_{0,m_1}\left(s,m_1^2,0\right)\right]\left(\Delta B0\left(s,m_3^2,0\right)+\Delta B0\left(s,m_4^2,0\right)-2\Delta B0\left(s,0,0\right)\right).$$

$$\Delta B0\left(s,m_a^2,m_b^2\right)=\frac{1}{s}\lambda \left(s,m_a^2,m_b^2\right)\Theta \left(s-{\left(m_a+m_b\right)}^2\right),$$

$$\Delta B_{0,m_1}\left(s,m_a^2,m_b^2\right)=\frac{m_a^2-m_b^2-s}{s\lambda \left(s,m_a^2,m_b^2\right)}\Theta \left(s-{\left(m_a+m_b\right)}^2\right).$$

$$\lambda (x,y,z)=\sqrt{x^2+y^2+z^2-2(xy+yz+zx)}$$

$$\frac{{\delta }^{(1)}{M}_f}{M_f}=\mathrm{Dp}\left[{\Sigma }_f^V\left(M_f^2\right)+{\Sigma }_f^S\left(M_f^2\right)\right],$$

$${\Sigma }_f\left(p^2\right)=\cancel{p}{\Sigma }_f^V\left(p^2\right)+\cancel{p}{\gamma }_5{\Sigma }_f^A\left(p^2\right)+M_f{\Sigma }_f^S\left(p^2\right)+M_f{\gamma }_5{\Sigma }_f^P\left(p^2\right).$$

$$\displaystyle\mathcal{L}_{\tilde{m}_{f}}^{\rm{bare}}=-\frac{1}{2}\left(\begin{array}[]{cc}\tilde{f}_{L}^{\dagger}&\tilde{f}_{R}^{\dagger}\end{array}\right)M_{L,R}^{2}\left(\begin{array}[]{c}\tilde{f}_{L}\\ \tilde{f}_{R}\end{array}\right),$$

Full text

1 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\leavevmode\nobreak\ \rm{GeV}$ with a total integrated luminosity of $2\leavevmode\nobreak\ ab^{-1}$ there are proposals for the $e^{+}\leavevmode\nobreak\ e^{-}$ ILC collider to precisely measure the couplings of the $125\leavevmode\nobreak\ \rm{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 ($\Lambda_{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 ($\Lambda_{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, $\Delta M_{h}\approx 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($\alpha\alpha_{s}$) can be found in [12], in particular the O($\alpha_{b}\alpha_{s}$) [13] and O($\alpha_{t}\alpha_{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($\alpha\alpha_{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 ($\overline{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-$\overline{DR}$ scheme and including corrections of $O(p^{2}\alpha_{t}\alpha_{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_{t}^{2}\alpha_{t}\alpha_{s}^{2})$. 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 $\overline{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($\alpha_{t}\alpha_{s}^{2}$) [28]. We have followed the Feynman diagrammatic (FD) procedure to obtain a renormalized correction in the $\overline{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.

2 Renormalization of the CP-Even Higgs Boson Masses

The Higgs sector of the real MSSM Lagrangian requires the definition of two doublets with opposite hyper-charges

[TABLE]

with complex components that are vevless scalar fields, $H_{j}^{0}=\phi_{j}^{0}+i\chi_{j}^{0}$, from the expansion around the vacuum expectation values $v_{1,2}$. In terms of $H_{1,2}$ doublets the expression of the Higgs sector of the bare Lagrangian is:

[TABLE]

$V(H_{1},H_{2})$ contains the soft SUSY breaking terms:

[TABLE]

The terms proportional to the mixing parameter $\left|\mu\right|$ come from the $F$-contribution to the SUSY Lagrangian, the terms with the EW gauge couplings ($g$, $g^{\prime}$) come from the $D$-contribution while $m_{H_{1}}^{2}$, $m_{H_{2}}^{2}$ and $b$ are the soft SUSY breaking parameters [30, 31]. The potential can be rotated to a basis where the quadratic terms in the real and imaginary components of $H_{j}^{0}$ are diagonalized through the transformations

[TABLE]

In this basis the Higgs sector has five physical Higgs bosons, three of them are neutral, the other two are charged. In lowest order these are the lightest ($h$) and heavy ($H$) CP-even Higgs bosons, the CP-odd Higgs boson ($A$), and the two charged Higgs bosons ($H^{\pm}$). On the other hand, we can rewrite the bare mass Lagrangian related to the neutral fields $H_{j}^{0}$ as

[TABLE]

where

[TABLE]

and $T_{0j}$ is the coefficient of the linear term in $\phi_{j}^{0}=ReH_{j}^{0}$. The angle $\alpha$ is a function of the parameters of the theory:

[TABLE]

We have used the short notation $s_{\theta}=sin(\theta)$ and $c_{\theta}=cos(\theta)$. 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\beta$, which is the ratio of the two vacuum expectation values. The masses of the five Higgs bosons particles $h$, $H$, $A$ and $H^{\pm}$ 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($\alpha_{t}\alpha_{s}^{2}$). Consequently, all the terms proportional to $M_{Z_{0}}$ are disregarded. In this limit the tree-level relation $\alpha=\beta$ 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\leavevmode\nobreak\ \rm{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\beta$. As the wave functions renormalization of the Higgs fields don’t make contributions, it is enough to rewrite (30) in terms of the renormalized parameters $M_{A}$ and $T_{j}$ according to

[TABLE]

iii) We have defined the vevs of the Higgs fields ($v_{1}$, $v_{2}$) as the minima of the full effective potential. The condition $\left\langle\Omega\left|H^{0}_{1,2}\right|\Omega\right\rangle=0$ is satisfied order by order in the perturbative expansion of the one-point Green function implying

[TABLE]

where $T_{j}^{(l)}$ 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

[TABLE]

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:

[TABLE]

where $h_{t}$ is the top Yukawa coupling,

[TABLE]

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

[TABLE]

which is equivalent to solve the determinantal equation

[TABLE]

The three-loop counter-terms derived from the mass Lagrangian (30) 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($\alpha_{t}\alpha_{s}^{2}$) we need the O($\alpha_{s}$)-contribution of the one-loop counter-terms 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($\alpha_{s}^{2}$). We have got all the one and two-loop counter-terms in the $\overline{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.

3 Technical Details to the Three-Loop Calculation of $\rm{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 $\alpha_{s}^{2}\alpha_{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 ($\sum$) and the tadpoles ($T$) for the fields $h$, $H$ and $A$ according with the equations (35) and (36). 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 $\gamma_{5}$ matrix in the Dirac traces. We have dealt $\gamma_{5}$ as an anti-commuting object, which is allowed in our computation because traces with an odd number of $\gamma_{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 $T3$ have very well known analytical expressions for any configuration of masses, while the three-loop functions $U4$, $U5$ and $U6$ 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 space-time 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($\varepsilon$) and O($\varepsilon^{2}$) in the Laurent expansion. It is no possible to evaluate these contributions with TVID because the higher-$\varepsilon$ terms of the real and imaginary parts of $B0$ and $C0$ and therefore of $U4$, $U5$ and $U6$ 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}\neq 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 $\overline{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/$\varepsilon$ divergence with coefficient $-1/2m_{1}^{2}$, independent of the values of $m_{2}$ and $m_{3}$. $Dm_{1}J$ has poles of one and second order, its 1/$\varepsilon^{2}$ divergence is 1/2, mass independent while its 1/$\varepsilon$ divergence is $1/2-ln(m_{1}^{2}/\mu_{r}^{2})$. The 1/$\varepsilon^{2}$ divergence coefficient of $J$ is $1/2(m_{1}^{2}+m_{2}^{2}+m_{3}^{2})$ while the 1/$\varepsilon$ divergence coefficient is

[TABLE]

For the cases where one or more masses $m_{j}^{2}$ vanish in the expression (48) one should take the zero order term of the Taylor expansion around $m_{j}^{2}=0$ to get the right expression. The same apply for $T3$ but choosing the external momentum $p^{2}$ equal to zero. Finally, we need the functions $F$ and $V$ with at most three independent mass scales. The function $F$ is finite, while $V\left[p^{2},m_{1}^{2},m_{2}^{2},m_{3}^{2},m_{4}^{2}\right]$ have a divergence 1/$\varepsilon^{2}$ with coefficient 1/2 and a 1/$\varepsilon$ divergence equal to $B0_{fin}\left[p^{2},m_{2}^{2},m_{4}^{2}\right]+1/2$, where $B0_{fin}$ 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\gamma_{E}\varepsilon)]$ 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.

4 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 (47). Our corrections depend upon 26 parameters: the renormalization scale $\mu_{r}$, the SM $\overline{DR}$ parameters $h_{t}$, $M_{t}$, $\alpha_{s}$ and the MSSM parameters $\mu$, $tan\beta$, $M_{A}$, $M_{\tilde{g}}$, $\theta_{t}$, $\tilde{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_{h}^{free}$ scenario, where $\mu_{r}$, $tan\beta$, $M_{A}$ and $A_{t}$ are left as free input parameters. ii) The $m_{h}^{max}$ and $m_{h}^{mod+}$ 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_{SUSY}<1.2\leavevmode\nobreak\ \rm{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_{SUSY}=1.2\leavevmode\nobreak\ \rm{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 super-symmetric 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_{h}^{free}$ scenario we fix the parameters $M_{\tilde{t}_{L}}=M_{\tilde{t}_{R}}=\tilde{m}_{q_{1,2}}=M_{SUSY}=1\leavevmode\nobreak\ \rm{TeV}$, where $q$ denotes any quark different than the top quark. We also set $M_{\tilde{g}}=1500\leavevmode\nobreak\ \rm{GeV}$, $\mu=200\leavevmode\nobreak\ \rm{GeV}$ and $M_{A}=1000\leavevmode\nobreak\ \rm{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\beta$ (Fig-3(b)) and on the renormalization scale $\mu_{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 ($\alpha_{s}\alpha_{t}$, $\alpha_{s}\alpha_{b}$, $\alpha_{t}^{2}$, $\alpha_{t}\alpha_{b}$, $\alpha_{b}^{2}$) in the rMSSM. The one-loop field-renormalization constants and the one-loop $tan\beta$ counter-term are set in the $\overline{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 (46). Besides, the dotted curve represents the three-loop predictions at order $O(M_{t}^{2}\alpha_{t}\alpha_{s}^{2})$ 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\beta$ for small values close to $tan\beta=3$, while for large values above $tan\beta=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 $\mu_{r}$. In Fig-4(a) the dependence is studied for three different values of $tan\beta$, namely $tan\beta=3,10$ and $20$. The three-loop corrections lead to a more stable dependence of $M_{h}$ with the renormalization scale $\mu_{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\beta$ in the $m_{h}^{free}$ scenario. Fig-4(b) shows the RG evolution for three values of soft breaking parameter $A_{t}$, $A_{t}=0.3\leavevmode\nobreak\ \rm{TeV},1\leavevmode\nobreak\ \rm{TeV}$ and $1.7\leavevmode\nobreak\ \rm{TeV}$. The evolution is more stable when the three-loop corrections are added and when the value of $A_{t}$ decreases, reducing the scale dependence by a factor of about 1.6 compared to the two-loop predictions.

In the $m_{h}^{max}$ and $m_{h}^{mod+}$ scenarios the renormalization scale is set to $\mu_{r}=m_{t}=173.2\leavevmode\nobreak\ \rm{GeV}$, where $m_{t}$ represents the combined Tevatron and LHC experimental value of the top quark mass [54]. Besides one has to fix $M_{\tilde{t}_{L}}=M_{\tilde{t}_{R}}=M_{SUSY}$, $A_{q}=0$ and $\tilde{m}_{q_{1,2}}=1500$ GeV. The parameter $A_{t}$ is fixed through the mixing term in the squark sector, $X_{t}=A_{t}-\mu cot\beta$, while $M_{A}$ and $tan\beta$ are left as free parameters. Within the $m_{h}^{max}$ scenario $X_{t}$ is chosen in order to maximize the value of $M_{h}$ for a given election of $\mu$ and $\tan\beta$. This occurs when $\left|X_{t}/M_{SUSY}\right|\approx 2$, where the radiative corrections in the FD calculation reach to the largest positive contribution. With this election of parameters the mass of the lightest CP-even Higgs boson is in agreement with the LHC Higgs boson signal just in a relatively small strip in the $M_{A}$-$tan\beta$ plane. A convenient way to enlarge the region of validity is to decrease the amount of mixing in the stop sector. The $m_{h}^{mod+}$ scenario is a modification of $m_{h}^{max}$ where this mixing represented by $\left|X_{t}/M_{SUSY}\right|$ is reduced. In detail we consider the input parameters shown in Table 1 for each scenario. The stop mixing angle $\theta_{t}$ and the stop masses $\tilde{m}_{t_{1}}$, $\tilde{m}_{t_{2}}$ are functions of the parameters specified in Table 1.

To draw the plots in Figure 5 the CP-odd Higgs boson mass $M_{A}$ is varied in the interval: $500\leavevmode\nobreak\ \rm{GeV}\leq M_{A}\leq 1\leavevmode\nobreak\ \rm{TeV}$ and we consider the values of $tan\beta=3$ and $tan\beta=10$. The numerical values of our three loop corrections to $M_{h}$ are reduced compared with the two-loop predictions showing a good behaviour of the perturbative expansion. In both scenarios the three-loop corrections give rise to a significant reduction of about $1\leavevmode\nobreak\ \rm{GeV}$ to the two-loop value of the Higgs boson mass. In addition the values of the corrections are consistent with the results obtained with H3m. The difference can amount a size of about 0.2 GeV in the $m_{h}^{max}$ limit and 0.05 GeV in $m_{h}^{mod+}$. The detailed analysis of the individual sources of uncertainty of the three-loop fixed order $\overline{DR}\>^{\prime}$ Higgs boson mass prediction of SoftSUSY, developed recently in [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.

5 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(\alpha_{t}\alpha_{s}^{2})$ 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_{h}^{free}$, $m_{h}^{max}$ and $m_{h}^{mod+}$ scenarios. We have considered scenarios with the SUSY scale $M_{SUSY}$ lower than $1.2\leavevmode\nobreak\ \rm{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\beta$, the soft-breaking parameter $A_{t}$, the renormalization scale $\mu_{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\leavevmode\nobreak\ \rm{GeV}$ compared with the two-loop results of FeynHiggs and a shift between $-0.4\leavevmode\nobreak\ \rm{and}\leavevmode\nobreak\ 0.4\leavevmode\nobreak\ \rm{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 $\overline{DR}$ three-loop results recently presented in [24] and implemented in FlexibleSUSY+Himalaya, the EFT calculations implemented in HSSUSY [57], the 2-loop + O($\alpha_{t}\alpha_{s}^{2}$) 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.

Acknowledgements. We would like to thank to Dr. S. Borowka for very useful discussions about the technical aspects of our three-loop computation. We thank to group of Prof. Dr. B. A. Kniehl for their hospitality and financial support during the stay of E. A. Reyes R. in DESY. Our computation was done using the facilities of the II. Institute for Theoretical Physics, Hamburg University. A. R. Fazio acknowledges the support received by the Erwin Schröedinger International Institute for Mathematical Physics, where the final stages of this work have been completed. This work is partially financially supported by the research grant N. 34729, Minimal Supersymmetric Standard Model Higgs Boson to High Accuracy of the call CONVOCATORIA NACIONAL DE PROYECTOS PARA EL FORTALECIMIENTO DE LA INVESTIGACION, CREACION E INNOVACION DE LA UNIVERSIDAD NACIONAL DE COLOMBIA 2016-2018.

APPENDICES

Appendix 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

[TABLE]

where

[TABLE]

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 $\nu_{j}$ are integers numbers, $\varepsilon=(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 $U4$ and the five-propagator function $U5$ can be represented in terms of one-dimensional integrals. In analogous way the six-propagator function $U6$ can be represented by a two-dimensional numerical integral. For the case of the $U4$ function the representation is

[TABLE]

where

[TABLE]

$U4_{div}$ contains the UV divergences of $U4$ while $\Delta B0$ and $\Delta B_{0,m_{j}}$ are the discontinuities of the scalar one-loop self-energy function, $B0$, and its mass derivative, $B_{0,m_{j}}=\frac{\partial}{\partial m_{j}^{2}}B0$, given by

[TABLE]

Here $\lambda(x,y,z)$ is the Källen function defined as

[TABLE]

and $\Theta$ is the Heaviside step function. The finite part of the vacuum integral $U4$ 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 $U5$ and $U6$ 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 $U4$ and $U5$ integrals and a precision of eight digits for the $U6$ vacuum integral.

Appendix B Renormalization of Non-Local Divergences

This appendix presents analytic expressions for the counter-terms needed in the sub-renormalization procedure discussed in Section 3. We have adopted a quite simple renormalization of the fermion masses, the gluino and the top quark masses are renormalized at one-loop in the $\overline{DR}$ scheme according with the condition:

[TABLE]

where the fermion self-energy was decomposed into a vector, an axial-vector, a scalar and a pseudo-scalar part as

[TABLE]

The function Dp$\left[\>f\>\right]$ takes the divergent part of the argument $f$. By its own side renormalization of the squark sector requires a more detail description. The necessary counter-terms emerge from the mass term of the bare squark Lagrangian, which is given by

[TABLE]

where

[TABLE]

for the squarks type up. An analogous expression for the squarks type down can be obtained just by changing $cot\beta\rightarrow tan\beta$ and the u-type soft SUSY breaking term $M_{\tilde{U}}^{2}$ by the d-type term $M_{\tilde{D}}^{2}$. The physical mass eigenstates $\tilde{m}_{f_{1}}$ and $\tilde{m}_{f_{2}}$ can be found by diagonalizing the mass matrix (66) through the orthogonal transformation

[TABLE]

In order to derive the counter-term expressions it is convenient to express the squark mass matrix in terms of the physical masses and the squark mixing angle $\theta_{\tilde{f}}$ as follows:

[TABLE]

The renormalization constants of the masses, the mixing angles and the fields are then defined via the transformations:

[TABLE]

where

[TABLE]

The bare Lagrangian (63) is therefore transformed to its renormalized version

[TABLE]

with

[TABLE]

Collecting the counter-terms generated from (93) and from the kinetic terms of the squark Lagrangian, we derive the expressions of the renormalized squark self-energies, which are given by

[TABLE]

Once again we have used the letter sigma with a hat to represent the renormalized self-energies. Keeping the terms up to two-loop order in the unrenormalized self-energy matrix

[TABLE]

as well as in $\Delta M_{1,2}^{2}$ it is possible to derived the components of (97) at one and two -loop level:

[TABLE]

[TABLE]

[TABLE]

As we are working at two-loop order in the $\overline{DR}$ scheme, the left-handed and right-handed components of the squark fields have the same renormalization constants, as a consequence $\delta^{(l)}Z_{\tilde{f}_{11}}=\delta^{(l)}Z_{\tilde{f}_{22}}$ and therefore we have used $\delta^{(l)}Z_{\tilde{f}_{21}}=0$ in the derivation of the above two-loop expressions. At one-loop level we have imposed the conditions

[TABLE]

leading to the renormalization constants:

[TABLE]

At two-loop order the renormalization conditions are derived by imposing

[TABLE]

that is to say, the finite part of the inverse propagator must be zero. Explicitly the two-loop conditions are:

[TABLE]

[TABLE]

The corresponding two-loop counter-terms in the $\overline{DR}$ scheme are:

[TABLE]

[TABLE]

[TABLE]

In the NLF limit the mixing angles $\theta_{\tilde{q}}$ with $q=u,d,b,c,s$ are equal to zero and therefore there are no mixing counter-terms. Thus, the pole equation (112) and the counter-term (114) don’t exist for the $q$-type squarks. Additionally, the equation (111) and the counter-terms (113) and (115) can be even further reduced just by putting $\delta^{(1)}\theta_{\tilde{q}}=0$.

Appendix C Renormalization Constants

In this appendix we are going to list the relevant counter-terms involved in our calculation. The one-loop counter-term of the top quark mass in the $\overline{DR}$ scheme is:

[TABLE]

The gluino mass counter-term at one-loop order is:

[TABLE]

In the NLF limit $M_{q}=0$ where $q=u,d,b,c,s$. From now on we are going to preserve the definitions of the indices $f$ and $q$.

The counter-terms of the squark sector are:

[TABLE]

The expression for $\delta^{(1)}\tilde{m}_{f_{2}}^{2}$ is obtained by interchanging $1\leftrightarrow 2$ in the above equation. Besides

[TABLE]

and

[TABLE]

The one-loop Passarino-Veltman functions $A0$ and $B0$ have the Laurent expansions:

[TABLE]

and

[TABLE]

where $\lambda(x,y,z)$ is the Källen function. For a complex number $z$ the logarithm is defined as $Log(z)=ln|z|+i\phi$, where $\phi\leavevmode\nobreak\ \in\leavevmode\nobreak\ (-\pi,\pi]$.

At two-loop level we need the $\overline{DR}$ counter-terms:

[TABLE]

[TABLE]

The expression for $\delta^{(2)}\tilde{m}_{t_{2}}$ can be derived from (124) by interchanging the indices $1\leftrightarrow 2$ and changing the sign of $\theta_{\tilde{t}}$. Furthermore in the $\overline{DR}$ scheme the two-loop counter-terms of the stop masses have a dependence on the $\varepsilon$-scalar mass $m_{\varepsilon}$. This dependence cancels in the final result of our three-loop correction to $M_{h}$ if we renormalize the $\varepsilon$-scalar mass in the $\overline{DR}$ scheme. It is enough to consider the one-loop counter-term

[TABLE]

Finally, we need the two-loop counter-term of the stop mixing angle:

[TABLE]