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

TL;DR

This paper calculates two-loop ${ m O}(\a_t^2)$ corrections to neutral Higgs boson masses in the CP-violating NMSSM, improving the precision of Higgs spectrum predictions in supersymmetric models.

Contribution

It provides the first implementation of ${ m O}(\a_t^2)$ corrections in the NMSSMCALC code, enhancing the accuracy of Higgs mass calculations in the CP-violating NMSSM.

Findings

Corrections add 4-7% to the Higgs mass with $ar{ ext{DR}}$ renormalization.

Corrections reduce the Higgs mass by 6-9% with OS renormalization.

Theoretical uncertainties are analyzed via scale and scheme variations.

Abstract

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

Tables

Table 1. Table 1 : Scan ranges for the NMSSM scan, all parameters varied independently between the given minimum and maximum values.

	$M_1$	$M_2$	$A_t$	$A_b$	$A_{\tau }$	$m_{{\tilde{Q}}_3}$	$m_{{\tilde{t}}_R}$	$m_{{\tilde{b}}_R}$	$m_{{\tilde{L}}_3}$	$m_{{\tilde{\tau }}_R}$	$M_{H^{\pm }}$	$A_{\kappa }$	${\mu }_{\text{eff}}$
	in TeV
min	0.4	0.4	-2.0	-2.0	-2.0	0.4	0.4	2.0	0.4	0.4	0.2	-2.0	0.2
max	1.0	1.0	2.0	2.0	2.0	3.0	3.0	3.0	3.0	3.0	1.0	2.0	0.3

