The O(αt+αλ+ακ)2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal{O}(\alpha _t+\alpha _\lambda +\alpha _\kappa )^2$$\end{document} correction to the ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} parameter and its effect on the W boson mass calculation in the complex NMSSM

We present the prediction of the electroweak ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} parameter and the W boson mass in the CP-violating Next-to-Minimal Supersymmetric extension of the Standard Model (NMSSM) at the two-loop order. The ρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} parameter is calculated at the full one-loop and leading and sub-leading two-loop order O(α+αtαs+αt+αλ+ακ2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(\alpha + \alpha _t\alpha _s + \left( \alpha _t+\alpha _\lambda +\alpha _\kappa \right) ^2)$$\end{document}. The new Δρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \rho $$\end{document} prediction is incorporated into a prediction of MW\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_W$$\end{document} via a full supersymmetric (SUSY) one-loop calculation of Δr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta r$$\end{document}. Furthermore, we include all known state-of-the-art SM higher-order corrections to Δr\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta r$$\end{document}. By comparing results for Δρ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \rho $$\end{document} obtained using on-shell (OS) and DR¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\text {DR}}$$\end{document} renormalization conditions in the top/stop sector, we find that the scheme uncertainty is reduced at one-loop order by 55%, at two-loop O(αsαt)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(\alpha _s\alpha _t)$$\end{document} by 22%, and at two-loop O(αt+ακ+αλ)2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(\alpha _t+\alpha _\kappa +\alpha _\lambda )^2$$\end{document} by 16%, respectively. The influence of the two-loop results on the MW\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_W$$\end{document} mass prediction is found to be sub-leading. The new calculation is made public in the computer program NMSSMCALC. We perform an extensive comparison in the W-mass, Higgs boson mass and the muon anomalous magnetic moment prediction between our calculation and three other publicly available tools and find very good agreement provided that the input parameters and renormalization scales are treated in the same way. Finally, we study the impact of the CP-violating phases on the W-mass prediction which is found to be smaller than the overall size of the SUSY corrections.


Introduction
The Standard Model (SM) of particle physics has seen a tremendous success story that certainly culminated in the discovery of the Higgs boson in 2012 by the Large Hadron Collider (LHC) experiments ATLAS [1] and CMS [2].The theory has been tested extensively for its internal consistency at the quantum level.In 2021, the combination of the measurements of the W boson mass has lead to a world average of M exp W = 80.379 ± 0.012 GeV [3].The SM predicts the W boson mass to be M SM, OS W = 80.353 ± 0.004 GeV [4] in the on-shell (OS) renormalization scheme using the stateof-the-art calculations [5], and to be M SM,MS W = 80.351 ± 0.003 GeV in the MS scheme [6,7] with the top mass being 172.76 GeV.This implies a deviation between the SM prediction of the W boson mass and the experimental value of less than 2σ standard deviation, a tension which has been existing between theory and experiment for a long time.In 2022, the CDF collaboration has reported a new result of the W boson mass, which amounts to [8] M CDF W = 80.4335 ± 0.0094 GeV. (1.1) The considerable shift of the central value and the small uncertainties of both the individual CDF measurement and the SM prediction lead to a discrepancy of more than 7σ .
Combining the new CDF result with the other measurements from LEP, Tevatron and the LHC leads to a new world average [9] of and a new deviation of order 6σ . 1 This result has caused a lot of attention in the particle physics community.Many calculations have been performed and analyzed in numerous models beyond the SM in order to resolve this tension.Of particular interest are (minimal) supersymmetric extensions of the SM which introduce fermionic/bosonic superpartners for each SM degree of freedom.In many extensions, these superpartners carry an odd R-parity while the SM particles carry an even R-parity.As a consequence of the imposed Rparity conservation, only an even number of superpartners can contribute to interaction vertices that also involve SM particles.Therefore, amplitudes with only SM-like fields on external legs receive contributions from superpartners at most starting from the one-loop order but not at tree level.
The new particle content beyond the SM predicted by supersymmetry such as e.g. the superpartners of the SM fermions may give significant loop contributions to the muon decay amplitude.Loop corrections to the muon decay are usually parametrised with the quantity r [11] defined through the matching of the Fermi theory and the high-energy theory.Subsequently, the experimentally well-known muon life-time allows to exploit the relation between r , the W/Z masses M W/Z , the Fermi constant G F and the fine-structure constant α to make a precise prediction for M W in terms of the other input parameters.Therefore, the combined W boson mass measurements can be used to constrain those parameters of the considered model that enter r at a given loop level.
Recent studies [4,7,12] in the minimal supersymmetric extension of the SM (MSSM) have shown that light and compressed spectra of electroweakinos can yield M MSSM W up to 80.376 GeV if taking into account the constraints from current LHC supersymmetric (SUSY) searches and limits on Dark Matter (DM) direct detection cross-sections and also satisfying the recent experimental result for the anomalous magnetic moment of the muon.In this work we consider exp W = 80.360 ± 0.016 GeVobtained from their advanced fitting technique at the Rencontres de Moriond conference [10].Thus, more efforts are needed to have better measurements of M W .
the prediction of M W within the next-to-minimal supersymmetric Standard Model (NMSSM) [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28].The model contains two Higgs doublet-and an additional complex singletsuperfield.The W boson mass in this model has been first studied in [29] which includes the full one-loop corrections to r .The authors of [30] improved the calculation by taking into account not only the one-loop but also the two-loop corrections of O(αα s ) [31,32] and O(α 2  t , α t α b , α 2 b ) [33,34] to the ρ parameter.It is important to stress that all contributions to ρ beyond one-loop order have been computed in the MSSM-limit so far.In this work, however, we explicitly calculate ρ including the full dependence on all NMSSMspecific parameters at the two-loop level.
Our group has contributed to the precise calculation of the Higgs boson masses in the real and complex NMSSM by computing the full one-loop [35,36], and the two-loop O(α t α s ) [37], O(α 2 t ) [38] as well as O((α t + α λ + α κ ) 2 ) [39] corrections.In these two-loop computations the W and Z boson self-energies had to be computed as a by-product at the corresponding orders.Therefore, these results can be used to further improve the prediction of the W boson mass in the NMSSM by including not only the two-loop QCD corrections but also the two-loop electroweak ones of O((α t + α λ + α κ ) 2 ), which contain specific NMSSM corrections related to the NMSSM superpotential parameters λ and κ.The ρ parameter and W boson mass calculations are implemented in the Fortran code NMSSMCALC [35][36][37][38][39][40] which is publicly available.The program computes the Higgs boson masses and mixing angles up to two-loop O((α t + α λ + α κ ) 2 + α t α s ), together with the Higgs boson decay widths and branching ratios taking into account the most up-to-date higher-order QCD corrections.The corrections to the trilinear Higgs self-couplings are included at full one-loop level [41] and at two-loop O(α t α s ) [42] and two-loop O(α 2 t ) [43].For the CP-violating NMSSM the calculation of the electric dipole moments (EDMs) has been implemented [44] to be checked against the experimental constraints.Recently, the electron and muon anomalous magnetic moments have been included in the code [45].The code has been extended to include the electroweak corrections to the NMSSM Higgs decays in NMSSMCALCEW [46][47][48].
This paper is organized as follows.In Sect.2, we briefly introduce the used notation to describe the complex NMSSM as well as the tree-level transformations to rotate from gauge to mass eigenstates.The calculation of the one-and twoloop corrections to ρ and the loop-corrected W boson mass will be presented in Sects.3 and 4, respectively.Furthermore, we discuss the renormalization schemes used to obtain UVfinite ρ and r .Section 5 is dedicated to the numerical analysis where we present the size of the loop corrections and discuss their behaviour as a function of the NMSSM-specific parameters.Finally, we give our conclusions in Sect.6.

