A New Tool for the study of the CP-violating NMSSM

Supersymmetric extensions of the Standard Model open up the possibility for new types of CP-violation. We consider the case of the Next-to-Minimal Supersymmetric Standard Model where, beyond the phases from the soft lagrangian, CP-violation could enter the Higgs sector directly at tree-level through complex parameters in the superpotential. We develop a series of Fortran subroutines, cast within the public tool NMSSMTools and allowing for a phenomenological analysis of the CP-violating NMSSM. This new tool performs the computation of the masses and couplings of the various new physics states in this model: leading corrections to the sparticle masses are included; the precision for the Higgs masses and couplings reaches the full one-loop and leading two-loop order. The two-body Higgs and top decays are also attended. We use the public tools HiggsBounds and HiggsSignals to test the Higgs sector. Additional subroutines check the viability of the sparticle spectrum in view of LEP-limits and constrain the phases of the model via a confrontation to the experimentally measured Electric Dipole Moments. These tools will be made publicly available in the near future. In this paper, we detail the workings of our code and illustrate its use via a comparison with existing results. We also consider some consequences of CP-violation for the NMSSM Higgs sector.


Introduction
After the discovery of a signal at a mass of about 125 GeV in the LHC Higgs searches [1,2], the question of the identification of the associated state(s) and the underlying physics remains open. While the general properties are consistent so far with those expected for the Higgs boson of the Standard Model (SM), a wide range of alternatives could equally well fit the experimental data. In particular, softly-broken supersymmetric (SUSY) extensions of the SM [3] count among the appealing options to solve the Hierarchy Problem [4] and allow for a smooth transition to higher energy physics (e.g. Grand Unification, neutrino physics or weakly-coupled dark matter). The scarceness of evidence for new physics effects in precision physics or direct searches should also be weighed by the considerations that SUSY-inspired models offer a SM-like decoupling regime, but also that complex mechanisms in the Higgs or SUSY sectors -see e.g. [5] -may account for this relative invisibility thus far.
The Next-to-Minimal Supersymmetric Standard Model (NMSSM), a singlet-extension of the simplest viable SUSY-inspired extension of the SM [6], has raised renewed interest ever since the Higgs discovery, notably due to its properties in the Higgs sector, e.g. allowing for an uplift of the mass of the SM-like Higgs related to F-terms or to the mixing of this state with a lighter singlet [7]. The original motivation for this singlet extension rests with the 'µ-problem' of the MSSM [8], which can be solved elegantly if this µ-term is generated dynamically, via a singlet vacuum expectation value (v.e.v.) [9]. Correspondingly, the Z 3 -conserving version of the NMSSM -allowing only cubic terms in the superpotential -is the most studied form of this model, while more general singlet couplings can be justified by higher-energy considerations -see e.g. [10,11]. Another usual feature in SUSY extensions of the SM is R-parity, which both constrains the possibility of baryon-number violation and provides a stable SUSY particle, hence a dark-matter candidate.
A troubling fact rests with the observation that several NMSSM parameters -especially in the Higgs sector -can take complex values, hence lead to CP-violation beyond that in the quark sector. On the one hand, CP-violation is known as a cosmological necessity for baryogenesis. On the other, it receives severe limits at the phenomenological level, from the non-observation of Electric Dipole Moments (EDM's; see e.g. [12]). In this paper, we aim at presenting a tool which allows to study the NMSSM with complex parameters within the framework of the public code NMSSMTools [13]. In a first step, we will focus on the Z 3 -conserving version although we plan on a generalization to tadpole and quadratic couplings of the singlet in the future.
The current version of NMSSMTools allows to perform several operations in connection with the spectrum of the CP-conserving NMSSM: in particular, it computes radiative corrections to the Higgs and SUSY spectrum, calculates the widths of Higgs decays or confronts the NMSSM parameter space to theoretical -e.g. vacuum stability -or phenomenological -e.g. Higgs searches, B-physics -limits. Several other tools aiming at the calculation of radiative corrections to the NMSSM spectrum have been developed in the past few years, e.g. NMSSMCALC [14] or SoftSUSY [15]. The latter focuses on the CP-conserving NMSSM, while NMSSMCALC allows for CP-violation but specializes in corrections to the Higgs spectrum. Other tools, less specialized in a given model, allow for similar manipulations, provided the implementation of a model-file as input: this applies to SPHENO [16] or FlexibleSUSY [17].
Our goal consists in generalizing NMSSMTools to the CP-violating case. While this task remains far from complete, the tools which we present here already allow for numerous operations: radiative corrections to the SUSY and the Higgs masses are implemented -so far, only the leading double-log corrections (beyond the full one-loop) are taken into account at two-loop order in the Higgs spectrum -; Higgs and top two-body decays are computed; phenomenological limits from LEP SUSY searches or Higgs physics are tested -the latter via an interface with the public tools HiggsBounds [18] and HiggsSignals [19] -; finally we designed a subroutine to estimate the EDM's. All these routines should become available on the NMSSMTools website [13] in the near future. This paper is intended to serve as a presentation of the calculations implemented in our tool, as well as a short illustration of its uses. In the following section, we will detail the characteristics of the model under study, the underlying assumptions and the tree-level spectrum. The third section will present the chain of subroutines that we designed and the operations which they carry out. Finally, we will consider phenomenological consequences and compare some of our results to the predictions of existing tools, before we conclude.

Model, Phase-counting and Tree-level
In this section, we present the details of the model under consideration, our notations as well as the spectrum at tree-level.

The CP-violating NMSSM
The NMSSM is a supersymmetry-inspired extension of the SM with soft SUSY-breaking terms. It differs from the minimal supersymmetric extension of the SM, the MSSM, in that it includes, in addition to the two Higgs SU (2) Ldoublet superfieldsĤ u andĤ d with opposite hypercharge ±1, a supplemental gauge-singlet chiral superfieldŜ. While the couplings of this singlet may take a more complex form in the general case, we will be considering only the R-parity and Z 3 -conserving NMSSM here, which is characterized by the following superpotential and SUSY-breaking terms: The 'matter' (super)fields 1 Q L , U c R , D c R , L L , E c R should be understood as summed over generations and the parameters within brackets should correspondingly be seen as (complex) matrices. '·' denotes the usual SU (2) L product.b,w α andg a stand for the U (1) Y , SU (2) L and SU (3) c gauginos, respectively. In the following g , g and g S will denote the corresponding gauge couplings and α S ≡ g 2 S /4π. While the Z 3 -conserving NMSSM offers the simplest solution to the µ-problem of the MSSM, the inclusion of Z 3 -violating terms can be justified from higher-energy considerations [10,11] and turns up as a phenomenological necessity in view of the domain-wall problem. Our restriction to the Z 3 -conserving lagrangian follows considerations of simplicity and our work shall be extended to the Z 3 -violating case in the near future.
The minimization of the scalar potential will generate Higgs vacuum expectation values (v.e.v.'s) so that we may write the Higgs (super)fields in terms of their (real and positive) v.e.v.'s s, v u , v d , and their charged and neutral components: 1 We will omit theˆdistinguishing the superfields from their scalar component, from now on.
The three 'dynamical' phases φ s , φ u and φ d add to the 'static' phases appearing in the lagrangian density (Eqs. 1,2). From now on, we will make the following replacements in our notations: The Yukawa matrices may be written in terms of (real and positive) matrices Y u , Y d , Y e , diagonal in flavour space, using unitary transformations: Redefining the quark and lepton (super)fields accordingly, and introducing the Cabibbo-Kobayashi-Maskawa (CKM) matrix V CKM ≡ X u † L X d L , the superpotential of Eq.1 now reads: Finally, we make the following assumptions to ensure minimal flavour violation in the sfermion sector: where m 2 Q is a diagonal (and, without loss of generality, real) matrix in flavour space. The approximation ' ' only holds for a matrix proportional to the identity, in the strict sense, but is viable, considering that the CKM matrix is hierarchical. Note that we will assume degeneracy for the first two generations of sfermions.
Eqs.7 and 8 fully characterize the model that we will be considering from now on -note that the three latter terms of Eq.7 as well as the second and fourth lines of Eq.8 are still implicitly summed over fermion generations. All the phases have been explicited and reduce, at this level, to four phases in the Higgs sector -ϕ λ , ϕ κ , ϕ 1 , ϕ 2 ; we will see that the minimization conditions further constrain these, as could be expected from the 'dynamical' nature of some phases -, three gaugino phases -φ M1 , φ M2 , φ M3 -, three sfermion phases per generation -ϕ Au , ϕ A d , ϕ Ae -and the CKM phase finally. Given that we will neglect the Yukawa couplings of the first two generations, only the sfermion phases of the third generation will intervene in practice.

The tree-level Higgs sector
The Higgs potential collects terms from the soft lagrangian (Eq.8), F-terms from the superpotential (Eq.7) and D-terms from the gauge interactions. We obtain: A κ e ıϕ2 S 3 + h.c.
The neutral part reduces to: At tree level, the Higgs v.e.v.'s are assumed to minimize this potential. A consequence is the cancellation of first derivatives with respect to the neutral Higgs fields at the minimum, which provides us with the minimization conditions: Here we see that the four phases of the Higgs sector are not independent but that, on the contrary, the minimization conditions relate ϕ 1 and ϕ 2 to ϕ λ − ϕ κ , the latter being the one and only 'observable' phase in the Higgs sector. Note that ϕ λ and ϕ κ intervene independently in other parts of the spectrum however. We will make an explicit use of the minimization conditions of Eq.11 in the following lines, replacing m 2 Hu , m 2 H d , m 2 S , A λ sin ϕ 1 and A κ sin ϕ 2 by their expressions in terms of the v.e.v.'s.
The terms of Eq.9, bilinear in the charged Higgs fields, define the 2 × 2 (hermitian) mass-matrix of the charged-Higgs states: which determines the charged Goldstone boson G ± = − sin β H ± u + cos β H ± d and the physical charged Higgs state H ± = cos β H ± u + sin β H ± d . Similarly, the terms bilinear in the neutral Higgs fields provide the 6×6 (symmetric) mass-matrix of the neutral 2 , with the notation meaning that fields are frozen to their v.e.v.'s. In As in the charged case, the neutral Goldstone boson can be singled out via a β-angle rotation G 0 ≡ − sin β a 0 u + cos β a 0 d . The remaining 5 × 5 symmetric block spanning the space (h 0 u , h 0 d , h 0 s , a 0 ≡ cos β a 0 u + sin β a 0 d , a 0 s ) may be diagonalized via an orthogonal matrix X H 0 : which defines the mass eigenstates: We will use the second notation which allows more clarity in the identification of the components. Additionally, we define X I iu ≡ cos β X I ia and X I id ≡ sin β X I ia . Note that the positivity of the squared Higgs-masses is a stability condition of the vacuum. Remember also that, at 0 th order in the electroweak v.e.v.'s, one can isolate the CP-even and CP-odd sectors and diagonalize their doublet subspaces via rotations of angle −β / β (the singlet states are then unmixed), which disentangles the 'light' (then fully massless) 'SM-like' doublet states from the 'heavy' states with approximate squared-mass vuv d (degenerate at this order with the charged state).

The supersymmetric spectrum at tree-level
The whole tree-level spectrum will be treated with further details in appendix B. Here we simply summarize, for the sake of notations, the basic ingredients concerning the treatment of masses and mixings of SUSY particles.

i) Gluinos
The gluinos are the fermionic partners of the gluons and, as such, form a color octet. Their bilinear terms originate in the soft lagrangian: −L soft −M 3 e ıφ M 3g aga . The mass statesG a , with mass M 3 (which we assume positive), then relate to the eigenstates of the SUSY vector superfieldg a asG a ≡ −ıe ı 2 φ M 3g a . The phase shift then affects the couplings of the gluinos to coloured matter.

ii) Charginos
The charginos are composed of the charged components of the electroweak gauginos and higgsinos. Their bilinear terms originate from both supersymmetry-conserving and violating terms and may be cast into the following form: We may diagonalize M χ −+ with the help of two unitary matrices U and V : )V . The mass eigenvalues may be assumed real and positive without any loss of generality and the mass eigenstates (i = 1, 2) relate to the gauge ones as: iii) Neutralinos The neutralinos are combinations of the neutral components of the electroweak gauginos and higgsinos. Their bilinear terms, resulting from both supersymmetry-conserving and violating terms, form a Majorana mass matrix: M χ 0 being symmetric, it can be diagonalized by a single unitary matrix N according to: . . , 5)N . Without loss of generality the eigenvalues m χ 0 i can be chosen real and positive (remember that N is complex) and the mass eigenstates relate to the gauge ones in the following fashion: iv) Sfermions The scalar partners of the SM fermions receive hermitian mass matrices. Due to our assumptions with respect to flavour violation, the three generations decouple. We keep a generic notation although only the Yukawa couplings of the third generation (u = t, d = b, e = τ ) will be treated as non-vanishing in practice: Each mass matrix M 2 F -F = U, D, N, E -can be diagonalized via a special-unitary matrix X F , according to: The positivity of the squared masses m 2 Fi is a stability condition of the vacuum. The mass eigenstates are then defined as: This completes this short presentation of the tree-level spectrum. More details are presented in appendix B, together with the Higgs couplings.

