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

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 MSSM is studied in a pure $\overline{\text{DR}}$ context. For this purpose, we implement the results of Kant et al. 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 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.


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 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 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 (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], Flexi-bleSUSY [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 twoloop 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 CP-even 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 threeloop 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 three-loop 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 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 α 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 summarizes all the contributions that enter the final Higgs mass prediction in FlexibleSUSY+Himalaya. Section 4 analyzes 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 Section 5. Technical details of Himalaya, its link to FlexibleSUSY, 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 stabilize 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) 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 minimizes 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 Section 1, a complete account of the top-and bottom-Yukawa effects to order α 2 s would require 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 renormalized 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) symbolizes 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 for M (αtα 2 s ) (mt) as well as 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 Figure 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 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 DR-renormalized 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 Σ S,L,R t (p 2 , Q) denote the scalar (superscript S), and the left-and right-handed parts (L, R) of the DR renormalized 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 bottomquark mass m b and the down-type VEV at the scale M Z as  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 diagonalizing the loop-corrected mass matrix 1 at the momenta p 2 = M 2 h and p 2 = M 2 H , respectively (M 2L and M 3L are evaluated at p 2 = 0). The one-loop correction M 1L (p 2 ) contains the full one-loop MSSM Higgs self energy and tadpole contributions, including electroweak corrections and the momentum dependence. The two-loop correction M 2L contains the known corrections of order [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 Section 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 Figure 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 Section 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  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 twoloop 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 Figure 2, but 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, Figure 3 shows the renormalization scale dependence of the one-, two-and three-loop Higgs pole mass for the scenario defined in Section 2.2 with tan β = 5 and X t = 0. The one-and two-loop    Figure 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.
calculations correspond to the original FlexibleSUSY. In the one-loop calculation the threshold corrections to α s and y t are set to zero, and in the two-loop calculation the oneloop 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 Section 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 β-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 MSSM within the interval [M S /2, 2M S ], while in HSSUSY it is varied in the Standard Model within the interval [M t /2, 2M 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 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.

Comparison with lower order and EFT results
In Figures 4-5 (M Z ) = 0.1184 and G F = 1.1663787 · 10 −5 GeV −2 . All DR soft-breaking mass parameters as well as the µ parameter of the super-potential in the 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 − µ/ tan β is defined in the DR scheme and left as a free parameter. The lightest CP-even Higgs pole mass is calculated at the scale Q = √ mt ,1 mt ,2 .       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 two-loop 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 2 S . The O v 2 /M 2 S corrections calculated in Ref. [66] have not been taken into account here. Consider first Figure 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 2 S , 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 Figure 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 on-shell 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 Section 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 fixedorder 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 ≈ 1 TeV, instead of a more sophisticated matching procedure along the lines of Ref. [31,56]. Nevertheless, the latter is clearly desirable through order α t α 2 s , 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 | 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 Section 2.1), its numerical α t α 2 s contribution does occasionally differ slightly from the one of H3m.
In Figure 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 Figure 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 one-loop 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., O (α t + α b )α s + (α t + α b ) 2 + α 2 τ + (α t + α b )α 2 s , further reduces the Higgs mass by up to 2 GeV, as shown by the red solid line. 4 The right panel of Figure 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 Figure 4, we observe that the higher-order terms lower the predicted Higgs mass and bring 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.
h HSSUSY Figure 7: Comparison of the lightest Higgs pole mass calculated at the one-, twoand three-loop level with FlexibleSUSY, FlexibleSUSY+Himalaya, H3m and HSSUSY as a function of the SUSY scale for the "heavy sfermions" scenario of Ref. [68]. The horizontal orange band shows the measured Higgs mass M h = (125.09 ± 0.32) GeV including its experimental uncertainty.  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 FlexibleSUSY. These three-loop contributions have been available in the public program H3m before, where they were combined with the onshell 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 Flexi-bleSUSY 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 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 SUSY models.

Acknowledgments
We would like to thank Luminita Mihaila, Matthias Steinhauser, and Nikolai Zerf for helpful comments on the manuscript, and valuable help in the comparison with H3m. Further thanks go to Pietro Slavich for his valuable comments, in particular for pointing out an inconsistency in Section 3.1 of the original manuscript. Alexander Bednyakov kindly provided the general two-loop SQCD corrections to the running top and bottom masses in the MSSM in Mathematica format. RVH would like to thank the theory group at NIKHEF, where part of this work was done, for their kind hospitality. AV would like to thank the Institute for Theoretical Physics (ITP) in Heidelberg for its warm hospitality. Financial support for this work was provided by DFG.

A Installation of Himalaya
Himalaya can be downloaded as compressed package from [47]. After the package has been extracted, Himalaya can be configured and compiled by running cd $HIMALAY_PATH mkdir build cd build cmake .. make where $HIMALAY_PATH is the path to the Himalaya directory. When the compilation has finished, the build directory will contain the Himalaya library libHimalaya.a. For convenience, a library named libDSZ.a is created in addition, which contains the two-loop O(α t α s ) corrections from Ref. [12].

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 MSSM models: • MSSMNoFVHimalaya: This model represents the MSSM without (s)fermion flavour violation, where tan β is fixed at the scale M Z and the other SUSY parameters are fixed at a user-defined input scale. The parameters µ and Bµ are fixed by the electroweak symmetry breaking conditions. The SUSY mass spectrum, including the Higgs pole masses, is calculated at the scale Q = √ mt ,1 mt ,2 , where mt ,i are the two 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 ). The file LesHouches.out.MSSMNoFVHimalaya will then contain the SUSY particle spectrum in SLHA format. Alternatively, the Mathematica interface of FlexibleSUSY can be used: math -run " < < \" models / MSSMNoFVHimalaya / run_MSSMNoFVHimalaya . m \"" For each model an example SLHA input file and an example Mathematica script can be found in models/<model>/. (M Z ) as described in Section 3.1. To achieve that in FlexibleSUSY, the global threshold correction loop order (EXTPAR [7]) must be set to 1 (or higher) and the specific threshold correction loop order for α s (3rd digit from the right in EXTPAR [24]) 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 (EXTPAR [7]) 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 EXTPAR [24]) to 2 in the SLHA input file: Here, <model> is the used FlexibleSUSY model from above, i.e. either MSSMNoFVHimalaya, MSSMNoFVatMGUTHimalaya or NUHMSSMNoFVHimalaya.

