1 Introduction

The discovery of the Higgs boson was the biggest success of Run-I of the Large Hadron Collider (LHC) [1, 2]. The mass of the Higgs is already known very precisely up to a few hundred MeV and its properties are in good agreement with the expectations from the Standard Model (SM). These observations now give strong constraints on any extension of the SM. Therefore, it is necessary to calculate these properties with increasing accuracy to close the gap between the experimental and theoretical uncertainty. In the context of supersymmetric models, most efforts have been put into a precise calculation of the Higgs mass in the minimal supersymmetric standard model (MSSM) assuming real parameters [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37]. In addition, for the next-to-minimal supersymmetric standard model (NMSSM) most calculations of the Higgs mass considered the CP conserving case [38,39,40,41].

Recently, however, there has been increasing interest in theories that extend these minimal models, and this has led to the present authors’ work [42, 43], building on that of [44,45,46,47], to extend those calculations and provide a public implementation for two-loop Higgs mass calculations in generic theories. Up to this point the corrections were only available for CP-even scalars in theories with real parameters; this paper discusses the extension to all neutral scalars with or without CP violation (CPV).

Indeed, the focus on the real versions of the MSSM and NMSSM to study the Higgs mass can hardly be motivated from first principles: there is no strong argument why the CP phases in the soft-breaking sector of SUSY models should be small—especially if SUSY breaking is transmitted via gravity. Moreover, SUSY models with CPV can have very interesting phenomenological aspects; see for instance [48,49,50,51,52,53,54,55,56,57,58]. Putting the phases to zero is often an assumption to circumvent conflicts with experimental limits and simplify calculations. Due to this the impact of CP phases on the Higgs mass has so far only been partially considered in both models. In the MSSM it was first studied using renormalisation group techniques [59,60,61], while diagrammatic calculations at the one-loop level [62] and certain two-loop level contributions [63,64,65] were performed much later. For the complex NMSSM the one-loop results [66] are so far only accompanied by two-loop corrections of \(O(\alpha _s \alpha _t)\) [67].

CPV in supersymmetric theories is constrained by several observations, principally meson mixing (in particular \(K^0\), \(B^0_{s}, B^0_d\) and \(D^0\) mesons) and decays; the electric dipole moment (EDM) of nucleons and electrons; and Higgs coupling measurements. Meson physics typically places extremely stringent bounds, but there is little overlap between those constraints and constraints on the Higgs mass/mixings as relevant for this work, because generational mixing is required—and even then, large enough generational mixing to have a sizeable effect on the Higgs mass at two loops may be unconstrained by flavour [68]. Furthermore, the measurement of the Higgs couplings at the LHC is rather insensitive to parity violation, with the parity-violating couplings still allowed to be of the same order as, if somewhat less than,Footnote 1 the Standard-Model-like ones [69], and so direct searches for additional Higgs bosons actually place more stringent constraints.

Therefore the most relevant constraint on the parameter space that we shall consider comes from electric dipole moments (see e.g. [70, 71]), in particular that of the electron \(\mathrm{d}_e\), which is constrained to be [72]

$$\begin{aligned} |\mathrm{d}_e| < 8.9 \times 10^{-29}\ e\ \mathrm {cm} = 4.5\times 10^{-15} e\ \mathrm {GeV}^{-1}. \end{aligned}$$
(1)

The typical value for electric or chromoelectric moments (CDMs) for fermions i of mass \(m_i\) and a common SUSY scale \(M_\mathrm{SUSY}\) is [71]

$$\begin{aligned} \kappa _i \equiv \frac{m_i}{16\pi ^2 M_{\mathrm{SUSY}}^2} =1.3 \times 10^{-25}\ {\mathrm{cm}} \times \frac{m_i}{\mathrm{MeV}} \times \left( \frac{{\mathrm{TeV}}}{M_{\mathrm{SUSY}}} \right) ^2 \end{aligned}$$
(2)

multiplied by a numerical factor, three Yukawa or gauge couplings, and the sine of a CP-violating phase. In the case of the electron dipole moment in the MSSM with only CP violation entering through the \(\mu \)-term with phase \(\varphi _\mu \) we have

$$\begin{aligned} \mathrm{d}_e/e \simeq \frac{5 g^2 }{24} \kappa _e \tan \beta \sin ( \varphi _\mu ) \end{aligned}$$
(3)

where \(\tan \beta \) is the ratio of up- to down-type Higgs vevs, and we therefore need a large suppression of the phase by roughly three orders of magnitude; however, if we just consider the sneutrino–chargino sector then we obtain in the same limit

$$\begin{aligned} |\mathrm{d}_e/e| \simeq \frac{4 g^2 }{3} \kappa _e \tan \beta \big |\sin ( \varphi _\mu + \varphi _{M_2} + \eta )\big |, \end{aligned}$$
(4)

which similarly constrains more CP-violating phases; here \(\varphi _{M_2}\) is the Wino mass phase and \(\eta \) is that of the up-type Higgs doublet (which will be defined in Sect. 3).

There is a constraint from the neutron electric dipole moments, a recent limit being [73] (see also [74,75,76])

$$\begin{aligned} |\mathrm{d}_n/e| \lesssim 3.0 \times 10^{-26}\ \mathrm {cm} . \end{aligned}$$
(5)

Here the calculation is more complicated, since it depends on the electric dipole moments of the light squarks, their chromoelectric moments \(\tilde{\mathrm{d}}_{u,d}\) (naively suppressed by a factor \(\frac{e}{4\pi }\)) and the theta-angle of QCD. While this has the power to restrict the gluino phase through a squark–gluino loop, either through a direct EDM since the squarks are charged, or through the chromoelectric moments, this is not relevant for our study because the bound is not sufficiently strong: it can easily be satisfied just by, for example, taking the first two generations of squarks to have masses of a few TeV.

There is an additional strong constraint from the mercury dipole moment [71, 77]:

$$\begin{aligned} |\mathrm{d}_\mathrm{Hg}/e| \simeq 10^{-29} \ \mathrm {cm}\ \times \left| \frac{\tilde{d}_u - \tilde{d}_d}{10^{-26}\ \mathrm {cm}} \right| \lesssim 7.4\times 10^{-30} \ \mathrm {cm} . \end{aligned}$$
(6)

If the first two generations of squarks are heavy then, again, they will not constrain the parameter space relevant for the Higgs mass.

Finally, we note that the Higgs sector itself contributes to electric dipole moments through two-loop diagrams: these consist of Barr–Zee diagrams which directly produce quark EDMs/CDMs via Higgs/quark loops and the Weinberg three-gluon operator (\(\mathcal {L} \supset \frac{\mathrm{d}_W}{6} f_{abc} \epsilon ^{\mu \nu \alpha \beta } G_{\alpha \beta }^a G_{\mu \rho }^b G^{c\ \rho }_\nu \) where \(f_{abc}\) are the SU(3) structure constants, \(G_{\mu \nu }^a\) is the field strength), which is also generated by quark/Higgs loops (both charged and neutral Higgs states can contribute).Footnote 2 The contribution from the latter to the neutron EDM is [71]

$$\begin{aligned} |\mathrm{d}_n/e| \sim 4 \times 10^{-26}\ \mathrm {cm} \times \big [ \mathrm{d}_W \times 10^{10} \times (\mathrm {GeV})^{2} \big ]. \end{aligned}$$
(7)

The typical size for the more model-independent of these contributions is (see e.g. [78, 79])

$$\begin{aligned} \tilde{\mathrm{d}}_q^\mathrm{Barr-Zee} \sim&\frac{m_q g_s \alpha _s}{(4\pi )^3 v^2} Z_{1h} Z_{1A} \sim 10^{-26} \ \mathrm {cm}\ \times Z_{1h}Z_{1A} \nonumber \\ \mathrm{d}_W \sim&\frac{g_s \alpha _s}{(4\pi )^3 v^2} Z_{1h} Z_{1A} \sim 10^{-9} \ (\mathrm {GeV})^{-2}\ \times Z_{1h} Z_{1A}, \end{aligned}$$
(8)