A short walk-through the code
In this section, we shall describe the operations which are conducted throughout our subroutines from the perspective of the phenomenology of the CP-violating NMSSM.

Interface with NMSSMTools
Before coming to the actual computations of our code, let us remind the reader that we embed it within the NMSSMTools package. We actually use the NMHDECAY routines to define its input. In particular, we do not alter the running of parameters -such as the Yukawa, gauge or soft couplings -, e.g. to the average scale of the squarks of third generation. We simply use the corresponding quantities as calculated by NMHDECAY as our input and introduce the complex phases at this level. This is sufficient in consideration of the order of precision which we aim at in the radiative corrections. A short subroutine init_CPV.f defines this interface and stores all the relevant quantities within commons of the code. The case of the parameters A λ and A κ is somewhat more subtle: given that the phases ϕ 1 and ϕ 2 are not free but, in our approach, determined by the minimization conditions of the potential (see Eq.11), we will only be using the quantities A λ cos ϕ 1 and A κ cos ϕ 2 as degrees of freedom in practice. Therefore, we identify the NMHDECAY input for A λ and A κ as ours for A λ cos ϕ 1 and A κ cos ϕ 2 . The wave-function scaling factors for the Higgs fields are also defined slightly differently from the original implementation in NMSSMTools, as we shall describe in section 3.3.1.
Given our discussion in section 2, the following eight phases are added as new degrees of freedom:

Supersymmetric spectrum
The first actual operations which are carried out in connection to the CP-violating NMSSM consist in the calculation of the masses of the supersymmetric matter content. Similarly to the evaluation by NMSSMTools in the CP-conserving case, we take into account the leading radiative corrections to the masses. In the following, we list the new subroutines and provide relevant information concerning the calculations which are performed.
i) mcha_CPV.f The purpose of this subroutine rests in diagonalizing the chargino mass matrix (Eq.16) according to Similarly to the corresponding implementation within NMSSMTools for the CP-conserving case, the entries of the mass matrix receive one-loop radiative corrections which are calculated in the approximation where mass and gauge eigenstates coincide. The corresponding effects are presented in section 4.2 of [20] -in the context the MSSM and still in the CP-conserving case. Small modifications appear in the CP-violating NMSSM, as gaugino and higgsino scalar couplings are rotated by phase factors of e −ıφ M i /2 and e −ıϕ λ /2 . Nevertheless, the factors of B 1 functions as well as the corrections involving gauge bosons are immune to this phase shift, so that only the scalar interactions resulting in a B 0 function -in the approximations of [20], this reduces to the Higgs / higgsino loops -are affected. Another difference with respect to ref. [20] originates from the presence of singlets and singlinos in the higgsino self-energies. A summary of these corrections is explicited in appendix C.1.
The following steps are essentially identical to their counterparts in the tree-level case, which is treated into details in appendix B.1.4: we define two special-unitary matrices U 0 and V 0 diagonalizing the hermitian matrices is then a diagonal matrix with, in general, non-real entries. We thus define the unitary matrices U and V via a phase-shift of U 0 and V 0 , where the phase of the lightest state is absorbed in U while that of the heavier one is absorbed in V : the resulting chargino masses are real and positive.
ii) mneu_CPV.f The case of the neutralinos follows the same principles as that of the charginos. The tree-level gaugino and higgsino masses are corrected in accordance with the one-loop effects presented in appendix C.1. We then diagonalize the complex symmetric neutralino mass matrix according to M χ 0 = N T diag(m χ 0 i , i = 1, . . . , 5)N . For that purpose, we consider the 10×10 real symmetric matrix Re Im −Im Re M χ 0 † M χ 0 , which can be diagonalized numerically by an orthogonal matrixÑ 0 . We then extract a special unitary matrix N 0 so that N * 0 M χ 0 N † 0 is diagonal. We finally absorb the remaining phases in a phase-shift of N 0 , which defines the real and positive neutralino masses as well as the mixing matrix N . Details are provided in appendix B.1.4.
iii) msferm_CPV.f We now turn to the sfermion masses. The hermitian tree-level mass matrices are diagonalized via special-unitary matrices X F , according to M 2 F = X F † diag(m 2 F1 , m 2 F2 )X F . We remind the reader that the parameters entering the matrices, e.g. the top and bottom Yukawa couplings, have been run to the average squark scale. The Yukawa couplings of the first two generation are neglected, so that the corresponding diagonalizing matrices are trivial. Details can be found in appendix B.1.3.
We then apply O(α S ) corrections to the squark squared masses (consistently with what was implemented in the original CP-conserving treatment in NMSSMTools). Gluons, gluinos as well as the quartic sfermion D-term contribute to the sfermion self energy at this order. CP-phases -here, φ M3 and ϕ A f -intervene in the gluinosfermion couplings leading to a B 0 function. A summary is proposed in appendix C.2.
Finally, we check the positivity of the sfermion squared masses, a vacuum-stability requirement.
iv) mgluino_CPV.f mgluino_CPV.f computes the gluino mass, including the O(α S ) radiative corrections, which are obtained in a similar manner to the discussion in section 4.1 of [20]. Relevant corrections include the gluon / gluino and the quark / squark loops. Complex phases again enter the couplings of gluinos to squarks. Details are provided in appendix C.3.

