Improved MSSM Higgs mass calculation using the 3-loop FlexibleEFTHiggs approach including $x_t$-resummation

We present an improved calculation of the light CP-even Higgs boson pole mass in the MSSM based on the FlexibleEFTHiggs hybrid method. The calculation resums large logarithms to all orders and includes power-suppressed terms at fixed order. It uses state-of-the-art 2- and 3-loop matching of the quartic Higgs coupling and renormalization group running up to 4-loop, resulting in a resummation of large logarithmic corrections up to N$^3$LL level. A conceptually novel ingredient is the expansion of the matching conditions in terms of high-scale MSSM parameters instead of SM parameters. In this way leading terms in the stop-mixing parameter are effectively resummed, leading to an improved numerical convergence of the perturbative expansion. Furthermore, the avoidance of double counting of loop corrections is more transparent than in other approaches and more independent of the high-scale model. We present numerical results and a detailed discussion of theoretical uncertainties for standard benchmark scenarios.


Introduction
Since the discovery of the Higgs boson [1,2], the Higgs boson mass M h = (125.10±0.14) GeV [3,4] has become a high-precision observable [5], which represents another useful tool to search for physics beyond the Standard Model (SM) and constrain the large zoo of proposed SM extensions, such as supersymmetric (SUSY) models. The latter are particularly interesting, as they require the existence of scalar fields and predict the quartic Higgs coupling and thus the Higgs boson mass. The precise prediction of the SM-like Higgs boson mass in the Minimal Supersymmetric Standard Model (MSSM), however, is a long-standing challenge, because in viable MSSM scenarios large radiative loop corrections of the order ∆m 2 h ∼ (100 GeV) 2 are required, resulting in a large truncation error of the perturbation series.
There are two main mechanisms which can generate such large loop corrections: (i) Large SUSY masses M S (in particular stop masses) lead to large logarithmic corrections of the form log (M S /v), where v represents the electroweak scale. (ii) A large mixing in the stop sector, governed by the parameter X t , leads to power corrections of the order (X t /M S ) n . Effective field theory (EFT) techniques are a well-known tool to perform a resummation of the large logarithmic corrections, thus effectively avoiding a large truncation error of the contributions from mechanism (i). Concerning the (X t /M S ) n power corrections, however, no similar resummation technique has been used so far. In the present work we present a technique to effectively resum leading terms in X t in the prediction of the light Higgs boson mass.
There are different approaches to calculate the Higgs boson mass in supersymmetric models, which can be roughly classified into fixed-order   1 , EFT [40][41][42][43][44][45][46][47][48][49][50][51], and hybrid [52][53][54][55][56][57][58][59][60][61] approaches, which combine the virtues of the former two. Fixed-order approaches truncate the perturbation series at a certain order in loops and couplings, neglecting in particular large logarithmic corrections arising at higher orders. Thus, when M S ≫ v, the fixed-order approaches usually suffer from a large uncertainty due to missing large higher-order corrections. EFT approaches, on the other hand, resum the large logarithmic corrections to all orders, but usually neglect terms of the order v 2 /M 2 S . As a consequence, EFT approaches become imprecise when M S ∼ v.
Hybrid approaches combine the virtues of fixed-order and EFT calculations: They resum large logarithmic corrections to all orders and include terms suppressed by v 2 /M 2 S at fixed order. A first variant of such a hybrid approach was presented in ref. [52] and implemented into FeynHiggs. This approach uses a "subtraction method", where the large logarithmic corrections are subtracted from a fixed-order calculation and are replaced by resummed logarithms, avoiding double counting. This method was refined in refs. [53,57] and applied in the context of the DR ′ scheme in ref. [60].
An alternative way to realize a hybrid approach was presented in refs. [54,56]. This socalled FlexibleEFTHiggs approach is an EFT calculation in which the matching condition is suitably modified such that terms suppressed by powers of v 2 /M 2 S are included in the quartic Higgs coupling. One advantage of this method is the structural simplicity of the matching condition. As a result, the method is well suited for automation and has thus been implemented into the generic spectrum generators FlexibleSUSY [56,62] and SARAH/SPheno [55]. A difficulty of the FlexibleEFTHiggs approach is to make sure that large logarithms cancel in the matching between the EFT and the UV model, as required. Indeed, avoiding double counting leads to significant complications in all hybrid calculations [56,57,63]. In this paper we present an extension of the FlexibleEFTHiggs hybrid approach with a matching of the quartic Higgs couplingλ beyond 1-loop level (next-to-leading order, NLO) and apply it to perform a state-of-the-art hybrid calculation of the light CP-even Higgs boson mass in the real MSSM. Thereby our calculation incorporates several conceptual changes and significant improvements: • We parametrize the matching calculation at the high-energy scale in terms of parameters of the UV model (i.e. the MSSM). This is in contrast to the usually chosen parametrization in terms of EFT parameters. Our "full-model parametrization" has several significant advantages. An important advantage is that the cancellation of large logarithmic corrections in the matching is more transparent. Furthermore, our parametrization allows for a computer algebraic implementation which is to a large extent independent of the chosen UV model. This fact enables the straightforward application to a large class of SUSY models. The detailed discussion of the different possible parametrizations is presented in section 3.
• In our application to the MSSM we include the state-of-the-art radiative corrections in the matching up to the 3-loop level at O(1ℓ+g 2 3 (y 4 t +y 4 b )+(y 2 t +y 2 b ) 3 +(y 2 t +y 2 τ ) 3 +g 4 3 y 4 t ) and perform renormalization-group running up to 4-loop level in QCD. As a result, our calculation reaches a precision of N 3 LO with a resummation of N 3 LL, comparable with the calculation presented in ref. [60]. The details of the matching of the MSSM to the SM are presented in section 4, and numerical results are shown in sections 7-8.
• The most important advantage of our new approach and the chosen full-model parametrization is the effective resummation of QCD-enhanced terms leading in the stop mixing parameter X t , which is presented in section 5. More specifically, we show that the highest power contributions at O(y 4 t g 2n 3 , y 2 t g 2 1,2 g 2n 3 ) for all n > 0 are captured by our procedure. As a result, the perturbation expansion of the Higgs boson mass in terms of the MSSM parameters stabilizes significantly for large X t , leading to a reduced theory uncertainty of the prediction.
We begin with a recap of the SM and the MSSM in section 2, introducing our conventions. section 3 gives a general overview of the implementation of the EFT approach, discussing in particular the role of the parametrization. Our new realization of the FlexibleEFTHiggs approach within a numerical code is discussed in section 4. In section 5 we show how our chosen parametrization in terms of MSSM parameters results in a resummation of highest power X t contributions as described above. After a study of the numerical results of our new calculation in section 7, we perform a thorough analysis of the remaining theory uncertainty of our calculation in section 8.

Definition of the Standard Model and the MSSM
In the following we will denote the Standard Model (SM) parameters, defined in the MS scheme, asP = {ĝ 1 ,ĝ 2 ,ĝ 3 ,ŷ t ,ŷ b ,ŷ τ ,λ,v}, (2.1) whereĝ 1 = 5/3ĝ Y andĝ Y ,ĝ 2 andĝ 3 denote the gauge couplings of the gauge groups U (1) Y , SU (2) L and SU (3) C , respectively. The Yukawa couplings of the top quark, bottom quark and tau lepton are denoted asŷ t ,ŷ b andŷ τ , respectively. The 1 st and 2 nd generation Yukawa couplings as well as CP-violation effects are neglected and we will set the CKM and PMNS matrices to unity. The quartic couplingλ of the SM Higgs field Φ is defined by the Higgs potential We decompose the Higgs field as where h is the SM Higgs particle,v ≡ √ 2⟨Φ⟩ is the Higgs vacuum expectation value (VEV) (i.e. the minimum of the SM effective potential) which satisfiesv = (−2μ 2 /λ) 1/2 ≈ 246 GeV at tree level and G 0,± are the SM Goldstone bosons. After spontaneous electroweak symmetry breaking, the MS masses for the top, bottom and tau fermion and for the heavy physical bosons are given bŷ For convenience we define in addition the following symbols: (2.6) We denote the corresponding relevant parameters of the (R-parity conserving) Minimal Supersymmetric Standard Model (MSSM), defined in the DR ′ scheme, as P = {g 1 , g 2 , g 3 , y t , y b , y τ , v}, where g 1 = 5/3 g Y and v = (v 2 u + v 2 d ) 1/2 , whereas v u and v d denote the VEVs of the upand down-type Higgs fields which represent the minimum of the effective potential in the MSSM, (2.8) If not stated otherwise, we define tan β = v u /v d . After the spontaneous electroweak symmetry breaking in the MSSM, the DR ′ masses for the top, bottom and tau fermion as well as the SM-like Higgs in the decoupling limit are given by (2.10) We neglect inter-generation sfermion mixing, so the DR ′ masses of the stops, sbottoms and staus are given by the eigenvalues of the mass matrices and m 2 e 3 denote the squared soft-breaking mass parameters of the left-and right-handed 3 rd generation squarks and sleptons and electroweak contributions from D-terms have been omitted. 2 The sfermion mixing parameters X t , X b and X τ are defined as where A f (f = t, b, τ ) are the trilinear Higgs-sfermion-sfermion couplings and µ is a MSSM superpotential parameter. For convenience we define in addition the following symbols: where mt i denotes the i th DR ′ stop mass.

Matching procedure in general
We begin by recalling a few basic aspects of the effective field theory approach to compute weak-scale observables, such as the pole mass of the Higgs boson, M h , in SUSY models in scenarios where the SUSY scale Λ is significantly larger than the weak scale. This will help later in characterizing our approach and in comparing it to other approaches.

