The Trilinear Higgs Self-Couplings at O ( α 2 t ) in the CP-Violating NMSSM

,


Introduction
The measurement of the trilinear Higgs self-coupling is one of the most important tasks at the LHC and future colliders [1]. It is the first step towards the experimental reconstruction of the Higgs potential and hence the direct experimental verification of the Higgs mechanism sui generis [2][3][4]. At the LHC, it is accessible in gluon fusion into Higgs pairs. In models beyond the Standard Model (SM) with extended Higgs sectors, the Higgs self-couplings are also involved in Higgs-to-Higgs decays. Through the Higgs potential, the trilinear Higgs self-coupling is related to the Higgs boson mass. While in the SM the Higgs mass is an ad hoc input parameter, in supersymmetric theories [5][6][7][8][9][10][11][12][13][14][15][16][17][18] it is derived from the parameters of the model. In the Minimal Supersymmetric extension of the SM (MSSM) [19][20][21][22] the Higgs quartic couplings are given in terms of the gauge couplings leading to an upper bound of the tree-level mass of the order of the Z boson mass so that considerable higher-order corrections are required to shift the SM-like Higgs boson mass to the observed value of 125.09 GeV [23]. In the Next-to-MSSM (NMSSM) [24][25][26][27][28][29][30][31][32][33][34][35] the situation is somewhat more relaxed due to the tree-level contribution stemming from the inclusion of the additional complex singlet field. In the last years, a lot of effort has been put to provide precise predictions for the Higgs mass values at higher loop level both in the MSSM and the NMSSM. For recent reviews, see [36,37]. For the trilinear Higgs self-couplings the corresponding corrections have not yet been provided at the same level of precision as for the masses. In the MSSM, the one-loop corrections to the effective trilinear couplings have been provided many years ago in [38][39][40]. The process-dependent corrections to heavy scalar MSSM Higgs decays into a lighter Higgs pair have been calculated in [41,42]. The two-loop O(α t α s ) SUSY-QCD corrections to the top/stop-loop induced corrections have been made available within the effective potential approach in [43]. In the NMSSM, we provided the full one-loop corrections for the CP-conserving NMSSM [44]. They are sizeable so that the inclusion of the two-loop corrections is mandatory to reduce the theoretical uncertainties due to missing higher-order corrections. Consequently, we subsequently calculated the two-loop O(α t α s ) corrections in the limit of vanishing external momenta in [45], in the CP-violating NMSSM. The full one-loop corrections to the Higgs-to-Higgs decays and other on-shell twobody decays were implemented in [46]. For corrections to the trilinear Higgs self-couplings in non-supersymmetric (non-SUSY) Higgs models, see for example Refs. [47][48][49][50][51][52] for one-loop and Refs. [53][54][55][56] for two-loop results, and Refs. [57][58][59][60][61][62][63][64][65][66] for the process-dependent Higgs-to-Higgs decays at one-loop level.
In this paper we continue our effort in increasing the precision for the predictions of the trilinear Higgs self-couplings in the context of the NMSSM. We calculate the two-loop corrections at O(α 2 t ) to the trilinear Higgs self-couplings of the complex NMSSM. They are obtained in the limit of zero external momenta and vanishing gauge couplings. We consistently apply the same renormalisation schemes as in our computation of the loop-corrections to the Higgs boson masses, based on a mixed on-shell-DR renormalisation in the Higgs sector and the possibility to choose between on-shell (OS) and DR renormalisation in the top/stop sector. Our corrections have been implemented in our Fortran code NMSSMCALC [67,68] and can be downloaded from the URL: https://www.itp.kit.edu/~maggie/NMSSMCALC/ The paper is organized as follows. In Sec. 2 we introduce the tree-level sectors of the NMSSM that are relevant for our calculation and set up our notation. In Sec. 3 we give the definitions for the loop-corrected effective trilinear Higgs self-couplings and for the loop corrections to the Higgs-to-Higgs decays. We specify the approximations that we apply and the renormalisation schemes that we use. In Sec. 4 we briefly present the set-up of our numerical analysis and the scan that we performed. Sections 5 to 7 are dedicated to our numerical analysis. In Sec. 5 we discuss the impact of our corrections on the effective trilinear Higgs self-couplings and on the Higgs-to-Higgs decays for two specific parameter points, in Sec. 6 we investigate these effects for our whole sample to get a more general picture. The effects of our corrections in the context of Higgs pair production are analysed in Sec. 7. Our conclusions are given in Sec. 8.

