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

Mark D. Goodsell, Florian Staub

► To cite this version:

Mark D. Goodsell, Florian Staub. The Higgs mass in the CP violating MSSM, NMSSM and beyond. European Physical Journal C: Particles and Fields, 2017, 77 (1), pp.46. 10.1140/epjc/s10052-016-4495-9. hal-01449811

HAL Id: hal-01449811

https://hal.sorbonne-universite.fr/hal-01449811

Submitted on 30 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Mark D. Goodsell¹,²,ᵃ, Florian Staub³

¹ UPMC Univ Paris 06, UMR 7589, Sorbonne Universités, LPTHE, 75005 Paris, France

² CNRS, UMR 7589, LPTHE, 75005 Paris, France

³ CERN, Theory Department, 1211 Geneva 23, Switzerland

Received: 18 May 2016 / Accepted: 31 October 2016

© The Author(s) 2017. This article is published with open access at Springerlink.com

**Abstract** 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 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 $O(\alpha_s\alpha_t)$ for the scale-invariant NMSSM are briefly analysed.

## 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–37]. In addition, for the next-to-minimal supersymmetric standard model (NMSSM) most calculations of the Higgs mass considered the CP conserving case [38–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–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–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–61], while diagrammatic calculations at the one-loop level [62] and certain two-loop level contributions [63–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_s^0$, $B_d^0$ 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

<footnote>ᵃ e-mail: goodsell@lpthe.jussieu.fr</footnote>
allowed to be of the same order as, if somewhat less than,<fnref n="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 $d_e$, which is constrained to be [72]

$$|d_e| < 8.9 \times 10^{-29} \text{ e cm} = 4.5 \times 10^{-15} \text{ e GeV}^{-1}. \quad (1)$$

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

