Higgs mass prediction in the MSSM at three-loop level in a pure $\overline{{\text {DR}}}$ context

Abstract

The impact of the three-loop effects of order $\alpha _t\alpha _s^2$ on the mass of the light CP-even Higgs boson in the ${\text {MSSM}}$ is studied in a pure $\overline{{\text {DR}}}$ context. For this purpose, we implement the results of Kant et al. (JHEP 08:104, 2010) into the C++ module Himalaya and link it to FlexibleSUSY, a Mathematica and C++ package to create spectrum generators for BSM models. The three-loop result is compared to the fixed-order two-loop calculations of the original FlexibleSUSY and of FeynHiggs, as well as to the result based on an EFT approach. Aside from the expected reduction of the renormalization scale dependence with respect to the lower-order results, we find that the three-loop contributions significantly reduce the difference from the EFT prediction in the TeV-region of the ${\text {SUSY}}$ scale ${M_S}$. Himalaya can be linked also to other two-loop $\overline{{\text {DR}}}$ codes, thus allowing for the elevation of these codes to the three-loop level.

1 Introduction

The measurement of the Higgs boson mass at the Large Hadron Collider (${\text {LHC}}$) represents a significant constraint on the viability of supersymmetric (${\text {SUSY}}$) models. Given a particular ${\text {SUSY}}$ model, the mass of the Standard Model-like Higgs boson is a prediction, which must be in agreement with the measured value of $(125.09 \pm 0.21 \pm 0.11)\,\text {GeV}$ [2]. Noteworthy, the experimental uncertainty on the measured Higgs mass has already reached the per-mille level. Theory predictions in ${\text {SUSY}}$ models, however, struggle to reach the same level of accuracy. The reason is that the Higgs mass receives large higher-order corrections, dominated by the top Yukawa and the strong gauge coupling [3,4,5]. Both of these two couplings are comparatively large, leading to a relatively slow convergence of the perturbative series. Furthermore, the scalar nature of the Higgs implies corrections proportional to the square of the top-quark mass, on top of the top-mass dependence due to the Yukawa coupling, which enters the loop corrections quadratically. On the other hand, corrections from ${\text {SUSY}}$ particles are only logarithmic in the ${\text {SUSY}}$ particle masses due to the assumption of only soft ${\text {SUSY}}$-breaking terms. If the ${\text {SUSY}}$ particles are not too far above the TeV scale [6, 7], the ${\text {SUSY}}$ Higgs mass can be obtained from a fixed-order calculation of the relevant one- and two-point functions with external Higgs fields. In this case, higher-order corrections up to the three-loop level are known in the Minimal Supersymmetric Standard Model (${\text {MSSM}}$) [1, 5, 8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23].

There are plenty of publicly available computer codes which calculate the Higgs pole mass(es) in the ${\text {MSSM}}$ at higher orders: CPsuperH [24,25,26], FeynHiggs [9, 27,28,29,30,31], FlexibleSUSY [32, 33], H3m [1, 20], ISASUSY [34], MhEFT [35], SARAH/SPheno [36,37,38,39,40,41,42], SOFTSUSY [43, 44], SuSpect [45] and SusyHD [46]. FeynHiggs adopts the on-shell scheme for the renormalization of the particle masses, while all other codes express their results in terms of $\overline{{\text {MS}}}$/$\overline{{\text {DR}}}$ parameters. All these schemes are formally equivalent up to higher orders in perturbation theory, of course. The numerical difference between the schemes is one of the sources of theoretical uncertainty on the Higgs mass prediction, however. All of these programs take into account one-loop corrections, most of them also leading two-loop corrections. H3m is the only one which includes three-loop corrections of order $\alpha _t\alpha _s^2$, where $\alpha _t$ is the squared top Yukawa and $\alpha _s$ is the strong coupling. It combines these terms with the on-shell two-loop result of FeynHiggs after transforming the $\mathcal {O}\!\left( \alpha _t\right)$ and $\mathcal {O}\!\left( \alpha _t\alpha _s\right)$ terms from there to the $\overline{{\text {DR}}}$ scheme.

Here we present an alternative implementation of the $\mathcal {O}\!\left( \alpha _t\alpha _s^2\right)$ contributions of Refs. [1, 20] for the light CP-even Higgs mass in the ${\text {MSSM}}$ into the framework of FlexibleSUSY  [32], referring to the combination as FlexibleSUSY+Himalaya in what follows. This allows us to study the effect of the three-loop contributions in a pure $\overline{{\text {DR}}}$ environment, i.e. without the trouble of combining the corrections with an on-shell calculation. The three-loop terms are provided in the form of a separate C++ package, named Himalaya, which one should be able to include in any other $\overline{{\text {DR}}}$ code without much effort. The Himalaya package and the dedicated version of FlexibleSUSY, which incorporates the three-loop contributions from Himalaya, can be downloaded from Refs. [47, 48], respectively. In this way, we hope to contribute to the on-going effort of improving the precision of the Higgs mass prediction in the ${\text {MSSM}}$.

In the present paper we study the impact of the three-loop corrections for low and high ${\text {SUSY}}$ scales and compare our results to the two-loop calculations of the public spectrum generators of FlexibleSUSY and FeynHiggs. By quantifying the size of the three-loop corrections, we also provide a measure for the theoretical uncertainty of the $\overline{{\text {DR}}}$ fixed-order calculation.

As will be shown below, the implementation of the $\alpha _t\alpha _s^2$ corrections also applies to the terms of order $\alpha _b\alpha _s^2$, where $\alpha _b$ is the bottom Yukawa coupling. Therefore, Himalaya will take such terms into account, and we will refer to the sum of top- and bottom-Yukawa induced supersymmetric ${\text {QCD}}$ (${\text {SQCD}}$) corrections as $\mathcal {O}\!\left( \alpha _t\alpha _s^2+\alpha _b\alpha _s^2\right)$ in what follows. However, it should be kept in mind that this does not include effects of order $\alpha _s^2\sqrt{\alpha _t\alpha _b}$, which arise from three-loop Higgs self energies involving both a top/stop and a bottom/sbottom triangle. The results of Himalaya are thus unreliable in the (rather exotic) case where $\alpha _t$ and $\alpha _b$ are comparable in magnitude.

The remainder of this paper is structured as follows. Section 2 describes the form in which the three-loop contributions of order $(\alpha _t+\alpha _b)\alpha _s^2$ are implemented in Himalaya. Its input parameters are to be provided in the $\overline{{\text {DR}}}$ scheme at the appropriate perturbative order. Section 3 details how this input is prepared in the framework of FlexibleSUSY. It also summarises all the contributions that enter the final Higgs mass prediction in FlexibleSUSY+Himalaya. Section 4 analyses the impact of various three-loop contributions on this prediction as well as the residual renormalization scale dependence, and it compares the results obtained with FlexibleSUSY+Himalaya to existing fixed-order and resummed results for the light Higgs mass. In particular, this includes a comparison to the original implementation of the three-loop effects in H3m. Our conclusions are presented in Sect. 5. Technical details of Himalaya, its link to FlexibleSUSY, and run options are collected in the appendix.

2 Higgs mass prediction at the three-loop level in the ${\text {MSSM}}$

The results for the three-loop $\alpha _t\alpha _s^2$ corrections to the Higgs mass in the ${\text {MSSM}}$ have been obtained in Refs. [1, 20] by a Feynman diagrammatic calculation of the relevant one- and two-point functions with external Higgs fields in the limit of vanishing external momenta. The dependence of these terms on the squark and gluino masses was approximated through asymptotic expansions, assuming various hierarchies among the masses of the ${\text {SUSY}}$ particles. For details of the calculation we refer to Refs. [1, 20].

2.1 Selection of the hierarchy

A particular set of parameters typically matches several of the hierarchies mentioned above. In order to select the most suitable one, Ref. [1] suggested a pragmatic approach, namely the comparison of the various asymptotic expansions to the exact expression at two-loop level. Himalaya also adopts this approach, but introduces a few refinements in order to further stabilise the hierarchy selection (see also Ref. [49]).

In a first step the Higgs pole mass $M_h$ is calculated at the two-loop level at order $\alpha _t\alpha _s$ using the result of Ref. [12] in the form of the associated FORTRAN code provided by the authors. We refer to this quantity as $M_h^\text {DSZ}$ in what follows. Subsequently, for all hierarchies i which fit the given mass spectrum, $M_h$ is calculated again using the expanded expressions of Ref. [1] at the two-loop level, resulting in $M_{h,i}$. In the original approach of Ref. [1], the hierarchy is selected as the value of i for which the difference

$$\begin{aligned} \delta ^{\mathrm {2L}}_i = \left| M^\text {DSZ}_{h} - M_{h,i}\right| \end{aligned}$$

