aainstitutetext: Faculty of Fundamental Sciences, PHENIKAA University, Hanoi
12116, Vietnambbinstitutetext: Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germanyccinstitutetext: Institute for Theoretical
Physics, Karlsruhe Institute of Technology, Wolfgang-Gaede-Str. 1, 76131
Karlsruhe, Germany

#  h
h KA-TP-17-2023
h DESY-23-112
The ${\cal O}(\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2}$ Correction to
the $\rho$ Parameter and its Effect on the W Boson Mass Calculation in the
Complex NMSSM

Thi Nhung Dao b Martin Gabelmann c Margarete Mühlleitner
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present the prediction of the electroweak $\rho$ 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 $\rho$ parameter is
calculated at the full one-loop and leading and sub-leading two-loop order
$\mathcal{O}(\alpha+\alpha_{t}\alpha_{s}+\left(\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa}\right)^{2})$.
The new $\Delta\rho$ prediction is incorporated into a prediction of $M_{W}$
via a full supersymmetric (SUSY) one-loop calculation of $\Delta r$.
Furthermore, we include all known state-of-the-art SM higher-order corrections
to $\Delta r$. By comparing results for $\Delta\rho$ obtained using on-shell
(OS) and $\overline{\text{DR}}$ renormalization conditions in the top/stop
sector, we find that the scheme uncertainty is reduced at one-loop order by
55%, at two-loop $\mathcal{O}(\alpha_{s}\alpha_{t})$ by 22%, and at two-loop
$\mathcal{O}(\alpha_{t}+\alpha_{\kappa}+\alpha_{\lambda})^{2}$ by 16%,
respectively. The influence of the two-loop results on the $M_{W}$ 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.

##  1 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 Aad:2012tfa and CMS
Chatrchyan:2012ufa . 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_{W}^{\text{exp}}=80.379\pm 0.012\,\mathrm{GeV}$ Zyla:2020zbs . The SM
predicts the $W$ boson mass to be $M_{W}^{{\text{SM}},\,\text{OS}}=80.353\pm
0.004\,\mathrm{GeV}$ Bagnaschi:2022qhb in the on-shell (OS) renormalization
scheme using the state-of-the-art calculations Awramik:2003rn , and to be
$M_{W}^{{\text{SM}},{\overline{\text{MS}}}}=80.351\pm 0.003\,\mathrm{GeV}$ in
the ${\overline{\text{MS}}}$ scheme Degrassi:2014sxa ; Athron:2022isz 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\sigma$
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 CDF:2022hxs