C Configuration options to calculate the Higgs mass at three-loop level with FlexibleSUSY
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:  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: At the Mathematica level we recommend to use: D Himalaya interface Input parameters. To calculate the three-loop corrections to the light CP-even Higgs pole mass at order O α t α 2 s + α b α 2 s with Himalaya, the set of DR parameters is needed, which is shown in the following code snippet. The parameters are stored in the struct Parameters which contains the following members: Here, the integer mdrFlag 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 i 0 /M DSZ h , where δ 2L i is defined in Eq. (1), and i 0 denotes the "optimal" hierarchy as determined by the procedure of Section 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 HierarchyObject offers a set of member functions which provide access to all intermediate results. These functions are summarized in Table 1. The selection method described in Section 2 is also applied to the (s)bottom contri- getExpUncertainty(int loops) Returns the uncertainty of the expansion at the given loop order (cf. Section 2.1).

getDMh(int loops)
Returns the Higgs mass matrix proportional to α t or α b at the given loop order. Note that at the two-loop level only corrections of order O(α t α s ) are considered. Its arguments are a HierarchyObject, the Higgs mass matrix massMatrix up to the loop order of interest, and three flags (oneLoopFlag, twoLoopFlag, threeLoopFlag) to define the desired loop orders. Using the member function calculateDMh, the returned HierarchyObject provides the user with the quantity δ conv