Basics of the effective field theory approach
In SUSY models with very heavy new particles of mass Λ ≫ v, an observable O can be expanded perturbatively in a three-fold way: in terms of loops (counted by a generic loopcounting parameter α), large logarithms of the large mass ratio L ≡ log(Λ/v) and a mass suppression factor v/Λ. For the particular case of a dimensionless observable O which has a tree-level contribution of order α 0 , the leading n-loop contribution (n ≥ 0) is typically of the form α n L n . Subleading/higher-order contributions have more powers of α, fewer powers of L and/or additional factors of v/Λ. Hence, one can write to all orders where the sum of the terms with n = 0 represent the tree-level contribution O 0ℓ and the coefficients c nlk are constants which may contain the parameters of the full model and logarithms of small mass ratios. An effective field theory calculation allows to include all terms at the m-th subleading log level, Since terms of all loop orders are contained in the sum of these terms, their inclusion is also called "resummation of logarithms". Usually the resummation is achieved by performing the following three steps (see figure 1):  1. Construct a Lagrangian of the effective theory and derive a relation between the running parameters of the full and the effective theory by a matching calculation at the m-loop level at the high scale Q match ≈ Λ. 3 2. Use (m + 1)-loop renormalization group running to evolve the parameters of the EFT from the scale Q match to the low-energy scale Q low . In this process the large logarithms are resummed to the N m LL order.
3. Match the parameters of the EFT at the scale Q low to observed quantities and compute the observable in question at m-loop level.
It is not only crucial to take into account all m-loop terms, but it is also important to consistently truncate the perturbation expansion at the m-loop order. In particular, it is imperative not to include any spurious (> m)-loop terms enhanced by large logarithms as this would spoil the correct resummation. On the other hand it is allowed to incorporate m-loop terms suppressed by powers of v/Λ in the m-loop matching, i.e. to take into account fixed order terms of the form α m L k v/Λ. In this way the computation of low-energy observables can be improved by power suppressed terms at fixed loop order. Note, however, that as long as only running of operators of mass dimension ≤ 4 is used, power-suppressed large logarithms of the form L k v/Λ are not resummed [54], i.e. terms of order α n+m L n+(k≤m) v/Λ with n ≥ 1 are not correctly predicted. However, in ref. [44] it was shown that this effect is negligible for the purpose of Higgs pole mass prediction in the relevant parameter space of the MSSM.

Parametrization of the matching relations
In the following we discuss different possibilities to perform the high-scale matching. Specifically, for a matching at some given loop order, one needs to consistently expand either in terms of the running parameter of the fundamental theory α or of the EFTα. In principle both options are correct and equivalent. However, once perturbation theory is truncated it matters whether truncation is done at the order (α) m or (α) m , because these two kinds of expansions differ by higher-order terms. We give a simple illustration using a 1-loop toy example which is similar to the case of the Higgs pole mass calculation. We suppose the exact matching condition is given by the equality where Γ is some Green function. In the full theory, the 1-loop expression is where ∆ γ and ∆ c are numerical coefficients. In the EFT, the 1-loop expression reads The coefficient ∆ γ of the large logarithm L is the same in both cases, because it must cancel in the matching condition. We assume that α andα are related at 1-loop level bŷ The matching condition can now be solved perturbatively forλ in terms of α orα. At tree-level one obtains tree-level :λ = α =α.
(3.7) At the 1-loop level one obtains in terms of α: full-model parametrization 1ℓ :λ = α + α 2 ∆ c (3.8) and in terms ofα: Both expressions (3.8) and (3.9) are valid possibilities for the 1-loop matching relations, but the results forλ differ by non-log-enhanced 2-loop terms. In fact, this difference could be used as an estimate of the theory uncertainty. For the prediction of the Higgs boson pole mass, the EFT parametrization is used in several calculations such as HSSUSY 4 [47,56], MhEFT [43] and SusyHD [42], although further parametrizations have been presented in refs. [23,40]. We note that in an algorithmic implementation of the full-model parametrization, the Green function Γ full may be evaluated numerically, while Γ eft needs to be analytically expanded in terms of α and truncated consistently at the 1-loop level. Hence, an analytic manipulation of Γ eft is needed. Conversely, an algorithmic implementation of the EFT parametrization would require an analytic expansion of Γ full in terms ofα and a consistent truncation of that expansion.
Finally, we note that one might be tempted to plug the respective 1-loop results (3.4)-(3.5) into the matching condition (3.3) to obtain and solve forλ, e.g. numerically. One would then obtain Here, a spurious log-enhanced 2-loop term is generated. If such an implementation were used, the resummation of subleading logarithms would be spoiled. A problem of this kind appeared in refs. [54,55] and a solution was first discussed in ref. [56].

Matching of the quartic Higgs coupling
In the following we will discuss the differences between the two parametrizations in the context of predicting the quartic Higgs couplingλ from a matching of the Standard Model to the MSSM.
EFT (SM) parametrization. In this parametrization the quartic Higgs couplingλ is expressed in terms of the MS-renormalized SM parameters {ĝ 1 ,ĝ 2 ,ĝ 3 ,ŷ t ,ŷ b ,ŷ τ ,v} at the matching scale Q match . In the scenario with degenerate SUSY mass parameters and Q match = M S , the 1-loop contribution toλ from stops is given by where x t = X t /M S is the dimensionless stop-mixing parameter in the DR ′ scheme.
where y t denotes the MSSM top Yukawa coupling in the DR ′ scheme.
With respect to the top Yukawa and strong gauge coupling, the difference between the EFT and the full-model parametrization (3.12) and (3.13) is of 2-loop order. This can be seen by equivalently reparametrizing eq. (3.13) in terms of the SM MS top Yukawa couplinĝ y t , which leads to ∆λ 1ℓ (3.14) Comparing the two versions of the threshold corrections (3.12) and (3.14) reveals several important points. We note first that by construction the 2-loop term on the r.h.s. of eq. (3.14) does not contain large logarithms, in agreement with the effective field theory paradigm. Clearly, the 2-loop difference between eqs. (3.12) and (3.14) could be used as a measure of the theory uncertainty of the 1-loop prediction ofλ at the matching scale. Finally note that this reparametrization generates a 2-loop x 5 t term on the r.h.s. of eq. (3.14). This term is correct, i.e. it appears in the explicit 2-loop calculation of ref. [41]. In section 5 we will show that this is not an accident; the full-model parametrization includes important terms correctly, which in the EFT parametrization would require higher-order calculations. It can thus be used to improve the precision of Higgs pole mass prediction in the effective field theory approach.
Automatization of the matching beyond 1-loop level. Besides the higher precision, the full-model parametrization may also be easier to implement in generic spectrum generators that use the FlexibleEFTHiggs approach [54][55][56]. In this approach the condition is numerically solved forλ at the matching scale. As discussed in ref. [56] and in section 3.2, care has to be taken to avoid the occurrence of spurious large logarithms of higher-order in the matching. A correct application of FlexibleEFTHiggs approach beyond the 1-loop level using the EFT parametrization requires an expansion of the full-model BSM Higgs self-energy Σ BSM ϕ (P ) in terms of the full-model parameters P and a following expansion of P in terms of the parameters of the EFT (here the SM)P , including a truncation at some fixed order inP . This expansion introduces "implicit" terms beyond 1-loop of the form EFT parametrization: where ∆P = P −P is the threshold correction of SM-like parameters expressed through SM parameters. Thus, the inclusion of derivatives of the BSM Higgs self-energy w.r.t. SM-like parameters becomes mandatory for the cancellation of large logarithms in the matching beyond 1-loop. The calculation of these derivatives requires some extra computational effort, which must be performed for each BSM model. The application of this approach to arbitrary BSM models thus requires some cost.
Within the full-model parametrization, the Higgs self-energy in the EFT, Σ SM h (P ), must be expanded in terms of the parameters of the EFT,P , which then must be expanded in terms of the parameters of the full model, P . As a result, 2-loop structures of the following form are generated full-model parametrization: where ∆P =P −P is the threshold correction of SM-like BSM parameters expressed through BSM parameters. Thus, only derivatives of the EFT Higgs self-energy are required. As long as the employed EFT does not change, these derivatives can be computed once and reused in the matching to arbitrary BSM models. With respect to computational effort and model independence, the full-model parametrization is thus advantageous. Due to the re-usability of the appearing structures and the improved treatment of x t (and tan β) discussed later, we propose to use the full-model parametrization instead of the EFT parametrization used in HSSUSY, SusyHD and the original FlexibleEFTHiggs implementation [54,56].

New FlexibleEFTHiggs matching procedure
In the following we apply the conclusions of the previous section to the matching of the SM to the MSSM and describe a new improved matching procedure of the FlexibleEFT-Higgs approach, which also allows to extend the approach beyond the 1-loop level without introducing spurious logarithms of higher order.