(1)

is minimal. However, we found that this criterion alone causes instabilities in the hierarchy selection in regions where several hierarchies lead to similar values of $\delta ^\mathrm {2L}_i$. We therefore refine the selection criterion by taking into account the quality of the convergence in the respective hierarchies, quantified by

$$\begin{aligned} \delta ^{\mathrm {conv}}_i = \sqrt{\sum _{j=1}^n\left( M_{h,i} - M^{(j)}_{h,i}\right) ^2}. \end{aligned}$$

(2)

While $M_{h,i}$ includes all available terms of the expansion in mass (and mass difference) ratios, in $M^{(j)}_{h}$ the highest terms of the expansion for the mass (and mass difference) ratio j are dropped. We then define the “best” hierarchy to be the one which minimises the quadratic mean of Eqs. (1) and (2),

$$\begin{aligned} \delta _i = \sqrt{\left( \delta _i^{\mathrm {2L}}\right) ^2 + \left( \delta _i^{\mathrm {conv}}\right) ^2}. \end{aligned}$$

(3)

The relevant analytical expressions for the three-loop terms of order $\alpha _t\alpha _s^2$ to the CP-even Higgs mass matrix in the various mass hierarchies are quite lengthy. However, they are accessible in Mathematica format in the framework of the publicly available program H3m. We have transformed these formulas into C++ format and implemented them into Himalaya.

The hierarchies defined in H3m equally apply to the top and the bottom sector of the ${\text {MSSM}}$, so that the results of Ref. [1] can also be used to evaluate the corrections of order $\alpha _b\alpha _s^2$ to the Higgs mass matrix. Indeed, Himalaya takes these corrections into account. However, as already pointed out in Sect. 1, a complete account of the top- and bottom-Yukawa effects to order $\alpha _s^2$ would require one to include the contribution of diagrams which involve both top/stop and bottom/sbottom loops at the same time. These were not considered in Ref. [1], and thus the Himalaya result should only be used in cases where such mixed $\sqrt{\alpha _t\alpha _b}$ terms can be neglected.

2.2 Modified $\overline{{\text {DR}}}$ scheme

By default, all the parameters of the calculation are renormalised in the $\overline{{\text {DR}}}$ scheme. However, in this scheme, one finds artificial “non-decoupling” effects [12], meaning that the two- and three-loop result for the Higgs mass depends quadratically on a ${\text {SUSY}}$ particle mass if this mass gets much larger than the others. Such terms are avoided by transforming the stop masses to a non-minimal scheme, named $\overline{{\text {MDR}}}$ (modified $\overline{{\text {DR}}}$) in Ref. [1], which mimics the virtue of the on-shell scheme of automatically decoupling the heavy particles.

If the user wishes to use this scheme rather than pure $\overline{{\text {DR}}}$, Himalaya writes the Higgs mass matrix as

$$\begin{aligned} \hat{\mathsf {M}}(\hat{m}_{\tilde{t}})&= \hat{\mathsf {M}}^\text {tree}+ \hat{\mathsf {M}}^{(\alpha _t)}(\hat{m}_{\tilde{t}}) + \hat{\mathsf {M}}^{(\alpha _t\alpha _s)}(\hat{m}_{\tilde{t}})\nonumber \\&\quad + \hat{\mathsf {M}}^{(\alpha _t\alpha _s^2)}(\hat{m}_{\tilde{t}}) + \cdots \nonumber \\&= \mathsf {M}^\text {tree}+ \mathsf {M}^{(\alpha _t)}(m_{\tilde{t}}) + \mathsf {M}^{(\alpha _t\alpha _s)}(m_{\tilde{t}})\nonumber \\&\quad + \delta \mathsf {M}(m_{\tilde{t}},\hat{m}_{\tilde{t}}) + \hat{\mathsf {M}}^{(\alpha _t\alpha _s^2)}(\hat{m}_{\tilde{t}}) + \cdots , \end{aligned}$$

(4)

where $\mathsf {M}$ and $\hat{\mathsf {M}}$ are the Higgs mass matrices in the $\overline{{\text {DR}}}$ and the $\overline{{\text {MDR}}}$ scheme, respectively, $\mathsf {M}^\text {tree}=\hat{\mathsf {M}}^\text {tree}$ is the tree-level expression, and the superscript ${}^{(x)}$ denotes the term of order $x\in \{\alpha _t,\alpha _s,\alpha _t\alpha _s,\ldots \}$. The ellipsis in Eq. (4) symbolises any terms that involve coupling constants other than $\alpha _t$ or $\alpha _s$, or higher orders of the latter. For brevity we suppress the stop mass indices “1” and “2” here. Himalaya provides the numerical results for $\hat{\textsf {M}}^{(\alpha _t\alpha _s^2)}(\hat{m}_{\tilde{t}})$ as well as

$$\begin{aligned} \delta \mathsf {M}(m_{\tilde{t}},\hat{m}_{\tilde{t}})&\equiv \left( \hat{\mathsf {M}}^{(\alpha _t)}(\hat{m}_{\tilde{t}}) + \hat{\mathsf {M}}^{(\alpha _t\alpha _s)}(\hat{m}_{\tilde{t}})\right) \nonumber \\&\quad -\left( \mathsf {M}^{(\alpha _t)}(m_{\tilde{t}}) +\mathsf {M}^{(\alpha _t\alpha _s)}(m_{\tilde{t}})\right) , \end{aligned}$$

(5)

where the $\overline{{\text {MDR}}}$ stop mass $\hat{m}_{\tilde{t}}$ is calculated from its $\overline{{\text {DR}}}$ value $m_{\tilde{t}}$ by the conversion formulas through $\mathcal {O}\!\left( \alpha _s^2\right)$, provided in Ref. [1]. Note that these conversion formulas depend on the underlying hierarchy, and may be different for $m_{\tilde{t}, 1}$ and $m_{\tilde{t}, 2}$.

Even if the result is requested in the $\overline{{\text {MDR}}}$ scheme, the output of Himalaya can thus be directly combined with pure $\overline{{\text {DR}}}$ results through $\mathcal {O}\!\left( \alpha _t\alpha _s\right)$ according to Eq. (4) in order to arrive at the mass matrix at order $\alpha _t\alpha _s^2$. Of course, one may also request the plain $\overline{{\text {DR}}}$ result from Himalaya, in which case it will simply return the numerical value for $\mathsf {M}^{(\alpha _t\alpha _s^2)}(m_{\tilde{t}})$, which can be directly added to any two-loop $\overline{{\text {DR}}}$ result.

In any case, the difference between the $\overline{{\text {DR}}}$ and $\overline{{\text {MDR}}}$ result is expected to be quite small unless the mass splitting between one of the stop masses and other, heavier, strongly interacting ${\text {SUSY}}$ particles becomes very large. As a practical example, in Fig. 1 we show the difference of the lightest Higgs mass at the three-loop level calculated in the $\overline{{\text {DR}}}$ and $\overline{{\text {MDR}}}$ scheme. All $\overline{{\text {DR}}}$ soft-breaking mass parameters, the $\mu$ parameter of the ${\text {MSSM}}$ super-potential, and the running CP-odd Higgs mass are set equal to ${M_S}$ here. The running trilinear couplings, except $A_t$, are chosen such that the sfermions do not mix. The $\overline{{\text {DR}}}$ stop mixing parameter $X_t = A_t - \mu /\tan \beta$ is left as a free parameter. For this scenario we find that the difference between the $\overline{{\text {DR}}}$ and $\overline{{\text {MDR}}}$ scheme is below $100\,\text {MeV}$ for different values of the stop mixing parameter.

Fig. 1

Difference between the lightest Higgs pole mass calculated in the $\overline{{\text {DR}}}$ scheme and the $\overline{{\text {MDR}}}$ scheme as a function of the ${\text {SUSY}}$ scale ${M_S}$ for $\tan \beta = 5$. In the left panel the soft-breaking stop and gluino mass parameters are set equal to ${M_S}$. In the right panel, we use $m_{\tilde{g}} = 2{M_S}$. We have cut off curves with non-zero $X_t$ around or below the TeV scale, where the $\overline{{\text {DR}}}$ CP-even Higgs mass becomes tachyonic at the electroweak scale

Note that, for all terms in the Higgs mass matrix except $\alpha _t$, $\alpha _t\alpha _s$, and $\alpha _t\alpha _s^2$, it is perturbatively equivalent to use either the $\overline{{\text {DR}}}$ or the $\overline{{\text {MDR}}}$ stop mass as defined above. Predominantly, this concerns the electroweak contributions as well as the terms of order $\alpha _t^2$. In this paper, we use the $\overline{{\text {DR}}}$ stop mass for these contributions.