The NMSSM at tree level
In this section, we give a short description of the complex NMSSM and introduce the notation used in this paper.We follow the same notation which has been used in [36][37][38][39][40].The superpotential of the complex NMSSM is given by (i, j = 1, 2) with the quark and leptonic superfields Q, Û , D, L, Ê, the Higgs doublet superfields Ĥd , Ĥu , the singlet superfield Ŝ and the totally antisymmetric tensor 12 = 12 = 1.Charge conjugated fields are denoted by the superscript c.Color and generation indices have been suppressed for the sake of clarity.The Yukawa couplings y u , y d and y e are assumed to be diagonal 3×3 matrices in flavour space.The parameters λ and κ are in general complex.The soft SUSY-breaking Lagrangian reads with the CP-violating phases ϕ u,s .The relation to the SM VEV v ≈ 246.22 GeV is given by and we define the mixing angle tan β as Thus, the expressions for the tree-level weak gauge boson masses are the same as in the SM, where g 1 , g 2 are the gauge couplings of the U (1) Y and SU (2) L gauge groups, respectively.These couplings can be written in terms of the electric charge e and the weak mixing angle θ W as where the short hand notation c x ≡ cos(x), s x ≡ sin(x), t x ≡ tan(x) is used in this paper.The effective μ parameter is given by (2.10) Besides the particles of the SM, gauge bosons, quarks, charged leptons, and three left-handed neutrino fields, the NMSSM particle spectrum features an extended Higgs sector and new SUSY particles, in particular: • The CP-even and CP-odd Higgs interaction states (h d,u,s , a u,d,s ) mix to form five CP-indefinite Higgs mass eigenstates h i (i = 1, . . ., 5), with their masses per convention ordered as m h 1 ≤ • • • ≤ m h 5 , and one neutral Goldstone boson G 0 .We use a two-fold rotation to rotate from the interaction to the mass eigenstates, (2.11) , and the rotation matrix N transforms the fields ( B, W3 , Hd , Hu , S) T into the mass eigenstates.• The two chargino mass eigenstates, are obtained from the rotation of the interaction states, given by the charged Higgsinos H ± d , H ± u and the charged gauginos W ± = W1 ∓ i W2 / √ 2, to the mass eigenstates.This is achieved by using a bi-unitary transformation with the two 2 × 2 unitary matrices V χ and U χ , • The scalar partners of the left-and right-handed up-type quarks are denoted by ũi L/R , of the down-type quarks by di L/R , and of the charged leptons by li L/R (i = 1, 2, 3).We do not include flavor mixing.Within each flavour the left-and right-handed scalar fermions with same electric charge mix.The sfermions are rotated to the mass eigenstates by a unitary matrix U f .• There are three scalar partners of the left-handed neutrinos, denoted as νi (i = 1, 2, 3) with the masses (2.15) For detailed discussions and derivations of the tree-level spectrum we refer to earlier works [36][37][38][39][40].

The one-and two-loop corrections to the ρ parameter
In this section we describe the calculation of the full one-loop and the O(α t α s ) and O((α t + α λ + α κ ) 2 ) two-loop corrections to the ρ parameter.The ρ parameter is defined by the neutral-to charged-current coupling ratio at zero external momentum transfer [49][50][51] where the neutral current coupling G NC (0) can be determined from, for example, the coefficient of the ν e νe → ν e νe scattering amplitude and the charged current coupling G CC (0) from the e νe → e νe scattering amplitude.In the NMSSM, at tree level we have therefore ρ (tree-level) = 1, as in the SM.Higher-order corrections to the ρ parameter lead to a deviation from unity, which can be written as [51][52][53] ρ = (1) ρ + (2) ρ + • • • , (3.19) where the superscript (n) indicates that the calculation is performed at the n-loop order. (n) ρ can be obtained by computing G CC (0) and G NC (0) taking into account only contributions related to the W and Z self-energies and expanding the ratio G NC (0)/G CC (0) according to the considered loop order.We get the following expressions for the one-loop and two-loop corrections (1) where the transverse part of the gauge boson self-energy ) is evaluated at zero external momentum.These expressions are the same for the SM and for the 2-Higgs Doublet Model (2HDM) presented in [54] and hold in general for every theory with ρ = 1 at tree level.In this paper, we compute the gauge boson self-energies at oneand two-loop level taking into account the NMSSM particles and couplings.These contributions also include the SM-like corrections as a subset.When investigating the size of new physics effects, it is convenient to subtract the SM-like corrections.To be consistent, the computation of the SM-like contributions is performed in the same way as the NMSSMspecific contributions.This strategy can be followed when considering ρ.In the investigation of the W boson mass it is, however, crucial to include existing SM-like higher-order corrections to r beyond the two-loop level (cf.Sect. 4 for a detailed discussion).

The one-loop result
The transverse parts of the gauge boson self-energies are computed with the full content of the NMSSM at zero external momentum without any further approximations.In order not to break SUSY at loop-level, we perform the calculation using dimensional reduction (DR) rather than minimal subtraction (MS).While the results for the individual gauge boson self-energies differ when using MS or DR, we find that the difference of the W and Z self-energies, entering ρ, is the same in both regularization schemes.We confirmed this explicitly at the one-and two-loop level.This is considered as a further consistency check of the results.
We can decompose the ( 1) ρ into four contributions given by the SM fermion, the Higgs, the squark and slepton, and the chargino and neutralino ones, (1) They are separately UV-finite.The SM fermion contributions are the same as in the SM.The contributions from the first two generations can be neglected, and the contribution from the third generation reads The contribution from the Higgs sector, i.e. from G ± , G 0 , h i and H ± , reads where the one-loop functions are given at the end of this section and with c 2W ≡ cos 2θ W , s 2W ≡ sin 2θ W .The contribution from the third generation squarks and sleptons is given by ) , (3.26) where n f = 3 for squarks and n f = 1 for sleptons, t τi / bi = −1/2, t ti = 1/2, Q fi denotes the electric charge of the sfermions.The contribution from the charginos and neutralinos can be cast into the form where In the above expressions, the one-loop integrals are defined as where the bar symbol indicates the dimensionless quantities, i.e. x = x/μ 2 R etc. with μ R being the renormalization scale.