$$\kappa_i \equiv \frac{m_i}{16\pi^2 M_{\text{SUSY}}^2} = 1.3 \times 10^{-25} \text{ cm} \times \frac{m_i}{\text{MeV}} \times \left( \frac{\text{TeV}}{M_{\text{SUSY}}} \right)^2 \quad (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

$$d_e/e \approx \frac{5g^2}{24} \kappa_e \tan \beta \sin(\varphi_\mu) \quad (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

$$|d_e/e| \approx \frac{4g^2}{3} \kappa_e \tan \beta |\sin(\varphi_\mu + \varphi_{M_2} + \eta)|, \quad (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–76])

$$|d_n/e| \lesssim 3.0 \times 10^{-26} \text{ cm}. \quad (5)$$

Here the calculation is more complicated, since it depends on the electric dipole moments of the light squarks, their chromoelectric moments $\tilde{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]:

$$|d_{\text{Hg}}/e| \approx 10^{-29} \text{ cm} \times \left| \frac{\tilde{d}_u - \tilde{d}_d}{10^{-26} \text{ cm}} \right| \lesssim 7.4 \times 10^{-30} \text{ cm}. \quad (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{d_W}{6} f_{abc} \epsilon^{\mu\nu\alpha\beta} G_{\alpha\beta}^a G_{\mu\rho}^b G_V^{c\rho}$ 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).<fnref n="2" /> The contribution from the latter to the neutron EDM is [71]

$$|d_n/e| \sim 4 \times 10^{-26} \text{ cm} \times [d_W \times 10^{10} \times (\text{GeV})^2]. \quad (7)$$

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

$$\begin{aligned} \tilde{d}_q^{\text{Barr-Zee}} &\sim \frac{m_q g_s \alpha_s}{(4\pi)^3 v^2} Z_{1h} Z_{1A} \sim 10^{-26} \text{ cm} \times Z_{1h} Z_{1A} \\ d_W &\sim \frac{g_s \alpha_s}{(4\pi)^3 v^2} Z_{1h} Z_{1A} \sim 10^{-9} (\text{GeV})^{-2} \times Z_{1h} Z_{1A}, \end{aligned} \quad (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

<footnote n="1">¹ The $0^-$ hypothesis for the Higgs boson is excluded, but couplings up to 0.83 times the Standard-Model values are still allowed for parity-violating couplings to Z-bosons of the form $SZ_{\mu\nu}\tilde{Z}^{\mu\nu}$ where S is a singlet, $Z_{\mu\nu}$ the Z field strength and $\tilde{Z}^{\mu\nu}$ the dual field strength.</footnote>

<footnote n="2">² There are also diagrams for both involving squarks and gluinos, in particular stops are important for the Weinberg operator. However, these are less significant once we take gluino/first generation squark masses above the current bounds.</footnote>
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 $\phi_\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 trilinear 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–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{DR}$ mass-squared 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{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|_{\text{gaugeless}}$ to solve the tadpole equations, and then use the parameters determined from these in our pure gaugeless tree-level potential $V_0|_{\text{gaugeless}}$ to determine the masses in our theory. Since we are effectively working in a false minimum the $\overline{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

$$ \Delta^{(2)} V_{\text{(given topology)}} = \left( \begin{array}{c} \text{Real product of couplings} \\ \times \left( \text{loop function of masses} \right) \end{array} \right) $$

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

$$ V_{\overline{F}\overline{F}S}^{(2)} = \frac{1}{4} y^{IJK} y^{I'J'K} M_{II'}^* M_{JJ'}^* f_{\overline{F}\overline{F}S}(m_I^2, m_J^2, m_k^2) + \text{c.c.} $$

The function $f_{\overline{F}\overline{F}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 $\text{Re}(\frac{1}{2}y^{IJK}y^{I'J'K}M_{II'}^*M_{JJ'}^*)f_{\overline{F}\overline{F}S}(m_I^2, m_J^2, m_k^2)$ 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

$$V_0 = \frac{1}{2} m_{ij}^2 S_i S_j + \tilde{V}_0 + \tilde{V}_0^D \quad (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{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 \\ \mathcal{M}_{0,ij}^2|_{\text{gaugeless}} &\equiv m_{0,ij}^2 + \partial_i \partial_j \tilde{V}_0. \end{aligned} \quad (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:

$$\Delta V = V_1(\mathcal{M}_{0,ij}^2) + V_2(\mathcal{M}_{0,ij}^2|_{\text{gaugeless}}). \quad (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, \\ \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) \\ &\quad + \partial_i \partial_j V_2(\mathcal{M}_{0,ij}^2|_{\text{gaugeless}}). \end{aligned} \quad (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) =$

For $p^2 = 0$ we can rewrite the above as

$$ \begin{aligned} \mathcal{M}_{ij}^2(0) &= \mathcal{M}_{0,ij}^2 + \Delta_1 \mathcal{M}_{ij}^2(\mathcal{M}_{0,ij}^2) \\ &\quad + \Delta_2 \mathcal{M}_{ij}^2(\mathcal{M}_{0,ij}^2|_{\text{gaugeless}}) \\ \Delta_\ell \mathcal{M}_{ij}^2 &\equiv -\frac{\delta_{ij}}{v_i} \partial_i V_\ell + \partial_i \partial_j V_\ell. \end{aligned} \quad (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

$$\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. \quad (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

$$\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}_0^D}{\partial S_i} \quad (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

$$\alpha_i \frac{\partial^2}{\partial S_i \partial S_j} \left( \frac{1}{2} \hat{m}_{ij}^2 S_i S_j + \tilde{V}_0 + \hat{V}_1 + V_2 \right) = - \frac{\partial \alpha_i}{\partial S_j} \frac{\partial \tilde{V}_0^D}{\partial S_i}. \quad (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 computation 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} \left( \frac{1}{2} m_{0,ij}^2 S_i S_j + \tilde{V}_0 \right) &= - \frac{\partial \alpha_i}{\partial S_j} \frac{\partial \tilde{V}_0^D}{\partial S_i}, \\ \hat{m}_{ij}^2 &= m_{0,ij}^2 - \frac{1}{v_i} \delta_{ij} (\partial_i \hat{V}_1 + \partial_j \hat{V}_2) \Big|_{S_i=v_i} \\ &\rightarrow \alpha_i (\Delta_2 \mathcal{M}_{ij}^2) = 0. \end{aligned} \quad (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 $\phi_u \equiv \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$ TeV. On the same plot we show an interpolation between the values in the CP-preserving MSSM for the values $\phi_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

$$m_h^{\text{interpolation}} = m_h^{\text{1 loop}}(\phi_u) + m_h^{\text{2 loops}}|_{\phi_u=0} - m_h^{\text{1 loop}}|_{\phi_u=0} + (\phi_u/\pi)^2 \left[ \begin{aligned} & m_h^{\text{2 loops}}|_{\phi_u=\pi} \\ & - m_h^{\text{1 loop}}|_{\phi_u=\pi} - m_h^{\text{2 loops}}|_{\phi_u=0} + m_h^{\text{1 loop}}|_{\phi_u=0} \end{aligned} \right]. \quad (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.

<figure><img src="image_3.png" /><figcaption>Fig. 1 Plot of the lightest Higgs boson mass in the MSSM, for the parameters given in Sect. 4.4 with $M_3 = 2000$ 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$</figcaption></figure>

## 3.2 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, 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,

<figure><img src="image_4.png" /><figcaption>Fig. 2 The calculated mass for the neutral Goldstone boson in the MSSM in Landau gauge as a function of $\text{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$</figcaption></figure>
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 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 $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 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(\alpha_s\alpha_t)$ were turned off.

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

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 it 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 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 $κ$.

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. 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 from Standard Model inputs given in appendix A of [88].

3. The finite shifts to $g_1, g_2$ and $v$ were put to zero.

4. We fixed a bug in the two-loop calculation which we found during our comparison.<fnref n="3" />

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( \frac{1}{\sqrt{2}} (v_d + \phi_d + i\sigma_d) \atop H_d^- \right), \\
H_u &\equiv e^{i\eta} \left( \frac{H_u^+}{\sqrt{2}} (v_u + \phi_u + i\sigma_u) \atop H_u^- \right), \\
S &\equiv e^{i\eta_S} \frac{1}{\sqrt{2}} (v_S + \phi_S + i\sigma_S),
\end{aligned} 
\quad (19) $$

where $\eta$ is used as input, and $\eta_S$ is calculated from the complex input of $\mu_{eff}$ and $\lambda$ via

$$ \eta_S = \arg(\mu_{eff}) - \arg(\lambda). \quad (20) $$

As default point we have chosen

$$ 
\begin{aligned}
\lambda &= 0.6, \quad \kappa = -0.3, \quad A_\kappa = 200 \text{ GeV}, \\
A_\lambda &= 1000 \text{ GeV}, \quad \tan\beta = 2.5, \quad \mu_{eff} = 250 \text{ GeV}, \\
M_1 &= M_2 = 1000 \text{ GeV}, \quad M_3 = 1500 \text{ GeV}, \\
A_t &= 1500 \text{ GeV}, \quad m_{\tilde{t}_L} = m_{\tilde{t}_R} = 1000 \text{ GeV}.
\end{aligned} 
\quad (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_{eff})$ and $im(\kappa)$ is shown. Finally, Fig. 5 depicts the dependence on $im(M_3)$, $im(A_t)$ and $\eta$.

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_{FFS}^{(2)}$ as discussed in Sect. 2.

Finally, while this section is intended to discuss validation of the code (so we chose parameter ranges that may

<footnote n="3">³ The expression for $\delta^{(2)}M_{H^+}^2$, which is used to express the shifts in the two-loop scalar mass matrix, contained a wrong prefactor for $\delta t_{\sigma_d}$, but this has now been corrected.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 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</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 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_{eff})$ (left) and $im(\kappa)$ (right). The blue lines correspond to SPheno, the red ones to NMSSMCALC</figcaption></figure>

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 $Im(M_3)$ is unconstrained, and we also found no constraint upon $Im(\kappa)$ for the values considered. However, as discussed in the introduction around Eq. (4), $|\eta|$ and $|Im(\mu_{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 $|Im(A_t)| \lesssim 1600$ GeV from the neutron EDM, and $|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 $\bar{D}\bar{R}$ 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 cor-
<figure><img src="image_7.png" /><figcaption>Fig. 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</figcaption></figure>

responding $\overline{\text{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_{\mu} \cdot H_{d}$ and $B_{\mu} \equiv e^{i\varphi_{B_{\mu}}} |B_{\mu}|$ these read

$$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}}) + \frac{1}{v_d}\frac{\partial \Delta V}{\partial \phi_d}, \qquad (22)$$

$$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}}) + \frac{1}{v_u}\frac{\partial \Delta V}{\partial \phi_u}, \qquad (23)$$

$$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}| + \frac{1}{v_d}\frac{\partial \Delta V}{\partial \sigma_u}. \qquad (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{\text{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 $\text{Im}(B_{\mu})$, taking $\mu$ and $\text{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

$$ (m_h^{\text{approx}})^2 = M_Z^2 c_{2\beta}^2 + \frac{3}{32\pi^2 v^2} \left[ \begin{aligned} & s_{2\theta}^2 (m_{\tilde{t}_1}^2 - m_{\tilde{t}_2}^2)^2 \\ & + 8m_t^2 \log \frac{m_{\tilde{t}_1}^2 m_{\tilde{t}_2}^2}{m_t^4} + \frac{s_{2\theta}^2}{2} (m_{\tilde{t}_1}^2 - m_{\tilde{t}_2}^2) \\ & \times \left( 8m_t^2 - s_{2\theta}^2 (m_{\tilde{t}_1}^2 + m_{\tilde{t}_2}^2) \right) \log \frac{m_{\tilde{t}_1}^2}{m_{\tilde{t}_2}^2} \end{aligned} \right] \quad (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

$$ s_{2\theta}^2 = \frac{2v^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}. \quad (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} \operatorname{Im} \left[ T_u^{3,3} e^{i\eta} (e^{i\eta} T_u^{3,3} s_\beta - y_t \mu^* c_\beta)^* \right] \\ & \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) \\ \frac{\partial \Delta V}{\partial \sigma_d} \supset & -\frac{3v}{16\pi^2} \operatorname{Im} \left[ \mu^* y_t (e^{i\eta} T_u^{3,3} s_\beta - y_t \mu^* c_\beta)^* \right] \\ & \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} \quad (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

$$ \operatorname{Im}(B_\mu) \sim \frac{3}{16\pi^2} |y_t \mu T_u^{3,3}| \sin \beta \sin(\varphi_u + \varphi_\mu). \quad (28) $$

So for $\tan \beta = 5$ (for example) and $|T_u^{3,3}| = |\mu| = 2000$ GeV, $y_t \sim 0.9$ and maximal CP violation in the combination of phases on the right-hand side we have $\operatorname{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 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. [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 $M_h^2$ and charged Higgs mass matrix $M_{H^\pm}^2$ in the gaugeless limit in the basis $\phi_d, \phi_u, \sigma_d, \sigma_u$ as

$$ \begin{aligned} M_h^2 &= \begin{pmatrix} -M_2^2(-\tan \beta) & 0_{2 \times 2} \\ 0_{2 \times 2} & M_2^2(\tan \beta) \end{pmatrix}, \\ M_{H^\pm}^2 &= M_2^2(\tan \beta), \end{aligned} \quad (29) $$

where

$$ M_2^2(x) = \begin{pmatrix} |B_\mu| x & |B_\mu| \\ |B_\mu| & |B_\mu|/x \end{pmatrix}. \quad (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.
<figure><img src="image_8.png" /><figcaption>Fig. 6 Plots of Higgs masses as the trilinear phase $\phi_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</figcaption></figure>

<figure><img src="image_9.png" /></figure>

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$

## 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\phi_A}$ and vary the phase for two points:

<table><thead><tr><th>tan &beta;</th><th>|A<sub>0</sub>|</th><th>m<sub>0</sub></th><th>M<sub>1/2</sub></th><th>m<sub>h<sub>1</sub></sub>|&phi;<sub>A</sub>=0</th><th>m<sub>h<sub>2</sub></sub>|&phi;<sub>A</sub>=0</th><th>m<sub>h<sub>3</sub></sub>|&phi;<sub>A</sub>=0</th></tr><tr><th></th><th>(GeV)</th><th>(GeV)</th><th>(GeV)</th><th>(GeV)</th><th>(GeV)</th><th>(GeV)</th></tr></thead><tbody><tr><td>5</td><td>4500</td><td>2000</td><td>2000</td><td>125.0</td><td>4118</td><td>4118</td></tr><tr><td>10</td><td>3300</td><td>1500</td><td>1000</td><td>125.5</td><td>2494</td><td>2495</td></tr></tbody></table>

Here we have given the values that we compute at two loops for the three Higgs neutral Higgs scalars with $\phi_A = 0$. In the figure we show the mass of the lightest Higgs at one and two loops as we vary $\phi_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 \beta$ values to maximise the visibility of the CP-violating effects, because in that way we can have a large difference between $\phi_A = 0$ and $\phi_A = \pi$. This also requires large values of the trilinear couplings, and the

<figure><img src="image_10.png" /><figcaption>Fig. 7 Electric dipole moment $d_e$ divided by electric charge $e$ as the trilinear phase $\phi_A$ is varied, with exclusion bands shown, for values of $\tan \beta = 5$ (red) and 10 (blue)</figcaption></figure>

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 $\phi_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 $\phi_A = 0$ and $\phi_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 trilinear 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 $\phi_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

$$ T_u^{3,3} = 0.9A_0e^{i\phi_u}, \quad T_d^{3,3} = 0.064A_0, \quad T_e^{3,3} = 0.05A_0, \\ A_0 = 3800 \text{ GeV}, \quad \mu = 3800 \text{ GeV}, \quad \tan \beta = 5, \\ \text{Re}(B_\mu) = 10^6 \text{ GeV}^2 \left( \frac{\tan \beta}{1 + \tan^2 \beta} \right) $$

$$ 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 \\ = (3 \times 10^3 \text{ GeV})^2, \quad i = 1, 2 $$

$$ m_{\tilde{L},33}^2 = m_{\tilde{e},33}^2 = m_{\tilde{d},33}^2 = (10^3 \text{ GeV})^2, \\ 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 $$

$$ M_1 = 200 \text{ GeV}, \quad M_2 = 500 \text{ GeV}, \quad M_3 = M_3^0 e^{i\phi_{M_3}}. \quad (31) $$

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 $\phi_u = \pi$ points have a smaller Higgs mass than for $\phi_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.

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 \text{ GeV}$, $v = 246 \text{ GeV}$ we find for $\tan \beta = 5$, $M_3^0 = 2 \text{TeV}$ the following calculated values for $m_h^{\text{approx}}$ from Eq. (25):

<table><thead><tr><th>φ<sub>u</sub></th><th>m<sub>τ<sub>1</sub></sub> (GeV)</th><th>m<sub>τ<sub>2</sub></sub> (GeV)</th><th>s<sup>2</sup><sub>2θ</sub></th><th>m<sub>t</sub> (GeV)</th><th>v (GeV)</th><th>m<sup>approx</sup><sub>h</sub> (GeV)</th></tr></thead><tbody><tr><td>0</td><td>1403</td><td>1716</td><td>0.936</td><td>173</td><td>246</td><td>123</td></tr><tr><td>π</td><td>1320</td><td>1781</td><td>0.970</td><td>173</td><td>246</td><td>129</td></tr></tbody></table>

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

<table><thead><tr><th>φ<sub>u</sub></th><th>φ<sub>M<sub>3</sub></sub></th><th>m<sub>τ<sub>1</sub></sub> (GeV)</th><th>m<sub>τ<sub>2</sub></sub> (GeV)</th><th>s<sup>2</sup><sub>2θ</sub></th><th>m<sup>DR</sup><sub>t</sub> (GeV)</th><th>v (GeV)</th><th>m<sup>approx</sup><sub>h</sub> (GeV)</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>1403</td><td>1716</td><td>0.936</td><td>151</td><td>244.2</td><td>101</td></tr><tr><td>π</td><td>0</td><td>1320</td><td>1781</td><td>0.970</td><td>145</td><td>243.6</td><td>88.3</td></tr><tr><td>0</td><td>π</td><td>1403</td><td>1716</td><td>0.936</td><td>147</td><td>244.6</td><td>95.6</td></tr><tr><td>π</td><td>π</td><td>1320</td><td>1781</td><td>0.970</td><td>152</td><td>243.0</td><td>100</td></tr></tbody></table>

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_{\tilde{H}_u}^2, m_{\tilde{H}_d}^2$, Im($B_\mu$)) is such that at tree level the heavy Higgs masses $m_{h_{2,3}}$ are at one TeV.

### 4.4.1 Variation of gluino phase

We expect that the gaugino phase $\phi_{M_3}$, 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 $\phi_{B_\mu}$, for $\phi_u = 0$ and $\phi_u = \pi$. The first observation is that the difference between one and two loops is strongly dependent on $\phi_{M_3}$; for $\phi_u = 0$ this changes between 4 and 7 GeV, and for $\phi_u = \pi$ between 2 and -5 GeV—an overall shift of 7 GeV in the latter case! While this point has been

Clearly the one-loop variation in the mass when we change $\phi_{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 $\phi_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 $\phi_u$ could have the opposite effect and increase the Higgs mass as per our naive expectation.

The effect of $\phi_u, \phi_{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 $\phi_u = 0, \pi$ and between one and two loops for $\phi_u = \pi$; this underlines the tuning involved in the on-shell
<figure><img src="image_11.png" /><figcaption>**Fig. 8** Plots as the gluino phase $\phi_{M_3}$ is varied for $M_3^0 = 2000$ GeV; left plots have $\phi_u = 0$, right plots have $\phi_u = \pi$. The top plots show the lightest Higgs mass at one (blue, dashed) and two (red) loops; the</figcaption></figure>

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 $\phi_{B_\mu}$ at one (blue, dashed) and two (red) loops

scheme when maintaining a constant heavy Higgs/charged Higgs mass. Finally, the variation of the phase $\phi_{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 $\phi_u$, then we obtain Figs. 9 (for $\phi_{M_3} = 0$) and 10 (for $\phi_{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 GeV variation in the difference of one- and two-loop results by itself.

Overall we find that $\phi_u$ has a markedly larger impact on the Higgs mass than $\phi_{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 $\phi_u$ is clearly not to compensate for the shift in the top coupling: in contrast
<figure><img src="image_12.png" /><figcaption>**Fig. 9** Left plot of the lightest Higgs mass at one (blue) and two (red) loops as $\phi_u$ is varied. The bands show the variation between $M_3 = 1500$ GeV and 2500 GeV. Right difference between the one- and two-loop lightest Higgs masses for three different values of $M_3$ marked in the plot, as $\phi_u$ is varied</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10** Same as Fig. 9 but with $\phi_{M_3} = \pi/2$</figcaption></figure>

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 GeV in the maximally CP-violating case of $\phi_{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 $O(\alpha_s\alpha_t)$

In Sect. 3, we concentrated on the impact of complex parameters on the one-loop corrections as well as the two-loop corrections $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 $O(\alpha_s\alpha_t)$ corrections. In the following, we

<figure><img src="image_14.png" /><figcaption>**Fig. 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 $O(\alpha_s\alpha_t)$, while the blue lines include all other corrections in the gaugeless limit approximation</figcaption></figure>

use the same parameter values as in Eq. (21) if not stated otherwise.<fnref n="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

<footnote n="4">⁴ Note also that therefore the discussion of CPV observables at the end of Sect. 3 applies.</footnote>
three different values of $|A_t|$. As expected, the overall differ-
ence between the full two-loop corrections and the approxi-
mation using $O(\alpha_s\alpha_t)$ grows with increasing $|A_t|$. However,
for given $|A_t|$, the difference shows only a very mild depen-
dence 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 λ: 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) cor-
rections 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/Pheno
can alter the SM-like Higgs mass by O(GeV).

<figure><img src="image_15.png" /><figcaption>**Fig. 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 $O(\alpha_s\alpha_t)$, while the *blue lines*</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13** On the *left* light singlet mass as a function of $\arg(\lambda)$ for the benchmark scenario TP3 of Ref. [94]. The *right side*, show the size of the two-loop corrections $O(\alpha_s\alpha_t)$ only (*dashed, red line*), and of the</figcaption></figure>

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
$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.

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 differ-

<figure><img src="image_17.png" /></figure>

include *all* other corrections in the gaugeless limit approximation. On
the *right*: the differences between the *red* and *blue lines*

<figure><img src="image_18.png" /></figure>

sum of all two-loop corrections (*full, blue line*). The *black dashed line* gives the mass at the one-loop level
<figure><img src="image_19.png" /><figcaption>**Fig. 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 $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 $O(\alpha_s\alpha_t)$ only (red), and the complete corrections</figcaption></figure>

(blue). Bottom, left the singlet fraction of second (red) and third (blue) scalar using two-loop corrections $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

ence 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(\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 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 cor-
rections 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.

**Acknowledgements** We thank Dao Thi Nhung for clarifying discussions concerning NMSSMCALC and Alexander Voigt for double checking one-loop results for the complex MSSM. MDG acknowledges support from Agence Nationale de Recherche Grant ANR-15-CE31-0002 "HiggsAutomator," and would like to thank Pietro Slavich for interesting discussions.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP³.

## References

1. CMS Collaboration, S. Chatrchyan et al., Phys. Lett. B 716, 30-61 (2012). arXiv:1207.7235

2. ATLAS Collaboration, G. Aad et al., Phys. Lett. B716, 1-29 (2012). arXiv:1207.7214

3. H.E. Haber, R. Hempfling, Phys. Rev. Lett. 66, 1815-1818 (1991)

4. J.R. Ellis, G. Ridolfi, F. Zwirner, Phys. Lett. B 257, 83-91 (1991)

5. Y. Okada, M. Yamaguchi, T. Yanagida, Prog. Theor. Phys. 85, 1-6 (1991)

6. Y. Okada, M. Yamaguchi, T. Yanagida, Phys. Lett. B 262, 54-58 (1991)

7. J.R. Ellis, G. Ridolfi, F. Zwirner, Phys. Lett. B 262, 477-484 (1991)

8. A. Brignole, Phys. Lett. B 281, 284-294 (1992)

9. P.H. Chankowski, S. Pokorski, J. Rosiek, Phys. Lett. B 274, 191-198 (1992)

10. A. Dabelstein, Z. Phys. C 67, 495-512 (1995). arXiv:hep-ph/9409375

11. D.M. Pierce, J.A. Bagger, K.T. Matchev, R.-J. Zhang, Nucl. Phys. B 491, 3-67 (1997). arXiv:hep-ph/9606211

12. R. Hempfling, A.H. Hoang, Phys. Lett. B 331, 99-106 (1994). arXiv:hep-ph/9401219

13. M.S. Carena, M. Quiros, C. Wagner, Nucl. Phys. B 461, 407-436 (1996). arXiv:hep-ph/9508343

14. S. Heinemeyer, W. Hollik, G. Weiglein, Phys. Rev. D 58, 091701 (1998). arXiv:hep-ph/9803277

15. R.-J. Zhang, Phys. Lett. B 447, 89-97 (1999). arXiv:hep-ph/9808299

16. S. Heinemeyer, W. Hollik, G. Weiglein, Eur. Phys. J. C 9, 343-366 (1999). arXiv:hep-ph/9812472

17. S. Heinemeyer, W. Hollik, G. Weiglein, Phys. Lett. B 455, 179-191 (1999). arXiv:hep-ph/9903404

18. J.R. Espinosa, R.-J. Zhang, JHEP 03, 026 (2000). arXiv:hep-ph/9912236

19. J.R. Espinosa, R.-J. Zhang, Nucl. Phys. B 586, 3-38 (2000). arXiv:hep-ph/0003246

20. A. Brignole, G. Degrassi, P. Slavich, F. Zwirner, Nucl. Phys. B 631, 195-218 (2002). arXiv:hep-ph/0112177

21. G. Degrassi, P. Slavich, F. Zwirner, Nucl. Phys. B 611, 403-422 (2001). arXiv:hep-ph/0105096

22. S.P. Martin, Phys. Rev. D 67, 095012 (2003). arXiv:hep-ph/0211366

23. A. Brignole, G. Degrassi, P. Slavich, F. Zwirner, Nucl. Phys. B 643, 79-92 (2002). arXiv:hep-ph/0206101

24. A. Dedes, P. Slavich, Nucl. Phys. B 657, 333-354 (2003). arXiv:hep-ph/0212132

25. S.P. Martin, Phys. Rev. D 66, 096001 (2002). arXiv:hep-ph/0206136

26. A. Dedes, G. Degrassi, P. Slavich, Nucl. Phys. B 672, 144-162 (2003). arXiv:hep-ph/0305127

27. S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Eur. Phys. J. C 39, 465-481 (2005). arXiv:hep-ph/0411114

28. K. Sasaki, M.S. Carena, C. Wagner, Nucl. Phys. B 381, 66-86 (1992)

29. M.S. Carena, J. Espinosa, M. Quiros, C. Wagner, Phys. Lett. B 355, 209-221 (1995). arXiv:hep-ph/9504316

30. H.E. Haber, R. Hempfling, A.H. Hoang, Z. Phys. C 75, 539-554 (1997). arXiv:hep-ph/9609331

31. J.R. Espinosa, I. Navarro, Nucl. Phys. B 615, 82-116 (2001). arXiv:hep-ph/0104047

32. S.P. Martin, Phys. Rev. D 75, 055005 (2007). arXiv:hep-ph/0701051

33. P. Kant, R. Harlander, L. Mihaila, M. Steinhauser, JHEP 1008, 104 (2010). arXiv:1005.5709

34. R. Harlander, P. Kant, L. Mihaila, M. Steinhauser, Phys. Rev. Lett. 100, 191602 (2008). arXiv:0803.0672

35. S.P. Martin, Phys. Rev. D 71, 016012 (2005). arXiv:hep-ph/0405022

36. S. Borowka, T. Hahn, S. Heinemeyer, G. Heinrich, W. Hollik, Eur. Phys. J. C 74(8), 2994 (2014). arXiv:1404.7074

37. G. Degrassi, S. Di Vita, P. Slavich, Eur. Phys. J. C 75(2), 61 (2015). arXiv:1410.3432

38. G. Degrassi, P. Slavich, Nucl. Phys. B 825, 119-150 (2010). arXiv:0907.4682

39. F. Staub, W. Porod, B. Herrmann, JHEP 1010, 040 (2010). arXiv:1007.4049

40. K. Ender, T. Graf, M. Muhleitner, H. Rzehak, Phys. Rev. D 85, 075024 (2012). arXiv:1111.4952

41. M.D. Goodsell, K. Nickel, F. Staub, Phys. Rev. D 91, 035021 (2015). arXiv:1411.4665

42. M.D. Goodsell, K. Nickel, F. Staub, Eur. Phys. J. C 75(1), 32 (2015). arXiv:1411.0675

43. M. Goodsell, K. Nickel, and F. Staub, Eur. Phys. J. C 75(6), 290 (2015). arXiv:1503.03098

44. S.P. Martin, Phys. Rev. D 65, 116003 (2002). arXiv:hep-ph/0111209

45. S.P. Martin, Phys. Rev. D 68, 075002 (2003). arXiv:hep-ph/0307101

46. S.P. Martin, Phys. Rev. D 70, 016005 (2004). arXiv:hep-ph/0312092

47. S.P. Martin, Phys. Rev. D 71, 116004 (2005). arXiv:hep-ph/0502168

48. M. Carena, J.R. Ellis, S. Mrenna, A. Pilaftsis, C.E.M. Wagner, Nucl. Phys. B 659, 145-178 (2003). arXiv:hep-ph/0211467

49. A. Bartl, H. Fraas, S. Hesselbach, K. Hohenwarter-Sodek, T. Kernerreiter, G. Moortgat-Pick, Eur. Phys. J. C 51, 149-161 (2007). arXiv:hep-ph/0608065

50. K.E. Williams, G. Weiglein, Phys. Lett. B 660, 217-227 (2008). arXiv:0710.5320

51. C.-W. Chiang, E. Senaha, JHEP 06, 019 (2008). arXiv:0804.1719
52. M. Terwort, O. Kittel, G. Moortgat-Pick, K. Rolbiecki, P. Schade, **1201**, 5272 (2012)

53. A. Bharucha, S. Heinemeyer, F. von der Pahlen, Eur. Phys. J. C **73**(11), 2629 (2013). arXiv:1307.4237

54. S. Munir, Phys. Rev. D **89**(9), 095013 (2014). arXiv:1310.8129

55. S. Heinemeyer, C. Schappacher, (2015). arXiv:1511.06002

56. S. Moretti, S. Munir, Adv. High Energy Phys. **2015**, 509847 (2015). arXiv:1505.00545

57. J. Berger, M. W. Cahill-Rowley, D. Ghosh, J. L. Hewett, A. Ismail, T.G. Rizzo, Phys. Rev. D **93**(3), 035017 (2016). arXiv:1510.08840

58. M. Carena, J. Ellis, J.S. Lee, A. Pilaftsis, C.E.M. Wagner, JHEP **02**, 123 (2016). arXiv:1512.00437

59. M.S. Carena, H. Haber, S. Heinemeyer, W. Hollik, C. Wagner et al., Nucl. Phys. B **580**, 29–57 (2000). arXiv:hep-ph/0001002

60. M.S. Carena, J.R. Ellis, A. Pilaftsis, C. Wagner, Nucl. Phys. B **586**, 92–140 (2000). arXiv:hep-ph/0003180

61. M.S. Carena, J.R. Ellis, A. Pilaftsis, C. Wagner, Nucl. Phys. B **625**, 345–371 (2002). arXiv:hep-ph/0111245

62. M. Frank, T. Hahn, S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, JHEP **02**, 047 (2007). arXiv:hep-ph/0611326

63. S. Heinemeyer, W. Hollik, H. Rzehak, G. Weiglein, Phys. Lett. B **652**, 300–309 (2007). arXiv:0705.0746

64. W. Hollik, S. Paßehr, Phys. Lett. B **733**, 144–150 (2014). arXiv:1401.8275

65. W. Hollik, S. Paßehr, JHEP **10**, 171 (2014). arXiv:1409.1687

66. T. Graf, R. Gröber, M. Mühlleitner, H. Rzehak, K. Walz, JHEP **1210**, 122 (2012). arXiv:1206.6806

67. M. Mühlleitner, D.T. Nhung, H. Rzehak, K. Walz, JHEP **05**, 128 (2015). arXiv:1412.0918

68. M.D. Goodsell, K. Nickel, F. Staub, Phys. Lett. B **758**, 18–25 (2016). arXiv:1511.01904

69. ATLAS collaboration, ATLAS-CONF-2015-008 (2015)

70. S. Abel, S. Khalil, O. Lebedev, Nucl. Phys. B **606**, 151–182 (2001). arXiv:hep-ph/0103320

71. M. Pospelov, A. Ritz, Ann. Phys. **318**, 119–169 (2005). arXiv:hep-ph/0504231

72. J. Acme, Baron et al. Science **343**, 269–272 (2014). arXiv:1310.7534

73. J. Pendlebury et al., Phys. Rev. D **92**(9), 092003 (2015). arXiv:1509.04411

74. C.A. Baker et al., Phys. Rev. Lett. **97**, 131801 (2006). arXiv:hep-ex/0602020

75. C.A. Baker et al., Phys. Rev. Lett. **98**, 149102 (2007). arXiv:0704.1354

76. A.P. Serebrov et al., Phys. Rev. C **92**(5), 055501 (2015)

77. B. Graner, Y. Chen, E.G. Lindahl, B.R. Heckel, Phys. Rev. Lett. **116**(16), 161601 (2016). arXiv:1601.04339

78. M. Jung, A. Pich, JHEP **04**, 076 (2014). arXiv:1308.6283

79. W. Dekens, J. de Vries, J. Bsaisou, W. Bernreuther, C. Hanhart, U.-G. MeiÄßner, A. Nogga, A. Wirzba, JHEP **07**, 069 (2014). arXiv:1404.6082

80. J. Baglio, R. Gröber, M. Mühlleitner, D.T. Nhung, H. Rzehak, M. Spira, J. Streicher, K. Walz, Comput. Phys. Commun. **185**(12), 3372–3391 (2014). arXiv:1312.4788

81. S.F. King, M. Mühlleitner, R. Nevzorov, K. Walz, Nucl. Phys. B **901**, 526–555 (2015). arXiv:1508.03255

82. F. Staub, (2008). arXiv:0806.0538

83. F. Staub, Comput. Phys. Commun. **181**, 1077–1086 (2010). arXiv:0909.2863

84. F. Staub, Comput. Phys. Commun. **182**, 808–833 (2011). arXiv:1002.0840

85. F. Staub, Comput. Phys. Commun. **184**, 1792–1809 (2013). arXiv:1207.0906

86. F. Staub, Comput. Phys. Commun. **185**, 1773–1790 (2014). arXiv:1309.7223

87. F. Staub, Adv. High Energy Phys. **2015**, 840780 (2015). arXiv:1503.04200

88. W. Porod, F. Staub, Comput. Phys. Commun. **183**, 2458–2469 (2012). arXiv:1104.1573

89. W. Porod, Comput. Phys. Commun. **153**, 275–315 (2003). arXiv:hep-ph/0301101

90. S.P. Martin, Phys. Rev. D **90**(1), 016013 (2014). arXiv:1406.2355

91. J. Elias-Miro, J. Espinosa, T. Konstandin, JHEP **1408**, 034 (2014). arXiv:1406.2652

92. S. Heinemeyer, W. Hollik, G. Weiglein, Comput. Phys. Commun. **124**, 76–89 (2000). arXiv:hep-ph/9812320

93. J. Lee, A. Pilaftsis, M.S. Carena, S. Choi, M. Drees et al., Comput. Phys. Commun. **156**, 283–317 (2004). arXiv:hep-ph/0307377

94. F. Staub, P. Athron, U. Ellwanger, R. Gröber, M. Mühlleitner, P. Slavich, A. Voigt, Comput. Phys. Commun. **202**, 113–130 (2016). arXiv:1507.05093

95. B.C. Allanach et al., Comput. Phys. Commun. **180**, 8–25 (2009). arXiv:0801.0045

96. P.Z. Skands et al., JHEP **07**, 036 (2004). arXiv:hep-ph/0311123