3 Implementation into FlexibleSUSY

3.1 Determination of the ${\text {MSSM}}$ $\overline{{\text {DR}}}$ parameters

FlexibleSUSY determines the running $\overline{{\text {DR}}}$ gauge and Yukawa couplings as well as the running vacuum expectation value of the ${\text {MSSM}}$ along the lines of Ref. [50] by setting the scale to the Z-boson pole mass $M_Z$. In this approach, the following Standard Model (${\text {SM}}$) input parameters are used:

$$\begin{aligned}&\alpha _{\text {em}}^{{\text {SM}} (5)}(M_Z), \alpha _s^{{\text {SM}} (5)}(M_Z), G_F, M_Z,\nonumber \\&M_e, M_\mu , M_\tau , m_{u,d,s}(2\,\text {GeV}), m_c^{{\text {SM}} (4),\overline{{\text {MS}}}}(m_c),\nonumber \\&\quad m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}(m_b), M_t, \end{aligned}$$

(6)

where $\alpha _{\text {em}}^{{\text {SM}} (5)}(M_Z)$ and $\alpha _s^{{\text {SM}} (5)}(M_Z)$ denote the electromagnetic and strong coupling constants in the $\overline{{\text {MS}}}$ scheme in the Standard Model with five active quark flavours, and $G_F$ is the Fermi constant. $M_e$, $M_\mu$, $M_\tau$, and $M_t$ denote the pole masses of the electron, muon, tau lepton, and top quark, respectively. The input masses of the up, down and strange quark are defined in the $\overline{{\text {MS}}}$ scheme at the scale $2\,\text {GeV}$. The charm and bottom quark masses are defined in the $\overline{{\text {MS}}}$ scheme at their scale in the Standard Model with four and five active quark flavours, respectively.

The ${\text {MSSM}}$ $\overline{{\text {DR}}}$ gauge couplings $g_1$, $g_2$ and $g_3$ are given in terms of the $\overline{{\text {DR}}}$ parameters $\alpha _{\text {em}}^{{\text {MSSM}}}(M_Z)$ and $\alpha _s^{{\text {MSSM}}}(M_Z)$ in the ${\text {MSSM}}$ as

$$\begin{aligned} g_1(M_Z)&= \sqrt{\frac{5}{3}} \frac{\sqrt{4\pi \alpha _{\text {em}}^{{\text {MSSM}}}(M_Z)}}{\cos \theta _w(M_Z)}, \end{aligned}$$

(7)

$$\begin{aligned} g_2(M_Z)&= \frac{\sqrt{4\pi \alpha _{\text {em}}^{{\text {MSSM}}}(M_Z)}}{\sin \theta _w(M_Z)}, \end{aligned}$$

(8)

$$\begin{aligned} g_3(M_Z)&= \sqrt{4\pi \alpha _s^{{\text {MSSM}}}(M_Z)} . \end{aligned}$$

(9)

The couplings $\alpha _{\text {em}}^{{\text {MSSM}}}(M_Z)$ and $\alpha _s^{{\text {MSSM}}}(M_Z)$ are calculated from the corresponding input parameters as

$$\begin{aligned} \alpha _{\text {em}}^{{\text {MSSM}}}(M_Z)&= \frac{\alpha _{\text {em}}^{{\text {SM}} (5)}(M_Z)}{1 - \Delta \alpha _{\text {em}}(M_Z)},\end{aligned}$$

(10)

$$\begin{aligned} \alpha _s^{{\text {MSSM}}}(M_Z)&= \frac{\alpha _s^{{\text {SM}} (5)}(M_Z)}{1 - \Delta \alpha _s(M_Z)}, \end{aligned}$$

(11)

where the threshold corrections $\Delta \alpha _i(M_Z)$ have the form

$$\begin{aligned} \Delta \alpha _{\text {em}}(M_Z)&= \frac{\alpha _{\text {em}}}{2\pi } \Bigg (\frac{1}{3}- \frac{16}{9} \log {\frac{m_{t}}{M_Z}}\nonumber \\&\quad - \frac{4}{9} \sum _{i=1}^6 \log {\frac{m_{\tilde{u}_i}}{M_Z}} - \frac{1}{9} \sum _{i=1}^6 \log {\frac{m_{\tilde{d}_i}}{M_Z}} \nonumber \\&\quad -\frac{4}{3} \sum _{i=1}^2 \log {\frac{m_{\tilde{\chi }^+_i}}{M_Z}} - \frac{1}{3} \sum _{i=1}^6 \log \frac{m_{\tilde{e}_i}}{M_Z} \nonumber \\&\quad - \frac{1}{3} \log {\frac{m_{H^+}}{M_Z}} \Bigg ), \end{aligned}$$

(12)

$$\begin{aligned} \Delta \alpha _s(M_Z)&= \frac{\alpha _s}{2\pi } \left[ \phantom {\left. - \frac{1}{6} \sum _{i=1}^6 \left( \log {\frac{m_{\tilde{u}_i}}{M_Z}} + \log {\frac{m_{\tilde{d}_i}}{M_Z}} \right) \right] } \frac{1}{2} - 2 \log {\frac{m_{\tilde{g}}}{M_Z}} - \frac{2}{3} \log {\frac{m_t}{M_Z}} \right. \nonumber \\&\quad \left. - \frac{1}{6} \sum _{i=1}^6 \left( \log {\frac{m_{\tilde{u}_i}}{M_Z}} + \log {\frac{m_{\tilde{d}_i}}{M_Z}} \right) \right] . \end{aligned}$$

(13)

The $\overline{{\text {DR}}}$ weak mixing angle in the ${\text {MSSM}}$, $\theta _w$, is determined at the scale $M_Z$ from the Fermi constant $G_F$ and the Z pole mass via the relation

$$\begin{aligned} \sin ^2\theta _w \cos ^2\theta _w = \frac{\pi \,\alpha _{\text {em}}^{{\text {MSSM}}}}{\sqrt{2} M_Z^2 G_F (1-\delta _r)}, \end{aligned}$$

(14)

where

$$\begin{aligned} \delta _r&= \hat{\rho }\frac{{{\mathrm{Re}}}\Sigma _{W,T}(0)}{M_W^2} - \frac{{{\mathrm{Re}}}\Sigma _{Z,T}(M_Z^2)}{M_Z^2} + \delta _{{\text {VB}}} + \delta _r^{(2)}, \end{aligned}$$

(15)

$$\begin{aligned} \hat{\rho }&= \frac{1}{1-\Delta \hat{\rho }},\nonumber \\ \Delta \hat{\rho }&= {{\mathrm{Re}}}\Biggl [ \frac{\Sigma _{Z,T}(M_Z^2)}{\hat{\rho }\,M_Z^2} - \frac{\Sigma _{W,T}(M_W^2)}{M_W^2}\Biggr ] + \Delta \hat{\rho }^{(2)} . \end{aligned}$$

(16)

Here, $\Sigma _{V,T}(p^2)$ denotes the transverse part of the $\overline{{\text {DR}}}$-renormalised one-loop self-energy of the vector boson V in the ${\text {MSSM}}$. The vertex and box contributions $\delta _{{\text {VB}}}$ as well as the two-loop contributions $\delta _r^{(2)}$ are taken from Ref. [50]. The $\overline{{\text {DR}}}$ vacuum expectation values of the up- and down-type Higgs doublets are calculated by

$$\begin{aligned} v_u(M_Z)&= \frac{2 m_Z(M_Z) \sin \beta (M_Z)}{\sqrt{3/5 g_1^2(M_Z) + g_2^2(M_Z)}}, \end{aligned}$$

(17)

$$\begin{aligned} v_d(M_Z)&= \frac{2 m_Z(M_Z) \cos \beta (M_Z)}{\sqrt{3/5 g_1^2(M_Z) + g_2^2(M_Z)}}, \end{aligned}$$

(18)

where $\tan \beta (M_Z)$ is an input parameter and $m_Z(M_Z)$ is the Z boson $\overline{{\text {DR}}}$ mass in the ${\text {MSSM}}$, which is calculated from the Z pole mass at the one-loop level as

$$\begin{aligned} m_Z^2(M_Z) = M_Z^2 + {{\mathrm{Re}}}\Sigma _{Z,T}(M_Z^2) . \end{aligned}$$

(19)