The two-loop O(α t α s ) corrections
For the O(α t α s ) corrections, the first term in Eq. (3.20) will vanish and hence (3.30) We take the results of α t α s ,T Z Z (0) and α t α s ,T W W (0) implemented in NMSSMCALC, which have been computed by our group in Ref. [42] in the complex NMSSM.For the detailed calculation we refer the reader to [42], here we summarize only the main features.We use the gaugeless limit, i.e. we set the electric charge and the W and Z boson masses to zero, The ratio In [42] we have discussed two renormalization schemes for these parameters, on-shell (OS) and DR renormalization.We keep these two options of renormalization schemes here, too.The expressions for the required counterterms were presented in [42].We have explicitly confirmed that the obtained result for α t α s NMSSM ρ is UV finite.For the SM-like contributions we reproduce the known result in the literature [55,56], where the top mass is renormalized using the OS renormalization scheme.Note that all our calculations have been performed in dimensional reduction for both the SM and the NMSSM, while the results in [55,56] were obtained using dimensional regularization.It should be stressed that α t α s ρ in the NMSSM is identical to the corresponding quantity in the MSSM which has been calculated in [31,32].
For the sake of a convenient notation we denote Since the electroweak sector contributes to the gauge boson self-energies at one-loop level, its contribution to ρ at twoloop level contains also the one-loop squared terms, hence The calculation of the transverse part of the gauge boson self-energies at the one-and two-loop order in the complex NMSSM has been presented in [39] and implemented in NMSSMCALC.The results have been obtained in the gaugeless limit, cf.Eq. (3.31).In this limit the Higgs-Goldstone couplings are non-zero at O(α λ + α κ ) while the Goldstone tree-level masses vanish which leads to intermediate infrared (IR) divergences in some of the two-loop diagrams that cancel in the sum of all self-energy diagrams.The top/stop, bottom/sbottom, Higgsino/singlino and Higgs sectors contribute already at one-loop level.The two-loop self-energies (2),T V V consist of contributions from the genuine two-loop diagrams and the counterterm inserted one-loop diagrams.Therefore the parameters of these sectors need to be renormalized at one-loop level to compute the two-loop selfenergies (2),T V V .For the parameters of the Higgsino/singlino and Higgs sectors we apply a mixed DR-OS renormalisation scheme while for the top/stop sector we apply an OS scheme or a DR scheme.All counterterms for the complex phases ϕ α (α = s, u, κ, λ) can be set to zero.The remaining input parameters together with the applied renormalization conditions are given by if Re A λ is chosen as independent input rather than M H ± .The tadpole counterterms are defined such that the minima of the Higgs potential do not change at higher order.We explicitly confirmed that α 2 EW ρ is UV finite and is free of IR divergences (cf.Ref. [39] for a detailed discussion on the cancellation of all IR divergences).
The corresponding result α 2 EW SM ρ within the SM, which is subtracted from the NMSSM result, is computed in the same way as described above, i.e. using OS conditions for v, sin θ W and OS/DR conditions for m t .In the SM only contributions from top and Higgs diagrams enter at O(α 2 EW ), hence EW SM ρ does not depend on the renormalization of the Higgs mass and the Higgs tadpole.Note that our expression of α 2 t SM ρ in Eq. (3.35) looks different from the one computed by Fleischer, Tarasov and Jegerlehner in [53] for the following reasons.In their computation, ρ at two-loop order is defined as where V V (0) gets contributions only from the genuine two-loop diagrams and the counterterm inserted diagrams, i.e. the one-loop-squared pieces from Eq. (3.20) are not contributing.This is because they used the Fermi constant G μ as an input parameter.Therefore, in [53] only the top mass needs to be renormalized in order to obtain a UV finite result.However, in the following, we use the VEV v (defined through M W , M Z and α(0)) as an input.Considering the one-loop expression for ρ in the SM, which involves both m t and v, the two-loop result needs to include the renormalization of m t and v. Since v is renormalized in the OS scheme, its counterterm contributes with a non-zero finite part of δv which is not present in the result of Eq. (3.38).The UV-finite part of δv also gives rise to non-vanishing contributions to the single pole of the counterterm inserted one-loop diagrams.This additional contribution is precisely canceled by the one-loopsquared term in Eq. (3.35).Our final results for α 2 t SM ρ are the same as the results of [53] The Higgs mass in the SM result SM ρ is set to be equal to the possibly non-zero tree-level SM-like Higgs boson mass in α 2 EW NMSSM ρ obtained in the gaugeless limit.This is a crucial difference to the MSSM, where the SM-like Higgs boson mass in the gaugeless limit always vanishes.

Calculation of the W boson mass in the OS scheme
The W boson mass can be computed from the following relation [11], where G μ is the Fermi constant, α ≡ α(0) = e 2 /4π is the fine-structure constant in the Thomson limit, and r includes all loop contributions to the amplitude of the μ → e νe ν μ decay after subtracting the Fermi-model-type QED correc-tion.By using the OS weak mixing angle s where the NMSSM r is taken from our implementation in NMSSMCALC.The quantity NMSSMCALC r , the subscript (n) denotes the n-loop order, depends also on M W .Therefore, Eq. (4.40) has to be evaluated iteratively.In the following subsections we describe in detail the NMSSM-specific oneand two-loop contributions as well as the higher-order SM corrections which are included in r .

One-loop corrections
The one-loop correction to r can be written as (1) where the origin of the individual terms in this expression is explained in the following.Using the OS renormalization scheme for the electric charge, its counterterm is given by where and sgn = 1 for the SM and sgn = −1 for the NMSSM accounts for different sign conventions in the covariant derivatives of the two models.The photon self-energy in Eq. (4.42) evaluated at vanishing external momentum contains large logarithmic contributions, O(ln(μ 2 R /m 2 q )), which depend on ratios of the light quark masses and the renormalization scale μ R , leading to numerical instabilities.To avoid this dependence, one writes for the light SM fermions where Re A A f (M 2 Z ) is computed perturbatively, A A rem (0) contains contributions from remaining charged particles of the model and α = α lepton + α (5) had .The contribution 123 (2023) 83:1079 α (5) had is determined from the dispersion relation using experimental data.We take α (5) had = 0.02768 from Ref. [57].The quantity α lepton = 0.03150 includes contributions up to three loops [58].The counterterm of s W is derived from the OS relation s where the Z and W boson mass counterterms δ M 2 Z and δ M 2 W are defined in the OS scheme, such that M W and M Z correspond to the real part of the complex pole of the W and Z propagator with a constant decay width W and Z , respectively.However, experimentally the resonances are parametrized using a Breit-Wigner line shape that features an energy-dependent rather than fixed decay width.To account for this, we convert the vector boson masses that correspond to the standard OS definition to the running-width masses M run W/Z = M W/Z + 2 W/Z /(2M run W/Z ) as described in [59].For Z we use the experimentally measured value while for the W boson width we use the approximate formula [60] neglecting possible beyond-the-SM (BSM) contributions.While the conversion from M run Z to M Z needs to be done only once before the iteration in Eq. (4.40), the conversion from M W to M run W needs to be done at the end of the iteration.For more details on the running/fixed-width mass conversion we refer to [30,60].
The wave function (WF) counterterms of the external leptons, δ Z l , and the triangle and box contributions, r , , are evaluated for vanishing lepton masses.The WF counterterms of the external leptons are defined in the OS scheme, where L ll (0) is the form-factor of the left-handed vector component of the lepton self-energy.The Fermi-model type QED correction appears in δ Z l , the triangle and the box diagrams.In order to remove this contribution, one should replace the photon propagator in the loop diagrams of the WF counterterms by and remove the box diagrams with virtual photons.The vertex diagrams with virtual photons are calculated as usual.This procedure has been first introduced in Refs.[11,61].Using this trick we recover the well-known relation between triangle, box and vertex corrections for the SM, derived e.g. in [62], Collecting all previous derivations leads to the known SM result which can be written as where α and (1) ρ contain the numerically dominant part of the corrections.Within the complex NMSSM, we have computed the one-loop contributions to (1) NMSSM r using the expression in Eq. (4.41) and using the resummation of the light fermions in the photon self-energy described above.We have checked that the resulting expression is UV-finite and renormalization scale independent.Contributions from squarks, sleptons, charginos and neutralinos form separate UV-finite subsets and can be studied independently.On the other hand, contributions from gauge and Higgs bosons must be combined in order to obtain UV-finite results.