The Tree-Level NMSSM
In order to set our notation we briefly introduce the two sectors relevant for the renormalisation, the Higgs and the top/stop sectors. While the computation of the O(α 2 t ) contributions to the Higgs self-energies and tadpoles involves neutralinos and charginos, they do not need to be renormalised as vertices and propagators with these particles only enter at the two-loop level. For the definition of the electroweakino masses and mixing angles in the gaugeless limit we refer to Ref. [69]. We work in the Z 3 symmetric NMSSM including CP violation. For the computation of the two-loop corrections to the trilinear Higgs self-couplings of the neutral Higgs bosons at O(α t α s + α 2 t ) we apply the gaugeless limit and follow the notation of our previous calculations [45,[69][70][71]. Note that in contrast to the MSSM, the trilinear and quartic Higgs couplings in the NMSSM do not vanish in the gaugeless limit, but involve the parameters λ, κ, A λ , A κ . The NMSSM superpotential is given by in terms of the quark and lepton superfieldsQ,Û ,D,L,Ê, the Higgs doublet superfieldsĤ d , H u and the singlet superfieldŜ. Charge conjugated fields are denoted by the superscript c.
We have suppressed color and generation indices for better readability. The symplectic product x · y = ij x i y j (i, j = 1, 2) is built with the anti-symmetric tensor 12 = 12 = 1. Working in the CP-violating NMSSM, the parameters λ, κ are in general complex. All y x (x = e, d, u) are taken to be real by rephasing the left-and right-handed Weyl-spinor fields as x L,R → x L,R e iϕ L,R . In our computation, the Yukawa couplings y x are assumed to be diagonal in flavour space, and we only include y t while all other Yukawa couplings are set to zero. The soft SUSY breaking Lagrangian is given by where again quark and lepton generation indices are suppressed. TheQ,ũ R ,d R andL,ẽ R denote the complex scalar components of the corresponding quark and lepton superfields. The soft SUSY breaking gaugino mass parameters M i (i = 1, 2, 3) of the bino, wino and gluino fieldsB, W i (i = 1, 2, 3) andG as well as the soft SUSY breaking trilinear couplings A x (x = λ, κ, u, d, e) are complex in the CP-violating NMSSM in contrast to the soft SUSY breaking mass parameters of the scalar fields, m 2 X (X = S, H d , H u ,Q,ũ R ,d R ,L,ẽ R ), which are real.

The Higgs Boson Sector
The tree-level Higgs boson potential is obtained from L soft, NMSSM , the F -terms of W NMSSM and the D-terms originating from the gauge sector, In the gaugeless limit, the U (1) Y and SU (2) L gauge couplings g 1 → 0 and g 2 → 0 while tan θ W = g 2 /g 1 is kept constant, with θ W being the weak mixing angle. This is equivalent to the limit of vanishing electric charge and tree-level vector boson masses, e, M W , M Z → 0, while keeping tan θ W constant.
The Higgs boson fields are expanded around their vacuum expectation values (VEVs) v u , v d , and v s as where ϕ u,s denote the CP-violating phases. We can replace the three VEVs by tan β, the SM VEV v and the effective µ parameter µ eff as Note that the MSSM limit is smoothly retraced by taking the limit λ, κ → 0, v s → ∞ and at the same time keeping µ eff and κ/λ constant. From the Higgs potential of Eq. (3) we obtain the tree-level tadpoles, the Higgs mass matrices and the trilinear Higgs self-couplings. For the tadpole coefficients we have Only five of the tadpoles are independent, since t au = t a d /t β . The neutral Higgs mass matrix in the interaction basis is obtained as and the charged one as (r, s = 1, 2) The trilinear couplings which need to be renormalised at two-loop level for the calculation of the O(α 2 t ) corrections in this paper, are obtained as The explicit expressions for the tadpoles and the squared mass matrices M φφ and M h + h − are given in Ref. [69] and those for the trilinear Higgs self-couplings can be found in the Appendix of Ref. [45]. The neutral Higgs mass eigenstates are obtained by a two-fold rotation that first separates the Goldstone component through the rotation R G (β n ), i.e. it transforms from the basis (h d , h u , h s , a d , a u , a s ) to (h d , h u , h s , a, a s , G 0 ), and afterwards rotates into the mass basis (h 1 , h 2 , h 3 , h 4 , h 5 , G 0 ) with the rotation matrix R, where the neutral Goldstone boson mass is equal to the Z boson mass, m G 0 = M Z , in 't Hooft Feynman gauge. It vanishes in the gaugeless limit. The charged Higgs fields are rotated to the mass eigenstates with a single rotation R G − (β c ), In the 't Hooft Feynman gauge the charged Goldstone boson mass is equal to the charged W boson mass, m G ± = M W , and vanishes in the gaugeless limit. At tree-level the rotation angles β n and β c coincide with β, β c = β n = β. They are distinguished here as β n and β c are mixing angles and do not need to obtain a counterterm which is not the case for β that arises from the ratio of the VEVs. It has to be renormalised and receives a non-vanishing counterterm.
After the renormalisation they are set equal to the tree-level value of β again. Our tree-level masses are denoted by small letters m, apart from the charged Higgs boson mass. When we talk about loop-corrected masses, they are denoted by capital M . In our renormalisation of the trilinear coupling we will adapt the same renormalisation conditions as those used in the two-loop corrections of the masses.
We apply the SUSY Les Houches Accord (SLHA) [72,73] and in accordance with this accord decompose the complex parameters A λ and A κ into their imaginary and real parts. While in our program code NMSSMCALC also λ and κ are read in in terms of their real and complex part in accordance with the SLHA, internally, we choose a different, more convenient, parametrisation. We decompose λ and κ into their absolute values and phases ϕ λ and ϕ κ . We note that the phases enter the tree-level Higgs mass matrix in two combinations together with ϕ u and ϕ s , where ϕ y is the only CP-violating phase at tree level in the Higgs sector. If ϕ y = 0, the CP-even components, h u , h d , h s , hence do not mix with the CP-odd ones, a d , a u , a s . We use the tadpole conditions to replace ImA λ,κ as well as m 2 H u,d ,S by the tadpole parameters t a d ,as and t h d,u,s , respectively, cf. Ref. [69] for details.
In NMSSMCALC, we have two possibilities to choose the set of input parameters in the Higgs sector, either or In the first choice the charged Higgs mass is an input parameter while in the second one we have ReA λ as an input parameter.

The Top/Stop Sector
For the calculation of the Higgs self-couplings at the order O(α 2 t ), the top/stop sector needs to be renormalised at O(α t ). The top mass and the top-quark Yukawa coupling are related as, with m t and y t being real in our convention. Applying the freedom of choice of the phases ϕ L , ϕ R of the left-and right-handed top-quark fields, we define ϕ L = −ϕ R = −ϕ u /2. Thereby the stop mass matrix in the (t L ,t R ) T basis in the gaugeless limit is given by where Ut denotes the rotation matrix for the left-and right-handed stop fieldst L,R into the mass eigenstatest 1,2 . We set the bottom quark mass to zero everywhere so that the right-handed sbottom states decouple and only left-handed sbottom states appear in the computation. In the stop sector the parameters to be renormalised at one-loop level are 3 The Loop-Corrected Couplings

Definition
The renormalised trilinear Higgs self-couplingλ ijk at two-loop order between the interaction states h i , h j and h k is given bŷ Here the indices i, j, k refer to the interaction basis (h d , h u , h s , a d , a u , a s ). We denote the trilinear tree-level Higgs self-coupling by λ ijk and the one-and two-loop corrections to it by ∆ (1) λ ijk and ∆ (2) λ ijk , respectively. The explicit expressions for the tree-level couplings in the interaction basis are given in App. A of [45].
Applying the description in [45], we define the so-called effective trilinear Higgs self-couplings as follows. Both one-loop and two-loop corrections are computed in the approximation of zero external momenta, more specifically: • In the one-loop corrections, ∆ (1) λ (for simplicity, here and in the following we drop the indices 'ijk' where they are not needed), we include only corrections at O(α t ). These are coming from the top/stop sector and are hence the dominant ones. They have been discussed in detail in [45].
• For the two-loop corrections, where the QCD corrections ∆ αtαs λ have been computed in [45]. In this paper the ∆ α 2 t λ corrections are calculated for the first time.
• After calculatingλ ijk in the interaction basis (h d , h u , h s , a d , a u , a s ) we rotate it first to the basis (h d , h u , h s , a, a s , G 0 ) to single out the couplings with the neutral Goldstone bosons asλ • To obtain the effective trilinear couplings in the mass eigenstate basis we use the loopcorrected rotation matrix R l,eff . This matrix diagonalizes the loop-corrected mass matrix evaluated in the approximation of vanishing external momentum, where the indices a, b, c refer to the mass basis (H 1 , H 2 , H 3 , H 4 , H 5 ) and n, m, q to the interaction basis (h d , h u , h s , a, a s ). The matrix R l,eff is a 5 × 5 matrix and returned as an SLHA output of NMSSMCALC in the block NMHMIXC. Note that we denote the loop-corrected Higgs boson masses by capital letters (H i ) and the tree-level ones by lower letters (h i ).
We will later also calculate the Higgs-to-Higgs decays where we have to ensure the proper onshell conditions of the external Higgs bosons. In this case the one-loop corrections ∆ (1) λ include the full electroweak corrections together with non-vanishing momentum effects. They have been computed by us in the context of the CP-conserving and CP-violating NMSSM in Ref. [44] and Ref. [45], respectively. The two-loop part ∆ (2) λ contains the O(α t α s ) and O(α 2 t ) part, computed in the zero-momentum approximation as described above. Throughout, at two-loop order we apply the gaugeless limit. In order to ensure the proper on-shell conditions of the Higgs bosons, to the maximum extent possible in the context of our calculation, the amplitude for the decay process H a → H b + H c is computed by including the effect from the finite wave-function renormalisation factor matrix Z which is defined by where with R being the matrix that rotates the interaction eigenstates (h d , h u , h s , a, a s , G 0 ) to the tree-level mass eigenstates (h 1 , h 2 , h 3 , h 4 , h 5 , G 0 ). The definition of the matrix Z can be found in Ref. [44] for the CP-conserving case and Ref. [46] for the CP-violating case. 1

One-and Two-Loop Corrections
To be consistent, we compute the one-and two-loop corrections to the trilinear Higgs selfcouplings in accordance with the corresponding one-and two-loop corrections to the Higgs boson masses. This means we use the same renormalisation conditions in the higher-order corrections to the trilinear couplings as the ones we used in our computation for the masses. For our mass calculations, the detailed presentation of the one-loop corrections can be found in [74,75] and of the two-loop corrections up to order O(α t α s ) in [70], to order O(α 2 t ) in [69] and to order O((α t + α λ + α κ ) 2 ) in [71], together with the corresponding renormalisation conditions and the explicit definitions of the counterterms. The one-loop corrections to the trilinear Higgs self-couplings in the real NMSSM have been presented in [44] and to two-loop order O(α t α s ) in the CP-violating NMSSM in [45]. Throughout our computations we apply a mixed on-shell (OS)-DR renormalisation scheme. In the two-loop corrections which require the renormalisation of the top/stop sector we provide the option to choose between OS and DR renormalisation. All details can be found in the respective papers. Here we focus on a minimal description and refer the reader for further information to this literature.
In case the charged Higgs mass is used as independent input the parameters related to the Higgs sector that need to be renormalised are given by 2 and by for Re(A λ ) as independent input. Note, that if we apply the gaugeless limit we do not need to renormalise the neutral and charged gauge boson masses, M Z and M W , and the electric coupling e. For the Higgs fields, which need to be renormalised as well, we choose DR conditions. The details of the renormalisation procedure and the counterterms are given in the above mentioned papers so that we do not repeat them here.
The one-loop corrections ∆ (1) λ of the trilinear Higgs self-couplings can be decomposed as where the first term denotes the unrenormalised part given by the genuine one-loop diagrams.
For the O(α t ) correction, they comprise the one-loop diagrams with top and stops running in the loops and we restrict ourselves to the gaugeless limit. For the trilinear couplings used in the Higgs-to-Higgs decays we include the complete one-loop corrections at non-vanishing gauge couplings. The explicit expressions for the order O(α t ) corrections to the trilinear self-couplings are given in App. B and the counterterm expressions ∆ (1) λ CT are given in App. C of Ref. [45].
The two-loop corrections ∆ (2) λ of the trilinear Higgs self-couplings are composed of The unrenormalised part ∆ (2) λ UR consists of the genuine two-loop diagrams contributing at order O(α t α s ) and O(α 2 t ). Some sample diagrams for the newly computed O(α 2 t ) are depicted in Fig. 1. In the approximation of zero external momenta all two-loop three-point functions can be written in terms of products of one-loop integrals or the two-loop tadpole integral. Their analytic expressions are given in the literature [76][77][78][79][80][81][82]. The counterterm contributions ∆ (2)  to generate the model file including the vertex counterterms. The file was then used in FeynArts 3.10 [89,90] to generate all required one-and two-loop Feynman diagrams for the calculation of the corrections to the trilinear Higgs self-couplings. The evaluation of the fermion traces and the tensor reduction of the one-and two-loop integrals and the amplitudes with the counterterminserted diagrams was done with the help of FeynCalc 9.2.0 [91,92] and its TARCER plugin [93]. We performed three independent calculations which all agreed. We also explicitly checked the ultraviolet (UV)-finiteness of the loop-corrected Higgs self-couplings.
The calculation of the trilinear Higgs self-couplings at one-and two-loop order as well as the Higgs-to-Higgs decays including these corrections, has been implemented in NMSSMCALC [67] both for the CP-conserving and the CP-violating NMSSM. The new NMSSMCALC version 5.1 can be downloaded from the URL: https://www.itp.kit.edu/~maggie/NMSSMCALC/ The input file inp.dat includes the option to choose between the different loop orders in the trilinear couplings and correspondingly the Higgs-to-Higgs-decay widths. The effective trilinear Higgs self-couplings as defined above are given out in the output file.

The Parameter Scan
For the numerical discussion of our results we used the data set that we had generated for Ref. [71] by performing a scan in the NMSSM parameter space and keeping only those data sets that are in accordance with the relevant experimental constraints. We briefly summarise them here for convenience of the reader. We ensured compatibility with experimental constraints from the Higgs data by using HiggsBounds 5.9.0 [94][95][96] and HiggsSignals 2.6.1 [97]. The required effective NMSSM Higgs boson couplings normalised to the corresponding SM values were generated with NMSSMCALC. For valid points, χ 2 computed by HiggsSignals-2.6.1 needs to be consistent with an SM χ 2 within 2σ. 3 For this analysis, we checked the sample again with the updated HiggsBounds 5.10.2 and HiggsSignals 2.6.2 4 and found that more than 90% of the points (and in particular the two benchmark points discussed below) are still in the allowed region. We required one of the neutral CP-even Higgs bosons, called h from now on, to behave as the SM-like Higgs boson and to have a mass in the range Table 1: Scan ranges for the random scan over the NMSSM parameter space, withX =b R ,L,τ when including the two-loop corrections at O((α t +α λ +α κ ) 2 +α t α s ) in the default mixed DR-OS scheme specified above and with OS renormalisation in the top/stop and charged Higgs boson sectors as well as an infrared mass regulator M R with M 2 R = 10 −3 µ 2 R to treat the Goldstone problem. For details, we refer to [71]. The SM input values have been chosen as [98,99] In accordance with the SLHA format the soft SUSY breaking masses and trilinear couplings are understood as DR parameters at the scale This is also the renormalisation scale that we use in the computation of the higher-order corrections. The scan ranges of our input parameters are given in Tab. 1. Note, that both λ and κ are required to remain below 0.7 in order to roughly ensure perturbativity below the GUT scale. Also λ, κ, µ eff and tan β are understood to be DR parameters at the scale M SUSY according to the SLHA format. For the scan we kept all CP-violating phases equal to zero. We neglected parameter points with any of the following mass configurations, With the first condition (i) we avoid large logarithms in our fixed-order calculation. The second condition (ii) omits degenerate mass configurations for which the two-loop part of the NMSSMCALC code is not yet optimised. The third condition (iii) takes into account model-independent lower limits for the lightest chargino and stop masses.

Investigation of Specific Benchmark Points
In the following, we present results for two benchmark points. One point is the benchmark point P2OS from our investigation of the Higgs mass corrections at O(α 2 λκ ) ≡ O((α t +α λ +α κ ) 2 +α t α s ) in [71]. The other point is the benchmark point BP10 of Ref. [100]. They have been chosen such that the SM-like Higgs boson mass complies with our required mass window Eq. (34) at O(α 2 λκ ) when we choose OS renormalisation in the top/stop sector. The charged Higgs mass here and in all other results presented in the following is renormalised OS. The first parameter point P2OS features a large singlet admixture to the h u -like mass and is defined by the following input parameters: Parameter Point P2OS: All complex phases are set to zero and the remaining input parameters are given by We apply the SLHA format in which µ eff is taken as input parameter. From this we compute v s by using Eq. (7) (ϕ s is set to zero).
Since the trilinear Higgs self-couplings and the mass values are closely related through the Higgs potential a discussion of the higher-order corrections to the trilinear Higgs self-couplings should be completed by the information on the Higgs mass corrections. In Table 2 we hence give the mass values obtained for P2OS at tree level, at one-loop order and at two-loop level at O(α t α s ), O(α t (α s + α t )) and the latest computed two-loop order O(α 2 λκ ) for OS renormalisation in the top/stop sector, and in round brackets those for DR renormalisation in the top/stop sector. Note that the numbers slightly changed compared to those given in [71] due to a bug in the VEV counterterm. The changes are in the sub percentage level.
In the table we also list in square brackets the main singlet/doublet and scalar/pseudoscalar component of each mass eigenstate. At O(α 2 λκ ) the lightest Higgs boson h 1 obtains a mass of around 125.3 GeV. Since it is h u -like it couples maximally to top quarks so that the LHC Higgs signal strengths are reproduced and it hence behaves SM-like. In the following plots we will always label the Higgs bosons according to their dominant admixture 5 , as this determines the Higgs coupling strengths and consequently the size of the loop corrections. This allows us to consistently compare and interpret the impact of the loop corrections.
The second parameter point BP10 features a resonantly enhanced Higgs pair production cross section in gluon fusion and is given by:     The impact of the loop corrections on the Higgs boson masses has been discussed extensively in [71]. Let us therefore here state only the main features. The h u -like tree-level Higgs mass value changes considerably when one-loop corrections are included, with a smaller change in the DR scheme, as in this scheme we already partly resum higher-order corrections. The relative O(α t α s ) corrections compared to the one-loop result are at the several per-cent level and move the obtained mass values in the two renormalisation schemes closer to each other, whereas the additional inclusion of the O(α 2 t ) corrections increases the difference again (in the OS scheme). The newest corrections at O(α 2 λκ ) move the two values a little bit closer again.

Impact on the Effective Trilinear Higgs Self-Coupling
In Fig. 2 (left) we present for the parameter point P2OS with the large singlet admixture the effective trilinear Higgs self-couplingλ eff 111 , as defined in Eq. (27), of the dominantly h u -like Higgs boson for OS (full) and DR (dashed) renormalisation in the top/stop sector as a function of the stop trilinear coupling A t . The dominantly h u -like Higgs boson here always is the lightest mass eigenstate. Note, that here and in the following the A t is always the DR value. 6 Shown are the results at one-loop order (black), two-loop O(α t α s ) (blue) and at the newly calculated two-loop O(α t (α s + α t )) (red). Note that the loop-corrected rotation matrix R l,eff for the rotation to the mass eigenstates is taken consistently at the respective loop order. We plot here only the variation of A t between −500 and +500 GeV. In this region the phenomenology at O(α t (α s +α t )) is in accordance with the LHC Higgs data. 7 The steep decrease of the full red curve towards negative A t values is due to the approach to the cross-over point where the singlet-like and doublet h u -like Higgs state interchange their roles with respect to their mass ordering. This point is located outside of the shown region in the plot, at A t = −900 GeV, cf. Fig. 3 in Ref. [71]. The interchange of the roles of the two lightest mass eigenstates takes place due to the large singlet admixture for this parameter point which induces the transition between the h s -and h u -like interaction state. For the trilinear couplingλ ijk in the interaction basis after singling out the Goldstone boson, cf. Eq. (26), this of course does not occur. The multiplication with the mixing matrices R l,eff causes the mixing of the singlet and doublet states and also mixes in higher orders, as we do not evaluate the mixing matrix multiplication strictly at the considered loop order.
In order to quantify the impact of the new additional corrections we define -for a given renormalisation scheme -the relative change in the trilinear coupling value when going from loop order α i to the loop order α i+1 , which includes the next level of corrections, as The relative corrections amount to ∆ αtαs one-loop = 14-19% in the OS scheme when we include the O(α t α s ) corrections beyond one-loop order. When we include the O(α t (α s + α t )) corrections in addition to the available two-loop O(α t α s ) corrections the relative change is smaller with ∆ αt(αs+αt) αtαs = 1.5-18%. As expected the relative change decreases with increasing higher-order in the corrections. Note, that the relative corrections can be positive or negative. In the DR scheme the corrections are much smaller. We have for the relative O(α t α s ) corrections compared to one-loop order 0.8-4%, and the new corrections change the coupling by about 1%. The reason is that the DR scheme already partly resums higher-order corrections.
The lower panel in Fig. 2 shows the relative change in the corrections to the trilinear Higgs self-coupling at fixed loop order when we change the renormalisation scheme in the top/stop sector, 8  The comparison of the results in the two different renormalisation schemes can be used to estimate the uncertainty on the trilinear Higgs self-coupling due to missing higher-order corrections. As expected, the renormalisation scheme dependence is reduced when more higher-order corrections are included. The renormalisation scheme dependence continuously decreases from one-loop order with 26-33%, to 9-13% at O(α t α s ) and to 0.01-8% at O(α t (α s + α t )).
In Fig. 2 (right) we show our results for the benchmark point BP10. For this point, the relative corrections in the OS scheme are slightly larger than for P2OS. We have ∆ αtαs one-loop = 17-24% and ∆ αt(αs+αt) αtαs = 9-17%. In the DR scheme the corrections are smaller, the relative O(α t α s ) corrections compared to one-loop order are of 4-7%, and the new corrections change the coupling by 0.4-2%. The renormalisation scheme dependence decreases from one-loop order with 26-39% to 0.7-1.4% at O(α t α s ). The scheme dependence slightly increases again at O(α t (α s + α t )) where it is 7-13%. This is a behaviour that we already observed in the loop corrections to the Higgs boson masses [69]. 9 In summary, for both benchmark points the inclusion of the new two-loop corrections has an impact of a few per cent and we find a renormalisation scheme dependence of typical two-loop order. The behaviour is similar to the one we found for the Higgs mass corrections. Fig. 3 we show for the parameter point P2OS the loop corrections to the effective trilinear Higgs self-couplingλ eff 111 of the dominantly h u -like Higgs boson as a function of the CP-violating phase ϕ At of A DR t . The colour and line codes are the same as in Fig. 2. In order to avoid too large singlet-doublet mixing effects we chose |A DR t | = 250 GeV. CP violation due to the phase ϕ At is a loop-induced effect. Since A t enters at one-loop level, the trilinear Higgs selfcoupling shows a dependence on the CP-violating phase, which for this parameter point turns out to be larger in the OS than in the DR renormalisation scheme. The almost flat dependence on the CP-violating phase of the OS curve at O(α t α s ) is due to accidental cancellations which we explicitly checked. At O(α t (α s +α t )), we see a stronger dependence again. In the DR scheme both two-loop orders show about the same dependence on the phase ϕ At .

Impact on the Higgs-to-Higgs Decays
We now turn to the impact of the computed higher-order corrections on the partial decay widths for Higgs-to-Higgs decays. The decay width for the Higgs decay h i into a Higgs pair h j h k is given by where β(x, y, z) = (x − y − z) 2 − 4yz is the two-body phase space function and the decay amplitude M h i →h j h k is calculated according to Eq. (28). We show in Fig. 4 (left) the partial decay width of the doublet-like CP-even Higgs boson h d into a pair of a SM-like Higgs boson (right) at one-loop order (black), two-loop O(α t α s ) (blue) and two-loop O(α t (α s + α t )) (red) in the OS (full) and DR scheme (dashed) as a function of A DR t . In the right plot the dashed blue and red lines nearly lie on top of each other. Middle: The relative correction defined analogously to Eq. (39), but for the partial decay width. Lower: The relative renormalisation scheme dependence defined analogously to Eq. (40), but for the partial decay width, for all three loop orders.
h u and a singlet-dominated Higgs h s , Γ(h d → h u h s ), at one-loop level and at two-loop O(α t α s ) and O(α t (α s + α t )) for P2OS, as a function of A t . 10 This decay is the largest of the Higgs-to-Higgs decays for this parameter point. For BP10 the largest one is given by h s → h u h u which we show in Fig. 4 (right). This is also the resonant contribution that increases the production process of an h u h u Higgs pair which we will discuss later. We include both the Z matrix of Eq. (29) and the Higgs mass values calculated at the corresponding same loop order as the one for which we calculate the higher-order corrections to the trilinear Higgs self-coupling. This of course also implies that the kinematical factor in the decay amplitude changes with the loop order. For both parameter points, we observe a reduction of both the relative correction and the renormalisation scheme dependence when we move from one-to two-loop order with the effect being less pronounced in the DR than in the OS scheme as the former already partly resums higher-order corrections.
For P2OS the relative corrections for the partial decay width in the OS scheme amount to more than 100% when including the O(α t α s ) corrections in addition to the one-loop corrections. The reason is the small one-loop decay width. In the DR scheme the relative corrections amount to 10-14%. The relative effect of the new two-loop corrections is much less as expected and reaches ∆ αt(αs+αt) αtαs = 7-15% (0.4%) in the OS (DR) scheme. The renormalisation scheme dependence decreases from O(71-90%) at one-loop level to a maximum of 13% at two-loop O(α t α s ) and at most 14% at two-loop O(α t (α s + α t )). For the benchmark point BP10 we find that for most A t values the relative corrections of our new two-loop corrections are smaller compared to the relative corrections when moving from one-to two-loop order. Note that we cut some of the lines in the middle plot of Fig. 4 (right) as here the relative corrections become artificially large due to comparatively very small widths at the previous loop order. The reduction in the renormalisation scheme dependence when moving from one-to two-loop order is less obvious as can be seen from the lower panel in Fig. 4 (right). The renormalisation scheme dependence becomes artificially large here where the DR result for the partial decay width is very small.
In both scenarios the partial decay widths can be become as large as about 1.7 GeV. In P2OS this leads to a branching ratio of about 12% at most, taking into account the dominant decay channels. In BP10 we get a maximum branching ratio of more than 70%.
Our results show that the higher-order corrections to the decay width have a substantial impact in particular when moving from one-loop to two-loop order. Furthermore, also at the two-loop level the inclusion of the O(α 2 t ) corrections on top of the available O(α t α s ) corrections leads to significant changes in the decay width both for P2OS and BP10. At the phenomenological level, the impact of the changes depends on the relative size of the Higgs-to-Higgs decay widths compared to the other decay widths. 11

Scatter Plots
After the investigation of two specific benchmark points we aim to get an overall picture of the corrections by investigating scatter plots. These plots contain all parameter scenarios that we obtained from our scan and that comply with the included constraints described above. Figure 5 (left) displays for all generated valid parameter scenarios the effective trilinear Higgs self-couplingλ eff huhuhu of the Higgs mass eigenstate that is dominantly h u -like, at O(α t (α s + α t )) in the OS renormalisation scheme of the top/stop sector as a function of A t ≡ A DR t . Note that, depending on the parameter point, this is not necessarily always the same Higgs mass eigenstate. The right plot shows the renormalisation scheme dependence at the one-and considered two-loop orders. We see the same trend as observed for the benchmark points. The scheme dependence at one-loop order is rather large, varying between about 20% and more than 50%. It is considerably reduced upon inclusion of the two-loop O(α t α s ) corrections where it ranges between about 1 and 5%. After the additional inclusion of the O(α 2 t ) corrections the scheme dependence increases again to values between 5 and 18% and reflects the necessity to include all corrections at a given loop order in order to make a reliable statement on the scheme dependence. The scheme dependence at two-loop order is well below the one at one-loop level as expected.

The Trilinear Higgs Self-Coupling
Turning to the values of the trilinear Higgs self-couplingλ eff huhuhu we find that it lies between 190 and 228 GeV. For the SM-coupling we have for M H = 125.09 GeV and v = √ 2G F ≈ 246.22 GeV. Taking into account the residual theoretical uncertainty at O(α t (α s + α t )) due to missing higher-order corrections, the NMSSM h u -like trilinear Higgs self-coupling complies with the one found in the SM. This means that taking into account the LHC Higgs data results which push the discovered Higgs boson very close to the SM expectation, we find that the NMSSM h u -like trilinear Higgs self-coupling is also very SM-like once all dominant higher-order corrections are taken into account. This has important implications for the cross section values of Higgs pair production (see our discussion in Sec. 7). Note also that these values forλ eff huhuhu lie well within the present experimental limits on the SM trilinear Higgs self-coupling which are between −0.4 and 6.3 times the SM value as reported by ATLAS [101] and between −1.7 and 8.7 times the SM value as found by CMS [102] (both assuming a SM-like top-Yukawa coupling). in the OS and of roughly 3-12% in the DR scheme for the trilinear Higgs self-coupling. As for the masses, we find here 8-19% in the OS and 3-8% in the DR scheme. Both corrections are correlated, larger corrections in the trilinear Higgs self-coupling correspond to larger corrections for the Higgs mass. As can be inferred from Fig. 6 (lower), the relative size of the additional O(α 2 t ) corrections amounts to about 5-27% in the OS scheme and roughly 1-6% in the DR scheme for the trilinear Higgs self-coupling which compares to 3-11% in the OS scheme and 0.1-1.1% in the DR scheme for the masses. The mass corrections are in general smaller than the corrections to the trilinear Higgs self-couplings, and the corrections in the OS scheme are generally larger than in the DR scheme which partly resums higher-order corrections.

Correlation between Trilinear Higgs Self-Coupling and Mass
Overall we find for our parameter points compatible with the applied constraints that the effective trilinear Higgs self-coupling values at O(α t (α s + α t )) are in general smaller in the DR

Higgs Pair Production
In this section we want to analyse what we can learn from our higher-order results to the trilinear Higgs self-coupling about the impact of the electroweak corrections on Higgs pair production. Higgs pair production gives access to the trilinear Higgs self-coupling and the measurement of the Higgs self-interactions provides the ultimate test [1][2][3][4] of the Higgs mechanism for the generation of particle masses. At the LHC the dominant Higgs pair production process is given by gluon fusion into Higgs pairs [1,103,104]. The loop-induced process is mediated by top-quark loops and by bottom-quark loops, the latter contributing at the percent level. Higherorder QCD corrections are important, increasing the cross section by roughly a factor two at next-to-leading order (NLO). A lot of effort is put in providing increasingly precise predictions. The first NLO results were presented in the large top-quark mass limit more than two-decades ago [105]. Full NLO QCD corrections including the top-quark mass dependence were finally made available in [106][107][108][109]. The next-to-next-to-leading order (NNLO) QCD corrections in the large m t limit were provided by [110], and the next-to-next-to-leading logarithmic corrections in this limit by [111,112]. Recently, the corrections due to the resummation of soft-gluon emission were provided up to next-to-next-to-next-to-leading logarithmic accuracy in [113]. The NNLO FT approx 12 result was presented in [114]. A combination of the usual renormalisation and factorization scale uncertainties with the uncertainties originating from the scheme and scale choice of the virtual top mass was given in [115]. In the NMSSM we have additional diagrams involving top and bottom squarks as well as the s-channel exchange of non-SM-like Higgs bosons, cf. Fig. 7. In [44,100] we computed the NLO QCD corrections in the heavy-top limit.
So far the complete electroweak (EW) corrections for gluon fusion into Higgs pairs are not yet available. While in the SM we can expect them to be of the order of a few percent by looking at the EW corrections to single Higgs production [116][117][118][119][120] this might not be the case in beyond-the-SM models where couplings can be enhanced compared to the SM or where light Higgs bosons could run in the loops. The computation of the EW corrections to Higgs pair production through gluon fusion is a major task and technical challenge, which requires the computation of massive two-loop integrals with several different mass scales. First steps have been taken recently within the SM. In [121] the top-Yukawa induced part of the EW corrections and their relation to the effective trilinear Higgs coupling have been provided and discussed. The subset of two-loop diagrams where the Higgs boson is exchanged between the virtual top quark lines has been calculated in the high-energy limit in [122].
In this work we use our effective loop-corrected Higgs self-couplings that make up part of the EW corrections to get some insights on their importance. For this we choose the parameter point BP10 where the di-Higgs cross section is dominated by the resonant production of two SM-like Higgs bosons via an intermediate heavy Higgs boson, as in this case the other diagrams will give a subleading contribution to the total cross section so that the missing EW higher-order corrections might be of less importance. In the triangle diagram involving the resonant heavy Higgs boson in the s-channel we still miss, however, the EW corrections to the top triangle.
For this benchmark point we compute the gluon fusion production cross section for a pair of SM-like Higgs bosons. We choose a c.m. energy of 14 TeV, we use the CT14 pdf set [123] and we set the top quark mass to m t = 172.74 GeV. By using the loop-corrected Higgs masses as inputs and the corresponding higher-order corrected Higgs mixing angles to compute the Yukawa couplings and the trilinear Higgs self-couplings that enter the process through tree-levellike formulae, we take into account the higher-order corrections to the input parameters. By additionally including the loop-corrected trilinear Higgs self-coupling computed in this paper we explicitly include higher-order corrections to the observable, namely the Higgs pair production process, though at an incomplete level as mentioned above.
In Tab. 4 we compare the di-Higgs cross sections for the case where only corrections to the input parameters are considered (called 'inp' in the following) and the case where we additionally include EW corrections to the process through the loop-corrected trilinear Higgs self-couplings (called 'proc'). For simplicity we focus on the case where we include the full 1-loop corrections both to the masses/mixing angles [75] and to the trilinear Higgs self-couplings [44,45] (named '1L1L') and on the case where we include the 2-loop O(α t (α s + α t )) corrections both to the masses/mixing angles [69] and the trilinear Higgs self-couplings, computed in this paper (named 'at2at2'). We furthermore list in the This value is also given in the table as the resonant enhancement of the cross section basically comes from the resonant H 2 production with subsequent decay into H 1 H 1 . The resonant production from the H 3 decay into H 1 H 1 plays only a minor role for this benchmark point. We provide all these values both for OS and DR renormalisation in the top/stop sector. Note that in all renormalisation schemes and at all considered loop levels the Yukawa coupling of the SM-like Higgs H 1 is practically SM-like (it differs by only 1% from the SM-value).
From the cross section values we first of all notice the resonant enhancement compared to the SM Higgs pair production cross section, which at tree-level amounts to 19.72 fb. When we compare the absolute values of the cross section we have to be careful as the changes not only come from the use of different trilinear Higgs self-couplings and renormalisation schemes, but also from the change in the kinematics as the Higgs mass values depend on the loop order and the renormalisation scheme. In the last column of Tab. 4 we give the relative change of the cross section with respect to the applied renormalisation scheme, ∆ ren σ ≡ |σ mt(DR) − σ mt(OS) | σ mt(DR) .
We observe that the inclusion of only the parameter corrections does not decrease the renormalisation scheme dependence when moving from one-loop order to two-loop O(α t (α s + α t )), on the contrary. This is not astonishing as the scheme dependence of the input parameters has to be compensated by the scheme dependence of the process-dependent corrections at the same loop '1L1L'  Table 4: BP10: Cross section values for the production of a SM-like Higgs pair H 1 H 1 for OS and DR renormalisation in the top/stop sector when using loop-corrected masses and mixing angles ('inp') and additionally loop-corrected effective trilinear Higgs self-couplings ('proc') at 1-loop order and at O(α t (α s + α t )), respectively. The corresponding normalized trilinear Higgs self-couplings are given as well as the relative change of the cross section with the applied renormalisation scheme (see the text for definition).
order. When these are included we observe a decrease in the renormalisation scheme dependence of the cross section at the same loop order when including higher and higher loop orders as expected in perturbation theory. Still the renormalisation scheme dependence with 14.6% at O(α t (α s + α t )) is significant. This gives a hint that the remaining electroweak corrections that we did not take into account in our approach might be significant. It will hence be important to provide the complete EW corrections to the cross section to be able to reduce the uncertainty in its prediction due to missing higher loop corrections.

Conclusions
In this paper we presented the O(α 2 t ) corrections to the trilinear Higgs self-couplings in the context of the CP-violating NMSSM. They are part of our ongoing program of increasing the precision in the predictions of the NMSSM Higgs potential parameters, the masses and the Higgs self-couplings. We find that the corrections to the effective trilinear Higgs self-couplings are in general larger than those to the Higgs boson masses. The relative corrections on top of the already existing two-loop corrections at O(α t α s ) are much smaller, however, than when moving from one-to two-loop order and indicate perturbative convergence. The remaining residual theoretical uncertainties due to missing higher-order corrections, estimated from the variation of the renormalisation scheme in the top/stop sector, range at the per-cent level and are also reduced compared to the one-loop results. In general, the results obtained in the DR scheme show a better convergence than in the OS scheme, which is to be expected as they already partly resum higher-order corrections. Within the theoretical uncertainties, the obtained loopcorrected trilinear Higgs self-coupling of the SM-like Higgs boson is in accordance with the result for the SM Higgs boson with same mass value. The impact of the higher-order corrections on the Higgs-to-Higgs decay widths is similar to the one on the effective Higgs self-couplings. We also investigated the effect of the inclusion of our loop-corrected effective Higgs self-couplings in the Higgs pair production process. The estimates of the theoretical uncertainty based on the variation of the renormalisation scheme indicate that the remaining missing electroweak corrections to the process may be significant.