Table 2. Table 2 : P1OS: Mass values in GeV and main components of the neutral Higgs bosons at tree-level, one-loop, two-loop 𝒪 ​ ( α t ​ α s ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 {\cal O}(\alpha_{t}\alpha_{s}) and at two-loop 𝒪 ​ ( α t ​ α s + α t 2 ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 superscript subscript 𝛼 𝑡 2 {\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2}) obtained by using OS renormalization in the top/stop sector.

	$H_1$	$H_2$	$H_3$	$H_4$	$H_5$
tree-level	74.29	91.43	704.12	895.91	897.83
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$
one-loop	86.58	135.0	700.03	895.83	897.83
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$
two-loop $𝒪\left({\alpha }_t{\alpha }_s\right)$	86.16	118.11	700.04	895.83	897.76
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$
two-loop $𝒪\left({\alpha }_t{\alpha }_s+{\alpha }_t^2\right)$	86.35	125.05	700.04	895.83	897.79
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$

Table 3. Table 3 : P1OS: Mass values in GeV and main components of the neutral Higgs bosons at tree-level, one-loop, two-loop 𝒪 ​ ( α t ​ α s ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 {\cal O}(\alpha_{t}\alpha_{s}) and at two-loop 𝒪 ​ ( α t ​ α s + α t 2 ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 superscript subscript 𝛼 𝑡 2 {\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2}) obtained using DR ¯ ¯ DR \overline{\text{DR}} renormalization in the top/stop sector.

	$H_1$	$H_2$	$H_3$	$H_4$	$H_5$
tree-level	74.29	91.43	704.12	895.91	897.83
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$
one-loop	85.93	112.77	700.05	895.79	897.71
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$
two-loop $𝒪\left({\alpha }_t{\alpha }_s\right)$	86.21	118.62	700.04	895.78	897.73
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$
two-loop $𝒪\left({\alpha }_t{\alpha }_s+{\alpha }_t^2\right)$	86.22	119.1	700.04	895.78	897.73
main component	$h_s$	$h_u$	$a_s$	$a$	$h_d$

Table 4. Table 4 : P2DR: Mass values in GeV and main components of the neutral Higgs bosons at tree-level, one-loop, two-loop 𝒪 ​ ( α t ​ α s ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 {\cal O}(\alpha_{t}\alpha_{s}) and at two-loop 𝒪 ​ ( α t ​ α s + α t 2 ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 superscript subscript 𝛼 𝑡 2 {\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2}) obtained by using DR ¯ ¯ DR \overline{\text{DR}} renormalization in the top/stop sector.

	$H_1$	$H_2$	$H_3$	$H_4$	$H_5$
tree-level	89.38	409.50	788.76	790.98	1828.56
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$
one-loop	120.86	407.68	788.64	791.01	1827.81
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$
two-loop $𝒪\left({\alpha }_t{\alpha }_s\right)$	124.58	407.69	788.65	791.0	1827.81
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$
two-loop $𝒪\left({\alpha }_t{\alpha }_s+{\alpha }_t^2\right)$	125.67	407.69	788.65	791.0	1827.81
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$

Table 5. Table 5 : P2DR: Mass values in GeV and main components of the neutral Higgs bosons at tree-level, one-loop, two-loop 𝒪 ​ ( α t ​ α s ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 {\cal O}(\alpha_{t}\alpha_{s}) and at two-loop 𝒪 ​ ( α t ​ α s + α t 2 ) 𝒪 subscript 𝛼 𝑡 subscript 𝛼 𝑠 superscript subscript 𝛼 𝑡 2 {\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2}) obtained using OS renormalization in the top/stop sector.

	$H_1$	$H_2$	$H_3$	$H_4$	$H_5$
tree-level	89.38	409.50	788.76	790.98	1828.56
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$
one-loop	142.91	407.74	788.62	790.9	1827.81
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$
two-loop $𝒪\left({\alpha }_t{\alpha }_s\right)$	123.92	407.71	788.57	790.91	1827.81
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$
two-loop $𝒪\left({\alpha }_t{\alpha }_s+{\alpha }_t^2\right)$	133.56	407.71	788.59	790.91	1827.81
main component	$h_u$	$a_s$	$h_d$	$a$	$h_s$

Table 6. Table 6 : Table adapted from Tab. 3 of [ 81 ] with the masses for the CP-even scalars (in GeV) for TP1–TP6 (defined in Tab. 2 of [ 81 ] ) when using the spectrum generators “out-of-the-box”. The values correspond to the two-loop results obtained by the different tools (the NMSSMTools value corresponds to one using the option leading to the most precise calculation implemented in NMSSMTools ). For NMSSMCALC the DR ¯ ¯ DR \overline{\text{DR}} scheme in the top/stop sector is applied in three different variants as defined in the text. The masses for the SM-like scalar are written in bold fonts, those for the singlet-like scalar in italics.

	TP1	TP2	TP3	TP4	TP5	TP6
	$h_1$
FlexibleSUSY	123.55	122.84	91.11	127.62	120.86	126.46
NMSSMCALC ’QCDrun+$𝒪(2,a)$’	120.34	118.57	90.88	126.37	120.32	123.45
NMSSMCALC ’EWrun+$𝒪(2,b)$’	123.58	121.51	90.99	127.38	120.82	124.89
NMSSMCALC ’gaugerun+$𝒪(2,b)$’	124.31	122.21	91.01	127.69	120.92	125.32
NMSSMTOOLS	123.52	121.83	90.78	127.30	119.31	126.63
SOFTSUSY	123.84	123.08	90.99	127.52	120.81	126.67
SPHENO	124.84	124.74	89.54	126.62	119.11	131.29
	$h_2$
FlexibleSUSY	1797.46	5951.36	126.58	143.11	125.08	700.80
NMSSMCALC ’QCDrun+$𝒪(2,a)$’	1797.45	5951.36	124.86	142.59	123.14	701.02
NMSSMCALC ’EWrun+$𝒪(2,b)$’	1797.45	5951.36	126.13	142.79	124.16	701.06
NMSSMCALC ’gaugerun+$𝒪(2,b)$’	1797.45	5951.36	126.51	142.84	124.51	701.06
NMSSMTOOLS	1797.46	5951.36	127.28	144.07	126.95	700.46
SOFTSUSY	1797.46	5951.36	126.59	143.02	125.12	701.01
SPHENO	1798.01	5951.35	126.77	144.04	125.61	689.30
	$h_3$
FlexibleSUSY	2758.96	6372.08	652.95	467.80	627.28	1369.53
NMSSMCALC ’QCDrun+$𝒪(2,a)$’	2756.70	6371.48	652.58	467.48	627.10	1368.08
NMSSMCALC ’EWrun+$𝒪(2,b)$’	2756.70	6368.58	652.70	467.73	627.16	1368.95
NMSSMCALC ’gaugerun+$𝒪(2,b)$’	2756.70	6368.47	652.67	467.72	627.15	1368.91
NMSSMTOOLS	2758.51	6345.72	651.03	466.38	623.79	1368.90
SOFTSUSY	2758.41	6370.3	652.78	467.73	627.14	1369.19
SPHENO	2757.11	6366.88	651.21	467.5	624.02	1363.02

Equations

$$W_{\text{NMSSM}}={ϵ}_{ij}[y_e{\hat{H}}_d^i{\hat{L}}^j{\hat{E}}^c+y_d{\hat{H}}_d^i{\hat{Q}}^j{\hat{D}}^c-y_u{\hat{H}}_u^i{\hat{Q}}^j{\hat{U}}^c]-{ϵ}_{ij}\lambda \hat{S}{\hat{H}}_d^i{\hat{H}}_u^j+\frac{1}{3}\kappa {\hat{S}}^3.$$

$${\mathcal{L}}_{\text{soft},\text{ NMSSM}}=$$

$$-m_{{\tilde{e}}_R}^2{\tilde{e}}_R^{\ast }{\tilde{e}}_R-({ϵ}_{ij}[y_e{A}_e{H}_d^i{\tilde{L}}^j{\tilde{e}}_R^{\ast }+y_d{A}_d{H}_d^i{\tilde{Q}}^j{\tilde{d}}_R^{\ast }-y_u{A}_u{H}_u^i{\tilde{Q}}^j{\tilde{u}}_R^{\ast }]+h.c.)$$

$$-\frac{1}{2}(M_1\tilde{B}\tilde{B}+M_2{\tilde{W}}_i{\tilde{W}}_i+M_3\tilde{G}\tilde{G}+h.c.)$$

$$-m_S^2|S{|}^2+({ϵ}_{ij}\lambda A_{\lambda }S{H}_d^i{H}_u^j-\frac{1}{3}\kappa A_{\kappa }{S}^3+h.c.).$$

$$V_H$$

$$H_d=\left(\begin{array}{c}\frac{1}{\sqrt{2}}(v_d+h_d+i{a}_d)\\ h_d^{-}\end{array}\right),H_u=e^{i{\phi }_u}\left(\begin{array}{c}{h}_u^{+}\\ \frac{1}{\sqrt{2}}(v_u+h_u+i{a}_u)\end{array}\right),S=\frac{e^{i{\phi }_s}}{\sqrt{2}}(v_s+h_s+i{a}_s).$$

$$\tan \beta =\frac{v_u}{v_d}.$$

$$t_L\to e^{-i{\phi }_u/2}{t}_L\text{mbox}and{t}_R\to e^{i{\phi }_u/2}{t}_R.$$

$$V_H$$

$$t_{\varphi }\equiv {\frac{\partial V_{\text{Higgs}}}{\partial \varphi }|}_{\text{Min.}}=0$$

$$\frac{t_{h_d}}{v{c}_{\beta }}$$

$$+\frac{1}{2}|\lambda {|}^2\left({s_{\beta }}^2{v}^2+v_s^2\right)$$

$$\frac{t_{h_u}}{v{s}_{\beta }}$$

$$+\frac{1}{2}|\lambda {|}^2\left({c_{\beta }}^2{v}^2+v_s^2\right)$$

$$\frac{t_{h_s}}{v_S}$$

$$+\frac{|\kappa |v_s(\text{Re}{A}_{\kappa }{c}_{{\phi }_{\omega }}-\text{Im}{A}_{\kappa }{s}_{{\phi }_{\omega }})}{\sqrt{2}}$$

$$\frac{t_{a_d}}{v{s}_{\beta }}$$

$$t_{a_u}$$

$$t_{a_s}$$

$$-\frac{|\kappa |v_s^2(\text{Im}{A}_{\kappa }{c}_{{\phi }_{\omega }}+\text{Re}{A}_{\kappa }{s}_{{\phi }_{\omega }})}{\sqrt{2}},$$

$${\phi }_y$$

$${\phi }_{\omega }$$

$${\mathcal{M}}_{h^{+}{h}^{-}}$$

$$+\left(\begin{array}{ccc}\frac{t_{h_d}-t_{a_d}\tan ({\phi }_{\omega }-{\phi }_y)}{v{c}_{\beta }}& & -\frac{t_{a_d}(\tan ({\phi }_{\omega }-{\phi }_y)+i)}{s_{\beta }v}\\ -\frac{t_{a_d}(\tan ({\phi }_{\omega }-{\phi }_y)-i)}{s_{\beta }v}& & \frac{s_{\beta }{t}_{h_u}-c_{\beta }{t}_{a_d}\tan ({\phi }_{\omega }-{\phi }_y)}{{s_{\beta }}^{2}v}\end{array}\right)$$

$$(h_d,h_u,h_s,a,a_s,G{)}^T$$

$$(h_1,h_2,h_3,h_4,h_5,G{)}^T$$

$$\text{diag}(m_{h_1}^2,m_{h_2}^2,m_{h_3}^2,m_{h_4}^2,m_{h_5}^2,0)$$

$${\mathcal{R}}^G=\left(\begin{array}{cc}\text{mathbbm}{1}_{3\times 3}& 0\\ 0& {\tilde{\mathcal{R}}}^G\end{array}\right),{\tilde{\mathcal{R}}}^G=\left(\begin{array}{ccc}{s}_{{\beta }_n}& c_{{\beta }_n}& 0\\ 0& 0& 1\\ c_{{\beta }_n}& -s_{{\beta }_n}& 0\end{array}\right).$$

$$\left(\begin{array}{c}{G}^{-}\\ H^{-}\end{array}\right)={\mathcal{R}}^{G^{-}}\left(\begin{array}{c}{h}_d^{-}\\ h_u^{-}\end{array}\right),\text{diag}(M_{H^{\pm }}^2,0)={\mathcal{R}}^{G^{-}}{\mathcal{M}}_{h^{+}{h}^{-}}({\mathcal{R}}^{G^{-}}{)}^T,$$

$${\mathcal{R}}^{G^{-}}=\left(\begin{array}{cc}-\cos {\beta }_c& \sin {\beta }_c\\ \sin {\beta }_c& \cos {\beta }_c\end{array}\right).$$

$$M_{H^{\pm }}^2=$$

$$+\frac{s_{\beta }\left(c_{\beta }{c}_{{\beta }_c}^2{t}_{h_u}+s_{\beta }{s}_{{\beta }_c}^2{t}_{h_d}\right)+c_{\beta -{\beta }_c}^2{t}_{a_d}\tan ({\phi }_y-{\phi }_w)}{c_{\beta }{s_{\beta }}^{2}v},$$

$${\mu }_{\text{eff}}=\frac{\lambda v_s{e}^{i{\phi }_s}}{\sqrt{2}},$$

$$\left\{t_{h_d},t_{h_u},t_{h_s},t_{a_d},t_{a_s},M_{H^{\pm }}^2,v,s_{{\theta }_W},e,\tan \beta ,|\lambda |,v_s,|\kappa |,\text{Re}{A}_{\kappa },{\phi }_{\lambda },{\phi }_{\kappa },{\phi }_u,{\phi }_s\right\},$$

$$\left\{t_{h_d},t_{h_u},t_{h_s},t_{a_d},t_{a_s},v,s_{{\theta }_W},e,\tan \beta ,|\lambda |,v_s,|\kappa |,\text{Re}{A}_{\lambda },\text{Re}{A}_{\kappa },{\phi }_{\lambda },{\phi }_{\kappa },{\phi }_u,{\phi }_s\right\}.$$

$${\mathcal{M}}_{\tilde{t}}=\left(\begin{array}{cc}{m}_{{\tilde{Q}}_3}^2+m_t^2& m_t\left(A_t^{\ast }{e}^{-i{\phi }_u}-\frac{{\mu }_{\text{eff}}}{\tan \beta }\right)\\ m_t\left(A_t{e}^{i{\phi }_u}-\frac{{\mu }_{\text{eff}}^{\ast }}{\tan \beta }\right)& m_{{\tilde{t}}_R}^2+m_t^2\end{array}\right),$$

$$({\tilde{t}}_1,{\tilde{t}}_2{)}^T$$

Full text

1 Introduction

Supersymmetric theories [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] belong to the best motivated and most intensively studied extensions of the Standard Model (SM). Supersymmetry (SUSY) between bosonic and fermionic degrees of freedom solves the hierarchy problem and relates the Higgs boson masses and self-couplings to the electroweak gauge couplings. The tree-level mass value of the SM-like Higgs boson is therefore bound to be of the order of the electroweak scale. In the Minimal Supersymmetric extension of the SM (MSSM) [15, 16, 17, 18] it is less than or equal to the $Z$-boson mass. This upper bound is lifted to higher values after the inclusion of the radiative corrections to the Higgs boson masses. In order to match the measured mass value of 125 GeV of the discovered Higgs boson large values of the stop masses and/or mixing are required. This challenges the model from the point of view of fine-tuning. The situation is more relaxed in the Next-to-Minimal SUSY extension (NMSSM) [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. The introduction of a complex singlet superfield coupling with the strength $\lambda$ to the two Higgs doublet superfields of the MSSM, adds new contributions to the quartic coupling so that the tree-level mass of the lightest CP-even MSSM-like Higgs boson is shifted to a higher value.

The Higgs boson mass is a crucial input parameter in all Higgs boson observables (see e.g. [31]) and plays an important role for the stability of the electroweak vacuum [32, 33, 34]. A precise measurement of the Higgs boson mass (the current experimental value is $125.09\pm 0.24$ GeV [35]) is hence indispensable. The experimental accuracy has to be matched by the precision of the theory predictions in order to fully exploit the experimental information to constrain the still viable parameter space of beyond-the-SM (BSM) extensions and to distinguish between new physics models in case of discovery.666For a discussion see e.g. Ref. [36] comparing various SUSY and non-SUSY extensions of the Higgs sector. In the recent years there has been a lot of activity, and is still ongoing, in order to improve the precision of the Higgs mass predictions within the NMSSM, both for the CP-conserving and the CP-violating case. In the CP-conserving NMSSM, after the computation of the leading one-loop contributions to the Higgs boson masses in [37, 38, 39, 40, 41, 42], the full one-loop corrections including the momentum dependence became available in the $\overline{\mbox{DR}}$ [43, 44] as well as a mixed $\overline{\text{DR}}$-OS scheme [45, 46].777We note that a full one-loop renormalization of the NMSSM has also been presented in [47, 48]. Two-loop corrections at the order ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{b}\alpha_{s})$ in the $\overline{\mbox{DR}}$ scheme have been provided in [43] in the approximation of zero external momentum. Two-loop corrections beyond these have been derived in [49], again exploiting the $\overline{\mbox{DR}}$ scheme. In the CP-violating case first results for the dominant one-loop corrections [50, 51, 52, 53, 54] have been followed by the full one-loop and logarithmically enhanced two-loop effects computed in the renormalization group approach [55]. Our group complemented these calculations by computing the full one-loop corrections in the Feynman diagrammatic approach [56] (see also [57], published recently). We subsequently calculated the ${\cal O}{(\alpha_{t}\alpha_{s})}$ corrections in the approximation of vanishing external momentum in the gaugeless limit [58].888We also provided the complete one-loop corrections to the trilinear Higgs self-couplings in the CP-conserving NMSSM [59] and extended these in [60] to the ${\cal O}(\alpha_{t}\alpha_{s})$ corrections in the CP-violating NMSSM. The adopted renormalization scheme is a mixture between $\overline{\text{DR}}$ and OS conditions with the possibility to choose between the $\overline{\text{DR}}$ or OS scheme for the renormalization of the top/stop sector. An automatized two-loop calculation of the Higgs boson masses in CP-violating SUSY theories, applying the $\overline{\text{DR}}$ scheme, was provided by [61].

The higher-order corrections to the NMSSM Higgs boson masses have been implemented in publicly available tools that partly also compute the Higgs boson decays. The program package NMSSMTools [62, 63, 64] calculates the masses and decay widths in the CP-conserving $\mathbb{Z}^{3}$ and can be interfaced with SOFTSUSY [65, 66], which includes the possibility of $\mathbb{Z}^{3}$ violation. Recently, a tool became available to support the extension of NMSSMTools to the CP-violating case [67, 57, 68]. The spectrum generation of different SUSY models, including the NMSSM, on the other hand, is possible through the interface of SARAH [69, 70, 71, 72, 49] with SPheno [73, 74]. In the same spirit, SARAH has been interfaced with the package FlexibleSUSY [75, 76]. All these programs include the Higgs mass corrections up to two-loop order, where in particular the two-loop corrections are obtained in the effective potential approach.999With the exception of the recent NMSSMTools extension to the complex NMSSM, which relies on a Feynman diagrammatic calculation in the OS scheme [57]. In FlexibleEFTHiggs [77], an effective field theory approach has been combined with a diagrammatic calculation to obtain the SM-like Higgs pole mass in various models, including the NMSSM. Our program package NMSSMCALC [78] computes the NMSSM Higgs boson masses and decay widths in the Feynman diagrammatic approach both for the CP-conserving and CP-violating case.101010One-loop corrected decay widths are also included in the code SloopS [47, 48, 79]. In the SARAH and SPheno framework a generic implementation of the two-body partial decays widths exists at the full one-loop level [80]. It incorporates the complete one-loop corrections at non-zero external momentum and the two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ corrections in the limit of vanishing external momentum in the gaugeless approximation. The renormalization is performed in a mixed $\overline{\text{DR}}$-OS scheme with the possibility to choose between the $\overline{\text{DR}}$ or OS scheme for the renormalization of the top/stop sector. For the CP-conserving NMSSM, a detailed comparison of the results of the various tools for the Higgs boson masses in the $\overline{\text{DR}}$ scheme was performed in [81]. In [82], a comparison of the Higgs boson masses and mixing matrices in the OS scheme up to ${\cal O}(\alpha_{t}\alpha_{s})$ was performed between the two codes NMSSMCALC and NMSSM-FeynHiggs.

With this paper we take another step in improving our predictions for the NMSSM Higgs boson masses. We compute the two-loop corrections to the neutral NMSSM Higgs boson masses in the Feynman diagrammatic approach at the order ${\cal O}(\alpha_{t}^{2})$ at vanishing external momentum in the gaugeless limit. Strictly speaking, we provide the contributions of ${\cal O}(m_{t}^{2}\alpha_{t}^{2})$ neglecting further terms beyond the approximation of the gaugeless limit and vanishing external momentum. For simplicity, we call our newly calculated contributions ${\cal O}(\alpha_{t}^{2})$ terms in the following. The calculation is performed both for the CP-conserving and the CP-violating case. We apply a mixed $\overline{\text{DR}}$-OS renormalization scheme. Our corrections are included in the program package NMSSMCALC.111111The program package can be downloaded from the url: https://www.itp.kit.edu/$\sim$maggie/NMSSMCALC/

The paper is organized as follows. In Section 2 we introduce the model and set our notation. Section 3 contains the detailed presentation of our calculation including the description of the renormalization and the checks that were performed. Section 4 is dedicated to the numerical presentation of our results. We conclude in Section 5. We include an extensive appendix that contains the presentation of the two-loop self-energy diagrams, the scalar one-loop integrals expanded up to ${\cal O}(\varepsilon)$, the details on the computation of the running $\overline{\mbox{DR}}$ top quark mass without and with the inclusion of the gauge coupling contributions, the required NMSSM renormalization group equations for the investigation of the scale dependence of our results, the neutral tree-level and counterterm mass matrices as well as the charged Higgs mass counterterms.

2 The CP-Violating NMSSM

In order to introduce the model and set our notation we briefly review the Higgs sector of the CP-violating NMSSM. For more details on the NMSSM, see the reviews [30, 29]. Since the two-loop diagrams appearing in our calculation also involve contributions from the stops and the electroweakinos, we briefly introduce them in this section as well.

2.1 The Higgs Sector

The framework for our computation is the CP-violating NMSSM with a scale-invariant superpotential and a discrete $\mathbb{Z}^{3}$ symmetry. The Higgs potential is derived from the superpotential, the soft SUSY breaking Lagrangian and the $D$-term contributions. The superpotential reads, in terms of the two Higgs doublet superfields $\hat{H}_{d}$ and $\hat{H}_{u}$, the singlet superfield $\hat{S}$, the quark and lepton superfields and their charge conjugates, denoted by the superscript $c$, $\hat{Q}$, $\hat{U}^{c}$, $\hat{D}^{c}$, $\hat{L}$, $\hat{E}^{c}$,

[TABLE]

Here $\epsilon_{ij}$ ($i,j=1,2$) is the totally antisymmetric tensor with $\epsilon_{12}=\epsilon^{12}=1$ and $i,j$ denoting the indices of the fundamental $SU(2)_{L}$ representation. Here and in the following, we sum over repeated indices. For simplicity, we neglect colour and generation indices. We neglect flavour mixing so that the Yukawa couplings $y_{e},y_{d}$ and $y_{u}$ are diagonal, and therefore complex phases can be reabsorbed by redefining the quark fields without effect on the physical meaning [83]. The dimensionless NMSSM-specific parameters $\lambda$ and $\kappa$ are complex in the CP-violating NMSSM. For the soft SUSY breaking Lagrangian we have in terms of the scalar component fields $H_{u}$, $H_{d}$ and $S$,

[TABLE]

The summation over all three quark and lepton generations is implicit. By $\tilde{Q}$ and $\tilde{L}$ we denote the complex scalar components of the corresponding quark and lepton superfields, so that we have for the first generation e.g. $\tilde{Q}=(\tilde{u}_{L},\tilde{d}_{L})^{T}$ and $\tilde{L}=(\tilde{\nu}_{L},\tilde{e}_{L})^{T}$. Working in the CP-violating case, the soft SUSY breaking trilinear couplings $A_{x}$ ($x=\lambda,\kappa,d,u,e$) and the gaugino mass parameters $M_{k}$ ($k=1,2,3$) of the bino, wino and gluino fields $\tilde{B}$, $\tilde{W}_{i}$ ($i=1,2,3$) and $\tilde{G}$, are complex. Application of the $R$-symmetry transformation allows to choose either $M_{1}$ or $M_{2}$ to be real. However, we keep both $M_{1}$ and $M_{2}$ complex in NMSSMCALC. The soft SUSY breaking mass parameters of the scalar fields, $m_{X}^{2}$ ($X=S,H_{d},H_{u},\tilde{Q},\tilde{u}_{R},\tilde{d}_{R},\tilde{L},\tilde{e}_{R}$), are real. The final Higgs potential at tree level reads

[TABLE]

where $g_{1}$ and $g_{2}$ denote the $U(1)_{Y}$ and $SU(2)_{L}$ gauge couplings, respectively. Two more CP-violating phases, $\varphi_{u}$ and $\varphi_{s}$, are introduced in the expansion of the two Higgs doublets and the singlet field about their vacuum expectation values (VEVs) $v_{d}$, $v_{u}$ and $v_{s}$,

[TABLE]

The VEV $v\approx 246$ GeV is related to $v_{u}$ and $v_{d}$ through $v^{2}=v_{d}^{2}+v_{u}^{2}$, and the ratio between the two VEVs is denoted by $\tan\beta$,

[TABLE]

The phase $\varphi_{u}$ enters the top quark mass term.121212For simplicity we focus on the top quark mass but the discussion holds for all up-type quark mass terms. However, we keep the coupling $y_{t}\exp(i\varphi_{u})$, appearing in the mass term, real by absorbing this phase into the left- and right-handed top quark fields through the replacements

[TABLE]

This affects all couplings involving one top quark. After inserting Eq. (2.4) into Eq. (2.1) the Higgs potential can be cast into the form

[TABLE]

with $\phi^{0}\equiv(h_{d},h_{u},h_{s},a_{d},a_{u},a_{s})^{T}$ and $\phi^{c}\equiv((h_{d}^{-})^{*},h_{u}^{+})^{T}$. The six tadpole coefficients $t_{\phi}$ ($\phi=h_{d},h_{u},h_{s},a_{d},a_{u},a_{s}$) are defined via

[TABLE]

at tree-level. This definition yields the following tree-level relations for the tadpole coefficients,

[TABLE]

with $M_{Z}$ being the $Z$ boson mass. Here and in the following we use the short-hand notations $c_{x}\equiv\cos x$ and $s_{x}\equiv\sin x$, and we express the complex parameters $A_{\lambda}$ and $A_{\kappa}$ by their imaginary and real parts in order to comply with the SUSY Les Houches Accord (SLHA) [84, 85] as well as $\lambda$ and $\kappa$ by their absolute values and their phases. The phases enter in two combinations together with $\varphi_{u}$ and $\varphi_{s}$,

[TABLE]

Since $t_{a_{u}}$ and $t_{a_{d}}$ are linearly dependent, only five out of the six tadpole equations yield linearly independent conditions in the Higgs sector. While the tadpole coefficients vanish at tree level due to the minimization conditions of the Higgs potential, they affect the higher-order corrections through the appearance of tadpole counterterms. We therefore keep the tadpole coefficients explicitly in all quantities that need to be renormalized and set them to zero only after the renormalization is performed. Note that at lowest order the two parameters $\text{Im}\,A_{\lambda}$ and $\text{Im}\,A_{\kappa}$ can be eliminated by using the minimisation conditions $t_{a_{d}}=0$ and $t_{a_{s}}=0$, and the soft-breaking mass parameters $m_{H_{d}}^{2}$, $m_{H_{u}}^{2}$ and $m_{S}^{2}$ can be eliminated by using $t_{h_{d}}=0$, $t_{h_{u}}=0$ and $t_{h_{s}}=0$.

The $6\times 6$ mass matrix for the neutral Higgs bosons is denoted by ${\mathcal{M}}_{\phi\phi}$ and the $2\times 2$ mass matrix for the charged Higgs bosons by ${\mathcal{M}}_{h^{+}h^{-}}$. The explicit expressions for the mass matrix ${\mathcal{M}}_{\phi\phi}$ can be found in App. F, while ${\mathcal{M}}_{h^{+}h^{-}}$ is given by

[TABLE]

where we explicitly keep the dependence on the tadpole parameters. Constant terms and trilinear and quartic interactions are summarized in $V_{H}^{\mbox{\scriptsize const}}$ and $V_{H}^{\phi^{3},\phi^{4}}$, respectively. The mass eigenstates $h_{i}$ ($i=1,...,5$) are obtained by rotating from the interaction to the mass basis with two consecutive rotations, where the first rotation ${\cal R}^{G}$ singles out the would-be Goldstone boson, and the second one, ${\cal R}$, performs the rotation to the mass eigenstates,

[TABLE]

with the diagonal mass matrix

[TABLE]

and

[TABLE]

The mass eigenstates $h_{i}$ are ordered by ascending mass, i.e. $m_{h_{1}}\leq...\leq m_{h_{5}}$. At tree level, the rotation angle $\beta_{n}$ is equal to the angle $\beta$ defining the ratio of the VEVs. However, we distinguish them here because in the subsequently applied renormalization procedure only $\beta$ receives a counterterm but not $\beta_{n}$. The charged Higgs and Goldstone boson $H^{-}$ and $G^{-}$, respectively, are obtained through the rotation

[TABLE]

with

[TABLE]

Like the mixing angle $\beta_{n}$, the rotation angle $\beta_{c}$ of the charged sector is considered to be already renormalized and hence does not receive a counterterm. At tree level $\beta_{c}=\beta_{n}=\beta$ holds. The relation between the charged Higgs mass and $\text{Re}\,A_{\lambda}$ is given by

[TABLE]

where $M_{W}$ is the $W$ boson mass. Note that we explicitly kept the difference between $\beta$ and $\beta_{c}$ as well as the tadpole parameters, all of which would vanish at tree level, in the formula for the charged Higgs boson mass since its counterterm at one- and two-loop order receives contributions from the counterterms of $\beta$ and of the tadpole parameters.

The MSSM limit of the complex NMSSM is obtained by $\lambda,\kappa\to 0$ (and constant ratio of $\kappa/\lambda$ for a smooth approximation of the limit) and keeping $A_{\lambda}$, $A_{\kappa}$ and the effective $\mu_{\text{eff}}$ parameter,

[TABLE]

fixed.

We choose the set of independent parameters entering the Higgs potential to be

[TABLE]

where $\theta_{W}$ denotes the weak mixing angle. Alternatively, we can choose to use $\text{Re}\,A_{\lambda}$ as input and calculate $M_{H^{\pm}}^{2}$ by means of Eq. (2.23). In this case, the set of independent parameters entering the Higgs potential is given by

[TABLE]

2.2 The Squark Sector

The top, bottom, stop and sbottom particles appear in the $\mathcal{O}(\alpha_{t}^{2})$ two-loop diagrams of the self-energies of the neutral and charged Higgs bosons. The two-loop corrections to the charged Higgs boson need to be computed for the definition of the two-loop counterterm of the squared charged Higgs mass, since this counterterm explicitly appears in the $\mathcal{O}(\alpha_{t}^{2})$ two-loop corrections to the masses of the neutral Higgs bosons. In the gaugeless approximation that we apply in our calculation, i.e. $e\to 0$ (or, in other words, $g_{1}\to 0$ and $g_{2}\to 0$), the stop mass matrix reads

[TABLE]

where $m_{t}$ denotes the top quark mass. The matrix is diagonalized by a unitary matrix $\mathcal{U}_{\tilde{t}}$, rotating the interaction states $\tilde{t}_{L}$ and $\tilde{t}_{R}$ to the mass eigenstates $\tilde{t}_{1}$ and $\tilde{t}_{2}$,

[TABLE]

Throughout our calculation of the $\mathcal{O}(\alpha_{t}^{2})$ two-loop corrections, we set the bottom quark mass to zero, i.e. $m_{b}=0$. Therefore, the left- and right-handed sbottom states will not mix and only the left-handed sbottom with a mass of $m_{\tilde{Q}_{3}}$ will contribute. With the parameters $\varphi_{u},\tan\beta$ and those defining $\mu_{\text{eff}}$ already appearing in the Higgs sector, we are hence left with the following set of independent parameters for the third generation quark/squark sector

[TABLE]

With this parameter choice for the mass matrix in the interaction basis, the rotation matrix $\mathcal{U}_{\tilde{t}}$ does not need to receive a counterterm. We follow here the same approach as in the Higgs sector, where we assumed the rotation matrices to be renormalized already.

2.3 The Chargino and Neutralino Sectors

The computation of the ${\cal O}(\alpha_{t}^{2})$ contributions to the Higgs self-energies and tadpoles also involves the neutralino and chargino sectors. In the gaugeless approximation, electroweak gauginos do not mix with the fermionic superpartners of the Higgs bosons and do not contribute to the $\mathcal{O}(\alpha_{t}^{2})$ correction. Only the Higgsinos enter the two-loop diagrams of the $\mathcal{O}(\alpha_{t}^{2})$ neutral and charged Higgs self-energies. The neutralino mass matrix in the Weyl spinor basis $\psi^{0}=(\tilde{B},\tilde{W}_{3},\tilde{H}^{0}_{d},\tilde{H}^{0}_{u},\tilde{S})^{T}$ simplifies to

[TABLE]

with the neutral $2\times 2$ gaugino mass matrix $M_{G}$ given by

[TABLE]

whose eigenstates decouple and do not contribute to the $\mathcal{O}(\alpha_{t}^{2})$ results in the gaugeless limit in our calculation. The neutral $3\times 3$ Higgsino mass matrix $M_{N}$ reads

[TABLE]

The symmetric matrix $M_{N}$ can be diagonalized by a $3\times 3$ matrix $N$, yielding

[TABLE]

where we denoted the Higgsino mass eigenstates by $\tilde{\chi}_{3,4,5}^{0}$. The neutral Higgsino mass eigenstates, expressed as a Majorana spinor ($i=3,4,5$)

[TABLE]

are obtained from

[TABLE]

They are ordered by ascending mass, i.e. $|m_{\tilde{\chi}^{0}_{3}}|\leq|m_{\tilde{\chi}^{0}_{4}}|\leq|m_{\tilde{\chi}^{0}_{5}}|$.

In the gaugeless approximation, the mass matrix of the charginos in the interaction basis $\psi_{R}^{-}=(\tilde{W}^{-},\tilde{H}^{-}_{d})^{T}$, respectively, $\psi_{L}^{+}=(\tilde{W}^{+},\tilde{H}^{+}_{u})^{T}$, simplifies to

[TABLE]

and is equal to the mass matrix $M_{\chi^{\pm}}$ in the basis of the mass eigenstates $\chi^{\pm}$. The fermionic superpartner of the charged Higgs boson has a mass of

[TABLE]

Here we reabsorbed the phase $\mu_{\text{eff}}$ into the mixing matrices, hence it implicitly appears in the charged Higgsino couplings. The other superpartner with mass $M_{2}$ decouples and does not give any contribution to our calculation at $\mathcal{O}(\alpha_{t}^{2})$.

Since vertices and propagators with charginos and neutralinos enter at the two-loop level only, they do not need to be renormalized. In the one-loop inserted counterterms the effective Higgsino-mass parameter $\mu_{\text{eff}}$ enters via the couplings of Higgs bosons to stops. The parameter is defined in terms of $\lambda$ and $v_{s}$ whose renormalization is specified in the Higgs sector already, and therefore no further renormalization conditions need to be specified.

3 The Higher-Order NMSSM Higgs Boson Masses

The loop-corrected Higgs boson masses are given by the real parts of the poles of the propagators of the Higgs bosons. They are obtained numerically as the zeroes of the determinant

[TABLE]

of the renormalized two-point correlation function

[TABLE]

The matrix part of the renormalized two-point correlation function is given by the tree-level masses $m_{h_{i}}$ of the neutral Higgs bosons and the renormalized self-energies $\hat{\Sigma}_{ij}(p^{2})$ for the $h_{i}\to h_{j}$ transition of the tree-level mass eigenstates $h_{i,j}$ at $p^{2}$. It reads

[TABLE]

and contains all contributions of one- and two-loop order, where for the latter only the $\mathcal{O}(\alpha_{t}^{2})$ and $\mathcal{O}(\alpha_{t}\alpha_{s})$ corrections in the gaugeless limit and at vanishing external momentum are taken into account.

The numerical recipe to extract the zeroes of the determinant proceeds along the same lines as in [45, 56], which is based on an iterative procedure. While such an iterative procedure automatically mixes different orders of perturbation theory and is thus not of strict $\mathcal{O}(\alpha_{t}^{2})$ any more in our present case131313Moreover, the iterative procedure can introduce intricate gauge dependences by mixing different orders of perturbation theory. Since our order $\mathcal{O}(\alpha_{t}^{2})$ contributions are calculated in the gaugeless limit, however, additional gauge dependences are only introduced beyond the considered two-loop corrections through the iterative procedure., it nevertheless gives an improvement of the numerical results when compared to a fixed-order calculation, as argued for the one-loop case in [86]. In the first iteration step for the calculation of the $n^{\text{th}}$ neutral Higgs boson mass, the square of the external momentum in the renormalized self-energies is chosen to be at its tree-level mass shell, i.e. $p^{2}=m_{h_{n}}^{2}$. The mass matrix part of the two-point correlation function, Eq. (3.49), is then diagonalized, yielding the $n^{\text{th}}$ eigenvalue as a first approximation to the squared mass of the $n^{\text{th}}$ neutral Higgs boson. In the next step of the iteration, this value is taken as the new OS value for $p^{2}$, and the matrix part of $\hat{\Gamma}$ is again diagonalized. This procedure is repeated until the $n^{\text{th}}$ eigenvalue changes by less than $10^{-9}$ between two consecutive steps of the iteration. This iterative algorithm is applied for the calculation of all five neutral Higgs boson masses.

The renormalized Higgs boson self-energies $\hat{\Sigma}_{ij}$ consist of one-loop and two-loop contributions labeled by the superscript $(1)$ and $(2)$, respectively,

[TABLE]

In Refs. [45, 56] we computed the complete one-loop corrections to the neutral NMSSM Higgs bosons in the CP-conserving and CP-violating NMSSM, respectively. For details, we refer to these papers. The renormalized two-loop self-energies $\hat{\Sigma}_{ij}(0)$, which are evaluated in the approximation of vanishing external momentum, i.e. $p^{2}=0$, comprise the ${\cal O}(\alpha_{t}\alpha_{s})$ corrections, which we computed in [58], and the ${\cal O}(\alpha_{t}^{2})$ contributions computed in this paper,

[TABLE]

In the following, we concentrate on the $\mathcal{O}(\alpha_{t}^{2})$ corrections and suppress the superscript $\alpha_{t}^{2}$. Products of one-loop terms of ${\cal O}(\alpha_{t})$ give also contributions to the renormalized Higgs self-energies at $\mathcal{O}(\alpha_{t}^{2})$. Therefore, if not stated explicitly otherwise, the following ${\cal O}(\alpha_{t})$ contributions are also evaluated in the gaugeless limit and at zero external momentum, just as the $\mathcal{O}(\alpha_{t}^{2})$ contributions141414Note, however, that in the pure one-loop corrections to the Higgs boson masses we do not apply these approximations.. Although the renormalized Higgs self-energies are evaluated in the approximation $p^{2}=0$, in the following formulae we keep the momentum dependence for the purpose of introducing our notation.

The renormalized one-loop Higgs self-energy for the transition $h_{i}\to h_{j}$ ($i,j=1,\ldots,5$) at ${\cal O}(\alpha_{t})$ can be decomposed as

[TABLE]

where $\Sigma^{\text{\tiny{(1)}}}_{ij}(p^{2})$ denotes the unrenormalized self-energy at ${\cal O}(\alpha_{t})$. The remaining part is the one-loop counterterm consisting of the wave function renormalization constant matrix $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\mathcal{Z}$ and the mass counterterm $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!{\cal M}_{hh}$. At two-loop $\mathcal{O}(\alpha_{t}^{2})$, the renormalized self-energy is given by

[TABLE]

where $\Sigma^{(2)}_{ij}(p^{2})$ is the unrenormalized two-loop self-energy which is evaluated at $p^{2}=0$ and $\left(\delta^{\text{\tiny(}\!\text{\tiny 2}\!\text{\tiny)}}\!M^{2}\right)_{ij}$ is the $\mathcal{O}(\alpha_{t}^{2})$ mass counterterm given by

[TABLE]

The Higgs field renormalization constant matrix is diagonal and can be expressed as

[TABLE]

in terms of the renormalization constants $\Delta^{(n)}Z_{H_{u},H_{d},S}$ for the doublet and singlet fields. The definition of $\Delta^{(n)}Z_{H_{u},H_{d},S}$ is given subsequently in Eq. (3.97) for $n=1$ and Eq. (3.74) for $n=2$, respectively. We note once more that in the above formulae we kept the momentum dependence for the purpose of defining the Higgs field renormalization constants. Our results are obtained, however, in the approximation of vanishing external momentum, *i.e. *for $p^{2}=0$.

We remind the reader that in the mass matrix ${\cal M}_{hh}$, and hence also in its counterterm matrix, we have dropped the Goldstone component, cf. Eq. (2.19) with the rotation matrix ${\cal R}$ defined in Eq. (2.18). Hence, we neglect higher-order corrections due to mixing of the Goldstone boson with the remaining neutral Higgs bosons, which we have verified to be numerically negligible. The mass matrix counterterms $\delta{\cal M}_{hh}^{(1,2)}$ implicitly contains the counterterms of the parameters that need to be renormalized. We explicitly specify the renormalization of these parameters and the Higgs wave function renormalization constants in Secs. 3.2 and 3.3.

3.1 The Unrenormalized Self-Energies of the Neutral Higgs Bosons

The unrenormalized self-energies $\Sigma^{(n)}_{ij}(p^{2})$ of the neutral Higgs bosons consist of all Feynman diagrammatic contributions that are relevant at $\mathcal{O}(\alpha_{t})$ and $\mathcal{O}(\alpha_{t}^{2})$, respectively. The one-loop diagrams contributing at $\mathcal{O}(\alpha_{t})$ are shown in Fig. 1, while the two-loop diagrams of $\mathcal{O}(\alpha_{t}^{2})$ are presented in Fig. 14 in App. A.1. The two-loop contributions consist of genuine two-loop diagrams as well as one-loop diagrams with insertions of counterterm contributions from the vertices and the particle masses relevant at $\mathcal{O}(\alpha_{t}^{2})$, denoted by a cross in Fig. 14. The definition of the counterterms necessary for the calculation of these counterterm-inserted diagrams is discussed in Sec. 3.2.

The calculation of the loop integrals is performed in the framework of dimensional reduction [87, 88] in $D=4-2\varepsilon$ dimensions, where $\varepsilon$ is the dimensional regulator. For the MSSM, it has been proven that for two-loop corrections of $\mathcal{O}(\alpha_{t}^{2})$, SUSY is preserved by dimensional reduction and no additional SUSY-restoring counterterms are needed [89]. Since our $\mathcal{O}(\alpha_{t}^{2})$ calculations in the NMSSM are formally two-loop corrections calculated in the MSSM limit, this proof is directly applicable to our two-loop $\mathcal{O}(\alpha_{t}^{2})$ corrections in the NMSSM as well, and no SUSY-restoring counterterms are needed. For the calculation of the genuine two-loop integrals, we use the analytic results of the two-loop integrals, evaluated up to $\mathcal{O}(\varepsilon^{0})$, as presented in [90]. Some of the genuine two-loop integrals as well as all integrals appearing in counterterm-induced Feynman diagrams can be represented as the product of two one-loop one-point and two-point functions which are defined in [91, 92]. For these, we derived and implemented their expansion up to $\mathcal{O}(\varepsilon)$, given in App. B.

The calculation of the renormalized two-loop self-energies is performed fully analytically. The implementation of the loop integrals at $\mathcal{O}(\varepsilon^{0})$ for the two-loop and at $\mathcal{O}(\varepsilon)$ for the one-loop case allows for an explicit check of the cancellation of all UV-divergent poles of $\mathcal{O}(\varepsilon^{-2})$ and $\mathcal{O}(\varepsilon^{-1})$, which was confirmed explicitly at $\mathcal{O}(\alpha_{t}^{2})$ in the MSSM limit of the NMSSM, i.e. for vanishing $\lambda$ and $\kappa$ while keeping $\mu_{\text{eff}}$ fixed at the value given in the input file. Note that the calculation of the $\mathcal{O}(\alpha_{t}^{2})$ two-loop contributions to the neutral Higgs masses is always evaluated in the MSSM limit to ensure UV finiteness, while all other contributions are not restricted to the MSSM limit.

3.2 The Renormalization of the Higgs Sector

Since we work in the gaugeless limit at $\mathcal{O}(\alpha_{t}^{2})$, the Higgs potential does not depend on $s_{\theta_{W}}$ and $e$ anymore so that we restrict ourselves to one of the following two new sets of independent input parameters entering the Higgs potential at $\mathcal{O}(\alpha_{t}^{2})$,

[TABLE]

or

[TABLE]

depending on whether $M_{H^{\pm}}^{2}$ or $\text{Re}\,A_{\lambda}$ is chosen as independent input. These are the parameters that need to be renormalized in order to obtain a UV-finite result for the mass corrections. For the renormalization we replace the parameters by the renormalized ones and their corresponding counterterms as

[TABLE]

where the superscript $(n)$ stands for the $n$-loop level.

For the consistent incorporation of the $\mathcal{O}(\alpha_{t}^{2})$ corrections with the previously computed one-loop corrections in [56] and two-loop corrections of ${\cal O}(\alpha_{t}\alpha_{s})$ in [58], we use here also the mixed $\overline{\text{DR}}$-OS renormalization scheme, i.e. the parameters are renormalized as

[TABLE]

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

[TABLE]

for $\text{Re}\,A_{\lambda}$ as independent input. The counterterm matrix for the neutral Higgs mass matrix $\mathcal{M}_{hh}$ is obtained by inserting the replacements from Eqs. (3.58)-(3.67) in the tree-level mass matrix and expanding order by order, yielding

[TABLE]

In App. G, we give the explicit expressions of $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\mathcal{M}_{hh}$ and $\delta^{\text{\tiny(}\!\text{\tiny 2}\!\text{\tiny)}}\!\mathcal{M}_{hh}$ in terms of all parameter counterterms.

Also the Higgs fields need to be renormalized. We introduce the renormalization constants for the doublet and singlet fields before rotating to the mass eigenstates as

[TABLE]

where we defined

[TABLE]

In the following, we discuss in detail the wave function renormalization constants and parameter counterterms.

Higgs wave function renormalization constants

In analogy to our one-loop and two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ calculation [45, 56, 58] the Higgs fields are renormalized through $\overline{\text{DR}}$ conditions. The renormalization conditions ($n=1,2$, $i=h_{d},h_{u},h_{s}$)

[TABLE]

yield

[TABLE]

where the subscript ’div’ indicates that we take the divergent part only. The wave function renormalization constants are defined through the Higgs self-energies in the gauge basis. They are obtained from the derivatives with respect to $p^{2}$ of the diagrams shown in Fig. 1 at one-loop and in Fig. 14 at two-loop level by choosing the mixing matrix of the neutral Higgs bosons to be diagonal. At ${\cal O}(\alpha_{t})$ we find

[TABLE]

and at $\mathcal{O}(\alpha_{t}^{2})$ we have, for the $\overline{\text{DR}}$ renormalization scheme of the top/stop sector,

[TABLE]

This result is in agreement with [93, 94] derived on the basis of the renormalization group. In the OS scheme of the top/stop sector, we find

[TABLE]

Note that the superscripts $\overline{\text{DR}}$ and OS on $\delta^{\text{\tiny(}\!\text{\tiny 2}\!\text{\tiny)}}\!Z_{H_{u}}$ refer to the renormalization of the top/stop sector. The expression for $dm_{t}^{\alpha_{t}}$ is given in Eq. (C.183) and for $dv^{\alpha_{t}}$ reads

[TABLE]

where

[TABLE]

The tadpole counterterms

The renormalization conditions for the tadpoles are chosen such that the minimum of the Higgs potential does not change at two-loop order, leading to ($\phi=h_{d},h_{u},h_{s},a_{d},a_{s}$)

[TABLE]

*where $t_{\phi}^{\text{\tiny{ (1)}}}$ and $t_{\phi}^{\text{\tiny{ (2)}}}$ represent the one- and two-loop tadpole contributions shown in Fig. 2 and Fig. 3 for $\mathcal{O}(\alpha_{t})$ and $\mathcal{O}(\alpha_{t}^{2})$, respectively.

The charged Higgs boson mass counterterm

If the charged Higgs boson mass is chosen as independent input, hence defined as OS parameter, we renormalize it in the OS scheme accordingly. In the approximation of vanishing external momentum, the OS counterterm of the charged Higgs mass at one-loop order is given by

[TABLE]

and the counterterm at two-loop level reads

[TABLE]

with

[TABLE]

where

[TABLE]

and $\Delta^{(2)}Z_{i}$ ($i=H_{u},H_{d},S$) has been defined in Eq. (3.74). The $\mathcal{O}(\alpha_{t})$ contributions to the unrenormalized self-energy of the charged Higgs boson are depicted in Fig. 4. The sum of all diagrams yields the following analytic expression of the unrenormalized one-loop self-energy of the charged Higgs boson at $\mathcal{O}(\alpha_{t})$,

[TABLE]

*with the one-loop scalar integrals $A_{0}$ and $B_{0}$ defined in App. B. Note that since the mass of the charged Higgs boson is calculated at vanishing external momentum the counterterm of the mass of the charged Higgs boson also involves the field renormalization constants. The two-loop $\mathcal{O}(\alpha_{t}^{2})$ contributions to the unrenormalized charged Higgs self-energies are depicted in Fig. 15 in App. A.2. Due to its lengthy structure, we do not display the analytic result of the corresponding unrenormalized two-loop self-energy explicitly here.

If $\text{Re}\,A_{\lambda}$ is chosen as independent input parameter instead of $M_{H^{\pm}}^{2}$, then the counterterms $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!M_{H^{\pm}}^{2}$ and $\delta^{\text{\tiny(}\!\text{\tiny 2}\!\text{\tiny)}}\!M_{H^{\pm}}^{2}$ of the charged Higgs boson mass at one- and two-loop order, respectively, are calculated as functions of all other counterterms by inserting the two-loop expansions of Eqs. (3.58) to (3.67) in the formula for the charged Higgs boson mass, Eq. (2.23). The explicit formulae of the counterterms are presented in App. H. The loop-corrected charged Higgs mass $(M_{H^{\pm}}^{2})^{(2)}$ is calculated iteratively by solving

[TABLE]

with the renormalized one-loop self-energy of the charged Higgs boson

[TABLE]

evaluated at the scale of the squared tree-level charged Higgs mass and the renormalized two-loop self-energies at $\mathcal{O}(\alpha_{t}\alpha_{s})$ and $\mathcal{O}(\alpha_{t}^{2})$, respectively,

[TABLE]

*evaluated in the approximation of zero external momentum.

The VEV counterterm

The OS renormalized one-loop counterterm $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!v/v$ of the VEV is given by

[TABLE]

with the gauge bosons renormalized OS. Their OS counterterms in the gaugeless approximation read

[TABLE]

Here, $\Sigma^{(1),T}_{V}(0)$ ($V=W,Z$*) is the transverse part of the unrenormalized one-loop vector boson self-energy evaluated at vanishing external momentum. The relevant diagrams at $\mathcal{O}(\alpha_{t})$ are depicted in Figs. 5 and 6. *

Note that while $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!M_{V}^{2}$ and $M_{V}^{2}$ are separately zero in the gaugeless limit, their ratio is non-zero for each gauge boson $V$ and hence contributes to our $\mathcal{O}(\alpha_{t}^{2})$ calculation. The explicit evaluation of the UV divergent part of $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!v/v$ shows that it is related to $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!Z_{H_{u}}$ as

[TABLE]

*as expected from Refs. [93, 94].

The $\tan\beta$ counterterm

The ratio of the vacuum expectation values of the Higgs doublets, $\tan\!\beta$, is renormalized in the $\overline{\text{DR}}$ scheme with the counterterm given by **[95, 96, 97, 98, 99, 100]**

[TABLE]

*where the subscript $(n)$ again indicates the loop level. Note that the last identity only holds at $\mathcal{O}(\alpha_{t})$ and $\mathcal{O}(\alpha_{t}^{2})$ in our approximation, but is not valid in general. *

The remaining $\overline{\mbox{DR}}$ counterterms

The counterterms of the remaining $\overline{\text{DR}}$ parameters $|\lambda|,|\kappa|,v_{s},\mbox{Re}A_{\kappa},\varphi_{\lambda},\varphi_{\kappa},\varphi_{u},$ and $\varphi_{s}$ have to cancel the left-over UV-divergent parts of the self-energies of the neutral Higgs bosons. For the $\mathcal{O}(\alpha_{t})$ one-loop counterterm of $|\lambda|$, we find

[TABLE]

while the two-loop counterterm $\delta^{\text{\tiny(}\!\text{\tiny 2}\!\text{\tiny)}}\!|\lambda|$ as well as the remaining one- and two-loop counterterms of $|\kappa|,v_{s},\mbox{Re}A_{\kappa},$ $\varphi_{u},\varphi_{s},\varphi_{\lambda}$ and $\varphi_{\kappa}$ are not needed to yield a UV-finite result of the neutral Higgs masses at $\mathcal{O}(\alpha_{t}^{2})$.

3.3 Renormalization of the Quark/Squark

Sector

The set of independent parameters to be renormalized in the third generation quark/squark sector is given in Eq. (2.30). The renormalization conditions of the corresponding counterterms

[TABLE]

*are specified in the following. Since all counterterms of the quark/squark sector necessary for the calculation of the $\mathcal{O}(\alpha_{t}^{2})$ corrections to the neutral Higgs boson masses are of one-loop order, we suppress the superscript $(1)$ of the counterterms in this subsection. *

We proceed along the same lines as in our ${\cal O}(\alpha_{t}\alpha_{s})$ calculation **[58]** where we have implemented in NMSSMCALC both the OS and $\overline{\text{DR}}$ renormalization scheme for the (s)quark sector. We follow the SLHA **[84, 85]** where the input top quark mass is understood to be the pole mass whereas the soft SUSY breaking masses and trilinear couplings are $\overline{\mbox{DR}}$ parameters at the renormalization scale $\mu_{R}=M_{\text{SUSY}}$. The translation between the two schemes has to be done consistently both in the counterterm part and at the level of the input parameters. Expanding the OS and $\overline{\mbox{DR}}$ counterterms of the parameters $X=m_{t},m_{\tilde{Q}_{3}},m_{\tilde{t}_{R}},A_{t}$ in terms of the dimensional regulator $\varepsilon$, we have

[TABLE]

Note that in our definition of the parameters of the OS scheme we did not take into account any terms proportional to $\varepsilon$, i.e. $\varepsilon\delta X_{\varepsilon}$. These terms, that could also be chosen to be included, cf. **[101, 102*]**, would manifest themselves as additional finite contributions originating from the counterterm inserted diagrams multiplying $1/\varepsilon$ terms from the one-loop functions with the $\varepsilon$ parts of the counterterms. We verified through explicit calculation that the contributions of finite terms arising through the inclusion of the $\mathcal{O}(\varepsilon)$ terms of the OS-defined counterterms cancel in the calculation of the $\mathcal{O}(\alpha_{t}^{2})$ corrections. Therefore, we neglect the $\mathcal{O}(\varepsilon)$ contributions from the counterterms and apply our thus defined OS scheme consistently throughout the whole calculation. *

In case of $\overline{\text{DR}}$ renormalization of the (s)quark sector the $\overline{\text{DR}}$ top quark mass $m_{t}^{\tiny{\overline{\mbox{DR}}}}$ has to be computed from the corresponding top pole mass, as described in App. C. If the OS scheme is chosen for the (s)quark sector, then the translation of the parameters $m_{\tilde{Q}_{3}},m_{\tilde{t}_{R}}$ and $A_{t}$ from the $\overline{\mbox{DR}}$ scheme to the OS scheme is performed by

[TABLE]

In the above equations, the finite counterterm parts have to be computed with OS input parameters, which we achieve by applying an iterative procedure. Note that we include both ${\cal O}(\alpha_{s})$ and ${\cal O}(\alpha_{t})$ corrections into these conversions in the numerical analysis. The OS conditions in the NMSSM (s)quark sector are the same as the ones in the complex MSSM presented in **[103, 104]**. We give here, for completeness, the expressions of the counterterms:

• •

Decomposing the unrenormalized top quark self-energy $\Sigma_{t}$ as

[TABLE]

in terms of the left- and right-handed projectors $P_{L/R}=(1\mp\gamma_{5})/2$, the top mass counterterm reads

[TABLE]

Here $\widetilde{\mbox{Re}}$ indicates that the real part is taken only of the one-loop function, but not of the parameters. The $\mathcal{O}(\alpha_{t})$ contributions to the top mass counterterm are depicted in Fig. 7.

• •

The counterterm of the trilinear stop coupling parameter reads

[TABLE]

where

[TABLE]

with $\Sigma_{\tilde{t}_{i}\tilde{t}_{j}}$ $(i,j=1,2)$ denoting the unrenormalized self-energy for the transition $\tilde{t}_{i}\to\tilde{t}_{j}$, whose contributions at $\mathcal{O}(\alpha_{t})$ are depicted in Fig. 8.

• •

The counterterms for the soft SUSY breaking left-handed squark and right-handed stop mass parameters are given by

[TABLE]

3.4 Tools and Checks

We performed two independent calculations of the two-loop NMSSM Higgs boson masses at $\mathcal{O}(\alpha_{t}^{2})$ and cross-checked the results against each other. In both cases, we used SARAH 4.14.0 **[105, 69, 70, 71]** to generate a model file for the complex NMSSM that was used in FeynArts 3.8 **[106, 107]** to generate all one- and two-loop Feynman diagrams necessary for the calculation of the mass corrections. For the simplification of expressions containing the Dirac matrices as well as for the calculation of fermion traces, we used FeynCalc 8.2.0 **[108, 109]**. The tensor reduction of the one- and two-loop integrals was also performed in FeynCalc, where for the two-loop case, we additionally used the TARCER **[110*]** extension of FeynCalc. *

As additional cross-checks for the calculation of the fermion traces, we used the Mathematica package FormTracer **[111]**. The fact that FormTracer treats the Dirac matrix $\gamma_{5}$ in the Larin scheme **[112]**, while FeynCalc 8.2.0 treats $\gamma_{5}$ per default in the ’naive scheme’ **[113]**, allowed us to compute the fermion traces in the framework of the two different $\gamma_{5}$ schemes. For the traces relevant for our $\mathcal{O}(\alpha_{t}^{2})$ computations, we explicitly verified that a change of the $\gamma_{5}$ scheme does not affect the calculation of the NMSSM Higgs boson masses at $\mathcal{O}(\alpha_{t}^{2})$.

We have compared the results of the self-energies at $\mathcal{O}(\alpha_{t}^{2})$ with the corresponding self-energies in FeynHiggs **[114, 115, 116, 117, 118, 119, 120, 121, 122]**. Since the $\mathcal{O}(m_{t}^{2}\alpha_{t}^{2})$ corrections are equivalent in the MSSM and the NMSSM with $\mu$ in the MSSM identified as $\mu^{\text{eff}}$ in the NMSSM and with further contributions of $\mathcal{O}(\alpha_{t}^{2})$ neglected, we found agreement after adapting the counterterms for the weak mixing angle and the vacuum expectation value to the renormalization scheme applied in FeynHiggs, see for example **[82]**, thus ensuring the same input values. In addition, we also compared to the on-shell and $\overline{\text{DR}}$ results for real parameters and the mass of the pseudoscalar Higgs boson $M_{A}$ as input of Ref. **[123]**, using the corresponding computer code, and found also complete agreement.

4 Numerical Analysis

4.1 The Parameter Scan

For our numerical analysis we perform a scan in the NMSSM parameter space in order to find scenarios that are compatible with the recent experimental constraints. We proceed as described in **[124, 125, 126]**, where also further details can be found. We vary $\tan\beta$, $\lambda$ and $\kappa$ in the ranges

[TABLE]

thus not exceeding the rough constraint

[TABLE]

which ensures perturbativity. Further parameter scan ranges are summarized in Tab. 1. We set

[TABLE]

and the mass parameters of the first and second generation sfermions are chosen as

[TABLE]

The soft SUSY breaking trilinear couplings of the first two generations are set equal to the corresponding values of the third generation. Following the SLHA format **[84, 85]**, the soft SUSY breaking masses and trilinear couplings are understood as $\overline{\text{DR}}$ parameters at the scale

[TABLE]

The SM input parameters have been chosen as **[127, 128]**

[TABLE]

The spectrum of the Higgs particles including the higher-order corrections presented in this work is calculated with the new NMSSMCALC version which includes the higher-order corrections of ${\cal O}(\alpha_{t}^{2})$ calculated in this paper. For the scan, $M_{H^{\pm}}$ has been used as input parameter, cf. Eq. (2.25).151515Note, however, that in NMSSMCALC we also have the option to set $A_{\lambda}$ as input parameter which is then renormalized in the $\overline{\mbox{DR}}$ scheme. NMSSMCALC also checks for the constraints from the electric dipole moments (EDMs) that become relevant for the CP-violating case **[129]**. 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

[TABLE]

Agreement with the Higgs exclusion limits from LEP, Tevatron and LHC is checked by using HiggsBounds 5.3.2 **[130, 131, 132]**, and with the Higgs rates by using HiggsSignals 2.2.3 **[133]**. We demand the total $\chi^{2}$ computed by HiggsSignals with our given effective coupling factors to be compatible with the total $\chi^{2}$ of the SM within 1$\sigma$*. The required input for HiggSignals is computed with NMSSMCALC. *

We also take into account the most relevant LHC exclusion bounds on the SUSY masses. These constrain the gluino mass and the lightest squark mass of the second generation to lie above 1.8 TeV, see **[134]**. The stop and sbottom masses in general have to be above 800 GeV **[135, 134]**, and the slepton masses above 400 GeV **[134]**.

We perform the scan in the limit of the CP-conserving NMSSM. For the numerical analysis we start from a valid parameter point and subsequently turn on various CP-violating phases. The thus obtained parameter points do not necessarily fulfill the constraints from the EDMs any more but nevertheless, we keep them for illustrative purposes. The strongest constraint originates from the electron EDM **[136]**. We check the EDMs of our parameter points against the experimental limit given by the ACME collaboration **[137]**.

4.2 Results

For our numerical analysis we have chosen two sample points among the parameter points compatible with all described constraints that we obtained from our scan. They both feature a SM-like Higgs boson with mass around 125 GeV at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ when - in one case - the top/stop sector is renormalized OS and - in the other case - the top/stop sector is $\overline{\mbox{DR}}$ renormalized. We call the former point P1OS, the latter P2DR. In the following we give the relevant input values for these two points. Note that we deliberately chose parameter points with not too large NMSSM-specific coupling values $\lambda$ and $\kappa$ as at two-loop order we so far do not include the corrections proportional to these couplings. This is left for future work. Since we include the complete set of one-loop corrections, however, our results for these parameter points should not be affected significantly by the missing corrections. This can also be inferred from the results given in **[81]**.

Parameter Point P1OS:

Besides the SM values defined above, the parameter point is given by the following soft SUSY breaking masses and trilinear couplings,

[TABLE]

with the CP phases given by

[TABLE]

The remaining input parameters have been set to161616The imaginary part of $A_{\kappa}$ is obtained from the tadpole condition.**

[TABLE]

As required by the SLHA, $\mu_{\text{eff}}$ is taken as input parameter, from which $v_{s}$ and $\varphi_{s}$ are obtained through Eq. (2.24). The parameters $\lambda,\kappa,A_{\kappa},\mu_{\text{eff}},\tan\beta$ as well as the soft SUSY breaking masses and trilinear couplings are understood as $\overline{\text{DR}}$ parameters at the scale $\mu_{R}=M_{\text{SUSY}}$171717For $\tan\beta$ this is only the case if it is read in from the block EXTPAR as done in NMSSMCALC. Otherwise it is the $\overline{\text{DR}}$ parameter at the scale $M_{Z}$.. The charged Higgs mass, however, is an OS parameter. The SUSY scale $M_{\text{SUSY}}$ is set to be

[TABLE]

In the following the subscript ’eff’ for $\mu$ is dropped and we use the expressions OS and $\overline{\text{DR}}$ in order to refer to the chosen renormalization conditions in the top/stop sector only, while all the other renormalization conditions remain unchanged.

Parameter Point P2DR:

Besides the SM values defined above, the parameter point is given by the following soft SUSY breaking masses and trilinear couplings,

[TABLE]

The remaining input parameters have been set to

[TABLE]

4.3 Results and Analysis Parameter Point P1OS

In Tab. 2 we summarize the values of the masses that we obtain for the chosen P1OS at tree level, at one-loop level, at two-loop level including only the ${\cal O}(\alpha_{t}\alpha_{s})$ corrections and at two-loop level including furthermore our newly calculated ${\cal O}(\alpha_{t}^{2})$ corrections. In Tab. 3 the results are given for the $\overline{\text{DR}}$ scheme in the top/stop sector. The tables also contain the information on the main singlet/doublet and scalar/pseudoscalar component of the respective mass eigenstate. The tree-level stop masses obtained within the OS and $\overline{\text{DR}}$ scheme are given by

[TABLE]

The $\overline{\text{DR}}$ top quark mass in our scenario amounts to $m_{t}^{\overline{\text{DR}}}=141.8$ GeV, and has been computed as described in App. C.

*The scenario features three heavy Higgs bosons with masses between 700 and 900 GeV. Additionally, we have a light Higgs boson with mass value below 125 GeV. For the correct interpretation of the importance of the loop corrections, the Higgs bosons with a similar admixture must be compared and not the ones according to their mass ordering. In the following plots we therefore label the Higgs bosons according to their dominant admixture and not by their mass ordering. However, in the scenario P1OS both orderings lead to the same result. The admixtures of the Higgs bosons are detailed in the Tables. While the dominance is very pronounced at loop level, at tree level, for $H_{1}$ and $H_{2}$ the $h_{s}$ and $h_{u}$ admixtures are approximately the same, respectively. For $H_{1}$, however, the $h_{s}$ admixture is larger, and for $H_{2}$ the $h_{u}$ admixture is larger. In order to comply with the Higgs rate measurements of the LHC the SM-like Higgs boson has to be $h_{u}$ dominated, as is the case for the 125 GeV Higgs boson in our scenario (at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ with OS renormalization in the top/stop sector). *

*Defining the absolute value of the relative change in the mass value when going from loop order $a$ to loop order $b$ including the next level of corrections, as $|m^{b}-m^{a}|/m^{a}$, we see that the lightest singlet-like Higgs boson $H_{1}$ receives rather large one-loop corrections of ${\cal O}(16\%)$. The two-loop corrections are below the per-cent level. The $h_{u}$-dominated Higgs boson $H_{2}$ receives important one-loop corrections of ${\cal O}(48\%)$ in the OS scheme, respectively ${\cal O}(23\%)$ in the $\overline{\text{DR}}$ scheme. The two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ corrections reduce the mass value by $12\%$ in the OS scheme and add 5% in the $\overline{\text{DR}}$ scheme so that the mass values in the two renormalization schemes move close to each other with values of $\sim 118$ GeV. The two-loop ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ correction adds another 6% to the ${\cal O}(\alpha_{t}\alpha_{s})$ result in the OS scheme while in the $\overline{\text{DR}}$ scheme this loop correction barely alters the mass value, so that at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ the mass values in the two renormalization schemes move further away from each other. *

In Fig. 9, we display the one-loop and two-loop corrected mass values $M_{h_{u}}$ of the $h_{u}$-like Higgs boson including at two-loop level the ${\cal O}(\alpha_{t}\alpha_{s})$ and the ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections, as a function of the $\overline{\text{DR}}$ parameters $A_{t}$ (left) and $\varphi_{A_{t}}$ (right) for both OS (full lines) and $\overline{\text{DR}}$ (dashed lines) renormalization in the top/stop sector. Starting from our initial parameter point P1OS, we vary $A_{t}$ in the left plot (keeping $\varphi_{A_{t}}$ unchanged) and for the right plot we vary $\varphi_{A_{t}}$ while $A_{t}$ is fixed at its initial absolute value. The point P1OS corresponds to the values at $A_{t}^{\overline{\text{DR}}}=-1922$ GeV (left) and vanishing CP-violating phase $\sin\varphi_{A_{t}}=0$, which corresponds to $\varphi_{A_{t}}=\pm\pi$ in the right plot. The two lower plots show the relative difference in the masses (at the same loop order, which is indicated by the colour of the lines) obtained using the two renormalization schemes in the top/stop sector,

[TABLE]

*as a function of $A_{t}$ and $\varphi_{A_{t}}$, respectively. We show the corrections for the $h_{u}$-like Higgs boson mass, as it is affected most by the ${\cal O}(\alpha_{t}\alpha_{s})$ and ${\cal O}(\alpha_{t}^{2})$ corrections. The plots confirm what the discussion of P1OS already revealed. For the $\overline{\text{DR}}$ renormalization scheme the two-loop corrections of ${\cal O}(\alpha_{t}\alpha_{s})$ do not change the one-loop result as much as in the OS scheme. This behaviour can be understood as follows. The $\overline{\text{DR}}$ renormalization in the top/stop sector requires the conversion of the input OS top-quark mass to the $\overline{\text{DR}}$ mass at the scale $\mu_{R}=M_{\text{SUSY}}$. This conversion is described in App. C. It includes the ${\cal O}(\alpha_{s}+\alpha_{t}+\alpha_{s}^{2})$ contributions in the conversion of the SM OS top-quark mass to the $\overline{\mbox{MS}}$ mass at $\mu_{R}=M_{Z}$ and the ${\cal O}(\alpha_{s}+\alpha_{t})+{\cal O}((\alpha_{s}+\alpha_{t})^{2})$ corrections in the renormalization group equations needed for the running from $\mu_{R}=M_{Z}$ to $\mu_{R}=M_{\text{SUSY}}$. In this way the one-loop mass in the $\overline{\text{DR}}$ renormalization scheme already includes higher-order corrections beyond the one-loop level. Furthermore, the ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections barely change the ${\cal O}(\alpha_{t}\alpha_{s})$ mass value for the $\overline{\text{DR}}$ scheme, whereas the OS renormalization leads to a further change of a few GeV, so that the relative difference in the mass values at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ with values of 5-6% is larger compared to the relative difference at ${\cal O}(\alpha_{t}\alpha_{s})$ with values below 4-5% down to almost 0%. As can be read off from the plots, the effect on the Higgs mass values due to the change of the renormalization scheme becomes more pronounced when moving from ${\cal O}(\alpha_{t}\alpha_{s})$ to ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$, increasing thus the estimate of the remaining theoretical uncertainty due to missing higher-order corrections from the change of the renormalization scheme. At first sight, this might be counter-intuitive. Both types of corrections are of two-loop order, however. A reduction in the theoretical uncertainty cannot necessarily be expected when further contributions at the same loop level are taken into account. With the loop order included in our conversion of the parameters, we estimate an uncertainty due to missing corrections of three-loop order ${\cal O}(\alpha_{t}^{2}\alpha_{s}+\alpha_{t}\alpha_{s}^{2})$ and of four-loop order ${\cal O}(\alpha_{t}^{3}\alpha_{s}+\alpha_{t}^{2}\alpha_{s}^{2}+\alpha_{s}^{3}\alpha_{t})$ contributions when we include the ${\cal O}(\alpha_{t}\alpha_{s})$ corrections. The inclusion of the ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections provides an estimate of the missing contributions as before but additionally also of terms of the order ${\cal O}(\alpha_{t}^{3})$ and ${\cal O}(\alpha_{t}^{4})$. This shows that care has to be taken in the estimate of the uncertainty due to missing higher-order corrections at a given loop-level when only parts of the loop contributions at this given order are included. An estimate based on the ${\cal O}(\alpha_{t}\alpha_{s})$ corrections alone would be more optimistic than after the inclusion of additionally the ${\cal O}(\alpha_{t}^{2})$ corrections. *

*Note also that the one-loop corrections in the OS scheme are almost symmetric with respect to the sign change of $A_{t}$ in contrast to the $\overline{\text{DR}}$ scheme. This behaviour results from the threshold effect in the conversion of the top OS to the $\overline{\text{DR}}$ mass, which depends on the sign of $A_{t}$. In contrast, the sign dependence in the conversion of the $\overline{\text{DR}}$ stop parameters to OS parameters for the OS scheme almost cancels out. The dependence of the loop-corrected mass values with varying $\varphi_{A_{t}}$ is on the one hand due to the genuine dependence on the phase and on the other hand due to the dependence of the stop mass values on the phase. The stronger dependence of the one-loop $\overline{\text{DR}}$ mass values compared to the higher-order and OS values again is due to the necessary conversion from the OS to the $\overline{\text{DR}}$ top mass value in this renormalization scheme. Finally, we want to remark that the non-zero $\sin\varphi_{A_{t}}$ may lead to scenarios that are not compatible with the EDMs any more.181818Actually, a check with NMSSMCALC showed that the EDM constraints are not fulfilled anymore for a non-zero phase |\varphi_{A_{t}}|\raisebox{-3.69899pt}{~{}\shortstack{>$\\ [-0.07cm]$\sim}}~{}0.08\pi (no additional CP violation from other phases). We still keep these scenarios in the plot for illustrative purposes. *

Figure 10 shows the one-loop and two-loop corrections of ${\cal O}(\alpha_{t}\alpha_{s})$ and ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$, respectively, to the $h_{u}$-like Higgs boson mass as a function of the phases $\varphi_{A_{t}}$ and $\varphi_{\mu}$ where the latter denotes the phase of $\mu_{\text{eff}}$. The left plot shows results obtained for the $\overline{\text{DR}}$ renormalization scheme in the top/stop sector, the right plot those for the OS scheme. Starting from the above defined parameter point, we turn on separately one of the two phases. For illustrative purposes we vary the phases also beyond values already excluded by experiment.191919The EDM constraints are not fulfilled any more for a non-zero phase |\varphi_{\mu}|\raisebox{-3.69899pt}{~{}\shortstack{>$\\ [-0.07cm]$\sim}}~{}(9.5\cdot 10^{-10})\pi (no additional CP violation from other phases). For the plots, we have varied $\varphi_{\mu}$ in such a way that the CP-violating phase, which appears already at tree level in the Higgs sector, i.e. $\varphi_{y}=\varphi_{\kappa}-\varphi_{\lambda}+2\varphi_{s}-\varphi_{u}$, remains zero. For this, the phases $\varphi_{\lambda}$ and $\varphi_{s}$ were varied at the same time, in particular we set $\varphi_{\lambda}=2\varphi_{s}=2/3\varphi_{\mu}$. The phases $\varphi_{\kappa}$ and $\varphi_{u}$ were kept zero ($|A_{\kappa}|$* is kept constant). We thus make sure that the main effect in the plots on the dependence on the phases originates from the loop corrections. *

*As can be inferred from the plots, for both schemes the shape of the variation of the two-loop corrections with the phases is very similar. The dependence on the phase $\varphi_{A_{t}}$ is stronger than the one on $\varphi_{\mu}$. Overall, the influence of the investigated complex phases on the loop corrections is quite small as we study purely radiatively induced CP violation here. The effect of the phases on the mass values is at most 3% at one-loop and below the percent level at two-loop order, and much smaller than the overall mass corrections at each loop level. *

The conversion of the OS top-quark mass to the $\overline{\text{DR}}$ mass in the $\overline{\text{DR}}$ renormalization scheme resums higher-order corrections into the fixed-order calculation, explaining the large difference between the one-loop $\overline{\text{DR}}$ and OS results, and the better convergence when going from the one- to the two-loop level in the $\overline{\text{DR}}$ scheme. Thus Fig. 10 shows that the absolute relative corrections between the two loop levels, defined, for a fixed renormalization scheme, as

[TABLE]

*with the superscripts {2, 1} referring to the two- and one-loop level, respectively, and the superscript $x$ to the two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ and ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections, are larger in the OS than in the $\overline{\text{DR}}$ renormalization scheme. For the ${\cal O}(\alpha_{t}\alpha_{s})$ corrections they amount to about 3-5% in the latter case, whereas the masses are reduced by 11-12% in the former case. The ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections make up for a relative correction between 3.7 and 5.5% in the $\overline{\text{DR}}$ scheme and between -6 and -7% in the OS scheme with respect to the one-loop case. Overall, the two-loop corrections reduce the one-loop masses in the OS scheme whereas they are positive in the $\overline{\text{DR}}$ scheme. As already discussed and commented on above, including the ${\cal O}(\alpha_{t}^{2})$ corrections worsens the convergence of the higher-order corrections. *

An estimate of the theoretical uncertainty due to the missing higher-order corrections can also be obtained from a change of the renormalization scale $\mu$. In Fig. 11 we depict the one- and two-loop corrected mass values of the $h_{u}$-like Higgs boson, again including at two-loop level the ${\cal O}(\alpha_{t}\alpha_{s})$ and the ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections, as a function of the renormalization scale in terms of the default renormalization scale $\mu_{R}\equiv\mu_{0}=M_{\text{SUSY}}$. The renormalization scale is varied between 1/2 and 2 times the value of the central scale $\mu_{0}$. The lower panel shows the absolute value of the relative change of the mass at the scale $\mu$ with respect to the value obtained for $\mu_{0}$,

[TABLE]

For this the loop order and type of conversion and RGE running is kept fix, as indicated by the type and color of each line. Note, that the consistent comparison of the results for different renormalization scales also requires the conversion of the input parameters to the new scale $\mu_{R}=\mu$ from the original set given at $\mu_{R}=\mu_{0}$. We perform this conversion for the top quark mass by applying the renormalization group equations and the procedure described in App.  C, including the ${\cal O}(\alpha_{s}+\alpha_{t}+\alpha_{s}^{2})$ contributions in the conversion of the SM OS top-quark mass to the $\overline{\mbox{MS}}$ mass at $\mu_{R}=M_{Z}$ and the ${\cal O}(\alpha_{s}+\alpha_{t}+(\alpha_{s}+\alpha_{t})^{2})$ corrections in the RGEs for the running values of $m_{t}$ and $\alpha_{s}$ from $\mu_{R}=M_{Z}$ to $\mu_{R}=\mu$. The corresponding results are denoted by ’EWrun’ in the plot. The relative change of the corrections obviously depends on the running of the parameters that is applied. The plot also includes the results ’QCDrun’, where we include only the ${\cal O}(\alpha_{s}+\alpha_{s}^{2})$ contributions of the conversion of the top mass at the scale $\mu_{R}=m_{t}$ and into RGEs, as done in Ref. **[58*]** where we computed the ${\cal O}(\alpha_{t}\alpha_{s})$ corrections. The inclusion of the EW contributions in the running is consistent with the ${\cal O}(\alpha_{t}^{2})$ contributions in our two-loop fixed order masses, leading (at all loop levels) to a better convergence of the higher-order results compared to those where only QCD running is included. The remaining NMSSM $\overline{\text{DR}}$ input parameters are evaluated from the SLHA default input scale $M_{\text{SUSY}}$ to the renormalization scale $\mu$ by applying the RGEs given in App. E and by using as input into the RGEs $\alpha_{s}$ and $y_{t}$ at the scale $\mu=M_{\text{SUSY}}$ as obtained by applying the procedure ’EWrun’ or ’QCDrun’, respectively, in the curve with the corresponding name. *

*For both schemes we observe that the inclusion of the two-loop corrections of ${\cal O}(\alpha_{t}^{2})$ worsens the scale dependence as compared to the two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ corrections alone. In the OS scheme this even leads to a larger scale dependence than at one-loop order.202020The one-loop scale dependence is, however, reduced in the OS scheme with respect to the $\overline{\text{DR}}$ scheme due to the larger number of running parameters in the latter scheme. In the $\overline{\text{DR}}$ scheme we encounter a scale dependence that is not present in the OS scheme through the conversion of the OS to the $\overline{\text{DR}}$ top-quark mass at one-loop order, rendering the dependence on the scale $\mu$ larger than in the OS scheme so that in the $\overline{\text{DR}}$ scheme the black one-loop curve lies above the curves showing the two-loop scale dependences. In the OS scheme, the running parameters are the $\overline{\text{DR}}$ input parameters whereas the top/stop parameters are treated OS in contrast to the $\overline{\text{DR}}$ scheme.212121To be precise, in the OS scheme the values of the running parameters are calculated in the same way as before for the $\overline{\text{DR}}$ scheme. The soft SUSY breaking parameters $A_{t}$, $m_{\tilde{t}_{R}}$ and $m_{\tilde{Q}_{3}}$, however, are evolved to the scale $\mu$ and then converted to the corresponding OS parameters. For the top quark mass the input OS value is used. Additionally, in both schemes $\alpha_{s}$ is running. The latter enters the one-loop result only via the running of the other $\overline{\text{DR}}$ parameters and directly only at two-loop level. At ${\cal O}(\alpha_{t}\alpha_{s})$ a cancellation between the running parameters and the fixed order results leads to the observed flat behaviour in the scale dependence. The inclusion of additional ${\cal O}(\alpha_{t}^{2})$ terms, however, leads to a larger scale dependence again. These observations show that care has to be taken with respect to the estimate of the uncertainty due to missing higher-order corrections based on a fixed order calculation that takes into account only partial corrections at the given loop order and the conclusions drawn from the scale dependence. Overall, as expected and discussed above, however, the two-loop corrections are smaller than the one-loop corrections and the renormalization scale dependence is reduced. *

*In order to further illustrate the interplay between the RGEs applied on the running parameters and the fixed order results we show in Fig. 12 the two-loop ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections to the SM-like Higgs boson mass with $\overline{\mbox{DR}}$ renormalization in the top/stop sector where the running $\overline{\mbox{DR}}$ top mass value is obtained as described in App. C and denoted by ’EWrun’ in the previous plot (red line ’w/o $g_{1},g_{2}$’). The black line (’full $g_{1},g_{2}$’) is obtained by including the $g_{1}$ and $g_{2}$ contributions in the conversion and in the running, as described in App. D. For the blue line (’partly $g_{1},g_{2}$’), we included the $g_{1}$ and $g_{2}$ contributions in the RGEs of $m_{t}$, $\alpha_{s}$ and $\lambda$ but not in the RGE for the VEV, and we did not include $g_{1}$ nor $g_{2}$ contributions in the conversion at the scale $M_{Z}$. As can be inferred from the plot the results without the $g_{1}$, $g_{2}$ contributions - that have also not been included in our fixed-order two-loop results - and those with their full inclusion are rather close to each other, differing by less than 1 GeV. The partial inclusion, however, drives the result away from the two previous ones by up to about 3 GeV. This is because here we do not include the running of the VEV and therefore neglect the large positive gauge contributions to the RGE of the VEV. This once again shows that care has to be taken in the conversion of the running parameters included in the higher-order corrections and the estimate of the uncertainty due to the missing higher-order contributions from the scale dependence. *

4.4 Results Parameter Point P2DR

For completeness, we give also the results for a parameter point that has been chosen such that the ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections, now with $\overline{\text{DR}}$ renormalization applied in the top/stop sector, yield a SM-like Higgs boson with mass around 125 GeV. We called this point P2DR and presented its input values above. Table 4 summarizes the mass values that we obtain for the $\overline{\text{DR}}$ scheme in the top/stop sector at tree level, at one-loop level and at two-loop level including only the ${\cal O}(\alpha_{t}\alpha_{s})$ and the ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ corrections, respectively. In Tab. 5 the results are given for the OS scheme in the top/stop sector. The tables also contain the information on the main singlet/doublet and scalar/pseudoscalar component of the respective mass eigenstate. The tree-level stop masses obtained within the $\overline{\text{DR}}$ and OS scheme are given by

[TABLE]

*The $\overline{\text{DR}}$ top mass in our scenario amounts to $m_{t}^{\overline{\text{DR}}}=146.64$ GeV. *

In this scenario, the lightest Higgs boson corresponds to the SM-like one with a mass around 125 GeV at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ for $\overline{\text{DR}}$ renormalization in the top/stop sector. The remaining masses are around 400 GeV for the singlet-like pseudoscalar, the masses of the CP-even and CP-odd MSSM-like Higgs bosons are around 790 GeV, and the heaviest Higgs boson has a mass of 1828 GeV and is mostly CP-even singlet-like. The one-loop corrections to the mass of the $h_{u}$-dominated SM-like Higgs boson $H_{1}$ are important and lead to a relative increase of the mass value by 35% in the $\overline{\text{DR}}$ scheme, and by 59% in the OS scheme. The two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ corrections add a relative correction of 3% in the $\overline{\text{DR}}$ scheme and reduce the OS one-loop mass by 13% in the OS scheme so that the mass values in both renormalization schemes are close at this loop order. The ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ lead to a relative increase of the mass value by 1% in the $\overline{\text{DR}}$ scheme, and the OS renormalization in the top/stop sector increases the two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ mass by about 8% when the ${\cal O}(\alpha_{t}^{2})$ corrections are included as well. Overall, the relative corrections are somewhat larger for P2DR than for P1OS. Altogether, however, we observe the same behaviour, namely a slightly worse convergence of the loop corrections after inclusion of the ${\cal O}(\alpha_{t}^{2})$ corrections. Since this is also a two-loop correction, however, a better convergence cannot necessarily be expected. Only the inclusion of all loop corrections at this given order and finally also the three-loop corrections can be expected to reduce the theoretical uncertainty on the mass values. We do not display the plots for P2DR corresponding to those shown for P1OS. Their inspection shows the same qualitative behaviour as for P1OS.

4.5 Discussion of the Renormalization Schemes

After having presented our results, we want to finish with a discussion about the renormalization schemes. Compared to our previous results of **[58]** we not only included additionally the two-loop ${\cal O}(\alpha_{t}^{2})$ corrections in our fixed order result to the previously computed ${\cal O}(\alpha_{t}\alpha_{s})$ corrections, we furthermore changed the $\overline{\mbox{DR}}$ renormalization of the top/stop sector. Thus we included also ${\cal O}(\alpha_{t})$ contributions in the conversion of the OS top quark mass to the $\overline{\text{DR}}$ mass at $\mu=M_{Z}$, and we included the resummation of $y_{t}$ in the RGEs. Furthermore, we applied a running of the VEV. We found that in particular the last three changes moved our results with $\overline{\text{DR}}$ renormalization in the top/stop sector close to those of the other codes applying solely $\overline{\text{DR}}$ renormalization, FlexibleSUSY **[75]**, NMSSMTools **[62, 63, 64]**, SOFTSUSY **[65, 66, 138]** and SPHENO **[73, 74]**. The results of these codes and NMSSMCALC (applying the $\overline{\text{DR}}$ scheme) have been compared in **[81]**. In Table 6 we show the two-loop results for the masses of the CP-even scalars for six different test points TP1,…,TP6 using all $\overline{\text{DR}}$ spectrum generators “out-of-the-box”. This table is adapted from Tab. 3 of **[81]222222For more details and the definition of the test points, we refer the reader to [81]. by adding further results for NMSSMCALC, for which we give three different values. We give the result denoted by ’QCDrun+${\cal O}(2,a)$’, with ’QCDrun’ as defined in Subsec. 4.3 for the discussion of Fig. 11, and ’$2,a$’ referring to the two-loop correction of ${\cal O}(\alpha_{t}\alpha_{s})$. It corresponds to our calculation of **[58]** and reproduces the values given for NMSSMCALC in **[81]232323The small differences in the mass values of the $h_{s}$-like Higgs boson given here with respect to those of [81] are due to a correction in its counterterm.. The numbers called ’EWrun+${\cal O}(2,b)$’ are those for our newly computed fixed order two-loop result of ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ after applying the changes in the $\overline{\text{DR}}$ renormalization, described at the beginning of this subsection (and in detail in App. C), and defined as ’EWrun’ in Subsec. 4.3. The third value, with the label ’gaugerun+${\cal O}(2,b)$*’, is obtained at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ after including the gauge contributions in the conversion and running as described in App. D. The NMSSMCALC values for the SM-like Higgs boson that were below those of the other codes by up to slightly more than 3 GeV have moved close to within 1 GeV mainly due to the adaption of our $\overline{\text{DR}}$ renormalization scheme. The exception is point TP6 which features a rather large value for $\lambda$ with $\lambda=1.6$. We are using SM RGEs in the running of $y_{t}$ and $\alpha_{s}$, while other codes (FlexibleSUSY, SOFTSUSY and SPHENO) use NMSSM RGEs for scales $\mu\geq M_{Z}$, so that contributions which become significant for large NMSSM-like couplings are not taken into account in the same way. *

Finally, we discuss our OS results where all top/stop parameters are OS. The OS values of the stop parameters $A_{t},m_{\tilde{Q}_{L}}$ and $m_{\tilde{t}_{R}}$ are obtained from the default $\overline{\text{DR}}$ input parameters through the Eqs. (3.110)-(3.112) without the inclusion of any resummation, whereas the top-quark mass is an OS input parameter and hence does not require a conversion. The large difference between the $\overline{\text{DR}}$ and OS results originates from the inclusion of the $y_{t}$ running in the $\overline{\text{DR}}$ scheme. To verify this statement, we propose a third option for the renormalization of the top/stop sector, calling it OS${}^{\text{mod}}$. In this scheme, we renormalize the three stop parameters OS as before, but use the $\overline{\text{DR}}$ top quark mass. The results are shown in Fig. 13. The figure has been generated by starting from the point P1OS and varying $A_{t}$. The plot shows the mass of the $h_{u}$ dominated SM-like Higgs boson for the $\overline{\text{DR}}$ scheme and the modified OS scheme OS${}^{\text{mod}}$, called for simplicity ’OS’ in the plot, at one-loop, two-loop ${\cal O}(\alpha_{t}\alpha_{s})$ and two-loop ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$. The lower plot depicts for each loop order the relative change $\Delta$ in the mass when changing the renormalization scheme from the $\overline{\text{DR}}$ to the OS${}^{\text{mod}}$* scheme, as defined in Eq. (4.137) with OS replaced by OS*${}^{\text{mod}}$. It is obvious that (at each loop order) the $\overline{\text{DR}}$ and OS curves have moved much closer compared to Fig. 9 with our original OS definition. The relative error due to missing higher-order corrections based on the change between the two renormalization schemes $\overline{\text{DR}}$ and OS${}^{\text{mod}}$* has shrinked considerably and now amounts to less than 1.1% at one-loop level, less than about 0.9% at ${\cal O}(\alpha_{t}\alpha_{s})$ and further improves at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ with less than 0.5%. We conclude our analysis with two statements. First of all, care has to be taken when estimating the remaining theoretical uncertainty. Depending on the applied renormalization scheme, the conclusion drawn can be quite different. Second, we found that the application of a $\overline{\text{DR}}$ top-quark mass and OS stop parameters moves the pure $\overline{\text{DR}}$ result and the ’OS’ result much closer than the use of the OS top-quark mass, as is not surprising.*

5 Conclusions

We computed the fixed order ${\cal O}(\alpha_{t}^{2})$ corrections to the neutral Higgs bosons of the CP-violating NMSSM in the gaugeless limit at vanishing external momentum, thus improving our previous results at ${\cal O}(\alpha_{t}\alpha_{s})$. We applied a mixed $\overline{\text{DR}}$-OS renormalization scheme for the NMSSM input parameters. For the top/stop sector which has to be renormalized at two-loop order, we apply either a $\overline{\text{DR}}$ or an OS definition. The two-loop corrections at ${\cal O}(\alpha_{t}\alpha_{s}+\alpha_{t}^{2})$ are found to amount to a few percent for the SM-like Higgs boson mass. In order to discuss the remaining theoretical uncertainty due to missing higher-order corrections we both vary the renormalization scheme of the top/stop sector and the renormalization scale of the $\overline{\text{DR}}$ parameters. The discussion shows that care has to be taken when drawing conclusions on the theory error. In particular, a modification of the original OS definition of the top/stop sector to the inclusion of a running top-quark instead of an OS mass considerably improves the convergence between the $\overline{\text{DR}}$ and the thus defined modified OS scheme OS${}^{\text{mod}}$. This calls for further improvements in the fixed-order calculation including higher loop orders.

Acknowledgments

We thank Philipp Basler for providing us with a set of NMSSM parameter points. We are grateful to Pietro Slavich for sending us the $\alpha_{t}^{2}$ correction in the real MSSM for comparison. We thank Martin Gabelmann, Michael Spira and Florian Staub for useful discussions. HR’s work is partially funded by the Danish National Research Foundation, grant number DNRF90. This research was supported in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 - TRR 257. TND’s work is funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.01-2017.78. RG acknowledges support of the ’Berliner Chancengleichheitsprogramm’.

Appendix

Appendix A Two-Loop Self-Energy Diagrams

A.1 Two-Loop Self-Energies of the Neutral Higgs Bosons

Figure 14 shows the two-loop self-energies of the neutral Higgs bosons, needed at ${\cal O}(\alpha_{t}^{2})$.

A.2 Two-Loop Self-Energies of the Charged Higgs Boson

Figure 15 displays the generic two-loop self-energies of the charged Higgs boson that contribute to the computation of the two-loop counterterm of the charged Higgs boson mass.

Appendix B Scalar One-Loop Integrals to $\mathcal{O}(\varepsilon)$

In the following, we present the scalar one-loop one- and two-point integrals expanded up to $\mathcal{O}(\varepsilon)$ which are used for our calculation of the two-loop corrected neutral Higgs boson masses.

B.1 Conventions

For logarithms, we use the short-hand notation

[TABLE]

where $m_{i}$ is an arbitrary parameter with mass dimension and $Q$ is the modified renormalization scale, which is related to the renormalization scale $\mu_{R}$ by

[TABLE]

where $\gamma_{E}$ is the Euler-Mascheroni constant. With $\text{Li}_{2}$ we denote the dilogarithm, or Spence’s function, defined for all $z\in\mathbb{C}$ as

[TABLE]

For taking several limits in the one-loop two-point integral, the following specific values are useful,

[TABLE]

*where $\zeta(2)$ denotes the Riemann zeta function $\zeta(s)$ evaluated at $s=2$. *

The solutions of the equation

[TABLE]

needed later for the evaluation of the scalar two-point integral **[139]**, where $p^{2}$ is the squared four-momentum, are denoted by $r_{1}$ and $r_{2}$. They can be cast into the form

[TABLE]

with $\epsilon>0$ and $\epsilon\ll m_{i}^{2}$, $\epsilon\ll p^{2}$. We choose $r_{1}$ to be the solution with the positive sign. The solutions satisfy the relations

[TABLE]

where $x$ and $y$ are dimensionless quantities defined as

[TABLE]

B.2 The Scalar One-Loop One-Point Integral $A_{0}$ at

$\mathcal{O}(\varepsilon)$

The scalar one-loop one-point integral is defined as **[91]**

[TABLE]

in $D=4-2\varepsilon$ dimensions. The expansion of the integral up to $\mathcal{O}(\varepsilon)$ reads

[TABLE]

For the special case $m^{2}=0$, the integral vanishes,

[TABLE]

B.3 The Scalar One-Loop Two-Point Integral $B_{0}$ at $\mathcal{O}(\varepsilon)$

The scalar one-loop two-point integral is defined as **[91]**

[TABLE]

in $D=4-2\varepsilon$ dimensions. The most general solution of $B_{0}$ up to $\mathcal{O}(\varepsilon)$ is given by **[139]**

[TABLE]

This solution is valid for arbitrary values of $p^{2}$, $m_{1}^{2}$ and $m_{2}^{2}$. However, for several limits of the parameters, the computation becomes numerically unstable due to the appearance of large logarithms or vanishing denominators. We therefore give explicit expressions for several limiting cases:

[TABLE]

The parameter region of the last integral is a special case which represents a pure divergence. It can be regularized in the UV regime in the same way as the most general solution of the $B_{0}$ integral. However, for $p^{2}=m_{i}^{2}=0$, it additionally contains a divergence in the infrared (IR) regime, which can be regularized in dimensional regularization with the regulator $\varepsilon_{\text{IR}}$. We want to emphasize that the results of our $\mathcal{O}(\alpha_{t}^{2})$ corrections to the neutral Higgs boson masses in the gaugeless approximation do not depend on the IR regulator.

Appendix C The Running $\overline{\text{DR}}$ Top Mass

In this section, we present for the NMSSM the computation of the $\overline{\text{DR}}$ top quark mass at the SUSY scale from the given top quark pole mass, as it is implemented in NMSSMCALC. We follow the strategy of **[58*]** where we assume that below the SUSY scale the running top quark mass is described solely by the SM, and the NMSSM contributions enter at the SUSY scale. The computation is done in the gaugeless approximation and performed in four steps as follows: *

Step 1: Translation within the SM of the OS parameters to the $\overline{\text{MS}}$ parameters at $\mu_{R}=M_{Z}$:

We include in this translation the one-loop and two-loop QCD corrections **[140, 141]** together with the one-loop EW corrections of ${\cal O}(\alpha_{t})$ in the gaugeless limit, in particular

[TABLE]

with

[TABLE]

where $n_{f}=5$. The one-loop two-point functions are evaluated at the renormalization scale $M_{Z}$. We denote the top quark pole mass by $m_{t}$ and the SM Higgs mass by $M_{H}$, with $M_{H}=125.09\,~{}\text{GeV}$ **[35]**. Our ${\cal O}(\alpha_{t})$ correction in the gaugeless limit is in agreement with the result presented in **[142]**. We use the following convention for the classical SM Higgs potential

[TABLE]

The renormalized SM VEV is defined in such a way that it minimizes the loop-corrected Higgs potential. The $\overline{\text{MS}}$ running VEV at the scale $\mu_{R}=M_{Z}$ is related to the on-shell VEV via

[TABLE]

where $\delta v$ is the finite part of the on-shell VEV counterterm,

[TABLE]

with the finite part of the $W$ and $Z$ boson mass counterterms

[TABLE]

being computed at one-loop level in the gaugeless approximation. Additionally, we need the running Higgs mass at the renormalization scale $M_{Z}$, which we obtain from

[TABLE]

with

[TABLE]

with the scalar one- and two-loop functions evaluated at the scale $M_{Z}$.

The running $\overline{\text{MS}}$ top quark Yukawa coupling $y_{t}$ and the Higgs quartic coupling $\lambda_{\text{SM}}$ at $M_{Z}$ are computed via the relations242424We computed $\delta m_{t}$, $\delta v$ and $\delta_{H}$ in the general $R_{\xi}$ gauge with the full one-loop EW correction and we confirmed that these formula give gauge-independent results for the running top quark Yukawa and quartic Higgs couplings, as expected.**

[TABLE]

where the tadpole at one-loop in the gaugeless approximation is given by

[TABLE]

*The $\overline{\text{MS}}$ running parameters $y_{t}(M_{Z}),\lambda_{\text{SM}}(M_{Z})$ and $v(M_{Z})$ together with $\alpha_{s}(M_{Z})$ are input parameters of the renormalization group equations (RGEs) for the running from $M_{Z}$ to the SUSY scale, which is described in the second step. *

Step 2: RGE running to the SUSY scale $\mu_{R}=M_{\text{SUSY}}$:

We use the two-loop RGEs for $g_{3},y_{t},\lambda_{\text{SM}}$ where $g_{3}$ is the strong gauge coupling with $\alpha_{s}=g_{3}^{2}/(4\pi)$ and $v$ in the gaugeless limit that are relevant for the running of the top quark mass to the SUSY scale, namely **[143, 144, 145, 146]**

[TABLE]

where $n_{f}=5$ for $\mu_{R}<m_{t}$ and $n_{f}=6$ for $\mu_{R}\geq m_{t}$. For the solution of the coupled system of differential equations we use the Runge-Kutta algorithm evaluated to fourth order **[147, 148]**. The $\overline{\mbox{MS}}$ top quark mass at the SUSY scale $\mu_{R}=M_{\text{SUSY}}$ is then obtained from

[TABLE]

Step 3: Conversion from $\overline{\text{MS}}$ to $\overline{\text{DR}}$:

Within the SM, the $\overline{\text{DR}}$ top quark mass is computed from the $\overline{\text{MS}}$ top quark mass at the SUSY scale by using the two-loop relation*[149, 150, 151]***

[TABLE]

Step 4: Adding the NMSSM contributions:

In the final step, the NMSSM $\overline{\text{DR}}$ top quark mass is calculated from the SM $\overline{\text{DR}}$ top quark mass by

[TABLE]

where $dm_{t}^{\alpha_{s}/\alpha_{t}}$ are the SUSY contributions relevant at $\mu_{R}=M_{\text{SUSY}}$ at order ${\cal O}(\alpha_{s})$ and ${\cal O}(\alpha_{t})$, respectively,

[TABLE]

where $m_{t}$ is the NMSSM $\overline{\text{DR}}$ running mass at the SUSY scale and $\varphi_{M_{3}}$ the phase of the gaugino parameter $M_{3}$. Note that we keep only the heavy Higgs boson and Higgsino contributions in $dm_{t}^{\alpha_{t}}$, since the light Higgs boson contribution is identified with the SM Higgs contribution in Eq. (C.165) and has already been taken into account at scale $M_{Z}$.

Appendix D The Running $\overline{\text{DR}}$ Top Mass with the Inclusion of the Gauge Couplings

*Since we want to investigate the effect of the gauge contributions on the loop-corrected Higgs masses, we specify in this section the computation of the running $\overline{\mbox{DR}}$ top quark mass when the gauge contributions are included. Due to their inclusion, gauge dependences appear at several places in the formulae. We therefore performed the computation in the general $R_{\xi}$ gauge and kept the explicit dependence on the gauge parameter $\xi$. In our numerical analysis, however, we use the t’ Hooft-Feynman gauge (i.e. $\xi=1$) to be consistent with the gauge that we use in our one-loop calculation of the Higgs boson masses. In the following, we present analogously to Section C the steps necessary to obtain the running ${\overline{\mbox{DR}}}$ top mass at the scale $M_{\text{SUSY}}$. Since the fourth step does not differ from the gaugeless approximation, we solely describe here the changes in the steps 1-3: *

Step 1: Translation within the SM of the OS parameters to the $\overline{\text{MS}}$ parameters at $\mu_{R}=M_{Z}$:

Including the gauge contribution, the $\overline{\text{MS}}$ top mass at $M_{Z}$ reads

[TABLE]

where the QCD contribution is presented in Eq. (C.164) and the full one-loop EW correction is given by

[TABLE]

where $\xi_{W,Z}$ is the gauge fixing parameter related to the $W$ and $Z$ boson, respectively. For the running $\overline{\text{MS}}$ VEV at the scale $\mu_{R}=M_{Z}$, we still use the relation

[TABLE]

where $\delta v$ including the gauge contribution reads

[TABLE]

Since we use the fine structure constant at the $Z$ boson mass $M_{Z}$, $\alpha=\alpha(M_{Z})$, as an input, the counterterm of the electric coupling is defined as

[TABLE]

where the transverse part of the photon self-energy $\Sigma_{T}^{AA}$ includes only the light fermion contributions. The analytical result for the counterterm of the electric coupling constant is given by

[TABLE]

*where we have set the light fermion masses to zero. The finite parts of the $W$ and $Z$ boson mass counterterms are *

[TABLE]

and

[TABLE]

For the running Higgs mass at the renormalization scale $M_{Z}$, we use

[TABLE]

with

[TABLE]

For the $\overline{\text{MS}}$ top Yukawa coupling and SM coupling $\lambda$ at $M_{Z}$ we again use the Eqs. (C.172) and (C.173) with the one-loop tadpole given by

[TABLE]

We also need the $\overline{\text{MS}}$ SM gauge couplings $g_{1},g_{2}$, with the $U(1)_{Y}$ gauge coupling $g_{1}=\sqrt{5/3}g_{Y}$, where $g_{2}$ and $g_{Y}$ are the $SU(2)_{L}$ gauge coupling and the $U(1)_{Y}$ gauge coupling at the GUT scale, respectively. We remind the reader that we have used the running fine structure constant $\alpha(M_{Z}),M_{Z},M_{W}$ as input parameters, so that the gauge couplings $g_{2}$ and $g_{Y}$ read

[TABLE]

where $c_{W}=M_{W}/M_{Z}$ and $s_{W}=\sqrt{1-M_{W}^{2}/M_{Z}^{2}}$ are OS parameters. As result the running $\overline{\text{MS}}$ couplings $g_{2}$ and $g_{Y}$ at the scale $M_{Z}$ are given by252525Note that the counterterms $\delta s_{W},\delta c_{W}$ are gauge-independent quantities.**

[TABLE]

where

[TABLE]

*with the analytic expressions for $\delta M_{W}$, $\delta M_{Z}$ and $\delta Z_{e}^{\alpha(M_{Z})}$ presented in Eq. (D.190), Eq. (D.191) and Eq. (D.189), respectively. *

Step 2: RGE running to the SUSY scale $\mu_{R}=M_{\text{SUSY}}$:

Here we use the two-loop RGEs for $g_{1},g_{2},g_{3},y_{t},\lambda$ and $v$ **[143, 144, 145, 146]**. The computation of the RGE for the VEV is done with the help of the SARAH package **[69, 70, 71, 72, 49]**. The contributions from the light fermions have been neglected here as well. We have

[TABLE]

Step 3: Conversion from $\overline{\text{MS}}$ to $\overline{\text{DR}}$:* *

With the inclusion of the electroweak gauge contributions the conversion from the $\overline{\text{MS}}$ to the $\overline{\text{DR}}$ top quark mass is modified to

[TABLE]

Appendix E NMSSM RGEs for Investigations of the Scale Dependence

In order to estimate the remaining theoretical uncertainty due to missing higher-order corrections, we investigate their scale dependence. For a consistent investigation also the input parameters have to be evaluated at the respective scales. According to the SLHA the SUSY parameters in the input file are $\overline{\text{DR}}$ parameters evaluated at the user-specified scale, where the SUSY scale is the default choice. These parameters need to be evolved to the different scales used to determine the scale dependence. For this, the RGEs of the relevant parameters are needed. In order to match the approximations used for the computation of the Higgs masses at order $\mathcal{O}(\alpha_{t}^{2})$ and $\mathcal{O}(\alpha_{t}\alpha_{s})$, respectively, we use the RGEs in the gaugeless limit and set the light Yukawa couplings to zero. We computed all RGEs with SARAH and cross-checked them with the results given in **[152, 27, 30]**, and with the results of **[94]** for the RGEs of $\tan\beta$ and $v_{s}$. The RGEs are given by

[TABLE]

Due to the supersymmetric non-renormalization theorem **[153, 154]** the RGEs of the superpotential parameters are proportional to themselves. Hence, the phases $\varphi_{\lambda}$ and $\varphi_{\kappa}$ do not need to be renormalized separately and we can write the RGEs for $\lambda$ and $\kappa$ in terms of their absolute values, see Eqs. (E.211) and (E.212). Instead, for the parameters $A_{\lambda}$ and $A_{\kappa}$ we solve the RGEs for the real part and imaginary part separately,

[TABLE]

For the soft SUSY breaking parameters entering the computation of the Higgs boson masses at one-loop level, we implemented the RGEs accordingly at one-loop order. They are given by

[TABLE]

Like for $A_{\lambda}$ and $A_{\kappa}$, we implemented the RGEs for the real and imaginary parts of $M_{i}$ ($i=1,2,3$) and $A_{j}$ ($j=t,b,\tau$) separately. In addition, we need the RGEs for $m_{H_{u}}^{2}$, $m_{H_{d}}^{2}$ and $m_{S}^{2}$ as they enter the RGEs of $m_{\tilde{Q}_{3}}^{2}$ and $m_{\tilde{t}_{R}}^{2}$. They read

[TABLE]

The RGEs are solved by means of the Runge-Kutta algorithm evaluated to fourth order. Note, that in the mixed $\overline{\text{DR}}$-OS scheme we first convert all the parameters to $\overline{\text{DR}}$ and subsequently solve the full system of RGEs.

Appendix F Tree-level Neutral Higgs Mass Matrix

The matrix elements of the neutral Higgs boson mass matrix in terms of the independent input parameters read

[TABLE]

The gaugeless approximation is obtained by setting $M_{W}=M_{Z}=0$ but keeping $v$ fixed.

Appendix G Neutral Higgs Counterterm Mass Matrix

In the following we present the explicit analytic form of the counterterm mass matrix defined in Eq. (3.70) in the basis $(h_{d},h_{u},h_{s},a,a_{s})$. At one-loop level, the matrix elements of the Higgs counterterm mass matrix are given by

[TABLE]

where the $(\Delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\mathcal{M}_{hh})_{ij}$ are obtained from Eqs. (G.244)-(G.258) by setting $n=1$. At two-loop level, the matrix elements of the Higgs counterterm mass matrix not only contain counterterms at two-loop order but also the product of two one-loop counterterms,

[TABLE]

where the $(\Delta^{\text{\tiny(}\!\text{\tiny 2}\!\text{\tiny)}}\!\mathcal{M}_{hh})_{ij}$ are given in Eqs. (G.244)-(G.258) with $n=2$ and $(\Delta^{{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}\!}{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}\!}}\mathcal{M}_{hh})_{ij}$ in Eqs. (G.259)-(G.273). In the following formulae for the counterterm mass matrix elements, we already applied the gaugeless approximation.

[TABLE]

In the following, we present the matrix elements of the Higgs counterterm mass matrix containing the product of two one-loop counterterms which contribute at two-loop order. All one-loop counterterms $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!v_{s},\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!|\kappa|,$ $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\text{Re}\,A_{\kappa},\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\varphi_{y},\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\varphi_{z}$ are zero.

[TABLE]

Appendix H Charged Higgs Boson Mass Counterterm

In the following we present the explicit analytic form of the counterterm of the charged Higgs boson mass, defined in Eqs. (3.92) and (3.93), which is needed if $\text{Re}\,A_{\lambda}$ is chosen as independent input parameter. At one-loop level, the counterterm of the charged Higgs boson mass is given by

[TABLE]

where the $\Delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!M_{H^{\pm}}^{2}$ is obtained from Eq. (H.276) by setting $n=1$. At two-loop level, the charged Higgs mass counterterm not only contains counterterms at two-loop order but also the product of two one-loop counterterms,

[TABLE]

where the $\Delta^{\text{\tiny(}\!\text{\tiny 2}\!\text{\tiny)}}\!M_{H^{\pm}}^{2}$ are obtained from Eq. (H.276) for $n=2$ and $\Delta^{{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}\!}{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}\!}}M_{H^{\pm}}^{2}$ is given in Eq. (H.277). In the following formulae for the mass counterterm, we already applied the gaugeless approximation and present only the terms that are relevant for the calculation at order $\mathcal{O}(\alpha_{t}^{2})$. The counterterm contribution $\Delta^{(n)}M_{H^{\pm}}^{2}$ reads

[TABLE]

Next, we give the charged Higgs mass counterterm containing the product of two one-loop counterterms which contribute at two-loop order at order $\mathcal{O}(\alpha_{t}^{2})$. All one-loop counterterms $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!v_{s},\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!|\kappa|,$ $\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\text{Re}\,A_{\kappa},\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\varphi_{y},\delta^{\text{\tiny(}\!\text{\tiny 1}\!\text{\tiny)}}\!\varphi_{z}$ are zero. The counterterm contribution reads

[TABLE]