FlexibleEFTHiggs matching conditions
The FlexibleEFTHiggs approach is based on the central matching condition where Herev is the minimum of the loop-corrected SM effective potential and Σ SM h , t SM h are the MS-renormalized self-energy and tadpole, respectively. Within the MSSM, the Higgs pole mass is related to MSSM tree-level parameters and DR ′ renormalized loop corrections as Here, the tree-level mass matrix (m 2 ϕ ) ij is parametrized such that the soft-breaking Higgsdoublet mass parameters m 2 H u,d are eliminated by employing the EWSB equations at the loop level. This elimination introduces the tadpoles t ϕ,i on the r.h.s. of eq. (4.5), which are of the same loop order as the momentum-dependent self-energy matrix Σ ϕ,ij (s) of the BSM model. For later convenience we introduce the abbreviations where m 2 h is the SM-like tree-level mass eigenvalue of the matrix (m 2 ϕ ) ij . Combining the previous expressions gives rise to the following relation for the SM quartic couplingλ: This is the master formula for the determination ofλ in the FlexibleEFTHiggs approach; in principle it could be evaluated exactly and at arbitrarily high orders. In particular, it could be evaluated either in the limit v/M S → 0 or by keeping power-suppressed terms of order v/M S . The first option would correspond to the pure EFT approach pursued e.g. in HSSUSY and SusyHD. The second option corresponds to the FlexibleEFTHiggs approach.
For an extensive discussion of this method we refer to refs. [54,56]. As exemplified in appendix A of ref. [54] and in appendix A of the present paper the two options indeed coincide analytically in the limit M S → ∞.
In the following we evaluate the master formula (4.8) according to the following prescription: • We use the FlexibleEFTHiggs hybrid method introduced in ref. [54], i.e. we evaluate eq. (4.8) as it stands, including power-suppressed terms of O(v 2 /M 2 S ) arising in the self-energies and tadpoles.
• Eq. (4.8) is evaluated in the full-model (MSSM) parametrization, which is rather easy to generalize to other SUSY models and allows for a resummation of leading x t and tan β contributions in the Yukawa couplings y t and y b as well as in the quartic couplingλ.
• The threshold correction forλ is calculated at N 3 LO with all 1-loop corrections, 2loop corrections in the gaugeless limit at O(g 2 to the pole mass of the Higgs boson can be found in ref. [19], but we don't include these corrections in the present study. and 3-loop corrections of O(g 4 3 y 4 t ). 6 By including the SM β-functions up to the 4-loop level, this matching allows for a resummation of N 3 LL at the considered order in the couplings. The final Higgs mass prediction will include the complete series of powersuppressed (v 2 /M 2 S ) n terms at 1-loop and 2-loop level at the given orders. However, 3-loop suppressed terms are not included in our calculation, because they are neither publicly available in the literature [37] nor implemented in the Himalaya library [33].
In addition to the master formula (4.8), the matching conditions for the other SM parameters {ĝ 1 ,ĝ 2 ,ĝ 3 ,ŷ t ,ŷ b ,ŷ τ ,v} are given by where M denotes the pole mass of the corresponding particle and Γ is a Green function. The symbols A µ and g µ a denote the QED and QCD gauge fields, respectively. Quarks are denoted as q and SM fermions with a non-vanishing electric charge are referred to as f .

Perturbative expansion of the matching conditions
In this section we perform the explicit perturbative expansion of the master formula (4.8) and the matching conditions (4.9) in the full-model parametrization. As a result we will obtain all building blocks necessary for the 3-loop Higgs pole mass prediction in the improved FlexibleEFTHiggs approach.
We start by performing the matching at tree-level. At this order one has Inserting eqs. (4.10) into the master formula (4.8) and into the matching conditions (4.9) one obtains the MS SM parameters expressed in terms of DR ′ MSSM parameters at tree level:λ (4.11e) For later convenience we denote the tree-level SM MS parameters on the l.h.s. of eqs. (4.11) asP 0ℓ . At the 1-loop level we obtain accordinglŷ where the 1-loop threshold corrections on the r.h.s. of eqs. (4.12) are expressed in terms of MSSM DR ′ parameters. In the pure EFT limit v → 0 the 1-loop threshold corrections can be found for example in refs. [23,40,41]. The explicit calculation of ∆y 1ℓ t and ∆λ 1ℓ beyond the pure EFT limit will be exemplified below. For brevity we denote in the following the 1-loop SM MS parameters on the l.h.s. of eqs. (4.12) generically asP 1ℓ . Similarly, the nloop SM MS parameters are denoted asP nℓ . Furthermore we denote the general threshold correction as ∆P =P −P 0ℓ and specify the notation of a generic n-loop threshold correction as which are expressed in terms of MSSM DR ′ parameters, where ∆P | α n denotes all contributions to ∆P of order O(α n ) . For the prediction of the SM-like Higgs pole mass in the MSSM with the improved FlexibleEFTHiggs approach up to the order O(g 2 3 (y 4 t + y 4 b ) + (y 2 t + y 2 b ) 3 + (y 2 t + y 2 τ ) 3 + g 4 3 y 4 t ), it is sufficient to determine all SM parameters at the 1-loop level, except forλ andŷ t , which must be determined at a higher order. For this reason we describe in the following in more detail the calculation of the threshold corrections ∆y nℓ t and ∆λ nℓ . In order to express these threshold corrections consistently in the full-model parametrization, an extra expansion of the loop corrections in terms of MSSM DR ′ parameters must be performed. We will refer to this procedure as "double loop expansion".
Expansion of the top quark pole mass matching condition. The 2-loop threshold correction for the top Yukawa coupling, ∆y 2ℓ t , can be obtained from the top quark pole mass matching condition eq. (4.9b) with denote the renormalized scalar, left-and right-handed components of the 1-loop top self-energy evaluated at momentum p =m t and p = m t in the SM and MSSM, respectively, without the QCD contributions. The SM self-energies are renormalized in the MS scheme and the MSSM self-energies are renormalized in the DR ′ scheme. In the degenerate SUSY mass limit the 1-and 2-loop QCD contributions are given by [64] ∆m QCD,1ℓ , which yields Expansion of the master formula. The n-loop threshold correction ∆λ nℓ is obtained from the master formula eq. (4.8). To derive the necessary building blocks to express ∆λ nℓ in the full-model parametrization, the following three expansions must be performed: • The prefactor 1/v 2 on the r.h.s. of eq. (4.8) must be expressed in terms of MSSM parameters, which yields 7 • The squared Higgs pole mass in the MSSM, (M MSSM h ) 2 , on the r.h.s. of eq. (4.8) is naturally expanded in terms of MSSM parameters as where s MSSM,nℓ h is obtained from eq. (4.7) with all corrections included up to n-loop level. 8 The 1-loop self-energy and tadpoles can be obtained from SARAH and we expand the self-energy as Note, that the last term on the r.h.s. of eq. (4.22) contributes pure 2-loop Yukawa terms and thus must be taken into account for a consistent Higgs pole mass prediction , together with the corresponding explicit 2-loop self-energy and tadpole contributions. The explicit 2-loop corrections to ∆s MSSM,2ℓ which we take in the gaugeless limit at the order which we take from the Himalaya library [33].
• The pure SM contributions ∆s SM h (p 2 ) from eq. (4.6) must also be expressed in terms of MSSM parameters, which we achieve by a double loop expansion: 7 Threshold contributions to v at 2-loop would induce 2-loop contributions in eq. (4.8) beyond the gaugeless limit. 8 Note that due to the non-linearity of the determinant (4.5), the contribution ∆s MSSM,nℓ h contains products of self-energies of lower loop order, which are, however, suppressed by factors of v 2 /M 2 S as discussed in refs. [57,68]. where (4.26c) The sum on the r.h.s. of eq. (4.26b) runs over P ∈ {p 2 ,ŷ t ,ŷ b ,ŷ τ ,ĝ 3 ,v}, where the ∆p 2 contribution accounts for the fact that the momentum inserted in the 1-loop SM Higgs mass correction of eq. (4.8) is evaluated at Hence the relation includes corrections at 1-loop in the gaugeless limit. The sum on the r.h.s. of eq. (4.26c) 9 Mixed derivatives and products of threshold corrections do not appear in eq. (4.26c), since they would contribute beyond the considered O(y 4 t g 4 3 ) in ∆λ. The explicit 2-loop corrections in the effective potential approach at O(v 2ŷ4 tĝ 2 3 ) are taken from refs. [56,69].
With these ingredients, the expansion of the master formula (4.8) is given bŷ Note that for consistency the parameter shifts ∆P , which contribute to ∆λ 2ℓ and ∆λ 3ℓ , must be evaluated in the gaugeless limit. Note also that the double loop expansion for the 3-loop threshold correction ∆λ 3ℓ is relatively simple, because only a small sub-set of implicit corrections contribute at order O(y 4 t g 4 3 ).

Explicit results and comparison to the literature
In this subsection we present several explicit results and analytic expressions for the threshold corrections at the 1-loop, 2-loop and 3-loop level. The first purpose is to demonstrate the internal consistency of the new way of setting up the threshold corrections by checking that all explicit large logarithms cancel in eqs. (4.18) and (4.28). A second purpose is to verify the correctness of the results by comparing to results presented in the literature, appropriately reparametrized. Finally, we also present several new analytic results at the 2-loop level.
All results in this subsection are provided in the EFT limit v 2 ≪ M 2 S and for the degenerate mass case and non-trivial stop mixing.
We start with the derivation of the 1-and 2-loop threshold corrections for the top Yukawa coupling, ∆y 1ℓ t and ∆y 2ℓ t , at O(y t g 2n 3 ) from eqs. (4.18). The exact 2-loop pole mass contribution in the MSSM (self-energy + momentum iteration) is obtained from refs. [66,67]. The 2-loop pole mass contribution in the SM is taken from ref. [64]. Using the notation as introduced in eq. (4.14) and for the case of degenerate SUSY mass parameters the corrections on the r.h.s. of eqs. (4.18a) and (4.18b) evaluate tô (4.29c) By squaring eq. (4.29a) we obtain perfect analytic agreement with the 2-loop threshold correction presented in ref. [48]. Note that the presented threshold corrections ∆y αs t and ∆y α 2 s t are linear in x t , which will be relevant for the determination ofλ in the full-model parametrization below.
For future use we record here the corresponding result in EFT parametrization in an analogous notation as in eq. (4.29):

(4.30c)
A noteworthy difference is the term ∝ g 4 3 x 2 t , which appears in the EFT parametrized expression in eq. (4.30c), but not in the full-model parametrization in eq. (4.29c). This term originates from an implicit (conversion) correction. Its origin can be traced back to eq. (4.18b), which contains a derivative term of the form when evaluated in the EFT parametrization. The 1-loop correction ∆y αs t and the 1-loop MSSM top quark self-energy Σ MSSM,αs t each contain terms ∝ x t , which results in a contribution ∝ x 2 t in eq. (4.31). As will be discussed below, this 2-loop x 2 t term is also present in the MSSM-parametrized threshold correction, where it is implicitly taken into account by the 1-loop threshold of eq. (4.29b).
∆λ at O(g 2 3 y 4 t ) We continue and show that eq. (4.28c) in that form leads to the known expression of ∆λ 2ℓ at O(g 2 3 y 4 t ), presented in ref. [41]. The calculation of this correction from eq. (4.28c) at 2-loop level requires the explicit 2-loop corrections of O(v 2 g 2 3 y 4 t ) for both the SM [69] and MSSM [10] Higgs pole mass. Since the threshold correction to the VEV ∆v does not contribute at this order, eq. (4.28c) simplifies for Q = M S to ∆λ 2ℓ where we have used ∆y 1ℓ t from eq. (4.29b). In contrast to the result presented in ref. [41], eq. (4.32c) is expressed in the full-model parametrization. To compare our result with eq. (36) of ref. [41], which is presented in the EFT parametrization, one has to express the SM top Yukawa couplingŷ t in ∆λ 1ℓ from eq. (3.12) in terms of the MSSM top Yukawa coupling y t . After truncation at the 2-loop order O(g 2 3 y 4 t ), the combined expression is identical to eq. (4.32c).
∆λ at O(y 6 t ) At 2-loop O(y 6 t ) we can perform an even stricter consistency test of eq. (4.28c), because the implicit corrections at this order have all non-trivial contributions from eqs. (4.22), (4.26b) and (4.28c). Since the explicit 2-loop contributions in the MSSM [68] and SM [42] are known, the only missing contributions in eq. (4.28c) are the implicit corrections. The 1-loop threshold corrections of O(α t ) to P ∈ {p 2 , y t , v} are obtained from ref. [40]. In appendix A we present all contributions and the explicit derivation in an expansion of ct β ≡ 1/ tan β up to the second order for Q = M S . Besides this expansion we have also proven the equivalence between our threshold correction and the one presented in ref. [42] for the general tan β dependence. Since the expression from ref. [42] does (by construction) not contain large logarithmic contributions, we have shown that all large logarithmic contributions vanish in eq. (4.28c) at O(y 6 t ). Our explicit result for the 2-loop threshold correction in the full-model parametrization reads: ∆λ 2ℓ where the analytical result is shown in eq. (A.8).
∆λ at O(y 4 t y 2 τ ) At O(y 4 t y 2 τ ) an interesting comparison of threshold contributions in both parametrizations can be made. In a fixed-order calculation there are no 2-loop diagrams which explicitly contribute to the squared Higgs pole mass M 2 h at O(v 2 y 4 t y 2 τ ). However, in a fixed-order calculation contributions of this order are induced by momentum iteration. In an EFT calculation based on EFT parametrization, the contributions of O(ŷ 4 tŷ 2 τ ) to ∆λ 2ℓ vanish, see the remarks in section 2 of ref. [44]. However, in our matching procedure, which is based on the full-model parametrization, terms of this order are explicitly generated. This can be seen as follows. Evaluating the implicit correction on the r.h.s.
Since the tau lepton and stau slepton do not explicitly contribute to the top quark selfenergy at 1-loop level, the threshold corrections ∆y ατ t and ∆v ατ are related as where the threshold correction ∆v ατ can be obtained from the wave function renormalization of the Higgs, analogously to ∆v αt out of ref. [40]. Eventually, in the degenerate mass case and Q = M S the threshold correction evaluates to ∆λ 2ℓ where ∆s MSSM,y 4 t y 2 τ has been computed similarly to eq. (A.2). We note that large logarithmic contributions within the terms of eq. (4.36a) cancel against each other. In order to validate our result from eq. (4.36b), we have computed the analogous expressions in EFT parametrization by performing a reparametrization of the 1-loop correction in eq. (3.12). We find that the EFT-parametrized threshold correction vanishes, which is in line with the remarks in ref. [44]. The discussion of the correction of O(y 2 t y 4 τ ) is analogous to O(y 4 t y 2 τ ).

Higgs mass loop corrections at
to the CP-even Higgs pole mass are known by ref. [17] and are included in several public codes to calculate the CP-even Higgs pole masses in a fixed-order calculation. In the SM, however, the corresponding 2-loop corrections are not available in a simple and explicit form in the literature to our knowledge. 10 In order to include the 2-loop threshold corrections of O((y 2 t + y 2 b ) 3 ) into the quartic Higgs couplingλ, our approach requires the corresponding fixed-order corrections to be available for both the MSSM and the SM separately. If the contributions in the SM would be omitted, wrong logarithmic enhanced terms of O(v 2ŷ4 ) would propagate into the expression for the Higgs pole mass. Thus, we calculate here the contributions to the 2-loop effective potential of the SM in the gaugeless limit for a non-vanishing bottom Yukawa coupling. The relevant diagram contributing to the mixed contributions of O(v 4ŷ4 Following the approach in ref. [15], we compute the 2-loop bubble diagrams V SM 2ℓ and expand the 1-loop effective potential around the MS-renormalized masses of the top and bottom quark, where (V SM 1ℓ ) ϵ represents the part of the 1-loop effective potential which is proportional to (4 − D)/2 = ϵ and ϕ is a background field. We have checked thatV SM 2ℓ | (ŷ 2 with the subtracted integralsÎ andĴ instead of I and J, which have been introduced in ref. [71]. After expanding the 2-loop integrals around the renormalized Goldstone mass parameter and taking only Yukawa coupling enhanced contributions into account, the finite result expressed in SM MS parameters readŝ The squared top and bottom quark mass parameters T =ŷ 2 t ϕ 2 /2 and B =ŷ 2 b ϕ 2 /2 are expressed in terms of the background field ϕ. In the MSSM, the corresponding SM-like contributions are included in eqs. (3.39) and (3.40) of ref. [13] and we have checked that the effective potential in eq. (4.38) reproduces these results when omitting the corrections from BSM Higgs bosons. The loop corrections to the Higgs pole mass are derived from the effective potential by differentiating w.r.t. to the background field ϕ as Because the implicit corrections of eq. (4.26b) are analogous to eq. (4.34b), the shift ∆s SM,2ℓ is completed by derivatives of eq. (4.38). Neglecting terms of order O(ŷ 4 bb 2 /t), the loop corrections read . We note that the first line in eq. (4.40) corresponds to the loop corrections of O(ŷ 6 tv 2 ), which can be found in ref. [42], and which we reproduce here. In the same way as prescribed in eq. (2.49) of ref. [70], we checked the renormalization scale invariance of the Higgs pole mass at withτ ≡m 2 τ and log(τ ) ≡ log(τ /Q 2 ).
In order to check the consistency of the 3-loop expression of eq. (4.28d), we derive the second term on the r.h.s., ∆s SM,3ℓ in the full-model parametrization and compare it to the result presented in ref. [48]. At 3-loop O(v 2 g 4 3 y 4 t ), the SM Higgs mass correction of eq. (4.26c) receives an explicit self-energy and tadpole contribution, which we take from ref. [70]. We determine the implicit (derivative) contributions of eq. (4.26c) using the 1-loop threshold correction ∆g αs 3 from ref. [72] and the 2-loop threshold correction ∆y α 2 s t form eq. (4.29c). 11 In the EFT limit and for degenerate SUSY 11 Note that in this section we ignore contributions of O(v 2 /M 2 S ) for brevity and for cross-checking against expressions from the literature. In our actual implementation of the corrections, presented later, we take all available terms of O(v 2 /M 2 S ) into account.

Resummation of leading squark mixing contributions
In the introductory example of section 3.2, differences between the full-model and EFT parametrization were discussed. It was shown that both approaches are equivalent up to higher-order terms, which contribute numerically to the difference of both approaches. From a technical point of view, implementing the full-model parametrized matching is easier to achieve; in this section we present a second, more important argument in favor of this approach: the resummation of higher-order contributions of the full-model squark mixing parameter x f ≡ X f /M S in the context of Higgs mass predictions.
The resummation is analogous to the resummation of large n-loop (tan β) n -corrections to m b of refs. [73][74][75][76], suitably generalized. 12 We begin here by recalling main features of the tan β-resummation in m b , rephrase it in the appropriate language and then present the generalization. More details and further generalizations will be presented elsewhere [77]. Here ∆m b is the loop correction between the b-quark pole mass and DR ′ running mass, The quantity ∆m b contains a 1-loop term of the order α s tan β, terms 12 At this point we want to note that the value of xt is bounded by the necessity of avoiding charge and color breaking minima [41]. For large mq 3 ≈ mũ 3 the absolute value of the dimensionless stop-mixing parameter |xt| is restricted to be less than 3, whereas tan β can be as large as 50-60. with lower orders in tan β and terms governed by other couplings, but no higher-loop terms of the orders given in the theorem. The theorem only holds for "unsuppressed" terms, i.e. terms not suppressed by powers of v/M S . The theorem can be equivalently formulated in the language of full-model versus EFT parametrization of the b-quark mass matching between the MSSM and the SM. Full-model parametrization means to express the SM MS bottom quark massm b as a perturbative series in the full-model MSSM parameters m b , α s , . . . , and truncating the series at some desired order. The full-model parametrized relation betweenm b and m b , truncated at order α s then reads 13m (5.1) Here the notation of eqs. (4.14) and (4.29a) has been used, and the subscript m b α ≤1 s refers to the full-model parametrization and the chosen truncation order. The dots denote terms irrelevant for the present discussion (containing less powers of tan β and/or power suppressions). In a numerical code, this equation is often numerically solved for m b , effectively giving (5. 2) The point of the theorem is that eq. (5.1) is 1-loop exact with respect to (α s tan β) n -terms and eq. (5.2) correctly takes into account ("resums") all terms of these orders. 14 On the other hand, in a calculation using the EFT parametrization, m b would be expressed as a perturbative series in terms ofm b ,α s , . . . , 15 truncated at some desired order. E.g. at orderα s , which contains the correct (α s tan β) n term only for n = 1 but misses all higher-order terms. Accordingly, we can evaluate the difference y f -matching and x f -resummation: An analogous resummation is possible for the α s x f -enhanced contributions to the Yukawa matching for all colored fermions. We note that the factor µ tan β in the previous theorem arises via X b = (A b − µ * tan β). Hence the above theorem generalizes to 13 As implicitly indicated by ref. [75], a stronger restriction can be formulated which forbids unsuppressed terms of O(α n s tan >1 β) in the threshold corrections to the bottom mass matching between the THDM and MSSM. 14 We note that all codes mentioned in the present paper resum the m b (αs tan β) n corrections in this way, even if they otherwise do not use full-model parametrization. 15 For the purpose of the present section the distinction between the SM parameterαs and the MSSM parameter αs is not relevant, since the two parameters differ by terms which do not depend on tan β or x f .
There are no unsuppressed contributions to the threshold correction ∆y f at O(α n s x >1 f ) for n ≥ 1 in full-model parametrization.
In the following we specialize to the case of the top quark. If we express the SM couplingŷ t as a perturbative series in y t , α s , . . . and truncate at the order α s , we obtain the full-model parametrized expressionŷ The theorem then implies that this relation is 1-loop exact, i.e. there are no unsuppressed higher-order corrections to ∆y αs In a numerical code such as FlexibleEFTHiggs, eq. (5.5) is solved numerically for y t , giving where the chiral symmetry ensures that for QCD corrections ∆y αs t ∝ y t , such that the denominator has the structure 1 + O(α s x t ) + · · · . This expression correctly "resums" in particular all terms of the orders O((α s x t ) n ). In the EFT parametrization, however, y t is expressed as a perturbative series inŷ t ,α s , . . . , truncated at some desired order. 16 The resulting expression is then correct up to that order but misses all higher-order terms of the form (α s x t ) n . Similar to eq. (5.4), we can express the difference between the two parametrizations as Here the subscript on the r.h.s. denotes the truncation of the EFT parametrization at n-loop level, and the dots denote terms irrelevant for the present discussion. The main point is that the difference contains terms of the orders (α s x t ) k with k > n, and all these terms are correctly contained in the full-model parametrization but missing in the EFT parametrization, and the coefficients a k can be analytically computed at all orders. As an example of this discussion, we refer to the analysis on the y t matching given in section 4.3: In (4.30c) the explicit 2-loop threshold correction expressed in SM parameters contains the term ∝ŷ t (α s x t ) 2 . According to eq. (5.7), the MSSM Yukawa coupling in eq. (4.30a) misses terms ∝ŷ t (α s x t ) ≥3 which, however, are implicitly taken into account in the 1-loop correction of eq. (4.29a).
λ-matching and x f -resummation: Now we turn to the resummation of O(g 2 3 x t )contributions in the matching ofλ in the full-model parametrization. In section 3.3 we already presented an example where the 1-loop matching ofλ and y t in full-model parametrization leads to a correct 2-loop term leading in x t . This example illustrates a more general property; however the resummation withinλ is slightly more complicated than the two cases discussed above. We first recall that the threshold correction forλ in full-model parametrization, truncated at 1-loop order, contains the terms where the c 0x are coefficients and the dots denote terms irrelevant for the present discussion. The terms leading in x t are thus of the general form y 2 t g 2 x x m t with g x ∈ {y t , g 1 , g 2 } and m ∈ {4, 2, 2}. The resummation of higher-order terms governed by (g 2 3 x t ) n relies on the resummation within the Yukawa coupling discussed before and on the following theorem for the explicit contributions to the threshold corrections: There are no unsuppressed contributions to ∆λ at O(y 2 t g 2 x g 2n 3 x >m t ) for n > 0 in full-model parametrization.
Again the theorem means that the full-model expression (5.8) is 1-loop exact with respect to the leading O(g 2 3 x t ) terms. The proof of the theorem is analogous to the proof in ref. [75] and relies on the large mass expansion and inspection of individual contributions. For details and generalizations we refer to ref. [77]. If the Yukawa coupling in eq. (5.8) is replaced bŷ y t via eq. (5.6), we obtain equations of the form with appropriately modified coefficientsĉ 0x from tree-level matching and higher-order coefficientsĉ nx . In a numerical code based on the full-model parametrization these higher-order terms are fully taken into account; the coefficients are also fully calculable analytically. All the explicit higher-order terms in eq. (5.9) are correct. In this sense the full-model parametrization resums these terms of O(ŷ 2 . On the other hand, if theλ-matching is done in EFT parametrization, ∆λ is expanded in terms of SM parameters and truncated at some desired order. In that case, leading O(ĝ 2 3 x t ) terms are only taken into account up to that order. The difference between the two versions of the threshold corrections can be written as All the explicitly given terms ∝ĉ kx on the r.h.s. of eq. (5.10) are taken into account correctly in the full-model parametrization but are missing in the EFT parametrization.
In other words all these higher-order terms in eq. (5.10) do not arise from explicit multiloop diagrams; instead they only arise via the reparametrization, and the values of the coefficientsĉ kx can be deduced from the 1-loop terms leading in x t . As an example, we reconsider the concrete calculation in eq. 3 ). The 1-loop threshold correction ∆λ 1ℓ at O(y 2 t g 2 1,2 ) originates from D-term contributions and involves a maximum x t -dependence of order x 2 t . Thus, our prediction for the leading x 3 t contribution at 2-loop mixed electroweak order is (4π) 6 All these terms and the corresponding terms of ≥ 4-loop order are automatically taken into account by the new FlexibleEFTHiggs calculation based on full-model parametrization.
6 Running and matching procedure at the electroweak scale In this section we describe the computations which are performed after the SM parameters are obtained at the high scale by the matching characterized in section 4. As illustrated in figure 1, the subsequent calculation involves two further steps. First the SM parameters are run down to the low-energy electroweak scale by solving the RGEs (subsection 6.1). Second, the low-scale parameters are related to input values of observables, and the final prediction for the Higgs boson mass is computed. The corresponding low-scale matching procedure is described in subsection 6.2.  [86,87], see also ref. [56].

Matching of SM couplings to observables
In order to express the prediction for the Higgs pole mass in terms of physical quantities, the running SM MS parameters have to be related to observables. There are eight SM MS parameters relevant for the Higgs mass prediction (cf. eq. (2.1)): Among these eight parameters,λ is fixed by the matching to the MSSM at the SUSY scale, while the other seven parameters are fixed by low-energy observables. Following the approach described in ref. [56], we fix these seven parameters at the scale Q = M Z by relating them to the seven low-energy quantities shown in table 1.  Matching procedure. Before continuing the discussion about the included loop corrections at the electroweak scale, we want to emphasize a qualitative difference of our low-scale matching procedure w.r.t. the procedure described in ref. [56]. In our new FlexibleEFT-Higgs matching approach we consider the full-model (MSSM) DR ′ parameters at the SUSY scale to be fundamental. This includes the SUSY parameters and the SM-like full-model parameters P from eq. (2.7). As a consequence, we eventually express all observables, including the low-energy observables from table 1 as well as the predicted Higgs pole mass in terms of the full-model DR ′ parameters. Technically, this is achieved by (i) converting the

Quantity Description
In contrast, in spectrum generators working in the EFT parametrization (e.g. HSSUSY, SusyHD or MhEFT), the low energy EFT MS parametersP (Q low ) (except forλ) are not determined from the full-model parameters, but rather they are directly extracted from the observables asP with some function h i denoting the calculation. The difference between the EFT parameterŝ P from both approaches (6.2) and (6.4) is of higher order. The different higher-order terms depend dominantly on the values ofP (Q low ). However, eqs. (6.2) and (6.4) can be modified by higher orders in such a way that the differences expressed in terms ofP (Q low ) vanishes, e.g. see the discussion below eq. (6.5). Even if the conditions (6.2) and (6.4) coincide in their inclusion of loop corrections, a remaining difference in the numerical value ofP (Q low ) may occur depending on the parametrization of the high-scale matching condition forλ. Though the difference, which originates from reparametrization effects at the high scale, depends only subdominantly on the SUSY parameters, e.g. tan β. 17 Matching conditions. As described in section 4, we aim for a prediction of the Higgs pole mass at N 3 LO with N 3 LL resummation in QCD. This precision requires to relate the SM MS top quark massm t to the pole mass M t up to O(α 3 s ) [88,89]. We follow the prescription presented in ref. [56] and express M t as (6.5) 17 Considering the difference of the threshold correction ∆λ between the EFT and the full-model parametrization, the reparametrization effects do explicitly depend on SUSY parameters. Inserted into the β functions of the SM, the reparametrization terms receive further loop suppressions. Therefore, the higherorder reparametrization effects which contain MSSM-specific parameters propagate into all EFT parameters at the electroweak scale with additional loop factors.
Note, that we have included the 4-loop QCD contribution from ref. [90] for later use in section 8. In eq. (6.5) M input t denotes the observed top quark pole mass, Q is the renormalization scale and Σ SM,1ℓ t,{L,R,S} (p 2 , Q) are the left-handed, right-handed and scalar parts of the 1-loop SM top quark self-energy without QCD contributions. The separate SM QCD contributions read [88][89][90] ∆m QCD,1ℓ t (Q) ≈ĝ witht =m 2 t and log(x) ≡ log(x/Q 2 ). As referred to before, the construction of eq. (6.5) contains the observed value of the pole mass M input t , which introduces higher-order terms in eq. (6.2) such that the higher-order difference to the low-scale constraint used in HSSUSY is minimized, see eq. (7) from ref. [56]. The Fermi constant G F is calculated similarly to ref. [72] as , where the 1-loop correction ∆α 1ℓ em (Q) is parametrized in terms of the 6 flavorα em (Q) as The loop corrections of ∆α nℓ s (Q) are defined in eqs. (13)-(15) of ref. [56]. The 1-loop MSrenormalized self-energies Σ 1ℓ τ (p, Q) and Σ 1ℓ Z (p, Q) of the τ lepton and Z boson, respectively, are evaluated at full 1-loop precision. In eq. (6.8e) Σ 1ℓ,heavy b (p, Q) denotes the 1-loop top quark and electroweak gauge boson contributions to the bottom quark self-energy as described in ref. [62]. The predictedm to be compared to the input valuem ). In principle, one might treat the Higgs mass similar to all other low-scale observables and use it as a constraint in the sense of eq. (6.3) to fix one parameter of the MSSM. However, in our application we choose the Higgs mass at the low scale to be the output of our calculation. Analogous to eq. (1) of ref. [56] we compute the Higgs pole mass in the SM as in eq. (4.3) but instead of using MSSM parameters we express it in terms of SM parameters at Q = M input t . The included corrections are described in section 4 and are of

Numerical results
In this section we present the numerical results for the light CP-even Higgs pole mass calculation in the real MSSM, based on the improved FlexibleEFTHiggs approach. We highlight the numerical impact of the new included threshold corrections and the parametrization scheme and describe in particular the x t -resummation.
If not stated otherwise, the dimensionful DR ′ -renormalized parameters of the MSSM Lagrangian are set to a common SUSY scale M S , For scenarios with non-trivial squark mixing we parametrize our results in terms of the dimensionless DR ′ parameter x t ≡ X t /M S . In our numerical discussion we choose the input values given in table 2. Effects from 1 st and 2 nd generation (s)fermions are omitted in our analysis.

Impact of higher orders, the new parametrization and the x t -resummation
We begin with the discussion of the impact of higher-order corrections in the matching and the impact of the new full-model parametrization and the resulting x t -resummation. Figure 3 shows the light CP-even Higgs pole mass calculated by different versions of Flex-ibleEFTHiggs: • FEFT 1ℓ/2ℓ/3ℓ: The new FlexibleEFTHiggs hybrid calculation developed in the present paper with full-model (MSSM) parametrization of the matching calculation and x t -resummation. The only difference among these calculations stems from the orders taken into account in the calculation ofλ, see eq.   We first discuss the two 1-loop versions (blue dotted and green dashed-dotted lines). They differ essentially by the full-model versus EFT parametrization of the high-scale matching. As discussed in section 5, in the full-model parametrization certain leading x t terms are correctly taken into account. Eq. (5.9) provides the general form of those terms. To exemplify the correctly included terms, consider the 1-loop threshold correction ∆λ 1ℓ in the MSSM parametrization. It contains terms of the order O(y 4 t x 4 t ) and mixed electroweak terms of the form O(y 2 t g 2 1,2 x 2 t ). Upon reparametrization of these contributions in terms of SM parameters, using in particular the 1-loop top Yukawa threshold correction of O(y t g 2 3 x t ), the terms shown in table 3 are generated. As discussed in section 5, these terms are not modified by genuine n-loop contributions. Hence, already the new FEFT 1ℓ calculation correctly takes into account the leading QCD (n + 1)-loop contributions of the order x In contrast, none of the terms in table 3 is correctly taken into account in the previous FEFT 1ℓ (SM para.) calculation. As a result of the x t -resummation, we see a dramatic shift between the two 1ℓ versions in figure 3.
• Without x t -resummation, the FEFT 1ℓ (SM para.) result (green dashed-dotted line) and the EFT 3ℓ calculation (green solid line) deviate up to ∼ 1.1 GeV for |x t | < 3.  . Note, that the FEFT 1ℓ calculation is based on the full-model parametrization, but the terms in this table are provided in terms of SM parameters, to compare with other SMparametrized calculations. Note further, that the terms in this table are already contained in the 1-loop calculation; further, higher-order terms resummed by FEFT 2ℓ/3ℓ are determined by eq. (5.9).
• With x t -resummation, this deviation between the FEFT 1ℓ and FEFT 3ℓ calculation decreases to less than 0.3 GeV, compare the blue dotted and the red solid lines.
For all values of x t and M S , the new FEFT 1ℓ calculation is far closer to the FEFT 3ℓ calculation; hence the convergence of perturbation theory is significantly improved. As a side remark we want to stress that in contrast to the correctly included contributions, the 1ℓ calculations in both parametrizations fail to capture the pure Yukawa   highest power in x t , respectively, which are implicitly contained by the 1-loop calculation in full-model parametrization (second column) and in EFT parametrization (third column). The terms in these columns have been obtained by inverting the 1-loop relation between y t andŷ t in full-model parametrization perturbatively up to the 3-loop level. The terms are compared to the known/unknown correct result in EFT parametrization in the last column. The contributions in the EFT parametrization are vanishing by construction. At 2-loop we see that the implicitly included term from the 1-loop correction in full-model parametrization lies in between the correct result and the analogous one in the EFT parametrization. Without further information from an explicit 3-loop (or higher) calculation there is no indication that the full-model parametrization worsens the convergence of the perturbative expansion with respect to these orders. The 3-loop term is studied numerically together with the reparametrization terms from the 2-loop correction ∆λ 2ℓ y 6 t in figure 7  and figure 8.
Second, we discuss the impact of the 2-loop and 3-loop threshold corrections in the new FlexibleEFTHiggs calculation. As can be seen in figure 3, the impact of the higherorder corrections is very small and below 0.3 GeV for all values of M S and |x t | < 3 in the shown scenarios. The main reason is again the x t -resummation: When going from 1ℓ to 2ℓ, the 1ℓ calculation already contains the leading 2-loop QCD x 5 t term of table 3, and the actual 2ℓ calculation only adds subleading x ≤4 t terms. Similarly, one can show that the 2ℓ calculation already correctly contains the leading 3-loop QCD x 5 t and x 6 t terms, and the actual 3ℓ calculation only adds subleading x ≤4 t terms. 18

Comparison to state-of-the-art calculations
In this subsection we compare our new improved FlexibleEFTHiggs calculation with the two state-of-the-art 3-loop fixed-order and EFT calculations from refs. [33,48]. Both of these calculations are also based on the FlexibleSUSY framework [56,62], which facilitates the comparison. In detail, these calculations are: • FO 3ℓ: This is the fixed-order calculation, which has been presented in ref. [33] (dashed-double-dotted magenta line in figure 4). It includes loop corrections to the Higgs pole mass in the full-model (MSSM) parametrization at full 1-loop level and 2and 3-loop corrections in the gaugeless limit at O(v 2 (g 2 3 (y 4 t + y 4 b ) + (y 2 t + y 2 b ) 3 + (y 2 t + y 2 τ ) 3 )) and O(v 2 (g 4 3 y 4 t )), respectively.
• EFT 3ℓ: This calculation is the pure EFT calculation from ref. [48], where a matching at the SUSY scale is performed in the EFT (SM) parametrization (dashed-dotted green line in figure 4). The threshold correction ∆λ includes the known 1-loop contributions from ref.      The right column of figure 4 shows the differences between the calculations in more detail. In the following we discuss these differences. We first focus on the differences between the new FlexibleEFTHiggs and the pure EFT calculation at SUSY scales above a few TeV. τ ) have an impact of ∆M h ≈ 10 −3 MeV and are thus negligible. The essential difference between the FEFT 3ℓ calculation (red line) and EFT 3ℓ (green dashed-dotted line) is the different parametrization of the matching in terms of either MSSM or SM parameters, and the resulting leading higher-order x t terms included in FEFT 3ℓ. The precise origin of this difference in the threshold correction ∆λ can be inferred from table 3. The FEFT 3ℓ calculation correctly includes all terms of the table, while the EFT 3ℓ calculation only includes the terms of the left column at most up to the 3-loop level, but neither includes the 4-loop term nor any term of the right column. As will be shown in the following subsection, the numerically dominant effect comes from the mixed QCD-EW 2-loop terms of the form (λ FEFT 3ℓ − λ EFT 3ℓ ) ⊃ŷ 4 tĝ 2 1,2ĝ 2 3 x 3 t . As shown in the middle row of figure 4, the numerical difference originating mainly from these terms remains around 200 MeV for M S = 100 TeV (and x t = − √ 6). 19 On the other hand, the lowest row of figure 4 shows that for fixed M S the numerical difference is below 200 MeV for |x t | ≤ 2 and M S = 3 TeV, but for larger |x t | the difference rises strongly.
Secondly, we focus on the comparison between the fixed-order and the FlexibleEFT-Higgs calculations for SUSY scales below around 1 TeV, where both calculations should be valid. By construction, both calculations include the same Higgs pole mass contributions of the orders O(1ℓ + v 2 (g 2 3 (y 4 t + y 4 b ) + (y 2 t + y 2 b ) 3 + (y 2 t + y 2 τ ) 3 ) + v 2 (g 4 3 y 4 t )), including terms suppressed by v 2 /M 2 S . However, they differ at other orders. Numerically, the difference is below 0.5 GeV for small x t and M S ≲ 500 GeV (see top row of figure 4), but the difference reaches around 1 GeV for large |x t | and small M S (see middle row of figure 4). The origins of these differences are the following: • Parametrization: In contrast to our hybrid approach, the determination of the DR ′ MSSM top quark mass m t in the fixed-order calculation consists of the following expanded version of the exact relation where ∆m 1ℓ,2ℓ t represent the 1-and 2-loop corrections to the DR ′ top quark mass as described in refs. [47,56]. Analogously to section 5, eq. (7.2) does not represent an all order resummation of terms in the top mass parameter of m t ⊃m t × (ĝ 2 3 x t ) n . Consequently, eq. (7.2) does not lead to an all order resummation of terms in the Higgs pole mass of the form for n > 2. Besides these non-resummed terms, our new FlexibleEFTHiggs hybrid

Further details on the comparison of hybrid and pure EFT calculations
In the lower-right panel of figure 4 one can see a deviation between the hybrid FEFT 3ℓ calculation and EFT 3ℓ for large |x t |. In the following we elaborate on the large-x t behavior in more detail. For the discussion it is sufficient to consider the 2-loop calculations. Figure 5 shows the Higgs pole mass of different 2-loop calculations w.r.t. the FEFT 2ℓ calculation (red solid line). The black dashed line corresponds to the same 2-loop calculation, where 2-loop threshold corrections proportional to powers of y b and/or y τ have been omitted. One finds that the difference between these lines is smaller than 50 MeV for the shown parameter scenario. The blue dotted line represents a modified calculation of the black dashed line, where the 2-loop threshold correction toλ has been replaced by the analytic 2-loop expressions from eqs. (4.32c) and (4.33) of the order O(g 2 3 y 4 t + y 6 t ), where terms of O(v 2 /M 2 S ) have been neglected. Thus, the difference between the blue dotted and the black dashed lines corresponds to the impact of some 2-loop higher-dimensional operators. The effect of these higher dimensional operators has been discussed in ref. [44], where it has been shown that they are of high relevance for large stop mixing. For small |x t | ≲ 3 and the shown value M S = 3 TeV, however, their effect is negligible. Note, that in figure 5 the hybrid 2-loop result is subtracted from each calculation. Hence, the blue dotted line represents the negative correction due to power suppressed terms. In contrast, figure 4 of ref. [44] shows the positive influence of higher dimensional operators. From the figure we draw the following conclusions: • The excellent agreement between the black dashed and the blue dotted lines for |x t | ≲ 3 confirms numerically the correctness of our automatized FlexibleEFTHiggs polemass matching procedure forλ at O(g 2 3 y 4 t + y 6 t ).
• For |x t | ≳ 3 the effect of the higher-dimensional 2-loop operators is in line with the numerical results of ref. [44].
For reference we also show in figure 5 the EFT 2ℓ calculation, represented by the green dashed-dotted line. One finds that EFT 2ℓ deviates numerically from FEFT 2ℓ for |x t | ≳ 1. This discrepancy can be explained by contributions originating from the different parametrization schemes. As motivated above, we categorize the higher-order corrections in two classes of terms; the ones which are incomplete in both approaches and the ones which are captured correctly in our full-model parametrization scheme, but not in the other EFT parametrization.
Concerning the higher-order terms correctly captured by our new FlexibleEFTHiggs hybrid calculation, we find the most dominant contribution to the numerical difference between the EFT 2ℓ prediction and FEFT 2ℓ to be the 2-loop mixed QCD-EW term from 3 ) are taken into account. One finds that this reparametrized calculation agrees well with the corresponding EFT 2ℓ calculation (green dashed-dotted line), which uses the same parametrization. The only difference between the blue dotted and the green dasheddotted line are power suppressed contributions in the Higgs mass at 1-loop, which become significant for |x t | ≳ 3, as discussed above. When adding the 2-loop leading x t mixed QCD-EW contribution from eq. (5.11) to the black dashed line, one obtains the red solid line. The so obtained result agrees very well with the MSSM-parametrized calculation (blue dotted line), which explains the dominant part of the deviation between the MSSMparametrized FEFT 2ℓ calculation and the EFT-parametrized EFT 2ℓ calculation. Thus, the numerical effect coming from the correct inclusion of highest power x t contributions in our new FlexibleEFTHiggs approach improves the precision for large |x t | in comparison to the calculation performed in the EFT parametrization.
Besides the higher-order terms correctly taken into account by our new FlexibleEFT-Higgs calculation, the threshold corrections ∆λ differ in both approaches by further terms, which are incomplete both in the full-model parametrization and in the EFT parametrization. Such incomplete higher-order terms are for example top Yukawa enhanced 3-loop  Figure 7: Impact of incomplete higher-order contributions to ∆λ from reparametrization in FlexibleEFTHiggs.
terms with high x t powers of the form The reparametrization of the 1-loop correction alone was discussed in table 4. The discussion here is extended by the gauge-less 2-loop contributions to ∆λ in MSSM parametrization. In figure 7 we show the numerical influence of such terms. When these (incomplete) higherorder terms are added coherently (green solid line), both contributions almost cancel up to a remaining effect of ∼ 150 MeV in the Higgs pole mass for |x t | < 3.5 and M S = 3 TeV. Thus, the numerical effect from the x t -resummation terms in ∆λ QCD-EW from figure 6 remains the dominant reparametrization effect. However, when the numerical effect of each incomplete higher-order term is drawn individually, the contributions have a higher impact on the Higgs pole mass, see figure 7. The magenta dashed-triple-dotted line corresponds to the effect of the terms of O(ŷ 6 tĝ 2 3 x 7 t ) and the blue dashed line corresponds to O(ŷ 8 t x ≤8 t +ŷ 6 tĝ 2 3 x <7 t ). There is a cancellation between these incomplete contributions, which should be kept in mind when using such terms as an uncertainty estimate of missing higher-order corrections. Using the maximum effect of all terms provides a more conservative estimate of the remaining uncertainty than the coherent sum.