Higgs masses and radiative corrections
The following series of subroutines aim at computing the Higgs masses and mixing, including full one-loop and leading two-loop corrections. Consistently with the original approach in NMHDECAY, we will consider the effective Higgs potential at the average scale of the squarks of the third generation -denoted as Q -, where the running parameters are thus evaluated.

Wave-function renormalization
Momentum-dependent radiative corrections can be included in two fashions within the effective potential evaluation: one may reject them to the end of the calculation, as 'pole-corrections', or one may take them into account -at least partially -into the effective lagrangian as corrections to the kinetic terms. The latter choice leads to wave-function renormalization factors. While the two methods are formally equivalent, they lead to slightly divergent results at the numerical level, as we will discuss later. Following the original approach in NMHDECAYpresented e.g. in appendix C of [21] or appendix C of [6] -, we decide to include the leading p 2 terms -where p stands for the external energy-momentum of the Higgs self-energies -, originating in fermion or gauge effects, into the kinetic term of the effective lagrangian. Nevertheless, since we aim at a full computation at one-loop, all the missing momentum-dependent parts will be added as pole-corrections (see below).
In the general case, the modified Higgs kinetic terms involve a hermitian (non-degenerate) matrix Z H (p 2 ) as follows (here and below S i denotes any Higgs field; we work in momentum space and omit the factor 1/2 which should appear if the considered field is real): The normal procedure then consists in rotating and scaling Z H via an invertible matrix O H in order to recover the identity -Z H (p 2 ) = O † H 1O H -, then considering the 'new' set of fields with standard kinetic termS i ≡ O H ij S j . Yet, Eqs.(C.1) of [21] or (C.9-11) of [6] show that a clever choice of the corrections included into Z H can make this procedure particularly simple, as Z H would turn out to be diagonal in the base of gauge-eigenstates. Restricting to neutral Higgs fields, one has (with δ Si,Sj denoting the Kronecker symbol): Indeed, considering the contributions of SM-fermions to Z H (N c = 3 is the colour factor; while using the generic notations u, d, e, we will be considering only the third generation fermions since we neglect the Yukawa couplings of the two first families), the deviations of the diagonal scaling factors from unity read 2 : Similarly, in the approximation where higgsinos and gauginos are simultaneously gauge and mass eigenstates (µ denotes the doublet higgsino mass; ms, the singlino mass): The last source of corrections to Z H is the gauge sector -note that we will be working in the Feynmann gauge. Yet, the corresponding contributions are not diagonal in the gauge eigenbase, but rotated by an angle β (or −β, depending on the CP-eigenvalue) in the doublet sector. Noticing however that tan β > 1 in practice, we may keep the sin 2 β term in the wave-function scaling while rejecting the remaining sin β cos β and cos 2 β terms for later treatment as pole-corrections. Then: Note that this choice in the gauge sector differs from the default treatment by NMSSMTools in the CP-conserving case (see Eq.(C.9-10) of [6]), where, moreover, pole corrections from the gauge sector are ignored.
Before setting Z Hu,H d ,S = 1 + δ SM ferm + δh ,g + δ gauge Z Hu,H d ,S , one is confronted to the remaining p 2dependence of these coefficients (via the loop functions B 0 ). In the ideal case, p 2 would match the Higgs squared masses. This, however, is impractical since several mass eigenvalues are present: keeping this p 2 dependence, hence working with p 2 -dependent fields and mass-matrices, and setting this implicit dependence separately to the corresponding Higgs squared mass after diagonalization of the mass matrix would be possible, yet problematic in a numerical evaluation of the mass matrices. The choice of [6] in the CP-conserving case rested in adding an artificial dependence of Z Hu,H d on ln(M 2 A /m 2 t ) -M A standing for the mass of the heavy doublet, m t approximating the SM-like Higgs mass -, so as to mimic the correct logarithmic dependence after rotation by an angle −β (approximating the tree-level diagonalizing rotation in the CP-even doublet sector): however an explicit rotation by the angle −β shows that this purpose is missed as only the light state receives the proper logarithmic factor; in the case of the heavy doublet, the factor is wrong so that the result does not really improve on neglecting the logarithms ln(M 2 A /m 2 t ) altogether. Therefore, we settle for the choice which consists in freezing the external momentum to a scale µ H = 125 GeV, allowing for a good precision in the characteristics of the SM-like Higgs state -the most sensitive to radiative corrections. Adequate corrections when the mass is far from this scale are rejected to the level of pole-corrections. A final difference with [6] comes from the implementation of the loop functions: we explicitly compute the full relevant B 0 's while [6] only included the leading logarithmic terms in case of large mass hierarchies.
A summary of the wave-function scaling factors is provided in appendix D.1. Consistently, the neutral higgs fields are rescaled as: so that all related quantities (e.g. the mass matrices) must be rescaled accordingly. In particular the Higgs v.e.v.'s: All these operations are carried out in the initialization subroutine init_CPV.f. In the charged-sector, the p 2 -dependent terms are typically different from those appearing in the neutral case. However, to keep β as the relevant rotation angle in the charged sector together with the v.e.v. rescaling of Eq.27, we will use the same wave-function scaling factors Z Hu,H d : This concludes the presentation of the general formalism and we may now describe the various contributions to the effective potential which are computed within our code.
ii) mhiggstree_CPV.f This subroutine simply defines the tree-level mass matrices at the scale Q according to Eqs.12 and 13. However, the corrected Higgs masses are not the only information that we want to extract from the effective potential: the Higgs-to-Higgs couplings are also encoded within this formalism. Therefore, and for reasons that will become clear when we implement the various radiative contributions to the potential, we wish to match the full effective potential onto the following and simpler one: This is a subset of the most general singlet + two doublet potential which one can write up to dimension 4 terms 5 . The gauge symmetry is observed. However the Z 3 -symmetry only holds up to terms quadratic in the doublet fields and is explicitly broken by the terms in the last line of Eq.33. This potential is meant as an expansion of Eq.29 in the doublet fields and as we mentioned before, there is no reason why the Z 3 -symmetry should hold in such an expansion. The characteristics of this potential are studied in appendix E and matching the tree-level expression of Eq.9 is straightforward (see appendix E.1).
iii) mhiggsloop_sferm_CPV.f With this subroutine, we start adding radiative corrections to the effective potential, here those arising from SMfermion and sfermion loops. These -particularly the contribution associated to the top -are known to convey the dominant radiative effect and lead e.g. to a substancial shift of the squared-mass of the SM-like Higgs boson.
The corresponding one-loop effects to the neutral Higgs mass matrix are particularly easy to include in the Coleman-Weinberg formalism of Eq.29, since the bilinear terms provide relatively simple matrices (refer to the appendices B.1.1 and B.1.3). The details of the corrections are developed in the appendices D.2.1 and D.2.3. Note that we recover Eqs.(C16-18) of [6] in the CP-conserving limit.
The situation is slightly more complex for the charged Higgs as well as for the Higgs-to-Higgs couplings: we then decide to expand the potential in terms of the doublet fields, up to quartic order H 4 and match the corresponding expansion onto the simplified potential of Eq.33. This amounts to an expansion in v/M SUSY , where M SUSY here stands for any sfermion mass. The sfermion contributions to the coefficients of Eq.33 are also provided in appendix D.2.3. Note that this alternative approach allows for a numerical cross-check with the corrections applied to the mass matrix of the neutral Higgs states with the method described in the previous paragraph.
In addition to these one-loop effects, the subroutine mhiggsloop_sferm_CPV.f also includes two-loop effects of order 6 O(Y 6 t,b , Y 4 t,b α S ) leading to a product of large logarithms in the fermion sector: given that we are working at the average scale of squark masses, the squarks are assumed to give subleading contributions. On the other hand, effects associated to SM fermions and gauge bosons will not introduce any additional dependence on the new physics phases. The corresponding effects are implemented in the approximations of [22] -see also Eq.(C.19) of [6] -, i.e. only the contributions to the quartic doublet parameters λ u and λ d of Eq.33 are included -note that contributions to M 2 u or M 2 d leave the analysis unaffected, while contributions to e.g. A ud can be absorbed in a shift of the tree-level term A λ , hence only drive a displacement in the parameter space. While these contributions are of two-loop order, they may still affect the mass of the SM-like Higgs state by several GeV, which is why we include them. Comparisons to more-elaborate two-loop calculations show that this approximation works well numerically (at the GeV level).
iv) mhiggsloop_inos_CPV.f The next subroutine implements the radiative effects associated to charginos and neutralinos. Sticking to the Coleman-Weinberg approach, we consider the 9 × 9 bilinear term associated with gauginos and higgsinos (refer to appendix B.1.4). Due to the large rank of this matrix, we exclusively employ the method which consists in expanding the potential and matching it to the simplified version of Eq.33. The corresponding results are collected in appendix D.2.4. Note that they differ from e.g. Eq.(C.22-24) of [6] where additional simplifying assumptions had been made.
v) mhiggsloop_gaugehiggs_CPV.f The contributions of the electroweak gauge bosons to the Higgs potential seem easy to include in the Landau gauge: see appendix D.2.2. Yet the drawbacks of the Landau gauge are felt in the Higgs sector, where one then has to handle massless Goldstone bosons. The associated infrared divergences are of course purely spurious and disappear once confronted to momentum-dependent corrections, as already noted in [23]. Still it remains a technical issue to manipulate with caution. Moreover, the strategy consisting in diagonalizing the field-dependent bilinear matrices, which we have been employing until here, becomes impractical, even in an expansion in terms of the doublet fields, due to the large number of parameters and operators involved in the Higgs bilinear terms. Instead, we decide to employ the concurrent strategy in Higgs-mass calculations, which simply consists in a direct diagrammatic evaluation of the Higgs self-energies and tadpoles generated by Higgs loops. Nevertheless, disentangling the Higgs and gauge contributions in this approach proves quite artificial so that we will lead the calculation for both types of particles appearing in the loop at the same time.
Explicit expressions for the gauge and Higgs one-loop contributions to the Higgs self-energies and tadpoles are summarized e.g. in [20] or [24] (with different conventions for the loop functions), in the context of the MSSM, and the NMSSM differs only in the definition of the couplings and the presence of the singlet fields, hence leads to a formally comparable result. We choose to work in the Feynmann gauge as it is then possible to set the external momentum to 0 without generating IR-divergent logarithms. Indeed, we still aim at computing, not only the corrections to the Higgs masses, but also to the Higgs-to-Higgs couplings. For this, we proceed in the following fashion: after the radiative corrections to the Higgs mass matrices are evaluated at zero momentum, we substract the pure-gauge contribution in the Landau gauge (for which we already know the potential from appendix D.2.2). The remaining 'Higgs' contributions to the mass matrices can then be identified with those that a Z 3 -conserving renormalizable potential would produce, allowing for a reconstruction of the corrections to the Z 3 -conserving parameters of the potential: this procedure is described for the CP-conserving case in [25] and is straightforwardly generalized to the CP-violating case. Of course, we then miss contributions to the Z 3 -violating parameters (the last line of Eq.33) but, as discussed in [25], these are subleading in the leading-logarithmic approach 7 . Further details can be found in appendix D.2.5.
This completes the list of radiative contributions implemented in the effective potential.