In order to calculate the Higgs pole mass in the $\overline{{\text {DR}}}$ scheme at the three-loop level $\mathcal {O}\!\left( \alpha _t\alpha _s^2+\alpha _b\alpha _s^2\right)$, the $\overline{{\text {DR}}}$ top and bottom Yukawa couplings must be extracted from the input parameters $M_t$ and $m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}(m_b)$ at the two-loop level at $\mathcal {O}\!\left( \alpha _s^2\right)$. In order to achieve that, we make use of the known two-loop ${\text {SQCD}}$ contributions to the top and bottom Yukawa couplings of Refs. [51,52,53,54], as described in the following: We calculate the $\overline{{\text {DR}}}$ Yukawa couplings $y_t$ at the scale $M_Z$ from the $\overline{{\text {DR}}}$ top mass $m_t$ and the $\overline{{\text {DR}}}$ up-type VEV $v_u$ as

$$\begin{aligned} y_t(M_Z) = \sqrt{2} \frac{m_t(M_Z)}{v_u(M_Z)} . \end{aligned}$$

(20)

In our approach, we relate the $\overline{{\text {DR}}}$ top mass to the top pole mass $M_t$ at the scale $M_Z$ as

$$\begin{aligned} \begin{aligned} m_t(M_Z)&= M_t + {{\mathrm{Re}}}\Sigma _{t}^S(M_t^2,M_Z) \\&\quad + M_t \Big [ {{\mathrm{Re}}}\Sigma _{t}^L(M_t^2,M_Z) + {{\mathrm{Re}}}\Sigma _{t}^R(M_t^2,M_Z) \\&\quad + \Delta m_t^{(1),{\text {SQCD}}}(M_Z) + \Delta m_t^{(2),{\text {SQCD}}}(M_Z) \Big ], \end{aligned} \end{aligned}$$

(21)

where the $\Sigma _{t}^{S,L,R}(p^2,Q)$ denote the scalar (superscript S), and the left- and right-handed parts (L, R) of the $\overline{{\text {DR}}}$ renormalised one-loop top self-energy without the gluon, stop, and gluino contributions, and $\Delta m_t^{(1),{\text {SQCD}}}$ and $\Delta m_t^{(2),{\text {SQCD}}}$ are the full one- and two-loop ${\text {SQCD}}$ corrections taken from Refs. [51, 52],

$$\begin{aligned} \Delta m_t^{(1),\text {SQCD}}&= - \frac{\alpha _s}{4 \pi } C_F \left[ \left( \frac{m_g m_{\tilde{t}_1}^2 s_{2\theta _t}}{m_t \left( m_{\tilde{t}_1}^2-m_g^2\right) }-\frac{m_g m_{\tilde{t}_2}^2 s_{2\theta _t}}{m_t \left( m_{\tilde{t}_2}^2-m_g^2\right) }+\frac{m_{\tilde{t}_1}^4}{2 \left( m_{\tilde{t}_1}^2-m_g^2\right) ^2}\right. \right. \nonumber \\&\quad \quad \left. \left. -\frac{m_{\tilde{t}_1}^2}{m_{\tilde{t}_1}^2-m_g^2} +\frac{m_{\tilde{t}_2}^4}{2 \left( m_{\tilde{t}_2}^2-m_g^2\right) ^2}-\frac{m_{\tilde{t}_2}^2}{m_{\tilde{t}_2}^2-m_g^2} +1\right) \log \frac{m_g^2}{Q^2} \right. \nonumber \\&\quad \quad \left. + \left( -\frac{m_g m_{\tilde{t}_1}^2 s_{2\theta _t}}{m_t \left( m_{\tilde{t}_1}^2-m_g^2\right) }-\frac{m_{\tilde{t}_1}^4}{2 \left( m_{\tilde{t}_1}^2-m_g^2\right) ^2}+\frac{m_{\tilde{t}_1}^2}{m_{\tilde{t}_1}^2-m_g^2}\right) \log \frac{m_{\tilde{t}_1}^2}{Q^2}\right. \nonumber \\&\quad \quad \left. + \left( \frac{m_g m_{\tilde{t}_2}^2 s_{2\theta _t}}{m_t \left( m_{\tilde{t}_2}^2-m_g^2\right) }-\frac{m_{\tilde{t}_2}^4}{2 \left( m_{\tilde{t}_2}^2-m_g^2\right) ^2}+\frac{m_{\tilde{t}_2}^2}{m_{\tilde{t}_2}^2-m_g^2}\right) \log \frac{m_{\tilde{t}_2}^2}{Q^2} \right. \nonumber \\&\quad \quad \left. +\frac{m_{\tilde{t}_1}^2}{2 \left( m_{\tilde{t}_1}^2-m_g^2\right) }+\frac{m_{\tilde{t}_2}^2}{2 \left( m_{\tilde{t}_2}^2-m_g^2\right) }-3 \log \frac{m_t^2}{Q^2}+\frac{7}{2} \right] , \end{aligned}$$

(22)

$$\begin{aligned}&\Delta m_t^{(2),\text {SQCD}} = \left( \Delta m_t^{(1),\text {SQCD}}\right) ^2 - \Delta m_t^{(2),\text {dec}} . \end{aligned}$$

(23)

In Eq. (22), it is $C_F = 4/3$ and $s_{2\theta _t} = \sin 2\theta _t$, with $\theta _t$ the stop mixing angle. The two-loop term $\Delta m_t^{(2),\text {dec}}$ is given in Ref. [51] for general stop, sbottom, and gluino masses.

The ${\text {MSSM}}$ $\overline{{\text {DR}}}$ bottom-quark Yukawa coupling $y_b$ is calculated from the $\overline{{\text {DR}}}$ bottom-quark mass $m_b$ and the down-type VEV at the scale $M_Z$ as

$$\begin{aligned} y_b(M_Z) = \sqrt{2} \frac{m_b(M_Z)}{v_d(M_Z)} . \end{aligned}$$

(24)

We obtain $m_b(M_Z)$ from the input $\overline{{\text {MS}}}$ mass $m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}(m_b)$ in the Standard Model with five active quark flavours by first evolving $m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}(m_b)$ to the scale $M_Z$, using the one-loop ${\text {QED}}$ and three-loop ${\text {QCD}}$ renormalization group equations (RGEs). Afterwards, $m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}(M_Z)$ is converted to the $\overline{{\text {DR}}}$ mass $m_b^{{\text {SM}} (5),\overline{{\text {DR}}}}(M_Z)$ by the relation

$$\begin{aligned} m_b^{{\text {SM}} (5),\overline{{\text {DR}}}}(M_Z)&= m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}(M_Z) \nonumber \\&\quad \times \left( 1 - \frac{\alpha _s}{3 \pi } + \frac{3 g_2^2}{128 \pi ^2} + \frac{13 g_Y^2}{1152 \pi ^2}\right) . \end{aligned}$$

(25)

Finally, the ${\text {MSSM}}$ $\overline{{\text {DR}}}$ bottom mass $m_b(M_Z)$ is obtained from $m_b^{{\text {SM}} (5),\overline{{\text {DR}}}}(M_Z)$ via

$$\begin{aligned} m_b(M_Z)= & {} \frac{m_b^{{\text {SM}} (5),\overline{{\text {DR}}}}(M_Z)}{1 + \Delta m_b^{(1)} + \Delta m_b^{(2)}},\end{aligned}$$

(26)

$$\begin{aligned} \Delta m_b^{(1)}= & {} -{{\mathrm{Re}}}\Sigma _{b}^S\left( \left( m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}\right) ^2,M_Z\right) \Bigg /m_b \nonumber \\&\quad - {{\mathrm{Re}}}\Sigma _{b}^L\left( \left( m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}\right) ^2,M_Z\right) \nonumber \\&\quad - {{\mathrm{Re}}}\Sigma _{b}^R\left( \left( m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}\right) ^2,M_Z\right) , \end{aligned}$$

(27)

$$\begin{aligned} \Delta m_b^{(2)}= & {} \Delta m_b^{(2),\text {dec}} - \frac{\alpha _s}{3\pi }\Delta m_b^{(1)}, \end{aligned}$$

(28)

where $\Sigma _{b}^{S,L,R}(p^2,Q)$ are the scalar, left- and right-handed parts of the $\overline{{\text {DR}}}$ renormalised one-loop bottom quark self-energy in the ${\text {MSSM}}$, in which all Standard Model particles, except the bottom quark, the top quark and the W, Z, and Higgs bosons, are omitted. In Eq. (28) $\Delta m_b^{(2),\text {dec}}$ denotes the two-loop decoupling relation of order $\mathcal {O}\!\left( \alpha _s^2\right)$ between the $\overline{{\text {MS}}}$ bottom mass $m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}$ and the $\overline{{\text {DR}}}$ bottom mass in the ${\text {MSSM}}$ calculated in Refs. [53, 54].