where \(Z_{ij}\) is the mixing matrix of the Higgs fields; \(Z_{1h}\) is the amount of the SM-like Higgs h that couples to quarks in the lightest eigenstate (excluding the Goldstone boson; i.e. the light field with couplings most like the Standard Model Higgs) and \(Z_{1A}\) is the amount of pseudoscalar. In addition, there are more model-dependent diagrams involving squarks/charginos instead of quarks/gluons; since we have argued that the electroweakino phases are already heavily constrained, and since we shall consider only squark phases in the stop sector, which already produce dipole moments at one loop, we do not expect these to be relevant. Hence we see that the values of the above diagrams are only weakly constraining given current bounds on Higgs mixing. Given that only one-loop EDMs are available in SARAH, a precise computation and application of the two-loop contributions would be beyond the scope of this work, and we shall restrict ourselves to CP phases for the stop trilinear and gluino mass in the MSSM, which we know from the above should not be excluded. On the other hand, all of the relevant observables are available for the NMSSM in the package NMSSMCALC [80, 81]) and we shall check the above expectations with results from that code.

In summary, in supersymmetric models some sources of CPV in the Higgs sector are required to be small by experiment, but several parameters are essentially unconstrained—which could have a strong impact on the masses of the neutral and charged scalars. In general, it is the electric dipole moment of the electron that will restrict the phases \(\varphi _\mu , \eta \) and the phases of the electroweakinos to be close to zero, so we will not consider their effect on the Higgs masses; on the other hand, we shall treat the phase of the gluino and trilinears in the third generation of squarks to be important free parameters, keeping the first two generations of squarks heavy.

The aim of this work is to present the possibility of calculating the Higgs mass and that of all other neutral scalars in a wide range of supersymmetric models with and without CPV to the same accuracy: an automatised, diagrammatic calculation of the Higgs mass covering CPV at the two-loop level is now available via the combination of the public codes SARAH [82,83,84,85,86,87] and SPheno [88, 89]. This functionality extends the automatised two-loop calculations for the real case presented in Refs. [42, 43]. In general, the calculations are done in the gaugeless limit and neglecting the dependence of the external momenta, i.e. they are competitive with the current state-of-the-art calculations for the complex MSSM, but extend any existing two-loop calculation for other SUSY models by important corrections beyond \(O(\alpha _s \alpha _t)\). We explain in Sect. 2 the underlying methodology used in the calculations and some technical subtleties of the new extension before we present in Sect. 3 the validation of the routines in the presence of complex parameters. In Sects. 4 and 5 we discuss some applications of these routines in the context of the MSSM and NMSSM, before we conclude in Sect. 6.

2 Methodology

The calculation of CP-violating corrections at two loops is now available in SARAH via the diagrammatic approach described in Ref. [43]. Indeed, no modifications are required to the expressions given in that paper. For the computation of masses for CP-odd scalars in CP-conserving theories the same routines also apply; since the formalism in Ref. [43] is given in terms of real scalars, and the CP-odd scalars are just CP-even scalars with different labels. However, once we extend our computations to these cases we find two potential subtleties associated with our method of avoiding the Goldstone Boson Catastrophe.

To remind the reader, this problem, highlighted and resolved in Refs. [22, 90, 91], arises either in the MSSM beyond the gaugeless limit, or in theories beyond the MSSM even in the gaugeless limit, in that the \(\overline{\mathrm {DR}}\) mass-squred of Goldstone bosons may be negative (and the Goldstone bosons themselves have non-zero couplings to the Higgs). The full on-shell mass of course being zero, the \(\overline{\mathrm {DR}}\) Goldstone boson mass parameter is thus of the same order as loop corrections and is small. The problem is that this parameter appears in the loop corrections to the tadpoles and masses of Higgs bosons (and other particles) and the solution for a mass calculation is to include momentum dependence.

However, since this is computationally onerous, our solution (described in Refs. [41, 42]) is to exploit the fact that we work in the gaugeless limit and, in our two-loop calculation, can therefore neglect corrections to the mass proportional to electroweak gauge couplings: we use the full potential \(V_0 + V_1 + V_2|_\mathrm{gaugeless}\) to solve the tadpole equations, and then use the parameters determined from these in our pure gaugeless tree-level potential \(V_0|_\mathrm{gaugeless}\) to determine the masses in our theory. Since we are effectively working in a false minimum the \(\overline{\mathrm {DR}}\) Goldstone masses entering in our two-loop calculation are non-zero and of the order of the electroweak boson masses. This has the effect of taming the problem for most parts of the parameter space.

The first subtlety related to this approach as concerns CP violation is that the Goldstone masses are typically tachyonic, and we retain only the real part of the loop functions. It is legitimate to ask whether the tadpole and self-energy diagrams that we compute really then correspond to the first and second derivatives of the two-loop potential once we introduce CP violation, since the complex parts of the couplings may in principle multiply a complex loop function. However, the terms in the potential, given in Ref. [44], all have the form

$$\begin{aligned} \Delta ^{(2)} V_\mathrm{(given\ topology)}= & {} \bigg ( \mathrm {Real\ product\ of\ couplings}\bigg ) \\&\times \bigg ( \mathrm {loop\ function\ of\ masses}\bigg ) \end{aligned}$$

with the exception of one contribution involving fermions and scalars, given (in terms of Weyl fermions labeled by indices in capitals with masses \(m_I^2\) as the eigenvalues of a mass matrix \(M_{IJ}\), and real scalars labeled by lower-case indices of mass \(m_k^2\), having Yukawa couplings \(y^{IJk}\)) by

$$\begin{aligned} V^{(2)}_{\overline{FF}S} = {1\over 4} y^{IJk} y^{I'J'k} M^*_{II'} M^*_{JJ'} f_{\overline{FF}S} (m^2_I, m^2_J, m^2_k) + \mathrm{c.c.}. \end{aligned}$$

The function \(f_{\overline{FF}S} \) is a two-loop function that is real for positive masses-squared but having complex values for negative ones. However, this should be understood as \(\mathrm {Re} ({1\over 2} y^{IJk} y^{I'J'k} M^*_{II'} M^*_{JJ'}) f_{\overline{FF}S} (m^2_I, m^2_J, m^2_k)\) and so falls into the same class as the other terms. Then, when we take the derivatives of the loop functions, since the masses are real, their derivatives with respect to real scalars are real, and we find that the imaginary part of the derivatives of the potential is always the same as the derivatives of the imaginary part, as we require for the consistency of our approach.

The second subtlety once we calculate the masses of CP-odd scalars, or when we have CP violation which mixes originally even and odd scalars, is that among our scalars we now have (would-be) Goldstone bosons. We must therefore ensure that the final Goldstone boson masses should vanish in the Landau gauge once we add the two-loop corrections to the tree and one-loop terms. To show that this is the case in our approach, let us write the tree-level potential as

$$\begin{aligned} V_0 \equiv \frac{1}{2} m_{ij}^2 S_i S_j + \tilde{V}_0 + \tilde{V}_0^D \end{aligned}$$
(9)

for scalars \(S_i\), where \(\tilde{V}_0^D\) is the gauge-coupling dependent part that vanishes in the gaugeless limit, and \(\tilde{V}_0\) comprises the remaining trilinear and quartic couplings. At tree level, the tadpole equations are used to determine some subset of the \(\overline{\mathrm {DR}}\) mass parameters; we shall denote as \(m_{0,ij}^2\) the solutions to the tadpole equations arising from the full potential \(V_0\) where we take the expectation values \(\langle S_i \rangle = v_i\) to be the solutions at all orders (so that by computing loops we correct the mass-squared parameters). We then compute the particle masses at tree level

$$\begin{aligned}&\mathcal {M}_{0,ij}^2 \equiv m_{0,ij}^2 + \partial _i \partial _j \tilde{V}_0 + \partial _i \partial _j\tilde{V}_0^D \nonumber \\&\mathcal {M}_{0,ij}^2\big |_\mathrm{gaugeless} \equiv m_{0,ij}^2 + \partial _i \partial _j \tilde{V}_0. \end{aligned}$$
(10)

Next we compute the corrections to the effective potential up to two-loop order, and one-loop self-energies using these tree-level masses:

$$\begin{aligned} \Delta V \equiv V_1 (\mathcal {M}_{0,ij}^2) + V_2 (\mathcal {M}_{0,ij}^2\big |_\mathrm{gaugeless}) . \end{aligned}$$
(11)

From these we solve the tadpole corrections so that

