Two-Loop ${\cal O}(\alpha_t^2)$ Corrections to the Neutral Higgs Boson Masses in the CP-Violating NMSSM

We present our calculation of the two-loop corrections of ${\cal O}(\alpha_t^2)$ to the neutral Higgs boson masses of the CP-violating Next-to-Minimal Supersymmetric extension of the Standard Model (NMSSM). The calculation is performed in the Feynman diagrammatic approach in the gaugeless limit at vanishing external momentum. We apply a mixed $\overline{\mathrm{DR}}$-on-shell (OS) renormalization scheme for the NMSSM input parameters. Furthermore, we exploit a $\overline{\mathrm{DR}}$ as well as an OS renormalization in the top/stop sector. The corrections are implemented in the Fortran code NMSSMCALC for the calculation of the Higgs spectrum both in the CP-conserving and CP-violating NMSSM. The code also provides the Higgs boson decays including the state-of-the-art higher-order corrections. The corrections computed in this work improve the already available corrections in NMSSMCALC which are the full one-loop corrections without any approximation and the two-loop ${\cal O}(\alpha_t \alpha_s)$ corrections in the gaugeless limit and at vanishing external momentum. Depending on the chosen parameter point, we find that the ${\cal O}(\alpha_t \alpha_s + \alpha_t^2)$ corrections add about 4-7% to the one-loop mass of the SM-like Higgs boson for $\overline{\mathrm{DR}}$ renormalization in the top/stop sector and they reduce the mass by about 6-9% if OS renormalization is applied. For an estimate of the theoretical uncertainty we vary the renormalization scale and change the renormalization scheme and show that care has to be taken in the corresponding interpretation.