$M_{W}^{\text{CDF}}=80.4335\pm 0.0094\,\mathrm{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\sigma$. Combining the new CDF result with the
other measurements from LEP, Tevatron and the LHC leads to a new world average
deBlas:2022hdk of

$M_{W}^{\text{exp}}=80.4133\pm 0.0080\,\mathrm{GeV},$ (1.2)

and a new deviation of order $6\sigma$. 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 $R-$parity 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 $\Delta r$ Sirlin:1980nh 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 $\Delta r$, the $W/Z$ masses $M_{W/Z}$, the Fermi constant $G_{F}$ and
the fine-structure constant $\alpha$ 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 $\Delta r$ at a given loop level.

Recent studies Athron:2022isz ; Bagnaschi:2022qhb ; Yang:2022gvz in the
Minimal Supersymmetric extension of the SM (MSSM) have shown that light and
compressed spectra of electroweakinos can yield $M_{W}^{\text{MSSM}}$ up to
$80.376\,\mathrm{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 the prediction
of $M_{W}$ within the Next-to-Minimal Supersymmetric Standard Model (NMSSM)
Fayet:1974pd ; Barbieri:1982eh ; Dine:1981rt ; Nilles:1982dy ; Frere:1983ag ;
Derendinger:1983bz ; Ellis:1988er ; Drees:1988fc ; Ellwanger:1993xa ;
Ellwanger:1995ru ; Ellwanger:1996gw ; Elliott:1994ht ; King:1995vk ;
Franke:1995tc ; Maniatis:2009re ; Ellwanger:2009dp . The model contains two
Higgs doublet- and an additional complex singlet-superfield. The $W$ boson
mass in this model has been first studied in Domingo:2011uf which includes the
full one-loop corrections to $\Delta r$. The authors of Stal:2015zca improved
the calculation by taking into account not only the one-loop but also the two-
loop corrections of ${\cal O}(\alpha\alpha_{s})$ Djouadi:1996pa ;
Djouadi:1998sq and ${\cal
O}(\alpha_{t}^{2},\alpha_{t}\alpha_{b},\alpha_{b}^{2})$ Heinemeyer:2002jq ;
Haestier:2005ja to the $\Delta\rho$ parameter into account. It is important to
stress that all contributions to $\Delta\rho$ beyond one-loop order have been
computed in the MSSM-limit so far. In this work, however, we explicitly
calculate $\Delta\rho$ including the full dependence on all NMSSM-specific
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 Ender:2011qh ;
Graf:2012hh , and the two-loop ${\cal O}(\alpha_{t}\alpha_{s})$
Muhlleitner:2014vsa , ${\cal O}(\alpha_{t}^{2})$ Dao:2019qaz as well as ${\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})$ Dao:2021khm
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 ${\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})$, which contain specific
NMSSM corrections related to the NMSSM superpotential parameters $\lambda$ and
$\kappa$. The $\rho$ parameter and $W$ boson mass calculations are implemented
in the Fortran code NMSSMCALC Baglio:2013iia ; Ender:2011qh ; Graf:2012hh ;
Muhlleitner:2014vsa ; Dao:2019qaz ; Dao:2021khm which is publicly available.
The program computes the Higgs boson masses and mixing angles up to two-loop
${\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2}+\alpha_{t}\alpha_{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
Nhung:2013lpa and at two-loop ${\cal O}(\alpha_{t}\alpha_{s})$
Muhlleitner:2015dua and two-loop ${\cal O}(\alpha_{t}^{2})$
Borschensky:2022pfc . For the CP-violating NMSSM the calculation of the
electric dipole moments (EDMs) has been implemented King:2015oxa to be checked
against the experimental constraints. Recently, the electron and muon
anomalous magnetic moments have been included in the code Dao:2021vqp . The
code has been extended to include the electroweak corrections to the NMSSM
Higgs decays in NMSSMCALCEW Baglio:2019nlc ; Dao:2019nxi ; Dao:2020dfb .

This paper is organized as follows. In Section 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 two-loop corrections to $\Delta\rho$ and the loop-corrected $W$
boson mass will be presented in Section 3 and Section 4, respectively.
Furthermore, we discuss the renormalization schemes used to obtain UV-finite
$\Delta\rho$ and $\Delta 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 Section 6.

##  2 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 Graf:2012hh ; Baglio:2013iia ; Muhlleitner:2014vsa ;
Dao:2019qaz ; Dao:2021khm . The superpotential of the complex NMSSM is given
by ($i,j=1,2$)

$\displaystyle\mathcal{W}_{{\text{NMSSM}}}=\epsilon_{ij}[y_{e}\hat{H}^{i}_{d}\hat{L}^{j}\hat{E}^{c}+y_{d}\hat{H}_{d}^{i}\hat{Q}^{j}\hat{D}^{c}-y_{u}\hat{H}^{i}_{u}\hat{Q}^{j}\hat{U}^{c}]-\epsilon_{ij}\lambda\hat{S}\hat{H}_{d}^{i}\hat{H}_{u}^{j}+\frac{1}{3}\kappa{\hat{S}}^{3}\;,$
(2.3)

with the quark and leptonic superfields $\hat{Q}$, $\hat{U}$, $\hat{D}$,
$\hat{L}$, $\hat{E}$, the Higgs doublet superfields $\hat{H}_{d}$,
$\hat{H}_{u}$, the singlet superfield $\hat{S}$ and the totally antisymmetric
tensor $\epsilon_{12}=\epsilon^{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$\times$3 matrices in flavour space. The parameters
$\lambda$ and $\kappa$ are in general complex. The soft SUSY-breaking
Lagrangian reads

$\displaystyle\mathcal{L}_{\rm soft,NMSSM}=$ $\displaystyle-
m_{H_{d}}^{2}H_{d}^{\dagger}H_{d}-m_{H_{u}}^{2}H_{u}^{\dagger}H_{u}-m_{\tilde{Q}}^{2}\tilde{Q}^{\dagger}\tilde{Q}-m_{\tilde{L}}^{2}\tilde{L}^{\dagger}\tilde{L}-m_{\tilde{u}_{R}}^{2}\tilde{u}_{R}^{*}\tilde{u}_{R}-m_{\tilde{d}_{R}}^{2}\tilde{d}_{R}^{*}\tilde{d}_{R}$
$\displaystyle-
m_{\tilde{e}_{R}}^{2}\tilde{e}_{R}^{*}\tilde{e}_{R}-(\epsilon_{ij}[y_{e}A_{e}H_{d}^{i}\tilde{L}^{j}\tilde{e}_{R}^{*}+y_{d}A_{d}H_{d}^{i}\tilde{Q}^{j}\tilde{d}_{R}^{*}-y_{u}A_{u}H_{u}^{i}\tilde{Q}^{j}\tilde{u}_{R}^{*}]+\mathrm{h.c.})$
$\displaystyle-\frac{1}{2}(M_{1}\tilde{B}\tilde{B}+M_{2}\tilde{W}_{j}\tilde{W}_{j}+M_{3}\tilde{G}\tilde{G}+\mathrm{h.c.})$
(2.4) $\displaystyle-m_{S}^{2}|S|^{2}+(\epsilon_{ij}\lambda
A_{\lambda}SH_{d}^{i}H_{u}^{j}-\frac{1}{3}\kappa
A_{\kappa}S^{3}+\mathrm{h.c.})\;,$

where $H_{u,d}$, $S$, $\tilde{Q}$, $\tilde{L}$ and $\tilde{x}_{R}$ ($x=e,u,d$)
are the scalar components of the respective superfields. The fermionic
superpartners $\tilde{B}$ (bino), $\tilde{W}_{1,2,3}$ (wino) and $\tilde{G}$
(gluino) of the $U(1)_{Y}$, $SU(2)_{L}$ and $SU(3)_{c}$ gauge bosons obtain
the soft SUSY-breaking gaugino masses $M_{1}$, $M_{2}$ and $M_{3}$,
respectively. The mass parameters $M_{1,2,3}$ and the soft SUSY-breaking
trilinear couplings $A_{x}$ ($x=\lambda,\kappa,u,d,e$) are in general complex,
while the mass parameters of the scalar fields are real.

After electroweak symmetry breaking, the Higgs boson fields are expanded
around their vacuum expectation values (VEVs) $v_{u}$, $v_{d}$, and $v_{s}$,
respectively,

$H_{d}=\begin{pmatrix}\frac{v_{d}+h_{d}+ia_{d}}{\sqrt{2}}\\\
h_{d}^{-}\end{pmatrix},\,\,H_{u}=e^{i\varphi_{u}}\begin{pmatrix}h_{u}^{+}\\\
\frac{v_{u}+h_{u}+ia_{u}}{\sqrt{2}}\end{pmatrix},\,\,S=\frac{e^{i\varphi_{s}}}{\sqrt{2}}(v_{s}+h_{s}+ia_{s})\,,$
(2.5)

with the CP-violating phases $\varphi_{u,s}$. The relation to the SM VEV
$v\approx 246.22$ GeV is given by

$v^{2}=v_{u}^{2}+v_{d}^{2}\,,$ (2.6)

and we define the mixing angle $\tan\beta$ as

$\tan\beta=\frac{v_{u}}{v_{d}}\;.$ (2.7)

Thus, the expressions for the tree-level weak gauge boson masses are the same
as in the SM,

$M_{W}^{2}=\frac{1}{4}g_{2}^{2}(v_{u}^{2}+v_{d}^{2}),\quad
M_{Z}^{2}=\frac{1}{4}(g_{1}^{2}+g_{2}^{2})(v_{u}^{2}+v_{d}^{2}),$ (2.8)

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 $\theta_{W}$ as

$g_{1}=\frac{e}{c_{\theta_{W}}},\quad g_{2}=\frac{e}{s_{\theta_{W}}},$ (2.9)

where the short hand notation
$c_{x}\equiv\cos(x),s_{x}\equiv\sin(x),t_{x}\equiv\tan(x)$ is used in this
paper. The effective $\mu$ parameter is given by

$\mu_{\text{eff}}=\frac{\lambda v_{s}e^{i\varphi_{s}}}{\sqrt{2}}.$ (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}}\leq\dots\leq 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,

$\displaystyle(h_{d},h_{u},h_{s},a,a_{s},G^{0})^{T}$ $\displaystyle=$
$\displaystyle\mathcal{R}^{G}(\beta)\,(h_{d},h_{u},h_{s},a_{d},a_{u},a_{s})^{T},$
(2.11) $\displaystyle(h_{1},h_{2},h_{3},h_{4},h_{5},G^{0})^{T}$
$\displaystyle=$
$\displaystyle\mathcal{R}^{H}\,(h_{d},h_{u},h_{s},a,a_{s},G^{0})^{T},$ (2.12)

where the first rotation matrix $\mathcal{R}^{G}$ with one rotation angle
$\beta$ singles out the neutral Goldstone boson and the second rotation matrix
$\mathcal{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}^{\pm},h_{u}^{\pm}$ build up the
charged Higgs bosons $H^{\pm}$ with mass $M_{H^{\pm}}$ and the charged
Goldstone bosons $G^{\pm}$.

  * •

The fermionic superpartners of the neutral Higgs bosons, the neutral higgsinos
$\tilde{H}_{u}$, $\tilde{H}_{d}$ and the singlino $\tilde{S}$, mix with the
neutral gauginos $\tilde{B}$ and $\tilde{W}_{3}$, resulting in five neutralino
mass eigenstates denoted as $\tilde{\chi}^{0}_{i}$, $(i=1,...,5)$. The mass
ordering of the $\tilde{\chi}^{0}_{i}$ is chosen as
$m_{\tilde{\chi}^{0}_{1}}\leq...\leq m_{\tilde{\chi}^{0}_{5}}$, and the
rotation matrix $N$ transforms the fields $(\tilde{B},\ \tilde{W}_{3},\
\tilde{H}_{d},\ \tilde{H}_{u},\tilde{S})^{T}$ into the mass eigenstates.

  * •

The two chargino mass eigenstates,

$\displaystyle\tilde{\chi}_{i}^{+}=\left(\begin{array}[]{c}\tilde{\chi}_{L_{i}}^{+}\\\\[2.84544pt]
{\tilde{\chi}_{R_{i}}^{-^{*}}}\end{array}\right)\;,\quad i=1,2\;,$ (2.15)

are obtained from the rotation of the interaction states, given by the charged
Higgsinos $\tilde{H}^{\pm}_{d}$, $\tilde{H}^{\pm}_{u}$ and the charged
gauginos $\tilde{W}^{\pm}=\left(\tilde{W}_{1}\mp
i\tilde{W}_{2}\right)/\sqrt{2}$, to the mass eigenstates. This is achieved by
using a bi-unitary transformation with the two $2\times 2$ unitary matrices
$V^{\chi}$ and $U^{\chi}$,

$\displaystyle\tilde{\chi}_{L}^{+}=V^{\chi}(\tilde{W}^{+},\tilde{H}^{+}_{u})^{T},\
\ \ \tilde{\chi}_{R}^{-}=U^{\chi}(\tilde{W}^{-},\tilde{H}^{-}_{d})^{T}.$
(2.16)

  * •

The scalar partners of the left- and right-handed up-type quarks are denoted
by $\tilde{u}_{L/R}^{i}$, of the down-type quarks by $\tilde{d}^{i}_{L/R}$,
and of the charged leptons by $\tilde{l}^{i}_{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^{\tilde{f}}$.

  * •

There are three scalar partners of the left-handed neutrinos, denoted as
$\tilde{\nu}_{i}$ ($i=1,2,3$) with the masses

$m^{2}_{\tilde{\nu}_{i}}=\frac{1}{2}M_{Z}^{2}c_{2\beta}+m_{\tilde{L}_{i}}^{2}\;.$
(2.17)

For detailed discussions and derivations of the tree-level spectrum we refer
to earlier works Graf:2012hh ; Baglio:2013iia ; Muhlleitner:2014vsa ;
Dao:2019qaz ; Dao:2021khm .

##  3 The One- and Two-Loop Corrections to the $\rho$ Parameter

In this section we describe the calculation of the full one-loop and the
${{\cal O}(\alpha_{t}\alpha_{s})}$ and ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$ two-loop corrections
to the $\rho$ parameter. The $\rho$ parameter is defined by the neutral- to
charged-current coupling ratio at zero external momentum transfer Ross:1975fq
; Veltman:1977kh ; vanderBij:1986hy

$\rho=\frac{G_{NC}(0)}{G_{CC}(0)},$ (3.18)

where the neutral current coupling $G_{NC}(0)$ can be determined from, for
example, the coefficient of the $\nu_{e}\bar{\nu}_{e}\to\nu_{e}\bar{\nu}_{e}$
scattering amplitude and the charged current coupling $G_{CC}(0)$ from the
$e\bar{\nu}_{e}\to e\bar{\nu}_{e}$ scattering amplitude. In the NMSSM, at tree
level we have

$G_{NC}(0)=\frac{g_{2}^{2}}{2M_{Z}^{2}c_{W}^{2}},\quad
G_{CC}(0)=\frac{g_{2}^{2}}{2M_{W}^{2}},$ (3.19)

therefore $\rho^{(\text{tree-level})}=1$, as in the SM. Higher-order
corrections to the $\rho$ parameter lead to a deviation from unity, which can
be written as vanderBij:1986hy ; Consoli:1989fg ; Fleischer:1993ub

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{1}{1-\Delta\rho}=1+\Delta^{\text{\tiny(}\\!\text{\tiny
1}\\!\text{\tiny)}}\\!\rho+((\Delta^{\text{\tiny(}\\!\text{\tiny
1}\\!\text{\tiny)}}\\!\rho)^{2}+\Delta^{\text{\tiny(}\\!\text{\tiny
2}\\!\text{\tiny)}}\\!\rho))+\cdots,$ (3.20) $\displaystyle\Delta\rho$
$\displaystyle=$ $\displaystyle\Delta^{\text{\tiny(}\\!\text{\tiny
1}\\!\text{\tiny)}}\\!\rho+\Delta^{\text{\tiny(}\\!\text{\tiny
2}\\!\text{\tiny)}}\\!\rho+\cdots,$ (3.21)

where the superscript $(n)$ indicates that the calculation is performed at the
$n$-loop order. $\Delta^{(n)}\rho$ 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

$\displaystyle\Delta^{\text{\tiny(}\\!\text{\tiny 1}\\!\text{\tiny)}}\\!\rho$
$\displaystyle=$
$\displaystyle\frac{\Sigma_{ZZ}^{(1),T}}{M_{Z}^{2}}-\frac{\Sigma_{WW}^{(1),T}}{M_{W}^{2}},$
(3.22) $\displaystyle\Delta^{\text{\tiny(}\\!\text{\tiny
2}\\!\text{\tiny)}}\\!\rho$ $\displaystyle=$
$\displaystyle-\frac{\Sigma_{ZZ}^{(1),T}}{M_{Z}^{2}}\left(\frac{\Sigma_{ZZ}^{(1),T}}{M_{Z}^{2}}-\frac{\Sigma_{WW}^{(1),T}}{M_{W}^{2}}\right)+\left(\frac{\Sigma_{ZZ}^{(2),T}}{M_{Z}^{2}}-\frac{\Sigma_{WW}^{(2),T}}{M_{W}^{2}}\right),$
(3.23)

where the transverse part of the gauge boson self-energy $\Sigma_{VV}^{(n),T}$
($V\equiv W,Z$) is evaluated at zero external momentum. These expressions are
the same for the SM and for the 2-Higgs Doublet Model (2HDM) presented in
Hessenberger:2016atw and hold in general for every theory with $\rho=1$ at
tree level. In this paper, we compute the gauge boson self-energies at one-
and 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 NMSSM-specific
contributions. This strategy can be followed when considering $\Delta\rho$. In
the investigation of the $W$ boson mass it is, however, crucial to include
existing SM-like higher-order corrections to $\Delta r$ beyond the two-loop
level (cf. Section 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 ($\overline{\text{DR}}$) rather than
minimal subtraction (${\overline{\text{MS}}}$). While the results for the
individual gauge boson self-energies differ when using
${\overline{\text{MS}}}$ or $\overline{\text{DR}}$, we find that the
difference of the $W$ and $Z$ self-energies, entering $\Delta\rho$, 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 $\Delta^{\text{\tiny(}\\!\text{\tiny
1}\\!\text{\tiny)}}\\!\rho$ into four contributions given by the SM fermion,
the Higgs, the squark and slepton, and the chargino and neutralino ones,

$\Delta^{\text{\tiny(}\\!\text{\tiny
1}\\!\text{\tiny)}}\\!\rho=\Delta^{f}\rho+\Delta^{H}\rho+\Delta^{\tilde{f}}\rho+\Delta^{\tilde{\chi}}\rho\,.$
(3.24)

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

$\Delta^{f}\rho=\frac{m_{\tau}^{2}+3m_{b}^{2}+3m_{t}^{2}}{16\pi^{2}v^{2}}+\frac{3m_{b}^{2}m_{t}^{2}\ln\frac{m_{t}^{2}}{m_{b}^{2}}}{8\pi^{2}v^{2}(m_{b}^{2}-m_{t}^{2})}\,.$
(3.25)

The contribution from the Higgs sector, i.e. from $G^{\pm},G^{0},h^{i}$ and
$H^{\pm}$, reads

$\displaystyle\Delta^{H}\rho=$
$\displaystyle\frac{1}{8v^{2}\pi^{2}}\bigg{[}s_{2\theta_{W}}^{2}A_{0}(M_{H^{\pm}}^{2})-2c_{2\theta_{W}}^{2}B_{00}(M_{H^{\pm}}^{2},M_{H^{\pm}}^{2})-\sum_{i,j=1}^{5}4g_{h_{i}h_{j}Z}^{2}B_{00}(m_{h_{i}}^{2},m_{h_{j}}^{2})$
$\displaystyle+\sum_{i=1}^{5}\left(8\lvert
g_{h_{i}H^{-}W^{+}}\rvert^{2}B_{00}(m_{h_{i}}^{2},M_{H^{\pm}}^{2})+8\lvert
g_{h_{i}G^{-}W^{+}}\rvert^{2}B_{00}(m_{h_{i}}^{2},0)\right)\bigg{]}\,,$ (3.26)

where the one-loop functions are given at the end of this section and

$\displaystyle g_{h_{i}H^{-}W^{+}}$ $\displaystyle=\frac{1}{2}(-{\cal
R}^{H}_{i4}-i(s_{\beta}{\cal R}^{H}_{i1}-c_{\beta}{\cal R}^{H}_{i2}))\,,\quad
g_{h_{i}G^{-}W^{+}}=\frac{1}{2}i(c_{\beta}{\cal R}^{H}_{i1}+s_{\beta}{\cal
R}^{H}_{i2})\,,\quad$ $\displaystyle g_{h_{i}h_{j}Z}$
$\displaystyle=\frac{1}{2}\left({\cal R}^{H}_{i4}(s_{\beta}{\cal
R}^{H}_{j1}-c_{\beta}{\cal R}^{H}_{j2})-{\cal R}^{H}_{j4}(s_{\beta}{\cal
R}^{H}_{i1}-c_{\beta}{\cal R}^{H}_{i2})\right)\;,$ (3.27)

with $c_{2W}\equiv\cos 2\theta_{W},s_{2W}\equiv\sin 2\theta_{W}$.

The contribution from the third generation squarks and sleptons is given by

$\displaystyle\Delta^{\tilde{f}}\rho=$
$\displaystyle\frac{1}{8v^{2}\pi^{2}}\bigg{[}-2\sum_{i=1,2}\sum_{\tilde{f}=\tilde{\tau},\tilde{t},\tilde{b}}n_{\tilde{f}}Q_{\tilde{f}_{i}}^{2}s_{\theta_{W}}^{2}\left[n_{\tilde{f}}(n_{\tilde{f}}\lvert
Q_{\tilde{f}_{i}}\rvert-2t_{\tilde{f}_{i}})\lvert
U_{i1}^{\tilde{f}}\rvert^{2}-2s_{\theta_{W}}^{2}\right]A_{0}(m_{\tilde{f}_{i}}^{2})$
$\displaystyle-2B_{00}(m_{\tilde{\nu}_{3}}^{2},m_{\tilde{\nu}_{3}}^{2})-2\sum_{i=1,2}\left[\lvert
U_{i1}^{\tilde{\tau}}\rvert^{2}-2s_{\theta_{W}}^{2}\right]^{2}B_{00}(m_{\tilde{\tau}_{i}}^{2},m_{\tilde{\tau}_{j}}^{2})$
$\displaystyle-\frac{2}{3}\sum_{i=1,2}\sum_{\tilde{f}=\tilde{t},\tilde{b}}n_{\tilde{f}}^{2}\left[\lvert
U_{i1}^{\tilde{f}}\rvert^{2}-2\lvert Q_{\tilde{f}_{i}}\rvert
s_{\theta_{W}}^{2}\right]^{2}B_{00}(m_{\tilde{f}_{i}}^{2},m_{\tilde{f}_{j}}^{2})$
$\displaystyle+4\sum_{i=1,2}\lvert
U_{i1}^{\tilde{\tau}}\rvert^{2}B_{00}(m_{\tilde{\tau}_{i}}^{2},m_{\tilde{\nu}_{3}}^{2})+4n_{\tilde{f}}\sum_{i,j=1,2}\lvert
U_{i1}^{\tilde{t}}\rvert^{2}\lvert
U_{i1}^{\tilde{b}}\rvert^{2}B_{00}(m_{\tilde{b}_{i}}^{2},m_{\tilde{t}_{j}}^{2})\bigg{]}\,,$
(3.28)

where $n_{\tilde{f}}=3$ for squarks and $n_{\tilde{f}}=1$ for sleptons,
$t_{\tilde{\tau}_{i}/\tilde{b}_{i}}=-1/2,t_{\tilde{t}_{i}}=1/2$,
$Q_{\tilde{f}_{i}}$ denotes the electric charge of the sfermions. The
contribution from the charginos and neutralinos can be cast into the form

$\displaystyle\Delta^{\tilde{\chi}}\rho=$
$\displaystyle\frac{1}{2v^{2}\pi^{2}}\bigg{\\{}\sum_{i,j=1}^{5}\bigg{[}\text{Re}\,{(g_{\tilde{\chi}^{0}_{i}\tilde{\chi}^{0}_{j}Z}^{L}g_{\tilde{\chi}^{0}_{j}\tilde{\chi}^{0}_{i}Z}^{R})}m_{\tilde{\chi}_{i}^{0}}m_{\tilde{\chi}_{j}^{0}}B_{0}(m^{2}_{\tilde{\chi}_{i}^{0}},m^{2}_{\tilde{\chi}_{j}^{0}})$
$\displaystyle-\frac{1}{2}(\lvert
g_{\tilde{\chi}^{0}_{i}\tilde{\chi}^{0}_{j}Z}^{L}\rvert^{2}+\lvert
g_{\tilde{\chi}^{0}_{i}\tilde{\chi}^{0}_{j}Z}^{R}\rvert^{2}){\cal
F}(m^{2}_{\tilde{\chi}_{i}^{0}},m^{2}_{\tilde{\chi}_{j}^{0}})\bigg{]}$
$\displaystyle+\sum_{i=1}^{2}\sum_{j=1}^{5}\,\bigg{[}-2\text{Re}\,{(g^{L}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{0}W}(g^{R}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{0}W})^{*})}m_{\tilde{\chi}_{i}^{\pm}}m_{\tilde{\chi}_{j}^{0}}B_{0}(m^{2}_{\tilde{\chi}_{i}^{\pm}},m^{2}_{\tilde{\chi}_{j}^{0}})$
$\displaystyle+(\lvert
g^{L}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{0}W}\rvert^{2}+\lvert
g^{R}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{0}W}\rvert^{2}){\cal
F}(m^{2}_{\tilde{\chi}_{i}^{\pm}},m^{2}_{\tilde{\chi}_{j}^{0}})\bigg{]}$
$\displaystyle+\sum_{i,j=1}^{2}\,\bigg{[}2\text{Re}\,{(g^{L}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{\pm}Z}(g^{R}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{\pm}Z})^{*})}m_{\tilde{\chi}_{i}^{\pm}}m_{\tilde{\chi}_{j}^{\pm}}B_{0}(m^{2}_{\tilde{\chi}_{i}^{\pm}},m^{2}_{\tilde{\chi}_{j}^{\pm}})$
$\displaystyle+(\lvert
g^{L}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{\pm}Z}\rvert^{2}+\lvert
g^{R}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{\pm}Z}\rvert^{2}){\cal
F}(m^{2}_{\tilde{\chi}_{i}^{\pm}},m^{2}_{\tilde{\chi}_{j}^{\pm}})\bigg{]}\bigg{\\}}\,,$
(3.29)

where