$$\begin{aligned} m_{ij}^2= & {} m_{0,ij}^2 - \frac{\delta _{ij}}{v_i} \partial _i \Delta V ,\nonumber \\ \mathcal {M}_{ij}^2 (p^2 )= & {} m_{ij}^2 + \partial _i \partial _j \tilde{V}_0 + \partial _i \partial _j\tilde{V}_0^D + \Pi _{1,ij} (p^2, \mathcal {M}_{0,ij}^2) \nonumber \\&+\, \partial _i \partial _jV_2 (\mathcal {M}_{0,ij}^2\big |_\mathrm{gaugeless}). \end{aligned}$$
(12)

Here \( \Pi _{1,ij}\) is the one-loop self-energy. The masses of the neutral scalars are then found as the eigenvalues of this matrix with \(p^2 = m^2\) via an iterative procedure; however, for the Goldstone bosons, we merely need to verify the presence of a null eigenvector for \(p^2 = 0\), when \( \Pi _{1,ij} (0, \mathcal {M}_{0,ij}^2) = \partial _i \partial _jV_1\). For \(p^2=0\) we can rewrite the above as

(13)

To prove that our procedure retains a massless Goldstone boson, we recall the standard proof: if a potential is invariant under a global symmetry where \(S_i \rightarrow S_i + \alpha _i\), then

$$\begin{aligned} \alpha _i \frac{\partial V}{\partial S_i}=0, \quad \frac{\partial \alpha _i}{\partial S_j} \frac{\partial V}{\partial S_i} + \alpha _i\frac{\partial ^2 V}{\partial S_i \partial S_j} =0. \end{aligned}$$
(14)

This is true order by order in perturbation theory. If we are at the minimum of the potential, then the first term in the second equation vanishes and we have a null eigenvector of the mass matrix given by \(\Delta _i\). However, for the two-loop calculation we are not working at the true minimum of the potential, nor are we using the same potential; instead our tree-level potential has \(\tilde{V}^D\) removed, and we solve the tadpole equations according to

$$\begin{aligned} \frac{\partial }{\partial S_i} \left( \frac{1}{2} \hat{m}_{ij}^2 S_i S_j + \tilde{V}_0 + \hat{V}_1 + V_2 \right) = - \frac{\partial \tilde{V}^D_0}{\partial S_i} \end{aligned}$$
(15)

where we have denoted by \(\hat{m}_{ij}^2, \hat{V}_1\) the masses and one-loop potential, to indicate that they are in the gaugeless limit; note, however, that we do not compute or require \(\hat{V}_1\). Therefore

$$\begin{aligned} \alpha _i \frac{\partial ^2 }{\partial S_i \partial S_j} \bigg ( \frac{1}{2} \hat{m}_{ij}^2 S_i S_j + \tilde{V}_0 + \hat{V}_1 + V_2 \bigg ) = - \frac{\partial \alpha _i}{\partial S_j} \frac{\partial \tilde{V}^D_0}{\partial S_i}. \end{aligned}$$
(16)

The term on the right-hand side of this equation is necessary to give a mass to the Goldstone boson at tree level to mitigate the Goldstone boson catastrophe. However, all that we take from this comptuation is the derivatives of the two-loop potential; since we solve for the solution in the same way at tree level and at two loops, and because the potential \( \frac{1}{2} m_{ij}^2 S_i S_j + \tilde{V}_0\) is invariant under the global symmetry (even if the minimum we choose is not), we find

$$\begin{aligned}&\alpha _i \frac{\partial ^2 }{\partial S_i \partial S_j} \bigg ( \frac{1}{2} m_{0,ij}^2 S_i S_j + \tilde{V}_0 \bigg ) = - \frac{\partial \alpha _i}{\partial S_j} \frac{\partial \tilde{V}^D_0}{\partial S_i}, \qquad \nonumber \\&\hat{m}_{ij}^2 = m_{0,ij}^2 - \frac{1}{v_i} \delta _{ij} ( \partial _i \hat{V}_1 + \partial _i V_2)\big |_{S_i = v_i} \nonumber \\&\quad \rightarrow \alpha _i \big (\Delta _2 \mathcal {M}_{ij}^2\big ) = 0. \end{aligned}$$
(17)

Since the one-loop computation used in actually calculating the scalar masses is performed in the minimum of the full potential it will automatically have massless Goldstone bosons; then \(\alpha _i \mathcal {M}_{ij}^2(0) = 0\) when including all corrections as required. It is a highly non-trivial check of our implementation that this should be true; we show this check of our code in the next section and find that it is satisfied to a high level of accuracy.

3 Validation

3.1 Comparison with CP-preserving case

The first verification of our new routines is to compare with the CP-conserving case. In Fig. 1 we show the one- and two-loop lightest Higgs masses obtained in our code as we vary the trilinear CP-violating phase \(\varphi _u \equiv \mathrm {arg}(T_u^{3,3})\) for a point in the complex MSSM; all definitions and other parameter values are given in Sect. 4.4, we take \(M_3 = 2\) \(\text {TeV}\). On the same plot we show an interpolation between the values in the CP-preserving MSSM for the values \(\varphi _u = 0,\pi \) corresponding to \(T_u^{3,3} = \pm |T_u^{3,3}\). Clearly the perfect agreement between the curves at the mid- and end-points indicates the agreement between the two codes. For the rest of the curve, the interpolation is

$$\begin{aligned} m_h^\mathrm{interpolation}= & {} m_h^\mathrm{1\ loop} (\varphi _u) + m_h^\mathrm{2\ loops}\big |_{\varphi _u = 0} - m_h^\mathrm{1\ loop}\big |_{\varphi _u = 0} \nonumber \\&+\, (\varphi _u/\pi )^2 \bigg [ m_h^\mathrm{2\ loops}\big |_{\varphi _u = \pi } \nonumber \\&-\, m_h^\mathrm{1\ loop}\big |_{\varphi _u = \pi } - m_h^\mathrm{2\ loops}\big |_{\varphi _u = 0} + m_h^\mathrm{1\ loop}\big |_{\varphi _u = 0}\bigg ].\nonumber \\ \end{aligned}$$
(18)

For this point, the two-loop corrections are clearly well modeled by a quadratic; we shall investigate more interesting cases in Sect. 4.4.

Fig. 1
figure 1

Plot of the lightest Higgs boson mass in the MSSM, for the parameters given in Sect. 4.4 with \(M_3 = 2000\) \(\text {GeV}\), as the phase of the trilinear soft terms is varied. In red is the result of the new CP-violating code; the blue dashed curve is a quadratic interpolation of the values \(\phi _u = 0, \pi \) from the Higgs mass calculated at two loops in the CP conserving routines, where the difference is added to the one-loop CP-violating case. As can be seen we find perfect agreement for the two routines at the values \(\phi _u = 0,\pm \pi \)

Full size image

3.2 Check of the Goldstone mass

Fig. 2
figure 2

The calculated mass for the neutral Goldstone boson in the MSSM in Landau gauge as a function of im(\(M_3\)). The green line corresponds to the correct calculation. The red lines show the impact of a possible inconsistency in the two-loop calculation: for the dashed line, the phase of \(T_u\) was only included in the calculation of the masses and rotation matrices entering the loop calculation, but dropped in the vertices. For the full line, the phase of the rotation matrix \(Z_U\) was put to zero. For the dotted line the phase of a single vertex in a single diagram involving gluino and (s)down (!) squarks was swapped. The other input parameters are \(m_0 = \text {re}(M_{1/2}) = 1\) TeV, \(A_0 = -2\) TeV, \(\tan \beta =10\), \(\text {sign}(\mu )>0\)

Full size image

As discussed in the previous section, there is an obvious but non-trivial check for the self-consistency of the entire loop calculation: the Goldstone boson mass has to be correct. Thus, choosing Landau gauge, the lightest eigenstate of the four neutral scalars must have zero mass. To obtain this in the complex case, a delicate cancellation of all phases appearing in the mass calculation of the fields in the loops, the phases in the vertices and the combination of the vertices in each diagram must happen. The impact of potential small inconsistencies in these calculations is demonstrated in Fig. 2 where we added by hand some mistakes: dropping imaginary parts of couplings only in the vertex calculation, neglecting the imaginary parts of the squark rotation matrices in the vertices, adding a wrong complex conjugation to a single vertex in a single diagram. For the last point to illustrate the delicacy of the cancellations we have chosen a diagram only involving down squarks and not up squarks, which would have given an even much larger effect. While these mistakes have no impact on the results in the real case, one sees that they immediately spoil the prediction for the neutral Goldstone mass as soon as CP violation is turned on, and this provides a sensitive check for the correctness of our results.