JHEP09(2021)193 1 Introduction
The Standard Model (SM) of particle physics belongs to the most successful theories ever tested. Despite its success the SM lacks explanations for a variety of open problems such as for example the nature of Dark Matter (DM) or the observed baryon-antibaryon asymmetry in the universe. Models based on supersymmetry (SUSY) [1][2][3][4][5][6][7][8][9][10][11][12][13][14] are promising beyond-the-SM (BSM) candidates that offer solutions to many problems the SM cannot address. In SUSY models the Higgs and gauge sectors are related so that in the Minimal Supersymmetric Standard Model (MSSM) [15][16][17][18] the tree-level mass m h of the lightest Higgs boson is bounded from above by the Z-boson mass m Z . Higher-order corrections, where those from the top/stop sector play the dominant role, can shift the upper bound to larger values. In order to reach the experimentally measured value of 125.09 GeV [19] for the SM-like Higgs boson mass, however, a large soft-SUSY breaking mass scale m SUSY and/or a large mixing in the stop sector is required so that naturalness arguments in favor of supersymmetry become questionable. In the Next-to-MSSM (NMSSM) [20][21][22][23][24][25][26][27][28][29][30][31] the situation is more relaxed as it features additional F -term contributions raising the tree-level Higgs mass to a higher value so that higher-order corrections can be smaller compared to the MSSM and still lead to the required Higgs mass value.
In the meantime, the Higgs boson mass has turned into a precision observable with an uncertainty of a few hundred MeV [19]. The precise knowledge of the Higgs boson mass is important as it is a crucial input parameter for all Higgs boson observables [32] and determines the stability of the electroweak vacuum [33][34][35]. Therefore, in order to make meaningful interpretations of the experimental results, the experimental accuracy has to be matched by the precision of the theory predictions. Only then sensible limits on the still allowed parameter space of the model can be derived from the experimental results and possibly distinguish between new physical models in case of discovery, cf. e.g. ref. [36]. Consequently, a tremendous effort has been put in the computation of the higher-order corrections to supersymmetric Higgs boson masses. These calculations can be grouped into three classes that are based on fixed-order (FO), effective field theory (EFT) or hybrid techniques that make use of FO as well as EFT results. A comprehensive and complete overview on the status of the higher-order calculations in the MSSM and NMSSM in the various approaches has been given in the recent review [37]. In the following, we briefly review the most relevant studies related to the computations and implementations of the higher-order corrections to the Higgs boson masses in the NMSSM with a Z 3 symmetry.
In the CP-conserving NMSSM, the leading one-loop contributions to the Higgs boson masses were presented in [28,[38][39][40][41][42][43][44] while the full one-loop corrections were provided in [45,46] for the DR scheme and in refs. [47,48] for a mixed DR-on-shell (OS) renormalisation scheme. Two-loop corrections at the order O(α t α s + α b α s ) in the DR scheme were obtained in the effective potential approach in [45]. The authors of [49] have provided the two-loop corrections beyond O(α t α s + α b α s ) in the gaugeless limit and in the DR scheme by differentiating numerically or analytically the generic two-loop effective potential presented in [50]. Including CP-violating phases, the leading one-loop corrections were provided in [51][52][53][54][55]. Subsequently, the full one-loop and logarithmically enhanced -1 -

JHEP09(2021)193
tive potential approach in the DR renormalization scheme and have been implemented in SARAH. However, what has been missing in the results is the two-loop OS effects from the electroweak VEV in the transition from the OS VEV to the DR VEV. In our present paper this effect is taken into account by using an OS renormalization scheme for the VEV. As common practice, we use the combination of the gaugeless limit (the weak mixing angle θ W is kept fixed and e → 0 which leads to massless Goldstone bosons) and the vanishing external momentum approximations. On the one hand, these approximations are simple and good in practice. On the other hand they can give rise to the appearance of infrared (IR) divergences which are present in Feynman diagrams with two and more Goldstone bosons in the internal lines. This problem is known as Goldstone Boson Catastrophe (GBC) [92][93][94][95][96][97][98][99]. A practical solution to this problem has been introduced in the computation of the Higgs boson masses for general renormalisable field theories by using a small external momentum expansion in the IR-divergent diagrams [99] and was implemented in SARAH [63]. In this paper, we address the issue of IR divergences in the CP-conserving and CP-violating NMSSM in the Feynman-diagrammatic approach going beyond the limit of vanishing/small external momentum expansion while using an OS scheme for the electroweak VEV, the charged Higgs boson as well as an OS or DR scheme in the top/stop sector and discuss our practical treatment in the code NMSSMCALC.
The paper is organised as follows. Section 2 introduces our notation as well as the NMSSM at tree level. Section 3 discusses the different renormalisation schemes and the derivation of all necessary one-and two-loop counterterms. In section 4 we describe our treatment of the GBC. The set-up of the calculation and the numerical analysis is given in section 5. Section 6 is dedicated to the numerical analysis of the newly calculated contributions. We conclude in section 7.

The NMSSM tree-level spectrum
We work in the Z 3 symmetric NMSSM. For the two-loop corrections of O((α t + α λ + α κ ) 2 ) we use the gaugeless limit and hence follow the same notation as in our previous calculations [59,61]. We describe here only the Higgs, higgsino and top/stop sectors which are relevant for our renormalisation procedure. For the purpose of introducing model parameters, we present here all terms appearing in the NMSSM superpotential with the quark and lepton superfieldsQ,Û ,D,L,Ê, and the Higgs doublet superfieldsĤ d , H u and the singlet superfieldŜ. Charge conjugated fields are denoted by the superscript c.
Color and generation indices have been suppressed. The symplectic product x·y = ij x i y j (i, j = 1, 2) is built with the anti-symmetric tensor 12 = 12 = 1. The parameters λ, κ are complex in general. For simplicity, the Yukawa couplings y x (x = e, d, u) are chosen to be diagonal matrices. This setting has a negligible effect on the Higgs mass calculation since only the Yukawa couplings of the third generation are important. Furthermore, we chose the convention that the y x (x = e, d, u) are real by rephasing the left and right-handed -3 -

JHEP09(2021)193
Weyl-spinor fields as x L,R → x L,R e iϕ L,R . The soft SUSY breaking Lagrangian reads where again the summation over quark and lepton generation indices is implicit. Thẽ Q,ũ R ,d R andL,ẽ R stand for the complex scalar components of the corresponding quark and lepton superfields. The soft SUSY breaking gaugino mass parameters M k (k = 1, 2, 3) of the bino, wino and gluino fieldsB,W l (l = 1, 2, 3) andG as well as the soft SUSY breaking trilinear couplings A x (x = λ, κ, u, d, e) are complex in the CPviolating NMSSM whereas the soft SUSY breaking mass parameters of the scalar fields, m 2 X (X = S, H d , H u ,Q,ũ R ,d R ,L,ẽ R ), are real.

The Higgs boson sector
The tree-level Higgs boson potential inferred from L soft, NMSSM and the F -terms of W NMSSM reads where we neglected the D-terms originating from the gauge sector as they vanish in the gaugeless limit and hence are not needed in the following. Note that gaugeless limit means that the U(1) Y and SU(2) L gauge couplings g 1 → 0 and g 2 → 0 while tan θ W = g 2 /g 1 is kept constant, where θ W is the weak mixing angle. This is equivalent to the limit of vanishing electric charge and tree-level vector boson masses, e, M W , M Z → 0, while keeping tan θ W constant. Expanding the Higgs boson fields around their VEVs v u , v d , and v s , respectively, yields with the CP-violating phases ϕ u,s . The three VEVs can be traded for tan β, the SM VEV v and the effective µ parameter µ eff as JHEP09(2021)193 The MSSM limit can be obtained by taking the limit λ, κ → 0, v s → ∞ while keeping µ eff and κ/λ constant. Using the Higgs potential given in eq. (2.3) we define tree-level tadpoles and mass terms in the broken phase, and the tadpole coefficients where only five of them are independent and t au = t a d /t β . The tadpoles vanish at tree level but affect the higher-order corrections. We keep them, however, for the renormalisation procedure, and set them to zero afterwards. The explicit expressions for the tadpoles and the squared mass matrices M φφ and M h + h − can be found in ref. [59]. The reference also contains a detailed discussion about the two-fold rotation of the where we kept the dependence on the tadpole parameters explicitly for m 2 G 0 for the later discussion of the cancellation of the IR divergences in section 4. It should be noted that M hh is only diagonal for vanishing tadpole parameters. For the charged Higgs fields a single rotation R G − (β c ) is used, where the charged Goldstone boson mass as well as all off-diagonal elements vanish for vanishing tadpoles. The rotation angles β n and β c coincide with β at tree level, β c = β n = β, which has been already applied in eqs. (2.12) and (2.14). We distinguish them since β n and β c as mixing angles do not need to obtain a counterterm while β arising from the ratio of VEVs has to be renormalised and receives a non-vanishing counterterm. After the renormalisation they are set equal to the tree-level value of β again. Note that we denote all masses apart from the charged Higgs boson mass M H ± by small letters m in order to indicate that they are tree-level -5 -

JHEP09(2021)193
masses while we will denote loop-corrected masses by capital M . As discussed later, M H ± will be renormalised on-shell so that there the distinction between tree-level mass and loop-corrected mass does not apply.
In accordance with the SUSY Les Houches Accord (SLHA) [100,101] and for purpose of renormalisation, we decompose the complex parameters A λ and A κ into their imaginary and real parts and λ and κ into their absolute values and phases ϕ λ and ϕ κ , respectively. 4 The phases enter the tree-level Higgs mass matrix in two combinations together with ϕ u and ϕ s , where ϕ y is the only CP-violating phase at tree level in the Higgs sector. In case of vanishing ϕ y , the CP-even components, h u , h d , h s , do not mix with the CP-odd ones, a d , a u , a s . We furthermore trade ImA λ,κ as well as m 2 H u,d ,S for the tadpole parameters t a d ,as and t h d,u,s , respectively, by using the tadpole conditions, cf. ref. [59] for details. In contrast to the MSSM, the trilinear and quartic Higgs couplings in the NMSSM do not vanish in the gaugeless limit, but involve λ, κ, A λ , A κ . The O((α t + α λ + α κ ) 2 ) corrections get contributions from two-loop diagrams containing these Higgs self-couplings. This causes the appearance of the GBC, which will be discussed in detail in section 4.
In NMSSMCALC, we have two possibilities to choose the set of input parameters in the Higgs sector: either In the first choice the charged Higgs mass is an input while ReA λ is an input in the second one. These parameters need to be renormalised at one-and two-loop level.

The squark sector
The relation between the top mass and the top quark Yukawa coupling is given by, (2.19) in which m t and y t are real by our convention. We use the freedom of choice of the phases ϕ L , ϕ R of the left-and right-handed top-quark fields and define ϕ L = −ϕ R = −ϕ u /2. As result, the stop mass matrix in the (t L ,t R ) T basis in the gaugeless limit is given by  where Ut rotates the left-and right-handed stop fieldst L,R into the mass eigenstatest 1,2 . Similar to our previous calculations [59,61], the bottom quark mass is set to zero everywhere, hence the right-handed sbottom states decouple and only left-handed sbottom states are involved in the computation. The parameters in the squark sector that need to be renormalised at one-loop level are

The electroweakino sector
The mass generation for the wino, bino, higgsino and singlino interaction states does not change significantly w.r.t. ref. [59] except that we do not assume λ = κ = 0. Therefore we only shortly repeat the used notation in this section. Since the gauged Weyl-fermions do not couple to any other particles in the gaugeless limit we only consider the 3×3 sub-matrix M χ 0 in the basis (H 0 d ,H 0 u ,S) T for the neutralinos, and the 1×1 matrix for the charginos, where the neutralino masses are ordered as |mχ0 | and U , V denote the 2 × 2 unitary matrices for the rotation from the gauge to the mass basis of the charginos. Note, that we absorbed the phase of µ eff into the chargino mixing matrix V so that the higgsino couplings entering the two-loop diagrams will depend on e iϕµ eff . In contrast to the previous O(α 2 t ) corrections, the singlino and mixed singlino-higgsino states now will also contribute in the two-loop diagrams. The vertex and propagator counterterms involving charginos and neutralinos enter one-loop counterterm inserted diagrams. We therefore need to renormalise them at one-loop level. However, all parameters in the electroweakino sector are also present in the Higgs sector so that we do not need further renormalisation conditions.

Renormalisation of the NMSSM Higgs bosons at the two-loop order
The loop corrected Higgs boson mass spectrum is obtained by iteratively solving 5 det 1 5×5 p 2 − M hh,5×5 +Σ hh (p 2 ) = 0 (3.1) 5 We have confirmed that the contributions from the Goldstone components are numerically negligible.
Thus we drop them in the final calculation.

JHEP09(2021)193
for the squared mass matrix M hh,5×5 = m 2 h i δ ij (i, j = 1, . . . , 5), where the m h i denote the tree-level masses of the tree-level mass eigenstates h i . The numerical recipe for the iterative solution is described in ref. [59]. The Σ hh ij ≡Σ ij stand for the renormalised self-energies for the h i → h j transition and contain the one-and two-loop contributions which are denoted by the superscripts (1) and (2), respectively, The one-loop renormalised Higgs self-energies have already been obtained with full momentum dependence in the CP-conserving and CP-violating NMSSM in refs. [47,57] to which we refer for further details. The two-loop renormalised Higgs self-energies consist of the O(α t α s ) corrections, which we computed in [58], and the O((α t + α λ + α κ ) 2 ) contributions computed in this paper, Note that the O(α t α s ) corrections are evaluated in the approximation of vanishing external momentum. In the O((α t + α λ + α κ ) 2 ) corrections, however, we can choose between including the finite momentum dependence or the Goldstone boson mass as regulator for the IR divergences (see section 4). We therefore keep the momentum dependence in the following formulae. We will drop the superscript (α t + α λ + α κ ) 2 on the self-energies for simplicity of the expressions. The neutral Higgs renormalised self-energies 6 are written as sum of the unrenormalised self-energies Σ ij and the counterterms at one-loop level aŝ (3.4) and at two-loop level In the above formulae, M hh and R are the tree-level Higgs mass matrix and the rotation matrix defined in eq. (2.10) and eq. (2.11). The Higgs mass counterterm matrix at n-loop level is denoted by δ (n) M hh and is obtained by replacing the parameters P i on which it depends by their renormalised quantities plus corresponding counterterms δP i up to nloop order, i.e. P i → P i + δ (1) P i + · · · + δ (n) P i , and expanding accordingly. Its explicit 6 For the inclusion of Goldstone components, i, j take the values 1 to 6 where h6 is identified with G 0 .

JHEP09(2021)193
expression can be found in the appendix G of ref. [59]. The Higgs field renormalisation constant matrix is given by where the renormalisation constants ∆ (n) Z Φ , Φ = H u , H d , S for the doublet and singlet fields will be given in section 3.1.
Similarly for the charged Higgs boson sector, the one-loop and two-loop renormalised self-energies for the transitions h + i → h + j with h + 1 ≡ G + and h + 2 ≡ H + are given by eq. (3.4) and eq. (3.5), respectively, with the following replacements where the charged Higgs mass matrix M h + h − and the rotation matrix R G − are defined in eq. (2.13) and the charged Higgs field renormalisation constant matrix is given by A list of all two-loop diagrams considered in this work to calculate the unrenormalised self-energies for the charged and neutral Higgs bosons is given in appendix C.2.

One-loop and two-loop counterterms
The Feynman integrals are calculated in dimensional regularisation in D = 4 − 2 dimensions. Therefore, intermediate results will contain ultraviolet (UV) divergences of the order O( −2 ) and O( −1 ). To render the renormalised self-energies UV-finite the relevant parameters need to be renormalised to either one-or two-loop order depending on the explicit dependence of the tree-level Higgs boson masses on the parameters. In our case, this means that the top/stop and the electroweakino sector need to be renormalised only at one-loop order while all other parameters are required up to O( −2 ). Note that the chargino and neutralino masses are derived quantities and depend on the (one-loop) counterterms of the input parameters. In the choice of the renormalisation schemes we follow our previous two-loop calculations [59,61]. Note that for parameters defined in the OS scheme, we also study the dependence on O( 1 )-terms in the corresponding one-loop counterterms which are potentially multiplied with O( −1 )-terms from loop-integrals and other one-loop counterterms, thereby generating additional finite contributions. In this section we give explicit expressions for all needed one-loop counterterms after applying our approximations. The remaining one-loop counterterms have already been computed in refs. [47,57] which worked out the renormalisation of the full NMSSM at the one-loop level in the DR, OS and the mixed DR-OS scheme.

The Higgs sector
In the Higgs sector we apply a mixed DR-OS renormalisation scheme. Working in the gaugeless limit at two-loop order the counterterm of the electric charge vanishes. Furthermore, the counterterms of the phases ϕ α (α = s, u, κ, λ) can be set to zero in order to obtain a UV-finite result. Since the charged Higgs mass M H ± can be traded for ReA λ and vice versa we have the following two possible sets of input parameters together with the applied renormalisation conditions, (3.11) in case M 2 H ± is used as independent input, or for ReA λ as independent input. All above listed parameters are renormalised at two-loop level except for the sine s θ W of the Weinberg angle θ W where only the non-vanishing oneloop counterterm contributes. The matrix-valued Higgs field renormalisation constants are needed up to two-loop level and are defined via DR conditions as explained in the following.
Higgs Boson wave-function renormalisation constants. The field renormalisation of the Higgs boson gauge eigenstates 7 (Φ = H u,d , S), with is carried out in the DR scheme. We obtain the counterterms in two equivalent ways. They are computed by either using Feynman diagrams or the renormalisation group equations (RGEs). For the former, they are given by the UV-divergent part of the derivative of the unrenormalised self-energies with respect to the momentum squared For the latter, they can be written as [102,103]

17)
7 All off-diagonal renormalisation constants have been verified to vanish at one-and two-loop order.