$\displaystyle g_{\tilde{\chi}^{0}_{i}\tilde{\chi}^{0}_{j}Z}^{L}$
$\displaystyle=-\frac{1}{2}(N_{i3}N_{j3}^{*}-N_{i4}N_{j4}^{*}),\quad$
$\displaystyle g_{\tilde{\chi}^{0}_{i}\tilde{\chi}^{0}_{j}Z}^{R}$
$\displaystyle=-g_{\tilde{\chi}^{0}_{j}\tilde{\chi}^{0}_{i}Z}^{L},$
$\displaystyle g^{L}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{\pm}Z}$
$\displaystyle=c_{\theta_{W}}^{2}V^{\chi}_{i1}V^{\chi*}_{j1}+\frac{1}{2}c_{2\theta_{W}}V^{\chi}_{i2}V^{\chi*}_{j2},\quad$
$\displaystyle g^{R}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{\pm}Z}$
$\displaystyle=c_{\theta_{W}}^{2}U^{\chi}_{j1}U^{\chi*}_{i1}+\frac{1}{2}c_{2\theta_{W}}U^{\chi}_{j2}U^{\chi*}_{i2},$
$\displaystyle g^{L}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{0}W}$
$\displaystyle=V^{\chi}_{i1}N_{j2}^{*}+\frac{1}{\sqrt{2}}V^{\chi}_{i2}N_{j3}^{*},\quad$
$\displaystyle g^{R}_{\tilde{\chi}_{i}^{\pm}\tilde{\chi}_{j}^{0}W}$
$\displaystyle=U^{\chi}_{i1}N_{j2}-\frac{1}{\sqrt{2}}U^{\chi}_{i2}N_{j4},$
$\displaystyle{\cal F}(x,y)$
$\displaystyle=A_{0}(y)-2B_{00}(x,y)+yB_{0}(x,y).$ (3.30)

In the above expressions, the one-loop integrals are defined as

$\displaystyle A_{0}(x)$ $\displaystyle=x(1-\ln(\bar{x}))$ $\displaystyle
B_{0}(x,y)$
$\displaystyle=\begin{cases}(1+\frac{-x\ln(\bar{x})+y\ln(\bar{y})}{x-y}),\quad&x\neq
y\\\ -\ln(\bar{x}),\quad&x=y\end{cases}$ $\displaystyle B_{00}(x,y)$
$\displaystyle=\begin{cases}\frac{3}{8}(x+y)+\frac{-x^{2}\ln(\bar{x})+y^{2}\ln(\bar{y})}{4(x-y)},\quad&x\neq
y\\\ \frac{x}{2}(1-\ln(\bar{x})),\quad&x=y\end{cases},$ (3.31)

where the bar symbol indicates the dimensionless quantities, i.e.
$\bar{x}=x/\mu_{R}^{2}$ etc. with $\mu_{R}$ being the renormalization scale.

The two-loop ${{\cal O}(\alpha_{t}\alpha_{s})}$ corrections
For the ${{\cal O}(\alpha_{t}\alpha_{s})}$ corrections, the first term in Eq.
(3.23) will vanish and hence

$\displaystyle\Delta^{\alpha_{t}\alpha_{s}}\rho$ $\displaystyle=$
$\displaystyle\frac{\Sigma_{ZZ}^{\alpha_{t}\alpha_{s},T}(0)}{M_{Z}^{2}}-\frac{\Sigma_{WW}^{\alpha_{t}\alpha_{s},T}(0)}{M_{W}^{2}}.$
(3.32)

We take the results of $\Sigma_{ZZ}^{\alpha_{t}\alpha_{s},T}(0)$ and
$\Sigma_{WW}^{\alpha_{t}\alpha_{s},T}(0)$ implemented in NMSSMCALC, which have
been computed by our group in Ref. Muhlleitner:2015dua in the complex NMSSM.
For the detailed calculation we refer the reader to Muhlleitner:2015dua , 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,

$e,M_{W},M_{Z}\to 0,\,\,\frac{M_{W}}{M_{Z}}=\text{const.}$ (3.33)

The ratio $\frac{\Sigma_{VV}^{\alpha_{t}\alpha_{s},T}(0)}{M_{V}^{2}}$
($V=W,Z$), however, is non-zero and proportional to $\alpha_{t}\alpha_{s}$. In
our calculation, we set the bottom quark mass to zero. The two-loop ${{\cal
O}(\alpha_{t}\alpha_{s})}$ corrections can be decomposed into the
contributions from the genuine two-loop diagrams (containing either a
gluon/gluino loop or a loop with a stop/sbottom quartic coupling) and the
contributions from the counterterm inserted one-loop diagrams (containing
either coupling-type counterterms or propagator-type counterterms). The set of
independent parameters entering the top/stop and bottom/sbottom sector, that
need to be renormalized at ${\cal O}(\alpha_{s})$ are

$m_{t},\;m_{\tilde{Q}_{3}},\;m_{\tilde{t}_{R}}\quad\mbox{and}\quad A_{t}\;.$
(3.34)

In Muhlleitner:2015dua we have discussed two renormalization schemes for these
parameters, on-shell (OS) and $\overline{\text{DR}}$ renormalization. We keep
these two options of renormalization schemes here, too. The expressions for
the required counterterms were presented in Muhlleitner:2015dua . We have
explicitly confirmed that the obtained result for
$\Delta^{\alpha_{t}\alpha_{s}}_{{\text{NMSSM}}}\rho$ is UV finite. For the SM-
like contributions we reproduce the known result in the literature
Djouadi:1987gn ; Djouadi:1987di ,

$\Delta^{\alpha_{t}\alpha_{s}}_{{\text{SM}}}\rho=-(1+\pi^{2}/3)\frac{\alpha_{s}m_{t}^{2}}{(8\pi^{3}v^{2})},$
(3.35)

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 Djouadi:1987gn ;
Djouadi:1987di were obtained using dimensional regularization. It should be
stressed that $\Delta^{\alpha_{t}\alpha_{s}}\rho$ in the NMSSM is identical to
the corresponding quantity in the MSSM which has been calculated in
Djouadi:1996pa ; Djouadi:1998sq .

