Higgs mass prediction in the MSSM at three-loop level in a pure DR¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{\text {DR}}}$$\end{document} context

The impact of the three-loop effects of order αtαs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _t\alpha _s^2$$\end{document} on the mass of the light CP-even Higgs boson in the MSSM\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {MSSM}}$$\end{document} is studied in a pure DR¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{\text {DR}}}$$\end{document} 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 SUSY\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {SUSY}}$$\end{document} scale MS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${M_S}$$\end{document}. Himalaya can be linked also to other two-loop DR¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{{\text {DR}}}$$\end{document} codes, thus allowing for the elevation of these codes to the three-loop level.


Introduction
The measurement of the Higgs boson mass at the Large Hadron Collider (LHC) represents a significant constraint on the viability of supersymmetric (SUSY) models. Given a particular 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 ± 0.21 ± 0.11) GeV [2]. Noteworthy, the experimental uncertainty on the measured Higgs mass has already reached the per-mille level. Theory predictions in 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 a e-mail: harlander@physik.rwth-aachen.de b e-mail: klappert@physik.rwth-aachen.de c e-mail: avoigt@physik.rwth-aachen.de 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 SUSY particles are only logarithmic in the SUSY particle masses due to the assumption of only soft SUSY-breaking terms. If the SUSY particles are not too far above the TeV scale [6,7], the SUSY Higgs mass can be obtained from a fixedorder 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 (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 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 MS/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 α t α 2 s , where α t is the squared top Yukawa and α s is the strong coupling. It combines these terms with the on-shell two-loop result of FeynHiggs after transforming the O(α t ) and O(α t α s ) terms from there to the DR scheme.
Here we present an alternative implementation of the O α t α 2 s contributions of Refs. [1,20] for the light CPeven Higgs mass in the 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 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 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 MSSM. In the present paper we study the impact of the threeloop corrections for low and high 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 DR fixed-order calculation.
As will be shown below, the implementation of the α t α 2 s corrections also applies to the terms of order α b α 2 s , where α 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 QCD (SQCD) corrections as O α t α 2 s + α b α 2 s in what follows. However, it should be kept in mind that this does not include effects of order α 2 s √ α t α b , which arise from threeloop 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 α t and α 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 (α t + α b )α 2 s are implemented in Himalaya. Its input parameters are to be provided in the 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 Flex-ibleSUSY, and run options are collected in the appendix.

Higgs mass prediction at the three-loop level in the MSSM
The results for the three-loop α t α 2 s corrections to the Higgs mass in the 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 SUSY particles. For details of the calculation we refer to Refs. [1,20].

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 α t α 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 DSZ h 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 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 δ 2L i . We therefore refine the selection criterion by taking into account the quality of the convergence in the respective hierarchies, quantified by ( While M h,i includes all available terms of the expansion in mass (and mass difference) ratios, in M 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), The relevant analytical expressions for the three-loop terms of order α t α 2 s 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 MSSM, so that the results of Ref. [1] can also be used to evaluate the corrections of order α b α 2 s 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 topand bottom-Yukawa effects to order α 2 s 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 √ α t α b terms can be neglected.

Modified DR scheme
By default, all the parameters of the calculation are renormalised in the 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 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 MDR (modified 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 DR, Himalaya writes the Higgs mass matrix aŝ where M andM are the Higgs mass matrices in the DR and the MDR scheme, respectively, M tree =M tree is the tree-level expression, and the superscript (x) denotes the term of order x ∈ {α t , α s , α t α s , . . .}. The ellipsis in Eq. (4) symbolises any terms that involve coupling constants other than α t or α s , or higher orders of the latter. For brevity we suppress the stop mass indices "1" and "2" here. Himalaya provides the numerical results forM where the MDR stop massmt is calculated from its DR value mt by the conversion formulas through O α 2 s , provided in Ref. [1]. Note that these conversion formulas depend on the underlying hierarchy, and may be different for mt ,1 and mt ,2 . Even if the result is requested in the MDR scheme, the output of Himalaya can thus be directly combined with pure DR results through O(α t α s ) according to Eq. (4) in order to arrive at the mass matrix at order α t α 2 s . Of course, one may also request the plain DR result from Himalaya, in which case it will simply return the numerical value for M (α t α 2 s ) (mt ), which can be directly added to any two-loop DR result.
In any case, the difference between the DR and MDR result is expected to be quite small unless the mass splitting between one of the stop masses and other, heavier, strongly interacting 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 DR and MDR scheme. All DR soft-breaking mass parameters, the μ parameter of the 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 DR stop mixing parameter X t = A t − μ/ tan β is left as a free parameter. For this scenario we find that the difference between the DR and MDR scheme is below 100 MeV for different values of the stop mixing parameter.
Note that, for all terms in the Higgs mass matrix except α t , α t α s , and α t α 2 s , it is perturbatively equivalent to use either the DR or the MDR stop mass as defined above. Predominantly, this concerns the electroweak contributions as well as the terms of order α 2 t . In this paper, we use the DR stop mass for these contributions.

Determination of the MSSM DR parameters
FlexibleSUSY determines the running DR gauge and Yukawa couplings as well as the running vacuum expectation value of the 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 (SM) input parameters are used: where α SM ( The couplings α MSSM em (M Z ) and α MSSM s (M Z ) are calculated from the corresponding input parameters as where the threshold corrections α i (M Z ) have the form The DR weak mixing angle in the MSSM, θ w , is determined at the scale M Z from the Fermi constant G F and the Z pole mass via the relation where Here, V,T ( p 2 ) denotes the transverse part of the DRrenormalised one-loop self-energy of the vector boson V in the MSSM. The vertex and box contributions δ VB as well as the two-loop contributions δ (2) r are taken from Ref. [50]. The DR vacuum expectation values of the up-and down-type Higgs doublets are calculated by where tan β(M Z ) is an input parameter and m Z (M Z ) is the Z boson DR mass in the MSSM, which is calculated from the Z pole mass at the one-loop level as In order to calculate the Higgs pole mass in the DR scheme at the three-loop level O α t α 2 s + α b α 2 s , the DR top and bottom Yukawa couplings must be extracted from the input parameters M t and m s . In order to achieve that, we make use of the known two-loop SQCD contributions to the top and bottom Yukawa couplings of Refs. [51][52][53][54], as described in the following: We calculate the DR Yukawa couplings y t at the scale M Z from the DR top mass m t and the DR up-type VEV v u as In our approach, we relate the DR top mass to the top pole mass M t at the scale M Z as where the S,L ,R t ( p 2 , Q) denote the scalar (superscript S), and the left-and right-handed parts (L , R) of the DR renormalised one-loop top self-energy without the gluon, stop, and gluino contributions, and m (1),SQCD t and m (2),SQCD t are the full one-and two-loop SQCD corrections taken from Refs. [51,52], In Eq. (22), it is C F = 4/3 and s 2θ t = sin 2θ t , with θ t the stop mixing angle. The two-loop term m (2),dec t is given in Ref. [51] for general stop, sbottom, and gluino masses.
The MSSM DR bottom-quark Yukawa coupling y b is calculated from the DR bottom-quark mass m b and the down-type VEV at the scale M Z as Finally, the MSSM DR bottom mass and the DR bottom mass in the MSSM calculated in Refs. [53,54].
Note that the matching of the SM to the MSSM leads to large logarithmic contributions in the MSSM DR parameters in the case of a heavy SUSY particle spectrum. These contributions can be resummed in a so-called EFT approach [31,33,46,55,56].

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 1 1 We do not distinguish between DR [12][13][14][15][16]. The three-loop correction M 3L incorporates the terms of order O α t α 2 s + α b α 2 s from the Himalaya package, as described in Sect. 2. In Eq. (29) all contributions are defined in the DR scheme by default. 2 The renormalization scale is chosen to be Q = √ mt ,1 mt ,2 and the DR parameters which enter Eq. (29) are evolved to that scale by using the three-loop RGEs of the MSSM [57,58]. Since the two CP-even Higgs pole masses are the output of the diagonalization of M but at the same time must be inserted into 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.

Size of three-loop contributions from different sources
In the DR calculation within FlexibleSUSY+Himalaya, there are three sources of contributions which affect the Higgs pole mass at order O α t α 2 s + α b α 2 s : The one-loop threshold correction O(α s ) to the strong coupling constant, the two-loop threshold correction O α 2 s 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 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.
First, we observe that the inclusion of the one-loop threshold correction to α s , Eq.(13), (blue dashed line) leads to a significant positive shift of the Higgs pole mass of around +2.5 GeV for M S ≈ 1 TeV. For larger 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 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 α s and the two-loop SQCD corrections to y t . The genuine three-loop contribution to the Higgs pole mass (black dotted line) is again positive and around +2 GeV for M S ≈ 1 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 ≈ √ 6) the size of the individual three-loop contributions changes by 1-2 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 α s and the SM 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 α t α 2 s contributions in the Higgs mass prediction is negative.

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 β = 5 and X t = 0. The one-and two-loop calculations correspond to the original Flex-ibleSUSY. In the one-loop calculation the threshold corrections to α s and y t are set to zero, and in the two-loop calculation the one-loop threshold corrections to α s and the two-loop QCD corrections to y t are taken into account. The three-loop result of FlexibleSUSY+Himalaya includes all three-loop contributions at (α t + α b )α 2 s discussed above, i.e. the one-loop threshold correction to α s , the full two-loop 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 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 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.     HSSUSY   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 β = 5 and X t = 0     [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 α s and the two-loop QCD contributions to y t , and it uses the threeloop RGEs of the MSSM [57,58]. FlexibleSUSY 1.7.4 (and SOFTSUSY) use the explicit two-loop Higgs pole mass contribution of order O α s (α t + α b ) + (α t + α b ) 2 + α 2 τ of Refs. [12][13][14][15][16].  . 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 O(α t (α t + α s )) when integrating out the SUSY particles at a common 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 twoloop level in the Standard Model at order O(α t (α t + α s )).
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 O v 2 /M S 2 corrections calculated in Ref. [66] have not been taken into account here.  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 ≈ 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 DR result loses its applicability once M S is larger than a few TeV, while the deviation between the non-resummed onshell result of FeynHiggs and HSSUSY increases more rapidly above M S ≈ 1 TeV. Note that the good agreement of FlexibleSUSY with HSSUSY above the few-TeV region is accidental, as shown in Ref. [33].
The effect of the three-loop α t α 2 s 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 fixedorder 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 fixedorder and the resummed approach at M S ≈ 1 TeV, instead of a more sophisticated matching procedure along the lines of Refs. [31,56]. Nevertheless, the latter is clearly desirable through order α t α 2 s , in particular in the light of the observa-  HSSUSY   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 SUSY scale M S for tan β = 5 and X t = − √ 6M S . The red band shows the size of the hierarchy selection criterion δ i . In the fixed-order calculations of FlexibleSUSY and FlexibleSUSY+Himalaya the Higgs mass becomes tachyonic for M S 350 GeV tions 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 | 3M 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 = − √ 6M 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 δ i as defined in Eq. (3), which is used to select the best fitting hierarchy. We find that δ i is comparable to the size of the kink, which indicates a reliable treatment of the hierarchy selection criterion.

Comparison with other three-loop results
The three-loop O α t α 2 s corrections to the light MSSM Higgs mass discussed in this paper were originally implemented in the Mathematica code H3m. We checked that the implementation of the α t and α t α s terms in Himalaya leads to the same numerical results as in H3m, if the same set of DR parameters is used as input. Since the α t α 2 s 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 α t α 2 s 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 O α t + α t α s + α t α 2 s , where the MSSM 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 DR top mass    at the 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., O α t + α t α s + α t α 2 s ; 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 DR parameters. The inclusion of all oneloop contributions to M h and the momentum iteration reduces the Higgs mass by 4-6 GeV, as shown by the red dashed line. Including all two-and three-loop corrections which are available in FlexibleSUSY+Himalaya, i.e., s , further reduces the Higgs mass by up to 2 GeV, as shown by the red solid line. 4 The right panel of Fig. 7 shows again our one-, two-, and three-loop predictions obtained with Flex-ibleSUSY, 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. 1 of Ref. [70], 6 we observe a reduction of M h towards higher loop orders, thus leading to the opposite conclusion of a heavy 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 GeV even at comparatively large SUSY scales.

Conclusions
We have presented the implementation Himalaya of the three-loop O α t α 2 s + α b α 2 s terms of Refs. [1,20] for the light CP-even Higgs mass in the MSSM, and its combination with the DR spectrum generator framework Flexible-SUSY. 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 DR environment. Despite the fact that the genuine O α t α 2 s 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 O α b α 2 s are found to be negligible in all scenarios studied here.
To indicate the remaining theory uncertainty due to higherorder 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 threeloop 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 highprecision determination of the Higgs mass in SUSY models.
See for more options. One can use to speed-up the compilation if CPU cores are available. When the compilation has finished, the MSSM spectrum generators can be run from the command line as The file LesHouches.out.MSSMNoFVHimalaya will then contain the SUSY particle spectrum in SLHA format. Alternatively, the Mathematica interface of Flexible-SUSY 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 ) must be set to 1 (or higher) and the specific threshold correction loop order for α 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.  [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:

Top and bottom Yukawa couplings
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 O α t α 2 s + α b α 2 s to the light CP-even Higgs mass in the 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 DR and 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 Flexible-SUSY's Mathematica interface Three-loop renormalization group equations Optionally, the known three-loop renormalization group equations can be used to evolve the MSSM DR parameters from M Z to M S [57,58]. To activate the three-loop RGEs, the β function loop order must be set to 3 in the SLHA input file: In the Mathematica interface of FlexibleSUSY the β 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:

Calculation of the three-loop corrections
All the functions which are required for the calculation of the threeloop 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 DR-(0) and the MDR-scheme (1). The 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 δ 2L where δ 2L i 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 DSZ h . 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 → b, so that only terms of order O(α b α s ) are considered in the comparison. By setting the Boolean parameter to ( ) the function returns the for the loop corrections proportional to α t (α 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 δ 2L i 0 /M DSZ h (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 δ 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 δ conv i 0 at two and three loops by default. Example For the benchmark point SPS2 one could estimate the uncertainty by calling