Pole corrections
The operations described in the previous lines have provided us with mass matrices for the Higgs states where radiative corrections from the potential (i.e. at zero external momentum) have been included. We will now detail how we account for momentum-dependent corrections. These calculations are conducted in the subroutine mhiggsloop_pole_CPV.f. First, let us remind the reader that the radiative effects associated with non vanishing external momentum have been partially encoded into the wave-function scaling factors of paragraph 3.3.1. It is necessary to rescale the Higgs mass-matrices in order to account for the re-scaling of the Higgs fields: While ideally the Higgs self energies Π S 0 should be evaluated at the pole masses, we approximate the latter by the DR masses. The full one-loop pole-corrections are applied. Shifts of the wavefunction scaling factors δZ Hu,H d ,S are know from Eqs.23-25. The shifts in the Higgs self energies are provided in appendix D.4.
This concludes our evaluation of the masses in the Higgs sector. We now wish to comment briefly on the precision achieved in this calculation. For this, it is instructive to consider the impact of the one-loop corrections with respect to the situation at tree-level. For mostly-doublet states, the leading effect is driven by the topquark loop and, respectively to the tree-level mass m H , can be quantified as ∼ , this amounts to a correction at the percent level for m H = O(TeV), but reaching a magnitude of ∼ 100% for m H = O(100 GeV): this accounts for the well-known sensitivity of the SM-like Higgs mass to radiative corrections. Contributions at the two-loop order will involve the strong coupling g S , or the top Yukawa coupling again, multiplying logarithms of a similar magnitude, so that the typical effect would easily amount to ∼ 30% of the one-loop contribution. Now, considering that we have included the leading double-logarithmic effects in the calculation, we can estimate a reduced uncertainty from higher orders, say at the level of ∼ 10% of the one-loop corrections. For a Higgs mass at ∼ 125 GeV, this still amounts to an uncertainty of several GeV. For a state at m H = O(TeV), this reduces to the permil level. The latter accuracy is treacherous however, as other sources of uncertainty appear e.g. in the determination of the couplings or neglected electroweak corrections entering the definition of the Higgs v.e.v.'s. In the outcome, the precision on the Higgs masses should not improve on O(1%) for heavy doublet states. Corrections to singlet states are typically smaller, since the associated couplings -λ, κ -are of order < ∼ O(0.1) and the hierarchies between Higgs bosons and higgsinos may not be as large as those between SM fermions and sfermions. However, when the singlets mix significantly with doublet states, they will correspondingly acquire part of the larger uncertainties on doublet masses.

Couplings, decays and constraints
After the previous subroutines are run, one has a complete set of corrected masses and rotation matrices at one's disposal. The following move consists in confronting this spectrum to physical processes.

Supersymmetric and Higgs couplings
The couplings of supersymmetric particles and Higgs bosons can take somewhat lengthy expressions. We thus design two subroutines, susycoup_CPV.f and higgscoup_CPV.f, in order to evaluate and store them within the code: Note that the rescaling of Higgs fields in Eqs.26 and 28 is also accounted for when computing the couplings of Higgs bosons.

Higgs and top decays
We then adapt the existing NMSSMTools subroutines decay.f and tdecay.f -respectively computing the Higgs and top two-body decays in the CP-conserving NMSSM -to the CP-violating case.
The subroutine hidecay_CPV.f calculates the Higgs widths and the dominant branching ratios. The following decay channels are considered: • decays into a pair of SM fermions: • decays into (on-shell) gauge bosons 8 : S 0 i → W W , ZZ, γγ, Zγ, gg; • decays into one Higgs and one gauge boson: • Higgs-to-Higgs decays: In the subroutine tdecay_CPV.f, we compute the following top decays: t → W + b, H + b,T χ 0 . As in the original CP-conserving version, leading QCD corrections have been taken into account.