Note that the matching of the ${\text {SM}}$ to the ${\text {MSSM}}$ leads to large logarithmic contributions in the ${\text {MSSM}}$ $\overline{{\text {DR}}}$ parameters in the case of a heavy ${\text {SUSY}}$ particle spectrum. These contributions can be resummed in a so-called EFT approach [31, 33, 46, 55, 56].

3.2 Calculation of the CP-even Higgs pole masses

FlexibleSUSY calculates the two CP-even Higgs pole masses $M_h$ and $M_H$ by diagonalising the loop-corrected mass matrix¹

$$\begin{aligned} \mathsf {M} = \mathsf {M}^\text {tree}+ \mathsf {M}^{1L}(p^2) + \mathsf {M}^{2L} + \mathsf {M}^{3L} \end{aligned}$$

(29)

at the momenta $p^2 = M_h^2$ and $p^2 = M_H^2$, respectively ($\mathsf {M}^{2L}$ and $\mathsf {M}^{3L}$ are evaluated at $p^2=0$). The one-loop correction $\mathsf {M}^{1L}(p^2)$ contains the full one-loop ${\text {MSSM}}$ Higgs self-energy and tadpole contributions, including electroweak corrections and the momentum dependence. The two-loop correction $\mathsf {M}^{2L}$ contains the known corrections of order $\mathcal {O}\!\left( \alpha _s(\alpha _t+ \alpha _b) + (\alpha _t+\alpha _b)^2 + \alpha _{\tau }^2\right)$ [12,13,14,15,16]. The three-loop correction $\mathsf {M}^{3L}$ incorporates the terms of order $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$ from the Himalaya package, as described in Sect. 2. In Eq. (29) all contributions are defined in the $\overline{{\text {DR}}}$ scheme by default.² The renormalization scale is chosen to be $Q = \sqrt{m_{\tilde{t}, 1}m_{\tilde{t}, 2}}$ and the $\overline{{\text {DR}}}$ parameters which enter Eq. (29) are evolved to that scale by using the three-loop RGEs of the ${\text {MSSM}}$  [57, 58]. Since the two CP-even Higgs pole masses are the output of the diagonalization of $\mathsf {M}$ but at the same time must be inserted into $\mathsf {M}^{1L}(p^2)$, an iteration over the momentum is performed for each mass eigenvalue until a fixed point for the Higgs masses is reached with sufficient precision.

4 Results

4.1 Size of three-loop contributions from different sources

In the $\overline{{\text {DR}}}$ calculation within FlexibleSUSY+Himalaya, there are three sources of contributions which affect the Higgs pole mass at order $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$: The one-loop threshold correction $\mathcal {O}\!\left( \alpha _s\right)$ to the strong coupling constant, the two-loop threshold correction $\mathcal {O}\!\left( \alpha _s^2\right)$ to the top and bottom Yukawa couplings, and the genuine three-loop contribution to the Higgs mass matrix. In Fig. 2, the impact of these three sources on the Higgs pole mass is shown relative to the two-loop calculation without these three corrections. The left panel shows the impact as a function of the ${\text {SUSY}}$ scale ${M_S}$, and the right panel as a function of the relative stop mixing parameter $X_t/{M_S}$ for the scenario defined in Sect. 2.2.

Fig. 2

Influence of different three-loop contributions to the Higgs pole mass. In the left panel we show the shift in the Higgs pole mass with respect to $M_h^{2L}(y_{t,b}^{1L},\alpha _s^{0L})$ for $\tan \beta = 5$ and $X_t = 0$ as a function of the ${\text {SUSY}}$ scale. In the right panel we fix $\tan \beta = 5$ and ${M_S}= 2\,\text {TeV}$ and vary the relative stop mixing parameter $X_t/{M_S}$

First, we observe that the inclusion of the one-loop threshold correction to $\alpha _s$, Eq. (13), (blue dashed line) leads to a significant positive shift of the Higgs pole mass of around $+ 2.5\,\text {GeV}$ for ${M_S}\approx 1\,\text {TeV}$. For larger ${\text {SUSY}}$ scales the shift increases logarithmically as is to be expected from the logarithmic terms on the r.h.s. of Eq. (13). The inclusion of the full two-loop ${\text {SQCD}}$ corrections to $y_t$ (green dash-dotted line) leads to a shift of similar magnitude, but in the opposite direction (the effect due to $y_b$ is negligible). Thus, there is a significant cancellation between the three-loop contributions from the one-loop threshold correction to $\alpha _s$ and the two-loop ${\text {SQCD}}$ corrections to $y_t$. The genuine three-loop contribution to the Higgs pole mass (black dotted line) is again positive and around $+2\,\text {GeV}$ for ${M_S}\approx 1\,\text {GeV}$. This is consistent with the findings of Ref. [1], of course. As a result, the sum of these three three-loop effects (red solid line) leads to a net positive shift of the Higgs mass relative to the two-loop result without all these corrections.

The size of the individual three-loop contributions depends on the stop mixing parameter $X_t/{M_S}$, as can be seen from the r.h.s. of Fig. 2: between minimal ($X_t/{M_S}= 0$) and maximal stop mixing ($X_t/{M_S}\approx \sqrt{6}$) the size of the individual three-loop contributions changes by 1–$2\,\text {GeV}$. For maximal (minimal) mixing, their impact is maximal (minimal). The direction of the shift is independent of $X_t/{M_S}$.

Note that the nominal two-loop result of the original FlexibleSUSY (i.e., without Himalaya) includes by default the one-loop threshold correction to $\alpha _s$ and the ${\text {SM}}$ ${\text {QCD}}$ two-loop contributions to the top Yukawa coupling [32, 33]. This means that the two-loop Higgs mass as evaluated by the original FlexibleSUSY already incorporates partial three-loop contributions. As a result, the two-loop result of the original FlexibleSUSY does not correspond to the zero-line in Fig. 2, but it is rather close to the blue dashed line. This implies that, compared to the two-loop result of the original FlexibleSUSY, the effect of the remaining $\alpha _t\alpha _s^2$ contributions in the Higgs mass prediction is negative.

4.2 Scale dependence of the three-loop Higgs pole mass

To estimate the size of the missing higher-order corrections, Fig. 3 shows the renormalization scale dependence of the one-, two- and three-loop Higgs pole mass for the scenario defined in Sect. 2.2 with $\tan \beta = 5$ and $X_t = 0$. The one- and two-loop calculations correspond to the original FlexibleSUSY. In the one-loop calculation the threshold corrections to $\alpha _s$ and $y_t$ are set to zero, and in the two-loop calculation the one-loop threshold corrections to $\alpha _s$ and the two-loop ${\text {QCD}}$ corrections to $y_t$ are taken into account. The three-loop result of FlexibleSUSY+Himalaya includes all three-loop contributions at $(\alpha _t+\alpha _b)\alpha _s^2$ discussed above, i.e. the one-loop threshold correction to $\alpha _s$, the full two-loop ${\text {SQCD}}$ corrections to $y_{t,b}$, and the genuine three-loop correction to the Higgs pole mass from Himalaya. In addition, the Higgs mass predicted at the two-loop level in the pure EFT calculation of HSSUSY is shown as the black dotted line, see Sect. 4.3. The bands show the corresponding variation of the Higgs pole mass when the renormalization scale is varied using the three-loop renormalization group equations [57,58,59,60,61,62,63] for all parameters except for the vacuum expectation values, where the $\beta$-functions are known only up to the two-loop level [64, 65]. In FlexibleSUSY and FlexibleSUSY+Himalaya, the renormalizaion scale is varied in the full ${\text {MSSM}}$ within the interval $[{M_S}/2,2{M_S}]$, while in HSSUSY it is varied in the Standard Model within the interval $[M_t/2,2 M_t]$, keeping the matching scale fixed at ${M_S}$. The plot shows that the successive inclusion of higher-order corrections reduces the scale dependence, as expected. In particular, the three-loop corrections to the Higgs mass reduce the scale dependence by around a factor two, compared to the two-loop calculation. The scale dependence of HSSUSY is almost independent of ${M_S}$, because scale variation is done within the ${\text {SM}}$ after integrating out all SUSY particles at ${M_S}$. Note that the variation of the renormalization scale only serves as an indicator of the theoretical uncertainty due to missing higher-order effects.

Fig. 3

Variation of the Higgs pole mass when the renormalization scale is varied by a factor two at which the Higgs pole mass is calculated, for $\tan \beta = 5$ and $X_t = 0$

Fig. 4

Comparison of Higgs mass predictions between two- and three-loop fixed-order programs and a two-loop EFT calculation as a function of the ${\text {SUSY}}$ scale for $\tan \beta = 5$ and $X_t = 0$. In the left panel the absolute Higgs pole mass and in the right panel the difference w.r.t. the three-loop calculation is shown (FS = FlexibleSUSY, FS+H = FlexibleSUSY+Himalaya, FH = FeynHiggs)