Combination with known higher-order corrections
In the SM, corrections lit.
SM r up to four-loop order have been computed and are available in the literature.Following the procedure developed for the MSSM [63] and implemented in FeynHiggs [64][65][66][67][68][69][70][71], we have combined all available higher-order SM corrections with the presented NMSSM contributions.For fixed n (n = 1, 2) we first subtract NMSSM r and re-add lit.SM r including all known-higher corrections from the literature, and Note that (n) SM r should be computed in the same way as for The SM corrections from the literature consist of the following terms lit.
SM r = (1) r For a complete list of references we refer the reader to [63].
In the following, we refer to previous works that provide analytical/numerical results which are used in Eq. (4.55).
The full two-loop QCD corrections (αα s ) r are taken from [72], partial three-loop QCD corrections (αα 2 s ) r are taken from [73].The three-loop corrections to the ρ parameter are taken from [74] and the four-loop QCD corrections (G 2 μ m 2 t α 3 s ) r from [75,76].For the full twoloop electroweak corrections (αα 2 s ) r we use the fitting formula presented in [77].For consistency, we need to use the running-width definition of the W and Z masses in the twoloop electroweak corrections, cf.footnote 7 in [30], while the fixed-width masses are used in the rest of the r calculation due to the employed OS scheme.

NMSSM-specific two-loop corrections
As discussed in the previous sections, corrections to r can be categorized into corrections arising from the renormalization of the electric charge, vertex and box diagrams and other corrections arising from the W and Z boson self-energies.For the first three categories, higher-order corrections in the (N)MSSM beyond the one-loop level are unknown and not considered in this work.This should be a good approximation for most phenomenologically viable scenarios since these corrections are expected to be numerically small.The dominant contributions arise from the W and Z boson selfenergies and can be parameterized in terms of the loop corrections to the ρ parameter.Thereby, the two-loop corrections to the ρ parameter, SUSY ρ, discussed in Sect.3, can be computed in an efficient way.However, (2)  NMSSM ρ also includes the SM-like corrections which are already taken into account in lit.
SM r .Hence, the SM-like contributions have to be subtracted from the ρ parameter in order to avoid double counting.We define the two-loop SUSY correction to ρ as

The final expression for
(2) SUSY r including the one and two-loop corrections, that is used together with Eq. (4.54) in Eq. (4.40) to compute the W mass in NMSSMCALC, reads where (1) SUSY r is defined via Eq.(4.53).To cross-check the implementation of the known SM corrections as well as to check the correct decoupling behaviour of the NMSSM-specific corrections, we verified that our NMSSM W -mass prediction approaches the one in the SM, as found in [30,63], once all NMSSM particle masses are chosen to be well above the electroweak scale.

Numerical analysis
In this section, we investigate the phenomenological impact of the O((α t + α λ + α κ ) 2 ) corrections on the electroweak ρ parameter and on the W boson mass.Detailed studies of the one-loop and the leading two-loop corrections can be found in [7,29,30].Furthermore, we estimate the uncertainty due to missing higher-order SUSY corrections to ρ by changing the renormalization scheme of the top/stop sector.Possible sources of uncertainties entering the M W prediction are discussed as well.For illustrative purposes, the results are shown for a set of two parameter points, one obtained from a simple scan, which will be described in the next paragraph, and one taken from the literature.For a full investigation of the viable parameter space where the NMSSM can simultaneously explain the (g − 2) μ data, the Higgs data, the W mass, and where the lightest neutralino is still a good lightest supersymmetric particle (LSP) and Dark Matter candidate, we refer the reader to the recent studies [78,79] for the NMSSM, and to [4] for the MSSM.Finally, we compare the W -mass prediction in NMSSMCALC with the one obtained from the public NMSSM spectrum generators FlexibleSUSY, NMSSMTools and SARAH/SPheno.

Setup of the parameter scan
In order to find parameter points which are not excluded by the measured properties of the 125 GeV Higgs boson and by the direct LHC searches for SUSY particles, we have performed a scan over the NMSSM parameter space in the following way.For a given set of input parameters, we use NMSSMCALC to calculate the mass spectrum of all SUSY particles at the one-loop order.In addition, the Higgs boson masses are obtained including the available two-loop corrections at O(α s α t + (α t + α λ + α κ ) 2 ) [39].The Higgs boson masses and mixing matrices are computed using OS renormalization in the top/stop sector.One of the neutral CP-even Higgs bosons is identified with the SM-like Higgs boson, and its mass is required to lie in the range 122 GeV ≤ m h ≤ 128 GeV. (5.58) In the following, m h will always denote the mass of the SMlike Higgs boson which not necessarily needs to be the mass m h 1 of the lightest scalar state.We further require that (1) the lightest SUSY particle is the lightest neutralino; (2) the masses of the electroweakinos and sleptons satisfy the lower bounds from LEP [57]; (3) the masses of the stops, lightest sbottom and the gluino satisfy the lower bounds taken from the 2022 Particle Data Group review [57], in particular (5.59) The Higgs decay widths and branching ratios including the state-of-the-art higher-order QCD corrections as well as the effective Higgs couplings, i.e. using the Higgs mixing angles obtained from the diagonalization of the loop corrected mass matrices, are obtained with NMSSMCALC, too.Having all important properties of the Higgs sector at hand, we use HiggsTools [80] which contains HiggsBounds-5 [81], to check if the parameter points pass all the exclusion limits from the searches at LEP, Tevatron and the LHC, and HiggsSignals-2 [82] to check if the points are consistent with the LHC data for a 125 GeV Higgs boson within 2σ .To constrain the SUSY fermionic and scalar sector we use SModelS-2 [83] to check whether a given scenario is excluded by the LHC searches for the electroweakinos and sleptons.For the input of SModelS-2, we have implemented the chargino and neutralino decay widths and branching ratios in a private version of NMSSMCALC.Furthermore, the production cross sections of all pairs of electroweakinos have been computed at leading order using a private implementation that was obtained with the help of FeynArts-3.11 [84] and FormCalc-9.8 [85].These tree-level cross sections are corrected by a common K -factor to account for the NLO QCD corrections.We chose K = 1.3 which was obtained with the help of Prospino2.0[86,87].Since there are many experimental constraints applied in the scan, using a uniform random scan over all input parameters can be very time and resource consuming.We therefore performed a Markov Chain Monte Carlo sampling using EasyScan_HEP-1.0[88] in order to efficiently find phenomenologically viable parameter regions.Table 1 summarizes the ranges applied in the parameter scan.
Note that all parameters are chosen to be real in the parameter scan.The SM input parameters, that are relevant for the calculation of ρ and the W mass, 2 are taken from [57] α 0 = 1/137.035999084,α MS s (M Z ) = 0.1179, 2 In the calculation of all Higgs masses α(M Z ) = 1/127.955and M W = 80.377 GeV are chosen as input parameters.18 GeV, m τ = 1.77686GeV, α (5)  had (M 2 Z ) = 0.02768, α lepton (M 2 Z ) = 0.03150. (5.60) For the following numerical analysis, we have chosen two parameter points to study the impact of the new corrections.The first parameter point, called P1, which passes all our constraints specified above, is given by the following input parameters, (5.61) The resulting spectrum of the Higgs boson masses at O(α s α t + (α t + α λ + α κ ) 2 ) is given in Table 2.The SM-like Higgs boson is dominated by the h u component and its mass is about 125.4 GeV, the remaining Higgs bosons are heavier.The spectrum of the electroweakinos is given in Table 3.For this parameter point, the lightest neutralino is the lightest SUSY particle.While the masses of the electroweakinos are rather light with masses below about 510 GeV, the masses of sleptons are rather heavy, larger than 380 GeV.The lightest sleptons are mainly composed of the right-handed selectron and smuon components.The special feature of this point is that the wino, bino and higgsino components mix strongly.Thus, one can not distinguish which neutralino is wino-, binoor higgsino-like, which is a region where experimental constraints are rather weak.In particular, the cross sections of the electroweakino pair production processes become small and therefore the parameter point escapes the LHC constraints from the electroweakino searches.