Phenomenological tests
We finally propose several tools to confront the CP-violating NMSSM spectrum to experimental constraints.
checkmin_CPV.f compares the value of the neutral effective potential at the electroweak symmetry-breaking minimum with that at other points, e.g. for vanishing v.e.v.'s. Loop effects from the SM fermions and gauge bosons are included explicitely in this evaluation, while other radiative effects are encoded within the approximate potential of Eq.33. We also vary the dynamical phases and check whether this generates a deeper minimum. Finally, the minimization conditions of Eq.11 and 30 are calculated explicitly, which allows e.g. to test the naturalness of the squared masses m 2 H u,d of the potential: they should remain of the order of the SUSY-breaking scale. In constsusypart_CPV.f, we generalize to the CP-violating case LEP limits on superparticle searches that were included in NMSSMTools for the CP-conserving case: • test on chargino, slepton, gluino and squark masses; • limits onT → blÑ ,T → cχ 0 ,B → bχ 0 ; • constraint on the invisible Z-width and neutralino-pair production.
HBNMSSM_CPV.f converts our spectrum into input for HiggsBounds [18] and HiggsSignals [19]. This allows to test the Higgs sector in view of LEP, TeVatron and LHC results via a call to the subroutines included within these public tools -note that NMSSMTools, HiggsBounds and HiggsSignals must be interfaced to make use of this subroutine. The chosen input mode is that employing effective couplings (see the documentation in [18]). Additional widths and branching ratios are taken from the results in hidecay_CPV.f and tdecay_CPV.f. In the following section, we will be using the current versions HiggsBounds_4.2.0 and HiggsSignals_1.3.1, which incorporate all the experimental results released till december 2014. The default uncertainty on the Higgs mass precision is set to 3 GeV and modeled as a gaussian distribution.
We also include an alternative set of tests for the Higgs sector, based on the original subroutines of NMSSMTools. These collect: • LEP_Higgs_CPV.f: LEP limits applying on neutral Higgs bosons produced in association with • TeVatron_CHiggs_CPV.f: TeVatron limits applying on a charged Higgs boson produced in top decays [27]; • bottomonium_CPV.f: test for a light mostly CP-odd Higgs in bottomonium decay -based on [28]; • LHC_Higgs_CPV.f: the inclusion of LHC limits on neutral or charged Higgs searches as well as the confrontation to the signals at ∼ 125 GeV -after [29] -are in progress.
Note that these routines will not be used in the next section, as we will employ the currently more complete set of tests performed by HiggsBounds and HiggsSignals. Finally, we design a subroutine EDM_CPV.f to estimate the electric dipole moments of the electron, the thallium atom, the neutron and the mercury atom. We essentially follow the summary in [30] -in the context of the MSSM; see also [31] for a recent work in the NMSSM. The supersymmetric one-loop effects are mediated by charginos, neutralinos or gluinos and sfermions. Moreover the two-loop diagrams of the Bar-Zee type -involving a fermion or sfermion loop connected to the quark / electron line by a Higgs and a photon propagator -are known to convey a sizable effect: these are particularly sensitive to the phases appearing in the Higgs sector. Other contributions, mediated by dimension 6 operators, are included as well. We estimate the associated uncertainties by adding linearly a 10% error on effects involving no coloured particles and a 30% error on contributions involving the coloured sector. Additional uncertainties associated to scale-running or hadronic parameters are also incorporated.

A few applications
We shall now make use of the subroutines which we have just presented and study phenomenological effects associated with the CP-violating NMSSM. This will be the opportunity to test our tool and compare its predictions with existing results.

CP-conserving limit
Setting all the phases to zero, it is possible to consider the CP-conserving case: in particular this allows to study how our results connect to the precision calculations implemented within NMSSMTools. Given that the input is common, discrepancies directly give an insight on the differences of treatment and the numerical magnitude of the corresponding effects.

i) Higgs spectrum
We shall first consider the Higgs masses. NMSSMTools provides three levels of precision in the inclusion of the radiative corrections to the Z 3 -conserving Higgs sector: • 'Precision 0': the default one -essentially following the procedure described in appendix C of [6] -confines to leading logarithmic order. Momentum-dependent effects are taken into account only to the extent of wavefunction renormalization (where the implementation is slightly different from ours: remember the discussion in section 3.3.1) and pole-corrections associated with the SM-fermion sector.
• 'Precision 2': a full one-loop + leading two-loop (to order O(Y 4 t,b α S )) implementation including momentumdependent effects. It follows the work of [32].
Formally, our implementation -full one-loop including momentum-dependent corrections + leading two-loop double logarithms of order O(Y 6 t,b , Y 4 t,b α S ) -should fall somewhere between these three procedures in terms of precision. We In the down right-hand plot, the singlet composition S 2 i3 of the two lightest CP-even states is displayed for precision 0 (green solid lines), precision 2 (blue solid lines) and for our calculation (dotted lines).
A t = −2 TeV, A b,τ = −1.5 TeV and we scan over λ ∈ [0, 0.65]. The Higgs masses are displayed in Fig.1 and 2: the results of NMSSMTools for precision '0' (greenish colors), '1' (pink colors) and '2' (bluish colors) are shown as solid lines while our calculation corresponds to the dots (yellow to red tones). We observe a significant variation of the masses corresponding to the mostly-singlet states while the doublet masses are grossly constant with varying λ. A typical NMSSM effect develops when singlet masses are close to doublet masses, as significant mixing may appear. In particular, when the singlet state is slightly lighter than the doublet one, the mixing tends to uplift the mass of the mostly-doublet Higgs. This is what occurs in this example for the CP-even sector in the upper range of λ. In Fig.1, we see that our results fit quite closely the predictions of the procedure with precision 2, while larger discrepancies appear with respect to precision 0, especially at large λ.  Fig.1. On the upper left-hand quadrant, we show the Higgs masses close to ∼ 125 GeV for precision 0 (green lines), precision 2 (blue lines) and for our implementations (red dots). The plot on the upper right-hand side compares our results (dots) for the 'heavy' masses with those of precision 0 (green lines). The same exercise is carried out in the lower left-hand corner for precision 2 (blue lines). In the lower right-hand quadrant, we alter our implementation of the Higgs mass corrections so that all p 2 -dependent terms are taken into account as pole-corrections only (so that the wave-scaling factors are set to 1): the results are displayed as khaki dots and compared to the masses obtained in the procedure of precision 2.
of the order of several GeV. However, it justifies the observation that the leading two-loop effects are captured by the simpler inclusion of double logarithmic terms.
Concerning the heavy mass states, we observe in Fig.2 -in the upper-right and lower-left hand quadrants -that our results are intermediary between the calculations with precision 0 and precision 2. However, we note that the leading difference with precision 2 originates in the implementation of the wave-function scaling factors. Indeed, if we set the 'Z-factors' to 1 and modify the pole-corrections accordingly, we observe that our result -corresponding to the khaki dots in the lower-right-hand corner of Fig.2 -then matches that with precision 2 somewhat more closely (at the permil level). It is quite easy to see how the discrepancy develops between these two procedures. For this, let us focus on the CP-odd doublet state, where we will neglect the mixing with the singlet. In the case where the wave-function scaling factors are set to 1, the squared mass of this state is -schematically: the effect of potential and pole corrections are encoded as δ pot,pole -obtained as (all the p 2 -dependent terms are treated as In the approach where the wave-function scaling factors are taken into account at the level of the kinetic terms, the Z-factors intervene in the calculation at several steps: first, for the scaling of the v.e.v.'s, which transforms the tree-level mass-matrix in the CP-odd doublet sector to: The scaling effect on the potential corrections can be neglected as being of higher order. Then comes the scaling of the mass-matrix: so that we can extract the DR squared mass for the physical state via a β-angle rotation. Finally, the Z-factors have to be substracted from the pole corrections, since they have been accounted for elsewhere. This provides: Expanding the Z-factors as Z · = 1+δZ · , we see that Eqs.36 and 39 differ by a factor 1− δZ S 2 + cos 2β 2 (δZ Hu −δZ H d ). This explains the mismatch, reaching the order of one-loop effects, that is O(1%) here. In particular, the steeper apparent slope with varying λ, in Fig.2 is largely driven by the Z S factor. In principle, the approach including the wave-function scaling factors is the most refined among the two methods, hence should be prefered. On the other hand, our choice of setting the Z-factors at a low-value of the external momentum, µ H = 125 GeV, is not optimized for heavy states. In any case, a 1% effect should not be taken too seriously in view of the various additional sources of uncertainty (parametric errors, running, etc.).
We then consider a second example with λ = 0.7, Fig.3. This region of the parameter space highlights another effect in the NMSSM Higgs sector, namely the large contribution of F-terms to the mass of the SM-like state for large λ and low tan β. Indeed, the low value of tan β, the low mass of the squarks of third generation and the moderate trilinear soft terms would result in a Higgs mass below M Z in the MSSM, making this regime incompatible with LEP limits and the LHC measurement. In the NMSSM however, we observe that the mass of the SM-like state remains above 120 GeV: this is a consequence of the specific tree-level contributions to the Higgs mass matrices, associated with λ. Comparison of our results with the masses obtained with NMSSMTools for precision settings 0 and 2 again show that our calculation is typically closer to precision 2, although the differences are larger than in the previous scan (about 1 GeV for the two light CP-even states, as can be observed on the plot on the lower left-hand corner). We also display the output of HiggsBounds and HiggsSignals for our results (plot on the lower right-hand side): HiggsBounds exclusions apply e.g. in the presence of very light Higgs-states with non-vanishing doublet composition. The χ 2 test of HiggsSignals provides values down to ∼ 75 -for comparison, we obtain ∼ 78 in the SM limit -when a light doublet is present close to ∼ 125 GeV.

ii) Higgsino and gaugino masses
Our implementation of the chargino, neutralino and gluino masses should prove very similar to the original subroutines within NMSSMTools in the CP-conserving limit. Nevertheless, small technical differences should be noted: • we take into account the Higgs-higgsino-singlino couplings which had been neglected in NMSSMTools: this results in additional corrections to the higgsino and singlino masses; • similarly, bino and winos are not assumed degenerate in the calculation of loop corrections to the higgsino masses; • all masses are chosen real and positive: this is possible since the diagonalizing matrices are complex. The convention in NMSSMTools consisted in keeping these matrices real, so that some masses could take negative values.
We consider the following region in the NMSSM parameter space: The masses of the higgsinos and gauginos are shown in Fig.4. The scan over µ eff drives a significant variation of the higgsino masses, while the gaugino masses remain essentially constant. Once again, the masses obtained with the original routine of NMSSMTools are depicted with a solid line, whereas our results appear as dots: the general features are identical. More quantitatively, the main deviation reaches ∼ 3% at the level of the neutralino masses: it originates from the corrections to the singlino mass, which were neglected in NMSSMTools.
For the rest of the spectrum, e.g. the sfermion masses, our calculation reduces, in the CP-conserving limit, to the original implementation within NMSSMTools. Therefore, we will not push the comparison in this limit any further.