4.3 Comparison with lower-order and EFT results

In Figs. 4, 5, we compare the three-loop calculation of FlexibleSUSY+Himalaya (red) with other ${\text {MSSM}}$ spectrum generators. As input we use $M_t = 173.34\,\text {GeV}$, $\alpha _{\text {em}}^{{\text {SM}} (5)}(M_Z) = 1/127.95$, $\alpha _s^{{\text {SM}} (5)}(M_Z) = 0.1184$ and $G_F = 1.1663787\cdot 10^{-5} \,\text {GeV}^{-2}$. All $\overline{{\text {DR}}}$ soft-breaking mass parameters as well as the $\mu$ parameter of the super-potential in the ${\text {MSSM}}$, and the running CP-odd Higgs mass are set equal to ${M_S}$. The running trilinear couplings, except for $A_t$, are chosen such that there is no sfermion mixing. The stop mixing parameter $X_t = A_t - \mu /\tan \beta$ is defined in the $\overline{{\text {DR}}}$ scheme and left as a free parameter. The lightest CP-even Higgs pole mass is calculated at the scale $Q = \sqrt{m_{\tilde{t}, 1}m_{\tilde{t}, 2}}$.

Fig. 5

Comparison of Higgs mass predictions between two- and three-loop fixed-order programs and a two-loop EFT calculation as a function of the relative stop mixing parameter $X_t/{M_S}$ for $\tan \beta = 5$ and ${M_S}= 2\,\text {TeV}$. In the left panel the absolute Higgs pole mass and in the right panel the difference w.r.t. the three-loop calculation is shown

FlexibleSUSY 1.7.4  The blue dashed line shows the original two-loop calculation with FlexibleSUSY 1.7.4 [32]. Note that, by construction of FlexibleSUSY, this result coincides exactly with the one of SOFTSUSY 3.5.1. As described above, it includes the one-loop threshold corrections to $\alpha _s$ and the two-loop ${\text {QCD}}$ contributions to $y_t$, and it uses the three-loop RGEs of the ${\text {MSSM}}$  [57, 58]. FlexibleSUSY 1.7.4 (and SOFTSUSY) use the explicit two-loop Higgs pole mass contribution of order $\mathcal {O}\!\left( \alpha _s(\alpha _t+ \alpha _b) + (\alpha _t+\alpha _b)^2 + \alpha _{\tau }^2\right)$ of Refs. [12,13,14,15,16].

HSSUSY 1.7.4  The black dotted line has been obtained using the pure two-loop effective field theory (EFT) calculation of HSSUSY  [48]. HSSUSY is a spectrum generator from the FlexibleSUSY suite, which implements the two-loop threshold correction for the quartic Higgs coupling of the Standard Model at $\mathcal {O}\!\left( \alpha _t(\alpha _t+ \alpha _s)\right)$ when integrating out the ${\text {SUSY}}$ particles at a common ${\text {SUSY}}$ scale [46, 55]. Renormalization group running is performed down to the top mass scale using the three-loop RGEs of the Standard Model [59,60,61,62,63] and, finally, the Higgs mass is calculated at the two-loop level in the Standard Model at order $\mathcal {O}\!\left( \alpha _t(\alpha _t+ \alpha _s)\right)$. In terms of the implemented corrections, HSSUSY is equivalent to SusyHD  [46], and resums large logarithms up to NNLL level while neglecting terms of order $v^2/{M_S}^2$. The $\mathcal {O}\!\left( v^2/{M_S}^2\right)$ corrections calculated in Ref. [66] have not been taken into account here.

FeynHiggs 2.13.0-beta  The green dash-dotted line shows the Higgs mass prediction using FeynHiggs 2.13.0-beta without large log resummation [9, 27,28,29,30,31].³ FeynHiggs 2.13.0-beta includes the two-loop contributions of order $\mathcal {O}\!\left( \alpha _t\alpha _s+ \alpha _b\alpha _s+ \alpha _t^2 + \alpha _t\alpha _b\right)$.

Consider first Fig. 4. The left panel shows the Higgs mass prediction as a function of ${M_S}$ according to three codes discussed above, together with the FlexibleSUSY+Himalaya result (solid red). The stop mixing parameter $X_t$ is set to zero. The right panel shows the difference of these curves to the latter. Note that the resummed result of HSSUSY neglects terms of order $v^2/{M_S}^2$, and thus forfeits reliability towards lower values of ${M_S}$. The deviation from the fixed-order curves below ${M_S}\approx 400$ GeV clearly underlines this. In contrast, the fixed-order results start to suffer from large logarithmic contributions toward large ${M_S}$, which on the other hand are properly resummed in the HSSUSY approach. From Fig. 4, we conclude that the fixed-order $\overline{{\text {DR}}}$ result loses its applicability once ${M_S}$ is larger than a few TeV, while the deviation between the non-resummed on-shell result of FeynHiggs and HSSUSY increases more rapidly above ${M_S}\approx 1$ TeV. Note that the good agreement of FlexibleSUSY with HSSUSY above the few-TeV region is accidental, as shown in Ref. [33].

Fig. 6

Comparison of the lightest Higgs pole mass calculated at the two- and three-loop level with FlexibleSUSY, FlexibleSUSY+Himalaya and HSSUSY as a function of the ${\text {SUSY}}$ scale ${M_S}$ for $\tan \beta = 5$ and $X_t = -\sqrt{6}{M_S}$. The red band shows the size of the hierarchy selection criterion $\delta _i$. In the fixed-order calculations of FlexibleSUSY and FlexibleSUSY+Himalaya the Higgs mass becomes tachyonic for ${M_S}\lesssim 350\,\text {GeV}$

Fig. 7

Comparison of the lightest Higgs pole mass calculated at the one-, two- and three-loop level with FlexibleSUSY, FlexibleSUSY+Himalaya, H3m and HSSUSY as a function of the ${\text {SUSY}}$ scale for the “heavy sfermions” scenario of Ref. [68]. The horizontal orange band shows the measured Higgs mass $M_h = (125.09 \pm 0.32)\,\text {GeV}$ including its experimental uncertainty

The effect of the three-loop $\alpha _t\alpha _s^2$ terms on the fixed-order result is negative, as discussed in Sect. 4.1, and amounts to a few hundred MeV in the region where the fixed-order approach is appropriate. They significantly improve the agreement between the fixed-order and the resummed prediction for $M_h$ in the intermediate region of ${M_S}$, where both approaches are expected to be reliable. Between ${M_S}$ of about 500 GeV and 5 TeV, our three-loop curve from FlexibleSUSY+Himalaya deviates from the HSSUSY result by less than 300 MeV. This corroborates the compatibility of the two approaches in the intermediate region. Considering the current estimate of the theoretical uncertainty in the Higgs mass prediction [28, 33, 46, 55, 67], our observation even legitimates a naive switching between the fixed-order and the resummed approach at ${M_S}\approx 1$ TeV, instead of a more sophisticated matching procedure along the lines of Refs. [31, 56]. Nevertheless, the latter is clearly desirable through order $\alpha _t\alpha _s^2$, in particular in the light of the observations for non-zero stop mixing to be discussed below, but has to be deferred to future work at this point.

Figure 5 shows the three-loop effects as a function of $X_t$, where the value of ${M_S}=2$ TeV is chosen to be inside the intermediate region. The figure shows that, for $|X_t|\lesssim 3{M_S}$, the qualitative features of the discussion above are largely independent of the mixing parameter, whereupon the quantitative differences between the fixed-order and the resummed results are typically larger for non-zero stop mixing. Figure 6 underlines this by setting $X_t=-\sqrt{6}{M_S}$ and varying ${M_S}$. The kink in the three-loop curve originates from a change of the optimal hierarchy chosen by Himalaya. The red band shows the uncertainty $\delta _i$ as defined in Eq. (3), which is used to select the best fitting hierarchy. We find that $\delta _i$ is comparable to the size of the kink, which indicates a reliable treatment of the hierarchy selection criterion.

4.4 Comparison with other three-loop results

The three-loop $\mathcal {O}\!\left( \alpha _t\alpha _s^2\right)$ corrections to the light ${\text {MSSM}}$ Higgs mass discussed in this paper were originally implemented in the Mathematica code H3m. We checked that the implementation of the $\alpha _t$ and $\alpha _t\alpha _s$ terms in Himalaya leads to the same numerical results as in H3m, if the same set of $\overline{{\text {DR}}}$ parameters is used as input. Since the $\alpha _t\alpha _s^2$ terms of Himalaya are derived from their implementation in H3m, it is not surprising that they also result in the same numerical value if the same set of input parameters is given and the same mass hierarchy is selected. But since Himalaya has a slightly more sophisticated way of choosing this hierarchy (see Sect. 2.1), its numerical $\alpha _t\alpha _s^2$ contribution does occasionally differ slightly from the one of H3m.