JHEP09(2021)193
where γ φφ is the anomalous dimension of the corresponding scalar field φ = h d , h u , h s , x = {y t , λ, κ} with y t = √ 2m t /(vs β ) and β (1) (x) is the one-loop beta function of the coupling x. The functions γ φφ and β(y t ) at one-and two-loop level can be obtained from either [102,103] or the package SARAH. Note that the RGE results are in the pure DR scheme which means that all parameters are renormalised in the DR scheme. In the following we will use the superscript DR on the wave-function renormalisation constants to indicate the pure DR scheme while we use the superscript OS for the scheme where y t and v are renormalised in the OS scheme. Our diagrammatic results in the pure DR scheme are in full agreement with the RGE results. At one-loop order, we find while the two-loop results yield, where k = 1/(4π) 2 . A closer look at eq. (3.19) shows that the scheme change from DR to OS in the top/stop sector introduces additional higher-order contributions via the one-loop field constant δ (1) Z Hu As noted above, the field constants denoted by the superscript OS are actually still DRrenormalised, i.e. only the UV-divergent parts are taken into account, but only refer to the additional UV-divergent sub-loop contributions from the top/stop sector. Note that we write the one-loop OS counterterm in the following form and similarly for all other one-loop OS counterterms. Solving the top mass counterterm eq. (3.67) for δ (1) y t and expanding eq. (3.25) to O( −1 ) we find and and where A 0 and B 0 denote the scalar one-loop one-point and two-point functions and B 1 the tensor one-loop two-point function [104] and log(x) = log x/µ 2 R with the renormalization scale µ R .

JHEP09(2021)193
In the Feynman diagrammatic approach we have computed δ (2) Z Φ in two different ways. In one computation we kept the full momentum dependence in the UV-divergent parts of all diagrams. We then evaluated each contribution with non-zero momentum and found the sum of all contributions being independent of p 2 . In another computation we took the limit p 2 → 0 right after taking the derivative. Here, the coefficients of intermediate results of the single poles feature logarithmic and quadratic IR divergences. As will be discussed later in section 4, a mass regulator can be introduced to deal with the IR divergences. We found full agreement with the finite p 2 -result and no dependence on the mass regulator in the sum of all Feynman diagrams when using the IR-save loop functions defined in appendix A, which gives us yet another possibility to verify if a mass-regularisation scheme is actually useful.
The VEV and the weak mixing angle counterterm. The VEV countertem in the OS scheme is given by where δ nm is the Kronecker delta. The weak mixing angle counterterm at one-loop order reads where M 2 W/Z are the squared vector-boson masses. The vector bosons are renormalised OS with the corresponding counterterms given by denoting the transverse part of the unrenormalised n-loop vector boson self-energy evaluated at zero external momentum. Note that whereas δ (n) M 2 V and M 2 V are separately zero in the gaugeless limit, their ratio entering the counterterms of the VEV and sin θ W is non-zero.
In the pure DR scheme, the one-loop counterterm δ (1) s DR θ W vanishes while the explicit evaluation of the UV-divergent part of the VEV counterterm is found to be This is in accordance with the relation given in refs. [102,103] which connects the counterterm of the VEV v i to the field renormalisation constant of the respective field H i . Exploiting v = v 2 u + v 2 d yields eq. (3.38). Note that since the Σ (2),T Z,W are evaluated at vanishing external momentum, we encounter intermediate IR divergences due to the appearance of massless Goldstone boson propagators. However, these divergences cancel in the sum of all two-loop self-energy diagrams. This will be discussed in detail in section 4.

Tadpole parameters.
Requiring the tree-level minimum of the potential to be the true minimum, higher-order tadpole contributions must be fully compensated by their counterterms, i.e. 8 In eqs. (3.11) and (3.12), we call this in slight abuse of the language an OS condition even though the tadpoles are strictly speaking not associated with any on-shell field. The full set of all two-loop tadpole diagrams considered in this work is given in appendix C.1.

Charged Higgs boson mass.
As mentioned earlier, we have the option to choose either ReA λ in the DR scheme as input parameter or M 2 H ± in the OS scheme. In the approximation of vanishing external momentum, the mixing between the charged Higgs boson and the charged Goldstone boson is negligible. If the external momentum squared is set equal to the charged Higgs boson mass squared one may have to consider this mixing effect, however. We follow our definition of the charged Higgs mass counterterms at one-and two-loop order given in ref. [59]. For convenience, we present here the most important formulae. At one-loop order, the charged Higgs boson mass counterterm in the OS scheme is given by while at two-loop order we have

JHEP09(2021)193
In case ReA λ is given as independent input parameter, the charged Higgs mass counterterms have to be obtained as functions of all other counterterms by inserting their respective loop expansions in the formula for the charged Higgs boson mass. For the explicit formulae of the counterterms and details on the calculation of the loop-corrected charged Higgs boson mass, we refer to [59].
Some of the two-loop charged Higgs boson self-energy diagrams suffer from IR divergences, but the sum of all contributions is indeed IR-finite and does not dependent on the regulator mass.
Ratio of the VEVs tan β. The parameter tan β is given by the ratio of the VEVs v u and v d . Using eq. (3.39) its counterterm can be related to the field renormalisation constants which are calculated in the DR scheme so that the one-loop DR counterterm is given by while the two-loop expansion yields where ∆ (2) Z H u,d were defined in eq. (3.14). As can be inferred from eq. (3.47) the two-loop counterterm of tan β also depends on δ (2) Z Hu which receives additional UV-divergent shifts when we change the renormalisation scheme of the top/stop sector. Accordingly, δ (2) tan β will be affected by such a scheme change.