CP-violating case
CP-violation could induce several phenomenological effects at colliders. The most immediate one would be the measurement of EDM's. The absence of any hint in corresponding searches thus places stringent limits on newphysics phases. Note however that, at one loop order, these effects are essentially driven by the gaugino phases. In other words, new-physics phases associated to the Higgs sector or the third generation sfermions enter the EDM's at the two-loop level only and are thus more loosely constrained. CP-violation could also intervene in rare flavour decays and oscillations, which are consistent so far with the SM-interpretation (where only the CKM phase is present): such effects have not been included in our study yet and we will not discuss them here.
i) CP-violating effects in the NMSSM Higgs spectrum CP-violation could enter the Higgs sector at tree-level, via a non-vanishing phase ϕ λ − ϕ κ , or at the loop-level, e.g. via the phases associated to the sfermions of third generation. As a first consequence, the neutral Higgs states would become scalar / pseudoscalar admixtures, which affects their couplings to SM particles: for doublet states, the pseudoscalar component does not couple to a ZZ or W + W − pair, so that the corresponding decay channels, as compared to the fermionic decays, are suppressed / enhanced with respect to the case of pure CP-even / CP-odd eigenstates. Other effects can be measured in the fermionic channels, provided, however, that the fermion masses are sufficiently large. Therefore the presence of CP-violation in the Higgs sector could be tested in precision analyses of the Higgs properties -for the observed or hypothetical new states. Note however that doublet Higgs states are typically shielded from CP-violating mixing -consider e.g. the zero-entries in the tree-level mass-matrix of Eq.13 -, so that only a very high degree of precision in the measurement of the branching ratios would be likely to detect the tiny -radiatively-generated -pseudoscalar component of a mostly CP-even state. Moreover, the current limits on Higgs searches tend to favour a sizable mass-hierarchy between the SM-like Higgs state and the approximately degenerate 'heavy-doublet' states. This makes the presence of a pseudoscalar doublet component within the observed Higgs state unlikely, as the mixing of this state with the 'heavy-doublet' pseudoscalar would be suppressed in proportion to the large mass gap. Another test would involve the two 'heavy-doublet' neutral states, which are generically close in mass, so that their mixing could be significant. Yet, the detection of CPviolation there will still require high-precision experiments (and the discovery of these states), due to a typically reduced production cross-section -with respect to a SM Higgs boson at the same mass; this is related to the mostly H d -nature of these states -as well as the opening of many less-controlled decays (e.g. towards new-physics states). In the NMSSM, another type of CP-violating mixing is allowed: a mostly CP-odd singlet may mix with the doublet CP-even states -provided λ and κ are large and ϕ λ − ϕ κ is non-vanishing -and this effect could be fairly important if these states are close in mass. In the following, we focus on the SM-like Higgs state at ∼ 125 GeV. Such a scenario is studied in Fig.5 for λ = 0.68, κ = 0.1, tan β = 2, µ eff = 635 GeV, M A = 1.5 TeV, mT ,B,τ = 0.5 TeV, mŨ ,D,Ẽ = 1.5 TeV, 2M 1 = M 2 = M 3 /3 = 0.5 TeV, A t,b,τ = −0.1 TeV. CP-violation is induced through variations of ϕ κ : note that this strategy is the safest in view of the EDM's, as non-vanishing ϕ λ produces direct CP-violation in the doublet higgsino sector (as well as in the sfermion sector). In the first column of Fig.5, A κ = −100 GeV, and the mostly CP-odd state is relatively far in mass (∼ 250 GeV): correspondingly, the mixing with the SM-like state does not reach 1%. The latter state has a somewhat low mass of ∼ 121 GeV which translates into a mediocre fit to the LHC-observed signals, hence a high χ 2 -value with HiggsSignals. In the second column, we take A κ = 0 GeV, so that the CP-odd singlet is close in mass to the SM-like state: at ϕ κ = 0, the singlet has a mass of about ∼ 115 GeV. Consequently, a significant mixing develops between the two light states as soon as ϕ κ = 0, the effect reaching the level of 30 to 40%. A consequence is the uplift in mass of the heavier SM-like state so that the associated signal gives an improved fit with the LHC data. The column on the right is obtained for A κ = 10 GeV: the CP-odd singlet is then somewhat lighter (∼ 100 GeV), so that the mixing effect at non-vanishing ϕ κ remains milder than in the previous case, yet generates an uplift of the mass of the SM-like state as well. It is to be noted that the mostly CP-odd singlet acquires a CP-even doublet component which reaches O(10%): the latter would generate a signal at the O(10%)-level as compared to a SM-like state at the same mass -indeed, the production cross-section at colliders is essentially mediated by the doublet components. For a state with mass ∼ 100 GeV, the corresponding signal could be consistent with the LEP ∼ 2.3 σ excess in Higgs searches with a bb final state [26], even though the state is dominantly CP-odd. Note that the two effects that we highlighted -uplift of the mass of the SM-like state via its mixing with the singlet and presence of a 'miniature' Higgs boson under 125 GeV -are well-known in the CP-conserving NMSSM [7], provided the auxiliary singlet is CP-even. CP-violation extends this possibility to CP-odd singlets. Further consequences appear on Fig.6 at the level of the branching fractions of the Higgs states -we display their values for the bb, cc and γγ final states -: similarly to the case where the SM-like Higgs boson mixes with a CP-even singlet, the proportions among doublet components h 0 u and h 0 d may fluctuate, displacing the branching ratios. However, the main effect in Fig.6 concerns the rates of the lighter singlet state which become dominated by CP-even-like channels -for fermionic final states, rates differ at the radiative level depending on the CP property 9 -, while the fluctuations of the branching fractions of the mostly CP-even doublet are dominated by the variations of the associated Higgs mass.
Disentangling this scenario -where a light mostly CP-odd singlet mixes with the SM-like Higgs boson -from the CP-conserving one -where the light singlet-like state is genuinely CP-even -is likely to prove very difficult. The reason rests with the observation that the singlets do not lend specific properties to the SM-like Higgs statethey simply reduce its total width and might alter its branching ratios at the percent level. Moreover their decays are essentially mediated by the doublet component which they acquire in the mixing, i.e. a CP-even one in both cases. Typical singlet decays -towards hypothetically lighter singlet states or singlinos -would not necessarily help to discriminate among CP-even and CP-odd mixing and would be problematic in terms of compatibility with the measured Higgs signals. Indeed, the standard rates would then be suppressed in proportion of the magnitude of the unconventional decays. While deviations of the rates of the observed Higgs state from the standard ones might be interpreted via such a mixing effect -should such deviations be detected at the LHC or a future linear collider -, it is questionable anyway whether the light singlet could be detected -possibly in Higgs-pair production: see e.g. [33] in the CP-even case.
At the outcome of this discussion, we see that, while the CP-violating effects involving singlets in the Higgs sector may be larger than in the pure doublet case, they are also more difficult to trace and could be mistaken for CP-conserving phenomena. For this reason, it is essential that CP-violation be tested in processes where the CP-properties are well-controlled, which brings us back to EDM's or rare flavour transitions. Spectral effects in the Higgs sector are unlikely to allow for discrimination with the CP-conserving case.
ii) Comparison of the Higgs mass predictions with the existing literature We will now compare some of our results with existing analyses in the literature, where CP-violation has been considered. Note that, contrarily to the comparison with the calculations in the CP-conserving NMSSMTools, one should not expect much more than a qualitative agreement. Indeed, diverging treatments in different tools, e.g. concerning the definition of the input -such as the choice of running Yukawa couplings or that of A κ versus A κ cos ϕ 2 -, are known to lead to sizable deviations, already in the CP-conserving case. The level of precision in radiative corrections is also to be considered.
NMSSMCALC [14] is a public tool computing the Higgs spectrum and decays in the Z 3 -conserving but possibly CP-violating NMSSM. The chosen approach is that of a diagrammatic calculation. The level of precision has recently been extended to include the dominant two-loop corrections [34].  Fig.6 of [34]: green for a scan over ϕ µ , red for a scan over ϕA t and blue for a scan over φ M3 . When the blue curve does not appear, the reason is that associated values are negligibly small. Note that several estimates are employed for the neutron EDM.
First, we focus on the results of [34] dealing with CP-violating effects, i.e. essentially Fig.6 and the surrounding text in that paper. If we blindly input the parameters given in section 4.1 of this reference into our framework 10 , the spectrum -not unexpectedly and already with the CP-conserving NMSSSMTools -turns out to be slightly different from the quoted one: in particular, the mostly CP-even and mostly CP-odd singlet states appear with masses ∼ 108 GeV and ∼ 113 GeV respectively. Yet, this discrepancy can be absorbed within a small shift of A κ cos ϕ 2 : using the value 203 GeV, we then recover states at ∼ 103 and ∼ 128 GeV so that the Higgs spectrum then largely coincides with the one provided in [34].
In any case, this manipulation has little effect on the properties of the mostly h 0 u -state, close to 125 GeV. Scanning over the phases ϕ At , ϕ µ -in the notations of [34], this means a scan over ϕ λ , keeping ϕ κ = ϕ λ so that CP-violation does not enter the Higgs sector at tree-level -and φ M3 , we obtain the plots of Fig.7. On the upper part, we observe that the general dependence of the 'SM-like' Higgs mass on ϕ At and ϕ µ is largely reminiscent in shape and magnitude of that observed in Fig.6 of [34]. In these two cases, CP-violation enters the Higgs sector via radiative corrections, where the leading effect is generated by the sfermion corrections. On the other hand, the mass obtained with our code is independent from φ M3 , while such a dependence already appears at one-loop in [34]. Note that one does not expect gluino corrections to the Higgs mass at one-loop order and it is thus not surprising that our implementation does not show any variation with φ M3 . The corresponding effect in [34] is explained there as an artifact of the top-Yukawa DR counterterm-fixing of higher order. Note also that the corresponding fluctuations, at the GeV level, are small compared to the uncertainty that one naively expects for the Higgs mass (a few GeV).
In addition, we show the values of the EDM's that we obtain in these scans. These have been normalized by the experimental upper bounds: ∼ 1 · 10 −28 e cm for the electron [35] -estimate from thorium monoxide experiment -, ∼ 1.3 · 10 −24 e cm for the Thallium atom [36], ∼ 3.1 · 10 −29 e cm for the Mercury EDM [37] and ∼ 3 · 10 −26 e cm for the neutron [38]. Note that only the central values are displayed, without error bands. The color code is the same as in Fig.6 of [34], i.e. green for the scan on ϕ µ , red for that on ϕ At and blue for the one over φ M3 (when the curve does not appear in the plot, this is because the corresponding values are negligibly small). We see that the scan over ϕ µ may generate tensions with the EDM's -mostly the electron EDM -when ϕ µ is not trivial. We now turn to the one-loop analysis proposed in [39]. We first consider the scenario presented in section 4.1.1 of this reference, where CP-violation intervenes in the Higgs sector at tree-level via the phase ϕ κ . Again, a qualitatively close spectrum can be recovered with little alteration of the input proposed in the reference and our results are displayed in Fig.8: while small differences appear, both the Higgs masses and the composition of the states agree reasonably well with those of [39]. The major source of deviation is associated to the use of different input -A κ in NMSSMCALC instead of A κ cos ϕ 2 in our case -, so that the comparison makes limited sense when ϕ 2 becomes large (i.e. for ϕ κ ∼ π/8). In the regime considered here, the CP-even and CP-odd singlet states are close Figure 9: Masses in the scenarios of section 4.1.2 (scan over ϕ µ ≡ ϕ λ = ϕ κ ) and section 4.1.3 (scan over ϕ At ) of [39].
in mass to the SM-like Higgs boson, so that the non-vanishing ϕ κ generates a substantial mixing of these three states.
[39] then considers the case where CP-violation is absent in the tree-level Higgs sector, but radiatively generated via phases in the supersymmetric spectrum. In the first case (section 4.1.2), the 'active' phase is ϕ λ but the condition ϕ κ = ϕ λ ensures that no CP-violation enters the Higgs potential at tree-level -we will recycle the previous notation ϕ µ for this scenario. In the second case, only the phase ϕ At is non-trivial. We display our results in Fig.9 and observe that they capture the effects depicted in Fig.5 and Fig.7 of [39].
Our code is thus able to reproduce the main qualitative features that were observed in the CP-violating case by NMSSMCALC analyses. We stress that a more quantitative study would have limited interest, as the divergent treatment of the input already generates discrepancies between the CP-conserving NMSSMTools and NMSSMCALC.