Uncertainty estimation
In this section we analyze missing higher-order contributions in our new FlexibleEFTHiggs approach in order to estimate the remaining theory uncertainty of our calculation. In accordance with refs. [41,42,47] we distinguish between missing higher-order contributions 20 Note that in order to investigate the complete reparametrization contributions of this order, the inclusion of 2-loop threshold corrections to ∆yt at O(y 5 t + g 2 3 y 3 t ) is required.
in the matching at the SUSY scale, which we denote as high-scale uncertainty, and missing loop corrections at the electroweak scale, denoted as low-scale uncertainty. Note, that since FlexibleEFTHiggs is a hybrid calculation, we do not assign an EFT uncertainty to our calculation from missing terms of O(v 2 /M 2 S ). 21

High-scale uncertainty
We begin our discussion by presenting our methods to estimate the high-scale uncertainty, i.e. the numerical impact of the missing higher-order corrections in the matching of the MSSM DR ′ to the SM MS parameters at the SUSY scale. We discuss three different approaches: the variation of the matching scale, implicit higher-order corrections from the double loop expansion and reparametrization terms.
Variation of the matching scale. A commonly applied strategy to estimate higherorder contributions is to vary the renormalization scale Q match at which the threshold corrections are computed. For reasons of comparability, we use the conventional range of Q match ∈ [M S /2, 2M S ] and take the maximum deviation from the value obtained at Q match = M S as an estimate The numerical variation of M h results from the fact that the matching corrections contain explicit dependencies of log Q 2 at fixed order, while the RGE running cancels those logarithms but also generates log Q 2 terms at higher orders. The quantity ∆M Q match h thus represents an estimate for these missing logarithmic higher-order terms. In particular, in the matching ofλ, the following 2-loop and 3-loop terms are generated: ∆λ Q match h ⊃ ∝ y 2 t g 2 3 g 2 1,2 + y 2 t g 2 1,2 + g 4 1,2 + y 2 t g 4 3 x 4 t + y 4 t g 2 3 + y 6 t + O(g 6 1,2 ).
The matching-scale variation thus provides an estimate of the theory uncertainty related to these terms, at least to their log Q 2 -dependent parts. We have omitted the specification of the powers of x t in most of the terms. The term of O(y 4 t g 4 3 x 4 t ) deserves special attention: In the degenerate mass case the Himalaya library up to version 3.0.1 does not provide the correct term in the Higgs mass correction at this order [33]. Since this is an important missing term of higher order in x t , but not of higher order in the couplings, we have verified that this missing term of this order has a non-vanishing log Q 2 dependence. Indeed, employing 2-loop β functions from ref. [96] on the 2-loop Higgs pole mass, derived from the effective potential of ref. [97], the renormalization scale dependence in the degenerate mass case is given by Our calculation of the Higgs mass does not include suppressed logarithms beyond the 2-loop gaugeless limit. In fact, in ref. [44] it has been demonstrated that their impact is very small ∆M h ≤ 20 MeV for the studied scenarios.
Thus, the matching scale variation in our calculation provides an estimate of the uncertainty originating from missing logarithmic terms at O(y 4 t g 4 3 x 4 t ) in particular. We'd like to point out a technical difficulty in this matching scale variation. The evolution of RGEs in the MSSM requires the numerical input values of MSSM DR ′ parameters as a boundary condition. However, in the MSSM two parameters cannot be fixed by the input; rather they have to be eliminated by imposing the two electroweak symmetry breaking conditions. Solving these so-called tadpole equations at the loop level will introduce logarithms which contain light masses. Hence, it is a legitimate question to ask whether such contributions spoil the automatized cancellation of large logarithms in the matching correction. In our calculation, the tadpole equations at the SUSY scale are solved for the dimensionful soft-breaking Higgs-doublet mass parameters m 2 Hu and m 2 H d . An explicit calculation up to leading 2-loop QCD order shows that large logarithms enter intoλ with a suppression of v 2 /M 2 S beyond the considered order. These contributions would be absent in a pure EFT calculation and they can be regarded as a power-suppressed contribution in a hybrid calculation.
Implicit corrections at higher order. In section 4 we discussed the expansion of the master formula (4.8) and explained how "explicit" contributions from genuine multi-loop diagrams are accompanied by "implicit" corrections in the double loop expansion, i.e. from the reparametrization of the SM self-energy in terms of MSSM parameters. These implicit corrections have the form of products of derivatives of the SM Higgs pole mass shift ∆s SM h times parameter shifts.
Hence, as another estimate of missing higher-order corrections, we compute further terms with such a structure at orders beyond the precision of the included threshold corrections and discard terms which contain logarithms of the form log(m t /Q). The resulting contributions take the form where ∆λ 2ℓ g 1 ,g 2 denotes terms which would arise in an actual 2-loop calculation beyond the gaugeless limit, and ∆λ 3ℓ g 3 ,yt contains terms which would arise in an actual 3-loop calculation in the gaugeless limit. The corresponding orders in couplings are The 3-loop gaugeless contributions contained in the generated terms on the r.h.s. of eq. (8.6) are of the order as indicated in the subscript. We can thus first define an estimate of the size of the missing 2-loop electroweak SUSY corrections as where M h (∆λ 3ℓ ) denotes the FEFT 3ℓ calculation. Next, we can define an estimate of the size of missing higher-order SUSY-QCD contributions as where we subtract the reparametrization terms from the FEFT 3ℓ calculation in order to reproduce the truncation of the EFT parametrization ofλ at O(1ℓ +ĝ 2 Up to a sign, the dashed-triple-dotted magenta line and the dashed blue line in figure 7 shows equivalently the numerical influence of the terms estimated by ∆M rep,g 3 ,3ℓ h and ∆M rep,g 3 yt,3ℓ h . At this point it is worthwhile to discuss the difference in the estimation of the uncertainty of an EFT-parametrized calculation at similar order, i.e. with a matching ofλ The reparametrization provides a way to estimate higher-order terms in this calculation, which are sensitive to high powers of x t . Furthermore, the uncertainty estimation should also cover terms of O(ŷ 8 t ,ŷ 6 tĝ 2 3 ), which are incomplete in both parametrizations, c.f. table 4. Consequently, if the discussed techniques are applied to construct higher-order terms for the uncertainty estimation of the EFT-parametrized calculation, we expect that they lead to very similar expressions for ∆λ ⊃ŷ 8 t x ≤8 t +ŷ 6 tĝ 2 3 x ≤7 t . Note, that in contrast to the full-model-parametrized calculation, the EFT-parametrized one would in addition have to estimate the size of the terms of ∆λ ⊃ŷ 2 which are implicitly captured in full-model parametrization. Thus, in EFT parametrization more higher order contributions would be needed to estimate the uncertainty for large |x t |.

Low-scale uncertainty
In this section we describe our method to estimate the low-scale uncertainty, i.e. the theory uncertainty from missing higher-order loop corrections in the matching to the SM input parameters at the electroweak scale. We consider two different approaches: the variation of the renormalization scale of the Higgs pole mass calculation and the variation of loop orders in the determination of the top Yukawa coupling.
Variation of the pole mass scale. First we discuss the variation of the renormalization scale at which the pole mass M h is computed in the SM. By default the scale Q pole = M t is chosen, which we vary by factor of two, This procedure estimates the impact of missing logarithmic higher-order corrections to the Higgs pole mass shift in the SM.
Variation of the loop order of threshold corrections at the low scale. As described in section 6, the relation between low-energy observables and MS-renormalized SM couplings contains corrections that can be switched off in the calculation without reducing the precision of the result for M h . As was shown in refs. [41,42,47,61], the dominant uncertainty obtained from this procedure is driven by the higher-order threshold correction in the relation between the top quark pole mass and the top Yukawa coupling. We define our estimation of missing threshold corrections at the electroweak scale in accordance with that reference as indicates that eq. (6.5) is evaluated at n-loop level. Since the consistent resummation of NNLL/N 3 LL logarithms requires an evaluation of eq. (6.5) at 2-/3-loop level, we estimate the uncertainty of the FEFT 2ℓ/3ℓ calculation by ∆M yt,2ℓ h and ∆M yt,3ℓ h , respectively.

Numerical size of individual uncertainties
In figure 8 we show the individual sizes of the uncertainty estimates discussed above for the parameter scenarios from figure 4. The two black lines correspond to uncertainties for FEFT 2ℓ and the other lines correspond to FEFT 3ℓ.

High-scale uncertainty
We start with a discussion of the high-scale uncertainty, shown in the left column of figure 8.
Estimate of missing 3-loop QCD and y t -enhanced contributions beyond O(y 4 t g 4 3 x ≤4 t ). In figure 8, the black and red solid lines represent the matching-scale uncertainties ∆M Q match h of the 2-and 3-loop FEFT calculation, respectively. The matching-scale uncertainty provides a global estimate of many kinds of terms, see eq. (8.2). The difference between the black and red solid lines corresponds to the inclusion of the known leading-QCD 3-loop contributions of O(y 4 t g 4 3 x ≤3 t ) to ∆λ. We find that this inclusion reduces the uncertainty very little, less than 0.2 GeV for all studied scenarios. In particular, since terms of the order O(y 4 t g 4 3 x 4 t ) are not known for all parameter scenarios, we expect a remaining uncertainty of significant size for large |x t | (see lower left panel of figure 8). Note, that ∆M Q match h is sensitive to terms of O(y 4 t g 4 3 x 4 t ) (c.f. eq. (8.3)) and thus includes an estimate of these missing terms.
To provide a direct estimate the size of the missing non-logarithmic 3-loop QCD and y t -enhanced contributions, we show as black and red dashed-dotted lines the uncertainties ∆M imp,g 3 yt,2ℓ h and ∆M imp,g 3 yt,3ℓ h for the 2-and 3-loop FEFT calculations, respectively. We find that these QCD uncertainties are very small already for the 2-loop calculation, ∆M imp,g 3 yt,2ℓ h ≲ 0.1 GeV. This is fully in line with the small difference between the 2-loop and 3-loop matching-scale uncertainty described above. The 3-loop QCD and y t -enhanced corrections missing in FEFT 3ℓ, ∆M imp,g 3 yt,3ℓ h , including terms with fewer powers of g 3 , are found to be negligible (red dashed-dotted line).
Taken together, all these results provide strong evidence that the contributions of leading QCD-type are already very well under control and inclusion of higher-order leading-QCD threshold corrections of O(y 4 t g 6 3 ) will not improve the precision of the calculation significantly.    predict that the uncertainty decreases at the same rate when going to higher M S . 22 This indicates that missing electroweak 2-loop terms contribute a theory uncertainty which is typically around 0.2-0.3 GeV, has a weak x t -dependence, and which is the dominant theory uncertainty for small |x t |.
Relevant higher-order contributions for large |x t |. For large |x t | ≳ 2 the matchingscale uncertainty ∆M Q match h is larger than for small x t . This cannot be attributed exclusively to the missing leading-QCD and 2-loop electroweak terms discussed so far. As discussed in section 8.1 and ref. [54], this is not unexpected because of the low powers of x t appearing in ∆M imp,g 3 yt,3ℓ h . On the other hand, the increased uncertainty for large x t is in line with the discussion of the impact of non-resummed large-x t contributions in section 7.3. In order to estimate missing terms with high x t -dependence, we employ the uncertainty estimates based on reparametrization terms. Indeed, reparametrization terms ∆M rep,g 3  For |x t | ≈ 3, the uncertainty estimate obtained from reparametrization becomes dominant. However, this is not specific to performing the calculation in full-model parametrization and it cannot be interpreted as an indication that the EFT parametrization would perform better with regard to missing contributions at O(ŷ 8 t +ŷ 6 tĝ 2 3 ). In fact, as discussed at the end of section 8.1, the estimation of the uncertainty for a calculation performed in SM parametrization would lead to a similar result at these orders. , is shown by the solid lines for FEFT 2ℓ (black solid line) and FEFT 3ℓ (red solid line). We find excellent agreement of the pole mass uncertainty of FEFT 2ℓ with the corresponding result shown in figure 3 of ref. [47]. 22 The middle-left plot in figure 8 shows a numerical instability at MS < 750 GeV in the FEFT 3ℓ calculation due to a hierarchy switch in Himalaya. The kink in the curves for ∆M Q match h in the lower left panel at xt ≈ −3.2 is due to a numerical artifact of our definition of the uncertainty. The irregularities at xt ≈ 2 in the middle-left plot are due to a numerical instability in our code, which is absent for lower values of tan β.

Low-scale uncertainty
Concerning the FEFT 3ℓ calculation we find a larger uncertainty of ∆M Q pole h than the corresponding FEFT 2ℓ calculation, which is surprising at first sight. The reason for this is the inclusion of the 3-loop Higgs pole mass shift in the SM of O(v 2ŷ4 tĝ 4 3 ) from ref. [70], which has the particular property that it increases the sensitivity of the Higgs pole mass on renormalization scale, if the scale is varied within Q pole ∈ [M t /2, 2M t ]. However, if the scale Q pole is varied within a larger range, the inclusion of this 3-loop correction leads to a significantly reduced dependence of the Higgs pole mass on Q pole . In order to keep our results comparable with the literature, we stick to the convention of using the Q pole ∈ [M t /2, 2M t ]. As a result, we find ∆M with the corresponding result from figure 3 of ref. [47], we find excellent agreement. Compared to the FEFT 2ℓ calculation, the FEFT 3ℓ calculation has a strong reduction of the uncertainty with ∆M yt,3ℓ h ≲ 0.2 GeV. This is the main source of the improved precision of our 3-loop calculation of M h in the studied scenarios.

Combined Uncertainty
In this subsection we combine the individual uncertainty estimates presented in the previous subsections to obtain a total uncertainty estimate of our new 2-loop and 3-loop Flexible-EFTHiggs calculations. Since the individual uncertainty estimates at the high-and lowenergy scales are sensitive to an overlap of higher-order terms, we define the following combined high-scale uncertainty, ∆M HS h , and low-scale uncertainty, ∆M LS h , for the FEFT nℓ calculation: h are sensitive to an overlap of higher-order contributions toλ that involve terms of O(y n t g m 3 ), we take their maximum in eq. (8.14). On the other hand, the electroweak contributions ∆M imp,g 1,2 ,2ℓ h are an independent subset of higher-order terms that involve electroweak gauge couplings, so we add it linearly to the other terms in eq. (8.14). To obtain the total combined uncertainty, ∆M h , of our calculations, we add the high-scale and low-scale uncertainties linearly,   For the degenerate SUSY mass scenarios defined in section 7, the results of our combined uncertainty estimates are shown in figure 9. The red solid line represents the Higgs pole mass M h obtained with the FEFT 3ℓ calculation and the red band in the lower sub-plots denotes the corresponding combined uncertainty ∆M h . The black dashed lines correspond to the FEFT 2ℓ calculation accordingly. The difference between the FEFT 3ℓ and 2ℓ calculations is of the order |M 3ℓ h − M 2ℓ h | ≲ 0.3 GeV. Compared to the 2-loop calculation, we find a more pronounced decrease of the uncertainty of the 3-loop calculation for large stop mixing |x t | ∼ 2 and M S ≳ 5 TeV. The dominant reduction of the total uncertainty of the 3-loop calculation is achieved in the low-scale uncertainty ∆M LS h , where ∆M yt,2ℓ h is the dominant uncertainty of the 2-loop calculation.
In general we find for the studied degenerate DR ′ SUSY mass parameter scenarios a combined uncertainty of the FEFT 3ℓ calculation of ∆M h ≲ 1 GeV for M S ≳ 1 TeV and |x t | ≲ 3. 23 This combined uncertainty becomes smaller for |x t | → 0 and larger M S , where it can reach ∆M h ∼ 0.5 GeV. These findings are compatible with the uncertainty estimates of refs. [60,61], where hybrid calculations with a comparable precision were studied. For large SUSY scales of M S ≳ 5 TeV we find that the remaining uncertainty of the FEFT 3ℓ calculation is dominated by the low-scale uncertainty induced by the determination of the top Yukawa coupling and the electroweak part of the high-scale uncertainty, which can be of similar size.

Conclusions
We have presented an extension of the FlexibleEFTHiggs method to calculate the SM-like Higgs pole mass in the MSSM. The method combines the virtues of an EFT and fixedorder calculation, resulting in a prediction that includes power-suppressed corrections and a resummation of large logarithms. We have applied our method to perform a state-of-theart calculation of the light CP-even Higgs pole mass in the MSSM, including corrections up to the 3-loop level and resummation of large logarithmic corrections up to N 3 LL.
The key of our extension is the generation of a consistent automatized pole mass matching procedure beyond the 1-loop level. The consistency of the FlexibleEFTHiggs method in this regard refers the cancellation of large logarithmic loop corrections and the inclusion of power-suppressed contributions in the matching to the EFT (assumed to be the Standard Model), thereby avoiding problems of double counting. Conceptually, this was achieved by a paradigm shift where the usually applied EFT-parametrized formulation of the high-scale matching was replaced by a parametrization in terms of full-model (MSSM) parameters. Technically, it required the inclusion of derivatives of the SM self energies and tadpoles in the multi-loop matching relations as described in section 4.
A thorough study of the new full-model parametrization shows that the new approach automatically resums leading contributions in the stop-mixing parameter x t , analogously to the well known tan β-resummation. This x t -resummation leads to significantly stabilized convergence of the perturbation series. For instance, in standard parameter scenarios such as in figure 3 and 4, the numerical impact of the known 2-loop (gaugeless) and 3-loop (leading QCD) threshold corrections is reduced to less than ∼ 0.3 GeV, compared to an impact of order 0.5-1.5 GeV in EFT-parametrized calculations.
Next, we have performed a detailed analysis of missing higher-order contributions of our 3-loop FlexibleEFTHiggs calculation. We have employed several different methods of uncertainty estimates, which have a complementary sensitivity to different types of missing higher-order contributions. Our analysis indicates that the remaining theory uncertainty of our calculation is dominated by (i) missing loop corrections to the top Yukawa coupling at the electroweak scale and (ii) missing electroweak 2-loop corrections to the quartic Higgs coupling at the SUSY scale, as shown in figure 8. Numerically, we find that the remaining theory uncertainty of our 3-loop FlexibleEFTHiggs calculation amounts to ∆M h ≲ 1 GeV 23 Note, that ∆M h is a measure of missing higher-order corrections in the relation between the predicted light CP-even Higgs pole mass and the DR ′ input parameters. As was stressed in ref. [61], there are additional uncertainties when the DR ′ input parameters are related to other physical observables.
for SUSY scales above 1 TeV and a stop-mixing of x t ≲ 3. This uncertainty is reduced to ∆M h ∼ 0.5 GeV for vanishing stop-mixing and/or SUSY scales of M S ≳ 10 TeV. Finally, we note that the resummation effects might be of high relevance for nonminimal supersymmetric extensions of the Standard Model, where the loop corrections to the Higgs mass are not known to the same order as in the MSSM. There, the matching correction in the full-model parametrization at NLO, for example, would result in a resummation of highest stop-mixing contributions of O(ŷ 2 t (ŷ 2 t +ĝ 2 1,2 )ĝ 2n 3 ) with n > 1, making resummation effects more advisable.

B Addendum
The following text in this section has been published by the authors as a separate addendum [98] to this publication. We include it here for convenience.

B.1 Introduction
We present the C++ program MSSMEFTHiggs3L, which implements the 3-loop Flexible-EFTHiggs state-of-the art calculation of M h in the real MSSM at N 3 LL and N 3 LO with x q resummation. The program is based on the FlexibleSUSY model NUHMSSMNoFVHimalaya and implements the matching and running described in our original publication, thus reproducing the results presented there. The program provides an easy-to-use SLHA interface for the MSSM input parameters and prints the value of M h as a single number to stdout.
We have structured the addendum as follows. In section B.2 we describe the technical details relevant for building the program. In section B.3 we discuss the user interface and relevant configuration options. Finally, we comment on the upcoming integration of the refined FlexibleEFTHiggs approach with full-model parametrization into the general FlexibleSUSY package.

B.2 Installation and usage of the stand-alone code
The MSSMEFTHiggs3L program can be downloaded as compressed package from https://flexiblesusy.hepforge.org/downloads/FlexibleEFTHiggs/MSSMEFTHiggs3L.tar.gz To build MSSMEFTHiggs3L, the boost C++ library, the Eigen3 library, the GNU Scientific Library and the Himalaya library [33] (version 4.0.0 or higher) are required. For installation instructions of the Himalaya library see e. g. ref. [33].
After the package has been extracted, it can be configured and compiled by running the following commands: $ ./ configure --enable -himalaya --enable -fflite \ --with -himalaya -incdir = $ { HIMALAYA_DIR }/ include \ --with -himalaya -libdir = $ { HIMALAYA_DIR }/ build \ --with -models = NUHMSSMNoFVHimalaya $ make The variable HIMALAYA_DIR contains the path to Himalaya root directory, required for the 3-loop pole-mass matching. Due to an improved numerical robustness, we recommend the configuration with the shipped 1-loop integral library FFLite. For more options see ./configure -h. After the compilation has finished, the program can be run with the shipped SLHA input file as follows: The multi-loop contributions entering the Higgs mass calculation are controlled by the configuration options in the FlexibleSUSY block of the SLHA input. A detailed documentation of the flags is given in ref. [56]. Here, we discuss the relevant options in the FlexibleSUSY block of the SLHA input, which controls the individual corrections of the Higgs pole mass calculation. Depending on the desired precision of the Higgs pole mass calculation, we present two configurations.