Superpotential parameters, soft-SUSY-breaking parameters and singlet VEV.
Due to SUSY-non-renormalisation theorems [105][106][107][108][109][110], the superpotential parameters are renormalised through the field renormalisation constants. Therefore, we have two possibilities to construct the one-loop counterterms: (i) using RGEs for the superpotential parameters [49,[73][74][75][76][77], (ii) using the calculated field renormalisation constants together with the non-renormalisation theorems. At one-loop order we verified that the two methods yield consistent results, resulting in the DR counterterms The DR counterterm of the singlet VEV can be obtained by using eq. (3.39), The DR counterterm of the soft SUSY breaking coupling ReA κ can be obtained from either the one-loop RGE or a diagrammatic calculation and reads

JHEP09(2021)193
Likewise, the DR counterterm for δ (1) ReA λ can be extracted from the RGEs, if ReA λ is chosen as input instead of the charged Higgs boson mass, Performing a one-loop counterterm expansion of the following expression 9 (3.53) in the gaugeless limit yields a relation between the DR counterterm of m 2 H ± and ReA λ which was used to cross-check the UV-pole of the charged Higgs boson self-energy.
Furthermore, eqs. (3.51) and (3.53) reveal the implicit dependence of the DR counterterm δ (1) ReA κ on OS defined parameters. This parametrisation of δ (1) ReA κ is useful to study the generation of additional two-loop contributions when performing a scheme change by expanding δ (1) At two-loop order the missing two-loop counterterms are constructed by demanding UV-finiteness in some components of the renormalised neutral self-energies, which were verified to also render all other components UV-finite. We found that the solutions to this system of equations are in agreement with the following expressions of the counterterms for |λ|, |κ| and v S in the pure DR scheme, where eq. (3.57) is in agreement with eq. (3.39) while eqs. (3.55) and (3.56) can also be derived from the NMSSM superpotential, eq. (2.1), using the SUSY-non-renormalisation theorem.
Even though the superpotential parameters are renormalised in the DR scheme, eq. (3.55) changes w.r.t. the single-pole when changing from DR to OS in v and/or the top/stop sector. This is due to the change in δ (2) Z Hu , cf. eq. (3.28). 9 The phases ϕw and ϕy have been defined in eqs. (2.15) and (2.16), respectively.

The higgsino sector
The chargino and neutralino masses are derived quantities that depend on the input parameters v, v S , β, λ and κ, cf. section 2.3. These parameters appear already in the Higgs sector, we therefore do not need additional renormalisation conditions. The chargino mass counterterm is found by expanding the tree-level mass in the gaugeless limit and reads while the neutralino mass counterterm matrix 10 is given by, where the non-vanishing components of the symmetric matrix δ (1) M χ 0 are They only enter in one-loop diagrams with a mass counterterm insertion. As a further crosscheck, we verified that these counterterms render the renormalised chargino/neutralino self-energies in the gaugeless limit UV-finite.

The squark sector
In the squark sector we apply both the OS and the DR renormalisation scheme for the following set of parameters at one-loop level, Expanding the OS and DR counterterms of the parameters X = m t , mQ 3 , mt R , A t in terms of the dimensional regulator ε, we have where we have kept the terms proportional to ε in the expansion of the OS one-loop counterterms, since we want to investigate the dependence of our results on these terms. We follow our definitions of the one-loop counterterms of the squark sector in ref. [59]. We therefore do not repeat the squark counterterm definitions here, but only present analytic JHEP09(2021)193 expressions 11 for these counterterms in the DR scheme. The one-loop counterterm for top quark mass reads where δ (1) v is obtained from eq. (3.38), δ (1) tan β is defined in eq. (3.46) and the counterterm for top quark Yukawa coupling reads The one-loop counterterm for the top-quark trilinear coupling is given by while the counterterms for mQ 3 and mt R read We have found full agreement between the UV-divergent parts of all top/stop counterterms of the OS and the DR scheme.

Independence of O( 1 ) counterterms
At O(α t α s ) [58,61,111] and O(α 2 t ) in the MSSM limit [59] it was shown that finite contributions generated by the O( 1 )-terms of the top/stop counterterms can be compensated by a finite shift in the two-loop wave-function renormalisation constant δ (2) Z Hu .
However, in the present O((α t + α κ + α λ ) 2 ) calculation additional finite contributions are generated in the renormalised self-energies by O( 1 )-terms of on-shell counterterms of v, the charged Higgs boson mass and the tadpoles. In order to keep the parameters tan β, λ, and ReA κ as pure DR parameters, which also corresponds to the aforementioned and applied RGEs, finite parts have to be taken into account in the respective counterterms. These finite parts can be derived from eqs.  11 These results are consistent with the results from the RGEs for yt obtained by SARAH [49,[73][74][75][76][77]].

JHEP09(2021)193
The superscript "OS" on the left-hand side refers to the OS-nature of the implicit parameters v, m t , t i and M 2 H ± . In addition, also the SM VEV counterterm receives a finite contribution, ∂ v δ (1) v| −1 δ (1) v| 1 , from the O( ) expansion at two-loop order. These finite parts in the counterterms cancel the O( 1 ) contributions from the one-loop counterterm insertion diagrams in the self-energies.
Therefore, the implementation of the calculation in NMSSMCALC for simplicity does not contain the O( ) contributions in the OS one-loop counterterms and does not include any finite terms in the two-loop DR counterterms. This leads to the same results for the renormalised self-energies but requires less computational resources.

Two-loop corrections in the gaugeless limit
We are working in the gaugeless limit where one neutral (G 0 ) and two charged Goldstone bosons (G ± ) are massless. In this limit, in contrast to the MSSM, the couplings between two Goldstone bosons with one/two neutral Higgs bosons, two charged Higgs bosons, one/two W and one/two Z bosons do not vanish in the NMSSM. In the computation of the tadpoles and self-energies at vanishing external momentum IR divergences appear due to the massless Goldstone bosons and the non-vanishing couplings. In this section we discuss our approaches to treat these IR divergences in the computation of the tadpoles, charged Higgs boson, W and Z boson self-energies and the neutral Higgs boson self-energies.
Before this study, there was only the code SARAH which implemented the two-loop corrections controlled by the NMSSM superpotential parameters λ and κ [63]. The code makes use of the two-loop effective potential for a general renormalisable theory computed in [50]. A solution of the IR divergences, which are also known as the Goldstone Boson Catastrophe, 12 has been presented in ref. [99]. In the effective potential approach, the tadpoles and Higgs self-energies are derived from the first and second field derivatives of the effective potential V eff H (x run ) where x run denotes MS or DR parameters in the model. The minimum of the effective potential is obtained by solving the tadpole equations order by order. At the minimum of the tree-level potential the running Goldstone boson squaredmass 13 (m run G ) 2 is zero, but is non-zero in general (it can be very small or negative [94]).
with µ R denoting the renormalisation scale. Therefore the tadpole equations contain terms proportional to log(m run G ) 2 which are divergent in the limit (m run G ) 2 → 0 and have an unphysical imaginary part if (m run G ) 2 < 0. The solution 14 proposed in [99] at two-loop 12 Originally, the term Goldstone Boson Catastrophe goes back to spurious IR divergences and imaginary parts encountered in the SM effective potential in the Landau gauge and its first derivative [112]. 13 For simplification, we do not distinguish between the masses of the neutral, positively and negatively charged Goldstone bosons here. In general, their running masses can be different. 14 It is close to the solutions worked out in refs. [94][95][96]113] for the SM and extended for the MSSM in ref. [97]. The Goldstone contributions are resummed by integrating out all heavy degrees of freedom when calculating the n − 1 loop-corrected Goldstone boson mass and using this effective Goldstone boson mass in the minimization of the n-loop effective potential. It was shown in ref. [98], that this procedure resums the spurious IR divergences to all orders in perturbation theory.
where the OS mass is set to be zero and Σ (1) GG denotes the one-loop unrenormalised selfenergy for the Goldstone boson component. Using the above relation directly in the expressions of the tadpoles and self-energies given in [114] a cancellation is found of the divergent terms log(m OS G ) 2 in the tadpoles at two-loop order. However, the divergences remain in some sets of neutral Higgs boson self-energy diagrams. These diagrams are then identified and calculated using a small momentum expansion.
In contrast to the previous study, we use the Feynman diagrammatic approach to directly calculate scalar one-and two-point functions at the two-loop order and do not use the available general expressions for the self-energies. Our intermediate results contain UV and IR divergences. The UV divergences are canceled by the counterterms of the parameters and fields introduced via the renormalisation procedure as discussed in the previous section. We are flexible in our choice of the renormalisation schemes for different parameters. In particular, we renormalise the tadpoles in the OS scheme. This means that the masses of all Goldstone bosons are always on-shell and zero at all orders in the perturbation theory. Using this OS scheme we find the full cancellation of the IR divergences in the tadpoles, charged Higgs, W and Z boson self-energies, but only a partial cancellation in the neutral Higgs boson self-energies. For the sets of diagrams that contain Goldstone bosons but do not have IR divergences we set the mass of Goldstone boson equal to zero and do not mention them furthermore in this section. We now discuss the sets of one-and/or two-loop diagrams where IR divergences appear. In appendix B, tables 5 and 6, we give the complete list of all IR-divergent two-loop tadpole and self-energy topologies together with the loop integrals causing the IR divergences. In the last column of these tables we specify whether the IR divergence is spurious, i.e. cancels against contributions of other diagrams, or whether a non-zero external momentum needs to be included for the regularisation.
We regulate all IR divergences in three different ways: (i) By using a mass regulator everywhere. This allows us to study which IR divergences are actually spurious. We expect a residual dependence on the mass regulator only in the subset of diagrams that require momentum regularisation for a physical result but which can in principle also be treated by mass regularisation, see section 4.1. (ii) By using a mass regulator only in the subset which features spurious IR divergences. In the remaining subset with the genuine IR divergences, we use non-zero external momentum and apply analytically known results for the small momentum expansion of the loop integrals. This is equivalent to the generalised effective potential approximation introduced in refs. [63,99] and is described in detail in section 4.2. (iii) By including the full external momentum dependence in the computation of all Feynman diagrams of O((α t + α λ + α κ ) 2 ) making use of TSIL [115]. This requires only the regularisation of the spurious IR divergences, see section 4.3. . Therefore, some of the involved two-loop integrals are identical to those in the expansion in the small Goldstone boson mass used in the resummation procedure of the generalised effective potential [97,99]. However, since we also want to regulate those divergences that do not cancel out in the sum of all diagrams, we require a larger set of expanded loop functions which is given in appendix A. Performing the expansion instead of simply setting the Goldstone mass equal to a finite mass value should reduce the dependence on the regulator mass further.

Infrared mass regulator
In order to investigate the cancellation of the IR divergences we divide the topologies that contain IR divergences into five sets, cf. tables 5 and 6 in the appendix. Set A includes topologies without external momentum flowing into the loop, i.e. the tadpole topologies 1-3 in table 5 and the self-energy topologies 8, 10 and 13 in table 6. It is evident that these must form IR-finite subsets. Set B contains the self-energy topologies 4, 7, 11 with two Goldstone bosons with the same momentum in table 6. This set B is also IR-finite. The self-energy topologies 4, 7 and 11 with two Goldstone bosons that couple with one external line belong to set C while the self-energy topologies with three Goldstone bosons belong to set E. Set D contains the self-energy topologies 5, 6, 9 and 12. Tadpoles contain only topologies of set A, while the charged Higgs self-energy contains topologies of the sets A and B. The W and Z boson self-energies contain topologies of the sets A, B, D, and E where each of these sets is separately IR-finite. The cancellation happens due to the contributions from the finite parts of the neutral and charged Goldstone boson mass counterterms δ (1) m 2 G 0 and δ (1) m 2 G ± which are related to the finite part of the one-loop counterterms as, The neutral Higgs boson self-energies contain all five sets of diagrams. The first derivative of the UV parts of the self-energies with respect to the momentum squared are used to compute the wave function renormalisation constants (WFRs) δ (n) Z Φ . We find that the WFRs are IR finite. This is consistent with [116,117] where it was shown that the WFRs do not depend on the regulator mass as long as appropriate counterterms are introduced, cf. [118] for a recent application. For the neutral Higgs boson self-energies we do not see such a cancellation in the sets D, and E but only in the sets A, B and C. Individual contributions to the sets D and E may contain divergences proportional to 1/M 2 R , however in the sum they are canceled out. The remaining IR divergence depends on the mass regulator as logM 2 R and log 2 M 2 R . A detailed study of the phenomenological impact of the residual M 2 R dependence on the Higgs boson mass corrections is performed in section 6.3.

