On-Shell neutral Higgs bosons in the NMSSM with complex parameters

The Next-to-Minimal Supersymmetric Standard model (NMSSM) appears as an interesting candidate for the interpretation of the Higgs-measurement at the LHC and as a rich framework embedding physics beyond the Standard Model. We consider the renormalization of the Higgs sector of this model in its $\mathcal{CP}$-violating version, and propose a renormalization scheme for the calculation of on-shell Higgs masses. Moreover, the connection between the physical states and the tree-level ones is no longer trivial at the radiative level: a proper description of the corresponding transition thus proves necessary in order to calculate Higgs production and decays at a consistent loop order. After discussing these formal aspects, we compare the results of our mass calculation to the output of existing tools. We also study the relevance of the on-shell transition-matrix in the example of the $h_i \to \tau^+ \tau^-$ width. We find deviations between our full prescription and popular approximations that can exceed $10\%$.

The Next-to-Minimal Supersymmetric Standard model (NMSSM) appears as an interesting candidate for the interpretation of the Higgs-measurement at the LHC and as a rich framework embedding physics beyond the Standard Model. We consider the renormalization of the Higgs sector of this model in its CP-violating version, and propose a renormalization scheme for the calculation of on-shell Higgs masses. Moreover, the connection between the physical states and the tree-level ones is no longer trivial at the radiative level: a proper description of the corresponding transition thus proves necessary in order to calculate Higgs production and decays at a consistent loop order. After discussing these formal aspects, we compare the results of our mass calculation to the output of existing tools. We also study the relevance of the on-shell transition-matrix in the example of the h i → τ + τ − width. We find deviations between our full prescription and popular approximations that can exceed 10%.

Introduction
Since the discovery of a Higgs-like particle with a mass around 125 GeV by the ATLAS and CMS experiments [1,2] at CERN, a lot of effort has been invested to reveal its nature as the particle responsible for electroweak symmetry breaking. While within the present experimental uncertainties the properties of the observed state are compatible with the predictions of the Standard Model (SM) [3] many other interpretations are possible as well, in particular as a Higgs boson of an extended Higgs sector.
One of the prime candidates for physics beyond the SM is softly-broken supersymmetry (SUSY), which doubles the particle degrees of freedom by predicting two scalar partners for each SM fermion, as well as fermionic partners for all bosons-for reviews see [4,5]. The Next-to-Minimal Supersymmetric Standard Model (NMSSM) [6,7] is a well-motivated extension of the SM. In particular it provides a solution for the "µ problem" [8] of the Minimal Supersymmetric Standard Model (MSSM), by naturally relating the µ parameter to a dynamical scale of the Higgs potential [9,10].
In this work, we specialize in the Z 3 -conserving version of the NMSSM, characterized by a scaleinvariant superpotential. The main effort of our project consists in analyzing radiative corrections in the Higgs sector of the CP-violating NMSSM. To serve this purpose, we elaborated a FeynArts [104,105] model file and a set of Mathematica routines for the evaluation of the Higgs masses and wavefunction normalization matrix at full one-loop order and beyond. These should serve as a basis for a future inclusion of the CP-violating NMSSM in the FeynHiggs [106][107][108][109][110][111][112] package-originally designed for precise calculations of the masses, decays, and other properties of the Higgs bosons in the CP-conserving or -violating MSSM. A first step in this direction is represented by Ref. [113], centering on the CP-conserving NMSSM. In the current paper, we expand this project further. We follow the general methodology of FeynHiggs, relying on a Feynman-diagrammatic calculation of radiative corrections, which employs FeynArts [104,105], FormCalc [114] and LoopTools [114]. Our chosen renormalization scheme differs somewhat from earlier proposals [88,113,115]. In particular, in our renormalization scheme, the electromagnetic coupling e-which is related to the fine-structure constant α = e 2 /(4π)-is defined in terms of the Fermi constant G F measured in muon decays.
In section 2, we shall introduce relevant notations and describe the renormalization procedure underpinning our model file for the CP-violating NMSSM. In this section we also describe our implementation of higher-order corrections in the Higgs sector. A numerical evaluation of our results follows in section 3, where we will validate our calculation by a comparison with public codes. We will also insist on the relevance of the field-renormalization matrix for a consistent evaluation of the Higgs decays at the one-loop level, before a short conclusion in section 4.

Higgs masses and mixing in the CP-violating NMSSM
After a few general remarks concerning our notations and conventions, we present the renormalization conditions that we employ in our calculation. There, we focus on effects beyond the MSSM in the Higgs and higgsino sectors, since we otherwise align with the conventions of FeynHiggs, described in [109]. Then we discuss how to formally extract the loop-corrected Higgs masses and the wave function normalization factors.

