Pole mass determination in presence of heavy particles

We investigate the determination of the Higgs-boson propagator poles in the MSSM. Based upon earlier works, we point out that in case of a large hierarchy between the electroweak scale and one or more SUSY masses a numerical determination with $\overline{\text{DR}}$ Higgs field renormalization induces higher order terms which would cancel in a more complete calculation. The origin of these terms is the momentum dependence of contributions involving at least one of the heavy particles. We present two different methods to avoid their appearance. In the first approach, the poles are determined by expanding around the one-loop solutions. In the second approach, a"heavy-OS"Higgs field renormalization is employed in order to absorb the momentum dependence of heavy contributions. We will find that the first approach leads to an unwanted behaviour of the Higgs boson mass predictions close to crossing points where two Higgs bosons that mix with each other are almost mass degenerate. These problems are avoided in the second approach, which became the default approach used in the public code FeynHiggs. Despite the discussion being very specific to the MSSM, the argumentation and the methods presented in this work are straightforwardly applicable to the determination of propagator poles in other models involving a large mass hierarchy.


Introduction
The mass measurement of the Higgs boson discovered at the Large Hadron Collider (LHC) by the ATLAS and CMS experiments fixed the last free parameter of the Standard Model (SM) [1][2][3]. It is, however, well known that beyond the SM (BSM) physics is needed to tackle several unresolved problems of both experimental and theoretical nature. One of the most common models for BSM physics is the Minimal Supersymmetric Standard Model (MSSM) [4,5], based upon the concept of supersymmetry (SUSY). Apart from adding a superpartner to each SM degree of freedom, it also introduces a second Higgs doublet leading to five physical Higgs bosons: At the tree-level, these are the CP-even h and H boson, the CP-odd A boson and the charged H ± bosons.
One of its unique features is the fact that the Higgs sector is highly predictive. Using one of the non-SM-like Higgs masses as an input, the masses of the remaining Higgs bosons, of which one plays the role of the Higgs boson discovered at the LHC, are calculable in terms of only a few relevant parameters. In consequence, the mass of the Higgs boson discovered at the LHC can be used as a precision observable to constrain the MSSM parameter space. In order to profit from the high experimental precision, it is crucial to take radiative corrections into account. In the most direct approach, corrections to the Higgs-boson propagators are evaluated by the means of a fixed-order diagrammatic calculation (see [6][7][8][9][10] for recent works and references therein). In this approach, large logarithms appear in the calculation if one or more of the SUSY particles become much heavier than the electroweak scale. These can spoil the convergence of the perturbative expansion. Therefore, effective field theory (EFT) techniques are used to resum these logarithms (see [11][12][13][14][15] for recent works and references therein). Since typically no higher dimensional operators are included in the effective Lagrangian, 1 EFT calculations miss terms which would be suppressed in case of a heavy SUSY scale. To obtain a precise predictions for low, intermediary and high SUSY scales, hybrid calculations combining the two approaches have been developed [14,[16][17][18][19][20][21].
In this paper, we will focus on the hybrid approach which is implemented into the publicly available code FeynHiggs [16,17,20,[22][23][24][25][26][27]. Its main idea is to add the logarithms resummed by means of an EFT calculation directly to the fixed-order result, namely the diagrammatic corrections to the Higgs-boson propagators. Appropriate subtraction terms ensure that double-counting of logarithms contained in both results is avoided.
The physical Higgs masses squared are then determined by finding the poles of the Higgs-boson propagator matrix as calculated in the fixed-order approach (supplemented by the logarithms resummed in the EFT approach) and taking their real part. This pole determination is the main topic of this paper. We discuss three different methods.
In the first method, used in FeynHiggs up to version 2.13.0, the poles are determined by finding them numerically employing a DR renormalization of the Higgs fields. As already found in [20], the momentum dependence of contributions involving heavy superpartners induces higher order terms which would cancel in a more complete calculation.
In [20], a solution for this issue was proposed in the limit of all non-SM-like Higgs bosons being heavy by determining the Higgs-boson propagator poles at a fixed order. As a second method, we will extend this approach to be applicable also for the case of light non-SM-like Higgs bosons. We will show, however, that the truncation at a fixed-order leads to unwanted discontinuities in the prediction of the Higgs boson masses close to crossing points, where two Higgs bosons that mix with each other are almost mass degenerate.
Therefore, we will present a third method not suffering from the problems of the other two methods. It resembles the first method since the poles are determined numerically. In contrast to the first method, the Higgs field renormalization is, however, not renormalized in the DR scheme. Whereas the Higgs field renormalization drops out of the calculation at every complete order of the calculation, it is used at higher incomplete orders to cancel the momentum dependence of contributions involving heavy superpartners. We refer to this renormalization scheme as "heavy-OS" scheme. A similar renormalization scheme has already been introduced in [14] in order to combine an EFT calculation suitable for light non-SM-like Higgs bosons with the fixed-order calculation implemented in FeynHiggs. In this paper, we extend this scheme to the case of scenarios involving complex CP-violating parameters and apply it to the pole determination. This method is employed by default in FeynHiggs since version 2.14.3. As we will show, its implementation was crucial for the precise prediction of the Higgs boson masses in some of the MSSM Higgs benchmark scenarios presented in [28].
The paper is structured as follows. In Section 2, we describe the MSSM Higgs sector and its renormalization. In Section 3, we shortly review the decoupling behaviour of heavy contributions in general theories (without being specific to the MSSM). The consequences for the Higgs pole determination in the MSSM are then discussed in Section 4. Numerical results are presented in Section 5.

The MSSM Higgs sector and its renormalization
In this Section, we present a short overview of the MSSM Higgs sector and its renormalization. We restrict our discussion to the one-and two-loop level. At the two-loop level, we neglect the electroweak gauge couplings and assume vanishing external momentum. In addition, we neglect all first and second generation as well as the tau Yukawa coupling. This corresponds to the accuracy level of the fixed-order corrections implemented as default in FeynHiggs [22][23][24][25][26][29][30][31][32][33][34][35][36][37][38][39][40][41]. 2 The discussion is especially focused on the Higgs field renormalization, which is of particular importance in this paper.

Higgs sector at the tree level
The Higgs sector of the MSSM consists out of two Higgs doublets H 1 and H 2 with opposite hypercharge. They are given by where v 1 and v 2 are the vacuum expectation values (vevs) of the doublets. The ratio of the two vevs is denoted by The mass matrices of the component fields are diagonalized by the means of unitary transformations yielding the mass eigenstates. The transformations can be written as , where Explicit expressions for these one-and two-loop counterterms can be found in App. A. As a consequence of Eq. (13), all field renormalization constants in the CP-even subsector, involving h and H, get contributions only from the real parts of Z H . The same holds true for the CP-odd subsector, involving A and G, and the diagonal field renormalization constants in the charged Higgs sector. In contrast, the field renormalization constants connecting the CPeven and the CP-odd subsectors (e.g. δZ hA ) receive contributions only from the imaginary parts of Z H . The off-diagonal field renormalization constants in the charged Higgs sector get contributions from both the real and imaginary parts of Z H .

Renormalization of tanβ
The vevs v 1 and v 2 are renormalized according to The appearance of the field renormalizaton constants follows from the renormalization prescription in Eq. (11).
The difference of the one-loop vev counterterms is UV finite [29,30]. Therefore, we set In the approximation of vanishing electroweak gauge couplings, which we apply at the twoloop level, also the difference of the two-loop vev counterterms is UV finite [43]. So, we set where the subscript "gl" is used to indicate the limit of vanishing electroweak gauge couplings, which is also often called gaugeless limit. Applying these simplifications, the resulting one-and two-loop counterterms of tan β, introduced via the renormalization transformation are then given by

Renormalization conditions for the Higgs fields
The UV divergent part of the field renormalization constants is fixed by demanding that the derivatives of the Higgs self-energies with respect to the external momentum are UV finite, where we used the subscript "div" to denote that only the UV divergent part is taken into account. We do not specify the momentum at which the derivatives of the self-energies are evaluated, since their UV divergent part does not depend on the external momentum. As noted already above, the off-diagonal field renormalization constants are not needed to render Green's functions with external Higgs bosons UV finite. In the addition to the UV divergent part, we also allow for a UV finite contribution to the field renormalization constants. This part will be specified in Section 4. Note that the specification of the field renormalization constants also fixes the renormalization and therefore definition of tan β.
In the approximation of vanishing electroweak gauge couplings and vanishing external momentum, as applied at the two-loop level, all two-loop field renormalization constants drop out of the calculation [39]. They also do not appear at higher-orders, since we work in the approximation of vanishing external momentum at the two-loop level. Therefore, we do not impose any condition on them in this work.

Renormalized self-energies
At the one-loop level, the renormalized self-energies are given in terms of the unrenormalized self-energies and counterterms bŷ where the various δ (1) m 2 denote the one-loop mass counterterms, which are given explicitly in App. A. The renormalized self-energies involving Goldstone bosons are obtained analogously. The renormalized two-loop self-energies are given bŷ where δ (2) m Z ij is a combination of mass counterterms and field renormalization constants. Explicit expressions are listed in App. A. Note that the unrenormalized two-loop self-energies contain subloop renormalization contributions.

Momentum dependence of heavy contributions
In this Section, we review the decoupling behaviour of heavy particles focusing especially on the momentum dependence induced by them. Our discussion is based on the proof of the decoupling theorem (see [44]) presented in [45]. The arguments in this Section are generally applicable and not restricted to the MSSM.
We investigate the kinematic behavior of a theory with one light boson with mass m and one heavy boson with mass M at energies p 2 , assuming 3 In order to study the momentum dependence of loop corrections involving the heavy particle, we consider the loop corrections to the unrenormalized n-point vertex function Γ (n) with external light fields only. The mass dimension of Γ (n) , depending on the number of external light bosons, is given by where we assume that the mass dimension of the overall Lagrangian is four. The square brackets are used to denote the mass dimension of the quantity within the brackets.
Next, we separate the n-point vertex function into a light and a heavy part, where the light part is assumed to contain all loop corrections involving only the light particle; the heavy part, to contain loop corrections involving at least one heavy particle. The heavy part can be expressed as where the sum runs over all involved loop diagrams. The involved loop integrals are labelled by I i 's are dimensionless prefactors (e.g. gauge couplings). 4 Consequently, the mass dimension of the integrals is the same as of the overall n-point vertex function, We want to know how this integral behaves in the limit of large M . Integrals can be expanded in the limit of an internal mass becoming heavy by using the technique of expansion by regions [46][47][48][49][50]. The leading contribution in this limit is given by the region in which all loop momenta are similar in size as the heavy mass M . In this region, all external momenta and light internal masses can be neglected. Consequently, the leading piece in the large mass limit can depend on neither the external momenta nor the light internal masses. By dimensional analysis, we obtain where the dots denote terms subleading in the M → ∞ limit. b

(n)
i,0 is a dimensionless coefficient, which does neither depend on m nor on p 2 . It can, however, contain logarithms involving the heavy mass and the renormalization scale µ R . 5 To investigate the momentum dependence, we take the derivative of the n-point vertex function with respect to the external momentum. The derivative lowers the mass dimension, with k ≥ 1. Consequently, we have d dp 2 and for the n-point vertex function This implies that the leading contribution in the M → ∞ limit is given by if there are more than two external particles (n > 2).
In case of only two external particles, there is a finite non-zero remainder if only one derivative is applied, whereas d dp 2 for k > 1. This means that in the limit M → ∞, the heavy parts of n-point vertex functions do not depend on the external momentum if there are more than two external particles, for n > 2. If there are only two external particles, the dependence on the external momentum is given by where we expanded around p 2 = m 2 . Note that we did not specify at which external momentum the term d dp 2 Γ The remaining dependence on the external momentum in the second line of Eq. (38) is not affecting physical observables. In fact, it can be removed by renormalizing the two-point vertex function and absorbing the remaining dependence on the external momentum into the field renormalization of the light external boson (see [45] for more details), The UV divergent part is chosen independently of the UV finite part such that the derivative of the two-point vertex function with respect to the external momentum is finite. In contrast to the OS scheme, in which the derivative of the full self-energy, including "light" and "heavy" contributions, is renormalized to zero, here only the "heavy" part is absorbed. Therefore, we will refer to this scheme as "heavy-OS scheme" in the following. In general, field renormalization constants drop out in the calculation of physical observable if all corrections of a specific order are completely taken into account. Consequently, the momentum dependence of heavy contributions vanishes in the limit of the heavy mass going to infinity -as expected from the decoupling theorem -at each complete order independently of the chosen renormalization scheme, since the result has to be independent of this choice. At higher incomplete orders, we have to impose Eq. (39) to ensure that the momentum dependence of heavy contributions is canceled out.
Following a very similar reasoning, also the renormalization scale dependence of the twopoint vertex function can be assessed. Analysing the degree of divergence after applying derivatives with respect to the external momentum shows: Leaving aside the intrinsic renormalization scale of the parameters entering the derivative of the two-point vertex function, also its explicit scale dependence can be absorbed into the Higgs field renormalization.

Higgs pole determination in the MSSM
Next, we will discuss how the decoupling behaviour of heavy particles, as discussed in Section 3, affects the determination of the Higgs-boson propagator poles in the MSSM. Note that the discussion, despite being specific to the MSSM here, is actually applicable to general models involving large mass hierarchies between two or more particles.
The Higgs-boson propagator poles, whose real parts are the squared physical masses, are obtained by solving the equations with the two-point vertex functionŝ We neglect mixing with the Goldstone and the Z boson, since these induce only subleading two-loop contributions. In order to keep the expressions short, we do not explicitly write down the contribution coming from the EFT calculation but assume that it is intrinsically contained in the self-energies.
In the following, we will discuss different methods for solving Eqs. (40).

Numerical pole determination with DR field renormalization
The most direct approach to determine the Higgs-boson propagator poles is to solve Eqs. (40) numerically. This can be done employing different schemes for the Higgs field renormalization. First, we discuss the case of DR Higgs field renormalization, which was used as default scheme in FeynHiggs until version 2.14.2. This means that no UV finite parts are added to the UV divergent parts as specified in Eqs. (22a).
In [20], this procedure was closely investigated in the limit of M A M Z . In this limit, h-H mixing effects in the CP-even Higgs sector are suppressed by powers of M A and the mass of the h boson can be obtained by solving the simpler equation up to corrections suppressed by M A . Solving this equation up to the two-loop order yields where the number in the superscript is used to indicate the loop order. Again, the subscript "gl" is used to indicate the approximation of vanishing electroweak gauge couplings at the two-loop level. In the last step, we split up the derivative of the hh self-energy into a part containing all contributions involving heavy particles (omitting terms suppressed by the heavy mass scale), labelled by the superscript "heavy" particles, and a part containing all contributions from light particles (including terms suppressed by a heavy mass scale), labelled by the superscript "light". Note that, in order to simplify the notation in this Section, the labels have a slightly different meaning in this Section as compared to Section 3. Here, the light contribution also includes terms suppressed by a large mass, whereas these terms are omitted in the heavy contribution. In [20], a scenario with all non-SM particles sharing one common mass scale was assumed. Therefore,Σ light hh contains contributions from SM particles as well as terms suppressed by a heavy non-SM mass.
The analysis in [20] showed that the term labelled by "heavy" in Eq. (43) is cancelled by parts of the subloop renormalization contained in the two-loop self-energy. The discussion in Section 3 clearly shows the reason for this cancellation: The momentum dependence of heavy contributions to two-point vertex functions, i.e. the term proportional toΣ (43), could be absorbed into a field renormalization constant and therefore drops out order by order in the perturbative expansion. 6 If the Higgs field renormalization is fixed in the DR scheme, it is therefore important for a proper cancellation to ensure that all corrections at the two-loop level are included at the same level of accuracy and that no higher order terms proportional toΣ Determining the Higgs-boson propagator poles numerically, however, does not allow to easily control the inclusion or exclusion of certain terms.
Concerning the corrections included in FeynHiggs, this means that in order to achieve a complete cancellation at the two-loop level only O(α 2 t , α t α b , α 2 b ) corrections should be taken into account for the termΣ hh (m 2 h ), since the two-loop self-energy is included up to this order.

Fixed-order pole determination
As an approach to solve this problem, it was suggested in [20] to simply employ Eq. (43) and to ensure that the last term in the equation is computed at the same level of accuracy as the two-loop self-energy. This simple procedure is, however, only applicable in case of high M A for which h-H mixing effects are negligible. Here, we work out a method based upon the same idea -solving Eqs. (40) at a fixed-order -applicable also for low M A . This method was used in the FeynHiggs versions 2.14.0 through 2.14.2.
A proper cancellation is guaranteed by including all two-loop corrections at the same level of accuracy. In case of low M A , we cannot neglect the off-diagonal elements of the Higgs two-point vertex function, since their inclusion is crucial for a precise prediction. Instead, we consider each element of the two-point vertex function separately.
The complex pole corresponding to the i-th boson (i = 1,2,3) is given at the one-loop level by where the superscript indicates the loop-order. We choose the tree-level masses to be and the one-loop self-energies to bê

=Σ
(1) Expanding each of the self-energies appearing in the two-point vertex function matrix (see Eq. (41a)) around the one-loop solution of Eq. (44) up to the two-loop level yields the expanded two-point vertex function matrixΓ i−exp hHA with the matrix elements given by with j, k = 1, 2, 3 and The subscript of the reducible termsΣ (1) Σ (1) , "gl", indicates that we calculate these in the gaugeless limit (i.e. vanishing electroweak gauge couplings corresponding to the approximation used at the two-loop level). Consequently, these terms are of O(α 2 t , α t α b , α 2 b ) and therefore of the same order as the corresponding parts of the two-loop self-energiesΣ (2) . Thereby, it is guaranteed that cancellations between the various two-loop contributions are correctly taken into account. The next step is to find the poles ofΓ i−exp hHA (p 2 ). We pick the pole with the largest overlap of the associated physical state with the tree-level mass eigenstate. This pole corresponds to the one-loop solution given in Eq. (44). The physical mass squared is then given by taking the real part of this pole. The mass of the H ± boson is determined accordingly. This procedure successfully avoids inducing higher order terms which would cancel in a more complete calculation.
However, as we will argue in Section 4.4 and explicitly show in Section 5, the approach presented in this Section yields unsatisfactory results close to crossings points where two mixing Higgs bosons are almost mass degenerate. This issue does not appear if the propagator poles are determined numerically.

Numerical pole determination with heavy-OS field renormalization
Last, we present a third approach to determine the Higgs-boson propagator poles avoiding the problems of the approaches presented above. As in Section 4.1, the Higgs-boson propagator poles are determined numerically. However, to avoid the problem with uncancelled terms beyond the order of the fixed-order calculation, we do not employ the DR scheme for the Higgs field renormalization. As shown in Section 3, the Higgs field renormalization can be used to absorb the momentum dependence of contributions from heavy particles. As discussed in Section 4.1, their momentum dependence is the origin of the uncancelled higher-order terms. Consequently, the issue is resolved by fixing the Higgs field renormalization according to the heavy-OS scheme (as presented in Section 3). This means that the one-loop Higgs field renormalization constants (introduced in Section 2) are fixed by Note that we did not specify at which external momentum the second term is evaluated, since it does not depend on the external momentum. In addition to absorbing the momentum dependence of the heavy contributions, we also add a term proportional to the Higgs anomalous dimension involving light fields only, 7 which is used to absorb the dependence of Higgs self-energies' derivatives on the renormalization scale µ R (see end of Section 3). In this way, we separate the scale dependence associated with the pole determination from the renormalization scale dependence of the fixed-order corrections. The corresponding reference scale at which the poles are determined is called Q pd . In order to avoid inducing large logarithms via the pole determination we choose it equal to the OS top mass, Q pd = M t . Explicit expressions for the one-loop renormalization constants in the gauge eigenstate basis are listed in App. C.
Since the fixed-order calculation used by default in FeynHiggs applies the approximation of vanishing external momentum at the two-loop level, no finite pieces have to be added to the two-loop field renormalization constants. Going beyond this approximation (see [6,10,42]) implies that also the two-loop field renormalization constants should be fixed in the heavy-OS scheme.
Having calculated the Higgs self-energies using the field renormalization scheme specified in Eq. (49), the Higgs-boson propagator poles are then determined by solving Eqs. (40) numerically.
Note that employing the full OS scheme, meaning that the full derivatives of the Higgs self-energies are renormalized to be zero, is not a good alternative. In the OS scheme, unphysical threshold effects coming from the "light" contributions are induced (see e.g. [51]).
It is important to note that the renormalization scheme choice made in Eq. (49) affects the definition of tan β (see Section 2). In particular tan β is not anymore a MSSM quantity but rather a Two-Higgs-Doublet Model (THDM) one. This is explained in detail in [14]. By choosing the field renormalization as in Eq. (49), tan β is defined at the scale Q pd , which is set to the OS top mass. In situations where the THDM scale M A is much larger than the top-quark mass, this definition is not well-suited. Furthermore, it would be preferable to have an MSSM tan β as input for a MSSM calculation. We present two methods to address these shortcomings in App. B.
In addition to the Higgs poles, also the Z-matrix relating the tree-level mass eigenstates and the external physical states in scattering amplitudes is calculated [26,29,[52][53][54][55]. This matrix depends on the choice of the Higgs field renormalization. For other calculations (e.g. decay rates or productions cross-sections), the heavy-OS scheme is, however, not needed to avoid uncancelled higher order terms, since these processes involve more than two external legs. In order to ease the use of the output Z-matrix, we therefore discuss how it can be transferred to the DR scheme, which is more widely used, in App. D.

Comparison of the different methods
To further illustrate the difference between the various approaches, we investigate the corrections to the mass of h boson up to the two-loop level.
In the first approach, numerical pole determination with DR Higgs field renormalization, they are given by In the fixed-order approach, we obtain We see that the fixed-order pole determination avoids including the terms identified in Section 4.1, which would cancel in a more complete calculation. Using numerical pole determination with OS pole determination, we get Note that in this Equation, the hat, used to mark renormalized quantities, refers to another renormalization scheme for the Higgs fields, namely the heavy-OS scheme, than in Eqs. (51) and (52), where the Higgs field renormalization constants are fixed employing the DR scheme. Again, terms which would cancel in a more complete calculation are not included. The difference between the numerical pole determination employing heavy-OS Higgs field renormalization and the fixed-order pole determination is given by 8 We see that the numerical pole determination employing heavy-OS renormalization includes additional terms. The term containing the light contribution to the derivative of the hh selfenergy and the hh self-energy is included fully and not only in the limit of vanishing electroweak gauge couplings. The electroweak next-to-leading logarithms which are contained in this term are of the same order as some of the logarithms resummed by the means of the EFT calculation, which are added on top of the fixed-order result. Therefore, this term should be included fully and not only in the limit of vanishing electroweak gauge couplings.
In addition to this differences at the two-loop level, there are also differences at higher orders. In the case of fixed-order pole determination, unphysical discontinuities in the prediction of the Higgs boson masses close to crossings points, where two mixing Higgs bosons are almost mass degenerate, appear (numerical results displaying this behaviour will be shown in Section 5). At these crossing points, a specific property of the Higgs-boson propagator matrix becomes important: Every complex pole is not only the solution of Eqs. (40) but also a pole of every element of the Higgs-boson propagator matrix [53]. This property is fulfilled in case of numerical pole determination. In case of the fixed-order pole determination, this property, however, holds only up to the level of the truncation, i.e., the two-loop level. The violation at higher orders is the origin of the discontinuities in the Higgs mass predictions.

Numerical results
In this Section, we will numerically compare the different methods to determine the Higgsboson propagator poles.
First, we investigate scenarios with M A = 200 GeV. All other soft-SUSY breaking masses are chosen to be equal to the common mass scale M SUSY , which we vary between 0.5 TeV and 10 TeV. The trilinear soft-breaking couplings are set equal to zero apart from A t , which is fixed by setting X OS t /M SUSY = 2. 9 All parameters are chosen to be real. Note that large parts of this parameter space are experimentally excluded due to the small M A value (see e.g. [28]). Here, we, however, only want to explore the numerical size of the differences between the various pole determination methods and do not perform a phenomenological analysis.
In the left plot of Fig. 1, we vary M SUSY setting tan β = 1.1 (this choice maximises the difference between the pole determination methods). Shown are the differences between 8 Note thatΣ (1) hh (m 2 h ) does not depend explicitly on the scheme chosen for the Higgs field renormalization. Assuming the same definition of tan β,Σ (1) hh (m 2 h ) is therefore equal for DR Higgs field renormalization and for heavy-OS Higgs field renormalization. 9 Here, we employ the OS scheme for the renormalization of the stop sector. In the right panel of Fig. 1, we show the dependence of the shifts between the various pole determination methods on tan β for M SUSY = 2 TeV. The lines correspond to those of the left plot of Fig. 1. For tan β 3, the shifts of M h are approximately constant and equal to each other. In the region tan β 3, the difference between the results for M h using numerical pole determination with heavy-OS Higgs field renormalization and fixedorder pole determination rises when lowering tan β. This behaviour reflects the decreasing importance of Higgs mixing effects if tan β is increased. Due to the suppression of the top Yukawa coupling of the H boson for high tan β, the shifts of M H quickly decrease if tan β is raised. The shifts of M H ± show no significant dependence on tan β.

Next, we investigate the M 125
H benchmark scenario proposed in [28]. It is defined such that the second lightest CP-even Higgs boson is SM-like and plays the role of the Higgs boson discovered at the LHC. Whereas we refer to [28] for the exact definition of the scenario, we want to point out the large µ/M SUSY value of ∼ 8 with M SUSY ∼ 750 GeV. We set There size is especially large here because of the large µ/M SUSY value. In most parts of the MSSM parameter space, where the lightest CP-even Higgs is SM-like, the shifts between different pole determination methods are much smaller.
In the left plot of Fig. 2, also results using the fixed-order pole determination method are shown (red). The curve for M H has discontinuities at M H ± ∼ 186.5 GeV and M H ± ∼ 196 GeV by ∼ 1 GeV and ∼ 4 GeV, respectively. We investigate these discontinuities more closely in the right plot of Fig. 2. As explained in Section 4.2, in the fixed-order pole determination approach each element of the inverse Higgs-boson propagator matrix is expanded around the one-loop solutions for the poles. In the case of real input parameters, these are HH (m 2 H ) (green curves). After determining the eigenvalues of these expanded matrices, the overlap of the corresponding eigenstates with the tree-level mass eigenstates are determined. 10 For the solid curves, the eigenstate overlaps mostly with the tree-level mass eigenstate corresponding to the one-loop solutions around which we expanded. These are the correctly chosen eigenvalues and correspond to the red curve in the left plot of Fig. 2. The dashed curves show the wrongly chosen eigenvalues whose eigenstate overlaps mostly with the tree-level mass eigenstate not corresponding to the one-loop solutions. We observe that the discontinuities visible in the left plot of Fig. 2 does not originate from numerical instabilities but are a consequence of the expansion of the propagator matrix at the two-loop level, which generates two different solution branches for the pole determination.
For M H ± 186 GeV and M H ± 197 GeV, the results using fixed-order pole determination and numerical pole determination with heavy-OS Higgs field renormalization are in good agreement. This confirms that both methods avoid inducing higher order terms which would cancel in a more complete calculation and which are included in the numerical pole determination employing DR Higgs field renormalization. In the crossing point region, the numerical pole determination with heavy-OS Higgs field renormalization yields no unphysical discontinuities and is therefore superior to the fixed-order pole determination method. Consequently, this method was used for the definition of not only the M 125 H scenario but also of the other benchmark scenarios presented in [28]. 11 After investigating scenarios with light non-SM Higgs boson masses, we turn to a highscale scenario where all non-SM masses are set equal to M SUSY (i.e., M A = M SUSY ), which is varied between 1 TeV and 100 TeV. Here, we fix the parameters of the stop sector in the DR scheme.
In the left plot of Fig. 3, we show the prediction for M h obtained using numerical pole determination with DR Higgs field renormalization (blue), fixed-order pole determination (red) and numerical pole determination with heavy-OS Higgs field renormalization (green). In addition, we display the pure EFT result obtained with FeynHiggs (orange). All curves are shown for X DR t /M SUSY = 0 (solid) and X DR t /M SUSY = √ 6 (dashed). In the considered range for M SUSY , terms suppressed by the SUSY scale are negligible and we expect to see a good agreement between the EFT prediction and the hybrid approach implemented in FeynHiggs. We observe that the remaining numerical difference between the EFT result and the hyprid approach is lowest if numerical pole determination with heavy-OS field renormalization is employed. In comparison to the fixed-order pole determination, especially the terms isolated in Eq. (54) improve the agreement between the EFT and the hybrid result. The still remaining difference (∼ 0.2 GeV for vanishing stop mixing and ∼ 0.5 GeV for X DR t /M SUSY = √ 6) can be explained by missing contributions proportional to the bottom Yukawa coupling in the EFT calculation and a different parameterization of non-logarithmic terms [20].
Finally, we investigate the dependence of the results on the renormalization scale µ R . As mentioned above, for the fixed-order calculation of FeynHiggs the OS scheme is employed by default. Only the field renormalization and thereby tan β are renormalized using the DR scheme. Consequently, the results for M h should only depend on the renormalization scale via the dependence of tan β.
As argued in Section 4, the numerical pole determination method using DR Higgs field renormalization leads to the inclusion of higher order terms depending on the renormalization scale µ R which would cancel in case of a more complete fixed-order calculation. This incomplete cancellation is reflected by an enhanced dependence on the renormalization scale.
We illustrate this behaviour in right plot of Fig. 3. In this Figure, all non-SM masses are chosen equal to M SUSY . We compare the results of the different pole determination methods for two different scale choices: µ R = M t (solid), the default choice in FeynHiggs, and µ R = M SUSY (dashed). We set tan β = 4. Looking at the curves produced using the numerical pole determination with DR Higgs field renormalization (blue), we see a small discrepancy between the two scale choices already for low SUSY scales. For rising M SUSY , the discrepancy increases to up to 12 GeV for M SUSY ∼ 100 TeV. In contrast, the curves using the fixed-order pole determination (red) lie within 1 GeV from each other indicating a small dependence on µ R over the whole considered M SUSY range. This remaining dependence can be explained by the renormalization scale dependence of tan β which is the only input parameter depending the renormalization scale. If we employ numerical pole determination with heavy-OS Higgs field renormalization, tan β is defined at the fixed scale Q pd = M t . Therefore, the corresponding curves (green) lie on top of each other. Note that also here in principle higher order terms including a dependence on the renormalization scale are induced. As explained in Section 4.3 this renormalization scale dependence is, however, absorbed into the field renormalization constants and thereby replaced by a dependence on the pole determination scale Q pd , which is kept fixed. The fact that the results using different renormalization scale choices are equal to each other confirms that this scale separation works as expected.
These plots show that using numerical pole determination with heavy-OS Higgs field renormalization allows to independently control the different scales of the calculation, i.e., the scales of the pole determination and the scale of tan β. This is an important step towards an improved estimate of the remaining theoretical uncertainties.

Conclusions
In this paper we investigated the determination of the Higgs-boson propagator poles in the MSSM. First, we gave a short overview of the MSSM Higgs sector and its renormalization. We especially focused on the Higgs field renormalization.
Afterwards, we discussed the behaviour of general n-point vertex functions if the masses of one or more internal particles approach infinity. These contributions have a dependence on the external momenta only if the number of external particles is equal two. In this case, the normalization of the external fields can be used to absorb the dependence.
Next, we focused on the determination of the Higgs-boson propagator poles in the MSSM. We showed that the most direct approach, determining the poles numerically employing a DR Higgs field renormalization, induces higher order terms which would cancel in a more complete calculation. The origin of these terms is the momentum dependence of diagrams involving heavy internal particles.
As an alternative approach to avoid this problem, we introduced a method to determine the Higgs-boson propagator poles at a fixed-order. This ensures that no unwanted higher order terms are induced. Close to crossing points, where two Higgs bosons that mix with each other are almost mass degenerate, this methods, however, suffers from unphysical discontinuities in the mass prediction.
Therefore, we introduced a third method which employs a numerical pole determination but fixes the Higgs field renormalization such that the momentum dependence of heavy contributions is absorbed. We call this scheme heavy-OS Higgs field renormalization. This method avoids the inclusion of unwanted higher order terms and also does not suffer from unphysical discontinuities close to crossings points.
In our numerical investigation, we compared the different methods in scenarios with light non-SM Higgs bosons finding differences of up to ∼ 3 GeV for M SUSY ∼ 10 TeV and tan β ∼ 1. In addition, we closely investigated the predictions for the M 125 H benchmark scenario introduced in [28], finding the third method to yield a smooth prediction approaching the results of the fixed-order pole determination away from the crossing point. Moreover, we found that the numerical pole determination with heavy-OS Higgs field renormalization brings the prediction of the combined fixed-order and EFT calculation closer to the prediction of the pure EFT calculation for very high SUSY scales. Last, we showed that the complete dependence of the overall result on the renormalization scale is caused by higher order terms induced through the pole determination. This is an important input for a re-evaluation of the remaining theoretical uncertainty, which we leave for future work.
We want to point out that despite the discussion being very specific to the MSSM, the argumentation and the methods presented can straightforwardly be applied to the determination of propagator poles in other models which involve a large mass hierarchy.

A One-and two-loop counterterms
In this Appendix, we list explicit expressions needed for the renormalization of the MSSM Higgs sector as described in Section 2. The expressions extend the ones listed e.g. in [14,26,39] by including off-diagonal field renormalization constants (see Section 2.2.1) and independent counterterms for tan β and the mixing angles (see discussion in App. C).

A.1 One-loop renormalization
The one-loop field renormalization constants of the neutral Higgs bosons in the mass basis are given by The corresponding constants in the charged Higgs sector read The field renormalization counterterms in the gaugeless limit can be obtained by setting α = β − π/2 in the expressions listed above. The one-loop mass counterterms of the neutral Higgs sector are given by and for the charged Higgs sector by

A.2 Two-loop renormalization
The relations between the two-loop field renormalization constants of the neutral Higgs sector in the mass basis and in the original gauge basis are given by where The corresponding expressions for the charged Higgs sector read The field renormalization counterterms in the gaugeless limit can be obtained by setting α = β − π/2. The two-loop mass counterterms in the neutral Higgs sector are given by and for the charged Higgs sector by Note that the expression are only valid in the limit of vanishing electroweak gauge couplings. Correspondingly, also the counterterm of α has to be inserted applying the limit of vanishing electroweak gauge couplings. The remaining two-loop mass counterterms, δ (2) m 2 A and δ (2) m 2 H ± , are fixed by employing OS conditions for either the mass of the A boson or the mass of the H ± boson (and by using δ (2) . In order to obtain the renormalized two-loop self-energies (see Eq. (24)) the following combinations of mass and field counterterms are needed in the neutral Higgs sector

B Reparametrization of tanβ
The heavy-OS Higgs field renormalization fixes tan β to be defined in the THDM at the scale M t (see discussion in Section 4). Here, we discuss how a different definition for tan β can be achieved.
The most direct solution is to accept that tan β is defined as stated above. If the input tan β is then defined in another scheme, we have to convert it to the scheme of the calculation, The obtained value is then inserted as new input parameter into the Higgs mass calculation.
In practice, a one-loop conversion is employed. This procedure is easy. It, however, has the disadvantage that the conversion induces terms beyond the order of the fixed-order conversion which should not be included if the result of the fixed-order calculation is combined with an infinite series of resummed logarithms. 12 To avoid the problem with induced higher-order terms, we can expand the result e.g. for M 2 h around the input tan β, δ (2) t β gl,fin = 0.
(71c) At the one-loop level, the following relations between the counterterms of tan β and the mixing angles hold, and at the two-loop level, δ (2) t βn gl = δ (2) t βc gl = δ (2) t β free,gl .
We find that both methods yield numerically very similar results. Only in scenarios with low tan β and a large hierarchy between Q in and M t , we observe non-negligible, but still small, numerical differences (e.g. for Q in = M SUSY = 20 TeV and tan β ∼ 2 the difference in M h is 0.4 GeV). These are induced through the higher-order terms generated in case of a conversion of the input tan β.

C Heavy-OS renormalization -explicit expressions
Here, we list explicit expressions for the finite parts of the Higgs field renormalization constants in the heavy-OS scheme. Similar expressions have already been presented in [14]. In there, however, all parameters were assumed to be real and all SUSY soft-breaking scalar masses as well as the gaugino masses were assumed to equal. For abbreviation, we introduce the constant where α is the fine structure constant and s W is the sine of the electroweak mixing angle. The contribution from up-type squarks is given by In addition, there is a contribution from non-SM Higgs bosons, In addition to the explicit dependence on the Higgs field renormalization constants, there also is an implicit dependence through the tan β counterterm which appears in δ (1) m 2 hH (see Eq. (57b)). This dependence vanishes if the Z-matrix is combined with renormalized n-point vertex functions in order to calculate physical observables. In order to get a DR-renormalized Z-matrix, which can easily be used in other calculations, we have to multiply the Z-matrix calculated employing heavy-OS Higgs field renormalization with the corresponding Higgs field renormalization matrix. In addition, we have to subtract the pieces proportional to the tan β counterterm in the off-diagonal entries. This procedure is applied in FeynHiggs.