3.3 Comparison with the known NMSSM corrections

The next check is to compare with existing results in the literature. For the MSSM with CP violation the codes FeynHiggs [92] and CPsuperH [93] exist. However, the two codes use another renormalisation scheme compared to SPheno. Therefore, already differences in the real case are present which are often larger than the expected effects from CP phases. Therefore, a quantitative comparison is not possible. The only other public code which supports CP violation is NMSSMCALC for the (scale-invariant) NMSSM [80]. NMSSMCALC makes use of mixed \(\overline{\mathrm {DR}}\)–OS renormalisation conditions for the computation of the Higgs masses, but the OS effects in the (s)top sector can be turned off. This option together with some modifications described in the following allows for a very precise comparison.

The Higgs mass calculation at one-loop level is performed in NMSSMCALC as in SPheno including the full momentum dependence and all possible contributions [40, 66]. At the two-loop level only the \(\mathcal{O}(\alpha _s \alpha _t)\) corrections are included [67]. The missing two-loop corrections will lead inevitably to a difference between SPheno and NMSSMCALC. Moreover, as has been discussed in detail in Ref. [94] the determination of the running \(\overline{\mathrm {DR}}\) parameters entering the Higgs mass calculation also differs between both codes. Therefore, to have a meaningful comparison between codes in the case of CPV, we made the following modifications:

  • SARAH 4.8.3 and SPheno 3.3.8:

    1. 1.

      All two-loop corrections but the ones \(\mathcal{O}(\alpha _s \alpha _t)\) were turned off.

    2. 2.

      The default input using SLHA-2 conventions [95] were changed to SLHA-1 conventions [96] for simpler comparison with NMSSMCALC.

    3. 3.

      The tadpole equations were modified to be solved for im\((A_\kappa )\) and im\((A_\lambda )\) instead of im\((T_\kappa )\) and im\((T_\lambda )\) at tree level: NMSSMCALC solves the tree-level tadpole equations to calculate im\((A_\kappa )\) and im\((A_\lambda )\), but it calculates the radiative shifts to im\((T_\kappa )\) and im\((T_\lambda )\), i.e. solves the loop-corrected tadpole equations with respect to other parameters than the tree-level ones. The SPheno code produced by SARAH always solves the tadpole equations at tree- and loop-level for the same parameters. This would have already given some difference at the one-loop level for specific complex phases, in particular for complex \(\kappa \).

    4. 4.

      The complex phases in the Yukawa couplings, which for instance appear via thresholds in the case of complex \(M_3\), were always put to zero because NMSSMCALC supports only real Yukawa couplings.

  • NMSSMCALC 2.0:

    1. 1.

      A flag to calculate only tree-level masses has been included.

    2. 2.

      The internal calculation of the running SM parameters has been overwritten. Instead, the values are now read in from the input file. This makes it possible to use exactly the same values as SPheno calculates from Standard Model inputs given in appendix A of [88].

    3. 3.

      The finite shifts to \(g_1\), \(g_2\) and v were put to zero.

    4. 4.

      We fixed a bug in the two-loop calculation which we found during our comparison.Footnote 3

Fig. 3
figure 3

Effect of the one-loop corrections (first row) and two-loop corrections (second row) on the masses of the three lightest scalars as a function of im(\(\lambda )\). The blue lines correspond to SPheno, the red ones to NMSSMCALC

Full size image

The conventions for the phases of the (complex) Higgs doublets \(H_u, H_d\) and singlet S, decomposed into real scalars \((\phi _i,\sigma _i)\) and vevs \(v_i\) where \(i=\{u,d,S\}\), are

$$\begin{aligned}&H_d \equiv \left( \begin{array}{c} \frac{1}{\sqrt{2}}(v_d + \phi _d + i \sigma _d) \\ H_d^- \end{array} \right) ,\qquad \nonumber \\&\quad H_u \equiv e^{i\eta } \left( \begin{array}{c} H_u^+ \\ \frac{1}{\sqrt{2}}( v_u + \phi _u + i \sigma _u) \end{array} \right) ,\qquad \nonumber \\&\quad S \equiv e^{i\eta _S} \frac{1}{\sqrt{2}}( v_S + \phi _S + i \sigma _S), \end{aligned}$$
(19)

where \(\eta \) is used as input, and \(\eta _S\) is calculated from the complex input of \(\mu _\mathrm{eff}\) and \(\lambda \) via

$$\begin{aligned} \eta _S = \text {arg}(\mu _\mathrm{eff}) - \text {arg}(\lambda ). \end{aligned}$$
(20)

As default point we have chosen

$$\begin{aligned} \lambda= & {} 0.6,\quad \kappa =-0.3\,,\quad \,A_\kappa =200~{\text {GeV}},\quad \nonumber \\ A_\lambda= & {} 1000~{\text {GeV}},\quad \tan \beta =2.5,\quad \mu _\mathrm{eff}=250~{\text {GeV}},\nonumber \\ M_1= & {} M_2 = 1000~{\text {GeV}},\quad M_3 = 1500~{\text {GeV}},\quad \nonumber \\ A_t= & {} 1500~{\text {GeV}},\quad m_{\tilde{t}_L} = m_{\tilde{t}_R} = 1000~{\text {GeV}}. \end{aligned}$$
(21)

All other sfermion soft masses were put to 1.5 TeV, and all other A-terms to zero.

In Fig. 3 we compare the radiative corrections to the three lightest scalars as a function of im(\(\lambda \)), while in Fig. 4 the impact of im(\(\mu _\mathrm{eff}\)) and im(\(\kappa \)) is shown. Finally, Fig. 5 depicts the dependence on im(\(M_3\)), im(\(A_t\)) and \(\eta \).

Fig. 4
figure 4

Effect of the one-loop corrections (first row) and two-loop corrections (second row) on the mass of the lightest scalars as a function of im(\(\mu _\mathrm{eff}\)) (left) and im(\(\kappa \)) (right). The blue lines correspond to SPheno, the red ones to NMSSMCALC

Full size image
Fig. 5
figure 5

Effect of the one-loop corrections (left) and two-loop corrections (right) on the mass of the lightest scalars as function of im(\(M_3\)) (first row), im(\(A_{t}\)) (second row) and \(\eta \) (third row). The blue lines correspond to SPheno, the red ones to NMSSMCALC

Full size image

We overall find very good agreement at the one- and two-loop level. The differences are at most \(O(10~\text {MeV})\), which corresponds to the agreement obtained in Ref. [94] for the real case and does not increase in the presence of very large phases. Even if it is only possible to compare the corrections \(O(\alpha _s \alpha _t)\) already the large majority of generic possible diagrams is covered. In particular all generic diagrams involving fermions are included, i.e. this confirms the correct treatment of the terms \(V^{(2)}_{\overline{FF}S}\) as discussed in Sect. 2.

Finally, while this section is intended to discuss validation of the code (so we chose parameter ranges that may not be physically interesting), we also used NMSSMCALC to compute the constraints on the CP-violating observables. For this point we found (as expected from the introduction) that \(\mathrm {Im}(M_3)\) is unconstrained, and we also found no constraint upon \(\mathrm {Im}(\kappa )\) for the values considered. However, as discussed in the introduction around Eq. (4), \(|\eta |\) and \(|\mathrm {Im}(\mu _\mathrm{eff})|\) are tightly constrained to \(\lesssim 0.007\) and \(\lesssim 2\) GeV, respectively, from one-loop processes. The two-loop processes become relevant to restrict \(|\mathrm {Im}(A_t)| \lesssim 1600\) GeV from the neutron EDM, and \(|\mathrm {Im}(\lambda )| \lesssim 0.07\) from the two-loop electron EDM.

4 MSSM

In SARAH, we express the Higgs doublets in the MSSM in the usual form shown in Eq. (19). Since we are considering CP violation, the \(\mu \)-term and holomorphic soft-breaking parameters are allowed to be complex, and in general all of the neutral scalars \(\phi _{u,d}, \sigma _{u,d}\) mix, with one state yielding the Goldstone boson of the Z. We use the pure \(\overline{\mathrm {DR}}\) renormalisation scheme, and once we use the measured values of the lepton and quark masses, electroweak gauge couplings and weak mixing angle, Z boson mass, to determine the corresponding \(\overline{\mathrm {DR}}\) quantities (Yukawa couplings, gauge couplings, Fermi constant), we still have a choice of parameters to be eliminated by the scalar tadpole equations. If we define the holomorphic Higgs mass term \(\mathcal {L} \supset - B_\mu H_u \cdot H_d\) and \(B_\mu \equiv e^{i\varphi _{B_\mu }} |B_\mu |\) these read