The two-loop ${{\cal O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$
corrections
For the sake of a convenient notation we denote

$\alpha_{\text{EW}}^{2}\equiv(\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2}.$
(3.36)

Since the electroweak sector contributes to the gauge boson self-energies at
one-loop level, its contribution to $\Delta\rho$ at two-loop level contains
also the one-loop squared terms, hence

$\displaystyle\Delta^{\alpha_{\text{EW}}^{2}}\rho$ $\displaystyle=$
$\displaystyle-\frac{\Sigma_{ZZ}^{(1),T}}{M_{Z}^{2}}\left(\frac{\Sigma_{ZZ}^{(1),T}}{M_{Z}^{2}}-\frac{\Sigma_{WW}^{(1),T}}{M_{W}^{2}}\right)+\left(\frac{\Sigma_{ZZ}^{(2),T}}{M_{Z}^{2}}-\frac{\Sigma_{WW}^{(2),T}}{M_{W}^{2}}\right).$
(3.37)

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 Dao:2021khm
and implemented in NMSSMCALC. The results have been obtained in the gaugeless
limit, cf. Eq. 3.33. In this limit the Higgs-Goldstone couplings are non-zero
at $\order{\alpha_{\lambda}+\alpha_{\kappa}}$ while the Goldstone tree-level
masses vanish which leads to intermediate infra-red (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 $\Sigma_{VV}^{(2),T}$
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
self-energies $\Sigma_{VV}^{(2),T}$. For the parameters of the
Higgsino/singlino and Higgs sectors we apply a mixed $\overline{\text{DR}}$-OS
renormalisation scheme while for the top/stop sector we apply an OS scheme or
a $\overline{\text{DR}}$ scheme. All counterterms for the complex phases
$\varphi_{\alpha}$ ($\alpha=s,u,\kappa,\lambda$) can be set to zero. The
remaining input parameters together with the applied renormalization
conditions are given by

$\displaystyle\underbrace{M_{H^{\pm}}^{2},v,s_{\theta_{W}},}_{\mbox{OS
scheme}}\underbrace{m_{t},m_{\tilde{Q}_{3}},m_{\tilde{t}_{R}},A_{t}}_{\mbox{OS}/\overline{\mbox{DR}}\mbox{
scheme}},\underbrace{\tan\beta,|\lambda|,v_{s},|\kappa|,\mbox{Re}A_{\kappa}}_{\overline{\mbox{DR}}\mbox{
scheme}}\,,$ (3.38)

in case $M_{H^{\pm}}^{2}$ is used as independent input, or

$\displaystyle\underbrace{v,s_{\theta_{W}},}_{\mbox{OS
scheme}}\underbrace{m_{t},m_{\tilde{Q}_{3}},m_{\tilde{t}_{R}},A_{t}}_{\mbox{OS}/\overline{\mbox{DR}}\mbox{
scheme}},\underbrace{\tan\beta,|\lambda|,v_{s},|\kappa|,\mbox{Re}A_{\lambda},\mbox{Re}A_{\kappa}}_{\overline{\mbox{DR}}\mbox{
scheme}}\,,$ (3.39)

if $\mbox{Re}A_{\lambda}$ is chosen as independent input rather than
$M_{H^{\pm}}$. The tadpole counterterms are defined such that the minima of
the Higgs potential do not change at higher order. We explicitly confirmed
that $\Delta^{\alpha_{\text{EW}}^{2}}\rho$ is UV finite and is free of IR
divergences (cf. Ref. Dao:2021khm for a detailed discussion on the
cancellation of all IR divergences).

The corresponding result $\Delta^{\alpha_{\text{EW}}^{2}}_{\text{SM}}\rho$
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\theta_{W}$
and OS/$\overline{\text{DR}}$ conditions for $m_{t}$. In the SM only
contributions from top and Higgs diagrams enter at
$\order{\alpha_{\text{EW}}^{2}}$, hence
$\Delta^{\alpha_{\text{EW}}^{2}}_{\text{SM}}\rho=\Delta^{\alpha_{t}^{2}}_{\text{SM}}\rho$.
It should be noted that $\Delta^{\alpha_{\text{EW}}^{2}}_{\text{SM}}\rho$ does
not depend on the renormalization of the Higgs mass and the Higgs tadpole.
Note that our expression of $\Delta^{\alpha_{t}^{2}}_{\text{SM}}\rho$ in Eq.
(3.37) looks different from the one computed by Fleischer, Tarasov and
Jegerlehner in Fleischer:1993ub for the following reasons. In their
computation, $\Delta\rho$ at two-loop order is defined as

$\Delta^{\alpha_{t}^{2}}_{\text{SM,
FTJ}}\rho=\left(\frac{\Sigma_{ZZ}^{\alpha_{t}^{2},T}(0)}{M_{Z}^{2}}-\frac{\Sigma_{WW}^{\alpha_{t}^{2},T}(0)}{M_{W}^{2}}\right),$
(3.40)

where $\Sigma_{VV}^{\alpha_{t}^{2},T}(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.23) are not contributing. This is because they
used the Fermi constant $G_{\mu}$ as an input parameter. Therefore, in
Fleischer:1993ub 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 $\alpha(0)$) as an input. Considering the one-
loop expression for $\Delta\rho$ 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 $\delta v$ which is not present in the result
of Eq. (3.40). The UV-finite part of $\delta 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-
loop-squared term in Eq. (3.37). Our final results for
$\Delta^{\alpha_{t}^{2}}_{\text{SM}}\rho$ are the same as the results of
Fleischer:1993ub provided that $v=1/\sqrt{\sqrt{2}G_{\mu}}$. The Higgs mass in
the SM result $\Delta^{\alpha_{\text{EW}}^{2}}_{\text{SM}}\rho$ is set to be
equal to the possibly non-zero tree-level SM-like Higgs boson mass in
$\Delta^{\alpha_{\text{EW}}^{2}}_{\text{NMSSM}}\rho$ 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.

##  4 Calculation of the $W$ Boson Mass in the OS Scheme

The $W$ boson mass can be computed from the following relation Sirlin:1980nh ,

$\frac{G_{\mu}}{\sqrt{2}}=\frac{\pi\alpha}{2M_{W}^{2}s_{W}^{2}}(1+\Delta r),$
(4.41)

where $G_{\mu}$ is the Fermi constant, $\alpha\equiv\alpha(0)=e^{2}/4\pi$ is
the fine-structure constant in the Thomson limit, and $\Delta r$ includes all
loop contributions to the amplitude of the $\mu\to e\bar{\nu}_{e}\nu_{\mu}$
decay after subtracting the Fermi-model-type QED correction. By using the OS
weak mixing angle $s_{W}^{2}=1-M_{W}^{2}/M_{Z}^{2}$, Eq. 4.41 can be solved
for $M_{W}$,

$M_{W}^{2}=\frac{M_{Z}^{2}}{2}\left\\{1+\sqrt{1-\frac{4\pi\alpha}{\sqrt{2}G_{\mu}M_{Z}^{2}}\left(1+\Delta^{(n)}_{\text{\tt
NMSSMCALC}}r\right)}\right\\},$ (4.42)

where the NMSSM $\Delta r$ is taken from our implementation in NMSSMCALC. The
quantity $\Delta^{(n)}_{\text{\tt NMSSMCALC}}r$, the subscript $(n)$ denotes
the $n$-loop order, depends also on $M_{W}$. Therefore, Eq. 4.42 has to be
evaluated iteratively. In the following subsections we describe in detail the
NMSSM-specific one- and two-loop contributions as well as the higher-order SM
corrections which are included in $\Delta r$.

###  4.1 One-loop Corrections

The one-loop correction to $\Delta r$ can be written as

$\displaystyle\Delta^{(1)}r$ $\displaystyle=\frac{\Sigma_{T}^{WW}(0)-\delta
M_{W}^{2}}{M_{W}^{2}}+2\delta Z_{e}-2\frac{\delta
s_{W}}{s_{W}}+\frac{1}{2}\left(\delta Z^{\mu}+\delta Z^{e}+\delta
Z^{\nu_{\mu}}+\delta Z^{\nu_{e}}\right)$ $\displaystyle+\Delta
r_{\triangle}+\Delta r_{\Box}\,,$ (4.43)

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

$\delta
Z_{e}=\frac{1}{2}\Pi^{AA}(0)-\text{sgn}\frac{s_{W}}{c_{W}}\frac{\Sigma_{T}^{AZ}(0)}{M_{Z}^{2}}\;,$
(4.44)

where

$\Pi^{AA}(0)\equiv\frac{\Sigma_{T}^{AA}(k^{2})}{k^{2}}\bigg{|}_{k^{2}=0}$
(4.45)

and $\text{sgn}=1$ for the SM and $\text{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.44) evaluated at vanishing external momentum
contains large logarithmic contributions,
$\order{\ln(\mu_{R}^{2}/m_{q}^{2})}$, which depend on ratios of the light
quark masses and the renormalization scale $\mu_{R}$, leading to numerical
instabilities. To avoid this dependence, one writes for the light SM fermions
($f=u,d,s,c,b,e,\mu,\tau$),

$\displaystyle\Pi^{AA}(0)$ $\displaystyle=$
$\displaystyle\Pi^{AA}_{f}(0)-\text{Re}\,\Pi^{AA}_{f}(M_{Z}^{2})+\text{Re}\,\Pi^{AA}_{f}(M_{Z}^{2})+\Pi^{AA}_{\text{rem}}(0)$
(4.46) $\displaystyle=$
$\displaystyle\Delta\alpha+\text{Re}\,\Pi^{AA}_{f}(M_{Z}^{2})+\Pi^{AA}_{\text{rem}}(0),$

where $\text{Re}\,\Pi^{AA}_{f}(M_{Z}^{2})$ is computed perturbatively,
$\Pi^{AA}_{\text{rem}}(0)$ contains contributions from remaining charged
particles of the model and
$\Delta\alpha=\Delta\alpha_{\text{lepton}}+\Delta\alpha_{\text{had}}^{(5)}$.
The contribution $\Delta\alpha_{\text{had}}^{(5)}$ is determined from the
dispersion relation using experimental data. We take
$\Delta\alpha_{\text{had}}^{(5)}=0.02768$ from Ref. ParticleDataGroup:2022pth
. The quantity $\Delta\alpha_{\text{lepton}}=0.03150$ includes contributions
up to three loops Steinhauser:1998rq . The counterterm of $s_{W}$ is derived
from the OS relation $s_{W}^{2}=1-M_{W}^{2}/M_{Z}^{2}$ and reads

$\delta s_{W}=\frac{c_{W}^{2}}{2s_{W}}\left[\frac{\delta
M_{Z}^{2}}{M_{Z}^{2}}-\frac{\delta M_{W}^{2}}{M_{W}^{2}}\right],$ (4.47)

where the $Z$ and $W$ boson mass counterterms $\delta M_{Z}^{2}$ and $\delta
M_{W}^{2}$ are defined in the OS scheme,

$\delta M_{Z}^{2}=\Sigma^{T}_{ZZ}(M_{Z}^{2})\;,\quad\delta
M_{W}^{2}=\Sigma^{T}_{WW}(M_{W}^{2})\;,$ (4.48)

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 $\Gamma_{W}$ and
$\Gamma_{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^{\text{run}}_{W/Z}=M_{W/Z}+\Gamma_{W/Z}^{2}/(2M_{W/Z}^{\text{run}})$ as
described in Bardin:1988xt . For $\Gamma_{Z}$ we use the experimentally
measured value while for the $W$ boson width we use the approximate formula
$\Gamma_{W}=3G_{\mu}M_{W}^{3}/(2\sqrt{2}\pi)(1+2\alpha_{s}/(3\pi))$
Freitas:2002ja neglecting possible beyond-the-SM (BSM) contributions. While
the conversion from $M_{Z}^{\text{run}}$ to $M_{Z}$ needs to be done only once
before the iteration in Eq. 4.42, the conversion from $M_{W}$ to
$M_{W}^{\text{run}}$ needs to be done at the end of the iteration. For more
details on the running/fixed-width mass conversion we refer to Freitas:2002ja
; Stal:2015zca .

The wave function (WF) counterterms of the external leptons, $\delta Z^{l}$,
and the triangle and box contributions, $\Delta r_{\triangle,\Box}$, are
evaluated for vanishing lepton masses. The WF counterterms of the external
leptons are defined in the OS scheme,

$\delta Z^{l}=-\Sigma_{ll}^{L}(0),\quad l=\mu,e,\nu_{\mu},\nu_{e}\;,$ (4.49)

where $\Sigma_{ll}^{L}(0)$ is the form-factor of the left-handed vector
component of the lepton self-energy. The Fermi-model type QED correction
appears in $\delta 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

$\frac{1}{k^{2}}\to\frac{1}{k^{2}-M_{W}^{2}},$ (4.50)

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. Sirlin:1980nh ; Sirlin:1977sv . Using this trick we
recover the well-known relation between triangle, box and vertex corrections
for the SM, derived e.g. in Denner:1991kt ,

$\displaystyle\delta_{vb}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}(\delta Z^{\mu}+\delta Z^{e}+\delta
Z^{\nu_{\mu}}+\delta
Z^{\nu_{e}})+\Delta_{\triangle}r+\Delta_{\Box}r-\frac{2}{c_{W}s_{W}}\frac{\Sigma^{T}_{AZ}(0)}{M_{Z}^{2}}\,,$
(4.51) $\displaystyle\delta_{vb}^{\text{SM}}$ $\displaystyle=$
$\displaystyle\frac{\alpha}{4\pi
s_{W}^{2}}\left(6+\frac{7-4s_{W}^{2}}{2s_{W}^{2}}\log c_{W}^{2}\right).$
(4.52)

Collecting all previous derivations leads to the known SM result

$\displaystyle\Delta^{(\alpha)}_{{\text{SM}}}r=$
$\displaystyle\Delta\alpha-\frac{c_{W}^{2}}{s^{2}_{W}}\left[\frac{\delta
M_{Z}^{2}}{M_{Z}^{2}}-\frac{\delta M_{W}^{2}}{M_{W}^{2}}\right]$ (4.53)
$\displaystyle+\text{Re}\,\Pi^{AA}_{f}(M_{Z}^{2})+\Pi^{AA}_{\text{rem}}(0)+\frac{2c_{W}}{s_{W}}\frac{\Sigma^{T}_{AZ}(0)}{M_{Z}^{2}}+\frac{\Sigma^{T}_{WW}(0)-\delta
M_{W}^{2}}{M_{W}^{2}}+\delta_{vb}^{{\text{SM}}},$

which can be written as

$\displaystyle\Delta^{(\alpha)}_{{\text{SM}}}r=\Delta\alpha-\frac{c_{W}^{2}}{s^{2}_{W}}\Delta^{(1)}\rho+\Delta
r_{\text{rem}},$ (4.54)

where $\Delta\alpha$ and $\Delta^{(1)}\rho$ contain the numerically dominant
part of the corrections.

Within the complex NMSSM, we have computed the one-loop contributions to
$\Delta^{(1)}_{{\text{NMSSM}}}r$ using the expression in Eq. (4.43) 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.

###  4.2 Combination with Known Higher-Order Corrections

In the SM, corrections $\Delta_{{\text{SM}}}^{\text{lit.}}r$ up to four-loop
order have been computed and are available in the literature. Following the
procedure developed for the MSSM Heinemeyer:2013dia and implemented in
FeynHiggs Bahl:2018qog ; Bahl:2017aev ; Bahl:2016brp ; Hahn:2013ria ;
Frank:2006yh ; Degrassi:2002fi ; Heinemeyer:1998np ; Heinemeyer:1998yj , we
have combined all available higher-order SM corrections with the presented
NMSSM contributions. For fixed $n$ ($n=1,2$) we first subtract
$\Delta^{\,(n)}_{{\text{SM}}}r$ from $\Delta^{(n)}_{{\text{NMSSM}}}r$ and re-
add $\Delta_{{\text{SM}}}^{\text{lit.}}r$ including all known-higher
corrections from the literature,

$\Delta^{(n)}_{{\text{SUSY}}}r=\Delta^{(n)}_{{\text{NMSSM}}}r-\Delta^{(n)}_{{\text{SM}}}r\,,$
(4.55)

and

$\Delta^{(n)}_{\text{\tt
NMSSMCALC}}r=\Delta^{\text{lit.}}_{{\text{SM}}}r+\Delta^{(n)}_{{\text{SUSY}}}r\,.$
(4.56)

Note that $\Delta^{(n)}_{{\text{SM}}}r$ should be computed in the same way as
for $\Delta^{(n)}_{{\text{NMSSM}}}r$.

The SM corrections from the literature consist of the following terms

$\displaystyle\Delta_{{\text{SM}}}^{\text{lit.}}r=$
$\displaystyle\Delta^{(1)}r+\Delta^{(\alpha\alpha_{s})}r+\Delta^{(\alpha\alpha_{s}^{2})}r+\Delta^{(\alpha^{2})}r$
$\displaystyle+\Delta^{(G_{\mu}^{2}m_{t}^{4}\alpha_{s})}r+\Delta^{(G_{\mu}^{2}m_{t}^{6})}r+\Delta^{(G_{\mu}^{2}m_{t}^{2}\alpha_{s}^{3})}r\,.$
(4.57)

For a complete list of references we refer the reader to Heinemeyer:2013dia .
In the following, we refer to previous works that provide analytical/numerical
results which are used in Eq. 4.57. The full two-loop QCD corrections
$\Delta^{(\alpha\alpha_{s})}r$ are taken from Halzen:1990je , partial three-
loop QCD corrections $\Delta^{(\alpha\alpha_{s}^{2})}r$ are taken from
Chetyrkin:1995js . The three-loop corrections
$\Delta^{(G_{\mu}^{2}m_{t}^{4}\alpha_{s})}r$,
$\Delta^{(G_{\mu}^{2}m_{t}^{6})}r$ to the $\rho$ parameter are taken from
Faisst:2003px and the four-loop QCD corrections
$\Delta^{(G_{\mu}^{2}m_{t}^{2}\alpha_{s}^{3})}r$ from Chetyrkin:2006bj ;
Boughezal:2006xk . For the full two-loop electroweak corrections
$\Delta^{(\alpha\alpha_{s}^{2})}r$ we use the fitting formula presented in
Awramik:2006uz . For consistency, we need to use the running-width definition
of the $W$ and $Z$ masses in the two-loop electroweak corrections, cf.
footnote 7 in Stal:2015zca , while the fixed-width masses are used in the rest
of the $\Delta r$ calculation due to the employed OS scheme.

###  4.3 NMSSM-Specific Two-Loop Corrections

As discussed in the previous sections, corrections to $\Delta 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
self-energies and can be parameterized in terms of the loop corrections to the
$\rho$ parameter. Thereby, the two-loop corrections to the $\rho$ parameter,
$\Delta^{(2)}_{{\text{SUSY}}}\rho$, discussed in Section 3, can be computed in
an efficient way. However, $\Delta^{(2)}_{{\text{NMSSM}}}\rho$ also includes
the SM-like corrections which are already taken into account in
$\Delta^{\text{lit.}}_{{\text{SM}}}r$. Hence, the SM-like contributions have
to be subtracted from the $\rho$ parameter in order to avoid double counting.
We define the two-loop SUSY correction to $\Delta\rho$ as

$\Delta^{(\alpha_{i}^{2})}_{{{\text{SUSY}}}}\rho=\Delta^{(\alpha_{i}^{2})}_{{\text{NMSSM}}}\rho-\Delta^{(\alpha_{i}^{2})}_{{\text{SM}}}\rho\,,$
(4.58)

with
$\alpha_{i}^{2}=\alpha_{t}\alpha_{s},\,\left(\alpha_{\lambda}+\alpha_{\kappa}+\alpha_{t}\right)^{2}$.
The final expression for $\Delta_{{\text{SUSY}}}^{(2)}r$ including the one and
two-loop corrections, that is used together with Eq. 4.56 in Eq. (4.42) to
compute the $W$ mass in NMSSMCALC, reads

$\Delta^{(2)}_{{\text{SUSY}}}r=\Delta^{(1)}_{{\text{SUSY}}}r-\frac{c_{W}^{2}}{s_{W}^{2}}(\Delta_{{\text{SUSY}}}^{(\alpha_{t}\alpha_{s})}\rho+\Delta_{{\text{SUSY}}}^{\left(\alpha_{\lambda}+\alpha_{\kappa}+\alpha_{t}\right)^{2}}\rho)\,,$
(4.59)

where $\Delta^{(1)}_{{\text{SUSY}}}r$ is defined via Eq. 4.55.

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 Heinemeyer:2013dia ; Stal:2015zca , once all NMSSM particle masses
are chosen to be well above the electroweak scale.

##  5 Numerical Analysis

In this section, we investigate the phenomenological impact of the ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$ corrections on the
electroweak $\rho$ parameter and on the $W$ boson mass. Detailed studies of
the one-loop and the leading two-loop corrections can be found in
Domingo:2011uf ; Stal:2015zca ; Athron:2022isz . Furthermore, we estimate the
uncertainty due to missing higher-order SUSY corrections to $\Delta\rho$ 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)_{\mu}$ 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 Tang:2022pxh ;
Domingo:2022pde for the NMSSM, and to Bagnaschi:2022qhb 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.

###  5.1 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 ${\cal
O}(\alpha_{s}\alpha_{t}+(\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})$
Dao:2021khm . 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

$\displaystyle 122\,\mathrm{GeV}\leq m_{h}\leq 128\,\mathrm{GeV}\;.$ (5.60)

In the following, $m_{h}$ will always denote the mass of the SM-like 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 ParticleDataGroup:2022pth ; 3) the masses of the
stops, lightest sbottom and the gluino satisfy the lower bounds taken from the
2022 Particle Data Group review ParticleDataGroup:2022pth , in particular

$\displaystyle M_{\chi_{1}^{\pm}}$ $\displaystyle>94\,~{}\text{GeV},$
$\displaystyle M_{\chi^{0}_{2}}$ $\displaystyle>62.4\,~{}\text{GeV},$
$\displaystyle M_{\chi^{0}_{3}}$ $\displaystyle>99.9\,~{}\text{GeV},$
$\displaystyle M_{\chi^{0}_{4}}$ $\displaystyle>116\,~{}\text{GeV},$
$\displaystyle M_{\tilde{e}_{1}}$ $\displaystyle>107\,~{}\text{GeV},$
$\displaystyle M_{\tilde{\mu}_{1}}$ $\displaystyle>94\,~{}\text{GeV},$
$\displaystyle M_{\tilde{\tau}_{1}}$ $\displaystyle>81.9\,~{}\text{GeV},$
$\displaystyle M_{\tilde{\nu}_{e/\mu/\tau}}$
$\displaystyle>94\,~{}\text{GeV},$ $\displaystyle M_{\tilde{b}_{1}}$
$\displaystyle>1270\,~{}\text{GeV},$ $\displaystyle M_{\tilde{t}_{1}}$
$\displaystyle>1310\,~{}\text{GeV},$ $\displaystyle M_{\tilde{g}}$
$\displaystyle>2300\,~{}\text{GeV}.$ (5.61)

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
Bahl:2022igd which contains HiggsBounds-5 Bechtle:2020pkv , to check if the
parameter points pass all the exclusion limits from the searches at LEP,
Tevatron and the LHC, and HiggsSignals-2 Bechtle:2020uwn to check if the
points are consistent with the LHC data for a 125 GeV Higgs boson within
2$\sigma$. To constrain the SUSY fermionic and scalar sector we use SModelS-2
Alguero:2021dig 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
Hahn:2000kx and FormCalc-9.8 Hahn:formcalc . 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
Beenakker:1996ed ; Beenakker:1999xh .

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 Shang:2023gfy in order to efficiently find phenomenologically
viable parameter regions. Table 1 summarizes the ranges applied in the
parameter scan.