Conventions and relations at the tree level
In the following, we consider the Z 3 -conserving version of the NMSSM and neglect flavor-mixing. The superpotential of the NMSSM (showing only one generation of fermions/sfermions) reads whereQ,û,d,L,ê,Ĥ 1 ,Ĥ 2 ,Ŝ denote the quark, lepton and Higgs superfields. The dot · stands for the SU (2) L -invariant product. The Yukawa couplings in Eq. (2.1) can be complex in general. However, their phases can be absorbed in a redefinition of the quark and lepton superfields. We may write the scalar fields inĤ 1 ,Ĥ 2 andŜ explicitly in terms of their (real and positive) vacuum expectation values (vevs), v 1 , v 2 and v s , respectively, as well as their CP-even, CP-odd, and charged components, φ i , χ i , and φ ± i , Here ξ 1 , ξ 2 and ξ s are the phases of the two Higgs doublets and the Higgs singlet, respectively. It is convenient to define the ratio tan β = v 2 /v 1 , the geometric mean of the doublet vevs v = v 2 1 + v 2 2 , as well as the sum ξ = ξ 1 + ξ 2 of the doublet phases. SinceŜ transforms as a singlet under the SM-gauge transformations, the D-terms of the scalar potential are unchanged with respect to the MSSM. On the other hand, as compared to the MSSM, additional dimensionless, complex parameters λ = |λ| e i φ λ and κ = |κ| e i φκ appear while the complex µ-term is absent. The latter is dynamically generated as an effective µ-term when the singlet field takes its vev, In the NMSSM the phases ξ and ξ s only appear in the combinations φ λ +ξ s +ξ and φ κ +3 ξ s , so that they could be absorbed in a re-definition of φ λ and φ κ . Nevertheless, we will keep the dependence on all phases of the Higgs sector explicitly, in order to allow for more flexibility on the choice of input.
Soft SUSY-breaking in the NMSSM is parametrized by the complex trilinear soft-breaking param- Expanding the Higgs potential in terms of the charged Higgs fields (φ (2.5) Here T = (T φ1 , T φ2 , T φs , T χ1 , T χ2 , T χs ) denotes the tadpole coefficients of the neutral Higgs fields, and M 2 and M 2 φ ± denote the mass matrices of the neutral and charged Higgs bosons, respectively. Since M 2 and M 2 φ ± are symmetric and hermitian matrices, respectively, we diagonalize them by an orthogonal (6 × 6) matrix U n and a unitary (2 × 2) matrix U c , respectively, These transformations define the five neutral Higgs boson mass eigenstates, h i , (i = 1, . . . , 5), and the (would-be) Goldstone boson G, as well as the charged Higgs and (would-be) Goldstone states, H ± and G ± , at the tree level, It is convenient to decompose U n into two matrices U G n and U 5 n , where U G n singularizes out the neutral Goldstone boson, In the CP-violating NMSSM the five fields h i are in general superpositions of the CP-even and -odd components φ i and χ j . In the special case of CP-conservation is restored in the Higgs sector at the tree level, and the neutral mass matrix M 2 becomes block-diagonal with two (3 × 3) sub-matrices for the CP-even and -odd entries. The five linearly independent tadpole coefficients are related to soft-breaking terms and combinations of phases as where the masses of the W and Z bosons are denoted by M W and M Z , respectively, and the phases combine to ζ 1 = ξ −2 ξ s +φ λ −φ κ , ζ 2 = ξ +ξ s +φ A λ +φ λ and ζ 3 = 3 ξ s +φ Aκ +φ κ . The expressions of Eq. (2.10) make plain that the tadpole coefficients can substitute the five parameters m 2 1 , m 2 2 , m 2 S , φ A λ and φ Aκ , so that the latter will not be regarded as free parameters in the following. Finally, the tadpole coefficients in the (tree-level) mass basis, T h = (T h1 , T h2 , T h3 , T h4 , T h5 , 0), where the zero denotes the vanishing tadpole coefficient of the Goldstone mode, are obtained by T h = U n T. The minimization of V H at the chosen Higgs vevs is guaranteed through the condition that all tadpole coefficients T vanish at the tree level.
The trilinear parameter |A λ | can be expressed in terms of the charged Higgs mass M H ± as Our renormalization scheme will also involve the fermionic superpartners of the Higgs bosons, known as the higgsinos. We thus introduce here the Dirac spinorsH ± of the charged higgsino fields, as well as the Majorana spinor of the singlinoS. In turn, these higgsino gauge eigenstates mix with the gauginos to form the mass states known as the neutralinos and charginos-see e. g. Eq. (11) in Ref. [113], where the NMSSM parameters λ, κ, M 1,2 and µ eff should be promoted to complex values. Yet these mass states will play no role in the discussion below.
In the present work, the radiative corrections to the Higgs sector are calculated in the diagrammatic approach. To this end, we first establish a list of the independent parameters appearing in the linear and bilinear terms of the Higgs potential in Eq. (2.5): where the electromagnetic coupling e is related to the fine structure constant α by e = √ 4 π α. Compared to Ref. [113], we use e as an independent parameter instead of the vev v. This choice allows to renormalize e to its value derived from the Fermi constant G F and does not require the reparametrization procedure employed in Ref. [113]. The difference between the renormalization employed here, and the renormalization and reparametrization procedure described in Ref. [113] is a sub-leading effect of two-loop order, however. Other proposals in the literature consist in fixing e from α (M Z ) [88,115].
To all real and complex independent parameters, g r and g c , respectively, that are given in Eq. (2.12), we apply the renormalization transformations The renormalization transformations for the Higgs, singlino and charged higgsino fields read H 1,2 → 1 + 1 2 δZ H1,2 H 1,2 , S → 1 + 1 2 δZ S S, (2.14a) with P L and P R denoting the left-and right-handed projectors, respectively. Since the singlino is a Majorana field, the corresponding wave-function counterterms δZ L S and δZ R S are complex conjugates of one another.

Renormalization conditions at the one-loop order
For the parameters M H ± , M W , M Z and tan β, which enter the one-loop calculation of the Higgs masses in the MSSM as well, we follow the renormalization prescription outlined in Ref. [109]: the on-shell renormalization scheme is employed for the gauge boson masses, M Z and M W , and the charged Higgs mass M H ± , while the parameter tan β is renormalized DR.
We apply the minimization conditions in order to fix the tadpole counterterms: where the T hi correspond to the one-loop contributions to the tadpole parameters. The counterterm of the electromagnetic coupling e is fixed by Here δZ Th e is the counterterm of the charge renormalization within the NMSSM according to the static (Thomson) limit. The quantities Π γγ (0) and Σ γZ T (0) are respectively the derivative of the transverse part of the photon self-energy and the transverse part of the photon-Z self-energy at zero momentum transfer. For the quantity ∆r NMSSM , relating the elementary charge to the Fermi constant G F measured in muon decays, we use the result of Ref. [127] (see also Ref. [128]). The numerical value for the electromagnetic coupling e in this parametrization is obtained from the Fermi constant in the usual way as e = 2 M W s w √ 2 G F . This choice differs from previous works, where either the charge renormalization condition was determined in terms of α(M Z ) [115], or instead v was renormalized DR and the result was subsequently reparametrized to use the value of e derived from the Fermi constant [113].
The remaining independent parameters and the field renormalization constants are renormalized DR. We present a detailed description of the DR renormalization conditions that we apply. The actual cancellation of UV-divergences, that we recover at the diagrammatic level, represents a non-trivial check for the validity of the FeynArts model-file employed for our calculation.
The DR field renormalization constants for the Higgs fields are obtained as where Σ (1) ii denotes the self-energy of field i at the one-loop order, and the subscript 'div' denotes the UV-divergent piece (along with the universal finite pieces that are associated in the DR scheme) of the quantity that it follows. The result does not depend on the momentum p 2 .
The field renormalization constants for the charged higgsino and the singlino fields are defined by the following conditions (the momentum p 2 again does not matter) where the self-energies of the fermion fields are decomposed into the left and right vector and the scalar contributions,

Σ
(1) The renormalization constants δ|λ|, δ|κ|, δφ λ , δφ κ , δξ, δξ s and δ|A κ | are fixed by DR conditions imposed on trilinear vertices involving scalar, CP-even Higgs fields φ 1,2,s , the singlinoS, and the charged higgsino fieldsH ± , in the interaction basis in analogy to the procedure outlined in [129]. The renormalization condition imposed on the renormalized three-point functionΓ ijk for three arbitrary fields i, j and k readŝ where Γ (0) ijk and Γ (1) ijk denote the vertex function at the tree-level and one-loop order, respectively, and δΓ ijk denotes the counterterm. The counterterm δΓ ijk is thus fixed by the divergent part of the vertex function. The renormalization constants of the independent parameters are subsequently fixed by linear relations to δΓ ijk .
• For δ|λ| and δφ λ we impose the DR renormalization condition of Eq. (2.20) on the vertices Γ • For δξ we impose the renormalization condition of Eq. (2.20) on the vertices Γ • We fix the renormalization constant δξ s for the phase ξ s by applying Eq. (2.20) on the vertices Γ (0) • The absolute value and phase of κ are renormalized by δ|κ| and δφ κ . We fix both renormalization constants by applying Eq. (2.20) on the vertex Γ • The parameter |µ eff | could be renormalized in the on-shell scheme for one of the charginos or neutralinos [130][131][132]. However, such schemes cannot be stabilized over the whole parameter space (due to mass-crossings). We thus prefer to apply the DR condition • In order to fix δ|A κ |, we impose the renormalization condition of Eq. (2.20) on the vertex where we used the one-loop relation δcos ζ 3 = − sin ζ 3 δζ 3 , and δv is not an independent counterterm, but a quantity depending on the counterterms to the electroweak parameters We performed various consistency checks of our model file at the one-loop order: • all the renormalized Higgs self-energies are UV-finite, for arbitrary values of the momentum, • all the vertex-diagram amplitudes of a Higgs state decaying to SM-particles or a pair of charginos/neutralinos are UV-finite, • the UV-divergences of the counterterms to gauge couplings, superpotential parameters or soft terms are consistent with the corresponding one-loop beta functions (see e. g. Refs. [6,133]), • in the CP-conserving limit, our parameters and couplings are identical to the findings of a previously developed model file [113], • in the MSSM limit, we have found agreement of the values of all our couplings with their counterparts in the complex MSSM, obtained with the model file of [134].
• we checked that φ λ + ξ + ξ s and φ κ + 3 ξ s were the only relevant combinations of the phases φ λ , φ κ , ξ and ξ s at the level of amplitudes, • finally, we checked explicitly, that the counterterms δφ λ , δφ κ , δξ and δξ s vanish when all NMSSM contributions are included, as pointed out in Ref. [115]. This can also be placed in the perspective of the β-functions [133,[135][136][137][138][139]: the phases from the superpotential parameters have no scale-dependence (at least up to two-loop order); since ξ and ξ s are spurious degrees of freedom, we could expect their counterterms to present the same vanishing behaviors as δφ λ and δφ κ .

Quark Yukawa couplings
The Yukawa couplings of the top and bottom quarks, Y t and Y b , have a sizable impact on radiative corrections to the Higgs masses. We present our prescriptions in this subsection. The top Yukawa coupling Y t = 2 √ 2 G F m t / sin β is defined by the on-shell top mass m t . For the bottom quark, we employ the running DR bottom-mass of the SM (containing one-loop QCD corrections), m b , at the scale m t [140]. Additionally, we subtract the possibly large tan βenhanced one-loop contributions to m b -induced by gaugino-squark and higgsino-squark loopsfrom the numerical definition of Y b at the tree level: Refs. [98,[140][141][142][143][144][145][146].

Higgs masses at higher orders
The masses of the Higgs bosons are obtained from the complex poles of the full propagator matrix. After rotating out the Goldstone mode 1 the inverse propagator matrix for the five Higgs fields h i is a (5 × 5) matrix that reads∆ −1 (2.28) } denotes the diagonalized mass matrix of the Higgs fields without the Goldstone at the tree level, andΣ hh denotes the matrix of the renormalized self-energy corrections of the neutral Higgs fields.
The five complex poles of the propagator are given by the values of the squared external momentum k 2 for which the determinant of the inverse propagator matrix vanishes, where we have explicitly stated the connection between the pole M 2 i , the Higgs mass M hi and the total width Γ hi for each Higgs field h i .
In order to account for the imaginary parts of the poles of the propagator matrix, we perform an expansion of the self-energies in terms of the imaginary part of the momentum, which is assumed to be small (also see section 4.3.5 of [147]), In this work, the renormalized self-energyΣ hh , is evaluated by taking into account the full contributions from the CP-violating NMSSM at one-loop order and, as an approximation, the MSSM-like contributions at two-loop order of O(α t α s ) [148] and O α 2 t [149,150] at vanishing external momentum as implemented in FeynHiggs. 2 We note that the two-loop O(α b α s ) contributions to the Higgs self-energies are not included in our calculation. Still, as we employ the running bottom mass in the definition of Y b enterinĝ Σ (1L) hh k 2 NMSSM , we expect that the missing two-loop piece is numerically subleading [151][152][153].

Wave function normalization factors: the matrix Z mix
In the Feynman-diagrammatic approach physical processes with external Higgs fields are defined in terms of the tree-level mass states h i . When higher-order contributions are considered, however, the tree-level mass states are not physical states. Indeed, radiative corrections induce additional mass and kinematic mixing among the fields h i , and the poles of the tree-level propagators do not coincide with M 2 i . A relation between the amplitudes with an external tree-level Higgs mass state and those with an external physical Higgs state is necessary (though this relation is trivial if the fields are renormalized on-shell). For example, for a Higgs decaying into two fermions f this relation is given by the LSZ reduction formula, Here the superscript 'phys' denotes the amplitude with an external physical field. The coefficients Z mix ij can be expressed explicitly in terms of the full propagator matrix (see Refs. [154][155][156][157] and also section 5.3 of Ref. [158]), as However, the analytical inversion becomes time-consuming in the case of a (5×5) propagator matrix. Additionally,∆ hh needs to be evaluated at (or close to) its singular points M 2 i , which can lead to numerical instabilities on the right-hand side of Eq. (2.33) (only the ratio of propagator matrix elements (∆ hh ) ij is finite). In order to avoid these issues we employed an equivalent formulation of the coefficients Z mix ij , which is outlined below. We consider the following effective Lagrangian for the tree-level mass states h (tree) and set Z mix = (Z mix ij ) as the transition matrix to the physical Higgs fields: In other terms, the h (loop) i should appear as on-shell fields with standard kinetic terms close to their mass-pole. Thus the coefficients Z mix ij should satisfy the 'eigenvalue' conditions 3 Once the roots of det ∆ −1 hh k 2 are known, the i-th line of Z mix is thus determined as the eigenvector of D hG −Σ hG for the eigenvalue M 2 hi . Finally, the normalization of the i-th line of Z mix is specified by the following condition on the kinetic term, such that the coefficients Z mix ij are uniquely specified (up to a sign without physical meaning) by Eq. (2.35). For the study of effects from the normalization of Z mix , it is convenient to define the (squared) norm |Z i | 2 of its rows, As we discussed, the components of the matrix Z mix establish the connection between the physical fields and the tree-level external legs. In the literature this matrix Z mix is often replaced by simplified versions neglecting the momentum dependence of the self-energies. With the aim of performing numerical comparisons in the following section, we introduce two such approximate definitions of the relation between loop-corrected and tree-level fields: • The first approach consists in freezing the momentum to k 2 = 0 in the self-energy of Eq. (2.28).
This assumption is known as the effective potential approximation. In this approach the inverse propagator matrix ∆ −1 hh (k 2 = 0), as given by Eq. (2.28), is diagonalized by a simple orthogonal matrix U 0 , which approximates Z mix . This procedure aims at more closely mimicking the actual values of the self-energies involved in the mass calculation. In this approach the inverse propagator is also diagonalized by an orthogonal matrix U m .
While these procedures capture the mixing effects induced by radiative corrections, at least partially, it is nevertheless obvious that they miss the normalization of the fields outlined in Eq. (2.37). This means in particular that the norm as defined in Eq. (2.38) will always be identical to 1 if the matrix Z mix is approximated by either U 0 or U m . We will discuss the impact of these approximations in the following section.

Numerical Analysis
In this section we present the results of our Higgs mass calculation and compare them with the output of public tools for several CP-violating scenarios. We also investigate the relevance of the matrix Z mix for transition amplitudes in the example of the one-loop corrected decays of one Higgs field into a tau/anti-tau pair, h i → τ + τ − . The choice for the top-quark mass is m t = 173.2 GeV. Throughout this section all DR parameters are defined at m t , and all stop-parameters are on-shell parameters.
From the point of view of the Higgs phenomenology we test the scenarios presented in this section with the full set of experimental constraints and signals implemented in the public tools HiggsBounds-4.3.1 [159][160][161][162][163][164] and HiggsSignals-1.3.1 [164,165].

Comparison with FeynHiggs in the MSSM-limit
In the limit of vanishing λ and κ, the singlet superfield decouples from the MSSM sector and one is left with an effective MSSM-the µ eff term persists as long as κ ∼ λ. We may then compare our results for the Higgs masses and the matrix Z mix in this limit to those of FeynHiggs-2.12.0. For this to be meaningful, we adjust the settings of FeynHiggs, so that they match the higher-order contributions and renormalization scheme of our NMSSM calculation. In particular, we impose that the one-loop field-renormalization constants and tan β are DR-renormalized and select full oneloop and leading two-loop MSSM contributions of O(α t α s + α 2 t ). We also require that FeynHiggs takes the tan β-enhanced contributions to the down-type Yukawa couplings into account. The corresponding FeynHiggs input flags read FHSetFlags[4,0,0,3,0,2,0,0,1,1].  We consider a region in the parameter space of the NMSSM with the following characteristics: λ = κ = 10 −5 , tan β = 10, m H ± = 500 GeV, µ eff = 250 GeV, A κ = −100 GeV; the sfermion soft masses are set to the universal value of 1.5 TeV and the sfermion trilinear couplings to a value of 0.5 TeV, with the exception of the third generation parameters |A t | = A b = 2.5 TeV; the gaugino masses are chosen as follows: 2M 1 = M 2 = M 3 /5 = 0.5 TeV. We then vary the phase φ At . A variation of φ At (or of any MSSM-like phase) in such a naive direction is of limited phenomenological interest, since in this case limits from Electric Dipole Moments (EDM) are violated almost as soon as CP, see e.g. Ref. [84]. In the following we dismiss this issue, however, and allow φ At to vary over its full range. Indeed, we are only interested in comparing our results with those of FeynHiggs. Due to the largely SM-like properties of the state with a mass close to 125 GeV, our scenario appears to retain characteristics that are compatible with the experimental data implemented in HiggsBounds and HiggsSignals, over the whole range of φ At . The results for the Higgs masses are displayed in Fig. 1. We observe a near perfect agreement between our results (solid curves) and those of FeynHiggs (squares) with differences of order MeV. This agreement is expected, since we closely follow the procedure for the renormalization and processing of the MSSM-like input of FeynHiggs. Moreover, due to the small values for λ and κ, deviations induced by genuine NMSSM effects remain negligible. The results for the elements of the matrix Z mix are displayed in Fig. 2. Again, we find a very good agreement between our results and FeynHiggs with differences below 1 for the modules.

Comparison in the CP-conserving limit
We now turn away from the MSSM limit. Our mass calculation can be confronted to the routines presented in [113] in the CP-conserving case. Both approaches employ an identical renormalization scheme in this limit, with the exception of the electroweak vev, which receives a DR renormalization in [113] while we parametrize v in terms of M W , M Z and e (see Eq. (2.27)). However, in [113] the input for v is obtained via a reparametrization from our scheme (the scheme using α(M Z ) as input is also considered), as explained in section 2.3 of that reference. Therefore, both mass predictions are directly comparable and the mismatch between them should be understood as an effect of two-loop electroweak order, due to the approximations in the reparametrization used by [113].
We consider the following region in the parameter space of the NMSSM: κ = λ/2, tan β = 10, M H ± = 1 TeV, µ eff = 125 GeV, A κ = −70 GeV; the soft masses are taken as in the previous subsection while the trilinear soft sfermion couplings are all set to 0.5 TeV, with the exception of A t = 1.2 TeV. We scan over λ from ∼ 0 (the MSSM-limit) to 0.5. The masses of the three lightest Higgs states are displayed in the plot on the left-hand side of Fig. 3. The dominantly SM-like state is the heaviest of the three (green curve). The two lighter states are dominantly singlet, CPeven (red curve) or CP-odd (blue curve). In the MSSM-limit, these three states are significantly lighter than 125 GeV: this results in an unsatisfactory phenomenological situation in view of the LHC measurements. With increasing λ, the CP-even singlet-doublet mixing uplifts the mass of the dominantly SM-like state, leading to phenomenologically viable characteristics-as tested with HiggsSignals-for λ ∼ 0.2. There, we observe that h 1 possesses a sizable doublet component and a mass M h1 100 GeV, so that this state could offer an interpretation of the local excess observed at LEP in e + e − → Z + (H → bb) searches [166].
We then focus on the comparison with the masses predicted by [113]. The plot on the left-hand side of Fig. 3 illustrates a general agreement between our calculation (solid curves) and the results of [113] (squares). On the right-hand side of Fig. 3, we display the mass differences between the two procedures, which are due to differences of two-loop order induced by the reparametrization used by [113]. We observe vanishing effects in the MSSM-limit while the mass differences eventually reach O(40 MeV) for λ 0.16. This can be understood in the following fashion: the leading effect originates in the Higgs mass matrix at the tree level, where an explicit dependence on v  On the left-hand side, we present our result (solid line) and that of [113] (squares). On the right-hand side, we plot the massdifferences between these two codes, due to the scheme applied to the electric coupling.
appears only 4 through terms of the form λ v and κ v (quadratically for the doublet and singlet mass entries, and linearly for the doublet-singlet mixing). These terms are processed differently in both approaches: in [113] v is regarded as an independent DR parameter, while in our calculation v is a dependent quantity that is expressed in terms of the independent parameters M W , M Z and e. While the reparametrization of [113] should restore the agreement between the two procedures, neglected effects of two-loop electroweak order in this reparametrization result in a small mismatch. Since the terms that convey this mismatch come with prefactors λ or κ, the difference vanishes in the MSSM limit (λ, κ → 0). Moreover, in the regime under consideration, where tan β 1, it is possible to understand why the mass of the CP-odd singlet (blue curve) is largely insensitive to the mismatch: terms ∝ (λ v) 2 in the CP-odd singlet mass entry are suppressed as 1/ tan β. Additionally, leading one-loop radiative corrections of O(α t ) induce further dependence on the processing of v. However, these corrections are suppressed for the points of Fig. 3, as the stops are relatively light.
On the whole, the numerical mismatch with the procedure of [113] is very minor, which places our current code in the direct continuity of this earlier work.

Comparison with NMSSMCALC
NMSSMCALC is particularly suitable for a comparison with our calculation, since its mixed DR/onshell renormalization scheme is relatively close to the one that we use. 5 Yet, we note several differences between the prescriptions implemented by NMSSMCALC and the procedure that we have outlined in section 2 (defining our "default" calculation). First, NMSSMCALC applies a renormalization scheme for the electric charge employing α(M Z ) as input-whereas we decided to define α via its relation to G F . Then the input parameters in the stop sector are defined in the DR scheme in NMSSMCALC-while we employ on-shell definitions. Additionally, we resum large tan β effects from our definition of the bottom Yukawa, contrarily to the Higgs-mass calculation of NMSSMCALC. Finally, NMSSMCALC includes only O(α t α s ) corrections at the two-loop order-where we consider 4 We remind the reader that both in [113] and in our calculation M 2 W and M 2 Z are chosen as independent, on-shell parameters. Therefore, the corresponding terms in the Higgs mass-matrix are not affected by the differences in the renormalization/reparametrization discussed here. 5 Note that it is somewhat more involved to compare our results quantitatively with RGE-based tools, as the input requires a conversion to the appropriate scheme (usually DR) and a running to the correct input scale [167]. For this reason, we shall confine our discussion to comparisons with NMSSMCALC, which shares closer characteristics with our approach. A similar comparison for real parameters has been presented in [168].
O α 2 t effects as well. However, the two-loop O(α t α s ) contributions of NMSSMCALC are exhaustive in the NMSSM (including corrections for the self-energies with at least one external singlet field)whereas ours are obtained in the MSSM approximation.
These observations mean that our mass-calculation is not directly (at least, not quantitatively) comparable to the predictions of NMSSMCALC, since, of the items listed above, the first few produce a deviation relative to the scheme, while the later ones generate a mismatch of higher orders. Consequently, several adjustments need to be performed in order to make a comparison meaningful and control the sources of deviations. Thus, NMSSMCALC has been adjusted in view of accepting on-shell input in the stop sector. 6 Moreover, we also establish a "modified" version of our routines that attempts to mimic the choices of NMSSMCALC-i. e. employing α(M Z ), discarding large-tan β effects for Y b and subtracting O α 2 t corrections-although we cannot currently include O(α t α s ) corrections beyond the MSSM, so that this effect should control the difference of our modified version with NMSSMCALC. Beyond this comparison with NMSSMCALC, we will also try to quantify the magnitude of the other higher-order effects that distinguish our "default" result from NMSSMCALC.
First, we consider the regime of the NMSSM with low tan β and large λ. This region in parameter space is well-known for maximizing the specific NMSSM tree-level contributions to the mass of the SM-like Higgs state as well as for stimulating singlet-doublet mixing effects and other genuine aspects of the NMSSM phenomenology. We employ the following parameters: λ = 0.7, |κ| = 0.1, tan β = 2, M H ± = 1170 GeV, µ eff = 500 GeV, A κ = −70 GeV; the soft masses are taken as in Fig. 1 with the exception of the squarks of the third generation, for which the soft masses and trilinear couplings are set to 500 GeV and 100 GeV, respectively. In the regime under consideration genuine NMSSM effects are indeed sufficient to produce a SM-like state in the observed mass-range without requiring large top/stop corrections. We vary the phase φ κ (we restrict to a range where the treelevel squared Higgs masses remain positive). We note that, contrarily to MSSM-like phases, the phases from the singlet sector are allowed a wide range of variation without conflicting with the measured EDM [83,169].
The results for the mass prediction are presented in Fig. 4. At vanishing φ κ the mass of the SM-like state (in blue) is somewhat low, m h2 ∼ 120 GeV, so that this point in parameter space has a very marginal agreement with the observed characteristics of the Higgs state. For non-vanishing φ κ , however, a CP-violating mixing with the lighter pseudoscalar singlet (in red) develops: this effect increases the mass of the light mostly CP-even state h 2 but affects its otherwise SM-like properties only in a subleading way. Consequently, we recover an excellent agreement with the LHC resultsas tested by HiggsSignals and HiggsBounds-for e. g. φ κ −0.11. Additionally, the dominantly CP-odd singlet h 1 then has a mass close to 100 GeV. As it acquires a doublet CP-even component via mixing, it could explain the LEP local excess in bb final states [166]. The mostly CP-even singlet h 3 (in green), with mass at ∼ 210 GeV plays no significant role. The masses of the heavier doublet-like fields h 4 and h 5 are approximately constant and close to M H ± .
In Fig. 4, we observe a good agreement between our results (solid lines), computed as described in section 2, and the predictions of NMSSMCALC (squares), although the corresponding masses are defined in different schemes and at different orders. For a more quantitative comparison, we turn to our "modified" scheme for the mass calculation. On the left-hand side of Fig. 5, we plot the deviation between the corresponding results and the predictions of NMSSMCALC for the three lightest Higgs states. We checked that the one-loop results are virtually identical, so that the differences between NMSSMCALC and our calculation are entirely controlled by two-loop effects. We observe typical deviations of order 0.5-1 GeV that should be interpreted as the impact of O(α t α s ) corrections beyond the MSSM-approximation. As could be expected, the masses of the mostly singlet states (red and green lines) tend to exhibit the largest effect, though the mass-predictions for the SM-  blue is associated to the essentially SM-like state h2 and green corresponds to the mostly CP-even singletlike state h3. We display our "default" result for the masses (solid curves) as well as the predictions of NMSSMCALC (squares).
h i correspond to the mass-differences (for each of the three lightest Higgs states) between the predictions of our "modified" scheme and NMSSMCALC: O(αtαs) should dominate these deviations. On the right-hand side, ∆M h i corresponds to the mass-shifts associated to O α 2 t contributions, which are calculated in our "default" scheme. The color of the curves match the convention of Fig. 4. like state may still differ by ∼ 0.5 GeV (for φ κ 0). The plot on the right-hand side of Fig. 5 depicts the magnitude of O α 2 t effects, which is quantified in our "default" scheme. Here again, the typical impact on the masses is of order 1 GeV. Expectedly, the masses of the almost pure singlet states (red curve at φ κ 0 or green curve) are insensitive to the corrections implemented in the MSSM-approximation. The mass of the mostly CP-odd singlet (red curve) is only affected when the corresponding state acquires a non-vanishing doublet component (φ κ = 0).
In Fig. 6, we compare the U 0 matrix elements that are delivered by NMSSMCALC (squares) with ours (solid line; NMSSMCALC does not provide Z mix ). The results show a satisfactory agreement also at this level.
Subsequently, we present our results in another region of the parameter space: λ = 0.2, |κ| = 0.6, tan β = 25, m H ± = 1 TeV, µ eff = 200 GeV, A κ = −750 GeV, the gaugino soft masses as well as the soft masses for the sfermions of first and second generations are chosen as before; for the third generation, the soft sfermion mass is set to 1.1 TeV; the trilinear soft terms are set to −2 TeV. With this choice of parameters, the singlet-like CP-even state and the heavy CP-even and CP-odd doublet-like states receive comparable masses of the order of 1 TeV. This results in a sizable mixing for the corresponding fields h 2 , h 3 and h 4 , which includes both singlet-doublet admixture as well as CP-violation (for non-vanishing φ κ ). The SM-like Higgs state has a mass close to ∼ 124 GeV on the whole range of φ κ , which leads to a good agreement with the Higgs properties measured at the LHC (as tested with HiggsSignals). The heaviest state h 5 has a mass of ∼ 1.2 TeV, which we will not comment further below.
On the left-hand side of Fig. 7, we show our prediction for the mass of the lightest (SM-like) Higgs state (full curve). The mass delivered by NMSSMCALC is represented by the squares at about 118 GeV, which is substantially smaller than ours (by ≈ 6 GeV). If we mimic the settings of NMSSMCALC (our "modified result", dotted curve), this discrepancy is considerably reduced. In fact, the difference between our full result and NMSSMCALC's is largely driven by the O α 2 t two-loop contributions, missing in NMSSMCALC. Again, both results are virtually identical at the one-loop order.
On the right-hand side of Fig. 7, we turn to the heavier states h 2 , h 3 and h 4 of this scenario. Our default results (full curves) are compatible with the predictions of NMSSMCALC (squares). The discrepancies are of order 1-3 GeV only, which should be considered from both the perspective of the different renormalization scheme of the electric coupling e and the different two-loop contributions. Actually, the mass predictions match almost exactly when comparing NMSSMCALC with our modified scheme (dotted curves). The corresponding deviations are shown on the left-hand side of Fig. 8 and fall in the range of 100 MeV. In this precise case, the difference between our results and NMSSMCALC is essentially driven by the resummation of large-tan β effects in the b-quark Yukawa coupling. On the right-hand side of Fig. 8, we quantify the associated mass-shift and find an impact of a few GeV.
Finally, we turn to the U 0 matrix elements for h 2 , h 3 and h 4 in Fig. 9. There, we observe sizable deviations between our default result (solid curves) and NMSSMCALC (squares), which, however, have no deep-reason to agree in view of the diverging options. If we keep in mind that the main difference between our full scheme and NMSSMCALC is controlled by the large-tan β corrections to Y b in this precise scenario, it is not surprising to observe large shifts, as mixing angles are indeed very sensitive to small deviations in the mass-matrix for states that are very close in mass. These differences largely vanish when we identify the output of NMSSMCALC with our modified results (dotted lines), which is better equipped for comparisons with this code.
In summary, the results of our mass calculation are largely compatible with the predictions of NMSSMCALC. Deviations with a magnitude of O(GeV) are indeed within the expected range if we allow for the different renormalization scheme of e and higher-order contributions. Such discrepancies tend to be reduced sizable when we modify our routines to adopt the assumptions of NMSSMCALC. The impact of two-loop O(α t α s ) contributions beyond the MSSM, two-loop O α 2 t and the resummation of large-tan β effects in the Yukawa couplings can be clearly identified from this comparison. As a final remark, let us emphasize that the uncertainty on the Higgs-mass calculation from unknown higher-order contributions and parametric errors may reach several GeV [108,167,168].

The matrix Z mix and the
The matrix Z mix is not an observable quantity in itself. It is a renormalization-scheme dependent object relating the tree-level mass states of the Higgs sector to the physical Higgs fields. For on-shell renormalized fields Z mix is trivial. In any other renormalization scheme, however, it is mandatory to include this transition to the physical fields for a proper description of external legs in Feynman diagrams at higher orders.
our modified two loop minus NMSSMCalc two loop two-loop shift by resummed m b Figure 8. Impact of the higher-order effects in the mass-predictions for h2, h3 and h4 in the scenario of Fig. 7. On the left-hand side, we show the deviation in mass between NMSSMCALC and our modified scheme (controlled by O(αtαs) corrections beyond the MSSM-approximation). The plot on the right-hand side illustrates the impact of large-tan β effects in Y b : the considered mass-difference is that induced in our default scheme by the ∆ b term. The colors follow the conventions of Fig. 7. The discontinuities at φκ = ± π 2 originate from the mass calculation in NMSSMCALC.
our default two loop our modified two loop NMSSMCalc two loop   A remarkable aspect of Z mix is that the eigenvectors that it contains do not preserve unitarity with respect to the tree-level fields. Instead, they satisfy the normalization condition given in Eq. (2.37). This is a feature that the approximations U 0 and U m are unable to capture (by construction). In a first step, we will show that the norms in Z mix can differ from 1 by a few percent in the scheme that we have described in section 2. Beyond the normalization of the fields, U 0 and U m also differ from Z mix in that they diagonalize the mass-matrix away from the poles of the propagator.
However, as we wrote above, Z mix is a scheme-dependent object and we should not pay excessive attention to its actual structure. In order to characterize its role in an observable quantity, we will consider the h i → τ + τ − decays at the one-loop level. We have chosen this particular channel as it is one of the main fermionic Higgs decays and proves technically simple to implement in a predictive way. Moreover, one-loop corrections are of purely electroweak nature-QCD contributions occur only at three-loop order and beyond-so that radiative corrections are expected to be moderate. This allows for a clean appreciation-free of large higher-order uncertainties-of the impact of the wave-function normalization matrix Z mix . Radiative corrections 7 are computed with our model file, except for the QED contributions, which are included according to the prescriptions of Refs. [170,171]. There, Z mix intervenes in the decay amplitudes of the physical fields according to Eq. (2.32) (we dismiss the superscript 'phys' throughout this section). We will show that the substitution of Z mix by the approximations U 0 and U m may lead to sizable deviations in certain regions of the the NMSSM parameter space. This result confirms the outcome of similar studies in the MSSM [109].
We turn to the following NMSSM input: the parameters are chosen as in Fig. 4, except for M H ± = 1.2 TeV and A κ = −100 GeV. We have plotted the masses of the three lighter Higgs fields as a function of φ κ in the plot on the left-hand side of Fig. 10. For vanishing φ κ the lightest Higgs state h 1 is SM-like and we checked with HiggsBounds and HiggsSignals that this point is consistent with the experimental data. The dominantly CP-odd, singlet-like state h 2 is only slightly heavier than the state h 1 in this case. For increasing values of φ κ the mixing of the states h 1 and h 2 tends to lower the mass of the SM-like state h 1 , which eventually becomes too light to accommodate the experimental data. The dominantly CP-even, singlet-like state h 3 has a near constant mass of ∼ 210 GeV for all depicted values of φ κ . The two heavier, CP-even and CP-odd doublet-like states have masses close to ∼ 1.2 TeV.
The results for the squared norms |Z i | 2 of the eigenvectors-see Eq. (2.38)-in this scenario are shown in the plot on the right-hand side of Fig. 10. We observe a departure from the value 1-which would correspond to a unitary transition, as modeled by the approximations U 0 and U m -by a few percent. The local extrema at φ κ 0 for |Z 1 | 2 (red curve) and |Z 2 | 2 (blue curve) are associated to the sudden disappearance of the mixing between the light CP-odd singlet and the SM-like states at φ κ = 0 (CP-conserving limit). The discontinuities of |Z 2 | 2 and |Z 3 | 2 (green curve) at φ κ ±0.5 and ±0.7 correspond to the crossing of decay thresholds (h 2 → W + W − , h 2 → 2 Z, h 3 → h 1 h 2 ). These "spikes" are associated to the singularities of the first derivatives of the loop functions involved in the determination of Z mix -the apparent singularities actually come with a finite height due to the imaginary parts of the poles. A proper description of these threshold regions would require that the interactions among the daughter particles (of the decays at threshold) are properly taken into account, which would result in e. g. interactions between the Higgs state and bound-states or s-waves of the daughter particles. This, however goes beyond the scope of the present work.
We now turn to the decay widths Γ(h i → τ + τ − ) in the scenario of Fig. 10. The widths are displayed in the left column of Fig. 11 in the exhaustive description of the Higgs external leg (i. e. employing Z mix ; solid curves), in the U m approximation (dashed lines) and in the U 0 approximation (dotted lines), for the five Higgs mass-eigenstates. We observe a sharp variation close to φ κ = 0 for the decays of h 1 and h 2 , both in the full and approximate descriptions. It is associated to the mixing that develops between the SM-like state h 1 and the dominantly CP-odd, singlet-like state h 2 : this effect transfers part of the doublet component of h 1 to h 2 , so that the second state acquires a nonvanishing coupling to SM fermions at the expense of the first. The sum of the decay widths for both these states remains approximately constant in the vicinity of φ κ = 0. The width Γ(h 3 → τ + τ − ) appears to be an order of magnitude smaller than the corresponding widths for h 1 and h 2 , an effect that is associated to the dominantly CP-even, singlet-like nature of h 3 . Still, Γ(h 3 → τ + τ − ) nearly doubles in the considered interval of φ κ , while the mass of h 3 is fairly stable: we can understand this fact in terms of the acquisition of a larger doublet component, which is channeled by the increased proximity of the masses of h 2 and h 3 . The widths of h 4 and h 5 are essentially constant with only small relative changes. The general φ κ -dependency of the approximated and the full results are very similar. Yet, a systematic shift can be observed, especially in the case of U 0 . This is consistent with the findings of similar studies in the context of the MSSM [109].
On the right-hand side of Fig. 11, we show the difference between the full and the approximate results, ∆Γ = Γ − Γ appr , normalized to the more accurate one obtained with Z mix . When Z mix is approximated by U 0 (dotted lines), the typical discrepancy averages 4%, although the deviation reaches beyond 10% in the case of the decays of the two lightest Higgs states in the vicinity of, but not at, φ κ 0. We stress that this interval close to φ κ = 0 corresponds to the phenomenologically relevant region from the perspective of the measured Higgs properties. The approximation of Z mix by U m tends to provide better estimates of the full result, though deviations reach up to ∼ 7%. For both approximations the largest deviations from the more complete result employing Z mix are intimately related to the proximity in mass of the SM-like and dominantly CP-odd, singlet-like states: as the approximations capture the dependence on the external momentum either partially (U m ) or not at all (U 0 ), the gap between the diagonal elements of the Higgs mass-matrix, hence the mixing between the two states, is not quantified properly. While this precise configuration might appear somewhat anecdotal, we wish to point out the popularity of NMSSM scenarios with a sizable singlet-doublet mixing. Dismissing this extreme case, the approximations of Z mix by U 0 , and to a lesser extent by U m , still generate errors of the order of a few percent at the level of the decay widths. In view of the precision of the measurements achievable at the LHC [3,[172][173][174][175], such discrepancies may appear of secondary importance. In the long run, however, if the Higgs couplings are studied more closely, e. g. at a linear collider [176][177][178][179], one would have to try and keep such sources of error to a minimum. mixing via Z mix mixing via U m mixing via U 0 Figure 11. In the left column, we show the decay widths Γ(hi → τ + τ − ) in the scenario of Fig. 10 for the five neutral Higgs states. The widths are computed at the one-loop level, and the mixing of the external, physical Higgs fields is expressed in terms of the matrix Z mix (solid), or approximated by the matrices U m (dashed) or U 0 (dotted). In the right column, the differences ∆Γ = Γ − Γappr between the widths obtained with Z mix and its approximate treatments U m (dashed) and U 0 (dotted) are depicted, normalized to the width Γ obtained with Z mix .

Conclusions
In this paper, we have discussed the renormalization of the NMSSM Higgs sector, including complex parameters. Radiative contributions to the Higgs self-energies have been included up to the leading two-loop MSSM-like effects of O α t α s + α 2 t . Beyond the calculation of on-shell Higgs masses in this scheme, we were interested in determining the transition matrix Z mix between the mass-and tree-level states. The latter plays an essential role in the proper description of external Higgs legs in physical processes at the radiative level.
Our predictions for the Higgs masses have been compared to the calculations of existing tools in several NMSSM scenarios. In the MSSM-limit of the model, we have recovered an excellent agreement with FeynHiggs. For non-vanishing λ and κ, we first compared our Higgs-mass prediction with the findings of a previous extension of FeynHiggs to the NMSSM in the case of real parameters, and found nearly identical results. Second, we compared our calculation in the case of complex parameters with NMSSMCalc and found values of the Higgs masses which are compatible, although small differences emerge as a result of different processing of the two-loop pieces, both for low and high tan β.
Finally, we investigated the impact of the transition matrix Z mix on the h i → τ + τ − width in a scenario with low tan β and large λ, where the SM-like Higgs state may have a sizable mixing with the CP-odd singlet. We compared the full one-loop calculation of the width-i. e. including Z mix -with the popular approximations U 0 and U m -which are determined for fixed, unphysical momenta. We found typical deviations at the percent level, although larger effects can develop in the presence of almost-degenerate states, especially in the U 0 approximation. Such precision effects will matter when the measurement of fermionic Higgs couplings reaches comparable accuracy.
In its current form, our mass-computing tool is contained within a Mathematica package. In time, our routines should be incorporated in an extension of FeynHiggs to the NMSSM.