$$\begin{aligned} 0= & {} m_{H_d}^2 + \frac{1}{2} c_{2\beta }^2 M_Z^2 + |\mu |^2 - t_\beta |B_\mu |\cos (\eta + \varphi _{B_\mu }) \nonumber \\&+ \frac{1}{v_d} \frac{\partial \Delta V}{\partial \phi _d}, \end{aligned}$$
(22)
$$\begin{aligned} 0= & {} m_{H_u}^2 -\frac{1}{2} c_{2\beta }^2 M_Z^2 + |\mu |^2 - \frac{|B_\mu |}{t_\beta }\cos (\eta + \varphi _{B_\mu }) \nonumber \\&+ \frac{1}{v_u} \frac{\partial \Delta V}{\partial \phi _u}, \end{aligned}$$
(23)
$$\begin{aligned} 0= & {} \sin (\eta + \varphi _{B_\mu }) |B_\mu | + \frac{1}{v_u} \frac{\partial \Delta V}{\partial \sigma _d} = \sin (\eta + \varphi _{B_\mu }) |B_\mu | \nonumber \\&+ \frac{1}{v_d} \frac{\partial \Delta V}{\partial \sigma _u}. \end{aligned}$$
(24)

The last two equations are not independent due to the gauge symmetries. Note also that \(\eta \) and \(\varphi _{B_\mu }\) are not independent: in the above they appear in the combination \(\eta + \varphi _{B_\mu }\) so at tree level \(\eta = - \varphi _{B_\mu }\).

In the on-shell scheme used by Refs. [62, 64, 65], the charged Higgs mass and \(\mu \) are taken as input parameters; this is equivalent to specifying \(|B_\mu |\) in the above equations and using the third tadpole equation to determine \(\eta + \varphi _{B_\mu }\). This has the advantage of using a physical input, but the disadvantage of disguising potentially large tuning in the underlying values, particularly for small \(\tan \beta \) where the loop corrections to the Higgs mass must be large and typically lead to large corrections to the other Higgs masses too, as we shall see in the following. In the \(\overline{\mathrm {DR}}\) scheme \(\eta \) (hence \(\eta + \varphi _{B_\mu }\)) is not a fundamental parameter, but rather something that should be derived. On the other hand, it appears in the Higgs couplings and would therefore be complicated to solve for, requiring a computationally expensive iterative procedure and problems with non-zero Goldstone boson masses (since we would be violating the condition that the shift is linear in a mass-squared parameter required in Sect. 2). Instead, we fix \(\eta \) and use the tadpole equations to determine \(B_\mu \). Since the tadpole corrections are typically small compared to \(B_\mu v_u\) this is a small adjustment and we can regard our choice of \(\eta \) as being, to a good approximation, equivalent to minus the phase of \(B_\mu \). For expediency due to the very strong constraints upon it we take \(\eta =0\); we must then regard this as a tuning between \(B_\mu \) and the other CP-violating phases in the theory for large values of \(\varphi _{B_\mu }\).

There then remain two options: one conventional choice, as in the CP-conserving case, is to solve the tadpole equations for \(|\mu |^2\) and \(B_\mu \); this is appropriate when we have GUT-scale boundary conditions where we expect the other soft masses to unify (and we shall use this choice in Sect. 4.3). For our study with low-energy boundary conditions in Sect. 4.4 we shall choose to solve for \(m_{H_d}^2, m_{H_u}^2\) and \(\mathrm {Im}(B_\mu )\), taking \(\mu \) and \(\mathrm {Re}(B_\mu )\) as input parameters. This has the advantage that the tree-level heavy-Higgs masses are simply fixed by \(\frac{|B_\mu |}{\sin \beta \cos \beta }\) so is closer to the on-shell interpretation; but as we shall see the loop corrections can be so large as to render a direct comparison of calculations in the two approaches impossible. In fact, this is further exacerbated since Refs. [64, 65] use the charged Higgs mass as the on-shell parameter, and we are only able (so far) to calculate the charged Higgs mass to one-loop order compared to two loops in that reference.

4.1 One-loop masses and tadpoles

The stop/top sector dominates the one-loop corrections as in the CP-conserving case. If we are only concerned with the lightest Higgs mass in the decoupling limit, then the leading corrections in the top Yukawa coupling \(y_t\) in the effective potential approximation are

$$\begin{aligned} (m_h^\mathrm{approx})^2= & {} M_Z^2 c_{2\beta }^2 + \frac{3}{32\pi ^2 v^2} \bigg [ s_{2\theta }^2 (m_{\tilde{t}_1}^2 - m_{\tilde{t}_2}^2)^2 \nonumber \\&+ 8 m_t^2 \log \frac{m_{\tilde{t}_1}^2m_{\tilde{t}_2}^2}{ m_t^4} + \frac{s_{2\theta }^2}{2} (m_{\tilde{t}_1}^2 - m_{\tilde{t}_2}^2)\nonumber \\&\times \, \bigg ( 8 m_t^2 - s_{2\theta }^2 (m_{\tilde{t}_1}^2 + m_{\tilde{t}_2}^2) \bigg )\log \frac{m_{\tilde{t}_1}^2}{m_{\tilde{t}_2}^2}\bigg ] \end{aligned}$$
(25)

where \(m_{\tilde{t}_{1,2}} \) are the masses of the stop eigenstates, \(m_t\) is the top mass and the square of the sine of twice the mixing angle \(\theta \) is defined as

$$\begin{aligned} s_{2\theta }^2 \equiv \frac{2 v^2 |e^{i\eta } T_{u}^{3,3} s_\beta - y_t \mu ^* c_\beta |^2}{(m_{\tilde{t}_1}^2 - m_{\tilde{t}_2}^2)^2}. \end{aligned}$$
(26)

Note that the stop mixing includes an additional phase (compared to the CP-even case) corresponding to the phase of \(e^{i\eta } T_{u}^{3,3} s_\beta - y_t \mu ^* c_\beta \), but we do not need that here. The important observation here is that if we define \(\mu \equiv e^{i \varphi _\mu }|\mu |, T_u^{3,3} \equiv e^{i\varphi _u}|T_u^{3,3}|\) then the phases only enter through the combination \(\cos (\eta + \varphi _u + \varphi _\mu )\), and the result should be, to leading approximation, even in that combination of phases. Since our two-loop calculations are performed in the effective potential approach in the gaugeless limit, then this should also be true at two loops.

The tadpole contribution from the stops is also important to our calculation, in particular the tadpole for the \(\sigma _u,\sigma _d\) fields. We find that the one-loop stop contribution to these tadpoles neglecting the gauge couplings is

$$\begin{aligned} \frac{\partial \Delta V}{\partial \sigma _u}&\supset&- \frac{3v}{16\pi ^2} \mathrm {Im}\bigg [T_u^{3,3} e^{i\eta } \big ( e^{i\eta } T_u^{3,3} s_\beta - y_t \mu ^* c_\beta \big )^* \bigg ] \nonumber \\&\times \left( \frac{{\mathbf A}_0(m_{\tilde{t}_1}^2) - {\mathbf A}_0(m_{\tilde{t}_2}^2)}{ m_{\tilde{t}_1}^2 - m_{\tilde{t}_2}^2}\right) \nonumber \\ \frac{\partial \Delta V}{\partial \sigma _d}&\supset&- \frac{3v}{16\pi ^2} \mathrm {Im}\bigg [\mu ^* y_t \big ( e^{i\eta } T_u^{3,3} s_\beta - y_t \mu ^* c_\beta \big )^* \bigg ] \nonumber \\&\times \left( \frac{{\mathbf A}_0(m_{\tilde{t}_1}^2) - {\mathbf A}_0(m_{\tilde{t}_2}^2)}{ m_{\tilde{t}_1}^2 - m_{\tilde{t}_2}^2} \right) \end{aligned}$$
(27)

where we define \({\mathbf A}_0 (x) \equiv -x(\log x/Q^2 -1)\), Q being the renormalisation scale. The function of the masses on the right is a slowly varying function with a typical value of order unity, so with \(\eta =0\) we have