| $t_{\beta}$ | $\lambda$ | $\kappa$ | $M_{1},M_{2}$ | $A_{t}$ | $m_{\tilde{Q}_{3}}$ | $m_{\tilde{t}_{R}}$ | $m_{\tilde{e}_{iR}},m_{\tilde{L}_{i}}$ | $M_{H^{\pm}}$ | $A_{\kappa}$ | $\lvert\mu_{\text{eff}}\rvert$
---|---|---|---|---|---|---|---|---|---|---|---
|  |  | in TeV
min | 1 | 0 | -1 | 0.1 | -5 | 1 | 1 | 0.1 | 0.6 | -3 | 0.1
max | 20 | 1 | 1 | 1 | 5 | 3 | 3 | 2 | 3 | 3 | 2
Table 1: Input parameters for the NMSSM scan. All parameters have been varied
independently between the given minimum and maximum values.

The SM input parameters, that are relevant for the calculation of $\Delta\rho$
and the $W$ mass111In the calculation of all Higgs masses
$\alpha(M_{Z})=1/127.955$ and $M_{W}=80.377~{}~{}\text{GeV}$ are chosen as
input parameters., are taken from ParticleDataGroup:2022pth

$\alpha_{0}$ | = | 1/137.035999084 , |  $\alpha^{\overline{\mbox{MS}}}_{s}(M_{Z})$ | = | 0.1179 ,
---|---|---|---|---|---
$G_{F}$ | = |  $1.1663788\cdot 10^{-5}$ GeV-2 , |  $M_{Z}$ | = | 91.1876 GeV ,
$m_{t}$ | = | 172.69 GeV , |  $m^{\overline{\mbox{MS}}}_{b}(m_{b}^{\overline{\mbox{MS}}})$ | = | 4.18 GeV ,
$m_{\tau}$ | = | 1.77686 GeV , |  $\Delta\alpha_{\text{had}}^{(5)}(M_{Z}^{2})$ | = | 0.02768 ,
$\Delta\alpha_{\text{lepton}}(M_{Z}^{2})$ | = | 0.03150 . |  |  |
(5.62)

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,

$\displaystyle\text{\tt P1}:\quad$ $\displaystyle
m_{\tilde{t}_{R}}=2002\,~{}\text{GeV}\,,\;m_{\tilde{Q}_{3}}=2803\,~{}\text{GeV}\,,\;m_{\tilde{b}_{R}}=2765\,~{}\text{GeV}\,,\;$
$\displaystyle
m_{\tilde{L}_{1,2}}=565\,~{}\text{GeV}\,,\;m_{\tilde{e}_{R},\tilde{\mu}_{R}}=374\,~{}\text{GeV}\,,\;$
$\displaystyle
m_{\tilde{L}_{3}}=575\,~{}\text{GeV}\,,\;m_{\tilde{\tau}_{R}}=981\,~{}\text{GeV}\,,\;$
$\displaystyle|A_{u,c,t}|=2532\,~{}\text{GeV}\,,\;|A_{d,s,b}|=1885\,~{}\text{GeV}\,,\;|A_{e,\mu,\tau}|=1170\,~{}\text{GeV}\,,\;$
$\displaystyle|M_{1}|=133\,~{}\text{GeV},\;|M_{2}|=166\,~{}\text{GeV}\,,\;|M_{3}|=2300\,~{}\text{GeV}\,,$
$\displaystyle\lambda=0.301\,~{}\text{GeV},\;\kappa=0.299\,~{}\text{GeV}\,,\;\tan\beta=4.42\,~{}\text{GeV}\,,$
$\displaystyle\mu_{\text{eff}}=254\,~{}\text{GeV},\;\text{Re}A_{\kappa}=-791\,~{}\text{GeV}\,,\;M_{H^{\pm}}=1090\,~{}\text{GeV}\,,$
$\displaystyle\varphi_{A_{e,\mu,\tau}}=0\,,\;\varphi_{A_{d,s,b}}=\pi\,,\;\varphi_{A_{u,c,t}}=\varphi_{M_{1}}=\varphi_{M_{2}}=\varphi_{M_{3}}=0\;.$

