The light MSSM Higgs boson mass for large $\tan\beta$ and complex input parameters

We discuss various improvements of the prediction for the light MSSM Higgs boson mass in the hybrid framework of the public code FeynHiggs, which combines fixed-order and effective field theory results. First, we discuss the resummation of logarithmic contributions proportional to the bottom-Yukawa coupling including two-loop $\Delta_b$ resummation. For large $\tan\beta$, these improvements can lead to large upward shifts of the Higgs mass compared to the existing fixed-order calculations. Second, we improve the implemented EFT calculation by fully taking into account the effect of $\mathcal{CP}$-violating phases. As a third improvement, we discuss the inclusion of partial N$^3$LL resummation. The presented improvements will be implemented into FeynHiggs.


Introduction
The discovery of a Higgs boson at the LHC [1,2] was an important step towards the understanding of the fundamental laws of Nature. The properties of the detected particle allow a sensitive test of the predictions of the Standard Model (SM) and of theories of physics beyond the SM (BSM). In particular, in the Minimal Supersymmetric extension of the SM (MSSM) [3,4], based upon the concept of supersymmetry (SUSY), the mass of the discovered boson is not a free parameter, as in the SM, but is predicted in terms of the model parameters.
While the SM-like Higgs mass in the MSSM is smaller or equal to the mass of the Z boson at the tree-level, large quantum corrections shift it upwards towards the experimentally measured value of M h ∼ 125 GeV. In order to allow the use of the SM-like Higgs mass as a precision constraint on the MSSM parameter space, the precise determination of these quantum corrections is crucial [5].
The quantum corrections can be calculated in different frameworks. In the most direct approach, quantum corrections to the Higgs self-energies are calculated diagrammatically in the full theory (for recent works see [6][7][8][9][10][11][12][13]). This approach has the advantage of capturing all corrections at a specific order in perturbation theory. If the scale of SUSY particles is, however, much larger than the electroweak scale, large logarithms emerge in the fixed-order corrections exacerbating the behaviour of the perturbative expansion. In such situations, effective field theory (EFT) techniques allow the resummation of large logarithmic corrections (for recent works see [14][15][16][17][18][19][20][21]). Without including higher-dimensional operators into the lowenergy EFT, terms suppressed by the SUSY scale are, however, missed in this approach. 1 Therefore, the accuracy of the EFT approach can be diminished if one or more SUSY scales are comparable to the electroweak scale. In order to obtain a precise prediction for the SMlike Higgs boson mass for low, intermediary as well as high SUSY scales, both approaches -the fixed-order and the EFT approach -can be combined. Such hybrid approaches have been developed in [5,17,[22][23][24][25][26][27][28][29][30].
In this paper, we focus on the hybrid approach implemented in the publicly available code FeynHiggs [22,23,26,[31][32][33][34][35][36]. We will discuss various improvements of the incorporated EFT calculation as well as their combination with the implemented fixed-order calculation: resummation of large logarithms proportional to the bottom Yukawa coupling (including twoloop ∆ b resummation [37][38][39]), extension of the EFT calculation fully taking into account the effects of complex input parameters as well as an inclusion of partial N 3 LL resummation.
This paper is structured as follows: In Section 2, we explain how the resummation of logarithms proportional to the bottom Yukawa coupling is incorporated into the hybrid framework. The extension of the EFT calculation for complex input parameters is discussed in Section 3. We explain the implementation of partial N 3 LL resummation in Section 4. Numerical results are presented in Section 5. In Section 6, we give our conclusions. In the Appendices, we provide more details regarding the calculation of the two-loop threshold corrections (see App. A), analytic expressions for the threshold corrections to the SM-Higgs self coupling and to the couplings of the Split-SUSY model (see App. B), more details regarding the dependence of ∆ b on CP-violating phases (see App. C), and explicit formulas for the one-and two-loop logarithmic corrections to M h proportional to the bottom-Yukawa coupling (see App. D).

Resummation of logarithmic bottom Yukawa contributions
In this Section, we describe the procedure for resumming logarithmic contributions controlled by the bottom Yukawa coupling in our hybrid framework. We first describe the employed fixed-order and EFT calculations separately. Then we discuss their combination. The resummation of the bottom Yukawa coupling for large tan β is discussed in Section 2.4.

Fixed-order calculation
The fixed-order part of the calculation consists of the full one-loop and O(α t α s , α 2 t ) corrections (α s = g 2 3 /(4π) with g 2 3 being the strong gauge coupling, and α t = y 2 t /(4π) with y t being the top Yukawa coupling) implemented into FeynHiggs [40][41][42]. In these corrections CP-violating phases are fully taken into account.
Also, two-loop corrections proportional to the bottom-Yukawa coupling have been calculated using different renormalisation schemes for the sbottom sector. In Section 5, we will compare the results of two different schemes.

Renormalisation scheme 1
The present public version of FeynHiggs (version 2.16.1) includes the O(α b α s , α b α t , α 2 b ) corrections derived in [43,44] (α b = y 2 b /(4π) with y b being the bottom Yukawa coupling). For these corrections the following renormalisation scheme is employed: The squark masses, m 2 q , and mixing angles θq, and the top quark mass are renormalised in the on-shell (OS) scheme, where Σq iqi is used to denote the respective scalar self-energy, Σ L t , Σ R t and Σ S t are the coefficients in the Lorentz decomposition of the unrenormalised top-quark self-energy, Additionally, Re denotes the real part, and M t is the OS top-quark mass. The corrections in [43,44] have been calculated assuming vanishing CP-violating phases. Moreover, they have been derived in the large tan β limit which implies that the bottomquark mass m b is put to zero if it is not multiplied with tan β. In this approach, the soft SUSY-breaking masses, the trilinear couplings of the stop and sbottom sector, as well as the top-quark mass are independent parameters. In contrast, the bottom-quark mass is treated as a dependent quantity and the expression for its counterterm is derived from the equation connecting the bottom-quark mass and the sbottom mixing angle, 2 where we introduced the abbreviations s γ = sin γ, c γ = cos γ, t γ = tan γ for a generic angle γ. Therefore, the one-loop counterterm for the bottom-quark mass in this scheme has the following form, where δµ is the counterterm for the Higgsino mass parameter µ. The actual bottom quark mass, which is used in the calculation, is then given by and δs 2θ b in the bracket of Eq. (6). Therefore, m b is scale independent at the one-loop level if the Higgsino mass parameter µ and t β are renormalised in the DR scheme as assumed throughout this work. In FeynHiggs, the associated renormalisation scale is by default set equal to M t .
Due to the SU (2) L gauge symmetry, the bilinear soft SUSY-breaking parameters mb L and mt L are equal to each other at the tree level. This relation is broken at the one-loop level. The counterterms for mt L and mb L read δm 2 t L = cos 2 θt δm 2 t 1 + sin 2 θt δm 2 t 2 + (m 2 t 2 − m 2 t 1 ) sin 2θt δθt − 2 m t δm t , (7a) and are in general not equal to each other. In the following, we will assume that m 2 t L is given as an input parameter. Then the renormalised soft sbottom mass m 2 b L is given by The trilinear soft SUSY-breaking parameter A b is fixed via the sbottom-sbottom-A-boson vertex function (see [43,44] for more details). 2 We use a different sign convention for µ in comparison to [43,44].

Renormalisation scheme 2
For our present study, we, however, do not make use of the O(α b α s , α b α t , α 2 b ) corrections already implemented in FeynHiggs. Instead, we employ the O(α b α s , α b α t , α 2 b ) corrections presented in [8,10] (see also [45,46]). These also include terms subleading in tan β, allow for easier control of the renormalisation scheme and take CP-violating phases fully into account. They will be part of an upcoming FeynHiggs release. For the present work, we evaluate them, however, externally and feed the numerical result back to FeynHiggs. 3 In this scheme the soft SUSY-breaking mass mb L is not treated as an independent parameter and set equal to mt L . This implies where δm 2 t L is given by Eq. (7a). The consequence of this relation is that only one of the sbottom masses can be set on-shell. As a matter of convention, the mass of the second sbottom is defined in the on-shell scheme, This allows one to treat the mass of the bottom quark as an independent parameter which is renormalised in the DR scheme, where Σ L b , Σ R b and Σ S b are defined in analogy to Eq. (2). The trilinear soft SUSY-breaking parameter A b is also defined in the DR scheme, where Ub is the mixing matrix of the sbottom sector, and Eqs. (9)-(13) fix the renormalisation conditions for all parameters of the sector. In both schemes described above it is assumed that the stop sector is renormalised using the OS scheme. We furthermore implemented a pure DR renormalisation of the stop and sbottom sector as an additional option.

EFT calculation
We build upon the existing EFT calculation in FeynHiggs [22,23,26,36]. At the sfermion mass scale, M SUSY , all sfermions as well as the non-SM-like Higgs bosons are integrated out. 4 Performing the renormalisation-group running to lower scales and passing two additional independent thresholds for electroweakinos (charginos and neutralinos) and the gluino, the SM is recovered as EFT. 5 The currently implemented EFT calculation resums leading and next-to-leading (LL and NLL) logarithms as well as next-to-next-to-leading logarithms (NNLL) in the limit of vanishing electroweak gauge couplings. So far, however, all corrections proportional to the bottom Yukawa coupling are set to zero in the EFT calculation.
For the incorporation of the bottom Yukawa contributions, our aim was to reach the same level of accuracy as for the other corrections. For implementing LL and NLL resummation, we include the bottom Yukawa contributions to the one-loop matching condition of the SM Higgs self-coupling, λ [14]. Also the one-and two-loop RGEs are extended by the RGE of the bottom-Yukawa coupling and by bottom-Yukawa contributions to the RGEs of the other couplings (see e.g. [47]).
To achieve resummation at the NNLL level, we derive the O(α b α s , α b α t , α 2 b ) threshold corrections for λ making use of the two-loop Higgs self-energy corrections obtained in [8,10,42]. Details are given in App. A. In the limit of vanishing CP-violating phases, we find agreement with the expressions derived in [16] using the effective-potential approach. We also add the bottom Yukawa contributions to the three-loop SM RGEs in the limit of vanishing electroweak gauge couplings [48][49][50][51][52] and to the calculation of the SM MS vev at the electroweak scale.
For the EFT calculation, all sbottom input parameters are defined in the DR scheme at the scale M SUSY . We choose to define tan β in the DR scheme at the scale M SUSY .

Combination in the hybrid approach
For the combination of the fixed-order and the EFT calculation, we follow the procedure described in [22,23,26,28]. For the self-energy of the SM-like Higgs boson, the result of the fixed-order calculation, Σ FO hh (p 2 ), and the EFT result, −2λ(M t )(v MS ) 2 (with v MS being the SM MS vev at the scale M t ), are summed. Subtraction terms are used to ensure that contributions included in both results are not counted twice, In order to ease this combination, we choose to define the sbottom input parameters in the same scheme in the fixed-order and the EFT calculations: We fix them in the DR scheme at the scale M SUSY . Also tan β and µ are fixed in the DR scheme at the scale M SUSY . A complication arises through the use of the DR bottom quark mass in the EFT as well as the fixed-order calculation. After adding both results as shown in Eq. (14), the Higgs pole masses are determined taking into account the momentum dependence of the fixed-order selfenergy. As discussed in detail in [28], this momentum dependence arises only from SM-type corrections as well as contributions suppressed by the SUSY scale. In order to match the result of a pure EFT calculation, in which the Higgs pole mass is determined in the SM, we have to ensure that the SM-type corrections are evaluated using only SM quantities. The DR bottom-quark mass, however, is an MSSM quantity. For this reason, we reparametrise the SM bottom-quark contributions to the Higgs self-energies in terms of the SM MS bottom quark mass at the scale M t .
In our implementation, this is achieved by subtracting the UV-finite O(α b ) self-energy, containing only the SM contributions and parametrised in terms of the MSSM bottom quark m DR,MSSM b (M SUSY ), and then adding back the same quantity but parametrised via m MS,SM b (M t ). The explicit expression for the one-loop SM Higgs boson self-energy renormalised in the MS scheme with the tadpoles renormalised to zero is given by where the superscript "fin" means that in this expression we take only the finite part of the one-loop scalar integral function B 0 (p 2 , m 2 1 , m 2 2 ), for which we use the definition given in [53]. Since the SM-like Higgs mass is determined via an iterative solution of the pole equation, the same procedure has to be applied to the derivative of Σ MS,SM hh (p 2 ) with respect to the external momentum. The O(α b α s , α t α b , α 2 b ) SM self-energies computed in the gaugeless limit are extracted from the code FlexibleSUSY [24,27,54] and Refs. [55,56].
In order to allow for an OS definition of the stop input parameters, a conversion of the OS parameters, used in the fixed-order calculation, to the DR scheme, used in the EFT calculation, is necessary. As argued in [5,22,23,26], for this conversion only one-loop logarithmic terms should be taken into account. Only in the conversion formula for the stop mixing parameter, X t , large logarithms appear. For the present study, we extend the formula given in [23] by including the bottom Yukawa contributions, where X b is the sbottom mixing parameter (

Determination of the MSSM bottom quark mass and Yukawa coupling
The DR bottom quark mass is then determined by It is well-known that the relation between the DR bottom quark mass and the SM MS bottom Yukawa coupling includes terms proportional to tan β. For large tan β, the leading tan βenhanced terms can be resummed as described in [43,44,[57][58][59][60][61][62]. Typically, this resummation is written in the form where ∆ b stands for the tan β-enhanced terms and b for terms not enhanced by tan β. 6 We employ a similar relation for the matching of the bottom Yukawa coupling, A similar procedure for the calculation of the MSSM bottom Yukawa coupling was adopted in [16]. There, however, non-enhanced terms, b , and the threshold correction of the vev, ∆v, were included into the definition of ∆ b . In our approach, we separate them to resum only tan β enhanced corrections to the bottom-Yukawa coupling in the same way as in [43,63]. This results only in a small numerical difference since the main contribution to h DR,MSSM b (M SUSY ) comes from ∆ b (see also the discussion in Section 5). In our implementation, we include full one-loop corrections to ∆ b . The quantity ∆v is calculated at the one-loop level in the gaugeless limit. In addition, we include the leading two-loop corrections to ∆ b . These two-loop corrections are based on the results from [37][38][39]. 7 We, however, perform an expansion of ∆ b (at the one-and two-loop level) for large M SUSY omitting terms of higher-order in O(v 2 /M 2 SUSY ). In addition, we adapt the renormalisation scheme to match our scheme. More precisely, in [37,39] the soft supersymmetry breaking parameters in the stop and sbottom sectors as well as the gluino mass are renormalised onshell. Moreover, all supersymmetric particles and the top quark are decoupled from the scale dependence of the strong coupling α s . This decoupling of the top quark and the on-shell renormalisation of the top sector induces large logarithms, log(M 2 SUSY /M 2 t ), implying that the formulas in [37][38][39] are not directly applicable in our framework. Since in our case the low-energy model is the SM with possibly light gluinos (and electroweakinos), we do not decouple the top quark and the gluino. Also, to be consistent with the other parts of our 6 We include also terms proportional to m 2 b tan 2 β into b and ∆v. We checked that they are numerically irrelevant. 7 Similar results have been obtained in [64,65]. Moreover, the authors of [66] derived subleading two-loop corrections, which are not taken into account in the present work.
EFT calculation we renormalise the gluino mass and the stop/sbottom masses in the DR scheme at the matching scale Q.
In the limit of all involved non-SM masses having the same value, we obtain Here, A t is the stop trilinear coupling (A t = X t + µ * / tan β), Q is the renormalisation scale, C A = 3, C F = 4 3 and T R = 1 2 . Formulas also valid for non-degenerate masses are distributed as ancillary files together with this paper.

EFT calculation for complex input parameters
In the fixed-order approach, the dependence on CP-violating phases is known at the oneand two-loop level [35,[40][41][42]67]. In the EFT framework, the phase dependence has so far only been considered in case of a low-energy Two-Higgs-Doublet-Model [20,68]. Here, we work out the dependence on CP-violating phases in the case of the SM (and the SM plus electroweakinos and/or gluinos) as EFT, for which so far only an interpolation of the result in case of real input parameters has been available [36].
We first discuss the case of the SM as low-energy EFT. Since the SM includes no phases (apart from the CKM phase, whose effect is negligible for the determination of the Higgs mass), CP-violating effects in the full MSSM enter only via threshold corrections to real parameters. At the one-loop level, the only contribution to the matching of the Higgs self-coupling with a non-vanishing phase dependence is the electroweakino contribution. It depends on the phases of the bino and wino soft-breaking masses, φ M 1 and φ M 2 , as well as of the Higgsino mass parameter, φ µ (explicit expressions are listed in App. B). This implies that at the one-loop level, there is no dependence on the phases of the squark sector (at least if the absolute values of the squark mixing parameters, |X q |, are kept constant).
The phases of the stop and sbottom sector along with the gluino phase, φ M 3 , however, enter the matching of the Higgs self-coupling at the two-loop level. Based upon the fixedorder results presented in [8,10,41,42], we extract the dependence of the two-loop threshold correction on these phases at O(α b,t α s , α 2 b , α b α t , α 2 t ) without assumptions on the internal masses (details are given in App. A). In case of real input parameters, we find full agreement with the results of [14,16,69]. By analysing the obtained expressions, it becomes clear how the expressions derived in [14,16,69] can be generalised to the case of complex input parameters: • O(α q α s ) where q = t, b: The expression for zero phases is a polynomial in X q . To get the expression for non-zero phases every odd power of X q has to be multiplied by cos(φ Xq − φ M 3 ), and X q has to replaced by | X q |.
• O(α 2 q ) where q = t, b: The expression for zero phases is a sum of monomials in the variables X q and Y q = X q + 2 µ * sin 2β of one of three types: the monomials which contain only even powers of X q , the ones which contain only even powers of Y q and the ones which contain both X q and Y q . The latter contain only even or only odd powers of X q and Y q at the same time. To get the expression for non-zero phases, every monomial which contains odd powers of X q and Y q has to be multiplied by cos(φ Xq − φ Yq ), and every X q and Y q has to be replaced by | X q | and | Y q |, respectively.
The generalisation of the O(α b α t ) expression from the CP-conserving case to the CPviolating case is slightly more complicated since different multiplicative factors arise. Full explicit expressions in the limit of all sfermions having the same mass are given in App. B.2. Fully general expressions can be found in ancillary files distributed alongside this paper.
If the low-energy theory is the SM plus electroweakinos, effective Higgs-Higgsino-gaugino couplings are induced. These are potentially complex-valued. An explicit matching calculation at the one-loop level, however, shows that their phase is zero even if one or more of the electroweakino phases in the MSSM are non-zero. Correspondingly, also the RGEs of the SM plus electroweakinos are not modified in the presence of non-zero phases. The phases, however, enter in the threshold corrections for the bottom and top Yukawa couplings as well as the Higgs self-coupling when integrating out the electroweakinos (full expressions are listed in App. B).
In addition to the phase dependencies discussed above, also the ∆ b corrections (see Section 2.4) depend on φ µ , φ M 1,2,3 and φ At . The phase dependence of the one-loop correction has been derived in [63,70,71]. The phase dependence of the two-loop correction, which we derived based upon the result of [37][38][39] (see Section 2.4), has, however, been unknown so far. We find that this dependence is the same as for the one-loop result. Namely, Eq. (23) has to be multiplied by cos(φ µ + φ M 3 ) and Eq. (24) has to be multiplied by cos(φ µ + φ At ).
This can be understood by looking at the explicit two-loop diagrams (see App. C). They fall into three categories: either a gluon, a gluino or a sbottom quark is added to the oneloop graph. If a gluon is added, the phase dependence of the one-loop graph is obviously not changed, since the two additionally appearing strong gauge couplings do not include a phase dependence. The same is true if a sbottom quark is coupled to the one-loop graph by a four-sfermion vertex. Working in the chiral basis, it is again obvious that this coupling does not induce an additional phase dependence. The case of adding a gluino is slightly more complicated. The two additional gluon-gluino-sbottom couplings do depend on the phase of the gluino mass parameter. Working again in the chiral basis, it is easy to see that one of these two additional couplings is a left-handed coupling and the other one is a righthanded coupling. The dependence on the gluino phase cancels between the left-handed and the right-handed coupling. More details and all relevant two-loop diagrams can be found in App. C.

N 3 LL resummation
Up to now, the EFT calculation implemented in FeynHiggs was restricted to full LL and NLL resummation as well as NNLL resummation in the limit of vanishing electroweak gauge couplings. In this Section, we discuss the implementation of N 3 LL resummation at O(α t α 2 s ) based upon the work presented in [18].
The following ingredients are needed in addition to the already implemented corrections for NNLL resummation: s ) Higgs self-energy corrections, • leading QCD corrections to the three-loop RGEs of the Higgs self-coupling, the strong gauge coupling as well as the top Yukawa coupling, • O(α 3 s ) extraction of the MS top Yukawa coupling at the electroweak scale, • O(α t α 2 s ) matching condition for the Higgs self-coupling between the SM and the MSSM.
The SM O(α t α 2 s ) corrections to the Higgs self-energy have been obtained in [56,72,73]; the necessary RGEs in [73,74]. Formulas for extracting the SM MS couplings at the three-loop level can be found in [47]. The O(α t α 2 s ) matching condition of the Higgs self-coupling was computed in [18] based on the O(α 2 t α 2 s ) fixed-order calculation performed in [9,75,76]. The result is implemented in the publicly available code Himalaya [9,18]. As discussed in [18], this calculation is based on an expansion of three-loop diagrams for certain mass hierarchies. Himalaya provides an uncertainty estimate for this truncation error.
We implemented all these corrections into the EFT calculation of FeynHiggs (the link to Himalaya has already been implemented for the work presented in [5]). By default, Himalaya uses the DR scheme [77] for the renormalisation of the squark input parameters [18]. Correspondingly, also the input parameters of FeynHiggs are defined in the DR scheme if N 3 LL resummation is activated. In case of complex input parameters, we interpolate the Himalaya result. 8 The inclusion of N 3 LL resummation in the EFT calculation can also be used within the hybrid approach. In this case we, however, require that also in the fixed-order calculation the parameters entering the three-loop threshold correction are renormalised in the DR scheme. The two-loop conversion, that would be necessary between OS parameters used in the fixed-order calculation and DR parameters used in the EFT calculation, is beyond the scope of the present paper. 9

Numerical results
In this Section, we discuss the numerical effects of the various improvements discussed above. 8 An interpolation in case of a complex M 3 is not possible, since the expressions implemented in Himalaya are not dependent on the sign of M 3 . 9 The necessary two-loop squark self-energy corrections have already been calculated in [78,79].

Resummation of logarithmic bottom Yukawa contributions
Here, we investigate the numerical effect of resumming logarithmic contributions proportional to the bottom Yukawa coupling. First, we concentrate on a scenario presented in Ref. [16]. Namely, we assume that all soft SUSY-breaking masses are equal to M SUSY = 1.5 TeV except the gluino mass which is fixed by M 3 = 2.5 TeV. The stop mixing parameter is set by X t = √ 6M SUSY , and the trilinear couplings of the third generation fermions are equal to each other, A b = A τ = A t . The Higgsino mass parameter, µ, is chosen to be equal to −1.5 TeV. Due to this choice of the signs of M 3 , X t and µ the MSSM bottom Yukawa coupling is enhanced by the one-loop threshold corrections proportional to the top Yukawa coupling and the strong coupling. As in Ref. [16] all the input parameters listed above and tan β are assumed to be DR parameters at the scale M SUSY . 10 In the left panel of Fig. 1 we present results for M h in dependence on tan β. In addition to showing results obtained with the calculation presented in this paper, we display results obtained using the most recent public version of FeynHiggs (version 2.16.1). Moreover, we show the result obtained in [16] for comparison. This result was obtained in a pure EFT framework using the code HSSUSY [24,27]. In the right panel of Fig. 1 we show the bottom mass, m b , which is used in the corresponding calculations. In the case of the red dashed curve it is the "OS" bottom mass, m b , defined by Eq. (6) and in case of the blue, red, orange and green solid lines it is m DR,MSSM b (M SUSY ) given by Eq. (20). In the first step of our numerical analysis, we focus on the various EFT results in the left panel of Fig. 1: the black dot-dashed line corresponds to the result obtained in [16] (red solid line in Fig. 2 of [16]). For this curve the full LL and NLL resummation of large logarithms is performed. In addition to that, NNLL logarithms are resummed to all orders in the gaugeless limit (i.e., the electroweak gauge couplings are neglected in the two-loop threshold corrections to λ). One-loop ∆ b resummation, including O(α s , α t ) corrections, is used in the one-loop threshold correction for the bottom Yukawa coupling. The red, blue, green and orange solid lines correspond to the results of our EFT calculation with different approximation levels used in the calculation of the bottom Yukawa threshold correction. We should note here that the results presented in [16] have been obtained using the SM MS top Yukawa coupling extracted at the N 3 LO level while we by default use the NNLO value. For a proper comparison with the results of [16], we adapted our calculation to use the same level of corrections (see also the discussion in Sections 4 and 5.3). We use this determination of the SM MS top Yukawa coupling for all curves of Fig. 1.
We observe a very good agreement between our EFT result using only O(α s , α t ) corrections in the calculation of ∆ b (solid red curve), which is the same level of accuracy as used in [16], and the result of [16] (black dot-dashed curve). The absolute difference between the two curves equals ∼ 0.04 GeV for t β = 15 and ∼ 0.7 GeV for t β = 42 where the curves have a very steep behavior. This difference comes mainly from the determination of the MSSM bottom Yukawa coupling at the scale M SUSY . In [16], the threshold correction for the vacuum expectation value, ∆v, and non-enhanced terms were included in the definition of ∆ b while we do not include them (see Eq. (21)). If we include them into ∆ b as in [16], the absolute difference between our calculation and the calculation presented in [16] shrinks down even further (∼ 0.2 GeV for t β = 42).
For the green solid curve in the left plot of Fig. 1, we take into account electroweak corrections in the calculation of ∆ b in addition to the O(α s , α t ) corrections used for the solid red curve. As a consequence of Eq. (62), this choice leads to a partial cancellation in the calculation of ∆ b and hence to a suppression of the MSSM bottom mass at the scale M SUSY as one can see on the right panel of Fig. 1 showing m DR,MSSM b (M SUSY ) in dependence of tan β. This in turn reduces the downward shift in the Higgs mass by the one-loop threshold corrections to the SM Higgs self-coupling, λ, that is proportional to the bottom Yukawa coupling.
The blue solid curve in the left plot of Fig. 1 shows the prediction for M h neglecting the electroweak one-loop contributions to ∆ b but including the leading two-loop QCD corrections to ∆ b . For our parameter choice, these corrections increase the absolute value of ∆ b by approximately 5%. Correspondingly, also the MSSM bottom mass is increased as can be seen in the right plot of Fig. 1. This results in a significant change of the resulting Higgs mass for tan β 40 where the dependence on tan β is very pronounced. The orange curves in the left plot of Fig. 1 correspond to the inclusion of all corrections to ∆ b mentioned above. For the considered parameter choice, the electroweak corrections to ∆ b are roughly three times larger by absolute value than the two-loop corrections to ∆ b . This explains why the orange and the green curves lie quite close to each other.
The orange dashed curve represents the result of the hybrid calculation of M h . Namely, we have merged the fixed-order result with the NNLL EFT calculation (see Section 2). The orange solid and dashed curves differ essentially by the inclusion of terms which are suppressed by the ratio v 2 /M 2 SUSY into the hybrid result. Since in our case M SUSY is chosen above the TeV scale, the size of these terms is, as expected, quite small. Therefore, the observed good agreement between the two methods serves as a consistency check of our hybrid calculation.
Finally, the red dashed curve shows the prediction for M h obtained by FeynHiggs-2.16.1 which we ran using the default flags as explained in [36]. 11 As only modification of this version, we have used the N 3 LO instead of the NNLO SM MS top Yukawa coupling to allow for a direct comparison to the result of [16]. We see that the agreement between all the seven curves is quite good for small values of tan β, but for tan β 30 the red dashed curve shows a steep fall-off, while for the other curves the large downward shift from b/b-sector corrections sets in only at higher values of tan β. The reason for this behaviour becomes clear when looking at the right panel of Fig. 1: the red dashed curve, which corresponds to m b defined in Eq. (6), increases much more rapidly for rising tan β than the other four lines. 12 This expression for the bottom mass is inserted in the leading one-loop fixed order result which gives rise to a large downward shift of M h [80][81][82], This term grows rapidly in absolute value with increasing tan β. A similar effect occurs for all other curves but there the dependence of the bottom mass on tan β is much milder. This is a consequence of our choice of the renormalisation scheme. Namely, the bottom mass used in our setup is the DR MSSM bottom mass calculated at the scale M SUSY . All the quantities entering the calculation of ∆ b and b are also DR MSSM quantities at this scale. The most important ones are the top Yukawa coupling α t and the strong Yukawa coupling α s (see Eq. (62) in App. B). Since their values decrease with increasing scale, 13 the ∆ b correction calculated in our approach is smaller than the corresponding correction in FeynHiggs-2.16.1. In this way our approach yields more stable results for large values of tan β and for regions of the MSSM parameter space where the signs of the products µM 3 and µA t are negative.
Next, we discuss the numerical effect induced by the resummation of logarithms proportional to the bottom Yukawa coupling. First of all, to have an idea how numerically important the effect is, it is instructive to have a look at the analytic one-and two-loop expressions which one can find in App. D. The bottom mass we use in our calculation, even though being potentially enhanced by ∆ b effects, is the smallest mass taken into account in our EFT calculation. The only way corrections containing the bottom mass may become sizeable is when these terms are additionally proportional to tan β. This is the case when m b is, for example, multiplied by X b , Y t , t β or 1/c β . We only find such enhancements in the two-loop next-to-leading logarithmic contributions when the stop mixing parameter, X t , is renormalised in the OS scheme (see Eqs. (83) and (84)). As a consequence, we expect the effect of the resummation to be small if we renormalise X t in the DR scheme.
This qualitative consideration turns out to be reflected in the numerical results as one can see in the left panel of Fig. 2. The red curve corresponds to the hybrid result including the effects of the bottom Yukawa coupling only at the one-loop level in the fixed order calculation.
Hybrid calculation including: Hybrid calculation including:  [8,10]. Finally, the blue curve also contains the resummation of LL, NLL and NNLL logarithms controlled by the bottom Yukawa coupling. The same color scheme also applies to the right panel of Fig. 2 and to both plots in Fig. 3. For these plots we have picked a MSSM scenario where all soft-breaking masses and µ are equal by absolute value to the common mass scale M SUSY . Moreover, we set A b = 2.5M SUSY and t β = 45. The bino, wino and gluino masses are chosen to be positive, M 1,2,3 > 0, while the Higgsino mass parameter is negative, µ < 0.
For the left plot of Fig. 2 we have chosen First, we note that the green and the blue curves agree very well with each other for low values of M SUSY . The difference between the two amounts to only ∼ 0.3 GeV for M SUSY = 700 GeV. In this region, where the scale separation is relatively small, the resummation of higher-order logarithmic contributions is expected to be subdominant. 14 However, the two curves lie quite close to each other for the whole range of scales, even for M SUSY as high as 10 5 GeV. This is in line with our qualitative analysis above: the logarithms containing bottom Yukawa coupling are numerically small if X DR t (M SUSY ) is used as an input parameter. On the other hand, a large shift of about 10 GeV between the red and the green curves 14 It is worth noting here that the different treatment of the vacuum expectation value in the fixed-order calculation, where v G F is used, and in the EFT calculation, where v DR,MSSM (M SUSY ) is used, induces differences at order O(m 2 b α(α t + α b )) which are beyond the accuracy level of our calculation and hence may be regarded as part of the uncertainties from unknown higher-order corrections. Moreover, the different treatment of the vacuum expectation value in the two-loop threshold corrections and in the respective pieces of the fixed-order calculation, inducing a difference at the three-loop order, contributes to the small shift between the result containing the two-loop fixed order result in the b/b-sector and the one including the resummation of higher-order logarithmic contributions. A similar effect was also discussed in [26].  15 Thus, in this scenario two-loop fixed-order corrections from the b/b-sector that go beyond the ∆ b contribution are numerically important. 16 With increasing M SUSY the bottom mass m DR,MSSM b (M SUSY ) decreases and the three curves get close to each other. In the plot on the right panel of Fig. 2 we fix M SUSY = 1.5 TeV and vary X DR t . One can see that for negative X DR t all three curves give roughly the same result. This is due to the fact that the contributions to ∆ b proportional to the strong coupling and to the top Yukawa coupling partially cancel each other. Correspondingly, the bottom mass does not acquire a significant enhancement for negative X DR t . The inclusion of the two-loop fixedorder corrections as well as the resummation of the logarithms has only a small effect in this case. On the contrary, for positive X t the top Yukawa and strong corrections to ∆ b add up and enhance the bottom mass. 17 In accordance with Eq. (25), this shifts the Higgs mass downwards at the one-loop level. This effect can be seen in the shape of the red curve: while for scenarios where contribution of the b/b-sector is numerically small (see e.g. Fig. 9 in [5]) the local maximum for M h at positive values of X DR t is typically several GeV higher than the one at negative values of X DR t , for the red curve in the right plot of Fig √ 6, giving rise to an upward shift of more than 3 GeV. Since there is only a moderate splitting between M SUSY and M t , the effect of the resummation of higher-order logarithmic contributions remains relatively small (blue curve). It amounts to a downward shift of less than 1 GeV for X DR t = √ 6.
In Fig. 3, we renormalise the stop sector in the OS scheme. In the left plot X OS t = 2 is chosen. This plot shares the same features at low M SUSY as the corresponding plot in Fig. 2. However, for large M SUSY the effect of the resummation becomes more prominent due to the presence of logarithmic terms of O(m 2 b m 4 t ) that are enhanced by t 2 β (see Eq. (83) below), In the considered scenario, the resummation gives rise to a downward shift of the Higgs mass, visible as the difference between the blue curve and the green curve, by ∼ 2 GeV for M SUSY = 10 TeV and by (M SUSY ) 7.5 GeV. 16 It is worth noting that due to the parameter choices, which enhance ∆ b , different levels of approximation in ∆ b yield very different results for M h in this scenario. 17 In particular, in the considered scenario and for X DR  Hybrid calculation including: As a final phenomenological application of our improved calculation, we consider the M 125,µ− h benchmark scenario recently defined in [83], accompanying the benchmark scenarios proposed in [84,85]. In this scenario the SUSY input parameters are fixed as For the SM parameters the ones recommended by the LHC-HXSWG [86]  The stop SUSY soft-breaking parameters are defined in the OS scheme. In [83], also the sbottom trilinear coupling is renormalised in the OS scheme. For better comparison with our previous results, we instead choose to fix A b and the sbottom masses in the DR scheme. 18 In addition, we define tan β at the scale M SUSY instead of at the scale M t , which was used in [83].
Note that for this scenario µ = −2 TeV is chosen implying relatively large ∆ b corrections which enhance the cross section times branching ratio for the heavy Higgs bosons decaying to a pair of bottom quarks. In addition, the ∆ b corrections also affect the prediction for the SM-like Higgs boson, which we will investigate here. The stop mass scale is equal to 1.5 TeV, so we do not expect the resummation of logarithms controlled by the bottom Yukawa coupling to have a major numerical impact in this case (see discussion above). On the other hand, as we have seen in Figs.1-3, large ∆ b corrections imply that the prediction for M h can be sensitive to the level of accuracy in the determination of the bottom mass which is used in the fixed-order corrections at the oneand the two-loop level.
In Fig. 4 we present, in the (M A , tan β) plane, the contour lines of the SM-like Higgs boson mass ranging from 122 GeV to 125 GeV. 19 We do not consider any of the experimental constraints described in detail in [83][84][85] and concentrate only on the prediction for the mass of the lightest Higgs boson of the MSSM. The red dashed and green solid lines correspond to two different computational setups. We calculated the red contours including only the leading one-loop corrections to ∆ b of O(α s , α t ) and evaluated the bottom-quark mass according to Renormalisation scheme 1 as described in Section 2. Apart from the different definition of some of the input parameters, as mentioned above, this corresponds to the default settings of FeynHiggs-2.16.1, which was used in [83] for the analysis of the benchmark scenario. The green lines show the prediction based on the improved calculation described in this paper. In comparison to the red contours, we also include electroweak one-loop as well as the leading two-loop corrections to ∆ b , evaluate the bottom-quark mass at the SUSY scale according to Eq. (20) and resum logarithms proportional to the bottom-Yukawa coupling.
We notice that in the region of small tan β both calculations agree with each other very well since in this region the corrections from the bottom/sbottom sector are negligible. In this region the Higgs mass grows with increasing tan β mainly due to the growth of the treelevel mass. With a further increase of tan β the Higgs mass starts to decrease due to large ∆ b corrections and the rapid increase of the DR bottom mass in the MSSM. This behaviour corresponds to the one that we observed in the left plot of Fig. 1. As discussed there, the mass of the SM-like Higgs computed using FeynHiggs-2.16.1 falls faster with increasing tan β than the mass computed using the calculation presented in the current paper due to the lower accuracy level in the calculation of ∆ b of the previous result. Consequently, the tan βregion in which the SM-like Higgs mass is compatible (taking into account the theoretical uncertainties) with the experimentally measured value is enlarged. The corresponding upper bound on tan β in this scenario is shifted from ∼ 28 to ∼ 33.

EFT calculation for complex input parameters
In this Section, we discuss the numerical effect of including the full phase dependence into the two-loop threshold corrections to the Higgs self-coupling. First, let us briefly review the method used in FeynHiggs to handle non-zero phases so far. The treatment of the two-loop corrections in the presence of complex parameters is controlled by the flag tlCplxApprox. When it equals 3, the fixed-order O(α t α s , α 2 t ) corrections including the full phase dependence are activated and combined with the fixed-order O(α b α s , α b α t , α 2 b ) corrections. Since the implementation of the latter corrections up to now is based on the results of [43,44], that were obtained for the case of real parameters, an interpolation in the phases is invoked for this part of the two-loop corrections. Specifically, an interpolation is carried out in FeynHiggs when the phases of µ, M 3 , X t or X b are non-zero. The user can choose between interpolation in A t or X t , and A b or X b . In the EFT part of the code the interpolation is always carried out in the following way. First, the RGEs are integrated numerically and the subtraction terms are calculated for all possible combinations of +|P | and −|P | (where P ∈ {µ, X t /A t , M 3 }). 20 After that, linear interpolation is performed on the obtained grid. In this Section, we choose to interpolate X t in the comparison with our new results when the phase of X t or A t is non-zero.
The phases of the above-mentioned parameters enter the hybrid calculation via threshold corrections to the Higgs self-coupling and via the subtraction terms. As we mentioned in Section 3, both of them depend only on the absolute value | X t | at the one-loop level, so the interpolation would give a correct result if only LL and NLL resummation were included and the interpolation was performed in X t . However, the two-loop threshold corrections to the Higgs self-coupling (and hence the two-loop non-logarithmic subtraction terms) do not depend just on the absolute value of X t . For example, the O(α t α s ) threshold correction also depends on the cosine of the phase difference, cos(φ Xt − φ M 3 ), and the formula for the O(α 2 t ) threshold correction depends on | Y t | and cos(φ Xt − φ Yt ). In comparison to the full formula, the application of interpolation leads to deviations at the next-to-next-to-leading logarithmic order. The phases also enter the expression for the two-loop threshold corrections of the bottom Yukawa coupling and ∆ b . First, we will, however, concentrate on MSSM scenarios in which the effect of the bottom Yukawa coupling on the Higgs mass is negligible and so we will not include any two-loop corrections of O(α b α s , α b α t , α 2 b ) for the results that are presented in Fig. 5 and Fig. 6.
In order to test our approach we first consider the same MSSM scenario as in Fig. 3 of [36]: all soft SUSY breaking masses and µ are equal to the common mass scale M SUSY = 2 TeV, tan β = 10 and X DR t (M SUSY ) = √ 6. We vary the phase of the gluino mass parameter M 3 in the interval [−π, +π] and assume all the other input parameters to be real. In this way, we test the phase dependence of the O(α t α s ) threshold correction.
In the left plot of Fig. 5, we show the comparison between the pure EFT prediction of FeynHiggs-2.16.1 (red line) and our new calculation including the full phase dependence (green line). First, we notice that the two methods give the same answer for φ M 3 = 0, ±π which is expected because in these cases M 3 is a real parameter. This serves as a cross-check for our implementation. Second, we see that the interpolation in this particular scenario is a fairly good approximation: the absolute difference between the two curves does not exceed ∼ 0.3 GeV. The largest deviations occur for φ M 3 ± π 4 and φ M 3 ± 3π 4 . Since the interpolation is only performed in one parameter, φ M 3 , the resulting curve consists of two straight lines.
In the case of the hybrid calculation (see right plot of Fig. 5), the phase dependence at the two-loop level is fully included in the fixed-order part of the calculation. However, the subtraction terms are interpolated in the same way as the EFT calculation. As a consequence of those phase-dependent contributions in the fixed-order part and the subtraction terms, the curve showing the interpolated hybrid calculation (red) has a different behaviour than the interpolated EFT calculation shown in the left plot. As in the left plot of Fig. 5, the overall difference between the full hybrid and the interpolated hybrid calculation does not exceed ∼ 0.3 GeV.
Next, we proceed with a scenario which is similar to the one described above, but we assume that X t and Y t are purely imaginary while keeping the same absolute value for | X t | = √ 6 as before. As one can see in the left plot of Fig. 6, where again the result of the pure EFT calculation varying the phase of M 3 is shown, the trilinear interpolation procedure The two results do not agree for φ M 3 = 0, ±π since X t and Y t are chosen purely imaginary, and therefore an interpolation is also carried out with respect to those phases.
As a next step, we investigate the effects of the phase dependence in the O(α 2 t ) threshold correction. To enhance the numerical value of this correction, we choose a low value for tan β, namely tan β = 3. This choice, however, suppresses the tree-level Higgs mass, so to obtain a predicted value around 125 GeV we have to choose in this scenario a heavy SUSY scale of M SUSY = 20 TeV. In order to isolate the effects of the phase dependence in the considered corrections, we fix the phase of the gluino mass parameter to be equal to the phase of X t . As a consequence of this choice, the phase dependence in the O(α t α s ) threshold correction vanishes. We also choose the Higgsino mass parameter to be positive, φ µ = 0.
The EFT prediction, varying the phase φ Xt = φ M 3 , is shown in the right plot of Fig. 6. Even though we have chosen a low value of tan β = 3 in order to enhance the impact of the O(α 2 t ) threshold correction, the overall phase dependence of the full result (green) is quite small. The difference between the Higgs mass calculated at φ Xt = 0 and φ Xt = π is only ∼ 0.05 GeV. Lowering tan β even further (and pushing M SUSY higher) does not lead to a stronger phase dependence. The behaviour of the interpolated result (red) is different. As in the case of Fig. 5, we see that the results of both methods coincide for φ Xt = 0, ±π since for these three points all parameters are real. For other values of φ Xt , however, the interpolation procedure underestimates the value of M h predicted based on the full expression by up to ∼ 0.5 GeV.
This large deviation can be understood by looking at Fig. 7 showing the M h prediction of the EFT calculation including the full phase dependence. The same scenario as in the right plot of Fig. 6 is used, but φ Xt and φ M 3 are varied independently. As visible in the plot, the contours are almost diagonal due to the small phase dependence of the O(α 2 t ) threshold corrections. The parabola-like shape of the interpolated result, as visible for the red curve in the right plot of Fig. 6, is a consequence of the bilinear interpolation in φ Xt and φ M 3 . For φ Xt = φ M 3 > 0, the Higgs mass values at (φ Xt , φ M 3 ) = (0, 0), (0, π), (π, 0), (π, π) enter the interpolation procedure. For the values (φ Xt , φ M 3 ) = (0, π), (π, 0) the phase dependence of the O(α t α s ) threshold correction is picked up resulting in the large phase dependence observed for the red curve in the right plot of Fig. 6. In the considered case, an interpolation in φ M 3 = φ Xt rather than in φ M 3 and φ Xt separately would improve the quality of the interpolation.
It should be noted that for the hybrid result the difference between the EFT result incorporating the full phase dependence and the one based on the interpolation, shown in the right plot of Fig. 6, is further enhanced because of the different treatment of the phase dependence in the fixed-order contribution and the subtraction terms. As a consequence, in this extreme scenario, the incomplete cancellation between the corresponding terms in the fixed-order part and the subtraction terms leads to an artificial enhancement of the deviation that can amount up to ∼ 2 GeV.
As a final topic in this Section, we analyse the interplay between the resummation of the logarithms proportional to the bottom Yukawa coupling and the inclusion of the full phase dependence into the EFT part of our hybrid calculation. As a starting point we go back to the scenario discussed in Section 5.1. Namely, we consider a single scale scenario, where all soft-breaking masses as well as the mass of the charged Higgs boson 21 are equal to 1.5 TeV, A DR b = 2.5M SUSY , the Higgsino mass parameter is negative, µ = −M SUSY , the bino and wino masses are chosen to be positive, M 1,2 > 0, and X OS t = 2. The phase of the gluino mass parameter is a free parameter, and we vary it in the interval from −π to +π. We examine The behaviour as a function of φ M 3 is significantly different for tan β = 45. The red solid curve starts to grow when φ M 3 increases starting from −π, resembling the red dashed line in shape. However, it reaches a maximum value at φ M 3 − π 3 . This is a consequence of the fact that the ∆ b correction becomes important in this region, leading to a steep increase of the MSSM bottom mass (see right plot of Fig. 8). Thus, the one-loop corrections involving the bottom mass (see Eq. (25)) become important, giving rise to a downward shift in M h . At φ M 3 = 0 the bottom mass reaches ∼ 5.8 GeV, and the Higgs mass prediction has a minimum at ∼ 123 GeV. The point φ M 3 = 0 in this plot corresponds to the point where M SUSY = 1.5 TeV in the left plot of Fig. 3. As in Fig. 3, we observe that the inclusion of the two-loop fixed-order corrections controlled by the bottom Yukawa coupling (the difference between the red and the green curves) has a very significant effect. The resummation of higher-order logarithmic contributions (the difference between the blue and the green curves) leads to a downward shift of ∼ 1 GeV for φ M 3 ± π 3 and of ∼ 1.2 GeV for φ M 3 = 0. The results displayed in Fig. 8 demonstrate that the (s)bottom sector contributions can have an important impact on the phase dependence. Similarly to Fig. 3, we again find that the resummation of logarithms proportional to the bottom-Yukawa coupling amount to an O(1 GeV) effect for large tan β if the OS scheme is used for the renormalization of the stop sector.

N 3 LL resummation
Here, we study the numerical effects of including N 3 LL resummation at leading order in the strong gauge coupling (see Section 4) into our hybrid framework. We study a simple singlescale scenario in which all non-SM masses are set to the common scale M SUSY . Furthermore, we set all trilinear soft SUSY-breaking couplings, except for A t , to zero. We define the stop parameters in the DR scheme at the scale M SUSY . We set tan β = 10.
In Fig. 9, we compare the results obtained using three different accuracy levels to each other: NNLL resummation with the SM top Yukawa coupling extracted at the two-loop level, NNLL resummation with the SM top Yukawa coupling extracted at the three-loop level and N 3 LL resummation, which also involves the SM top-Yukawa coupling extracted at the threeloop level. In the upper plots, the different results are shown as a function of M SUSY . In the upper left plot, the three results (blue, red and green lines) are shown for vanishing stop mixing (solid lines) and for X t = − √ 6 (dashed lines). For vanishing stop mixing, all three results are in good agreement with each other for low M SUSY . If M SUSY is raised, there is, however, an increasing difference between the NNLL result (with the two-loop level SM top Yukawa coupling) and the two results involving the three-loop level SM top-Yukawa coupling of up to ∼ 1 GeV for M SUSY ∼ 100 TeV. This shift is almost completely caused by including the three-loop corrections to the extraction of the SM top Yukawa coupling, since the NNLL result with the SM top Yukawa coupling extracted at the three-loop level and the N 3 LL result are in very good agreement also for M SUSY ∼ 100 TeV. Also for X t = − √ 6, the NNLL result with the SM top Yukawa coupling extracted at the three-loop level and the N 3 LL result are in good agreement across the considered M SUSY range (within ∼ 0.3 GeV). This difference is displayed by the red curve in the upper right plot of Fig. 9. The NNLL result with the SM top Yukawa coupling extracted at the two-loop level deviates from the other two results by ∼ 0.7 GeV, as shown by the blue curve in the top right plot of Fig. 9. In this plot, furthermore the estimate of the uncertainty associated with the truncation error in the calculation of the O(α t α 2 s ) threshold correction for the Higgs self-coupling, obtained mainly by switching between different parameterizations of the top-Yukawa coupling (see [18] for more details), is shown as a green band. We find that this estimate is of the same size as the shift induced by including the O(α t α 2 s ) threshold correction. In the lower plots of Fig. 9, the same quantities as in the upper plots are shown, but M SUSY is set to 5 TeV and X t is varied. The shifts between the various results are only mildly dependent on X t (varying X t leads to shifts of up to 0.4 GeV). This dependence tan β = 10, would be stronger for lower M SUSY values. The estimate of uncertainty associated with the truncation error, however, shows a strong dependence on X t . Whereas it is negligible for −1 X t 1, it increases to up to 0.7 GeV for | X t | ∼ 3.5. As shown by the red curve in the lower right plot of Fig. 9, the difference between the NNLL result with the SM top-Yukawa coupling extracted at the three-loop level and the N 3 LL result is rather small except for large negative values of X t . Where this difference exceeds the level of 0.2 GeV, it is smaller than the estimated uncertainty of the truncation error.
As expected, the results for the N 3 LL resummation are in very good agreement with the results of [18]. We observe that the main part of the shift induced by including N 3 LL resummation is caused by taking into account the three-loop corrections to the extraction of the SM MS top Yukawa coupling from the measured top mass. The shift caused by including the O(α t α 2 s ) threshold correction for the Higgs self-coupling is smaller and also associated with a rather large uncertainty for large | X t | values. For small | X t | values, the shift induced by including the O(α t α 2 s ) threshold correction for the Higgs self-coupling is found to be very small. Therefore, we choose in our implementation to use the result obtained using NNLL resummation with the SM top Yukawa coupling extracted at the three-loop level as default result until the uncertainty in the calculation of the O(α t α 2 s ) threshold correction is further reduced by incorporating additional higher-order contributions.

Conclusions
In this paper, we have presented an improved prediction for the lightest Higgs boson mass in the MSSM in scenarios with large tan β, complex input parameters and large M SUSY . Our calculation builds on results that are contained in the publicly available code FeynHiggs.
The first improvement concerning scenarios with large tan β includes the change of the renormalisation scheme for the bottom mass with respect to the present implementation in FeynHiggs: instead of treating the bottom mass as a derived quantity, in the scheme used in our calculation it is as an independent parameter, renormalised in the DR scheme in the full MSSM at scale M SUSY . The scheme that we have adopted yields numerically more stable results and turned out to be better suited for the combination with the EFT calculation. In the calculation of the DR bottom mass, we have taken into account higherorder corrections enhanced by tan β by means of a resummation of the quantity ∆ b . We have incorporated full one-loop corrections to ∆ b . Moreover, we have adapted the leading two-loop QCD corrections to ∆ b obtained in [37][38][39] such that they are suitable for the framework of our calculation. The inclusion of this correction is formally a three-loop effect. While this correction is numerically not relevant for large parts of the parameter space, it can become sizeable for scenarios with large tan β.
Moreover, we have included one-and two-loop threshold corrections to the SM Higgs self-coupling proportional to the bottom Yukawa coupling well as the corresponding RGE contributions up to the three-loop level. This allows resummation up to the next-to-nextto-leading-order. In contrast to the resummation of the logarithms proportional to the top Yukawa coupling or electroweak couplings, here the one-and two-loop leading logarithms are numerically negligible due to the smallness of the bottom mass. However, at the two-loop level for the case where the stop sector is renormalised in the OS scheme the next-to-leading logarithms become parametrically enhanced for large tan β. In this case, the resummation can become numerically relevant for large tan β and large M SUSY .
Secondly, we used the two-loop fixed order results presented in Refs. [8,10] to derive twoloop threshold corrections to the SM Higgs self-coupling for the matching between the SM and the MSSM that are valid for the general case of complex input parameters. This enabled us to perform the EFT calculations for the case of complex parameters. We compared the results including the full phase dependence to the results obtained by the use of the interpolation routine that has been adopted in FeynHiggs up to now. For the pure EFT calculation, we have found the interpolation procedure to perform well in scenarios with only one non-zero phase. In scenarios with more than one non-zero phase, we observed deviations in the prediction for M h of up to 1 GeV. For the hybrid result the incorporation of the full phase dependence of the EFT part of the calculation yields another important improvement. Up to now the corresponding contributions in the fixed-order result (containing the full phase dependence) and the subtraction terms (based on the interpolated EFT contributions) were treated differently, which could lead to an incomplete cancellation between the two types of contributions. This can lead to numerical deviations of up to 2 GeV compared to our improved result where the treatment of the phase dependence is the same in all parts of the calculation. We furthermore analysed the interplay between the resummation of the logarithms proportional to the bottom Yukawa coupling and the inclusion of the full phase dependence into the EFT part of the hybrid calculation. We have found that the impact of phase variations on the prediction for M h can be modified very significantly through the contributions of the b/b sector.
Finally, we combined the publicly available code Himalaya with FeynHiggs in order to obtain a prediction for M h including N 3 LL resummation at leading order in the strong gauge coupling. A similar analysis was performed in [29], and we find a very good agreement with the results presented in that paper. The overall effect of the N 3 LL resummation is 1 GeV, and it only weakly depends on M SUSY . We have found that employing the extraction of the SM top Yukawa coupling at the three-loop level within the existing NNLL hybrid calculation yields a result that approximates the N 3 LL resummation well in view of the remaining theoretical uncertainties of the N 3 LL contribution.
The improvements described in this paper will be implemented into an upcoming version of the public code FeynHiggs.

A Derivation of two-loop threshold corrections
In this Appendix, we derive the two-loop threshold corrections for the SM Higgs self-coupling for the matching between the SM and the MSSM based on the fixed-order calculations presented in [8,10,42]. We fully take into account the dependence on CP-violating phases.
The threshold corrections to the quartic coupling λ can be obtained via the matching of the four-point vertex function involving the SM Higgs boson as external particle. Here, however, we follow a different approach. Since in the SM the running mass of the lightest Higgs boson is related to its quartic coupling via the threshold corrections to λ can be obtained via the threshold correction to the running Higgs mass m h . 22 Below we outline the method and derive the general formulas for the oneand two-loop threshold corrections to λ in the gaugeless limit. Similar methods can be found in [14,30].
In the limit M A M t , the SM-like Higgs pole mass in the MSSM up to the two-loop level is given by 23 where m h is the MSSM tree-level mass, and the prime indicates the derivative with respect to the external momentum squared. In the gaugeless limit and the decoupling limit (M A M t ), m h = 0 can be inserted. All parameters entering the self-energies Σ Below the matching scale Q, the effective field theory is the SM. We write the matching condition for the SM running Higgs mass m 2 h as a loop expansion, where the ellipsis denotes three-loop terms and higher. Since m h in Eq. (28) equals zero in the considered approximation, we have m h,tree = 0. The pole mass in the SM can then be obtained via the solution of the pole equation where Σ MS,SM hh is the SM Higgs boson self-energy renormalised in the MS scheme with the tadpoles renormalised to zero. Since the SM is treated as an effective field theory, its parameters are related to the corresponding parameters in the MSSM. This relation can be 22 This method is not sufficient to obtain the threshold corrections for all quartic couplings if the EFT below M SUSY is the Two-Higgs-Doublet-Model [87]. 23 In general, only the real part of each term in the sum on the right hand side of Eq. (28) should be considered, since the Higgs self-energies have imaginary parts arising from the contributions of the particles which are lighter than the SM-like Higgs. Since in the MSSM the mass of the SM-like Higgs is close to the electroweak scale, in the usually considered scenarios these imaginary parts arise only from SM particles and, therefore, cancel out in the matching procedure. schematically written as follows, where P is a coupling constant, a running quark mass, or the vacuum expectation value. Inserting this relation into the self-energy Σ MS,SM hh (M 2 h ) induces a shift at one order higher in the loop expansion, The first term on the right-hand side of this equation represents the self-energy which has the same analytic form as Σ MS,SM hh but with all MS SM coupling constants and masses being replaced with their DR MSSM counterparts. Therefore, the two self-energies are equal at the one-loop level (the different regularisation does not lead to a different result in this case). The difference between the two self-energies is encoded in the quantity Σ SM,shifts hh which is of two-loop order and higher.
The renormalised self-energy of the SM-like Higgs boson in the full MSSM can be split into parts, where the SM part contains contributions from the diagrams with only SM particles and the non-SM part (indicated as "n/SM") originates from the diagrams with at least one non-SM particle. At the one-loop level, the following identity holds, where in this equation the symbols " " and " " are used to denote the SM part of the MSSM Higgs self-energy renormalised in the DR scheme and the SM Higgs self-energy renormalised in the MS scheme, respectively. This equation means that the SM contributions in the full MSSM self-energy computed in the DR scheme, Σ SM hh , have the same analytic form as the self-energy computed in the SM in the MS scheme, Σ SM hh . In Eq. (34) the replacement rule on the right-hand side, P SM → P MSSM , implies that all MS SM parameters in the SM selfenergy have to be replaced by their DR MSSM counterparts without additional shifts. At the two-loop level, a relation analogous to Eq. (34) holds for the SM-type corrections to the Higgs mass proportional to the Yukawa couplings, i.e. the corrections of O(α 2 t , α t α b , α 2 b ) without including any parameter shifts in the one-loop self-energy. For the mixed Yukawa-QCD corrections of O(α t α s , α b α s ), the two different choices of the regularisation scheme (dimensional regularisation in the case of the SM and dimensional reduction in the case of the MSSM) lead to different expressions for the SM part of the self-energy. This can already be anticipated since the running MS and DR quark masses are not equal to each other at the one-loop level (see e.g. [81]). By using TwoCalc [88,89] and the scripts described in [67] we explicitly checked the following relation, and the analogous relation for the O(α b α s ) self-energies. In Eq. (35), 24 As explained above, the one-loop reparametrisation of the couplings and masses in the oneloop SM self-energy induces shifts at the two-loop order, where P are all SM parameters which enter Σ Here we have exploited the fact that in the gaugeless limit Σ SM,(1) hh scales as ∝ 1/v 2 . In [26] it was shown that in the heavy SUSY limit the non-SM part of the relative shift in the vacuum expectation value can be expressed via the non-SM part of the Higgs self-energy derivative, 24 Here we do not specify the renormalisation scheme for the strong coupling and for the top mass since a change of the renormalisation scheme is of three-loop order. 25 Note that the one-loop shift ∆ (1) v in this equation is the same quantity as ∆v in Eq. (18) in Section 2.4 for Q = M SUSY . 26 For brevity, we will omit arguments of the self-energies in the rest of this section, implying that they always equal zero.
Using this relation, Eq. (40) can be rewritten as follows, At the matching scale Q, the predictions for the physical Higgs mass, M h , in the SM and the MSSM have to be equal order by order, By equating the one-loop pieces in Eq. (44) and taking into account Eq. (28) and Eq. (43) we get, By equating the two-loop pieces in Eq. (44) and expanding the one-loop self-energies of the Higgs boson in the full MSSM according to Eq. (33), we get The running Higgs-boson mass m 2 h can be related to the threshold corrections to the quartic coupling λ via the relation To express the one-and two-loop corrections in terms of the MSSM coupling constants we have to perform the shift of the vacuum expectation value in Eq. (48), By solving this equation at the one-and two-loop levels, we obtain the expressions for the matching coefficients for the quartic coupling, As already mentioned, in the expressions above all couplings are DR MSSM couplings at the scale Q. Another option is to parametrise these threshold corrections in terms of the MS SM couplings at Q. In this work, we will use the MS top-Yukawa coupling in the SM and the DR MSSM bottom-Yukawa coupling to parametrise the one-and two-loop threshold corrections.
To express the two-loop threshold corrections in terms of the SM MS top-Yukawa coupling we have to reparametrise the top mass and the vacuum expectation value in the one-loop O(α t ) threshold correction. This generates the following two-loop terms, which has to be added to Eq. (50b).
We have evaluated the non-SM part of the one-loop Higgs boson self-energy with the help of FeynArts [90][91][92] and FormCalc [93] and then expanded in the limit mt L , mt R , mb R m t , m b . The explicit expressions for them (and their derivatives) in the gaugeless limit and for the case mt L = mt R = mb From these expressions and Eq. (50a) it is clear how the one-loop threshold corrections to λ computed in [16,69] can be generalised to the case of the MSSM with complex parameters. These corrections are polynomials in the squark mixing parameter X q . To obtain the expression in the MSSM with complex parameters X q has to replaced by | X q |.
The two-loop self-energies were taken from [8,10,42] and expanded in the limit without any additional assumptions on the internal masses and the phases of X t , X b , µ and M 3 . The two-loop SM self-energies in the MS scheme were taken from [55,56] and extracted from the code FlexibleSUSY [24,27,54]. Finally, the one-loop shifts ∆ (1) m t and ∆ (1) m b have been computed using FeynArts and FormCalc. The resulting two-loop formulas for the threshold corrections are presented in App. B.2.
B Threshold corrections for the case of non-vanishing CP-violating phases

B.1 One-loop threshold corrections
If the sfermions and heavy Higgses are integrated out from the MSSM, effective Higgsgaugino-Higgsino couplings,g 1u,1d,2u,2d , are generated (for their exact definition see e.g. [69]). In principle, they can be complex. An explicit calculation of their matching conditions at the SUSY scale M SUSY = √ mt L mt R , however, shows that they remain real if the sfermions and heavy Higgses are integrated out. All other couplings of the EFT below the SUSY scale are also real-valued. 27 The only exception are the mass parameters of the EWinos for the case of light EWinos. The phases of these parameters become relevant if the EWinos are integrated out at the EWino mass scale, M χ = |M 2 ||µ|, and the SM is recovered as EFT. 28 The threshold corrections of the top and bottom Yukawa couplings originate only from the corrections to the external Higgs leg. They read with In the expressions above, κ = 1/(4π) 2 is used to indicate the loop order, and ∆ WFR are the one-loop corrections to the external Higgs legs. Setting the phases to zero, we recover the result presented in Ref. [69]. The loop functions f and g are defined in the Appendix of Ref. [69]. Similarly, also the matching condition of the Higgs self-coupling is modified, The CKM phase is neglected here. 28 Here, we assume that the absolute values of the EWino mass parameters M 1 , M 2 and µ are of the same order of magnitude. − 4g 1ug1dg2ug2d − 2 g 2 1u +g 2 2d g 2 1d +g 2 2u ln |µ| 2 M 2 r 1 r 2 r 1 + r 2 f 8 (r 1 , r 2 ) + 5 2 h 6 (r 1 , r 2 ) The loop functions f i are defined in the Appendix of Ref. [69]. The loop functions h i are defined by h 2 (r) = 6 11(1 − r 2 ) 3 2 + 3r 2 − 4r 4 − r 6 + r 2 (4 + 5r 2 − r 4 ) ln r 2 , h 6 (r 1 , r 2 ) = 6 5(1 − r 2 1 ) 2 (1 − r 2 2 ) 2 (r 2 1 − r 2 2 ) − (1 − r 2 1 )(1 − r 2 2 )(r 4 2 − r 2 1 r 2 2 − r 4 1 + r 4 1 r 2 2 ) + (1 − r 2 2 ) 2 r 6 1 ln r 2 1 − (1 − r 2 1 ) 2 r 6 2 ln r 2 2 . (60f) In the limit of r, r 1 , r 2 → 1 all of the loop functions h i approach 1. Setting all phases to zero, we again recover the expression given in Ref. [69]. The corresponding expressions for the EWino contribution to the matching between the SM and the MSSM can be obtained by replacing the effective Higgs-Higgsino-gaugino couplings g 1u,1d,2u,2d in Eq. (59) using their tree-level matching conditions [94], (61a) g 2d = g cos β, (61b) g 1u = g sin β, (61c) g 2u = g sin β. (61d) The functions F 5,6,8,9 (x) are defined in Appendix A of Ref. [69]. Q is the renormalisation scale, which we set equal to M SUSY . We neglect electroweak contributions to ∆v.
In addition, we give here the result for the one-loop threshold correction for the top Yukawa coupling, which can be used to reexpress the two-loop threshold corrections of the Higgs self-coupling in terms of the MSSM top Yukawa coupling (see Section B.2). It is given by Setting all phases to zero, we again recover the expression given in Ref. [69].

B.2 Two-loop threshold corrections
Here, we list the two-loop threshold corrections to the Higgs self-coupling in the limit where all involved non-SM particles except for EWinos and gluinos have the same mass, The functions f 1,2,3 (x) are the same as the functions f 1,2,3 (x) from Ref. [95], . We use a different notation for them than in Ref. [95], since we have already used the notation f i (x) for the functions in the threshold correction to the quartic coupling above in Eq. (59). 29 The constant K is Fully general expressions for the two-loop threshold corrections to the SM Higgs self-coupling can be found in ancillary files distributed alongside this paper.
C Dependence of ∆ 2l b on CP-violating phases Before deriving the phase dependence of the two-loop correction to ∆ b , we first consider the one-loop correction. The O(α s , α t ) one-loop ∆ b correction in the heavy SUSY limit can be obtained by evaluating the diagrams in Fig. 10. 30 These diagrams include the incoming and the outgoing b quarks with different chirality. The diagrams which involve quarks with the same chirality are subleading with respect to their powers of tan β and do not contribute to ∆ b [63]. The diagrams are drawn using the chiral basis in which the "left" and the "right" squarks propagate and the off-diagonal mass term is interpreted as an additional interaction (denoted as ) which flips the chirality quantum number of the squark. In the limit M SUSY v only the diagrams with a single mass insertion contribute. This can be seen in the following way. The left diagram in Fig. 10 is proportional to 29 We note that these functions are not independent from each other. For example, using the identities for the Spence function Li 2 (x 2 ) one can show that f 3 (x) = (1 − 2x 2 − 2x 4 ) f 5 (x). However, we decided to stick to the notations of Ref. [95]. So, we expressed our result in terms of f 1,2,3 and added two more functions for better readability of the results. 30 All diagrams in this Section were produced with Axodraw [96]. where C 0 is a Passarino-Veltman function corresponding to the scalar vertex function with three external legs. If all soft SUSY-breaking masses are equal to M SUSY , the expression in Eq. (70) reduces to The diagram with two mass insertions, which is only possible if the two external quarks have the same chirality, is proportional to where D 0 is a Passarino-Veltman function corresponding to the scalar vertex function with four external legs. If all soft SUSY-breaking masses are equal to M SUSY this diagrams behaves like We see that it is suppressed by an additional factor m b /M SUSY compared to the diagram with one insertion and is therefore subleading. This is consistent with the fact that only diagrams with external legs of different chirality contribute to ∆ b . Clearly, diagrams with more insertions will be suppressed by additional factors of m b /M SUSY . Following similar arguments, one can show that diagrams similar to the right diagram in Fig. 10 with more than one mass insertion are suppressed by powers of m t /M SUSY . The same kind of argument applies to higher-order corrections to ∆ b [38] as can be proven by using the Kinoshita-Lee-Nauenberg theorem [97,98]. Namely, the diagrams contributing to the two-loop quantity ∆ b of order O(α 2 s , α s α t ) contain only one mass insertion. The phases of the complex parameters A t , µ and M 3 enter the diagrams of Fig. 10 through the mass insertion and through (b L b * L g), and (b L t * R H + 2 ) vertices. In particular, the mass insertion in the left diagram in Fig. 10 yields the phase factor ∝ e +iφµ while both vertices contain e +i φ M 3 2 . The overall diagram is then ∝ e i(φ M 3 +φµ) . The result for the analogous diagram with incoming b R and outgoing b L leads to the phase factor ∝ e −i(φ M 3 +φµ) . The overall phase dependence is then ∝ cos(φ M 3 + φ µ ) (see Eq. (62) in App. B.1). The mass insertion in the right diagram in Fig. 10 gives the phase factor ∝ e +iφ A t , while the (b L t * R H + 2 ) and (b R t L H − 2 ) vertices contain the entries of the chargino mixing matrices V * 22 and U * 22 . In the gaugeless limit, they are proportional to the phase factors ∝ e +i φµ 2 . The overall diagram (together with its complex conjugated) is ∝ cos(φ At + φ µ ) (see Eq. (62) in App. B.1). The two-loop diagrams contributing to the quantity ∆ b at O(α 2 s ) can be split into three categories: either a gluon, a sbottom or a gluino is added to the one-loop O(α s ) graph. Examples of the corresponding diagrams are depicted in Figs. 11-13. The particles which are added to the one-loop graph are highlighted with red color. 31 Following the argumentats given at the end of Section 3, we can conclude that in the case of the MSSM with complex parameters the two-loop ∆ b of O(α 2 s ) given in Eq. (23) has to be multiplied by cos(φ µ + φ M 3 ). The same reasoning can be applied to the O(α s α t ) ∆ b corrections: Eq. (24) has to be multiplied by cos(φ µ + φ At ). 32

D Leading and next-to-leading logarithms
In this part of the Appendix we present analytic expressions for the leading (LL) and nextto-leading (NLL) logarithms proportional to the bottom Yukawa coupling, which appear at the one and two-loop order in the calculation of the lightest Higgs boson mass in the MSSM. They are derived in the special case of 31 The rightmost diagram in Fig. 12 cannot be reduced to the left diagram in Fig. 10. This fact, however, does not change our arguments that we present here. 32 Corresponding two-loop diagrams can be found in Appendix C of [38].  In the following expressions, κ = 1/(4π) 2 , L is the logarithm of the ratio of M SUSY and M t , m t is MS SM top mass, m t ≡ m MS,SM t (m t ), m b is the DR MSSM bottom mass at the scale M SUSY defined in Eq. (20), v is the SM Higgs vacuum expectation value v ≡ v G F = (2 √ 2G F ) −1/2 174 GeV. The ratio of the vacuum expectation values of the two Higgs doublets, tan β, is renormalised in the DR scheme at the scale M SUSY . g 3 is the strong gauge coupling, g 2 3 = 4πα s . We do not specify the renormalisation prescription for it, since it appears in the expressions for M h starting at the two-loop level. Therefore, a change of the renormalisation scheme for g 3 is a three-loop effect. The parameter X b is renormalised in the DR scheme at the scale M SUSY , while X t is fixed either in the OS scheme or in the DR scheme at the scale M SUSY .
The logarithmic terms can be derived in one of the two following ways. The first method is an approximate perturbative solution of the renormalisation group equation (RGE) for the quartic coupling λ dλ d log Q 2 = κβ where Q is a renormalisation scale, β λ , β (2) λ are the one-and two-loop contributions to the beta function, and the ellipsis encodes higher-order terms. Truncating the result at the second order, one obtains [82,99] λ(M t ) = λ(M SUSY ) − β The pole Higgs mass is then calculated from the pole equation [26] ( The second method is an iterative solution of the Higgs boson pole mass equation in the MSSM which in the decoupling limit M A M Z reads [26] where the prime stands for the derivative of the self-energy with respect to the momentum squared, and m 2 h is the lightest Higgs boson mass at the tree level. We have checked analytically that the logarithmic terms are the same in (M h ) 2 FO and (M h ) 2 EFT . This is an important cross-check of our calculation. 33 At the one-loop level, we obtain the following contribution proportional to the bottom quark mass, At the two-loop level the leading logarithmic terms read If X t is renormalised in the OS scheme, the sub-leading logarithms read If X t is renormalised in the DR scheme at the scale M SUSY , we obtain We see that the choice of the renormalisation scheme for X t affects only the logarithmic terms proportional to m 2 b m 4 t .