$$\begin{aligned} \mathrm {Im} (B_\mu ) \sim&\frac{3}{16\pi ^2} | y_t \mu T_u^{3,3} | \sin \beta \sin (\varphi _u + \varphi _\mu ) . \end{aligned}$$
(28)

So for \(\tan \beta = 5\) (for example) and \(|T_u^{3,3}| = |\mu | = 2000\ {\text {GeV}}, y_t \sim 0.9\) and maximal CP violation in the combination of phases on the right-hand side we have \(\mathrm {Im} (B_\mu ) \sim (270\ {\text {GeV}})^2 \). For a purely imaginary \(B_\mu \) this would correspond to tree-level charged/heavy Higgs mass of \(600\ {\text {GeV}}\); hence charged Higgs masses below \(600\ {\text {GeV}}\)—or, more realistically, somewhat heavier—invoke additional fine-tuning and are difficult to impose in the \(\overline{\mathrm {DR}}\) scheme. In particular, this contributes to the fact that we cannot in any way reliably compare the results of our code to the benchmark scenarios of Refs. [64, 65], which involve lighter charged Higgs masses and larger trilinear couplings than we have quoted here.

4.2 Alternative approach to the Higgs sector

The MSSM is a special case as regards the two-loop mass computations in the gaugeless limit: there is no Goldstone boson catastrophe. To understand this, note that the tree-level Higgs potential consists only of the mass terms; the quartic couplings being given by the gauge couplings that are turned off. Hence the scalar masses are independent of the scalar expectation values in this limit. Now, the potential itself contains no divergences when taking the Goldstone boson masses to zero—the singularities only appear in derivatives of the potential with respect to the Goldstone boson masses—and so the derivatives of the two-loop potential with respect to the scalar expectation values are finite. Hence we are free to consistently use the tree-level solution of the tadpole equations and mass matrices for the Higgs sector in the gaugeless limit, as was done for the calculations in the CP-conserving MSSM in Refs. [20, 21].

We can then write the neutral scalar mass matrix \(\mathcal {M}^2_{h} \) and charged Higgs mass matrix \(\mathcal {M}^2_{H^{\pm }}\) in the gaugeless limit in the basis \(\phi _d, \phi _u, \sigma _d, \sigma _u\) as

$$\begin{aligned} \mathcal {M}^2_{h}= & {} \left( \begin{array}{cc} -\mathcal {M}_2^2 (-\tan \beta ) &{} 0_{2\times 2} \\ 0_{2\times 2} &{} \mathcal {M}_2^2 (\tan \beta ) \end{array} \right) ,\quad \quad \nonumber \\ \mathcal {M}^2_{H^{\pm }}= & {} \mathcal {M}_2^2 (\tan \beta ), \end{aligned}$$
(29)

where

$$\begin{aligned} \mathcal {M}_2^2 (x) \equiv&\quad \left( \begin{array}{cc} |B_\mu | x &{} |B_\mu | \\ |B_\mu | &{} \quad |B_\mu |/x \end{array} \right) . \end{aligned}$$
(30)

Using these mass matrices in the diagrammatic two-loop routines thus neatly avoids tachyonic masses in the two-loop functions, although it does not significantly affect the result.

Fig. 6
figure 6

Plots of Higgs masses as the trilinear phase \(\varphi _A\) is varied. On the left we show the lightest Higgs mass at one (dot-dashed) and two (solid) loops. On the right we show the masses of the heavier eigenstates, which differ for each model by less than one GeV and so the distinction is not visible in the plot; the lower, blue curve is for \(\tan \beta =10\) while the upper, red curve is for \(\tan \beta = 5\)

Full size image

4.3 Two-loop Higgs mass with GUT boundary conditions

If we take CMSSM boundary conditions, that is, a unified A-term parameter \(A_0\), scalar mass parameter \(m_0\) and gaugino masses \(M_{1/2}\), we should solve for \(|\mu |^2\) and \(B_\mu \) (both real and imaginary parts) at low energies (fixing the phase of \(\mu \) as a choice). Then due to the strong constraints from electric dipole moments the phases of \(\eta , \mu , m_{12}\) are constrained to be very small (of order \(10^{-3} \div 10^{-2}\)) and thus not interesting parameters for the Higgs mass; we are only left with the phase of \(A_0\). However, without performing a detailed scan to search for tuned corners of the parameter space, typically points that match LHC constraints on squarks and gluinos, alongside reproducing the correct Higgs mass, will tend to have large values of \(B_\mu \) and thus heavy additional Higgses, showing little CP-violating effects.

We illustrate this in Fig. 6, where we define \(A_0 \equiv -|A_0|e^{i\varphi _{A}}\) and vary the phase for two points:

\(\tan \beta \)

\(|A_0|\) (GeV)

\(m_0\) (GeV)

\(M_{1/2}\) (GeV)

\(m_{h_1}|_{\varphi _A =0}\) (GeV)

\(m_{h_{2}}|_{\varphi _A =0}\) (GeV)

\(m_{h_{3}}|_{\varphi _A =0}\) (GeV)

5

4500

2000

2000

125.0

4118

4118

10

3300

1500

1000

125.5

2494

2495

Here we have given the values that we compute at two loops for the three Higgs neutral Higgs scalars with \(\varphi _A = 0\). In the figure we show the mass of the lightest Higgs at one and two loops as we vary \(\varphi _A\); we also show the heavier Higgs masses whose differences are less than one GeV, either between one and two loops or between the second and third eigenstates.

Fig. 7
figure 7

Electric dipole moment \(d_e\) divided by electric charge e as the trilinear phase \(\varphi _A\) is varied, with exclusion bands shown, for values of \(\tan \beta = 5\) (red) and 10 (blue)

Full size image

We selected small \(\tan \beta \) values to maximise the visibility of the CP-violating effects, because in that way we can have a large difference between \(\varphi _A =0\) and \(\varphi _A = \pi \). This also requires large values of the trilinear couplings, and the contribution from stops to the Higgs mass must be large to obtain the correct value of the Higgs mass for at least some value of \(\varphi _A\). However, this means that for typical points of the parameter space \(B_\mu \) and \(\mu \) will also be large, leading to heavy additional Higgses with small corrections to their masses at two loops. These states then have little impact on the light Higgs mass calculation and so the net effect is still little variation (of about 1 GeV) in the two-loop contribution to the Higgs mass between \(\varphi _A =0 \) and \(\varphi _A = \pi \); almost all of the variation shown in the two plots of Fig. 6 is due to the one-loop effects.

We finally note that the effect of large phases in the trilinears leads to gaugino phases through RGE running, and this in turn has a significant effect on the electron EDM \(d_e\); we show the variation of this in Fig. 7 and note that a small region near \(\varphi _A = \pm \frac{\pi }{2}\) is already excluded for our \(\tan \beta = 10\) scenario (recall that the contributions are proportional to \(\tan \beta \) in e.g. Eq. (4)). Thus we conclude that for CMSSM boundary conditions the CP-violating phases are well constrained and the specific two-loop corrections to the Higgs masses are less important. In the next subsection we shall consider more general low-energy boundary conditions.

4.4 Two-loop shifts with SUSY-scale boundary conditions

If we take our boundary conditions at low energies (i.e. at the masses of the squarks and gauginos) then, without a particular bias for the conditions at high energies, we are free given the constraints to consider a phase of the trilinear terms and the gluino. Here we shall restrict our attention to the stop trilinear term \(T_u^{3,3}\) and maximise the effect of CP violation in the stop sector since this is typically the source of the largest contributions at two loops. Thus again we shall consider a small \(\tan \beta \) scenario; we take as our parameter values

$$\begin{aligned}&T_u^{3,3} = 0.9 A_0 e^{i\varphi _u},\quad T_d^{3,3} = 0.064 A_0 ,\quad T_e^{3,3} = 0.05 A_0 , \nonumber \\&A_0 = 3800\ {\text {GeV}},\quad \mu =3800\ {\text {GeV}}, \quad \tan \beta = 5,\quad \nonumber \\&\mathrm {Re}(B_\mu ) = 10^6\ {\text {GeV}}^2 \left( \frac{\tan \beta }{1 + \tan ^2 \beta }\right) \nonumber \\&m_{\tilde{L},ii}^2=m_{\tilde{e},ii}^2= m_{\tilde{q},ii}^2 = m_{\tilde{u},ii}^2= m_{\tilde{d},ii}^2\nonumber \\&\quad = (3\times 10^3\ {\text {GeV}})^2, \quad i=1,2 \nonumber \\&m_{\tilde{L},33}^2=m_{\tilde{e},33}^2= m_{\tilde{d},33}^2 = (10^3\ {\text {GeV}})^2,\quad \nonumber \\&\quad m_{\tilde{q},33}^2 = (1.6\times 10^3\ {\text {GeV}})^2,\quad m_{\tilde{u},33}^2= (1.52\times 10^3\ {\text {GeV}})^2 \nonumber \\&M_1 =200\ {\text {GeV}}, \quad M_2=500\ {\text {GeV}},\quad M_3 = M_3^0 e^{i\varphi _{M_3}}. \end{aligned}$$
(31)