The resulting spectrum of the Higgs boson masses at ${\cal
O}(\alpha_{s}\alpha_{t}+(\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{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\,\mathrm{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-, bino- or 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.

${H_{1}}$ | ${H_{2}}$ | ${H_{3}}$ | ${H_{4}}$ | ${H_{5}}$
---|---|---|---|---
125.4 | 230.6 | 770.5 | 1088.0 | 1090.1
$h_{u}$ | $h_{s}$ | $a_{s}$ | $a$ | $h_{d}$
Table 2: The parameter point P1: Higgs masses in GeV and the main components are shown. ${\tilde{\chi}_{1}^{0}}$ | ${\tilde{\chi}_{2}^{0}}$ | ${\tilde{\chi}_{3}^{0}}$ | ${\tilde{\chi}_{4}^{0}}$ | ${\tilde{\chi}_{5}^{0}}$ | ${\tilde{\chi}_{1}^{+}}$ | ${\tilde{\chi}_{2}^{+}}$
---|---|---|---|---|---|---
113.9 | 145.07 | 261.95 | 295.74 | 509.49 | 132.93 | 294.98
Table 3: Electroweakino masses in GeV for the parameter point P1.

###  5.2 Results for the $\rho$ Parameter

In the following, we investigate the prediction for the $\rho$ parameter
starting from the parameter point P1 as a function of the NMSSM-specific
superpotential parameters $\lambda$ and $\kappa$. In order to avoid negative
squared tree-level masses, we simultaneously vary both $\lambda$ and $\kappa$
around their values $\lambda_{0}=0.301$ and $\kappa_{0}=0.299$ such that the
ratios $\lambda/\lambda_{0}$ and $\kappa/\kappa_{0}$, respectively, are kept
equal. Furthermore, we vary
$A_{\kappa}=A_{\kappa_{0}}\left(1-(\lambda-\lambda_{0})/4\right)$ to avoid a
tachyonic tree-level Higgs mass even for very large values of $\lambda$. All
other parameters of P1 are kept fixed. The upper panel of Fig. 1 (left) shows
the $\Delta\rho$ parameter computed at one-loop (black), ${{\cal
O}(\alpha_{t}\alpha_{s})}$ (blue) and ${{\cal O}(\alpha_{\text{new}}^{2})}$
(red) as a function of $\sqrt{\lambda^{2}+\kappa^{2}}$ (lower axis) and the
obtained Higgs mass including all available two-loop corrections (upper axis).
We here introduced for better readability the notation

$\displaystyle\alpha_{\text{new}}^{2}$
$\displaystyle\equiv\alpha_{t}\alpha_{s}+(\alpha_{\lambda}+\alpha_{\kappa}+\alpha_{t})^{2}\;,$
(5.64a) $\displaystyle s_{\lambda\kappa}$
$\displaystyle\equiv\sqrt{\lambda^{2}+\kappa^{2}}\;.$ (5.64b)

Since the prediction of the $\rho$ parameter strongly depends on the top quark
mass, it is shown using both the OS (full lines) and the
$\overline{\text{DR}}$ (dashed lines) renormalization schemes in the top/stop
sector. The $W$ boson mass entering the $\rho$ parameter prediction has been
chosen to be the value for $M_{W}$ obtained at ${{\cal
O}(\alpha_{\text{new}}^{2})}$, cf. Section 5.3.

Figure 1: Upper left: The $\rho$ parameter at full one-loop order (black),
two-loop ${{\cal O}(\alpha_{t}\alpha_{s})}$ (blue) and two-loop ${{\cal
O}(\alpha_{\text{new}}^{2})}$ (red) as function of
$\sqrt{\lambda^{2}+\kappa^{2}}$. The dashed lines are for the
$\overline{\text{DR}}$ scheme in the top/stop sector while the full lines are
for the OS scheme. Lower left: Renormalization scheme dependence
$\Delta^{\rho}_{\text{ren}}$ of the $\rho$ parameter at full one-loop order
(black), at two-loop ${{\cal O}(\alpha_{t}\alpha_{s})}$ (blue) and ${{\cal
O}(\alpha_{\text{new}}^{2})}$ (red). Right: The part of the Higgs- and EW-ino
mass spectrum causing larger deviations.

Note that points with $s_{\lambda\kappa}>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 $\lambda,\kappa$ would be required in order to obtain sizeable corrections.
We therefore plot our results for values up to $s_{\lambda\kappa}=3$. In the
OS scheme, the ${{\cal O}(\alpha_{t}\alpha_{s})}$ and ${{\cal
O}(\alpha_{\text{new}}^{2})}$ corrections are both negative. In the region
that is free of low-energy Landau poles, the ${{\cal
O}(\alpha_{t}\alpha_{s})}$ corrections reduce the one-loop result by 10% while
the ${{\cal O}(\alpha_{\text{new}}^{2})}$ corrections reduce the ${{\cal
O}(\alpha_{t}\alpha_{s})}$ results by 4%. In contrast, the ${{\cal
O}(\alpha_{t}\alpha_{s})}$ and ${{\cal O}(\alpha_{\text{new}}^{2})}$
corrections are both positive in the $\overline{\text{DR}}$ scheme. The
${{\cal O}(\alpha_{t}\alpha_{s})}$ corrections increase the one-loop value by
9.8% while the ${{\cal O}(\alpha_{\text{new}}^{2})}$ corrections add 0.6% on
top of the ${{\cal O}(\alpha_{t}\alpha_{s})}$ results. In the region
$s_{\lambda\kappa}>2$, the ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$ corrections in the two
renormalization schemes start to deviate differently from the ${{\cal
O}(\alpha_{t}\alpha_{s})}$ corrections. The magnitude of the ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$ corrections becomes
smaller in the OS scheme, but larger in the $\overline{\text{DR}}$ scheme.

We observe an increase in $\Delta\rho$ with increasing $s_{\lambda\kappa}$
which first is rather weak and becomes stronger for very large values of
$s_{\lambda\kappa}$. 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_{\lambda\kappa}$, while the masses of their corresponding EW-
ino states as function of $s_{\lambda\kappa}$ are shown in the lower plot.

We conclude this section with a discussion about numerical instabilities that
can appear in the ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$ corrections to
$\Delta\rho$, especially for small values of $\lambda$ and $\kappa$. The
corrections to the $\Delta\rho$ 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
$\Delta\rho$, 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 ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{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
$\mathcal{O}(\alpha_{\lambda}+\alpha_{\kappa})$ 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 $\rho$ parameter prediction in the first place, the $M_{W}$
prediction which discussed in Section 5.3, is very sensitive to the value of
$\Delta\rho$. 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 $\Delta\rho$.
Therefore, the $M_{W}$ prediction in NMSSMCALC at ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$ can only be used for
parameter points that do not suffer from large numerical instabilities. If no
convergence is found at ${{\cal
O}((\alpha_{t}+\alpha_{\lambda}+\alpha_{\kappa})^{2})}$ in $M_{W}$, the
program automatically falls back to the ${{\cal O}(\alpha_{t}\alpha_{s})}$
predictions for both the $W$-mass and $\Delta\rho$. Furthermore, some
parameter points may feature tachyonic tree-level masses in the gaugeless
limit, while the actual tree-level masses are positive. While Ref.
Domingo:2021jdg 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. Domingo:2021jdg 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 diagrams222If 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..

In Fig. 1 we observe good convergence for $s_{\lambda\kappa}>0.1$ for $M_{W}$
at ${{\cal O}(\alpha_{\text{new}}^{2})}$, cf. Section 5.3. For
$s_{\lambda\kappa}<0.1$ the $M_{W}$ prediction does not converge and therefore
the $\rho$ parameter prediction is used at ${{\cal O}(\alpha_{t}\alpha_{s})}$
which explains the jump of the red line onto the blue line.

Uncertainty Estimate for $\Delta\rho$

In order to estimate the uncertainty due to missing higher-order corrections
to the $\rho$ parameter, we define the renormalization scheme dependence of
the $\rho$ parameter at a given loop order as

$\displaystyle\Delta_{\text{ren}}^{\rho}=\frac{\Delta\rho^{\text{OS}}-\Delta\rho^{\overline{\text{DR}}}}{\Delta\rho^{{\overline{\text{DR}}}}}.$
(5.65)

The lower panel of Fig. 1 (left) shows the resulting
$\Delta_{\text{ren}}^{\rho}$ obtained at the three considered loop orders as a
function of $s_{\lambda\kappa}$. We observe a renormalization scheme
dependence of up to $50\%,22\%$ and $16\%$ for the one-loop, ${{\cal
O}(\alpha_{t}\alpha_{s})}$ and ${{\cal O}(\alpha_{\text{new}}^{2})}$ results,
respectively. Therefore, including the two-loop QCD and EW corrections can
significantly reduce the theory uncertainty of the $\rho$ parameter. For the
comparison with the SM we present in Table 4 the $\rho$ parameter computed in
the SM at the corresponding one-loop order and two-loop ${\cal
O}(\alpha_{t}\alpha_{s})$ and ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$
using $\overline{\text{DR}}$ and OS renormalisation conditions in the top/stop
sector. In the SM, the ${{\cal O}(\alpha_{t}\alpha_{s})}$ and ${\cal
O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections are negative both in the
OS and the $\overline{\text{DR}}$ scheme. The renormalisation scheme
dependence $\Delta_{\text{ren}}^{\rho}$ in the SM is $55\%,\,45\%$ and $42\%$
at one-loop order, ${{\cal O}(\alpha_{t}\alpha_{s})}$ and ${\cal
O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$, respectively, which is significantly
larger than the corresponding results obtained in the NMSSM. Therefore, the
SUSY QCD contributions to the $\rho$ parameter seem to play an important role
in the reduction of the scheme dependence of the $\rho$ 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
$\rho$ parameter. In particular for the SM prediction of $M_{W}$ a full two-
loop prediction, beyond corrections to $\Delta\rho$, was found to yield only
an uncertainty of $\order{3-7\,\mathrm{MeV}}$ Degrassi:2014sxa implying
cancellations of the scheme dependence between $\Delta\rho$ and other
quantities entering $\Delta r$. For a discussion of the $M_{W}$ uncertainty in
the NMSSM, see the next section.

| one-loop | ${{\cal O}(\alpha_{t}\alpha_{s})}$ | ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$
---|---|---|---
$\overline{\text{DR}}$ | 5.734 | 5.478 | 5.295
OS | 8.918 | 7.950 | 7.509
$\Delta_{\text{ren}}^{\rho}$ [%] | 55 | 45 | 42
Table 4: The SM $\rho$ parameter multiplied by a factor of $10^{-3}$, computed
with $M_{W}=80.36\,~{}\text{GeV}$ and $m_{h}=125\,~{}\text{GeV}$ at one-loop
and two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ and ${\cal
O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$.

###  5.3 Results for the $W$ Boson Mass

In this section we discuss the prediction of $M_{W}$ at one-loop order, two-
loop ${{\cal O}(\alpha_{t}\alpha_{s})}$ and ${{\cal
O}(\alpha_{\text{new}}^{2})}$. It is important to stress that a consistent
inclusion of the known higher-order SM corrections to $\Delta r$ beyond two-
loop order, cf. Section 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 $M_{W}$:

$\displaystyle\Delta^{\text{SM}}_{m_{h}}$
$\displaystyle=M_{W}^{\text{NMSSM}}-M_{W}^{\text{SM}}(m_{h})\quad$ (5.66a)
$\displaystyle\text{and}\quad\Delta^{\alpha_{i}}_{\alpha_{j}}$
$\displaystyle=M_{W}^{(\alpha_{i})}-M_{W}^{(\alpha_{j})}.$ (5.66b)

The quantity $\Delta^{\text{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_{W}^{\text{SM}}(m_{h})$ varies with $\lambda$
and $\kappa$ as the Higgs mass prediction in the NMSSM changes. Since the SM
higher-order corrections drop out in $\Delta^{\text{SM}}_{m_{h}}$, it can be
used as a measure for the size of the genuine SUSY corrections to $M_{W}$. The
quantity $\Delta^{\alpha_{i}}_{\alpha_{j}}$ determines the size of specific
higher-order SUSY correction $\alpha_{i}$ w.r.t. the next-lowest order
$\alpha_{j}$.

Figure 2 (upper left) shows $\Delta^{\text{SM}}_{m_{h}}$ (left axis) as a
function of $s_{\lambda\kappa}=\sqrt{\lambda^{2}+\kappa^{2}}$ starting from
the parameter point P1 at one-loop (black solid) and two-loop ${{\cal
O}(\alpha_{t}\alpha_{s})}$ (blue solid) and ${{\cal
O}(\alpha_{\text{new}}^{2})}$ (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_{W}^{\text{SM}}(m_{h})$ (green dot-dashed, right axis). The lower panel
shows $\Delta^{\alpha_{t}\alpha_{s}}_{\text{one-loop}}$ (blue) and
$\Delta^{\alpha_{\text{new}}^{2}}_{\alpha_{t}\alpha_{s}}$ (red). We observe
that the NMSSM-specific corrections range between about 10-20 MeV, mostly
dominated by the one-loop corrections. 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_{\lambda\kappa}$.
Compared to the ${\cal O}(\alpha_{t}\alpha_{s})$ result the ${{\cal
O}(\alpha_{\text{new}}^{2})}$ corrections range between -0.2 MeV and +1.5 MeV
in the shown range of $s_{\lambda\kappa}$.

To get a better understanding of the individual contributions to the $W$-mass
prediction we plot the values of $\Delta r$ obtained after the $M_{W}$
iteration has converged. In Fig. 2 (right) the blue solid line shows the total
result of $\Delta r$ obtained with NMSSMCALC including all available
corrections. The green dotted line shows the one-loop SM-contributions which
also contain the $\Delta\alpha$ 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 $\Delta 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_{\lambda\kappa}\gtrsim 0.6-0.9$. Note that in Fig. 2
(right) we chose a log-scale for $|\Delta r|>10^{-5}$ and a linear scale
otherwise which sets the focus on the $s_{\lambda\kappa}$-dependence of the
two-loop SUSY corrections.

Figure 2: Upper left: the difference $\Delta^{\text{SM}}_{m_{h}}$ (left axis)
defined in Eq. 5.66a at one-loop (black), two-loop ${{\cal
O}(\alpha_{t}\alpha_{s})}$ (blue) and two-loop ${{\cal
O}(\alpha_{\text{new}}^{2})}$ (red) as a function of
$\sqrt{\lambda^{2}+\kappa^{2}}$, and $M_{W}$ in the SM (green dash-dotted,
right axis) using the loop-corrected NMSSM Higgs boson mass prediction (upper
axis). Lower left: the difference $\Delta^{\alpha_{t}\alpha_{s}}_{\text{one-
loop}}$ (blue) and $\Delta^{\alpha_{\text{new}}^{2}}_{\alpha_{t}\alpha_{s}}$
(red) defined via Eq. 5.66b. Right: Individual contributions to $\Delta r$
predicted by NMSSMCALC. For $|\Delta r|>10^{-5}$ a log-scale is chosen and a
linear scale otherwise.

For the other parameter points which pass the constraints of our scan
described in Section 5.1 we observe similar features. The ${{\cal
O}(\alpha_{t}\alpha_{s})}$ and ${{\cal O}(\alpha_{\text{new}}^{2})}$
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}$, $\Delta\alpha_{\text{had}}^{(5)}(M_{Z}^{2})$, $\alpha_{s}$ within
their $1\sigma$ PDG values ParticleDataGroup:2022pth . 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 Section 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 Bagnaschi:2022qhb ; Athron:2022isz .

|  $\Delta M_{W}$ in MeV
---|---
| P1 | BP3 | SM
$M_{Z}\pm 1\sigma$ | 1.5 | 1.5 | 1.5
$m_{t}\pm 1\,\mathrm{GeV}$ | 5.6 | 5.8 | 5.6
$\Delta\alpha_{\text{had}}^{(5)}(M_{Z}^{2})\pm 1\sigma$ | 1.27 | 1.28 | 1.27
$\alpha_{s}\pm 1\sigma$ | 0.5 | 0.6 | 0.5
$m_{h}\pm 1\,\mathrm{GeV}$ | 0.47 | 0.47 | 0.47
Table 5: The variation in the $W$ mass prediction around the central value
when varying $M_{Z}$, $\Delta\alpha_{\text{had}}^{(5)}(M_{Z}^{2})$,
$\alpha_{s}$ within $1\sigma$ of their central values and $m_{t}$ within
$172.69\pm 1\,\mathrm{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 $\pm 1\,\mathrm{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 Awramik:2003rn ; Degrassi:2014sxa ;
Bagnaschi:2022qhb and yield about 4 MeV and 3 MeV in the OS and
${\overline{\text{MS}}}$ calculation, respectively. A comparison between OS
and ${\overline{\text{MS}}}$ result, however, suggest an uncertainty of about
6 MeV Degrassi:2014sxa . Regarding the missing higher-order SUSY corrections,
one can expect them to not be significantly larger than the computed ${{\cal
O}(\alpha_{t}\alpha_{s})}$ and ${{\cal O}(\alpha_{\text{new}}^{2})}$
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
possibility 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.

###  5.4 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 Athron:2014yba ; Athron:2017fvs and
SARAH Staub:2009bi ; Staub:2010jh ; Staub:2012pb ; Staub:2013tta 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 Benakli:2022gjn while the one of FlexibleSUSY was applied
to the MSSM, MRSSM and a singlet extended SM Athron:2022isz . 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 Domingo:2011uf to compute the $W$-mass in the general NMSSM as well as the
$\mathbb{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 Section 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_{\mu}$ are known to have a connection to large corrections to
$M_{W}$ Bagnaschi:2022qhb . Since $a_{\mu}$ is of increasing recent
theoretical and experimental interest Muong-2:2021ojo ; Athron:2021evk ;
Dao:2022rui ; Tang:2022pxh , we also include a comparison of the $a_{\mu}$
prediction between the various codes.

We start by discussing the incorporation of the SM higher-order corrections in
the different codes. NMSSMCALC implements the results of
$\Delta_{\text{SM}}^{\text{lit.}}r$ to a large degree analytically (cf.
Section 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
${\overline{\text{MS}}}$/$\overline{\text{DR}}$ scheme and therefore rely on
the SM ${\overline{\text{MS}}}$ fit formula for $M_{W}^{\text{SM}}$ provided
in Degrassi:2014sxa :

$M_{W}^{{\begin{subarray}{c}{\tt FlexibleSUSY}\\\ {\tt
SARAH/SPheno}\end{subarray}}}=\sqrt{M_{W}^{\text{SM
fit.}^{2}}(m_{h},m_{t},\alpha,\alpha_{s})\left[1+\frac{s_{W}^{2}}{c_{W}^{2}-s_{W}^{2}}\Delta
r_{\text{SUSY}}^{(n)}\right]},\,$ (5.67)

where $M_{W}^{\text{SM fit.}}$ 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
$\Delta_{\text{SM}}^{\text{lit.}}r$ on the value of $M_{W}$, which is
correctly taken into account in NMSSMCALC via Eq. 4.42, is lost when using fit
formulas for $M_{W}$. In NMSSMTools the $M_{W}$ dependence is partially
restored by determining $M_{W}^{\text{SM fit.}}$ from the fit formula given in
Awramik:2003rn , inverting Eq. 4.42 for $\Delta_{\text{SM}}r(M_{W}^{\text{SM
fit.}})$ and adding the $M_{W}$ dependence using a further fit formula,

$\displaystyle\Delta^{\text{{\tt NMSSMTools}}}r(M_{W})=\,$
$\displaystyle\Delta_{\text{SUSY}}\,r(M_{W})+\Delta_{\text{SM}}\,r(M_{W}^{\text{SM
fit.}})+\sum_{n=1}^{3}a_{n}\left(M_{W}-M_{W}^{\text{SM fit.}}\right)^{n}$
(5.68a) $\displaystyle\approx$
$\displaystyle\,\Delta_{\text{SUSY}}\,r(M_{W})+\Delta_{\text{SM}}^{\text{lit.}}r(M_{W})\,.$
(5.68b)

We now discuss the different treatments of SUSY input parameters. All (SM and
BSM) quantities entering Eq. 5.67 are defined in the
${\overline{\text{MS}}}$/$\overline{\text{DR}}$ scheme. In FlexibleSUSY, Eq.
5.67 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_{\text{SUSY}}^{2}=m_{\tilde{t}_{R}}m_{\tilde{Q}_{3}}$. 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_{\mu}$ at $M_{\text{SUSY}}$ rather than $M_{Z}$. Likewise, ${\tt
SARAH/SPheno}$ computes $a_{\mu}$ by default at $M_{Z}$ which we changed to
$M_{\text{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 $\overline{\text{DR}}$
parameters defined at $M_{\text{SUSY}}$. In particular, the parameter
$\text{Re}A_{\lambda}$ is used as input parameter in the
$\overline{\text{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 $\overline{\text{DR}}$ scheme. The $M_{W}$
prediction in NMSSMCALC, however, is still carried out in the OS scheme for
the SM-sector described in Section 4. NMSSMTools interprets $\tan\beta$ per
default to be defined at $M_{Z}$ rather than $M_{\text{SUSY}}$. Thus, we run
the $\tan\beta(M_{\text{SUSY}})$ used in the other codes to $M_{Z}$ (using
two-loop RGEs generated with SARAH) and use $\tan\beta(M_{Z})$ in NMSSMTools.
Furthermore, many of the involved codes also compute loop-corrected masses to
sfermions and electroweakinos which in turn are used in the $M_{W}$
prediction. We also disabled such calculations in 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
Staub:2015aea ; Slavich:2020zjv . 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_{h}^{\text{(loop)}}\approx 125\,\mathrm{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.2\,\mathrm{GeV}$ 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.42, Eq. 5.67 and Eq. 5.68, yield the precise SM result in the limit of
heavy superpartners (see the discussion in Athron:2022isz for more details on
the correct decoupling behaviour). We explicitly verified that all codes still
feature a proper decoupling behaviour after we applied the changes discussed
above.

In order to compare the $m_{h}$, $M_{W}$ and $a_{\mu}$ 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 Section 5.1 and features rather light
electroweakinos. The second, BP3, was introduced in Domingo:2022pde 333Note
that BP3 has no preferred features compared to the other parameter points
(BP1, BP2, BP4) in Domingo:2022pde but was chosen by matter of taste.. This
parameter point is characterized by large one-loop corrections to the $W$
boson mass due to very light sleptons with masses of
$\order{100\,\mathrm{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 Domingo:2022pde
using the notation introduced in Section 2,

$\displaystyle\text{\tt BP3}:\quad$ $\displaystyle
m_{\tilde{t}_{R}}=2144\,~{}\text{GeV}\,,\;m_{\tilde{Q}_{3}}=1112\,~{}\text{GeV}\,,\;m_{\tilde{b}_{R}}=1539\,~{}\text{GeV}\,,\;$
$\displaystyle
m_{\tilde{L}_{1,2}}=131.9\,~{}\text{GeV}\,,\;m_{\tilde{e}_{R},\tilde{\mu}_{R}}=103.6\,~{}\text{GeV}\,,\;$
$\displaystyle
m_{\tilde{L}_{3}}=205.2\,~{}\text{GeV}\,,\;m_{\tilde{\tau}_{R}}=238.6\,~{}\text{GeV}\,,\;|A_{u,c,t}|=3971.2\,~{}\text{GeV}\,,\;$
$\displaystyle|A_{d,s,b}|=1210.3\,~{}\text{GeV}\,,\;|A_{e,\mu}|=3643\,~{}\text{GeV}\,,\;|A_{\tau}|=2052.4\,~{}\text{GeV}\,,\;$
$\displaystyle|M_{1}|=178.3\,~{}\text{GeV},\;|M_{2}|=128.6\,~{}\text{GeV}\,,\;|M_{3}|=1757.6\,~{}\text{GeV}\,,$
$\displaystyle\lambda=0.1229\,~{}\text{GeV},\;\kappa=0.0128\,~{}\text{GeV}\,,\;\tan\beta=8.7199\,~{}\text{GeV}\,,$
$\displaystyle\mu_{\text{eff}}=212\,~{}\text{GeV},\;\text{Re}A_{\kappa}=-10.48.\,~{}\text{GeV}\,,\;\text{Re}A_{\lambda}=2245\,~{}\text{GeV}\,,$
$\displaystyle\varphi_{A_{e,\mu,\tau}}=0\,,\;\varphi_{A_{u,c,t}}=\pi\,,\;\varphi_{A_{d,s,b}}=\varphi_{M_{1}}=\varphi_{M_{2}}=\varphi_{M_{3}}=0\;.$

In Table 6 we compare the prediction for $m_{h}$, $M_{W}$ and $a_{\mu}$
obtained with all four codes for P1 and BP3. Additionally, we list the
uncertainty estimates for $M_{W}$ and $a_{\mu}$ 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. Section 5.3), the uncertainty due to the use of the fit
formula in Eq. 5.68 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
$\gtrsim 1\,\mathrm{GeV}$ Dao:2019qaz ; Slavich:2020zjv ; Dao:2021khm . Since
the parameter point P1 is now interpreted with stop masses and stop-Higgs
trilinear couplings defined in the $\overline{\text{DR}}$ scheme, rather than
OS, the obtained Higgs boson mass for P1 is no longer at 125 GeV but around
$m_{h}\approx 119\,\mathrm{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_{\mu}$ 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 uncertainty of $4-7\,\mathrm{MeV}$ which
is about the size of the $M_{W}$ uncertainty of the SM prediction, cf. Section
5.3. The result for $a_{\mu}$ 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_{\mu}$ is calculated between
$M_{\text{SUSY}}/2$ and $2M_{\text{SUSY}}$. This is due to the scale-choice
used in the calculation of $a_{\mu}$ which is fixed in NMSSMCALC to
$\mu^{2}=m_{\tilde{L}_{2}}m_{\tilde{R}_{2}}$ 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_{\mu}\approx 3\times 10^{-10}$ when using the lightest neutralino mass as
the renormalisation scale.

|  | FlexibleSUSY | NMSSMCALC | NMSSMTools | SARAH/SPheno
---|---|---|---|---|---
P1 |  $m_{h}$ [GeV] | 119.77 | 119.19 | 118.61 | 118.95
|  $M_{W}$ [MeV] | 80366.3 | 80365.7 |  80370.8$\pm 23$ | 80366.2
| $a_{\mu}\times 10^{9}$ |  0.29$\pm 0.01$ | 0.256 |  0.329$\pm 0.03$ | 0.33
BP3 |  $m_{h}$ [GeV] | 125.60 | 125.63 | 124.63 | 123.97
|  $M_{W}$ [MeV] | 80396.9 | 80400.0 |  80404.2$\pm 22$ | 80401.3
| $a_{\mu}\times 10^{9}$ |  2.98$\pm 0.45$ | 2.89 |  3.19$\pm 0.34$ | 3.70
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 $\overline{\text{DR}}$ scheme (which is why the Higgs mass is no longer at
125 GeV). Figure 3: Comparison between NMSSMCALC (red solid), NMSSMTools (blue
dashed), SARAH/SPheno (magenta dash-dotted) and FlexibleSUSY (black dotted)
for the parameter point P1 in the $\overline{\text{DR}}$ scheme as a function
of $\sqrt{\lambda^{2}+\kappa^{2}}$. Left: $W$-mass prediction. Right:
prediction of the SM-like Higgs boson mass $m_{h_{1}}$.

In the following we compare the dependence of the $M_{W}$ prediction onto the
superpotential parameters $\lambda$ and $\kappa$ 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 ${\tt FlexibleSUSY}$ and ${\tt
SARAH/SPheno}$ agree almost perfectly and only start to slightly deviate for
very large values of $s_{\lambda\kappa}\gtrsim 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 SM-like Higgs boson mass $m_{h_{2}}$, the mass of the
lightest CP-even 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_{\lambda\kappa}$ which is rather flat in the
shown range of $s_{\lambda\kappa}\leq 0.3$. For larger values of
$s_{\lambda\kappa}$ 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_{\lambda\kappa}$. 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 SM-
uncertainty. The prediction obtained with NMSSMTools, however, seems to behave
much flatter for large values of $s_{\lambda\kappa}$. 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 Figs. 3 and 4, if we
remove the $M_{W}$-restoring fit function in Eq. 5.68 from its prediction.
Thus, we suspect that the fit coefficients in Eq. 5.68 in NMSSMTools are
outdated. Finally, in the MSSM limit $s_{\lambda\kappa}\to 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 Heinemeyer:2013dia . 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.

Figure 4: Comparison between NMSSMCALC (red solid), NMSSMTools (blue dashed),
SARAH/SPheno (magenta dash-dotted) and FlexibleSUSY (black dotted) for the
parameter point BP3 as a function of $\sqrt{\lambda^{2}+\kappa^{2}}$. Left:
$W$-mass prediction. Right (up to down): prediction of the SM-like Higgs boson
mass $m_{h_{2}}$, the lightest CP-even Higgs mass $M_{h_{1}}$ and lightest CP-
odd mass $m_{a_{1}}$, respectively.

###  5.5 CP-Violating Effects in the $M_{W}$ Prediction

In this section we study the effect of CP-violating beyond-the-SM phases on
the $M_{W}$ calculation. We consider the parameter point P1, which was
initially defined in the CP-conserving 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 ${\cal O}(\alpha_{\text{new}}^{2})$ 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
$\Delta^{\text{SM}}_{m_{h}}$ defined in Eq. 5.66b.

Figure 5: Left panel: $M_{W}$ at two-loop ${\cal O}(\alpha_{\text{new}}^{2})$
as a function of the complex phases $\varphi_{M_{1}}$ (black),
$\varphi_{M_{2}}$ (blue), $\varphi_{M_{3}}$ (red) and $\varphi_{A_{t}}$
(orange) for the parameter point P1. Right panel: same as the left panel but
showing the difference between the NMSSM and SM $M_{W}$ prediction that has
been obtained using NMSSMCALC and the value of the SM-like Higgs mass
prediction for the particular values of $\varphi_{i}$, cf. Eq. 5.66b.

We find that the phase of $\varphi_{M_{2}}$ has the largest impact on $M_{W}$
for this parameter point of about 2 MeV while the overall SUSY corrections
$\Delta^{\text{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 $\varphi_{A_{t}}$ and $\varphi_{M_{3}}$ only impact $M_{W}$ at the sub-
MeV level, which is in agreement with the findings in Heinemeyer:2006px 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}$.

##  6 Conclusions

In this paper, a consistent inclusion of the two-loop ${{\cal
O}(\alpha_{t}\alpha_{s})}$ and ${{\cal O}(\alpha_{\text{new}}^{2})}$
corrections to the $\rho$ parameter and its application 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/$\overline{\text{DR}}$
scheme and conveniently allows to switch between OS and $\overline{\text{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 $\rho$ parameter prediction due to missing higher-orders. We showed that
the ${{\cal O}(\alpha_{t}\alpha_{s})}$ and ${{\cal
O}(\alpha_{\text{new}}^{2})}$ corrections to the $\rho$ 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 ${{\cal O}(\alpha_{t}\alpha_{s})}$ and ${{\cal
O}(\alpha_{\text{new}}^{2})}$ 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_{\mu}$, 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.

## Acknowledgements

The authors thank G. Weiglein for various discussions about the SM higher-
order corrections. The research of MM was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762
- TRR 257. T.N.D is funded by the Vietnam National Foundation for Science and
Technology Development (NAFOSTED) under grant number 103.01-2020.17. M.G.
acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy — EXC 2121 “Quantum
Universe” — 390833306. This work has been partially funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) — 491245950.

## References

  * (1) ATLAS Collaboration collaboration, _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys.Lett._ B716 (2012) 1 [1207.7214].
  * (2) CMS Collaboration collaboration, _Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC_ , _Phys.Lett._ B716 (2012) 30 [1207.7235].
  * (3) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (4) E. Bagnaschi, M. Chakraborti, S. Heinemeyer, I. Saha and G. Weiglein, _Interdependence of the new “MUON G-2” result and the W-boson mass_ , _Eur. Phys. J. C_ 82 (2022) 474 [2203.15710].
  * (5) M. Awramik, M. Czakon, A. Freitas and G. Weiglein, _Precise prediction for the W boson mass in the standard model_ , _Phys. Rev. D_ 69 (2004) 053006 [hep-ph/0311148].
  * (6) G. Degrassi, P. Gambino and P.P. Giardino, _The $m_{\scriptscriptstyle W}-m_{\scriptscriptstyle Z}$ interdependence in the Standard Model: a new scrutiny_, _JHEP_ 05 (2015) 154 [1411.7040].
  * (7) P. Athron, M. Bach, D.H.J. Jacob, W. Kotlarski, D. Stöckinger and A. Voigt, _Precise calculation of the W boson pole mass beyond the standard model with FlexibleSUSY_ , _Phys. Rev. D_ 106 (2022) 095023 [2204.05285].
  * (8) CDF collaboration, _High-precision measurement of the W boson mass with the CDF II detector_ , _Science_ 376 (2022) 170.
  * (9) J. de Blas, M. Pierini, L. Reina and L. Silvestrini, _Impact of the Recent Measurements of the Top-Quark and W-Boson Masses on Electroweak Precision Fits_ , _Phys. Rev. Lett._ 129 (2022) 271801 [2204.04204].
  * (10) A. Sirlin, _Radiative Corrections in the SU(2)-L x U(1) Theory: A Simple Renormalization Framework_ , _Phys. Rev. D_ 22 (1980) 971.
  * (11) J.M. Yang and Y. Zhang, _Low energy SUSY confronted with new measurements of W-boson mass and muon g-2_ , _Sci. Bull._ 67 (2022) 1430 [2204.04202].
  * (12) P. Fayet, _Supergauge Invariant Extension of the Higgs Mechanism and a Model for the electron and Its Neutrino_ , _Nucl.Phys._ B90 (1975) 104.
  * (13) R. Barbieri, S. Ferrara and C.A. Savoy, _Gauge Models with Spontaneously Broken Local Supersymmetry_ , _Phys.Lett._ B119 (1982) 343.
  * (14) M. Dine, W. Fischler and M. Srednicki, _A Simple Solution to the Strong CP Problem with a Harmless Axion_ , _Phys.Lett._ B104 (1981) 199.
  * (15) H.P. Nilles, M. Srednicki and D. Wyler, _Weak Interaction Breakdown Induced by Supergravity_ , _Phys.Lett._ B120 (1983) 346.
  * (16) J. Frere, D. Jones and S. Raby, _Fermion Masses and Induction of the Weak Scale by Supergravity_ , _Nucl.Phys._ B222 (1983) 11.
  * (17) J. Derendinger and C.A. Savoy, _Quantum Effects and SU(2) x U(1) Breaking in Supergravity Gauge Theories_ , _Nucl.Phys._ B237 (1984) 307.
  * (18) J.R. Ellis, J. Gunion, H.E. Haber, L. Roszkowski and F. Zwirner, _Higgs Bosons in a Nonminimal Supersymmetric Model_ , _Phys.Rev._ D39 (1989) 844.
  * (19) M. Drees, _Supersymmetric Models with Extended Higgs Sector_ , _Int.J.Mod.Phys._ A4 (1989) 3635.
  * (20) U. Ellwanger, M. Rausch de Traubenberg and C.A. Savoy, _Particle spectrum in supersymmetric models with a gauge singlet_ , _Phys.Lett._ B315 (1993) 331 [hep-ph/9307322].
  * (21) U. Ellwanger, M. Rausch de Traubenberg and C.A. Savoy, _Higgs phenomenology of the supersymmetric model with a gauge singlet_ , _Z.Phys._ C67 (1995) 665 [hep-ph/9502206].
  * (22) U. Ellwanger, M. Rausch de Traubenberg and C.A. Savoy, _Phenomenology of supersymmetric models with a singlet_ , _Nucl.Phys._ B492 (1997) 21 [hep-ph/9611251].
  * (23) T. Elliott, S. King and P. White, _Unification constraints in the next-to-minimal supersymmetric standard model_ , _Phys.Lett._ B351 (1995) 213 [hep-ph/9406303].
  * (24) S. King and P. White, _Resolving the constrained minimal and next-to-minimal supersymmetric standard models_ , _Phys.Rev._ D52 (1995) 4183 [hep-ph/9505326].
  * (25) F. Franke and H. Fraas, _Neutralinos and Higgs bosons in the next-to-minimal supersymmetric standard model_ , _Int.J.Mod.Phys._ A12 (1997) 479 [hep-ph/9512366].
  * (26) M. Maniatis, _The Next-to-Minimal Supersymmetric extension of the Standard Model reviewed_ , _Int.J.Mod.Phys._ A25 (2010) 3505 [0906.0777].
  * (27) U. Ellwanger, C. Hugonie and A.M. Teixeira, _The Next-to-Minimal Supersymmetric Standard Model_ , _Phys.Rept._ 496 (2010) 1 [0910.1785].
  * (28) F. Domingo and T. Lenz, _W mass and Leptonic Z-decays in the NMSSM_ , _JHEP_ 07 (2011) 101 [1101.4758].
  * (29) O. Stål, G. Weiglein and L. Zeune, _Improved prediction for the mass of the W boson in the NMSSM_ , _JHEP_ 09 (2015) 158 [1506.07465].
  * (30) A. Djouadi, P. Gambino, S. Heinemeyer, W. Hollik, C. Junger and G. Weiglein, _Supersymmetric contributions to electroweak precision observables: QCD corrections_ , _Phys. Rev. Lett._ 78 (1997) 3626 [hep-ph/9612363].
  * (31) A. Djouadi, P. Gambino, S. Heinemeyer, W. Hollik, C. Junger and G. Weiglein, _Leading QCD corrections to scalar quark contributions to electroweak precision observables_ , _Phys. Rev. D_ 57 (1998) 4179 [hep-ph/9710438].
  * (32) S. Heinemeyer and G. Weiglein, _Leading electroweak two loop corrections to precision observables in the MSSM_ , _JHEP_ 10 (2002) 072 [hep-ph/0209305].
  * (33) J. Haestier, S. Heinemeyer, D. Stockinger and G. Weiglein, _Electroweak precision observables: Two-loop Yukawa corrections of supersymmetric particles_ , _JHEP_ 12 (2005) 027 [hep-ph/0508139].
  * (34) K. Ender, T. Graf, M. Muhlleitner and H. Rzehak, _Analysis of the NMSSM Higgs Boson Masses at One-Loop Level_ , _Phys.Rev._ D85 (2012) 075024 [1111.4952].
  * (35) T. Graf, R. Grober, M. Muhlleitner, H. Rzehak and K. Walz, _Higgs Boson Masses in the Complex NMSSM at One-Loop Level_ , _JHEP_ 10 (2012) 122 [1206.6806].
  * (36) M. Muhlleitner, D.T. Nhung, H. Rzehak and K. Walz, _Two-loop contributions of the order $\mathcal{O}\left({\alpha}_{t}{\alpha}_{s}\right)$ to the masses of the Higgs bosons in the CP-violating NMSSM_, _JHEP_ 1505 (2015) 128 [1412.0918].
  * (37) T. Dao, R. Gröber, M. Krause, M. Mühlleitner and H. Rzehak, _Two-loop $\mathcal{O}$ ( ${\alpha}_{t}^{2}$ ) corrections to the neutral Higgs boson masses in the CP-violating NMSSM_, _JHEP_ 08 (2019) 114 [1903.11358].
  * (38) T.N. Dao, M. Gabelmann, M. Mühlleitner and H. Rzehak, _Two-loop $\mathcal{O}$(($\alpha$t \+ $\alpha$λ \+ $\alpha$κ)2) corrections to the Higgs boson masses in the CP-violating NMSSM_, _JHEP_ 09 (2021) 193 [2106.06990].
  * (39) J. Baglio, R. Grober, M. Muhlleitner, D. Nhung, H. Rzehak et al., _NMSSMCALC: A Program Package for the Calculation of Loop-Corrected Higgs Boson Masses and Decay Widths in the (Complex) NMSSM_ , _Comput.Phys.Commun._ 185 (2014) 3372 [1312.4788].
  * (40) D.T. Nhung, M. Muhlleitner, J. Streicher and K. Walz, _Higher Order Corrections to the Trilinear Higgs Self-Couplings in the Real NMSSM_ , _JHEP_ 1311 (2013) 181 [1306.3926].
  * (41) M. Mühlleitner, D.T. Nhung and H. Ziesche, _The order $\mathcal{O}\left({\alpha}_{t}{\alpha}_{s}\right)$ corrections to the trilinear Higgs self-couplings in the complex NMSSM_, _JHEP_ 12 (2015) 034 [1506.03321].
  * (42) C. Borschensky, T.N. Dao, M. Gabelmann, M. Mühlleitner and H. Rzehak, _The trilinear Higgs self-couplings at $\mathcal{O}(\alpha_{t}^{2})$ in the CP-violating NMSSM_, _Eur. Phys. J. C_ 83 (2023) 118 [2210.02104].
  * (43) S.F. King, M. Muhlleitner, R. Nevzorov and K. Walz, _Exploring the CP-violating NMSSM: EDM Constraints and Phenomenology_ , _Nucl. Phys. B_ 901 (2015) 526 [1508.03255].
  * (44) T.N. Dao, M. Mühlleitner and A.V. Phan, _Loop-corrected Higgs masses in the NMSSM with inverse seesaw mechanism_ , _Eur. Phys. J. C_ 82 (2022) 667 [2108.10088].
  * (45) J. Baglio, T.N. Dao and M. Mühlleitner, _One-Loop Corrections to the Two-Body Decays of the Neutral Higgs Bosons in the Complex NMSSM_ , _Eur. Phys. J. C_ 80 (2020) 960 [1907.12060].
  * (46) T.N. Dao, L. Fritz, M. Krause, M. Mühlleitner and S. Patel, _Gauge dependences of higher-order corrections to NMSSM Higgs boson masses and the charged Higgs Decay ${H^{\pm}\rightarrow W^{\pm}h_{i}}$_, _Eur. Phys. J. C_ 80 (2020) 292 [1911.07197].
  * (47) T.N. Dao, M. Muhlleitner, S. Patel and K. Sakurai, _One-loop Corrections to the Two-Body Decays of the Charged Higgs Bosons in the Real and Complex NMSSM_ , _Eur. Phys. J. C_ 81 (2021) 340 [2012.14889].
  * (48) D.A. Ross and M.J.G. Veltman, _Neutral Currents in Neutrino Experiments_ , _Nucl. Phys. B_ 95 (1975) 135.
  * (49) M.J.G. Veltman, _Limit on Mass Differences in the Weinberg Model_ , _Nucl. Phys. B_ 123 (1977) 89.
  * (50) J.J. van der Bij and F. Hoogeveen, _Two Loop Correction to Weak Interaction Parameters Due to a Heavy Fermion Doublet_ , _Nucl. Phys. B_ 283 (1987) 477.
  * (51) M. Consoli, W. Hollik and F. Jegerlehner, _The Effect of the Top Quark on the M(W)-M(Z) Interdependence and Possible Decoupling of Heavy Fermions from Low-Energy Physics_ , _Phys. Lett. B_ 227 (1989) 167.
  * (52) J. Fleischer, O.V. Tarasov and F. Jegerlehner, _Two loop heavy top corrections to the rho parameter: A Simple formula valid for arbitrary Higgs mass_ , _Phys. Lett. B_ 319 (1993) 249.
  * (53) S. Hessenberger and W. Hollik, _Two-loop corrections to the $\rho$ parameter in Two-Higgs-Doublet Models_, _Eur. Phys. J. C_ 77 (2017) 178 [1607.04610].
  * (54) A. Djouadi and C. Verzegnassi, _Virtual Very Heavy Top Effects in LEP / SLC Precision Measurements_ , _Phys. Lett. B_ 195 (1987) 265.
  * (55) A. Djouadi, _O(alpha alpha-s) Vacuum Polarization Functions of the Standard Model Gauge Bosons_ , _Nuovo Cim. A_ 100 (1988) 357.
  * (56) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * (57) M. Steinhauser, _Leptonic contribution to the effective electromagnetic coupling constant up to three loops_ , _Phys. Lett. B_ 429 (1998) 158 [hep-ph/9803313].
  * (58) D.Y. Bardin, A. Leike, T. Riemann and M. Sachwitz, _Energy Dependent Width Effects in e+ e- Annihilation Near the Z Boson Pole_ , _Phys. Lett. B_ 206 (1988) 539.
  * (59) A. Freitas, W. Hollik, W. Walter and G. Weiglein, _Electroweak two loop corrections to the $M_{W}-M_{Z}$ mass correlation in the standard model_, _Nucl. Phys. B_ 632 (2002) 189 [hep-ph/0202131].
  * (60) A. Sirlin, _Current Algebra Formulation of Radiative Corrections in Gauge Theories and the Universality of the Weak Interactions_ , _Rev. Mod. Phys._ 50 (1978) 573.
  * (61) A. Denner, _Techniques for calculation of electroweak radiative corrections at the one loop level and results for W physics at LEP-200_ , _Fortsch. Phys._ 41 (1993) 307 [0709.1075].
  * (62) S. Heinemeyer, W. Hollik, G. Weiglein and L. Zeune, _Implications of LHC search results on the W boson mass prediction in the MSSM_ , _JHEP_ 12 (2013) 084 [1311.1663].
  * (63) H. Bahl, T. Hahn, S. Heinemeyer, W. Hollik, S. Paßehr, H. Rzehak et al., _Precision calculations in the MSSM Higgs-boson sector with FeynHiggs 2.14_ , _Comput. Phys. Commun._ 249 (2020) 107099 [1811.09073].
  * (64) H. Bahl, S. Heinemeyer, W. Hollik and G. Weiglein, _Reconciling EFT and hybrid calculations of the light MSSM Higgs-boson mass_ , _Eur. Phys. J. C_ 78 (2018) 57 [1706.00346].
  * (65) H. Bahl and W. Hollik, _Precise prediction for the light MSSM Higgs boson mass combining effective field theory and fixed-order calculations_ , _Eur. Phys. J. C_ 76 (2016) 499 [1608.01880].
  * (66) T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak and G. Weiglein, _High-Precision Predictions for the Light CP -Even Higgs Boson Mass of the Minimal Supersymmetric Standard Model_ , _Phys. Rev. Lett._ 112 (2014) 141801 [1312.4937].
  * (67) M. Frank, T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak and G. Weiglein, _The Higgs Boson Masses and Mixings of the Complex MSSM in the Feynman-Diagrammatic Approach_ , _JHEP_ 02 (2007) 047 [hep-ph/0611326].
  * (68) G. Degrassi, S. Heinemeyer, W. Hollik, P. Slavich and G. Weiglein, _Towards high precision predictions for the MSSM Higgs sector_ , _Eur. Phys. J. C_ 28 (2003) 133 [hep-ph/0212020].
  * (69) S. Heinemeyer, W. Hollik and G. Weiglein, _The Masses of the neutral CP - even Higgs bosons in the MSSM: Accurate analysis at the two loop level_ , _Eur. Phys. J. C_ 9 (1999) 343 [hep-ph/9812472].
  * (70) S. Heinemeyer, W. Hollik and G. Weiglein, _FeynHiggs: A Program for the calculation of the masses of the neutral CP even Higgs bosons in the MSSM_ , _Comput. Phys. Commun._ 124 (2000) 76 [hep-ph/9812320].
  * (71) F. Halzen and B.A. Kniehl, _$\Delta$ r beyond one loop_, _Nucl. Phys. B_ 353 (1991) 567.
  * (72) K.G. Chetyrkin, J.H. Kuhn and M. Steinhauser, _QCD corrections from top quark to relations between electroweak parameters to order alpha-s**2_ , _Phys. Rev. Lett._ 75 (1995) 3394 [hep-ph/9504413].
  * (73) M. Faisst, J.H. Kuhn, T. Seidensticker and O. Veretin, _Three loop top quark contributions to the rho parameter_ , _Nucl. Phys. B_ 665 (2003) 649 [hep-ph/0302275].
  * (74) K.G. Chetyrkin, M. Faisst, J.H. Kuhn, P. Maierhofer and C. Sturm, _Four-Loop QCD Corrections to the Rho Parameter_ , _Phys. Rev. Lett._ 97 (2006) 102003 [hep-ph/0605201].
  * (75) R. Boughezal and M. Czakon, _Single scale tadpoles and O(G(F m(t)**2 alpha(s)**3)) corrections to the rho parameter_ , _Nucl. Phys. B_ 755 (2006) 221 [hep-ph/0606232].
  * (76) M. Awramik, M. Czakon and A. Freitas, _Electroweak two-loop corrections to the effective weak mixing angle_ , _JHEP_ 11 (2006) 048 [hep-ph/0608099].
  * (77) T.-P. Tang, M. Abdughani, L. Feng, Y.-L.S. Tsai, J. Wu and Y.-Z. Fan, _NMSSM neutralino dark matter for CDF II W-boson mass and muon g $-$ 2 and the promising prospect of direct detection_, _Sci. China Phys. Mech. Astron._ 66 (2023) 239512 [2204.04356].
  * (78) F. Domingo, U. Ellwanger and C. Hugonie, _$M_{W}$ , dark matter and $a_{\mu}$ in the NMSSM_, _Eur. Phys. J. C_ 82 (2022) 1074 [2209.03863].
  * (79) H. Bahl, T. Biekötter, S. Heinemeyer, C. Li, S. Paasch, G. Weiglein et al., _HiggsTools: BSM scalar phenomenology with new versions of HiggsBounds and HiggsSignals_ , _Comput. Phys. Commun._ 291 (2023) 108803 [2210.09332].
  * (80) P. Bechtle, D. Dercks, S. Heinemeyer, T. Klingl, T. Stefaniak, G. Weiglein et al., _HiggsBounds-5: Testing Higgs Sectors in the LHC 13 TeV Era_ , _Eur. Phys. J. C_ 80 (2020) 1211 [2006.06007].
  * (81) P. Bechtle, S. Heinemeyer, T. Klingl, T. Stefaniak, G. Weiglein and J. Wittbrodt, _HiggsSignals-2: Probing new physics with precision Higgs measurements in the LHC 13 TeV era_ , _Eur. Phys. J. C_ 81 (2021) 145 [2012.09197].
  * (82) G. Alguero, J. Heisig, C.K. Khosa, S. Kraml, S. Kulkarni, A. Lessa et al., _Constraining new physics with SModelS version 2_ , _JHEP_ 08 (2022) 068 [2112.00769].
  * (83) T. Hahn, _Generating Feynman diagrams and amplitudes with FeynArts 3_ , _Comput.Phys.Commun._ 140 (2001) 418 [hep-ph/0012260].
  * (84) T. Hahn and M. Perez-Victoria, _Automatized one loop calculations in four-dimensions and D-dimensions_ , _Comput. Phys. Commun._ 118 (1999) 153 [hep-ph/9807565].
  * (85) W. Beenakker, R. Hopker and M. Spira, _PROSPINO: A Program for the production of supersymmetric particles in next-to-leading order QCD_ , hep-ph/9611232.
  * (86) W. Beenakker, M. Klasen, M. Kramer, T. Plehn, M. Spira and P.M. Zerwas, _The Production of charginos / neutralinos and sleptons at hadron colliders_ , _Phys. Rev. Lett._ 83 (1999) 3780 [hep-ph/9906298].
  * (87) L. Shang and Y. Zhang, _EasyScan_HEP: a tool for connecting programs to scan the parameter space of physics models_ , 2304.03636.
  * (88) F. Domingo and S. Paßehr, _Curing tachyonic tree-level syndrome in NMSSM light-singlet scenarios_ , _Eur. Phys. J. C_ 82 (2022) 98 [2109.00585].
  * (89) P. Athron, J.-h. Park, D. Stöckinger and A. Voigt, _FlexibleSUSY—A spectrum generator generator for supersymmetric models_ , _Comput. Phys. Commun._ 190 (2015) 139 [1406.2319].
  * (90) P. Athron, M. Bach, D. Harries, T. Kwasnitza, J.-h. Park, D. Stöckinger et al., _FlexibleSUSY 2.0: Extensions to investigate the phenomenology of SUSY and non-SUSY models_ , _Comput. Phys. Commun._ 230 (2018) 145 [1710.03760].
  * (91) F. Staub, _From Superpotential to Model Files for FeynArts and CalcHep/CompHep_ , _Comput.Phys.Commun._ 181 (2010) 1077 [0909.2863].
  * (92) F. Staub, _Automatic Calculation of supersymmetric Renormalization Group Equations and Self Energies_ , _Comput.Phys.Commun._ 182 (2011) 808 [1002.0840].
  * (93) F. Staub, _SARAH 3.2: Dirac Gauginos, UFO output, and more_ , _Computer Physics Communications_ 184 (2013) pp. 1792 [1207.0906].
  * (94) F. Staub, _SARAH 4: A tool for (not only SUSY) model builders_ , _Comput.Phys.Commun._ 185 (2014) 1773 [1309.7223].
  * (95) K. Benakli, M. Goodsell, W. Ke and P. Slavich, _W boson mass in minimal Dirac gaugino scenarios_ , _Eur. Phys. J. C_ 83 (2023) 43 [2208.05867].
  * (96) Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * (97) P. Athron, C. Balazs, A. Cherchiglia, D.H.J. Jacob, D. Stöckinger, H. Stöckinger-Kim et al., _Two-loop prediction of the anomalous magnetic moment of the muon in the Two-Higgs Doublet Model with GM2Calc 2_ , _Eur. Phys. J. C_ 82 (2022) 229 [2110.13238].
  * (98) T.N. Dao, D.N. Le and M. Mühlleitner, _Leptonic anomalous magnetic and electric dipole moments in the CP-violating NMSSM with and without inverse seesaw mechanism_ , _Eur. Phys. J. C_ 82 (2022) 954 [2207.12618].
  * (99) F. Staub, P. Athron, U. Ellwanger, R. Gröber, M. Mühlleitner, P. Slavich et al., _Higgs mass predictions of public NMSSM spectrum generators_ , _Comput. Phys. Commun._ 202 (2016) 113 [1507.05093].
  * (100) P. Slavich et al., _Higgs-mass predictions in the MSSM and beyond_ , _Eur. Phys. J. C_ 81 (2021) 450 [2012.15629].
  * (101) S. Heinemeyer, W. Hollik, D. Stockinger, A.M. Weber and G. Weiglein, _Precise prediction for M(W) in the MSSM_ , _JHEP_ 08 (2006) 052 [hep-ph/0604147].