The Order $\mathcal{O}(\alpha_t\alpha_s)$ Corrections to the Trilinear Higgs Self-Couplings in the Complex NMSSM

A consistent interpretation of the Higgs data requires the same precision in the Higgs boson masses and in the trilinear Higgs self-couplings, which are related through their common origin from the Higgs potential. In this work we provide the two-loop corrections at order ${\cal O}(\alpha_t \alpha_s)$ in the approximation of vanishing external momenta to the trilinear Higgs self-couplings in the CP-violating Next-to-Minimal Supersymmetric extension of the Standard Model (NMSSM). In the top/stop sector two different renormalization schemes have been implemented, the OS and the $\overline{\text{DR}}$ scheme. The two-loop corrections to the self-couplings are of the order of 10\% in the investigated scenarios. The theoretical error, estimated both from the variation of the renormalization scale and from the change of the top/stop sector renormalization scheme, has been shown to be reduced due to the inclusion of the two-loop corrections.


Introduction
While the discovery of the Higgs boson by the LHC experiments ATLAS [1] and CMS [2] certainly marked a milestone for particle physics, it also triggered a change of paradigm: The Higgs particle, formerly target of experimental research, has become a tool in the quest for our understanding of nature. Although the Standard Model (SM) of particle physics has been tested at the quantum level and the discovered scalar particle behaves SM-like [3] there are experimental and theoretical arguments to assume it to be a low-energy effective theory of a more fundamental theory appearing at some high scale. In the absence of any direct observation of new states the study of the Higgs boson and its properties may reveal the existence of beyond the SM (BSM) physics. In particular, the discovered particle could be the SM-like Higgs boson of the enlarged Higgs sector of a supersymmetric extension of the SM. Supersymmetric (SUSY) theories [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] require the introduction of at least two complex Higgs doublets in order to give masses to up-and down-type quarks and ensure an anomaly-free theory. This minimal setup is extended by a complex singlet superfield in the Next-to-Minimal Supersymmetric extension of the SM (NMSSM) [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. After electroweak symmetry breaking (EWSB) the NMSSM Higgs sector features seven Higgs bosons, which in the CP-conserving case are three neutral CP-even, two neutral CP-odd and two charged Higgs bosons. In contrast to the Minimal Supersymmetric extension (MSSM) [35][36][37][38] in the NMSSM CP-violation can occur in the Higgs sector already at tree level. The additional sources of CP-violation in SUSY theories are interesting not only because they clearly mark physics beyond the SM, but also because CP-violation is an important ingredient for successful baryogenesis [39]. From a phenomenological point of view it entails a plethora of interesting new physics (NP) scenarios not excluded by experiment yet.
In order to study NP extensions, to properly interpret the experimental data and to be able to distinguish different BSM realizations, from the theory side we need as precise predictions as possible not only for experimental observables 1 but also for the parameters of the theory under investigation. In the Higgs sector these are in particular the Higgs boson masses and couplings. In the recent years there has been quite some progress in the computation of the higher order corrections to the Higgs boson masses of both the CP-conserving and CP-violating NMSSM. Thus in the CP-conserving NMSSM after the computation of the leading one-loop (s)top and (s)bottom contributions [41][42][43][44][45] and the chargino, neutralino as well as scalar oneloop contributions at leading logarithmic accuracy [46], the full one-loop contributions in the DR renormalization scheme have first been provided in [47] and subsequently in [48]. In [47] also the order O(α t α s + α b α s ) corrections in the approximation of zero external momenta have been given, and recently, first corrections beyond order O(α t α s + α b α s ) have been published in [49,50]. Our group has calculated the full one-loop corrections in the Feynman diagrammatic approach in a mixed DR-on-shell and in a pure on-shell renormalization scheme [51]. In the CPviolating NMSSM the contributions to the mass corrections from the third generation squark sector, from the charged particle loops and from gauge boson contributions have been computed in the effective potential approach at one loop-level in Refs. [52][53][54][55][56]. The full one-loop and logarithmically enhanced two-loop effects in the renormalization group approach have subsequently been given [57]. We have contributed with the calculation of the full one-loop corrections in the Feynman diagrammatic approach [58] and recently provided the two-loop corrections to the neutral NMSSM Higgs boson masses in the Feynman diagrammatic approach for zero external momenta at the order O(α t α s ) based on a mixed DR-on-shell renormalization scheme [59].
Several codes have been published for the evaluation of the NMSSM mass spectrum from a user-defined input at a user-defined scale. The Fortran package NMSSMTools [60][61][62] computes the masses and decay widths in the CP-conserving Z 3 -invariant NMSSM and can be interfaced with SOFTSUSY [63,64], which provides the mass spectrum for a CP-conserving NMSSM, also including the possibility of Z 3 violation. Recently, it has been extended to include also the CP-violating NMSSM [65]. The spectrum of different SUSY models, including the NMSSM, can be generated by interfacing SPheno [66,67] with SARAH [49,[68][69][70][71]. This is also the case for the recently published package FlexibleSUSY [72,73], when interfaced with SARAH. All these codes include the Higgs mass corrections up to two-loop order, obtained in the effective potential approach. The program package NMSSMCALC [74,75] on the other hand, which calculates the NMSSM Higgs masses and decay widths in the CP-conserving and CP-violating NMSSM, provides the one-loop corrections and the O(α t α s ) corrections in the full Feynman diagrammatic approach, where the latter are obtained in the approximation of vanishing external momenta.
The Higgs self-couplings are intimately related to the Higgs boson masses via the Higgs potential. For a consistent description therefore not only the Higgs boson masses have to be provided at highest possible precision, but also the Higgs self-couplings need to be evaluated at the same level of accuracy. The trilinear Higgs self-coupling enters the Higgs-to-Higgs decay widths. These can become sizable in NMSSM Higgs sectors with light Higgs states in the spectrum [76][77][78], and via the total width these decays sensitively alter the branching ratios of these states. Also Higgs pair production processes are affected by the size of the trilinear Higgs self-couplings [79][80][81]. Their determination marks a further step in our understanding of the Higgs sector of EWSB [82][83][84]. We have provided the one-loop corrections to the trilinear Higgs self-couplings for the CP-conserving NMSSM [79]. They have been calculated in the Feynman diagrammatic approach with non-vanishing external momenta. The renormalization scheme that has been applied is a mixture of On-Shell (OS) and DR conditions. In this paper we present, in the framework of the CP-violating NMSSM, our computation of the dominant two-loop corrections due to top/stop loops to the trilinear Higgs self-couplings of the neutral NMSSM Higgs bosons. In addition, we give explicit formulae for the leading one-loop corrections at order O(α t ). We use the Feynman diagrammatic approach in the approximation of zero external momenta and furthermore work in the gaugeless limit. We find that the determination of the two-loop corrections reduces the error on the trilinear Higgs self-coupling due to unknown higher order corrections and hence contributes to the effort of providing precise predictions for NMSSM parameters and hence observables. We have furthermore expanded for this paper the full one-loop corrections with full momentum dependence to include CP-violating effects.
The outline of the paper is as follows. In section 2 we set our notation, introduce the NMSSM Higgs sector and present the determination of the loop-corrected effective trilinear Higgs selfcouplings. Section 3 is then dedicated to the numerical analysis. We discuss the effects of the loop corrections on the trilinear Higgs self-couplings and the implications for Higgs-to-Higgs decays. Section 4 contains our conclusions.

The effective trilinear Higgs self-couplings in the NMSSM
In this section we present the details of the calculation of the effective trilinear Higgs selfcouplings at order O(α t ) and at order O(α t α s ). We closely follow the convention and notation of our paper on the Higgs mass corrections in the complex NMSSM at order O(α t α s ) [59]. We therefore repeat here only the most important definitions relevant for our calculation. We work in the framework of the complex NMSSM with a scale invariant superpotential and a discrete Z 3 symmetry. In terms of the two Higgs doublets H d and H u , and the scalar singlet S, the Higgs potential reads, The indices of the fundamental representation of SU (2) L are denoted by i, j = 1, 2, and ij is the totally antisymmetric tensor with 12 = 12 = 1. The dimensionless parameters λ and κ and the soft SUSY breaking trilinear couplings A λ and A κ can in general be complex. The U (1) Y and SU (2) L gauge couplings are given by g 1 and g 2 , respectively. In order to obtain the Higgs boson masses, trilinear and quartic Higgs self-couplings from the Higgs potential, the Higgs doublets and the singlet field are replaced by the expansions about their vacuum expectation values (VEVs), v d , v u and v s , where two additional phases, ϕ u and ϕ s , have been introduced. Note, that in order to keep the Yukawa coupling neutral we absorb the phase ϕ u into the left-and right-handed top fields, which of course affects all couplings involving only one top quark [59].
We work in the approximation of zero external momenta and call the thus derived selfcouplings 'effective' self-couplings. The (loop-corrected) self-couplings are automatically real in this approach. In the interaction basis, the effective trilinear Higgs self-couplings at order O(α t α s ) can be cast into the form with φ = (h d , h u , h s , a d , a u , a s ) and i, j, k = 1, . . . , 6. The first term represents the tree-level trilinear couplings, which can directly be derived from the tree-level Higgs potential Eq. (2.1) by taking the derivative Explicit expressions for these couplings can be found in Appendix A. The second and third terms denote the one-and two-loop corrections to the Higgs self-couplings. They can be obtained by either taking the derivative of the corresponding loop-corrected effective potential or by using the Feynman diagrammatic approach in the approximation of zero external momenta. At one-loop level we use both methods and find that the results obtained in these two different approaches agree as expected. However, at two-loop level, for the sake of automatization of our codes we solely employ the Feynman diagrammatic approach. 2 Therefore only the latter is described in the following.
In order to obtain the effective trilinear couplings in the mass eigenstate basis, the selfcouplings in the interaction basis have to be rotated to the mass basis by applying the rotation matrix R (l) . In detail, we have, where l = 1, 2 stands for the loop order and Φ for the loop-corrected mass eigenstates. These are denoted by upper case H and ordered by ascending mass with The neutral Goldstone boson G has been singled out. Note in particular, that the mass eigenstates are no CP eigenstates any more since we work in the CP-violating NMSSM. In order to be as precise as possible in the computation of the loop-corrected trilinear Higgs self-couplings in the mass eigenstate basis, we employ the most precise rotation matrix R (l) that is available. This means, that we rotate to the mass eigenstates H i (i = 1, ..., 5) at two-loop order. The loopcorrected rotation matrix R (l) is computed by the Fortran package NMSSMCALC [74,75] where the zero momentum approximation is employed, so that the matrix is unitary. In particular, the rotation matrix includes the complete electroweak (EW) corrections at one-loop order and the order O(α t α s ) corrections at two-loop level. For more details see [51,58,59].
The rotation matrix R (l) can be decomposed in the rotation matrix R, that rotates the interaction eigenstates to the tree-level mass eigenstates Φ (0) ≡ (h 1 , h 2 , h 3 , h 4 , h 5 , G), singling out the Goldstone boson G, and in the finite wave-function renormalization factor Z, cf. [51,79], With this definition, the loop-corrected effective trilinear couplings between the Higgs bosons in the 2-loop mass eigenstate basis are hence given by (i, i , j, j , k, k = 1, ..., 5) where the couplings in the tree-level mass eigenstates Γ h i h j h k are obtained from Eq. (2.3) by rotation with R.

The order O(α t ) corrections
In this subsection we present the one-loop corrections at order O(α t ). Due to the large top quark Yukawa coupling, at one-loop level the corrections from the top/stop sector are the dominant corrections to the Higgs boson masses and self-couplings. This is in particular true for the SMlike Higgs boson. The latter must be dominantly h u -like, inducing via the top loop a sufficiently large coupling to the gluons, so that its rates are in accordance with the measured signal rates of the discovered Higgs boson, which at the LHC is dominantly produced through gluon fusion. The restriction to the order O(α t ) corrections with large top/stop masses in the loops furthermore ensures the approximation of zero external momenta to be reliable. This approximation breaks down if the masses of the particles running in the loops are small. 3 For the numerical analysis presented in section 3 we took care to choose scenarios where all possibly involved loop particles are sufficiently heavy so that not only the approximation of zero external momenta works well but also the order O(α) corrections do not play a significant role. The full EW and the order O(α t ) corrections differ by less than 4% for the chosen scenarios as we explicitly verified.
In the following we give the analytic formulae for the one-loop order O(α t ) corrections to the trilinear Higgs self-couplings in the interaction basis at vanishing external momenta. These formulae are compact enough to be easily implemented in computer codes. In order to extract only the O(α t ) and later on the O(α t α s ) corrections we neglect all D-term contributions to the Higgs potential and to the stop mixing matrix i.e. we work in the gaugeless limit, where the electric coupling e and the W and Z boson masses M W and M Z are taken to be zero but the vacuum expectation value v and the weak mixing angle θ W are kept finite. In this approximation, the stop mass matrix reads where m t denotes the top quark mass and the effective higgsino mixing parameter and the ratio of the two VEVs v u and v d , have been introduced. The soft SUSY breaking masses mQ 3 and mt R are real, whereas the trilinear coupling A t ≡ |A t | exp(iϕ At ) is in general complex. The matrix is diagonalized by a unitary matrix Ut, rotating the interaction statest L andt R to the mass eigenstatest 1 andt 2 , The order O(α t ) corrections to the trilinear Higgs self-couplings in the interaction basis are decomposed as The first term denotes the unrenormalized part arising from the one-loop diagrams with tops and stops running in the loops. The explicit expressions for ∆ (1) λ UR φ i φ j φ k are given in Appendix B. The contributions from the parameter counterterms are collected in the second part ∆ (1) λ CT φ i φ j φ k . Their explicit expressions in terms of the counterterms, defined in the following, are given in Appendix C.
For the order O(α t ) and and the order O(α t α s ) corrections, we need to renormalize the following set of parameters [59] where t φ , φ = h d , h u , h s , a d , a s denote the five independent tadpoles, M H ± stands for the mass of the charged Higgs boson and v ≈ 246 GeV is given by In order to renormalize the parameters, they are replaced by the renormalized ones and the corresponding counterterms as follows: Here the superscript (1) denotes the counterterms of O(α t ) and the superscript (2) the counterterms of O(α t α s ).
In addition to the parameter renormalization, also the wave function renormalization of the Higgs fields is needed in order to obtain a UV finite result. At O(α t ) and O(α t α s ), only the Higgs doublet H u has a non-vanishing wave function renormalization counterterm [59], which is introduced as The parameters are renormalized in a mixed OS-DR renormalization scheme as described in [59].
In this scheme part of the parameters, that are directly related to "physical" quantities, are renormalized on-shell, and the remaining parameters are defined via DR conditions, as While it is debatable if the tadpole parameters can be called physical quantities, their introduction is motivated by physical interpretation, so that in slight abuse of the language we call their renormalization conditions on-shell. For the wavefunction renormalization of the Higgs fields, the DR scheme is employed. Note that this procedure is applied for both the order O(α t ) and the order O(α t α s ) corrections. We do not repeat the renormalization conditions here, since they are introduced in detail in [59]. We give, however, the explicit expressions for the counterterms of order O(α t ). For the OS renormalization constants at order O(α t ) we find in D = 4 − 2 dimensions: And for the DR renormalization constants we get: Here φ x = ϕ u + ϕ s + ϕ λ + ϕ At , with ϕ λ and ϕ At being the complex phase of λ, and accordingly of A t , has been introduced. Furthermore we use c β ≡ cos β etc. The functions A 0 m 2 and B 0 p 2 , m 2 1 , m 2 2 denote the scalar one-point and two-point functions, respectively, in the convention of [85], and Q is the renormalization scale.

The order O(α t α s ) corrections
In order to obtain the order O(α t α s ) corrections we use the Feynman diagrammatic approach in the approximation of zero external momenta. These corrections are composed of The first part consists of the contributions from genuine two-loop diagrams. These must contain either a gluon or gluino or a four-stop coupling in order to give a contribution of order O(α t α s ).
Some sample diagrams are presented in Fig. 1. 5 In the approximation of zero external momenta all two-loop three-point functions can be reduced to the product of two one-loop tadpoles and to the two-loop one-point integral which are presented analytically in the literature [86][87][88][89][90][91][92].  The second part ∆ (2) λ CT1L φ i φ j φ k denotes the contributions arising from the one-loop diagrams with top quarks and stops as loop particles and with one insertion of a counterterm of order O(α s ) from the top/stop sector. Some representative diagrams for this set are depicted in Fig. 2. The parameters of the top/stop and bottom/sbottom sectors are renormalized at order O(α s ). The bottom quarks are treated as massless, so that the left-and right-handed sbottom states do not mix and only the left-handed sbottom with a mass of mQ 3 contributes. We choose the set of independent parameters entering the top/stop and bottom/sbottom sector, that we renormalize, to be given by Note that A t is in general complex. We renormalize these parameters both in the DR and in the OS scheme. The definition of their counterterms can be found in [59] 6 . According to the SUSY Les Houches Accord (SLHA) [93,94] convention, mQ 3 , mt R and A t are given as DR parameters. When we choose the OS scheme these parameters need finite shifts for the conversion into OS parameters. In the DR scheme on the other hand, the given top pole mass must be translated into a DR mass. These translations according to our conventions are described in detail in [59].
The third part consists of contributions arising from the order O(α t α s ) counterterms. The explicit expressions of ∆ (2) (2) tan β, δ (2) |λ| and δ (2) Z Hu are the same as in the one-loop case after replacing the one-loop by the two-loop counterterms. The formulae are given in Appendix C. For the exact definitions of the two-loop counterterms, we refer the reader to [59].
Our results have been obtained in two independent calculations. For the generation of the amplitudes we have employed FeynArts [95,96] using in one calculation a model file created by SARAH [68][69][70]97] and in the other calculation a model file based on the one presented in [98] which has been extended by our group to the case of the NMSSM. The contraction of the Dirac matrices was done with FeynCalc [99]. The reduction to master integrals was performed using the program TARCER [100], which is based on a reduction algorithm developed by Tarasov [101,102] and which is included in FeynCalc. We have applied dimensional reduction [103,104] in the manipulation of the Dirac algebra and in the tensor reduction. In our calculation no γ 5 terms appear that require a special treatment in D dimensions, so that we take γ 5 to be antisymmetric with all other Dirac matrices. The cancellation of the single pole and double poles has been checked. The results of the two computations are in full agreement. We furthermore compared our results in the limit of the real MSSM with Ref. [105] where the two-loop O(α t α s ) corrections to the MSSM Higgs self-couplings were given, and we found agreement between the two computations.

Numerical analysis 3.1 Scenarios
For the numerical analysis of the impact of the higher order corrections on the Higgs selfcouplings we made sure to choose scenarios that comply with the experimental constraints. In order to find viable scenarios we performed a scan in the NMSSM parameter space. We checked the scenarios for their accordance with the LHC Higgs data by using the programs HiggsBounds [106][107][108] and HiggsSignals [109]. The programs require as inputs the effective couplings of the Higgs bosons, normalized to the corresponding SM values, as well as the masses, the widths and the branching ratios of the Higgs bosons. These have been obtained for the SM and NMSSM Higgs bosons from the Fortran code NMSSMCALC [74,75]. A remark is in order for the loop-induced Higgs couplings to gluons and photons. The effective NMSSM Higgs boson coupling to the gluons normalized to the corresponding coupling of a SM Higgs boson with same mass is obtained by taking the ratio of the partial widths for the Higgs decays into gluons in the NMSSM and the SM, respectively. The QCD corrections up to next-to-next-to-next-to leading order in the limit of heavy quarks [110][111][112][113][114][115][116][117][118][119] and squarks [120,121] are included. As the EW corrections are unknown for the NMSSM Higgs boson decays, they are consistently neglected also in the SM decay width. The loop-mediated effective Higgs coupling to the photons has been obtained analogously. Here the NLO QCD corrections to quark and squark loops including the full mass dependence for the quarks [113,[122][123][124][125][126][127] and squarks [128] are taken into account. The EW corrections, which are unknown for the SUSY case, are neglected also in the SM.
For the numerical analysis of the corrections to the Higgs self-couplings, we have chosen two parameter sets that fulfill the above constraints. For both scenarios we use the SM input parameters [129,130] The remaining input parameters are given by And the remaining input parameters are set as follows, |λ| = 0.629 , |κ| = 0.208 , |A κ | = 179.7 GeV , |µ eff | = 173.7 GeV , We follow the SLHA format, which requires µ eff as input parameter. The values for v s and ϕ s can then be obtained by using Eq. (2.9). In the SLHA format, the parameters λ, κ, A κ , µ eff , tan β as well as the soft SUSY breaking masses and trilinear couplings are understood as DR parameters at the scale µ R = M s 7 , whereas the charged Higgs mass is an OS parameter. We set the SUSY scale M s to  The resulting supersymmetric particle spectrum from the thus chosen parameter values is in accordance with present LHC searches for SUSY particles [131][132][133][134][135][136][137][138][139][140][141][142][143][144][145]. Note, that in the following we will drop the subscript 'eff' for µ. Furthermore, whenever we will use the expressions OS and DR these refer to the renormalization in the top/stop sector.

Results for the loop-corrected self-couplings
In this and the following subsection we discuss the impact of the order O(α t α s ) corrections on the trilinear Higgs self-couplings and on Higgs-to-Higgs decay widths. We start by discussing the results for the parameter set called scenario 1 in the previous subsection. The masses of the Higgs bosons and their main composition in terms of singlet/doublet and scalar/pseudoscalar components at tree level, one-loop and two-loop order are summarized in Table 1   For definiteness, with respect to the mass corrections one-loop means here and in the following that we include the full EW corrections at non-vanishing external momenta, while at two-loop level the order O(α t α s ) corrections are computed at vanishing external momenta. As can be inferred from the table, the masses of the three lightest scalars are substantially different, so that mixing effects due to CP-violation for non-vanishing phases cannot be expected to be significant. The reason for choosing this scenario are higher order corrections to the trilinear Higgs self-coupling of the SM-like Higgs boson which are rather important for this parameter point. This boson is given by the state with the largest h u component and a mass value around 125 GeV. 8 At tree level it is the lightest Higgs boson H 1 that is mainly h u -like, and its mass thus receives large corrections which are dominantly stemming from the top/stop sector. The large corrections shift the H 1 mass above the one of H 2 so that the two Higgs bosons interchange their roles, as they are ordered by ascending mass. At one-and two-loop level it is therefore the second lightest Higgs boson, which is h u -like. For convenience, we denote in the following the mass eigenstate that is dominantly h u -like, by h. Furthermore, when we perform comparisons in the interaction basis at different loop-levels, these will be done for the h u state.
In Fig. 3 we show the dependence of the one-and two-loop corrections to the Higgs selfcoupling on the DR parameter A t in the two different renormalization schemes applied in the stop sector. The one-loop corrections have been obtained at order O(α t ) for vanishing external momenta. We explicitly verified, that the differences between the one-loop result in this approximation and the one including the full one-loop corrections for non-vanishing momenta at the threshold 9 are below 4% for the investigated parameter points. Two-loop corrections always refer to the order O(α t α s ) corrections at vanishing external momenta. The left plot of Fig. 3 shows the corrections to the self-coupling of h u in the interaction basis, λ huhuhu .  in both renormalization schemes. In the OS scheme the one-loop correction increases it by 140% while it is decreased by 24% to two-loop order. In the DR scheme the increase is of 74% going from tree-to one-loop order supplemented by another increase of 9% when adding the two-loop corrections. The reason, why the one-and two-loop corrections differ much more in the OS scheme than in the DR scheme can be understood as follows. In the DR scheme the top quark mass, which according to the SLHA accord is an OS parameter, has to be converted to the DR value. Thereby, the finite counterterm to the top mass, which in the OS scheme is included at two-loop level, is already induced at one-loop level in the value of the DR mass. In this way some corrections of order O(α t α s ), which in the OS scheme only appear at the two-loop level, are moved to the one-loop level, cf. also [59].
The lower panels of Fig. 3  where H both refers to the h u dominated mass eigenstate h, and to the h u interaction eigenstate. This value gives a rough estimate of the theoretical error in the Higgs self-coupling due to the unknown higher order corrections. In the interaction eigenstate it is of order O(50%) at one-loop level, decreasing to roughly 4% at two-loop level. In the mass eigenstate it is about 5% higher at both loop orders. The inclusion of the two-loop corrections hence substantially decreases the theoretical uncertainty. Figure 4 shows the same as Fig. 3 but now as a function of the phase ϕ At . All other CPviolating phases have been kept to zero. The figure shows that the dependence of the loop corrections on the phase is almost negligible, as expected for radiatively induced CP-violation. The size of the loop corrections and the remaining theoretical uncertainty are of the same order as for the variation of A t .
We now turn to the discussion of scenario 2. The masses and dominant composition of the mass eigenstates at tree level, one-and two-loop order are summarized in  one-loop level and from one-to two-loop level. In contrast to scenario 1 the masses of H 2 and H 3 are now much closer together, in particular after inclusion of the two-loop corrections. We therefore expect CP-violating effects to be more important here. The H 2 state is identified with the discovered Higgs boson. The stop masses are again rather heavy with OS : mt 1 = 1145 GeV , mt 2 = 1421 GeV , DR : mt 1 = 1126 GeV , mt 2 = 1387 GeV . (3.46) In Fig. 5 we show the dependence of the Higgs self-coupling of the h u state in the interaction basis (left) and of the h u -like mass eigenstate h (right) at one-(dashed) and two-loop order (full) for DR renormalization in the top/stop sector as a function of the phases ϕ M 3 , ϕ At and ϕ µ . For illustrative purposes we have varied the phases in rather large ranges although they might already be excluded by experiment. We start from our original CP-conserving scenario and turn on the phases one by one. Note, that ϕ µ has been varied such that the CP-violating phase ϕ y = ϕ κ − ϕ λ + 2ϕ s − ϕ u , that appears already at tree level in the Higgs sector, remains zero, i.e. ϕ λ and ϕ s were varied at the same time as ϕ λ = 2ϕ s = 2/3ϕ µ , while ϕ κ and ϕ u are kept zero. As expected, the loop-corrected couplings show a somewhat larger dependence on ϕ At than in scenario 1, in particular in the mass eigenstate basis. Defining as for the two-loop corrected self-coupling. Note, that the one-loop corrected self-couplings show a dependence on the phase of M 3 , although the genuine diagrammatic gluino corrections only appear at two-loop level. This dependence enters through the conversion of the OS top quark mass to the DR mass. Overall, the dependence of the loop corrected self-couplings on the CP-violating phases is smaller in the interaction states than in the mass eigenstates, which are obtained by rotating the interaction states with the mixing elements obtained from the loop corrected masses, that also depend on the CP-violating phases.
The lower panels show the relative corrections, defined at order n = 1, 2 as (3.49) In the interaction basis they are of order ∼ 70 − 80% for the one-loop corrections relative to the tree-level coupling and are somewhat larger than the corresponding values in the mass eigenstate basis, which are of order ∼ 50−60%. For the two-loop coupling relative to the one-loop coupling the corrections are significantly reduced to about 5 − 8% in both the interaction and the mass eigenstate basis. The two-loop corrections hence considerably reduce the theoretical uncertainty.
In order to further study the theoretical uncertainty, we show in Fig. 6 for scenario 1 (left) and scenario 2 (right) the scale variation of the trilinear Higgs self-coupling in the mass eigenstate h at one-and at two-loop order. We have varied the renormalization scale µ R between 1/3 and 3 times the central scale µ 0 = M s . The scale variation affects the DR parameters entering the calculation. In the absence of an implementation of the 2-loop renormalization group equations (RGE) for the complex NMSSM, which is devoted to future work, we obtain the parameters at the different scales by exploiting the relation between DR and OS parameters, as explained in Appendix D. This should approximate the results obtained from the RGE running rather well, in case the scale is not varied in a too large range. Since the scale variation provides only a rough estimate of the error made by neglecting higher order corrections this approach is sufficient for our purpose. As can be inferred from the figures, in scenario 1 the one-loop coupling is altered by up to 7% compared to its value at the central scale in the investigated range. This reduces to 2-5% at two-loop order. In scenario 2 the corresponding numbers at one-loop order are 3.5% compared to up to 2.5% for the two-loop coupling. As expected, the scale dependence reduces when going from one-to two-loop order. Note, however, that these numbers should not be taken as estimate for the residual theoretical uncertainty.

Phenomenological implications
We now turn to the discussion of the phenomenological implications due to the loop-corrected Higgs self-couplings. Higgs self-couplings are involved in Higgs-to-Higgs decays and in Higgs pair production processes. At the LHC, pair production dominantly proceeds through gluon fusion. This process, however, includes EW corrections beyond those approximated by the loopcorrected effective trilinear couplings. As they are not available at present we will not discuss Higgs pair production further and concentrate on Higgs-to-Higgs decays.
The decay width for the Higgs-to-Higgs decay H i → H j H k including the two-loop corrections to the Higgs self-coupling is obtained from where f = 2 for identical final state particles and f = 1 otherwise. The decay amplitude is denoted by M H i →H j H k and λ = (x − y − z) 2 − 4yz is the two-body phase space function. In case of CP-violation all Higgs-to-Higgs decays between the five neutral Higgs bosons are possible, if kinematically allowed, so that i, j, k = 1, ..., 5. In the CP-conserving case, however, only the trilinear Higgs couplings between three CP-even or one CP-even and two CP-odd Higgs bosons are non-vanishing. The matrix element is given by Here Γ h i h j h k is the loop corrected trilinear Higgs self-coupling in the tree-level mass eigenstate basis, where the Goldstone boson has been singled out. We here include at one-loop level the full electroweak corrections [79] at p 2 = 0, where we set the 4-momenta of the external Higgs particles equal to the respective loop-corrected Higgs mass values p 2 H i,j,k = M 2 H i,j,k as obtained with NMSSMCALC [51,58,59,75], and at two-loop the order O(α t α s ) corrections at p 2 = 0. 10 The proper on-shell conditions of the external Higgs bosons as required in the decay process are ensured by rotating the tree-level mass eigenstates h i ,j ,k to the loop corrected mass eigenstates H i,j,k with the matrix Z, cf. [51,79]. In this calculation we include at one-loop order the full electroweak corrections at non-vanishing external momenta. At two-loop order as usual the order O(α t α s ) corrections which are available only at p 2 = 0 are taken into account. 11 The δM mix H i →H j H k accounts for the contributions stemming from the mixing of the CP-odd components of the external Higgs bosons with the Goldstone and with the Z boson, respectively. These contributions, which are evaluated by setting the external momenta to the tree-level masses in order to maintain gauge invariance, are small already at one-loop order compared to the remaining contributions to the decay amplitude, as has been shown in [79]. We hence do not include the two-loop contributions, that can safely be expected to be negligible.
In the plots below we show apart from the two-loop corrected decay widths also the ones at one-loop order. The only change required to adapt formula (3.51) to this case is in Γ h i h j h k where solely the one-loop corrections to the vertex functions together with the corresponding counterterms are included. In particular we also use the two-loop corrected mass eigenstates and mixing matrix elements for the external particles (apart from the mixing contribution with the Goldstone and Z boson of course).
The scenario which the following discussion is based on and which has been checked to be compatible with the constraints from the LHC Higgs data, is given by: The remaining input parameters are given by This results in the Higgs mass spectrum given in Table 3 at tree, one-and two-loop level together with the main singlet/doublet and scalar/pseudoscalar components at each loop level. In Fig. 7   parameter point of scenario 3 with vanishing phase ϕ µ = 0. The phase is then varied in the range −π...π, such that at tree level the CP-violating phase ϕ y in the Higgs sector vanishes. As expected the dependence of the decay width on the CP-violating phase induced through the loop corrections is small, remaining below the per-cent level. For ϕ µ = 0 the tree-level decay width in the OS scheme is 0.171 GeV and 0.186 GeV in the DR scheme. 12 In the latter the one-loop corrections increase the decay width by 6.5% and the two-loop corrections add another 2.0% on top of that. In the OS we find a 21% increase at one-loop and a 7% decrease at two-loop so that at two-loop order the results in the two renormalization schemes approach each other, as can also be inferred from the lower left plot, which shows that the dependence on the renormalization scheme decreases from one-to two-loop level. Depending on the renormalization scheme different input parameters of the top/stop sector have to be converted to match the required renormalization scheme. The corresponding induced two-loop corrections into the one-loop corrections lead to a different dependence on the CP-violating phase of the two schemes at one-loop level. At twoloop level this difference in the phase dependence is then almost washed out and the scheme dependence is about 4.58% independent of ϕ µ . For the computation of the branching ratio of the decay, shown in Fig. 7 (right), we replace in the program package NMSSMCALC the tree-level decay widths Γ(H i → H j H k ) with our loop-corrected ones. The branching ratio, which with O(3%) is very small shows the same trend as the decay width with respect to the loop corrections.
The non-vanishing CP-violating phase induces through the higher order corrections CP mixing in the Higgs mass eigenstates, such that otherwise not allowed decays of e.g. the CP-odd doublet-like H 3 into a pair of SM-like H 2 bosons are possible. The branching ratio remains, however, tiny, reaching at most 0.58 per mille for |ϕ µ | ≈ π/2 in our scenario.

Conclusions
The search for New Physics and the proper interpretation of the experimental data requires from the theoretical side precise predictions of parameters and observables. In this work we computed the two-loop corrections to the trilinear Higgs self-couplings in the CP-violating NMSSM. Originating from the Higgs potential the Higgs boson masses and self-couplings are related to each other. For a consistent interpretation therefore the level of accuracy of the self-couplings has to match the one of the masses, that have been provided previously up to the two-loop level. Here, the two-loop corrections to the self-couplings have been calculated at the same order O(α t α s ) at vanishing external momenta. We have allowed for two renormalization schemes in the top/stop sector, namely OS and DR renormalization. Depending on the scenario and the renormalization scheme, the two-loop corrections are of the order of 5-10% relative to the one-loop couplings, compared to up to 80% for the one-loop corrections relative to the tree-level values. The investigation of the remaining theoretical uncertainty performed by varying the renormalization scheme of the top/stop sector or by changing the renormalization scale confirmed that the theoretical error is reduced through the inclusion of the two-loop corrections. As expected the dependence on the CP-violating phase due to radiatively induced CP-violation is small and of the order of a few percent.
The trilinear self-couplings are relevant in Higgs pair production processes and Higgs-to-Higgs decays, which now become accessible at run 2 of the LHC. While Higgs pair production requires the inclusion of further higher order corrections beyond the loop-corrected Higgs self-couplings provided in this work, the inclusion of the radiatively corrected trilinear Higgs self-coupling improves the prediction for the Higgs-to-Higgs decay rates and related branching ratios. For the investigated scenario and decay we find that the two-loop corrections alter the decay width by 2% (7%) in the DR (OS) scheme with respect to the one-loop level, which is to be compared to about 7% (21%) when going from tree-to one-loop level. The dependence on the renormalization scheme is reduced from ∼ 4.8 − 5.4% in the investigated range of the phase ϕ µ at one-loop level to ∼ 4.5% at two-loop level. The behaviour in the branching ratio is similar to the one of the decay width.
In summary, the inclusion of the two-loop corrections at O(α t α s ) in the approximation of vanishing external momenta in the trilinear Higgs self-couplings of the CP-violating NMSSM Higgs sector is necessary to match the available precision in the Higgs masses and to allow for a consistent interpretation of the Higgs data. Being of the order of 10% they have been shown to further reduce the theoretical error due to missing higher order corrections.
ijk (i, j, k = 1, ..., 6), with the same correspondences as introduced in Appendix A, can be cast into the form , (B.58) and the non-vanishing couplings read For the parameters p that are renormalized at one-loop order only, this relation simplifies to p pureDR (µ 1 ) − p pureDR (µ 2 ) = a 1 ln µ 2 1 µ 2 2 .