Here the prefactors of the trilinears are approximate values for the Yukawa couplings to allow simpler comparison between our points and those used in other work such as [64, 65]. The choice of the real part of \(B_\mu \) (here we solve the tadpole equations for \(m_{H_u}^2, m_{H_d}^2,\mathrm {Im}(B_\mu )\)) is such that at tree level the heavy Higgs masses \(m_{h_{2,3}}\) are at one \(\text {TeV}\).

4.4.1 Variation of gluino phase

We expect that the gaugino phase \(\varphi _{M_3}\), only entering the Higgs mass calculation at two loops in the \(\overline{\mathrm {DR}}\) scheme, should be an important parameter for our results. We show in Fig. 8 the effect that it has on the three neutral Higgs masses and the parameter \(\varphi _{B_\mu }\), for \(\varphi _{u} = 0\) and \(\varphi _u = \pi \). The first observation is that the difference between one and two loops is strongly dependent on \(\varphi _{M_3}\); for \(\varphi _u = 0\) this changes between 4 and 7 \(\text {GeV}\), and for \(\varphi _u = \pi \) between 2 and \(-5 \) \(\text {GeV}\)—an overall shift of 7 \(\text {GeV}\) in the latter case! While this point has been chosen to show a large variation, it underlines the importance of the two-loop corrections.

Looking more closely, we find that the effect at two loops is partly to compensate for the variation at one loop, giving a more constant value for all Higgs fields. We also have the potentially counter-intuitive result that the \(\varphi _u = \pi \) points have a smaller Higgs mass than for \(\varphi _u = 0\), when we might expect that when the A-terms are aligned with the \(\mu \)-terms we should have more mixing and thus larger masses.

Fig. 8
figure 8

Plots as the gluino phase \(\varphi _{M_3}\) is varied for \(M_3^0 = 2000\) \(\text {GeV}\); left plots have \(\varphi _u = 0\), right plots have \(\varphi _u = \pi \). The top plots show the lightest Higgs mass at one (blue, dashed) and two (red) loops; the middle plots show the next lightest (green) and heaviest (red) at one (dashed) and two (solid) loops. The bottom plots show the calculated variation of the phase \(\varphi _{B_\mu }\) at one (blue, dashed) and two (red) loops

Full size image

Both these observations have a simple explanation, which illustrates the need for the two-loop routines. In the points that we have chosen with near maximal mixing, the soft terms are nearly degenerate, and so changes in the trilinear terms make only a small difference to the mixing angles. If we take the tree-level values of the stop masses and mixing (as we are required to do) and naively take \(m_t = 173\) GeV, \(v=246\) GeV we find for \(\tan \beta = 5, M_3^0 = 2{\text {TeV}}\) the following calculated values for \(m_h^\mathrm{approx} \) from Eq. (25):

$$\begin{aligned} \begin{array}{||c|c|c|c|c|c||c||}\hline \varphi _u &{} m_{\tilde{t}_1}\ ({\text {GeV}})&{} m_{\tilde{t}_2}\ ({\text {GeV}}) &{} s_{2\theta }^2 &{} m_t\ ({\text {GeV}}) &{} v\ ({\text {GeV}}) &{} m_h^\mathrm{approx}\ ({\text {GeV}})\\ \hline 0 &{} 1403 &{}1716 &{} 0.936 &{} 173 &{} 246 &{} 123 \\ \pi &{} 1320 &{} 1781 &{} 0.970 &{} 173 &{} 246 &{} 129 \\ \hline \end{array}\nonumber \\ \end{aligned}$$
(32)

However, if we use the values actually calculated in SPheno we find

$$\begin{aligned} \begin{array}{||c|c|c|c|c|c|c||c||}\hline \varphi _u &{} \varphi _{M_3} &{} m_{\tilde{t}_1}\ ({\text {GeV}})&{} m_{\tilde{t}_2}\ ({\text {GeV}}) &{} s_{2\theta }^2 &{} m_t^{{\overline{\mathrm {DR}}}}\ ({\text {GeV}}) &{} v\ ({\text {GeV}}) &{} m_h^\mathrm{approx}\ ({\text {GeV}}) \\ \hline 0 &{} 0 &{}1403 &{}1716 &{} 0.936 &{} 151 &{} 244.2 &{} 101 \\ \pi &{} 0&{} 1320 &{} 1781 &{} 0.970 &{} 145 &{} 243.6 &{} 88.3 \\ 0 &{} \pi &{} 1403 &{}1716 &{} 0.936 &{} 147 &{} 244.6 &{} 95.6 \\ \pi &{} \pi &{} 1320 &{} 1781 &{} 0.970 &{} 152 &{} 243.0 &{} 100 \\ \hline \end{array} \end{aligned}$$
(33)

Clearly the one-loop variation in the mass when we change \(\varphi _{M_3}\) can only come from the change in the Yukawa coupling; the two-loop shifts then compensate for this (since the top–stop–gluino diagrams partly correspond to a self-energy correction to a top loop), which could presumably be more clearly seen in an on-shell scheme. We also see that when we vary \(\varphi _u\) the shift in the top Yukawa has a much larger effect than the change in the mixing angle (since this can only be small when the mixing is already large). Therefore this observation is particular to the large mixing case; if the mixing were smaller, then potentially the variation of \(\varphi _u\) could have the opposite effect and increase the Higgs mass as per our naive expectation.

The effect of \(\varphi _u, \varphi _{M_3}\) on the heavy Higgses is very similar to the light Higgs: the loops compensate for the variation of the top Yukawa. But we note the enormous variation in their masses between \(\varphi _u = 0, \pi \) and between one and two loops for \(\varphi _{u} = \pi \); this underlines the tuning involved in the on-shell scheme when maintaining a constant heavy Higgs/charged Higgs mass. Finally, the variation of the phase \(\varphi _{B_\mu }\) is significant enough so that, if we were fixing the phase of \(B_\mu \) and solving the tadpole equations for \(\eta \), for most of the parameter space there would be large electric dipole moments; instead we find throughout that \(|d_e/e| < 10^{-30}\).

4.4.2 Variation of trilinear phase

If we instead fix the gluino phase and vary \(\varphi _u\), then we obtain Figs. 9 (for \(\varphi _{M_3} = 0\)) and 10 (for \(\varphi _{M_3} = \pi /2\)). In those figures we show the variation of the light Higgs mass for three different values of \(M_3^0\): the absolute value of the gluino mass clearly has a significant effect on the shift in the Higgs mass, of up to 2 \(\text {GeV}\) variation in the difference of one- and two-loop results by itself.

Fig. 9
figure 9

Left plot of the lightest Higgs mass at one (blue) and two (red) loops as \(\varphi _u\) is varied. The bands show the variation between \(M_3 = 1500\ {\text {GeV}}\) and \(2500\ {\text {GeV}}\). Right difference between the one- and two-loop lightest Higgs masses for three different values of \(M_3\) marked in the plot, as \(\varphi _u\) is varied

Full size image
Fig. 10
figure 10

Same as Fig. 9 but with \(\varphi _{M_3} = \pi /2\)

Full size image

Overall we find that \(\varphi _u\) has a markedly larger impact on the Higgs mass than \(\varphi _{M_3}\); as we discussed in the previous subsection this is largely due to the shift in the top Yukawa coupling rather than the change in mixing of the stops. However, the effect of the two-loop corrections as we vary \(\varphi _u\) is clearly not to compensate for the shift in the top coupling: in contrast to the previous subsection we see large variations of \(\Delta m_h\) (this being the difference between two- and one-loop lightest Higgs masses) between \(-2\) and 6 \(\text {GeV}\) in the maximally CP-violating case of \(\varphi _{M_3} = \pi /2\); in that case we also see a significant asymmetry in the plot, which resembles in form Fig. 6 of Ref. [65] (indeed our parameter choices are deliberately similar), even though that calculation is in the on-shell scheme and as we have already remarked cannot be quantitatively compared. Even more than the previous subsection, this therefore underlines the importance of the two-loop corrections to obtaining a reliable calculation of the Higgs mass.