Conclusions
We have presented a series of Fortran tools extending NMSSMTools to the CP-violating case. Radiative corrections to the supersymmetric and Higgs masses are computed at one-loop order. Dominant two-loop effects to the Higgs masses are also included in the double-log approximation. Additionally, Higgs couplings and decays, as well as top two-body decays and EDM's are implemented and allow for phenomenological tests of the spectra. We have shown that our code compares competitively with existing results, both in the CP-conserving and CP-violating cases. The new tools will be made publicly available on the NMSSMTools website [13] in the near future.
We also highlighted a scenario made possible by CP-violation, where the SM-like Higgs would mix with a mostly CP-odd singlet state. The consequences on the Higgs phenomenology are similar to the CP-conserving mixing with a light CP-even singlet so that both scenarii should prove difficult to discriminate, unless genuine CP-violating effects -e.g. in EDM's or flavour physics -are discovered simultaneously.
Finally, we would like to close this discussion with some details concerning the future developments which we plan to consider. First, an extension of our tools including Z 3 -violating terms should raise little difficulty. Then, flavour constraints are relevant in the CP-violating NMSSM and we intend to design phenomenological tests accordingly. Finally, the dominant two-loop corrections to the Higgs masses will be calculated in a more quantitative way.

A Reference functions
We only consider the finite part of the loop integrals:

B The tree-level masses and couplings
This appendix provides the reader with a detailed presentation of the tree-level spectrum and couplings of the CP-violating, minimal-flavour-violating, R-parity and Z 3 conserving NMSSM.

B.1 Tree-level masses
Here we derive the tree-level bilinear terms of the lagrangian. For a later application to the Higgs couplings as well as to the loop-corrections in the Coleman-Weinberg effective potential, we will try to keep a full dependence in the Higgs scalar fields S, To evaluate the masses, one of course simply needs to replace these fields by their v.e.v.'s.

B.1.1 SM fermions
The Higgs-fermion potential reads: Focussing on the third generation (and neglecting off-diagonal CKM elements), we may cast under matrix form: from which we derive the squared-mass matrices: Replacing the Higgs fields by their v.e.v.'s, one obtains diagonal matrices M 2 q 3 and M 2 l 3 , with the usual relations:

B.1.2 Electroweak gauge bosons
From the Higgs kinetic terms, one obtains the Higgs-gauge potential (where we omit the derivative Higgs couplings): After the fields are rotated to the mass-states, we derive: This leads to the usual gauge-boson masses:

B.1.3 Sfermions
The Higgs-sfermion potential originates from soft, F and D terms: The bilinear sfermion terms can be cast under matrix form: with the matrix blocks: Moving to the v.e.v.'s, the matrices become block diagonal -each block being associated to a given electric charge of the sfermion fields. Under our Minimal Flavour Violation hypothesis the various generations also decouple so that we are left with 2 × 2 (hermitian) mass-matrices M 2 F . Those can be diagonalized via unitary matrices X F , according to: The mass-states are given by F i = X F iL F L + X F iR F c * R (where, in our notation 1 ↔ L and 2 ↔ R).

B.1.4 Charginos and neutralinos
The gaugino-higgsino-Higgs potential may also be cast under matrix form: The 2 × 2 chargino mass-matrix may be diagonalized via two unitary matrices U and V : , U and V , we consider the hermitian matrices: The choice of phasesφ U ,φ V is a priori arbitrary. We decide to determine them by the requirement that m χ ± 1 and m χ ± 2 , obtained in the matrix product U * M χ −+ V † , are real and positive. The associated mass-states are then: The 5 × 5 neutralino mass-matrix is symmetric, hence is diagonalizable via a single unitary matrix N : . . , 5)N . As before, we first consider the hermitian matrix

