The Higgs mass in the CP violating MSSM, NMSSM, and beyond

We discuss the automatised calculation of the Higgs mass in renormalisable supersymmetric models with complex parameters at the two-loop level. Our setup is based on the public codes SARAH and SPheno, which can now compute the two-loop corrections to masses of all neutral scalars in such theories. The generic ansatz for these calculations and the handling of the `Goldstone Boson catastrophe' is described. It is shown that we find perfect agreement with other existing two-loop calculations performed in the DR-bar scheme. We also use the functionality to derive results for the MSSM and NMSSM not available before: the Higgs mass in the constrained version of the complex MSSM, and the impact of CP phases in the two-loop corrections beyond order alpha-strong alpha-top for the scale invariant NMSSM are briefly analysed.

the Higgs mass at two loops may be unconstrained by flavour [66]. 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 1 , the Standard-Model-like ones [67], 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, in particular that of the electron de, which is constrained to be [68] |de| < 8.9 × 10 −29 e cm = 4.5 × 10 −15 e GeV −1 . (1) The typical value for electric or chromoelectric moments for fermions i of mass m i and a common SUSY scale M SUSY is [69] 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 µ-term we have de/e 5g 2 24 κe tan β sin(ϕµ) (3) and we therefore need a large suppression of the total angle ϕµ by roughly three orders of magnitude; however, if we just consider the sneutrino-chargino sector then we obtain in the same limit |de/e| which similiarly constrains more CP-violating phases.
Here the calculation is more complicated, since it depends on the electric dipole moments of the light squarks, their chromoelectric momentsd u,d (naively suppressed by a factor e 4π ), 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 be easily 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 [69]: |d Hg /e| 10 −28 cm × d u −d d 10 −25 cm 2 × 10 −28 cm. (6) If the first two generations of squarks are heavy then, again, this will not constrain the parameter space relevant for the Higgs mass. 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 ϕµ, η 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 [74][75][76][77][78][79] and SPheno [80,81]. This functionality extends the automatised two-loop calculations for the real case presented in Refs. [43,44]. 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(αsα t ). We explain in sec. 2 the underlying methodology used in the calculations and some technical subtleties of the new extension before we present in sec. 3 the validation of the routines in the presence of complex parameters. In secs. 4 and 5 we discuss some applications of these routines in the context of the MSSM and NMSSM, before we conlcude in sec. 6.

Methodology
The calculation of CP-violating corrections at two loops is now available in SARAH via the diagrammatic approach described in Ref. [44]. 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. [44] 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,82,83] arises either in the MSSM beyond the gaugeless limit, or in theories beyond the MSSM even in the gaugeless limit, in that the DR mass of goldstone bosons have indeterminate. The full on-shell mass of course being zero, the 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. [42,43]) 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 | gaugeless to solve the tadpole equations, and then use the parameters determined from these in our pure gaugeless tree-level potential V 0 | gaugeless to determine the masses in our theory. Since we are effectively working in a false minimum the DR goldstone masses entering in our two-loop calculation are non-zero and of order 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. [45], all have the form ∆ (2) V (given topology) = Real product of couplings × loop function of masses with the exception of one contribution involving fermions and scalars, given by However, this should be understood as Re( 1 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 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 for scalars S i , whereṼ D 0 is the gauge-coupling dependent part that vanishes in the gaugeless limit. At tree level, the tadpole equations are used to determine some subset of the DR mass parameters; let us take them here to be defined to be the diagonal terms: and we then compute the particle masses at tree level We then compute the potential and one-loop self-energies using these tree-level masses: From these we solve the tadpole corrections so that Here Π 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 Π 1,ij (0, M 2 0,ij ) = ∂ i ∂ j V 1 . For p 2 = 0 we can rewrite the above as 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 → S i + α i , then 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 ∆ 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Ṽ D removed, and we solve the tadpole equations according to where here we have denoted bym 2 ij ,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 requireV 1 . Therefore 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 1 2 m 2 ij S i S j +Ṽ 0 is invariant under the global symmetry (even if the minimum we choose is not) then we find 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 α i M 2 ij (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.

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 In red is the result of the new CP-violating code; the blue dashed curve is a quadratic interpolation of the values φu = 0, π 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 φu = 0, ±π.

Check of the Goldstone mass
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, added 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 his provides a sensitive check for the correctness of our results.

Comparison with the known NMSSM corrections
The next check is to compare with existing results in literature. For the MSSM with CP violation the codes FeynHiggs [84] and CPsuperH [85] exist. However, both 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. NMSSMCALC makes use of mixed DR-OS renormalization 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 allow 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 [41,64]. At the two-loop level only the O(α S α t ) corrections are included [65]. The missing two-loop corrections will lead inevitable to a difference between SPheno and NMSSMCALC. Moreover, as has been discussed in detail in Ref. [86] the determination of the running 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. All two-loop corrections but the ones O(α S α t ) were turned off.
3. The tadpole equations were modified to be solved for im(Aκ) and im(A λ ) instead of im(Tκ) and im(T λ ) at tree-level: NMSSMCALC solves the tree-level tadpole equations to calculate im(Aκ) and im(A λ ), but calculates the radiative shifts to im(Tκ) and im(T λ ), 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 solvs 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 κ. 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 becuase NMSSMCALC supports only real Yukawa couplings. -NMSSMCALC 2.0: 1. A flag to calculate only tree-level masses has been included.
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.
3. The finite shifts to g 1 , g 2 and v were put to zero to have a pure DR renormalisation.
4. We fixed a bug in the two-loop calculation which we found during our comparison. 2 The conventions for the phases of the Higgs fields are where η is used as input, and η S is calculated from the complex input of µ eff and λ via As default point we have chosen 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 function of im(λ), while in Fig. 4 the impact of the im(µ eff ), and im(κ) is shown. Finally, Fig. 5 depicts the dependence on im(M 3 ), im(A t ) and η. We find an overall very good agremment at the one-and two-loop level. The differences are at most O(10 MeV), which corresponds to the agreement obtained in Ref. [86] 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(αsα 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 F F S as discussed in sec. 2.

MSSM
In SARAH, we express the Higgs doublets in the MSSM in the usual form shown in eq. (18). Since we are considering CP-violation, the µ-term and holomorphic soft-breaking parameters are allowed to be complex, and in general all of the neutral scalars φ u,d , σ u,d mix, with one state yielding the Goldstone boson of the Z.
We use the pure 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 2 The expression for δ (2) M 2 H + , which is used to express the shifts in the two-loop scalar mass matrix, contained a wrong prefactor for δtσ d .  DR quantities (Yukawa couplings, gauge couplings, electroweak expectation value), we still have a choice of parameters to be eliminated by the scalar tadpole equations. If we define Bµ ≡ e iϕ Bµ |Bµ| these read The last two equations are not independent due to the gauge symmetries. Note also that η and ϕ Bµ are not independent: in the above they appear in the combination η + ϕ Bµ so at tree level η = −ϕ Bµ .
In the on shell scheme used by Refs. [37,62,63], the charged Higgs mass and µ are taken as input parameters; this is equivalent to specifying |Bµ| in the above equations and using the third tadpole equation to determine η + ϕ Bµ . 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 β 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 DR scheme η (hence η + ϕ Bµ ) 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 sec. 2. Instead, we fix η and use the tadpole equations to determine Bµ. Since the tadpole corrections are typically small compared to Bµvu this is a small adjustment and we can regard our choice of η as being, to a good approximation, equivalent to minus the phase of Bµ. For expediency due to the very strong constraints upon it we take η = 0; we must then regard this as a tuning between Bµ and the other CP-violating phases in the theory for large values of ϕ Bµ . There then remain two options: one conventional choice, as in the CP-conserving case, is to solve the tadpole equations for |µ| 2 and Bµ; 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 subsection 4.3). For our study with low-energy boundary conditions in subsection 4.4 we shall choose to solve for m 2 H d , m 2 Hu and Im(Bµ), taking µ and Re(Bµ) as input parameters. This has the advantage that the tree-level heavy-Higgs masses are simply fixed by |Bµ| sin β cos β 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 comparision of calculations in the two approaches impossible. In fact, this is further exacerbated since Refs. [37,63] 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 their two.

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 where mt 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 θ is defined as Note that the stop mixing includes an additional phase (compared to the CP-even case) corresponding to the phase of e iη T 3,3 u s β − y t µ * c β , but we do not need that here. The important observation here is that if we define µ ≡ e iϕµ |µ|, T 3,3 u ≡ e iϕu |T 3,3 u | then the phases only enter through the combination cos(η + ϕu + ϕµ), 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 σu, σ d fields. We find that the one-loop stop contribution to these tadpoles neglecting the gauge couplings is where we define A 0 (x) ≡ −x(log x/Q 2 − 1), Q being the renormalisation scale. The function of masses on the right is a slowly-varying function with typical value of order unity, so with η = 0 we have Im(Bµ) ∼ 3 16π 2 |y t µT 3,3 u | sin β sin(ϕu + ϕµ).
So for tan β = 5 (for example) and |T 3,3 u | = |µ| = 2000 GeV, y t ∼ 0.9 and maximal CP-violation in the combination of phases on the right hand side we have Im(Bµ) ∼ (270 GeV) 2 . For a purely imaginary Bµ this would correspond to tree-level charged/heavy Higgs mass of 600 GeV; hence charged Higgs masses below 600 GeV -or, more realistically, somewhat heavier -invoke additional fine-tuning and are difficult to impose in the 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. [37,63], which involve lighter charged Higgs masses and larger trilinear couplings than we have quoted here.

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 singluarities 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 M 2 h and charged Higgs mass matrix M 2 H ± in the gaugeless limit in the basis φ d , φu, σ d , σu as where 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.

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 |µ| 2 and Bµ (both real and imaginary parts) at low energies (fixing the phase of µ as a choice). Then due to the strong constraints from electric dipole moments the phases of η, µ, m 12 are constrained to be very small (of order 10 −3 ÷ 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µ and thus heavy additional Higgses, showing little CP-violating effects. We illustrate this in Fig. 6, where we define A 0 ≡ −|A 0 |e iϕ A and vary the phase for two points: Here we have given the values that we compute at two loops for the three Higgs neutral Higgs scalars with ϕ A = 0. In the figure we show the mass of the lightest Higgs at one and two loops as we vary ϕ 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. We selected small tan β values to maximise the visibility of the CP-violating effects, because in that way we can have a large difference between ϕ A = 0 and ϕ A = π. 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 ϕ A . However, this means that for typical points of the parameter space Bµ and µ 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 ϕ A = 0 and ϕ A = π; almost all of the variation shown in the both 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 de; we show the variation of this in Fig. 7 and note that a small region near ϕ A = ± π 2 is already excluded for our tan β = 10 scenario (recall that the contributions are proportional to tan β 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.

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 3,3 u 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 β scenario; we take as our parameter values 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 works such as [37,63]. The choice of the real part of Bµ (here we solve the tadpole equations for m 2 Hu , m 2 H d , Im(Bµ)) is such that at tree level the heavy Higgs masses m h2,3 are at one TeV.

Variation of gluino phase
We expect that the gaugino phase ϕ M3 , only entering the Higgs mass calculation at two loops in the 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 ϕ Bµ , for ϕu = 0 and ϕu = π. The first observation is that the difference between one and two loops is strongly dependent on ϕ M3 ; for ϕu = 0 this changes between 4 and 7 GeV, and for ϕu = π between 2 and −5 GeV -an overall shift of 7 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 ϕu = π points have a smaller Higgs mass than for ϕu = 0, when we might expect that when the A-terms are aligned with the µ-terms we should have more mixing and thus larger masses.
Both these observations have a simple explanation, that 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 β = 5, M 0 3 = 2TeV the following calculated values for m approx h from eq. (24): ϕu mt However, if we use the values actually calculated in SPheno we find: Clearly the one-loop variation in the mass when we change ϕ M3 can only come from the change in the Yukawa coupling; the two loop shifts are 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.
What we also see is that when we vary ϕ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 ϕu could have the opposite effect and increase the Higgs mass as per our naive expectation. The effect of ϕu, ϕ M3 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 ϕu = 0, π, and between one and two loops for ϕu = π; 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 ϕ Bµ is significant enough that, if we were fixing the phase of Bµ and solving the tadpole equations for η, for most of the parameter space there would be large electric dipole moments; instead we find throughout that |de/e| < 10 −30 .
In those figures we show the variation of the light Higgs mass for three different values of M 0 3 : the absolute value of the gluino mass clearly has a significant effect on the shift in the Higgs mass, of up to 2 GeVvariation in the difference of one-and two-loop results by itself. Overall we find that ϕu has a markedly larger impact on the Higgs mass than ϕ M3 ; 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 ϕu is clearly not to compensate for the shift in the top coupling: in contrast to the previous subsection we see large variations of ∆m h (being the difference between two-and one-loop lightest Higgs masses) between −2 and 6 GeV in the maximally CPviolating case of ϕ M3 = π/2; in that case we also see a significant asymmetry in the plot, which resembles in form figure 6 of Ref. [37] (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 quantitively compared. Even more than the previous subsection, this therefore underlines the importance of the two loop corrections to obtaining a realiable calculation of the Higgs mass.  In sec. 3, we concentrated on the impact of complex parameters on the one-loop corrections as well as the two-loop corrections O(α S α 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 O(α S α t ) corrections. In the following, we use the same parameter values as in eq. (20) if not stated otherwise.
We start with the phase of A t and show the change in the SM-like Higgs mass as 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 O(α S α t ) grows with increasing |A t |. However, for given |A t |, the difference shows only a very mild depdence on the phase of A t . Thus, at least for this parameter point the main sensitive of arg(A t ) appears in the (α S α t ) corrections. This is different for other phases like the one of λ: we show in Fig. 12 the SM-like Higgs mass as function of arg(λ) for three different values of |λ|. 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 O(α S α t ) corrections give the impression that the Higgs mass is reduced for a large phase of λ, the full calculation shows exactly the opposite. The Higgs mass actually increases for this point with increasing arg(λ). As consequence, the Higgs mass is underestimated in the real case by the O(α S α t ) corrections by about 1.3-1.6 GeV for all three vales of |λ|, while for arg(λ) = ± 0.4 the discrepancy increases to 1.6-2.6 GeV. Thus, for singlet scenarios with large λ and CP violation, the additional corrections now available with SARAH/SPheno can alter the SM-like Higgs mass by O(GeV). Moreover, leaving the discussion of the mass of the SM-like Higgs for a moment, we find even bigger effects of arg(λ) on the mass of a light singlet. This is depicted in Fig. 13 where we have taken as basis the benchmark point TP3 of Ref. [86], but with a large phase arg(λ)=0.40-0.43. In this range, the mass of the light singlet shows a large sensitivity to the phase of λ. We find here that the two-loop corrections O(α S α 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(λ)=-0.4 they alter the singlet mass by -2 GeV, while for arg(λ)=-0.43 they even cause a shift of more than -4 GeV. We also briefly give an example for the effect of arg(κ) by picking the very last point proposed in Tab. 3 of Ref. [57]. 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 κ 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 O(α S α 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.

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 summarising 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 DR scheme for applications to studying GUT models, gauge-mediation, etc. On the other hand, we do find a 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(αsα 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 λ in the full two-loop corrections deviates clearly from the impression one has when considering only the αsα 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 be easily 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.