Partial momentum dependence
A solution of the GBC not involving a mass regulator and without the need of the time consuming numerical evaluation of all loop integrals including external momentum is given -21 -

JHEP09(2021)193
by the generalised effective potential approximation [99]. Here only the subset of diagrams of figure 12 with a residual dependence on M 2 R is computed at p 2 = 0. Additionally, the number of independent mass scales is significantly reduced for this particular subset making it possible to evaluate it analytically either with exact p 2 dependence or approximately around p 2 = 0.
As mentioned in the previous section, there exist residual IR divergences in the sets D and E of the neutral Higgs boson self-energies. These IR divergences can be avoided by using non-zero external momentum. The obvious and fastest way is to use a small momentum expansion only in these sets which is similar to the implementation of the effective potential approximation in SARAH [63]. For the sake of comparison we also apply this method in our code. For sets D and E we calculate Feynman diagrams with full momentum dependence and use the expansion of the loop integral around p 2 = 0. In this evaluation the masses of the Goldstone bosons are set to zero. In the expansion we removed terms of O(p 2 ) but kept terms proportional to 1/p 2 , log(p 2 ) and terms independent of p 2 . The necessary special cases for the loop integrals have already been worked out in refs. [63,99,115,[119][120][121][122].
Even when using finite external momentum, the individual diagrams in set C, D and E still feature an IR divergence O(logM 2 R ) originating from the integrals V(x, 0, z, u) and C(x, 0, 0) which cancels in the sum of all contributions. The cancellation can be obtained by making use of the identity The integralṼ (0, z, u) is IR-finite for p 2 = 0 (since it scales with log p 2 ) and has been calculated in ref. [63] 15,16 Therefore, the choice of regulating C(x, 0, 0) is not important as there is no dependence on this function in the final result. The small momentum approximation, however, breaks down latest near the various thresholds involved in the diagrams, making the full momentum dependence necessary for a reliable result. At one-loop order, we find that a strict expansion of the massless scalar two-point integral around p 2 = 0 gives the same dependence on the momentum-regulator, B(0, 0)| p 2 ≈0 ∝ − ln p 2 , as the dependence on the mass regulator when starting with p 2 = 0 and expanding around a small Goldstone boson mass This might raise the question whether the two expansions are also connected at the two-loop level. However, we only expand the IR-divergent two-loop integrals around p 2 = 0 but use the exact analytic result for the one-loop integral (without expanding it). This leads to additional constant and log p 2 terms which are also present in the full-momentum calculation.

Full momentum dependence
In this approach the full momentum dependence is taken into account in all two-loop diagrams of the O((α t + α λ + α κ ) 2 ) corrections in the gaugeless limit. This has not been 15 Note that our notation slightly differs from appendix A.1.1 in ref. [63]. Their B(x, y ) and P(z, u) corresponds to our C(x, y, y) and − B(z, u)| p 2 =0 . 16 There is a sign mistake in the corresponding identity in ref. [63] eq. (A.5). However, the implementation in SARAH has the correct sign.

JHEP09(2021)193
done in the literature for the NMSSM before our study. The momentum dependence for the O(α t α s ) corrections has been studied in refs. [66,111,123,124] within the MSSM using differential equations and sector decomposition for the numerical evaluation of the loop integrals. We expect that the momentum dependence in the NMSSM is comparable to the one of the MSSM O(α t α s ) corrections which were found to be at most about one GeV for the loop-corrected SM-like Higgs boson mass compared to the zero momentum approximation.
Applying the integral basis used by TSIL the generalisation of our framework to nonzero external momentum is straightforward. There is only one class of diagrams, figure 11 (h), which requires more care as it contains additional 1/p 2 -terms that do not allow to numerically take the limit p 2 → 0. Therefore, we set p 2 = 0 before invoking TARCER in this particular diagram when using the zero/partial-momentum approximation and assume arbitrary p 2 when reducing the integral for the full-momentum calculation. Using p 2 = 0 in all two-loop diagrams requires the inclusion of additional wave-function renormalisation constants in eq. (3.5). We have checked that this indeed restores UV-finiteness in the full momentum approach. Further modifications to the two-loop counterterms are not necessary. However, the use of an OS scheme for the charged Higgs boson mass would in principle also generate additional finite shifts originating from the momentum dependence of the charged Higgs boson self-energy. For simplicity, we do not include momentum dependence in the calculation of the OS charged Higgs mass counterterm and the VEV counterterm.
The Higgs boson masses are obtained by iteratively solving for the pole of the (twoloop) propagator, eq. (3.1). However, so far the p 2 -dependence was only taken into account at the one-loop level in NMSSMCALC such that the result of the two-loop self-energies had to be calculated only once. The inclusion of external momentum in the new two-loop self-energy corrections would require them to take part in the iterative procedure and thus slow down the overall runtime by several orders of magnitude. Therefore, we chose a fixed value of p 2 = (m 2 h i + m 2 h j )/2, where m h i,j are the tree-level Higgs boson masses, for the calculation of the new self-energy correctionsΣ (2) ij (p 2 ). This approximation would in principle require a more detailed study to estimate its numerical impact. However, in ref. [63] the same situation was studied for the SM with the result that the Higgs boson mass prediction varies only by a few MeV if the external momentum is varied by several orders of magnitude.
Note that results for the NMSSM Higgs mass corrections including the momentumdependent contributions while at the same time applying the gaugeless limit should be taken with care, cf. also the discussion in [37]. While in the MSSM the momentum-dependent and the electroweak gauge contributions are of similar size, in the NMSSM there are additional F -term contributions of O(λ, κ) which could yield an additional enhancement of the momentum-dependent corrections if λ, κ ∝ O(1), especially through the mixing of heavy and light Higgs components in the self-energies. Therefore, the new momentum-dependent results should be taken with care as their contributions could be either comparable or sub-dominant to the missing electroweak gauge coupling contributions depending on the considered parameter point. We leave the inclusion of the gauge-dependent contributions for future work.

Tools, checks and NMSSMCALC release
We performed two independent calculations to derive the here presented new two-loop corrections to the NMSSM Higgs boson masses and cross-checked the results against each other. Both of them used SARAH 4.14.3 [49,[73][74][75][76][77] to generate the model file including the vertex counterterms. This file was used in FeynArts 3.1 [125,126] to generate all required one-and two-loop Feynman diagrams for the calculation of the mass corrections. We used FeynCalc 9.2.0 [127,128] for the evaluation of the fermion traces and the tensor reduction of the one-and two-loop integrals and the amplitudes with the counterterminserted diagrams. For the reduction to the two-loop master integrals including the full momentum dependence we additionally used TARCER 2.0 [129], a patched version that comes with FeynCalc. We use the loop integrals defined in TSIL [115]. They are the basis integrals of TARCER extended by a few convenient functions.
The two implementations differ in the way the tensor reduction is performed, namely: (i) setting p 2 = 0 before the reduction and writing the result in terms of one-and twoloop tadpole integrals and (ii) using general p 2 -dependence during the whole procedure. While method (i) is only able to regulate the GBC using a mass regulator as described in section 4.1, method (ii) is more flexible and also able to include partial as well as full external momentum dependence as described in sections 4.2 and 4.3. In addition to the consistency checks regarding UV-finiteness discussed in the previous sections, the two implementations have been cross-checked against each other in the limit p 2 = 0.
In section 6 we also investigate the p 2 dependence of our results. The computation at non-zero p 2 , however, significantly increases the complexity and runtime of the code due to the dependence on the external library TSIL. Even though TSIL was specifically designed for the evaluation of two-loop self-energy integrals, the runtime of one parameter point with a naive Fortran implementation can be of O(hours-days) because of the large amount of different mass scales and diagrams entering at the considered order. Therefore, the second implementation is not part of the public NMSSMCALC release but only consists of private Mathematica notebooks that make heavy use of caching and parallelisation in order to speed up the computation. Non-zero p 2 results can be provided on request. Taking the computation with full momentum dependence as a reference result, we will perform comparisons with the results obtained when using a mass regulator by choosing different sample values of the regulator mass. We observed that the phenomenological impact of using a mass regulator rather than partial external momentum is negligible compared to the overall size of the new two-loop corrections in the considered parameter space for the specific value of the regulator mass given by M 2 R = 10 −3 µ 2 0 . Therefore, the use of a mass regulator is considered to be a good compromise between accuracy and fast numerical performance.
The updated version of NMSSMCALC including the new two-loop corrections to the NMSSM Higgs boson masses in the CP-conserving and CP-violating NMSSM can be downloaded from the url: https://www.itp.kit.edu/~maggie/NMSSMCALC/ -24 -

JHEP09(2021)193
On this webpage we give a description of the new files that have been included. In the input file inp.dat, besides the option to choose the computation of the new two-loop corrections in the block MODSEL we added a new block REGFACTOR that allows to choose the size of the regulator mass which by default is set to M 2 R = 10 −3 µ 2 0 (cf. the discussion in section 6.3). Note, that our new two-loop computation becomes numerically unstable for nearly degenerate Higgs mass values. Therefore, in NMSSMCALC we automatically switch to the computation of the O(α t (α s + α t )) corrections for |m H i − m H i±1 | ≤ 10 −3 GeV where the m H i,i±1 denote the tree-level Higgs masses in the gaugeless limit. This is also done in case the user has chosen in the input file to compute the O((α t + α λ + α κ ) 2 + α t α s )) corrections. In the output file the actually computed loop-order will be stated.

The parameter scan
For the numerical discussion of our results we performed a scan in the NMSSM parameter space in order to obtain parameter scenarios that are in accordance with the most recent experimental constraints. We checked the parameter points of our random scan against compatibility with experimental constraints from the Higgs data by using HiggsBounds 5.9.0 [130][131][132] and HiggsSignals 2.6.1 [133]. The required effective NMSSM Higgs boson couplings normalised to the corresponding SM values have been obtained with the Fortran code NMSSMCALC [83]. One of the neutral CP-even Higgs bosons, called h from now on, is required to behave as the SM-like Higgs boson and have a mass in the range We follow the SLHA format [100] in which the soft SUSY breaking masses and trilinear couplings are understood as DR parameters at the scale This is also the renormalisation scale that we use in the computation of the higher-order corrections. In table 1 we summarize the ranges applied in the parameter scan. In order to roughly ensure perturbativity below the GUT scale we require that both λ and κ remain below 0.7. According to the SLHA format, also λ, κ, µ eff and tan β are understood to be DR parameters at the scale M SUSY . Note that in the scan we kept all CP-violating phases equal to zero. For the investigation of the impact of CP violation we will then turn on individual phases. We retain scan points with a χ 2 computed by HiggsSignals-2.6.1 that is consistent with an SM χ 2 within 2σ. 17 We omit parameter points with any of the following mass configurations, The first constraint ensures that no large logarithms appear that would jeopardize the validity of a fixed-order calculation. The second condition (ii) excludes degenerate mass configurations for which the two-loop part of the NMSSMCALC code is not yet optimised. 18 The third condition takes into account model-independent lower limits for the lightest chargino and stop masses.

Results
In the subsequent numerical analysis we show scatter plots summarising the overall behaviour of our new results and perform specific investigations for two sample parameter points. The first point, P1OS, has been chosen among our allowed parameter points and is defined as follows: Parameter point P1OS. All complex phases are set to zero and the remaining input parameters are given by In accordance with the SLHA format µ eff is taken as input parameter, from which v s and ϕ s can be computed using eq. (2.7). We call this point P1OS in order to mark that the SM-like Higgs boson mass value around 125 GeV is obtained for the OS renormalisation in the top/stop sector. Table 2 summarizes the mass values that we obtain for P1OS at tree level, at one-loop order and at two-loop level including the previously available O(α t α s ) and O(α t (α s + α t )) corrections and finally the corrections including our new results, the O((α t + α λ + α κ ) 2 + α t α s ) corrections. From now on, we denote these by α 2 new , i.e. the size of the loop corrections, in the following plots we label the Higgs bosons according to their dominant admixture and not by their mass ordering unless stated otherwise. 19 In this way we make sure to consistently compare and interpret the impact of the loop corrections.
Defining the absolute value of the relative correction to the mass value when going successively from loop order a in the tables to loop order b in the next row by |m b −m a |/m a , we see that the h u -like mass changes considerably upon inclusion of the one-loop corrections, by 53% in the OS and by 31% in the DR scheme. Since in the DR scheme we already partly resum higher-order corrections, the relative O(α t α s ) correction compared to the one-loop result is only 4% while in the OS scheme we have a reduction by 11%, moving the obtained mass values in the two schemes close to each other. The inclusion of the O(α 2 t ) correction in addition increases the discrepancy again, as already observed in our publication [59]. We have a relative correction of 5% in the OS and of close to 0% in the DR scheme. Our newly calculated corrections move the two values a little bit closer again. The DR result is increased by a very small amount and the OS value is slightly reduced. Still the absolute difference between the two results amounts to about 5 GeV. As expected the overall size of the two-loop corrections is much smaller than the one-loop corrections and amounts to a few percent. This behaviour is also reflected in figure 1 that we will discuss in the next section.
While the parameter point P1OS is characterized by a small singlet admixture to the h u -like Higgs mass we also present results for a parameter point P2OS which features large singlet admixture to the h u -like mass in order to investigate the impact of our newly computed corrections. It is defined by   is always h u -like (h s -like). The partial resummation of higher-order corrections through DR renormalisation in the top/stop sector implies smaller one-and two-loop corrections as compared to the OS scheme and thereby less sensitivity to possible singlet admixture effects in the mass corrections. The impact of our newly computed corrections with respect to the already available two-loop corrections at O(α t (α t + α s )) is less than 1% for the two considered parameter points. However, we will show in the following subsections that the corrections can be enhanced for large values of λ and κ.

Impact of the new two-loop corrections
In the following we discuss in more detail the impact of our newly computed two-loop corrections, both for the point P1OS with small singlet admixture and for the point P2OS with large singlet admixture.
with α 2 i = α t (α s + α t ) and α 2 i = α t α s , respectively. Here and in the following, loopcorrected Higgs mass values are always denoted by capital M . We find that in the DR scheme the relative impact of our new corrections with respect to the previous two-loop orders is about the same and varies between 0 and 7% for λ = 0 to 1.5. In the OS scheme the relative impact with respect to the O(α t (α s + α t )) corrections is in the same range as in the DR scheme while it varies between 6% and more than 10% for ∆  figure 1 (right) we furthermore infer that the corrections are asymmetric with respect to the sign of A OS t . The lower panels in figure 1 show the relative change in the mass corrections at fixed loop order when switching the renormalisation scheme in the top/stop sector, (6.5) The comparison of the results in the two different renormalisation schemes gives one ingredient for the estimate of the uncertainty on the Higgs mass values due to missing higher-order corrections. In the whole plotted λ range at O(α t α s ) the impact is less than 1%, while it increases to values between about 3 and more than 5% upon inclusion of the O(α t (α s +α t )) and slightly less in the new O(α 2 new ) corrections, respectively. Also for the plotted A OS t values the renormalisation scheme dependence is larger for these loop orders (with values between 2.5 and more than 5%) than for O(α t α s ) except for large negative A OS t values. In general one needs to be careful in drawing conclusions on the remaining theoretical uncertainty at a given loop order as long as not all existing contributions at the investigated loop order are included. Since the scheme dependence induced by the top/stop sector is also not significantly reduced upon inclusion of the O(α t (α λ + α κ )) corrections, the 3-loop corrections of O(α t α 2 s ) or even beyond might be required. However, in the rest of this section we show that there are cases with an h u -like Higgs boson at 125 GeV, where the scheme dependence is significantly reduced.
Large singlet admixture. We now turn to the impact of our corrections for the benchmark point P2OS which is characterized by a large singlet admixture to the h u -like Higgs state. In figure 2 we show the absolute mass values as a function of λ (upper) and the dependence on the renormalisation scheme (lower) for P2OS. The notation is the same as in figure 1 (left). Like in the case with small singlet admixture, all two-loop corrections are close to each other in the DR scheme for λ ≤ 1, for λ ≥ 1 the new corrections start to deviate from the previous ones reaching a relative correction ∆ In the OS scheme the impact is slightly more pronounced. We find non-zero relative corrections ∆ α 2 new αt(αs+αt) in the OS scheme starting for λ ≥ 0.5 increasing to up to 5% for λ = 2. The impact of the new corrections with respect to O(α t α s ), however, varies from 6% at λ = 0 to zero for λ around 1.2 and up to 5% again at λ = 2. As for the renormalisation scheme dependence, cf.   is reduced by 1-2% w.r.t. the O(α t (α s + α t )) result. This shows that an estimate of the residual uncertainty due to missing higher-order corrections based on the renormalisation scheme variation in the top/stop sector cannot be made without taking into account the computation of the complete two-loop corrections since the top/stop sector contributes with mixed contributions such as O(α t (α λ + α κ )).  lines signals the following. Transparent lines denote the A OS t range where the experimental Higgs signal constraints are not fulfilled any more. Full lines indicate the parameter region where the h u -like Higgs boson is the lightest Higgs state in the spectrum, dashed lines correspond to the second-lightest Higgs boson being h u -like. We have already seen in table 3 that the nature of the Higgs mass eigenstate can change depending on the loop order. We also see dips in the plots. Here the singlet-doublet admixture of the two lightest Higgs states becomes large inducing nearly same values for their respective coupling values to the SM particles. The comparison of the coupling values for the various two-loop orders (comparison of the red, blue and black lines) clearly shows that the Higgs couplings and hence the Higgs boson phenomenology is strongly affected by the order of included loop corrections. The comparison of the full and transparent regions shows how the allowed parameter range is impacted by the loop order. This underlines the need of precision calculations in order to be able to delineate the underlying parameter range through the measurement of the Higgs properties.
Whole parameter sample. In figures 4 and 5 we investigate the overall impact of our corrections by looking at the whole parameter sample. In figure 4 we show for all allowed parameter points obtained in our scan the relative sizes ∆ for h s (a s ) (not shown here) we find for most cases, when going from O(α t (α s + α t )) to O(α 2 new ), a relative increase in the mass values below 0.4% (0.1%) with some outliers up to 6% (3%).
In figure 5 (left) we show the relative change ∆ ren , see eq. (6.5), of our O(α 2 new ) corrections to the h u like mass M hu for all allowed parameter points as a function of the NMSSM-specific combination √ λ 2 + κ 2 where the color code indicates the value of |A OS t |. The change of the renormalisation scheme requires a conversion of the involved top/stop sector parameters, so that ∆ ren clearly depends on the value of |A OS t |. This is also observed in the plot, where the largest effects from the change of the renormalisation scheme in the top/stop sector are found for large |A OS t | values. The smallest renormalisation scheme dependence is obtained for small |A OS t | and large √ λ 2 + κ 2 . In this regime the impact of the new corrections from the Higgs-and electroweakino sectors becomes more pronounced while the top/stop sector contributes less. Therefore, the renormalisation scheme dependence introduced by the top/stop sector is reduced further. Overall the renormalisation scheme dependence at O(α 2 new ) is larger than 5% if |A OS t | > 2 TeV. This is not surprising since we only consider points with mt 2 < 2 TeV in the scan. Corrections of the order O(A t /mt i ) (and higher powers) become large and introduce a large scheme dependence if |A t | > mt i .  However, for |A OS t | 1 − 2 TeV the scheme dependence is under good control (i.e. smaller than the overall size of the new two-loop corrections) in the NMSSM-specific parameter region √ λ 2 + κ 2 > 0.6. The right plot of figure 5 shows the renormalisation scheme dependence ∆ ren at O(α 2 new ) for the singlet-like scalar h s (green) and pseudoscalar a s (blue) for all of the considered scenarios where they are lighter than h u . A tendency of increasing maximum values for the renormalisation scheme dependence with rising NMSSM-specific couplings can be inferred from the plot. Since the states are mostly h s -like, one would assume that they are less affected by higher-order corrections involving α t in the new corrections. However, we find renormalisation scheme dependences of up to about 5% for very large λ and κ. This shows that corrections of O(α t (α λ + α κ )) can indeed be important for singlet-like states if λ and κ become large. Overall the renormalisation scheme dependence is of typical sizes expected at two-loop order.

Renormalisation scale dependence
We now turn to the discussion of the renormalisation scale dependence of the two-loop corrected Higgs boson mass which can be taken as a rough estimate of the uncertainty due to missing higher-order corrections. The RGEs implemented in NMSSMCALC are used as described in appendix E of [59]. We define the scale dependence of the loop-corrected Higgs mass M h at a given loop order as the relative change of the mass value at the scale µ with respect to our default scale µ 0 = M SUSY , hence (blue) with DR renormalisation in the top/stop sector as a function of a variation of the renormalisation scale µ between one half and twice the default scale µ 0 . The lower plots show the relative scale dependence ∆ scale . For both points the scale dependence is rather small and remains below 3% for all three two-loop corrections. For P1OS we observe the smallest dependence for O(α t α s ). It increases after including the O(α t (α s +α t )) corrections and is reduced again with our new corrections O(α 2 new ). For P2OS we find a similar relation between the different contributions. It has to be noted here that the renormalisation group equations applied in the generation of the plot include all two-loop contributions while for consistency at the three different two-loop orders only the respective contribution corresponding to the given included loop corrections should be taken into account in the renormalisation group equations so that the comparison of the three curves should be taken with caution. Still, overall we see that the inclusion of the two-loop corrections leads to rather small remaining scale dependences. However, the inclusion of the new corrections from the Higgs-and electroweakino sector does only lead to a minor reduction of the scale dependence. Therefore, we argue that the largest uncertainty comes from the top/stop sector requiring either higher orders or the resummation of large logarithms.

Numerical comparison of the three regulation schemes
In section 4 we discussed three regulation schemes to cure the Goldstone boson catastrophe, which we compare in this section for the parameter point P2OS. with the results of the mass computation where the IR divergence is regulated by a mass regulator, cf. section 4.1.
Here, we consider four cases for the regulator mass M 2 R , namely R = M 2 R /µ 2 0 = 10 −5 , 10 −3 , 1 and 10 3 . Note that we included the full momentum dependence in all diagrams of the considered order whereas the partial momentum approximation and the regulator mass is only applied in the GBC subset which is always proportional to λ and κ. As expected partial momentum inclusion approximates the full momentum dependence best, as well as  the regulator mass result for R = 10 −3 . For regulator masses departing more and more from 10 −3 µ 2 0 the difference increases. The difference for all approximations shows a strong dependence on λ but agrees for λ = 0 as expected since the IR effects are related to λ. The largest deviation from M full−p 2 h 1 amounts to ∼ 1.5 GeV in the mass regulated scheme with R = 10 3 . In contrast, the maximum deviation remains below 150 MeV for R = 10 −3 for this parameter point and λ varying between 0 and 2. We note that R = 10 −3 corresponds to a regulator mass of O(30 − 60 GeV) for typical values of the renormalization scale of O(1−2 TeV). This is in good agreement with the findings of [63] which have best agreement between their mass-and momentum-regulated results for a mass-regulator of O(100 GeV). The kinks in the plot are not an artifact of the numerical integration but are due to the strong dependence of the neutral Higgs mixing matrix on λ and the resulting change of certain sub-dominant admixtures when departing from the region allowed by the current collider constraints (grey shaded region).
The result that the full-and partial-momentum corrections do not deviate by more than 1-200 MeV shows that our incomplete calculation of the full-momentum corrections has to be taken with care as these corrections are of similar size as the momentum corrections of O(α s α t ) one might expect comparing to results from the MSSM [66,111,123,124] (which are not calculated in this work) even for very large values of λ. Therefore, the calculation of the momentum-dependent corrections at the order O(α s α t ) is still an open task to be done.
In order to get a more general picture we compare the results in the various approximations for all of our allowed points. Figure 8 (left) shows the relative difference ∆ partial−p 2 R=10 −3 (M hu ) between the O(α 2 new )-corrected M hu value obtained in the partial momentum approximation and in the mass-regulated computation with a regulator mass squared of 10 −3 µ 2 0 for all allowed parameter points. Note that these two different treatments only af--37 -

JHEP09(2021)193
fect the IR-divergent diagrams. The results are shown as a function of the NMSSM-specific parameter √ λ 2 + κ 2 and the color code indicates tan β. The maximum relative difference increases with √ λ 2 + κ 2 and saturates at values close to 1 permille. This behaviour is expected, since the IR-regulated diagrams -and therefore the dependence on the IR regulator mass -are always proportional to λ and κ. The right plot displays the relative difference in M hu when including full momentum dependence in the O((α t + α λ + α κ ) 2 ) contributions and the mass-regulated result for R = 10 −3 when considering a subset of one thousand random points. 20 Note that the full momentum dependence is also taken into account in all IR-finite diagrams of the order O((α t + α λ + α κ ) 2 ). The maximum deviation shows only a weak dependence on √ λ 2 + κ 2 and reaches at most 2 permille for large values of √ λ 2 + κ 2 . The behaviour shows that the full-momentum corrections at O(α 2 t ) (i.e. for λ, κ → 0 in the right plot of figure 8) and those at O((α t + α λ + α κ ) 2 ) for finite λ are equally important. For the singlet-like masses of our parameter sample we found similar results. This allows us to conclude that the regulation of the GBC with a regulator mass of R = 10 −3 is a good compromise between accuracy and computational costs 21 with a difference that remains in the subpercentage range. Also R is not too small in this case to lead to numerically instable results. Comparing figure 8 (left) with figure 4 (left), we conclude that the error made for R = 10 −3 compared to the result obtained with partial momentum dependence is always around one order of magnitude smaller than the overall size of the new two-loop corrections. Therefore, all parameter samples that we present in this paper have been obtained with a regulator mass of R = 10 −3 , unless stated otherwise. The value R = 10 −3 is also the default setting that we have implemented in our new version of NMSSMCALC that has been extended to the here presented new two-loop corrections.

CP-violating phases
In this section, we discuss the influence of the CP-violating phases on the loop corrections to the Higgs masses. In figure 9 (upper) we show by choosing the parameter point P2OS as starting point the loop-corrected mass M hu of the h u -like Higgs boson 22 as a function of the CP-violating phases ϕ λ (full) and ϕ At (dashed) at O(α 2 new ) (red), O(α t (α s + α t ) (black), and O(α t α s ) (blue) for OS (left) and DR (right) renormalisation in the top/stop sector. In the phase variation only one phase is varied at a time. The phase ϕ λ is varied such that the CP-violating phase ϕ y (eq. (2.15)) appearing already at tree level is kept zero, more specifically ϕ λ = 2ϕ s = 2/3ϕ µ eff. and ϕ κ = ϕ u = 0. We thereby ensure to study only radiatively induced CP-violating effects. Otherwise all CP-violating phases are kept zero. Note also that for illustrative reasons, the phases are varied in ranges beyond their allowed validity by the EDM constraints. 23 The lower inserts quantify the effect ∆ α 2 new α 2 i of the newly calculated loop corrections with respect to the O(α t (α s + α t )) (black) and 20 This was done to save computational resources. We confirmed that this subset still contains all important features of the original sample. 21 The computation time considerably increases when the full momentum dependence is included. 22 The kinks in the plot appear at parameter configurations where the hu-like and hs-like Higgs bosons swap their mass ordering and are therefore nearly degenerate. 23 Actually, non-zero ϕ λ are excluded by the EDM constraints for these parameter points while ϕA t is allowed in the whole range.  the O(α t α s ) (blue) corrections. Both in OS and DR renormalisation all three two-loop corrections show the same behaviour with respect to a variation of ϕ At . This can also be inferred from the lower panels where the dashed lines are almost flat. The variation of the loop corrections with a change of the phases is more pronounced in the OS scheme than in the DR scheme, also the behaviour of the new corrections differs more from the other two two-loop orders. Still the curves for the relative deviation ∆

Conclusions and outlook
We have computed the O((α t + α λ + α κ ) 2 ) two-loop corrections to the Higgs boson masses of the CP-violating NMSSM in the Feynman-diagrammatic approach in the gaugeless limit and at vanishing external momentum. While these limits give a good approximation and simplify computations significantly, they may induce infrared divergences in two-loop Feynman diagrams with multiple massless Goldstone bosons. We have shown that using OS conditions for tadpoles at one-loop order makes the two-loop tadpoles, the charged Higgs, the W and the Z boson self-energies at vanishing external momentum at O((α t +α λ +α κ ) 2 ) IR finite, however, IR divergences remain in the neutral Higgs boson self-energies. For the treatment of the IR divergences we have followed three different approaches -the introduction of a regulator mass, the application of a small momentum expansion, and the inclusion of the full momentum dependence in all Feynman diagrams of O((α t +α λ +α κ ) 2 ). By comparing the three methods, we found that the regulator mass approach reproduces the momentum-dependent results well for squared regulator masses that amount to a permille of the renormalisation scale squared. Due to the robustness of this approach we have implemented the new corrections using this value for the regulator mass as default in the published version of NMSSMCALC. In order to quantify the impact of our newly computed corrections we have performed a scan in the NMSSM parameter range and kept only those points that are compatible with experimental constraints. We found that our corrections increase with λ and κ as expected. For λ and κ values compatible with perturbativity below the GUT scale, the corrections are less than 3% relative to the already available O(α t (α t + α s )) corrections. Our new corrections reduce slightly the theoretical uncertainties due to missing higher-order corrections that we estimated by changing the renormalization scheme in the top/stop sector and by varying the renormalization scale. We have also shown that the impact of the new corrections on the Higgs mixings, which manifest themselves in the couplings between the Higgs bosons and the SM particles, is significant and strongly affects the compatibility with the Higgs data. The impact of the CP-violating phases on the new corrections has been found to be small. With our calculation we further improve the precision on the NMSSM Higgs boson masses and mixings. The next steps to be taken are the inclusion of the full gauge dependence and of non-zero momentum in the computation of the two-loop corrections.

A Mass regulated one-and two-loop functions
This appendix complements the analytically known one-and two-loop self-energy integrals with a set of IR-regularised functions. We implicitly assume vanishing external momentum unless stated otherwise. In this limit, all two-(one-) loop integrals can be written in terms of two-(one-)loop tadpole integrals I(x, y, z) (A(x)). The analytical solution of I has been studied e.g. in refs. [50,136] in detail. While I(x, y, z) and A(x) themselves are IR-finite, their derivatives w.r.t. squared masses become divergent in the IR regime. In this regime, we replace vanishing scalar masses with a mass regulator 0 → M 2 R in the loop functions and expand around the small regulator mass while keeping all terms of order O(log n≤2 M 2 R ) and O(M −n≤0 R ).
The notation and definition of the integrals closely follow those of refs. [115,119]. In addition, we introduce the scalar three-point integral C in terms of the two-point integral B, allowing us to keep track on the spurious IRdivergences C(x, 0, 0) which are not connected to the GBC but cancelled between counterterm-inserted diagrams and genuine two-loop diagrams involving the V-integral introduced later. The required IR-save one-loop functions are Diagrams that factorise into products of one-loop functions can contribute with finite terms like e.g.
where A denotes the one-point function. However, these diagrams would vanish when including the full momentum dependence. Therefore, we always set the one-point function A(0) = 0 before we start with the expansion in a small regulator mass.
We continue with the IR-regulated two-loop functions. For convenience, we define the following abbreviations for derivatives of the tadpole integral I:

JHEP09(2021)193
For the NMSSM, we need the following special cases (A.14) The function T is IR-finite, has been introduced in ref. [119] and is identical to the function R SS (x, y) used in refs. [97,99]. For completeness, we recall here only the expressions needed in the regularisation procedure: With this set of functions we can define all remaining two-loop functions in an IR-regulated way. We start with the UV-divergent U-integral at vanishing external momentum, Therefore we need to regulate We have verified that the UV-IR mixing terms O(logM 2 R / ) cancel in the sub-loop renormalisation with counterterm-inserted diagrams that involve a vertex counterterm and an IR-divergent B(0, 0) integral.
The V-integral V(x, y, z, u) = − ∂ ∂y U(x, y, z, u) (A. 21) can also contribute to the UV-IR mixing terms because its single pole can be written as

JHEP09(2021)193
It precisely cancels against counterterm-inserted diagrams involving a Higgs boson mass counterterm and a three-point function. Furthermore, we need the following special cases of vanishing arguments, The UV-finite two-loop master integral M at vanishing external momentum requires IRregularisation for the following cases, while the finite part of the U-integral is not needed as we renormalise all wave functions in the DR scheme.
Keeping the M, U, V and S integrals even in the zero-momentum approximation (and expressing them in terms of the I integrals during the numerical evaluation) makes the transition to finite/partial external momentum dependence straightforward as we simply replace them by the functions computed with TSIL.  As a closing remark we want to mention that only the O(logM 2 R ) and O(log 2 M 2 R ) divergences are actually physical as they correspond to O( −1 IR ) and O( −2 IR ) terms in dimensional regularisation. As a further cross-check, we explicitly checked that indeed all O(M −n≤−2 R )-terms cancel exactly in the sum of all two-loop diagrams.

B IR-divergent topologies
In this appendix we give all IR-divergent cases for the two-loop tadpole and self-energy topologies ind table 5 and table 6. The first column shows the topology with labels on each generic propagator. The second column lists special cases of vanishing masses in the propagators that lead to IR-divergent loop functions. The third column collects the various sets discussed in section 4.1. The fourth column lists the IR-divergent loop functions that appear in all possible field-insertions after applying the TARCER's algorithm. The last column indicates whether the IR-divergence cancels in the final result or if it requires the inclusion of external momentum.
While the tadpole diagrams cannot be treated with external momentum, the selfenergy diagrams indeed need external momentum in a few cases. In addition, there are cases (such as for instance m 3 = m 4 = 0 in topology 4), where loop integrals do not require momentum-regularisation. In these cases, the IR-divergence was found to cancel against other diagrams that are connected by the BPHZ theorem [137] (such as topologies 4, 7 and 11). Therefore, the subset of topologies 4,7 and 11 with at least one massive Higgs in the outer loop forms an IR-finite set. Similarly to the tadpole diagrams, the topologies 8, 10 and 13 form an IR-finite subset as well. After these considerations, only the diagrams in figure 12 are regularised by momentum.   Figure 10. Generic two-loop tadpole diagrams considered in this work where S denotes scalars and F denotes fermions. Diagrams with red propagators become IR divergent for massless Goldstone bosons.

C.1 Tadpoles
In figure 10 we show all types of tadpole diagrams that enter at our considered loop order.
In the summation over internal fields we consider only those two-loop diagrams that obey the constraint where n Φ is the number of internal propagators of the field type Φ. Diagrams with one loop are taken into account if they contain at least one Higgs boson, electroweakino, stop or top field. In addition, we also consider the diagrams given at O(α 2 t ), cf. ref. [59] figures 14 and 15, but also taking into account terms of the order O(α t α λ ).
The two-loop tadpole diagrams with red colored propagators shown in figure 10 suffer from IR divergences for Φ = G 0 , G ± . Using the mass regularised loop functions from appendix A we find that the sum of all contributions is indeed IR-finite and does not depend on the regulator mass.

C.2 Self-energies
In figure 11 we show all generic self-energy diagrams appearing in the calculation of the O((α t + α λ + α κ ) 2 )-corrected Higgs boson masses. The external fields (s i ) can be the neutral (h i ) or charged (h − 2 ) Higgs bosons as well as massive vector bosons (Z or W − ). The summation over the internal scalars and fermions running in the loops is performed with the same selection principle as for the tadpoles, eq. (C.1). The diagrams (b), (d), (p) and (q) pick up an additional factor of two to account for their mirrored versions.

C.3 Momentum regulated diagrams
For the sake of simplicity, we also list the full set of Feynman diagrams in figure 12 which feature a residual dependence on the IR mass regulator (i.e. the IR divergences in all other diagrams do cancel). Figure 12 can also be inferred from the information given in table 6. si sj S2 S1 (s) Figure 11. All generic two-loop self-energy diagrams with internal scalars S and fermions F . The external fields can be neutral/charged scalars as well as vector bosons (for the latter some diagrams such as (j) and (m) are not present due to gauge invariance and the p 2 = 0 approximation). The internal fields are selected according to eq. (C.1). -47 -