Gluinos
The gluons of course remain massless. Concerning their supersymmetric partners, the gluino bilinear terms read: so that we define the mass statesG a ≡ −ıe ı 2 φ M 3g a , with mass M 3 .

B.1.6 Higgs sector
The tree-level Higgs potential is given in Eq.9.

i) Minimization Conditions
First derivatives of the potential must vanish at the minimum, which provides: So that one can express certain parameters in terms of the v.e.v.'s: ii) Charged Higgs The 2 × 2 charged-Higgs bilinear terms read: Obviously, with the Goldstone boson G ± = − sin βH ± u + cos βH ± d and the charged Higgs state H ± = cos βH ± u + sin βH ± d . We will denote the corresponding rotation matrix as follows:

iv) Charged-Neutral Higgs terms
For completeness we indicate the bilinear terms mixing charged and neutral Higgs states (note that M 2

B.2 Tree-level Higgs couplings
Having presented the spectrum and our conventions, we may now turn to the Higgs couplings.

B.2.1 Higgs-SM fermions
Employing the Dirac-fermion notation, the Higgs couplings to SM fermions may be cast in the following form (with the usual left-and right-handed projectors P L,R ): with the (non-vanishing) values of g Sf f L :

B.2.2 Higgs-gauge
The situation is unchanged with respect to the CP-conserving case:

B.2.5 Higgs-to-Higgs couplings
From the tree-level potential of Eq.10, one may derive the trilinear and quartic Higgs couplings: where:

B.3.3 Chargino and Neutralino gauge couplings
Using the notation: the chargino and neutralino gauge couplings may be written: C Radiative corrections to the supersymmetric spectrum

C.1 Electroweak gauginos and higgsinos
We follow the approach of [20] and consider the loops involving sfermions / fermions, higgsinos or gauginos / Higgs and gauge bosons in the self energies of the gauginos and higgsinos, under the assumption that the gauge eigenstates are approximately mass states. Taking the complex phases φ Mi and ϕ λ,κ into account, we find the following corrections to the gaugino and higgsino masses: We took over the notations of [20] to designate the approximate masses of the particles in the loops; note that µ ≡ λs stands for the doublet higgsino mass, ms ≡ 2κs for the singlino one and m h 0 S ,a 0 S for the singlet (pseudo)scalar masses.

C.2 Sfermions
O(α S ) corrections to the squark masses are generated by gluon / squark and gluino / quark loops. Another source are squark self-couplings, as these receive a contribution from the SU (3) c D-term. We use the following expressions to correct the squark squared masses m 2 Q : We recover the results of [40] in the CP-conserving limit.

C.3 Gluino
We follow [20] to include the O(α S ) corrections to the gluino mass: these involve the gluon /gluino and the quark / squark loops. The latter depend on squark and gluino phases via the quark / squark / gluino couplings. We obtain:

D Radiative corrections to the Higgs spectrum D.1 Wave-function renormalization
We summarize the discussion of section 3.3.1. Remember that µ H = 125 GeV replaces the external momentum. The squared-bilinear matrices of the fermions of third generation have been provided in Eq.42, 43. One observes that they split into (at-most) 2 × 2 blocks corresponding to the left-handed quark fields, the right-handed ones, the left-handed leptons and the right-handed one. The following eigenvalues can be derived: Similarly, one can derive the corrections to the trilinear Higgs couplings:

D.2.3 Sfermions
Considering the sfermion mass-matrices of Eq.49, 50, 51, 52, 53, 54 and focussing on the neutral-Higgs dependence first, one obtains decoupling 2 × 2 blocks -1 × 1 in the case of the sneutrinos -so that eigenvalues may be expressed as m 2 Corresponding contributions to the neutral Higgs mass-matrix thus read (we denote as E ijk the coefficients coming from the tadpole equations -see Eq.32): with: Additionally, the sneutrino mass reads: . To derive the corrections to the charged-Higgs masses and Higgs-to-Higgs couplings, one is confronted to the task of diagonalizing the matrix-system of Eq.49, 50, 51, 52, 53, 54. This can be performed perturbatively, as an expansion in the Higgs-doublet fields, which amounts to a series in v MSUSY . We confine to a precision of order O v 2 at the level of the masses, which means that we compute the potential up to terms of H 4 -order (H standing for any Higgs-doublet field) and freeze singlet fields to their v.e.v. s for terms of H 4 -order (they are kept explicitly for terms of lower order in the expansion). The ensuing corrections to the Higgs potential can be matched onto Eq.33 and we may then use the results of section E.3, e.g. for the charged-Higgs mass (Eq.67) or the Higgs couplings. For squarks of each generation (note that we neglect the Yukawa couplings of the two first families): For the sleptons:

D.2.4 Charginos and neutralinos
To include the chargino and neutralino contributions to the effective Higgs potential, we turn exclusively to the method that we have just presented in the case of the sfermions. In other words, we diagonalize the matrix system of Eq.56 perturbatively, in an expansion of doublet Higgs-fields and we match the ensuing potential to the form of Eq.33. One can then work out the contributions to the Higgs mass-matrix and couplings. Note, however, that instead of diagonalizing directly the 9×9 (squared) bilinear matrix of Eq.56, it is easier to consider the dependence on neutral Higgs fields only (that is replacing charged fields by 0), as the corresponding matrix then splits into various blocks. All the couplings of Eq.33 can be identified from the neutral potential, with the exception of λ 3,4 , which only appear in terms of the sum λ 3 + λ 4 . It is a straightforward task, however, to compute λ 4 in a second step, from the charged couplings of the full potential. Another remark accompanies the observation that, as higgsino and singlino masses depend on the singlet Higgs field, the coefficients in our matching procedure depend on the singlet fields as well (those would correspond to operators of dimension > 4), which leads to additional (but straightforward) terms, with respect to the results of section E.3. This S-dependence can be neglected for terms of order H 4 , as keeping it would produce terms of higher order in v 2 M 2

D.2.5 Higgs-to-Higgs contributions
Instead of diagonalizing the Higgs bilinear terms, we compute the Higgs self-energy (Π) and tadpole (T ) diagrams mediated by gauge and Higgs particles in the Feynmann gauge. We then set the external momentum to zero to determine the potential contributions to the Higgs mass matrices and substract the pure gauge effects in the Landau gauge -from the results in appendix D.2.2. Finally we identify the corrections to the Z 3 -conserving parameters of Eq.33.
i) Pure-gauge contributions to the Higgs self-energies and tadpoles For external neutral Higgs: For an external charged Higgs (we consider only the physical state): ii) Higgs / gauge diagrams The neutral self energy receives contributions from hybrid Higgs (Goldstone) / vector diagrams: Π SV S 0 i S 0 j (p 2 ) = 1 32π 2 g 2 BSV (p, MW , MW ) + g 2 + g 2 2 BSV (p, MZ , MZ ) × sin 2 βX R iu X R ju + cos 2 βX R id X R jd + sin β cos β(X R iu X R jd + X R id X R ju ) + g 2 BSV (p, MW , MW ) sin 2 βX I iu X I ju + cos 2 βX I id X I jd − sin β cos β(X I iu X I jd + X I id X I ju ) + g 2 BSV (p, m H ± , MW ) cos 2 β(X R iu X R ju + X I iu X I ju ) + sin 2 β(X R id X R jd + X I id X I jd ) + sin β cos β(X R iu X R jd + X R id X R ju − X I iu X I jd − X I id X I ju ) For the charged Higgs self-energy: × cos 2 β(X R 2 ku + X I 2 ku ) + sin 2 β(X R 2 kd + X I 2 kd ) − sin β cos β(X R ku X R kd − X I ku X I kd ) iii) Pure Higgs loops The Higgs sfermion couplings are given in appendix B.2.3.

E Simplified effective potential
In this appendix, we study the simplified effective Higgs potential of Eq.33, or more precisely the following and slightly modified version: This simplified potential is meant as an expansion of the effective potential -see Eq.29 -up to quartic order in the doublet fields. It slightly differs from Eq.33 in that the Z 3 -symmetry has been explicitly restored in the terms of the last line. Note that this way of restoring the Z 3 -symmetry is just the simplest educated guess, while any additional factor f (|S| 2 , S 3 , S * 3 ) could intervene. Therefore, the factors of S/s appearing in the last line are just chosen as such because they will provide improved results numerically. Formally however, the associated corrections will remain of subleading order in the expansion in the doublet v.e.v.'s.

E.4 Trilinear Higgs couplings
We omit the quartic couplings.