In Fig. 7 we compare our results to the three-loop calculation presented in Ref. [68], assuming the input parameters for the “heavy sfermions” scenario defined in detail in the example folder of Ref. [69]. In the left panel the blue circles show the H3m result, including only the terms of $\mathcal {O}\!\left( \alpha _t+ \alpha _t\alpha _s+ \alpha _t\alpha _s^2\right)$, where the ${\text {MSSM}}$ $\overline{{\text {DR}}}$ top mass is calculated using the “running and decoupling” procedure described in Ref. [68]. The black crosses show the same result, except that the $\overline{{\text {DR}}}$ top mass at the ${\text {SUSY}}$ scale is taken from the spectrum generator FlexibleSUSY+Himalaya. We can reproduce the latter result with FlexibleSUSY+Himalaya if we take the same terms into account, i.e., $\mathcal {O}\!\left( \alpha _t+ \alpha _t\alpha _s+ \alpha _t\alpha _s^2\right)$; see the dotted red line in Fig. 7. The small differences between the two results are due to the fact that H3m works with on-shell electroweak parameters, while FlexibleSUSY+Himalaya uses $\overline{{\text {DR}}}$ parameters. The inclusion of all one-loop contributions to $M_h$ and the momentum iteration reduces the Higgs mass by 4–$6\,\text {GeV}$, as shown by the red dashed line. Including all two- and three-loop corrections which are available in FlexibleSUSY+Himalaya, i.e., $\mathcal {O}\!\left( (\alpha _t+\alpha _b)\alpha _s+ (\alpha _t+ \alpha _b)^2 + \alpha _{\tau }^2 +(\alpha _t+\alpha _b)\alpha _s^2\right)$, further reduces the Higgs mass by up to $2\,\text {GeV}$, as shown by the red solid line.⁴ The right panel of Fig. 7 shows again our one-, two-, and three-loop predictions obtained with FlexibleSUSY, FlexibleSUSY+Himalaya, as well as the EFT result of HSSUSY. Similar to Fig. 4, we observe that the higher-order terms lower the predicted Higgs mass and render it closer to the resummed result. A detailed comparison of FlexibleSUSY+Himalaya to a result where H3m is combined with the lower-order results of FeynHiggs is beyond the scope of this paper and left to a future publication.

Fig. 8

Comparison of the lightest Higgs pole mass calculated at the one-, two- and three-loop level with FlexibleSUSY, FlexibleSUSY+Himalaya and HSSUSY as a function of the lightest stop pole mass for the benchmark point of Fig. 1 of Ref. [70]. The horizontal orange band shows the measured Higgs mass $M_h = (125.09 \pm 0.32)\,\text {GeV}$ including its experimental uncertainty. The bands around the calculated Higgs mass values show the parametric uncertainty from $M_t = (173.34 \pm 0.98)\,\text {GeV}$ and $\alpha _s^{{\text {SM}} (5)}(M_Z) = 0.1184 \pm 0.0006$

Figure 8 shows the lightest ${\text {MSSM}}$ Higgs mass as obtained by FlexibleSUSY at one- and two-loop level, the FlexibleSUSY+Himalaya result, as well as the EFT prediction obtained with HSSUSY. The ${\text {MSSM}}$ parameters are defined in the $\overline{{\text {DR}}}$ scheme and are chosen in the style of Ref. [70]:⁵ The soft-breaking mass parameters of the left- and right-handed stops are set equal at the ${\text {SUSY}}$ scale ${M_S}= \sqrt{m_{\tilde{t}, 1}m_{\tilde{t}, 2}}$, i.e. $m_{\tilde{t}_L}({M_S}) = m_{\tilde{t}_R}({M_S})$. All other soft-breaking sfermion mass parameters are set to $m_{\tilde{f}}({M_S}) = m_{\tilde{t}_{L,R}}({M_S}) + 1\,\text {TeV}$. Stop mixing is disabled, $X_t({M_S}) = 0$, and the remaining trilinear couplings are set to zero at the scale ${M_S}$. The gaugino mass parameters, the super-potential $\mu$ parameter and the CP-odd $\overline{{\text {DR}}}$ Higgs mass are set to $M_1({M_S}) = M_2({M_S}) = M_3({M_S}) = 1.5\,\text {TeV}$, $\mu ({M_S}) = 200\,\text {GeV}$ and $m_A({M_S}) = {M_S}$, respectively, and we fix $\tan \beta (M_Z) = 20$. As opposed to the results shown in Fig. 1 of Ref. [70],⁶ we observe a reduction of $M_h$ towards higher loop orders, thus leading to the opposite conclusion of a heavy ${\text {SUSY}}$ spectrum in this scenario, given the current experimental value for the Higgs mass. Reassuringly, the higher-order corrections move the fixed-order result closer to the resummed result, leading to agreement between the two at the level of about $1\,\text {GeV}$ even at comparatively large ${\text {SUSY}}$ scales.

5 Conclusions

We have presented the implementation Himalaya of the three-loop $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$ terms of Refs. [1, 20] for the light CP-even Higgs mass in the ${\text {MSSM}}$, and its combination with the $\overline{{\text {DR}}}$ spectrum generator framework FlexibleSUSY. These three-loop contributions have been available in the public program H3m before, where they were combined with the on-shell calculation of FeynHiggs. With the implementation into FlexibleSUSY presented here, we were able to study the size of the three-loop contributions within a pure $\overline{{\text {DR}}}$ environment. Despite the fact that the genuine $\mathcal {O}\!\left( \alpha _t\alpha _s^2\right)$ corrections are positive [1], the combination with the two-loop decoupling terms in the top Yukawa coupling lead to an overall reduction of the Higgs mass prediction relative to the “original” two-loop FlexibleSUSY result by about 2 GeV, depending on the value of the stop masses and the stop mixing. This moves the fixed-order prediction for the Higgs mass significantly closer to the result obtained from a pure EFT calculation in the region where both approaches are expected to give sensible results. Contributions of order $\mathcal {O}\!\left( \alpha _b\alpha _s^2\right)$ are found to be negligible in all scenarios studied here.

To indicate the remaining theory uncertainty due to higher-order effects, we have varied the renormalization scale which enters the calculation by a factor two. The results show that the inclusion of the three-loop contributions reduces the scale uncertainty of the Higgs mass by around a factor two, compared to a calculation without the genuine three-loop effects. We conclude that our implementation leads to an improved CP-even Higgs mass prediction relative to the two-loop results. Our implementation of the three-loop terms should be useful also for other groups that aim at a high-precision determination of the Higgs mass in ${\text {SUSY}}$ models.

Appendices

Appendix A: Installation of Himalaya

Himalaya can be downloaded as a compressed package from [47]. After the package has been extracted, Himalaya can be configured and compiled by running

where ${\$}{\texttt {HIMALAY\_PATH}}$ is the path to the Himalaya directory. When the compilation has finished, the build directory will contain the Himalaya library

. For convenience, a library named

is created in addition, which contains the two-loop $\mathcal {O}\!\left( \alpha _t\alpha _s\right)$ corrections from Ref. [12].

Appendix B: Installation of FlexibleSUSY with Himalaya

We provide a dedicated version of FlexibleSUSY 1.7.4, which uses Himalaya to calculate the Higgs pole mass at the three-loop level. This package contains three pre-generated ${\text {MSSM}}$ models:

• MSSMNoFVHimalaya  This model represents the ${\text {MSSM}}$ without (s)fermion flavour violation, where $\tan \beta$ is fixed at the scale $M_Z$ and the other ${\text {SUSY}}$ parameters are fixed at a user-defined input scale. The parameters $\mu$ and $B\mu$ are fixed by the electroweak symmetry breaking conditions. The ${\text {SUSY}}$ mass spectrum, including the Higgs pole masses, is calculated at the scale $Q = \sqrt{m_{\tilde{t}, 1} m_{\tilde{t}, 2}}$, where $m_{\tilde{t}, i}$ are the two $\overline{{\text {DR}}}$ stop masses.

• MSSMNoFVatMGUTHimalaya  This is the same model as the MSSMNoFVHimalaya, except that the input scale is the GUT scale $M_X$, defined to be the scale where $g_1(M_X) = g_2(M_X)$.

• NUHMSSMNoFVHimalaya  This is the same model as the MSSMNoFVHimalaya, except that the soft-breaking Higgs mass parameters $m_{H_u}^2$ and $m_{H_d}^2$ are fixed by the electroweak symmetry breaking conditions.