Results for the ρ parameter
In the following, we investigate the prediction for the ρ parameter starting from the parameter point P1 as a function of the NMSSM-specific superpotential parameters λ and κ.In order to avoid negative squared tree-level masses, we simultaneously vary both λ and κ around their values λ 0 = 0.301 and κ 0 = 0.299 such that the ratios λ/λ 0 and κ/κ 0 , respectively, are kept equal.Furthermore, we vary ) to avoid a tachyonic tree-level Higgs mass even for very large values of λ.All other parame-Table 1 Input parameters for the NMSSM scan.All parameters have been varied independently between the given minimum and maximum values (5.62b) Since the prediction of the ρ parameter strongly depends on the top quark mass, it is shown using both the OS (full lines) and the DR (dashed lines) renormalization schemes in the top/stop sector.The W boson mass entering the ρ parameter prediction has been chosen to be the value for M W obtained at O(α 2 new ), cf.Sect.5.3.Note that points with s λκ > 0.7 generally are in danger to violate perturbative unitarity and/or to run into a Landau pole close to their input scale.Nonetheless, it is possible to go to larger values in some regions of the parameter space (giving up the requirement that the NMSSM is UV-complete).From this point of view, it is interesting to study which values of λ, κ would be required in order to obtain sizeable corrections.We therefore plot our results for values up to s λκ = 3.In the OS scheme, the O(α t α s ) and O(α 2 new ) corrections are both negative.In the region that is free of low-energy Landau poles, the O(α t α s ) corrections reduce the one-loop result by 10% while the O(α 2 new ) corrections reduce the O(α t α s ) results by 4%.In contrast, the O(α t α s ) and O(α 2 new ) corrections are both positive in the DR scheme.The O(α t α s ) corrections increase the one-loop value by 9.8% while the O(α 2 new ) corrections add 0.6% on top of the O(α t α s ) results.In the region s λκ > 2, the O((α t + α λ + α κ ) 2 ) corrections in the two renormalization schemes start to deviate differently from the O(α t α s ) corrections.The magnitude of the O((α t + α λ + α κ ) 2 ) corrections becomes smaller in the OS scheme, but larger in the DR scheme.
We observe an increase in ρ with increasing s λκ which first is rather weak and becomes stronger for very large values of s λκ .This behaviour is correlated with an increase of the SU (2) mass splittings between the neutral and charged Higgs bosons on the one side and the neutral and charged electroweakinos on the other side.This can be inferred from the right plots, where we show in the upper plot the dependence of the charged Higgs mass and the heavy CP-even/odd Higgs masses, which have a dominant h d and a d component, respectively, as a function of s λκ , while the masses of their corresponding EW-ino states as function of s λκ are shown in the lower plot.
We conclude this section with a discussion about numerical instabilities that can appear in the O((α t + α λ + α κ ) 2 ) corrections to ρ, especially for small values of λ and κ.The corrections to the ρ parameter are composed of W and Z boson self-energies which individually can receive very large higher-order corrections.However, in the difference of the self-energies, entering ρ, large cancellations of many orders of magnitude can appear.In some cases the size of this cancellation may exceed the numerical precision of the program which is currently limited to double precision.The O((α t + α λ + α κ ) 2 ) corrections are particularly sensitive to this kind of instabilities since the tree-level masses, which enter the two-loop self-energies, are calculated at O(α λ + α κ ) in the gaugeless limit which can lead to very small but non-zero tree-level masses that enter the calculation of the two-loop self-energies.While these instabilities might not be apparent in the ρ parameter prediction in the first place, the M W prediction which discussed in Sect.5.3, is very sensitive to the value of ρ.In particular the VEV (and all parameters and couplings derived from it) depends on the floating value of M W during the iteration which can easily amplify numerical instabilities coming from ρ.Therefore, the M W prediction in NMSSMCALC at O((α t + α λ + α κ ) 2 ) can only be used for parameter points that do not suffer from large numerical instabilities.If no convergence is found at (2023) 83:1079   2 ) in M W , the program automatically falls back to the O(α t α s ) predictions for both the W -mass and ρ.Furthermore, some parameter points may feature tachyonic tree-level masses in the gaugeless limit, while the actual treelevel masses are positive.While Ref. [89] proposed a possible solution to a similar problem at the one-loop order, our calculation is mostly affected by tachyonic states required for partial two-loop corrections.The generalisation of the method developed in Ref. [89] is beyond the scope of this work and left for future work.Therefore, if the program encounters negative squared tree-level masses in the preparation of the masses entering at a given loop-order, it also automatically falls-back to the next-lowest order that is not involving tachyonic tree-level masses in the Feynman diagrams. 3 In Fig. 1 we observe good convergence for s λκ > 0.1 for M W at O(α 2 new ), cf.Sect.5.3.For s λκ < 0.1 the M W prediction does not converge and therefore the ρ parameter prediction is used at O(α t α s ) which explains the jump of the red line onto the blue line.

Uncertainty estimate for ρ
In order to estimate the uncertainty due to missing higherorder corrections to the ρ parameter, we define the renormalization scheme dependence of the ρ parameter at a given loop order as 3 If the full tree-level masses, which enter the full one-loop mass prediction, are found to be tachyonic, the program exits with an error message.The lower panel of Fig. 1 (left) shows the resulting ρ ren obtained at the three considered loop orders as a function of s λκ .We observe a renormalization scheme dependence of up to 50%, 22% and 16% for the one-loop, O(α t α s ) and O(α 2 new ) results, respectively.Therefore, including the twoloop QCD and EW corrections can significantly reduce the theory uncertainty of the ρ parameter.For the comparison with the SM we present in Table 4 the ρ parameter computed in the SM at the corresponding one-loop order and two-loop O(α t α s ) and O(α t α s + α 2 t ) using DR and OS renormalisation conditions in the top/stop sector.In the SM, the O(α t α s ) and O(α t α s + α 2 t ) corrections are negative both in the OS and the DR scheme.The renormalisation scheme dependence ρ ren in the SM is 55%, 45% and 42% at one-loop order, O(α t α s ) and O(α t α s + α 2 t ), respectively, which is significantly larger than the corresponding results obtained in the NMSSM.Therefore, the SUSY QCD contributions to the ρ parameter seem to play an important role in the reduction of the scheme dependence of the ρ parameter.We want to stress that it is not possible to draw conclusions about the scheme uncertainty of the M W prediction from the scheme uncertainty of the ρ parameter.In particular for the SM prediction of M W a full two-loop prediction, beyond corrections to ρ, was found to yield only an uncertainty of O(3 − 7 MeV ) [6] implying cancellations of the scheme dependence between ρ and other quantities entering r .For a discussion of the M W uncertainty in the NMSSM, see the next section.