5 Two-loop CPV effects in the NMSSM beyond \(\mathcal {O}(\alpha _s\alpha _t)\)

Fig. 11
figure 11

Size of the two-loop corrections as a function of arg(\(A_t\)) for different values of \(|A_t|\): 1.5 TeV (dashed), 2.0 TeV (dotted), 2.5 TeV (full). The red lines include only corrections \(\mathcal {O}(\alpha _s \alpha _t)\), while the blue lines include all other corrections in the gaugeless limit approximation

Full size image
Fig. 12
figure 12

On the left size of the two-loop corrections as function of arg(\(\lambda \)) for different values of \(|\lambda |\): 0.4 (dashed), 0.6 (dotted), 0.7 (full). The red lines include only corrections \(\mathcal {O}(\alpha _s \alpha _t)\), while the blue lines include all other corrections in the gaugeless limit approximation. On the right: the differences between the red and blue lines

Full size image
Fig. 13
figure 13

On the left light singlet mass as a function of arg(\(\lambda \)) for the benchmark secenario TP3 of Ref. [94]. The right side, show the size of the two-loop corrections \(\mathcal {O}(\alpha _s\alpha _t)\) only (dashed, red line), and of the sum of all two-loop corrections (full, blue line). The black dashed line gives the mass at the one-loop level

Full size image

In Sect. 3, we concentrated on the impact of complex parameters on the one-loop corrections as well as the two-loop corrections \(\mathcal {O}(\alpha _s\alpha _t)\) in the complex NMSSM. However, with the combination SARAH/SPheno one can immediately go beyond this: all non-vanishing two-loop corrections in the gaugeless limit are included automatically. Therefore, we can check how well the entire impact of the complex phases is covered by the \(\mathcal {O}(\alpha _s\alpha _t)\) corrections. In the following, we use the same parameter values as in Eq. (21) if not stated otherwise.Footnote 4

We start with the phase of \(A_t\) and show the change in the SM-like Higgs mass as a function of arg(\(A_t\)) in Fig. 11 for three different values of \(|A_t|\). As expected, the overall difference between the full two-loop corrections and the approximation using \(\mathcal {O}(\alpha _s\alpha _t)\) grows with increasing \(|A_t|\). However, for given \(|A_t|\), the difference shows only a very mild dependence on the phase of \(A_t\). Thus, at least for this parameter point the main sensitivity of arg(\(A_t\)) appears in the \((\alpha _s\alpha _t)\) corrections.

This is different for other phases like the one of \(\lambda \): we show in Fig. 12 the SM-like Higgs mass as function of arg(\(\lambda \)) for three different values of \(|\lambda |\). Here, we see not only a visible shift between the two different two-loop calculations, but also the dependence on the phase is very different. While the \(\mathcal {O}(\alpha _s\alpha _t)\) corrections give the impression that the Higgs mass is reduced for a large phase of \(\lambda \), the full calculation shows exactly the opposite. The Higgs mass actually increases for this point with increasing arg(\(\lambda \)). As consequence, the Higgs mass is underestimated in the real case by the \(\mathcal {O}(\alpha _s\alpha _t)\) corrections by about 1.3–1.6 GeV for all three vales of \(|\lambda |\), while for arg(\(\lambda \)) = ±0.4 the discrepancy increases to 1.6–2.6 GeV. Thus, for singlet scenarios with large \(\lambda \) and CP violation, the additional corrections now available with SARAH/SPheno can alter the SM-like Higgs mass by \(O({\text {GeV}})\).

Moreover, leaving the discussion of the mass of the SM-like Higgs for a moment, we find even bigger effects of arg(\(\lambda \)) on the mass of a light singlet. This is depicted in Fig. 13 where we have taken as a basis the benchmark point TP3 of Ref. [94], but with a large phase arg(\(\lambda \))=0.40–0.43. In this range, the mass of the light singlet shows a large sensitivity to the phase of \(\lambda \). We find here that the two-loop corrections \(\mathcal {O}(\alpha _s\alpha _t)\) shift the singlet mass between \(-0.5\) and −1.0 GeV for this range, while the full two-loop corrections are about four times as large: for arg(\(\lambda )=-0.4\) they alter the singlet mass by −2 GeV, while for arg(\(\lambda )=-0.43\) they even cause a shift of more than −4 GeV.

Fig. 14
figure 14

Impact of arg(\(\kappa \)) for a scenario with two Higgs states close to 125 GeV as proposed in Ref. [56]. Top, left The mass of the second and third physical scalar at one loop (dashed), and including two-loop corrections \(\mathcal {O}(\alpha _s \alpha _t)\) (dotted), as well as all two-loop corrections as calculated by SARAH/SPheno (full line). Top, right difference between the scalars when using \(\mathcal {O}(\alpha _s \alpha _t)\) only (red), and the complete corrections (blue). Bottom, left the singlet fraction of second (red) and third (blue) scalar using two-loop corrections \(\mathcal {O}(\alpha _s \alpha _t)\) only (dashed), and the complete two-loop corrections in the gaugeless limit (full). Bottom, right the CP-even fraction of the two states. The colour code is the same as on the left

Full size image

We also briefly give an example for the effect of arg(\(\kappa \)) by picking the very last point proposed in Table 3 of Ref. [56]. The particular feature of this point is that it has two scalars close to the desired mass of 125 GeV when choosing a phase of \(\kappa \) of about 0.52. We show in Fig. 14 the sensitivity of the properties of these two scalars to changing this phase. First, we find that the actual phase at which the two states are closest in mass is only slightly different between the full two-loop calculation and the one including only corrections involving the strong interaction. However, the minimal difference in mass which we find for the full calculation is smaller than 1.0 GeV, while with the incomplete calculation it is not possible to come closer than 1.6 GeV when keeping all other parameters but the phase fixed. An even more important effect can be seen when considering the character of the two states: the singlet admixture of the doublet state is quite dependent on the used two-loop calculation. One finds for instance that the full calculation has a smaller mixing when moving away from the cross-over point than the \(\mathcal {O}(\alpha _s\alpha _t)\) predicts. Also close to the cross-over point we find that the mixing between the CP-even and -odd states is different between both calculations and the CP-odd component of the singlet would be underestimated by up to 10% when not including all necessary two-loop corrections.

6 Conclusion

We have presented the possibility of calculating the two-loop corrections to real scalar masses in SUSY models with CP violation using the public packages SARAH and SPheno. After summarizing the generic approach used in these calculations, we showed the self-consistency of all results and the perfect agreement with corrections implemented in NMSSMCALC for the scale-invariant NMSSM. We demonstrated for selected examples in the MSSM and NMSSM how interesting physical results can be obtained easily with the available functionality, and that the variations of the corrections with the CP-violating phases can be large. We discussed the different options for the complex MSSM to fix CP phases by the tadpole equations, and the bias which enters the calculation by doing that. In the case of the MSSM, the only equivalent calculations have been performed in the on-shell scheme, rendering the results difficult to compare precisely; this underlines the utility of these results in SARAH since the majority of spectrum generators deliberately use the \(\overline{\mathrm {DR}}\) scheme for applications to studying GUT models, gauge-mediation, etc. On the other hand, we do find pleasing qualitative agreement of our results.

Afterwards, for the complex NMSSM we have briefly analysed the effect of CP phases in the two-loop corrections beyond \(O(\alpha _s\alpha _t)\), which have not previously been available in any scheme. It was shown that the dependence on the phases of \(M_3\) and \(A_t\) is included to a large extent in the corrections involving the strong coupling. However, it turned out that for instance the effect of the phase of \(\lambda \) in the full two-loop corrections deviates clearly from the impression one has when considering only the \(\alpha _s\alpha _t\) corrections.

Finally, we want to stress again that it is now possible with the demonstrated approach to obtain the Higgs masses with very high accuracy not only for the CP-violating MSSM and NMSSM. A large variety of other SUSY models can now easily be studied in the presence of significant CP phases without the problem of having a large theoretical uncertainty in the mass predictions. We hope that this gives a new impetus to interesting phenomenological studies of SUSY models with CP violation.