The package FlexibleSUSY-1.7.4-Himalaya.tar.gz can be downloaded from Ref. [48]. To extract the package at the command line, run

After the extraction, FlexibleSUSY must be configured and compiled by running

See

for more options. One can use

to speed-up the compilation if

CPU cores are available. When the compilation has finished, the ${\text {MSSM}}$ spectrum generators can be run from the command line as

The file LesHouches.out.MSSMNoFVHimalaya will then contain the ${\text {SUSY}}$ particle spectrum in SLHA format. Alternatively, the Mathematica interface of FlexibleSUSY can be used:

For each model an example SLHA input file and an example Mathematica script can be found in

.

Appendix C: Configuration options to calculate the Higgs mass at three-loop level with FlexibleSUSY

To calculate the CP-even Higgs pole masses at order $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$ at the scale $Q={M_S}$, the top and bottom Yukawa couplings $y_t({M_S})$ and $y_b({M_S})$ as well as the strong coupling constant $\alpha _s({M_S})$ must be extracted from the input parameters at the appropriate loop level.

Strong coupling constant  To calculate $M_h$ at the three-loop level at $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$ correctly, $\alpha _s({M_S})$ must be extracted at the one-loop level from the input value $\alpha _s^{{\text {SM}} (5)}(M_Z)$ as described in Sect. 3.1. To achieve that in FlexibleSUSY, the global threshold correction loop order (

) must be set to 1 (or higher) and the specific threshold correction loop order for $\alpha _s$ (3rd digit from the right in

must be set to 1 (or higher) in the SLHA input file. See the next paragraph for an example.

Top and bottom Yukawa couplings  FlexibleSUSY by default determines $y_t(M_Z)$ from the top pole mass at the full one-loop level including two-loop Standard Model ${\text {QCD}}$ corrections; see Ref. [32]. The bottom Yukawa coupling $y_b(M_Z)$ is determined at the full one-loop level from the running bottom quark mass in the Standard Model with five active quark flavours, $m_b^{{\text {SM}} (5),\overline{{\text {MS}}}}(m_b)$, where $\tan \beta$-enhanced higher-order corrections are resummed. Both calculations are not sufficient for the calculation of $M_h$ at the three-loop level at $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$, because strong two-loop corrections from ${\text {SUSY}}$ particles would be missing. For this reason, the complete two-loop strong corrections to the top and bottom Yukawa couplings of Refs. [51,52,53,54] have been implemented into FlexibleSUSY. They must be activated by setting the global threshold correction loop (

) order to 2 and by setting the threshold correction loop order for $y_t$ and $y_b$ (7th and 8th digit from the right in

) to 2 in the SLHA input file:

In the Mathematica interface of FlexibleSUSY these two settings are controlled using the

and

symbols:

Here,

is the used FlexibleSUSY model from above, i.e. either MSSMNoFVHimalaya, MSSMNoFVat MGUTHimalaya or NUHMSSMNoFVHimalaya.

Three-loop corrections to the CP-even Higgs mass  To use the three-loop corrections of order $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$ to the light CP-even Higgs mass in the ${\text {MSSM}}$ from Refs. [1, 20], the pole mass and EWSB loop orders must be set to 3 in the SLHA input file. In addition, the individual three-loop corrections should be switched on, by setting the flags 26 and 27 to 1. The user can select between the $\overline{{\text {DR}}}$ and $\overline{{\text {MDR}}}$ scheme for the three-loop corrections by setting the flag 25 to 0 or 1, respectively:

In the Mathematica interface of FlexibleSUSY the pole mass and EWSB loop orders are controlled using the

and

symbols, respectively. The individual three-loop corrections can be switched on/off by using the

and

symbols. The renormalization scheme is controlled by

. The above shown SLHA input settings read in FlexibleSUSY ’s Mathematica interface

Three-loop renormalization group equations  Optionally, the known three-loop renormalization group equations can be used to evolve the ${\text {MSSM}}$ $\overline{{\text {DR}}}$ parameters from $M_Z$ to ${M_S}$ [57, 58]. To activate the three-loop RGEs, the $\beta$ function loop order must be set to 3 in the SLHA input file:

In the Mathematica interface of FlexibleSUSY the $\beta$ function loop order is controlled using the

symbol:

Recommended configuration options for FlexibleSUSY+Himalaya   We recommend to run FlexibleSUSY+Himalaya with the following SLHA configuration options:

At the Mathematica level we recommend to use:

Appendix D: Himalaya interface

Input parameters  To calculate the three-loop corrections to the light CP-even Higgs pole mass at order $\mathcal {O}\!\left( \alpha _t\alpha _s^2+ \alpha _b\alpha _s^2\right)$ with Himalaya, the set of $\overline{{\text {DR}}}$ parameters is needed, which is shown in the following code snippet. The parameters are stored in the

which contains the following members:

Table 1 Description of the member functions of the

class

All these parameters are given at the scale stored in the

variable, which is typically the ${\text {SUSY}}$ scale. The input values of the stop/sbottom masses and their associated mixing angle are optional, so their default value is set to

(

). If no input is provided, the $\overline{{\text {DR}}}$ stop masses will be calculated by diagonalising the stop mass matrix,

$$\begin{aligned} \mathcal {M}_{\tilde{t}} = \begin{pmatrix} \big (m_{\tilde{Q}}^2\big )_{33} + m_t^2 + g_{t}M_Z^2c_{2\beta } &{} \tilde{X}_{t} \\ \tilde{X}_{t} &{} \big (m_{\tilde{u}}^2\big )_{33} + m_t^2 + Q_ts_W^2M_Z^2c_{2\beta } \end{pmatrix}. \end{aligned}$$

(30)

Here, $(m_{\tilde{Q}})_{33}$ is the left third generation scalar quark mass parameter, $g_{t} = 1/2 - Q_ts_W^2$, $\tilde{X}_t = m_t(A_t - \mu \cot \beta )$, $(m_{\tilde{u}})_{33}$ the right scalar top mass parameter, $Q_t = 2/3$, $s_W$ the sine of the weak mixing angle and $c_{2\beta }=\cos (2\beta )$. The sbottom mass matrix is obtained by replacing $t\rightarrow b$ and $\tilde{u}\rightarrow \tilde{d}$ in (30) with $g_b=-(1/2+Q_bs_W^2)$, $\tilde{X}_b = m_b(A_b - \mu \tan \beta )$ and $Q_b = -1/3$.

Calculation of the three-loop corrections  All the functions which are required for the calculation of the three-loop corrections are implemented as methods of the class

.

In the context of Himalaya, the procedure described in Sect. 2 is implemented by the member function

.

Here, the integer

is optional and can be used to switch between the $\overline{{\text {DR}}}$- (0) and the $\overline{{\text {MDR}}}$-scheme (1). The $\overline{{\text {DR}}}$-scheme is chosen as default. The returned object holds all information of the hierarchy selection process, such as the best fitting hierarchy, or the relative error $\delta _{i_0}^{\mathrm {2L}}/M^{\mathrm {DSZ}}_{h}$, where $\delta _i^\mathrm {2L}$ is defined in Eq. (1), and $i_0$ denotes the “optimal” hierarchy as determined by the procedure of Sect. 2.1. The latter represents a lower limit on the expected accuracy of the expansion by comparison to the exact two-loop result $M_h^\mathrm {DSZ}$. In addition to that, the

offers a set of member functions which provide access to all intermediate results. These functions are summarised in Table 1. The selection method described in Sect. 2 is also applied to the (s)bottom contributions by replacing $t \rightarrow b$, so that only terms of order $\mathcal {O}\!\left( \alpha _b\alpha _s\right)$ are considered in the comparison. By setting the Boolean parameter

to

(

) the

function returns the

for the loop corrections proportional to $\alpha _t$ ($\alpha _b$).

Example  Function calls for the benchmark point SPS2:

Estimation of the uncertainty of the expansion  In addition to the relative error of the hierarchy choice $\delta _{i_0}^\mathrm {2L}/M_h^\text {DSZ}$ (see above), we provide a member function which returns a measure for the quality of convergence of the expansion at a given loop order, given by $\delta ^\text {conv}_{i_0}$ defined in Eq. (2), where again $i_0$ labels the “optimal” hierarchy. It can be called with

Its arguments are a

, the Higgs mass matrix

up to the loop order of interest, and three flags (

,

,

) to define the desired loop orders. Using the member function

, the returned

provides the user with the quantity $\delta _{i_0}^\text {conv}$ at two and three loops by default.

Example  For the benchmark point SPS2 one could estimate the uncertainty by calling