Results for the W boson mass
In this section we discuss the prediction of M W at one-loop order, two-loop O(α t α s ) and O(α 2 new ).It is important to stress that a consistent inclusion of the known higher-order SM corrections to r beyond two-loop order, cf.Sect.4.2, is only possible if the SM sector entering the SUSY corrections is renormalized in the same renormalization scheme as in the SM calculation that has been implemented in NMSSMCALC.Therefore, we exclusively chose the OS scheme in the top/stop sector entering the prediction of M W .
We define two quantities to investigate the behaviour of and W . (5.64b) The quantity SM m h defines the difference between the NMSSM W -mass prediction at a given order and the SM prediction, including all SM higher-order corrections, evaluated with the Higgs boson pole-mass m h that is predicted by the NMSSM for the considered scenario at the two-loop level.Therefore, also M SM W (m h ) varies with λ and κ as the Higgs mass prediction in the NMSSM changes.Since the SM higher-order corrections drop out in SM m h , it can be used as a measure for the size of the genuine SUSY corrections to M W .The quantity α i α j determines the size of specific higher-order SUSY correction α i w.r.t. the next-lowest order α j .
Figure 2 (upper left) shows SM m h (left axis) as a function of s λκ = √ λ 2 + κ 2 starting from the parameter point P1 at one-loop (black solid) and two-loop O(α t α s ) (blue solid) and O(α 2 new ) (red solid).Note that we have used the same procedure for the variation of the parameters described in the previous section.The upper axis shows the obtained value for m h , including all available two-loop corrections, that is used in the prediction of M SM W (m h ) (green dot-dashed, right axis).The lower panel shows α t α s one-loop (blue) and 2 new α t α s (red).We observe that the NMSSM-specific corrections range between about 10-20 MeV, mostly dominated by the one-loop correc-tions.The QCD corrections to the W boson mass are negative compared to the one-loop prediction and subtract about 0.2 MeV from the one-loop result, independent of the value of s λκ .Compared to the O(α t α s ) result the O(α 2 new ) corrections range between −0.2 MeV and +1.5 MeV in the shown range of s λκ .
To get a better understanding of the individual contributions to the W -mass prediction we plot the values of r obtained after the M W iteration has converged.In Fig. 2 (right) the blue solid line shows the total result of r obtained with NMSSMCALC including all available corrections.The green dotted line shows the one-loop SM-contributions which also contain the α contributions that are numerically most significant.The green dashed line shows the size of the higher-order SM results taken from the literature which are the next-to-largest contribution to r .The third-largest contribution is the one-loop SUSY contribution (black dash-dotted) which is negative (hence a positive shift to M W ) followed by the -also negative -EW contributions (red dash-dotted).The SUSY QCD corrections are positive and numerically in competition with the EW corrections for s λκ 0.6 − 0.9.Note that in Fig. 2 (right) we chose a logscale for | r | > 10 −5 and a linear scale otherwise which sets the focus on the s λκ -dependence of the two-loop SUSY corrections.
For the other parameter points which pass the constraints of our scan described in Sect.5.1 we observe similar features.The O(α t α s ) and O(α 2 new ) corrections to M W are rather small and their sign w.r.t. to the one-loop corrections can change, depending on the considered parameter point.

Uncertainty estimate for M W
In this paragraph we comment on the uncertainties that contribute to the M W calculation in NMSSMCALC.We focus on two sources of uncertainties: (i) parametric uncertainties that are introduced through the dependence on experimental input parameters and (ii) theory uncertainties due to the missing higher-order corrections.To estimate (i) we vary the SM input parameters M Z , α (5) had (M 2 Z ), α s within their 1σ PDG values [57].The top quark mass is varied by 1 GeV.Table 5 lists the maximal differences in the W mass prediction compared to the result obtained at the central values for the parameter points P1, BP3 (which is introduced in Sect.5.4) and the SM prediction that includes all known higher-order corrections.Furthermore, we also vary the result of the loop-corrected Higgs boson mass by 1 GeV to account for the theoretical uncertainty in the Higgs mass prediction which is indirectly influencing the prediction of the W boson mass.We observe that the SUSY prediction does not introduce a significantly larger parametric uncertainty compared to the SM prediction.The values for the latter are in good agreement with those found in [4,7].
Fig. 2 Upper left: the difference SM m h (left axis) defined in Eq. ( 5.64a) at one-loop (black), two-loop O(α t α s ) (blue) and two-loop O(α 2 new ) (red) as a function of √ λ 2 + κ 2 , and M W in the SM (green dash-dotted, right axis) using the loop-corrected NMSSM Higgs boson mass predic-tion (upper axis).Lower left: the difference αt αs one-loop (blue) and  5 The variation in the W mass prediction around the central value when varying M Z , α (5)  had (M 2 Z ), α s within 1σ of their central values and m t within 172.69 ± 1 GeV, for P1 and BP3.The Higgs mass m h used in the SM-prediction is chosen to be the loop-corrected Higgs mass of the respective parameter point and varied by ±1 GeV Combining all uncertainties in Table 5 quadratically yields a total parametric uncertainty of about 6 MeV.To estimate the size of the missing of higher-order corrections, we first divide them into SM-like and SUSY corrections and discuss their individual uncertainties.Uncertainty estimates for the SM corrections have been studied in [4][5][6] and yield about 4 MeV and 3 MeV in the OS and MS calculation, respectively.A comparison between OS and MS result, however, suggest an uncertainty of about 6 MeV [6].Regarding the missing higher-order SUSY corrections, one can expect them to not be significantly larger than the computed O(α t α s ) and O(α 2 new ) corrections in large parts of the parameter space.For the parameter points of our scan which pass the applied constraints, the maximal corrections for M W are about 4 MeV and 2 MeV, respectively.However, there remains the possi-bility that two-loop SUSY corrections proportional to the electroweak gauge couplings, which are unknown so far, could be enhanced in cases where SU (2) states have large mass splittings.Therefore, the SUSY uncertainty should be at least about 4 MeV large.

Comparison with previous M W results
Among publicly available tools, not only NMSSMCALC is able to predict the W -mass in the NMSSM at high precision.In light of the CDF measurement, the spectrum-generator generators FlexibleSUSY [90,91] and SARAH [92][93][94][95] have been updated to be able to predict M W in a wide class of BSM models.So far, the updated implementation in SARAH/SPheno has been used to study M W in Dirac gaugino models [96] while the one of FlexibleSUSY was applied to the MSSM, MRSSM and a singlet extended SM [7].However, they are in principle also capable to generate spectrum-generators that allow to study M W in the NMSSM.Furthermore, the program NMSSMTools was extended in [29] to compute the W -mass in the general NMSSM as well as the Z 3 -symmetric NMSSM described by Eq. (2.3).In the following, we briefly review the main ingredients for the M W prediction implemented in FlexibleSUSY, NMSSMTools and SARAH/SPheno while focusing on treatments that are different from the NMSSMCALC implementation described in Sect. 4. Following the discussion of the four different M W calculations in the NMSSM, we numerically compare the prediction obtained for two concrete benchmark points.Higher-order corrections to the muon anomalous magnetic moment a μ are known to have a connection to large corrections to M W [4]. Since a μ is of increasing recent theoretical and experimental interest [78,[97][98][99], we also include a comparison of the a μ prediction between the various codes.
We start by discussing the incorporation of the SM higherorder corrections in the different codes.NMSSMCALC implements the results of lit.
SM r to a large degree analytically (cf.Sect.4.2), while the other three codes are based on fit formulas for the M W prediction within the SM.In contrast to NMSSMCALC, FlexibleSUSY and SARAH/SPheno compute all BSM corrections to M W in the MS/DR scheme and therefore rely on the SM MS fit formula for M SM W provided in [6]: SUSY , (5.65) where M SM fit.
W is a numerical fit that incorporates the SM higher-order corrections as a function of the SM input parameters.It is important to stress that the implicit dependence of lit.
SM r on the value of M W , which is correctly taken into account in NMSSMCALC via Eq.(4.40), is lost when using fit formulas for M W .In NMSSMTools the M W dependence is partially restored by determining M SM fit.W from the fit formula given in [5], inverting Eq.
SM r (M W ). (5.66b) We now discuss the different treatments of SUSY input parameters.All (SM and BSM) quantities entering Eq. (5.65) are defined in the MS/DR scheme.In FlexibleSUSY, Eq. (5.65) is evaluated with all running parameters at M Z while in SARAH/SPheno it is evaluated using parameters defined at the SUSY input scale.Thus, SARAH/SPheno is closer to the approach of NMSSMCALC and NMSSMTools which compute M W using the running SUSY input parameters that are given at the SUSY input scale M 2 SUSY = m tR m Q3 .In order to account for this systematic difference in the M W calculation, we modified the spectrum-generator generated by FlexibleSUSY to compute M W and a μ at M SUSY rather than M Z .Likewise, SARAH/SPheno computes a μ by default at M Z which we changed to M SUSY .While NMSSMCALC is in principle able to renormalize parts of the SUSY sector OS, a consistent comparison among the different tools is easiest when interpreting all SUSY parameters as DR parameters defined at M SUSY .In particular, the parameter Re A λ is used as input parameter in the DR scheme, in contrast to the previous sections.Likewise, the renormalization of the top/stop sector for the calculation of the Higgs boson masses is performed in the DR scheme.The M W prediction in NMSSMCALC, however, is still carried out in the OS scheme for the SM-sector described in Sect. 4. NMSSMTools interprets tan β per default to be defined at M Z rather than M SUSY .Thus, we run the tan β(M SUSY ) used in the other codes to M Z (using two-loop RGEs generated with SARAH) and use tan β(M Z ) in NMSSMTools.Furthermore, many of the involved codes also compute loopcorrected masses to sfermions and electroweakinos which in turn are used in the M W prediction.We also disabled such calculationsin all programs as far as possible.Another important ingredient for the W mass prediction in supersymmetric theories is the Higgs mass prediction.For a detailed discussion of the different ingredients to the NMSSM Higgs mass prediction as well as comparisons between the various spectrum generators we refer to [100,101].In the context of the W -mass prediction, a common approach is to use the loop-corrected Higgs mass, which is around 125 GeV in phenomenologically viable scenarios, in the SM-part of the M W calculation.This ensures that, for a parameter point with m (loop) h ≈ 125 GeV, the NMSSM in the decoupling limit predicts the same numerical value for M W as the SM.This is the approach implemented in FlexibleSUSY, NMSSMCALC and SARAH/SPheno.NMSSMTools uses a fixed value of m h = 125.2GeV in the SM fit formula.Furthermore, the SM-like Higgs boson is not necessarily the lightest scalar in the spectrum since the singlet-like states can in principle be lighter.For this reason, NMSSMCALC automatically determines the SM-like Higgs boson (based on the structure of the mixing matrix) which is to be used in the SM part of the calculation.In case of FlexibleSUSY, this information can be given via the SLHA input file by the user.To our best knowledge, SARAH/SPheno does not have a mechanism to determine the SM-like Higgs boson in the M W calculation but always assumes it to be the lightest scalar state.Since we perform the comparison between the different programs using a parameter point which has a singlet state lighter than 125 GeV, we modified the SARAH generated SPheno code such that is also takes the index of the SM-like Higgs boson as additional input in the SLHA input file.
By construction the different methodologies to determine M W described by Eq. (4.40), Eq. (5.65) and Eq.(5.66), yield the precise SM result in the limit of heavy superpartners (see the discussion in [7] for more details on the correct decoupling behaviour).We explicitly verified that all codes still (2023) 83:1079 feature a proper decoupling behaviour after we applied the changes discussed above.
In order to compare the m h , M W and a μ prediction numerically between the four codes FlexibleSUSY 2.7.1, NMSSMCALC 5.2.0, NMSSMTools 6.0.0 and SARAH 4.15.1, we chose two parameter benchmark points.The first point, P1, was already discussed in Eq. (5.61) and features rather light electroweakinos.The second, BP3, was introduced in [79] 4 .This parameter point is characterized by large one-loop corrections to the W boson mass due to very light sleptons with masses of O(100 GeV).Another interesting feature of this point is that the singlet-like CP-even and -odd Higgs bosons have masses lighter than 50 GeV.Therefore, both chosen parameter points give the opportunity to compare the four considered M W codes under rather extreme conditions.For convenience, we reprint the input parameters for BP3 from [79] using the notation introduced in Sect.2, (5.67) In Table 6 we compare the prediction for m h , M W and a μ obtained with all four codes for P1 and BP3.Additionally, we list the uncertainty estimates for M W and a μ that are returned by the programs NMSSMTools and FlexibleSUSY.The uncertainty for M W in NMSSMTools consists of a parametric uncertainty which is dominated by the top quark mass (cf.Sect.5.3), the uncertainty due to the use of the fit formula in Eq. (5.66) and a SUSY uncertainty of about 5 MeV.
The uncertainty for m h due to missing higher-orders can be estimated to be at least about 1 GeV [38,39,101].Since the parameter point P1 is now interpreted with stop masses and stop-Higgs trilinear couplings defined in the DR scheme, rather than OS, the obtained Higgs boson mass for P1 is no longer at 125 GeV but around m h ≈ 119 GeV.Despite the fact that there are many differences in the treatment of the parameters between the four programs, the obtained results for the loop-corrected SM-like Higgs boson mass, for the anomalous magnetic moment a μ and for M W are overall in good agreement.In particular the M W prediction is in agreement between all four codes even if we only impose an uncer-tainty of 4-7 MeV which is about the size of the M W uncertainty of the SM prediction, cf.Sect.5.3.The result for a μ obtained with NMSSMCALC for P1 is far outside the claimed uncertainty of NMSSMTools and FlexibleSUSY, which were obtained by changing the scale at which a μ is calculated between M SUSY /2 and 2M SUSY .This is due to the scale-choice used in the calculation of a μ which is fixed in NMSSMCALC to μ 2 = m L2 m R2 while in the other codes it is chosen dynamically to be the smallest mass of all sleptons and electroweakinos.For P1, the lightest SUSY particle is an electroweakino while for BP3 it is a slepton.Thus, we have better agreement for BP3 than for P1.We explicitly verified that NMSSMCALC predicts a similar number for P1 of a μ ≈ 3 × 10 −10 when using the lightest neutralino mass as the renormalisation scale.
In the following we compare the dependence of the M W prediction onto the superpotential parameters λ and κ between the different codes.In Fig. 3 (left) we plot the M W prediction obtained with the four codes.In the right plot we show the obtained value for m h to validate if similar Higgs mass values are used in the respective SM higher-order results.We observe that M W obtained with FlexibleSUSY and SARAH/SPheno agree almost perfectly and only start to slightly deviate for very large values of s λκ 0.6.The level of agreement in M W , however, always need to be seen in the light of the m h prediction.In case of FlexibleSUSY and SARAH/SPheno the Higgs mass predictions differ by about 0.8-1 MeV which means that their SM-prediction actually differs by 0.2-0.5 MeV such that the perfect agreement in Fig. 3 (left) seems accidental.Furthermore, they never differ from the NMSSMCALC prediction for M W by more than 0.63 MeV.NMSSMTools, however, seems to always predict a W mass which is about 2 MeV larger even though its Higgs mass prediction is also relatively close to the other codes.We perform a similar analysis for the parameter point BP3 in Fig. 4. For this point we plot, in addition to the SMlike Higgs boson mass m h 2 , the mass of the lightest CPeven and -odd Higgs boson, m h 1 and m a 1 , respectively.For all three shown scalar masses we find good agreement and a similar behaviour in terms of s λκ which is rather flat in the shown range of s λκ ≤ 0.3.For larger values of s λκ the singlet-like CP-even state becomes tachyonic at tree-level.The M W prediction of FlexibleSUSY, NMSSMCALC and SARAH/SPheno shows almost exactly the same behaviour when varying s λκ .The NMSSMCALC M W prediction differs with the one of FlexibleSUSY (SARAH/SPheno) by at most 1.7 MeV (3.1 MeV) which is smaller than the SMuncertainty.The prediction obtained with NMSSMTools, however, seems to behave much flatter for large values of s λκ .This is likely because NMSSMTools seems to use the loop-corrected scalar masses in all parts of the M W calculation.Furthermore, we find that the M W predicted by NMSSMTools agrees much better with the other codes in Table 6 Comparison of the prediction for the SM-like Higgs boson mass, the W boson mass and the muon anomalous magnetic moment using FlexibleSUSY, NMSSMCALC, NMSSMTools, and SARAH/SPheno.The benchmark point P1, which was obtained from a scan in NMSSMCALC using the OS scheme, is now interpreted in the DR scheme (which is why the Higgs mass is no longer at 125 GeV)  Figs. 3 and 4, if we remove the M W -restoring fit function in Eq. (5.66) from its prediction.Thus, we suspect that the fit coefficients in Eq. (5.66) in NMSSMTools are outdated.Finally, in the MSSM limit s λκ → 0, we are also able to also compare with the code FeynHiggs 2.19.0 which calculates m h and M W in the MSSM rather than the NMSSM.For a detailed description of the M W calculation in FeynHiggs we refer to [63].In Fig. 4 we find that FeynHiggs yields the smallest M W prediction which is, however, still in good agreement with the other codes given the SM uncertainty alone.In particular the difference to the NMSSMCALC prediction is about 5.7 MeV.

CP-violating effects in the M W prediction
In this section we study the effect of CP-violating beyondthe-SM phases on the M W calculation.We consider the parameter point P1, which was initially defined in the CPconserving NMSSM, and study its dependence on the phases of A t , M 1 , M 2 , and M 3 .Note, that this investigation is for illustrative purpose and we hence do not check the validity of the phases w.r.t. the EDM constraints.Figure 5 (left) shows the resulting prediction for M W at two-loop O(α 2 new ) if the phases are varied individually for the parameter point P1.In the right plot we show the difference of the M W prediction to the SM prediction SM m h defined in Eq. (5.64b).We find that the phase of ϕ M 2 has the largest impact on M W for this parameter point of about 2 MeV while the overall SUSY corrections SM m h are at most 12 MeV which is due to the very light electroweakino masses.The stop quarks and gluinos have already been decoupled from the M W prediction as they are heavier than 1 TeV.Thus, also the phases ϕ A t and ϕ M 3 only impact M W at the sub-MeV level, which is in agreement with the findings in [102] for the MSSM.We furthermore observe that the phase dependence is dominating in the one-loop corrections while the two-loop corrections only lead to a change in the phase dependence at the sub-MeV level.In conclusion, the overall phase dependence is smaller than the size of the total shift of the SUSY corrections to M W .

Conclusions
In this paper, a consistent inclusion of the two-loop O(α t α s ) and O(α 2 new ) corrections to the ρ parameter and its applica- tion in the calculation of the W boson mass has been presented in the context of the complex NMSSM.Both corrections have been computed by our group in the previous computation of the loop-corrected Higgs boson masses at the corresponding orders.These two calculations use the gaugeless limit and the zero external momentum approximation.The renormalization features a mixed OS/DR scheme and conveniently allows to switch between OS and DR conditions used in the top/stop sector and for the charged Higgs boson mass.A scheme change in the top/stop sector is used to estimate the uncertainty in the ρ parameter prediction due to missing higher-orders.We showed that the O(α t α s ) and O(α 2 new ) corrections to the ρ parameter are significant and can help to reduce the theory uncertainty.After subtracting the SM corrections we add them back in the evaluation of the W mass including all known higher-order SM corrections in the OS scheme.We show that the effects arising from O(α t α s ) and O(α 2 new ) for M W are of the order of a few MeV which is smaller than the parametric uncertainty of the top mass and is of similar size as the missing higher-order SM corrections.We have implemented all corrections in the new version of the Fortran code NMSSMCALC which takes into account the most up-to-date higher-order corrections for the CP-violating NMSSM in the computation of the Higgs boson masses, Higgs boson decay widths and branching ratios, the muon anomalous magnetic moment a μ , electric dipole moments, and the W boson mass.Finally, we have performed a detailed comparison of the W boson mass, Higgs boson mass, and muon anomalous magnetic moment prediction between NMSSMCALC and the public spectrum generators FlexibleSUSY, SARAH/SPheno and NMSSMTools.We found good agreement between all calculations once different treatments of the most-important input parameters and renormalization scales are taken into account.

Fig. 4 Fig. 5
Fig. 4 Comparison between NMSSMCALC (red solid), NMSSMTools (blue dashed), SARAH/SPheno (magenta dash-dotted) and FlexibleSUSY (black dotted) for the parameter point BP3 as .12)where the first rotation matrix R G with one rotation angle β singles out the neutral Goldstone boson and the second rotation matrix R H rotates the five interaction states (h d , h u , h s , a, a s ) to the five mass eigenstates (h 1 , h 2 , h 3 , h 4 , h 5 ).•The charged Higgs interaction states h ± d , h ± u build up the charged Higgs bosons H ± with mass M H ± and the charged Goldstone bosons G ± .•The fermionic superpartners of the neutral Higgs bosons, the neutral higgsinos Hu , Hd and the singlino S, mix with the neutral gauginos B and W3 , resulting in five neutralino mass eigenstates denoted as χ 0 i , (i = 1, . . ., 5).The mass ordering of the χ 0 i is chosen as m χ 0 1 ≤ ... ≤ 123 (2023) 83:1079

Table 2
The parameter point P1: Higgs masses in GeV and the main components are shown

Table 3
Electroweakino masses in GeV for the parameter point P1 ters of P1 are kept fixed.The upper panel of Fig.1(left) shows the ρ parameter computed at one-loop (black), O(α t α s ) (blue) and O(α 2 new ) (red) as a function of t α s + (α λ + α κ + α t ) 2 , (5.62a) Upper left: the ρ parameter at full one-loop order (black), two-loop O(α t α s ) (blue) and two-loop O(α 2 new ) (red) as function of √ λ 2 + κ 2 .The dashed lines are for the DR scheme in the top/stop sector while the full lines are for the OS scheme.Lower left: renormal-ization scheme dependence ρ ren of the ρ parameter at full one-loop order (black), at two-loop O(α t α s ) (blue) and O(α 2 new ) (red).Right: the part of the Higgs-and EW-ino mass spectrum causing larger deviations

Table 4
The SM ρ parameter multiplied by a factor of 10 −3 , computed with M W = 80.36 GeV and m h = 125 GeV at one-loop and two-loop O(α t α s ) and O(α t α s + α 2 t )