One-Loop Corrections to the Two-Body Decays of the Neutral Higgs Bosons in the Complex NMSSM

Since no direct signs of new physics have been observed so far indirect searches in the Higgs sector have become increasingly important. With the discovered Higgs boson behaving very Standard Model (SM)-like, however, indirect new physics manifestations are in general expected to be small. On the theory side, this makes precision predictions for the Higgs parameters and observables indispensable. In this paper, we provide in the framework of the CP-violating Next-to-Minimal Supersymmetric extension of the SM (NMSSM) the complete next-to-leading order (SUSY-)electroweak corrections to the neutralHiggs boson decays that are on-shell and non-loop induced. Together with the also provided SUSY-QCD corrections to colored final states, they are implemented in the Fortran code NMSSMCALC which already includes the state-of-the art QCD corrections. The new code is called NMSSMCALCEW. This way we provide the NMSSM Higgs boson decays and branching ratios at presently highest possible precision and thereby contribute to the endeavor of searching for New Physics at present and future colliders.


Introduction
The discovery of a scalar particle by the LHC experiments ATLAS [1] and CMS [2] and the subsequent investigation of its properties revealed a Higgs boson that behaves very Standard Model (SM)-like. Also years after its discovery there are no evidences for new physics from direct searches. In this situation the precise investigation of the Higgs sector plays an important role. Indirect effects of physics beyond the SM (BSM) might show up in the properties of the discovered Higgs boson. With a mass of 125.09 GeV [3] it does not exclude the possibility for the Higgs boson of a supersymmetric (SUSY) extension of the SM, like the minimal (MSSM) or the next-to-minimal (NMSSM) ones. Supersymmetry certainly belongs to the best motivated and most intensively studied BSM extensions, and the NMSSM, with a Higgs sector consisting of seven Higgs bosons arising after electroweak symmetry breaking (EWSB) from the two doublet and singlet fields of the Higgs sector, provides a rich phenomenology [4,5]. The experimental limits strongly restrict possible new physics effects in the Higgs sector and call for precision in the theory predictions for the Higgs boson observables. This is also necessary in order to be able to distinguish between new physics extensions in case of discovery.
In this paper, we concentrate on the NMSSM Higgs boson decays. While the (SUSY-)QCD corrections can be taken over from the MSSM case with the appropriate modifications and a minimum of effort, this is not the case for the electroweak (EW) corrections. In the recent years there has been some progress on this subject. In the CP-conserving NMSSM, members of our group computed the next-to-leading order (NLO) SUSY-EW and SUSY-QCD corrections to the decays of CP-odd NMSSM Higgs bosons into stop pairs and found that the both the EW and the SUSY-QCD corrections are significant and can be of opposite sign [6]. The authors of [7,8] provided in the framework of the CP-conserving NMSSM its full one-loop renormalization and the two-body Higgs decays at one-loop order in the on-shell (OS) renormalization scheme. A generic calculation of the two-body partial decays widths at full one-loop level was provided in [9] in the DR scheme. In [10] the full one-loop corrections for the neutral CP-violating NMSSM Higgs bosons were calculated to their decays into fermions and gauge bosons and combined with the leading QCD corrections. For the Higgs-to-Higgs decays, we provided in previous papers the complete one-loop [11] and the order O(α t α s ) two-loop [12] corrections in the CP-conserving and CP-violating NMSSM, respectively.
In this work, we compute, in the framework of the CP-violating NMSSM, the complete next-to-leading order (SUSY-)electroweak corrections to the neutral NMSSM Higgs boson decays into all tree-level induced SM final states, i.e. into fermion and massive gauge boson pairs, but also into non-SM pairs, namely gauge and Higgs boson final states, chargino and neutralino pairs, and into squarks. Where applicable we combine our corrections with the already available (SUSY-)QCD corrections. We furthermore include the complete one-loop corrections to the decays into Higgs boson pairs, cf. Refs. [11,12]. For the loop-induced decays into gluon and photon pairs as well as a photon and a Z boson no corrections are provided as they would be of two-loop order. For the first time, we present the one-loop corrections to the electroweakino, stop and sbottom masses in the context of the CP-violating NMSSM, by applying both OS and DR schemes. We have implemented our corrections in our original code NMSSMCALC [13], which calculates, based on a mixed OS-DR scheme, the NMSSM Higgs mass corrections and decays in both the CP-conserving and CP-violating case. This way we provide the NMSSM Higgs boson decays and branching ratios at presently highest possible precision including the state-of-the-art (SUSY-)QCD and the computed (SUSY-)EW corrections. In the EW higher-order corrections we not only include the NLO vertex corrections but also take into account the proper on-shell conditions of the external decaying Higgs bosons up to two-loop order O(α t α s + α 2 t ). This is the order up to which the mass corrections for the NMSSM Higgs bosons both in the CP-conserving [14] and CP-violating case [15][16][17] have been implemented in NMSSMCALC. The new program is called NMSSMCALCEW can be obtained at the url: https://www.itp.kit.edu/∼maggie/NMSSMCALCEW/. Here also a detailed description of the program and its structure are given, instructions on how to compile and run it as well as information on modifications, which is constantly updated. A brief description of the code is given in Appendix B.
The paper is organized as follows. In Section 2 we introduce the NMSSM sectors at tree level, that are relevant for our computation, and set our notation before we move on to the NMSSM at one-loop level in Section 3. We here describe the renormalization of the Higgs, chargino/neutralino, and squark sectors as well as the loop corrections to the Higgs boson masses and mixings, to the neutralino and chargino masses, and finally to the squark masses and their mixings. Section 4 is devoted to the detailed presentation of our calculation of the oneloop corrections to the neutral non-loop induced Higgs boson decays into on-shell final states, namely into fermion pairs, massive gauge boson pairs, final states with one gauge and one Higgs boson, neutralino and chargino pairs, and squark final states. In Section 5 we present the numerical analysis of the one-loop corrections to the Higgs boson branching ratios into SM and SUSY final states, where we discuss in particular the size of the newly implemented corrections to both the branching ratios and to the electroweakino and third generation squark masses. Our conclusions are given in Section 6. Explicit expressions of the counterterm couplings for the decays of neutral Higgs bosons into a squark pair are displayed in Appendix A.

The NMSSM at Tree Level
We are working in the complex NMSSM with a preserved Z 3 symmetry. The Lagrangian of the NMSSM can be divided into the Lagrangian of the MSSM and the additional part coming from the NMSSM. For convenience of the reader and to set our notation, we give here the parts of the Lagrangian that are relevant for our calculations. For the Higgs sector we need the NMSSM Higgs potential. It is derived from the NMSSM superpotential W NMSSM , the corresponding soft SUSY-breaking terms, and the D-term contributions. With the Higgs doublet superfieldŝ H u andĤ d coupling to the up-and down-type quark superfields, respectively, and the singlet superfieldŜ, we have for the NMSSM superpotential in terms of the left-handed quark and lepton superfield doubletsQ andL and the right-handed up-type, down-type, and electron-type superfield singletsÛ ,D, andÊ, respectively. Charge conjugation is denoted by the superscript c, and color and generation indices have been omitted. The NMSSM superpotential contains the coupling κ of the self-interaction of the new singlet superfield and the coupling λ for theŜ interaction with the two Higgs doublet superfields. Both couplings are complex. The quark and lepton Yukawa couplings y d , y u , and y e are in general complex. However, in case of no generation mixing, as assumed in this paper, the phases of the Yukawa couplings can be absorbed through a redefinition of the quark fields, so that the phases can be chosen arbitrarily without changing the physical meaning [18]. The soft SUSY-breaking NMSSM Lagrangian in terms of the component fields H u , H d and S reads 3) It contains two more complex parameters specific to the NMSSM, the soft SUSY-breaking trilinear couplings A λ and A κ . The soft SUSY-breaking MSSM contribution can be cast into the form where the SM-type and SUSY fields corresponding to a superfield (denoted with a hat) are represented by a letter without and with a tilde, respectively. The indicesQ (L) of the soft SUSY-breaking masses denote, exemplary for the first generation, the left-handed quark (lepton) doublet component fields of the corresponding quark and lepton superfields, andũ R ,d R ,ẽ R the right-handed component fields for the up-type and down-type quarks, and charged leptons, respectively,. The trilinear couplings A u , A d and A e of the up-type and down-type quarks and charged leptons are in general complex, whereas the soft SUSY-breaking mass terms m 2 x (x = S, H u , H d ,Q,ũ R ,d R ,L,ẽ R ) are real. The soft SUSY-breaking mass parameters of the gauginos, M 1 , M 2 , M 3 , for the bino, the winos, and the gluinos,B,W i (i = 1, 2, 3), andG, corresponding to the weak hypercharge U (1), the weak isospin SU (2), and the color SU (3) symmetry, are in general complex. The R-symmetry can be exploited to choose either M 1 or M 2 to be real. In this paper we keep both M 1 and M 2 complex.
Expanding the scalar Higgs fields about their vacuum expectation values (VEVs) v u , v d , and v s , two further phases, ϕ u and ϕ s , are introduced which describe the phase differences between the VEVs, (v s + h s + ia s ) .

(2.5)
For vanishing phases, the fields h i and a i with i = d, u, s correspond to the CP-even and CP-odd part of the neutral entries of H d , H u and S. The charged components are denoted by h ± d,u . In this paper, we set the phases of the Yukawa couplings to zero. We furthermore re-phase the left-and right-handed up-quark fields as u L → e −iϕu u L and u R → e iϕu u R , so that the quark and lepton mass terms yield real masses.
After electroweak symmetry breaking (EWSB) the six Higgs interaction states mix and in the basis φ = (h d , h u , h s , a d , a u , a s ) the mass term is given by The mass matrix M φφ is obtained from the second derivative of the Higgs potential with respect to the Higgs fields in the vacuum. The explicit expression of the mass matrix M φφ can be found in Ref. [15]. The transformation into mass eigenstates at tree level can be performed with orthogonal matrices R, R G , where the matrix R G is used first to single out the Goldstone boson G. The tree-level Higgs mass eigenstates are denoted by the small letter h, and their masses are ordered as m h 1 ≤ ... ≤ m h 5 .
The mass term for the charged components of the Higgs doublets in is given by where M W is the mass of the W boson, θ W the electroweak mixing angle, e the electric charge and ϕ λ , ϕ κ the complex phases of λ and κ, respectively. The angle β is defined as Here and in the following we use the short hand notation c x = cos x, s x = sin x and t x = tan x. Diagonalizing this mass matrix by a rotation matrix with the angle β c , for which at Born level β c = β, one obtains the charged Higgs boson mass as The charged Goldstone boson G ± , on the other hand, is massless.
The fermionic superpartners of the neutral Higgs bosons,H 0 d ,H 0 u ,S, and of the neutral gauge bosons,B,W 3 , mix, and in the Weyl spinor basis ψ 0 = (B,W 3 ,H 0 d ,H 0 u ,S) T the neutralino mass matrix M N is given by after EWSB, where M Z is the Z boson mass. The symmetric neutralino mass matrix can be diagonalized by a 5 × 5 matrix N , yielding diag(mχ0 where, in terms of the Pauli matrix σ 2 , The fermionic superpartners of the charged Higgs and gauge bosons are given in terms of the the mass term for these spinors is of the form (2.18) The chargino mass matrix M C can be diagonalized with the help of two unitary 2 × 2 matrices, U and V , yielding The left-handed and the right-handed part of the mass eigenstates arẽ 20) respectively, with the mass eigenstates (i = 1, 2) written as Dirac spinors. In summary, the bilinear terms in the chargino and neutralino mass eigenstates are given by where the left-and right-handed projectors are defined as P L/R = (1 ∓ γ 5 ) /2 and i, j = 1, 2 and k, l = 1, . . . , 5.
The scalar partners of the left-and right-handed quarks are denoted byq L andq R , respectively. The mixing matrix for the top squark is given by while the bottom squark matrix reads The mass eigenstates are obtained by diagonalizing these squark matrices with the unitary transformations with the usual convention mq 1 ≤ mq 2 .
3 The NMSSM at One-Loop Level

The Higgs Sector
For the Higgs sector we follow the mixed on-shell OS-DR renormalization scheme described and applied in Refs. [14][15][16][17]. We do not repeat all details here but quote the most important formulae. There are eighteen parameters entering the Higgs sector at tree level, Note that for the sake of convenience we decompose the complex trilinear couplings A λ and A κ into a real and imaginary part in contrast to Ref. [15] where the absolute values and complex phases were used. It was found in Ref. [15] that the four complex phases ϕ s , ϕ u , ϕ λ and ϕ κ do not need to be renormalized. We verified this statement and will discard them in our renormalization procedure.
In our introduction of the NMSSM Higgs sector in Sec. 2 we have already replaced the U (1) and SU (2) gauge couplings g and g and the VEV v by the physical observables M W , M Z and e. It is convenient to further convert, where possible, the input parameters of the set Eq. (3.27) to parameters that can be interpreted more easily in terms of physical quantities. Thus we trade the three soft SUSY-breaking mass parameters m 2 H d , m 2 Hu , m 2 S as well as Im A λ and Im A κ for the tadpole parameters t φ (φ = h d , h u , h s , a d , a s ). These coefficients of the terms of the Higgs potential V Higgs are linear in the Higgs boson fields and have to vanish, in order to ensure the minimum at non-vanishing VEVs v u , v d , v s , It is debatable whether the tadpole parameters can be called physical quantities, but certainly their introduction is motivated by physical interpretation. In the same way, in a slight abuse of the language, we will call the renormalization conditions for the tadpole parameters on-shell. With the new set of input parameters, we allow for two possible renormalization schemes in our Higgs mass calculation. The difference between the two schemes relates to the treatment of the charged Higgs mass. In the first scheme the charged Higgs mass is an OS input parameter, , tan β, v s , |λ|, |κ|, Re A κ DR , (3.29) while in the second scheme Re A λ is an input parameter renormalized as a DR parameter and the charged Higgs mass is a derived quantity, For the definition of the one-loop OS counterterms we refer the reader to Ref. [15]. The counterterms of the DR parameters do not contribute to the final physical results but they are kept to ensure UV finiteness. We define these counterterms solely in the Higgs sector by requiring that all renormalized self-energies of the Higgs bosons be finite. This is different with respect to the definition in Ref. [15] where the chargino and neutralino sectors were used. The numerical results between the two definitions are identical, however. We renormalize the Higgs fields in the DR scheme as described in Ref. [15] at one-loop level, and in Refs. [16,17] at two-loop order O(α t α s ) and O(α 2 t ), respectively.

The Chargino and Neutralino Sector
The chargino and neutralino sectors are described by fourteen real parameters: Since the first ten of these already appear in the Higgs sector, there remains to define the renormalization conditions for the four parameters There are no physical renormalization conditions to fix the counterterms of the phases ϕ M 1 , ϕ M 2 . It has been found in Ref. [15] that the complex phases of M 1 and M 2 do not need to be renormalized. We verified this statement in our computation. 2 In addition, we have to renormalize the chargino and neutralino fields in order to obtain finite self-energies. In the literature there exist two descriptions for the introduction of wave function renormalization (WFR) constants. In the Espriu-Manzano-Talavera (EMT) description two independent renormalization constants were introduced for incoming and outgoing fermions [19][20][21]. Thanks to more degrees of freedom one can keep contributions arising from absorptive parts of the loop integrals and eliminate completely the mixing self-energies thereby fulfilling the standard OS conditions. However, the hermicity of the renormalized Lagrangian is not satisfied anymore. In the Denner description [22], one WFR constant was used instead. It preserves the hermicity constraint, but the absorptive part of the loop integral must be eliminated. We want to investigate the effect of the absorptive part and therefore apply both descriptions. In the following we will derive the formulae in the EMT method. From these formulae, one can easily obtain the ones in the Denner description. The bare parameters and fields are replaced by the renormalized ones and the corresponding counterterms as Note that we do not need to renormalize the rotation matrices U, V and N because their counterterms are redundant. They always appear in combination with WFR constants. One can therefore always redefine the WFR constants to absorb the counterterms of the rotation matrices, as shown in Ref. [23]. In general, the renormalized self-energiesΣ of the fermions can be cast into the following form, cf. Ref. [22], where Σ ij (p) denotes the unrenormalized self-energy of the transitionχ + i →χ + j , i, j = 1, 2, for the charginos andχ 0 i →χ 0 j , i, j = 1, · · · , 5, for the neutralinos. For the charginos, the tree-level mass matrix M treẽ χ and its counterterm δM treẽ χ are given by and for the neutralinos by In the following, we will discuss the OS conditions for the general fermion fieldsχ i having the tree-level masses mχ i . The renormalized fermion propagator matrix is given bŷ where the renormalized one-particle irreducible (1PI) two-point functionsΓ are related to the renormalized self-energies asΓ The propagator matrix has complex poles at M 2 where Γχ i denotes the decay widths. In the OS scheme we require that the tree-level masses are equal to the physical masses, the mixing terms are vanishing at the poles and that the residues of the propagators are unity 3 , In addition, we require the chiral structure to vanish in the OS limit 4 , Applying the decomposition in Eq. (3.38) and the tree-level relations one obtains the mass counterterms δmχ i = Re (δM treẽ χ ) ii , with the off-diagonal wave function renormalization constants, δZ L/R,ij , δZ L/R,ij , 5,6 ReΣij and ReΓii, respectively. Re means that one takes only the real part of the loop integrals but leaves the couplings unaffected. 4 If we use the Denner OS conditions in Ref. [22] then these relations are automatically satisfied for the real parts only. 5 In the Denner description Σ L/R/Ls/Rs ij are replaced by ReΣ L/R/Ls/Rs ij . 6 Note that the δZ are obtained from δZ by interchanging the m 2 χ of the electroweakinos, but not the mχ. and the diagonal wave function renormalization constants 7 Our results coincide with those given in Ref. [19]. The above wave function renormalization constants have been chosen such that for all fermions the mixing terms are cancelled and the correct propagators are produced at the tree-level mass values. Note that in case we do not have enough parameters to renormalize all fermions on-shell, only some of these fermions satisfy OS conditions, i.e. their tree-level masses are equal to the physical masses. The remaining fermions have loop-corrected masses. This is the case for the electroweakino sector. Given the fact that we have already fixed the renormalization scheme for the Higgs sector, only the two gaugino masses M 1 and M 2 remain to be renormalized, while we have 7 masses (5 neutralino and 2 chargino masses), so that only two of them can be set OS. The remaining 5 particles receive loop-corrections to their masses. At loop level the mixing between fermions is in general not vanishing any more, and the residues of the propagators are not unity. These effects should be taken into account if the loop corrections to the masses are large. This is not the case for the renormalization schemes chosen here, therefore we neglect these effects.
In the chargino and neutralino sector hence only the two counterterms δ|M 1 | and δ|M 2 | remain to be determined. There are 21 different ways to choose two out of the 7 masses for the OS conditions. We will consider here two different schemes. In the first scheme (OS1), we require the masses of the wino-like charginoχ + i and the bino-like neutralinoχ 0 k to be OS. The bino-like neutralino is sensitive to M 1 while the wino-like chargino is sensitive to M 2 . Note that we do not choose the chargino and neutralino by referring to a fixed index order since they may not be sensitive to M 1 or M 2 . This can then lead to numerical instability, as was found in the MSSM [23][24][25][26] and in [7,8] for the NMSSM. We denote the tree-level masses for the neutralinos (charginos) by a small letter mχ0(+) The counterterms for M 1 and M 2 are then given by (3.64) (3.65) In the second scheme (OS2), we use the masses of the bino-like neutralino, denoted byχ 0 k , and of the wino-like neutralino, denoted byχ 0 l , as inputs. The renormalization conditions for their OS masses are given byΣχ This results in the two solutions for the counterterms δ|M 1 | and δ|M 2 |, (3.68) (3.71) For the field renormalization constants of the charginos and neutralinos, we impose the OS conditions for the tree-level masses. Besides the two OS schemes, we will also adopt the DR renormalization scheme for M 1 and M 2 , while for the field renormalization constants we use the OS conditions.

The Squark Sector
We consider here only the third-generation squarks. The results for the first-and secondgeneration squarks are obtained analogously. There are seven parameters to be renormalized in this sector, where A t , A b are complex and the mass terms are real. We denote their corresponding counterterms as and define the squark-mass counterterm matrices as with Mq given in Eq. (2.23) for the stops and in Eq. (2.24) for the sbottoms. The renormalization of the remaining parameters appearing in the squark mass matrices has been specified in the renormalization of the Higgs sector, more specifically see Refs. [14][15][16][17].
We have to renormalize the squark fields in order to make the squark self-energies finite. Here we use both the EMT and the Denner description. For the EMT description we have to introduce two separate WFR constants, one for the particle and one for the anti-particle. We introduce the squark WFR constants for the particles and anti-particles in the mass eigenstate basis as 8 The renormalized self-energies in the mass eigenstate basis are given by (i, j = 1, 2) where we denote by Σq ij the unrenormalized self-energies for theq * i →q * j transition. 9 In the following, we give the OS counterterms. The DR counterterms are then easily obtained by taking only the divergent parts of the corresponding OS counterterms.
Applying the decomposition of the fermionic self-energies as given in Eq. (3.38), the mass counterterm in the OS scheme for the top and bottom quark, respectively, is given by (q = t, b) 10 The OS conditions for the scalar renormalized self-energies are (i, j = 1, 2) 11 In the Denner description, we have δZq = δZ † q . 9 Note that in the real NMSSM the unrenormalized self-energies for theq * i →q * j transition and for theqi →qj transition are identical. They are different, however, in the complex case. 10 In the Denner description, Re is replaced by Re. 11 In the Denner description, ReΣq ii is replaced by ReΣq ii in Eqs. (3.78) and ( where There remain two parameters from the bottom squark sector to be determined, . We choose the OS scheme where the bottom squark with the dominant contribution from the righthanded sbottom, which we denote byb i R , is OS and the mixing between the two bottom squark states vanishes. The three counterterms δm 2 b R , Re δA b , Im δA b are then obtained by solving three linear equations is the same in both the EMT and the Denner description. We use Re in the definition of δYq (q = t, b) so that δAq = (δA * q ) * . The contribution from the imaginary part of the loop integrals is then moved into δZq ij and δZq ij . 13 In the Denner description, where we have introduced the notation δYb for later use. The other bottom squark mass gets loop corrections. 14 Its loop-corrected mass Mb j is obtained by solving iteratively the following equation 15 M We stop the iteration when the difference between two consecutive solutions is less than 10 −5 .
The OS wave function renormalization constants for the squarks are given by 16 where the δYq are given in Eq. (3.86) and Eq. (3.92) for q = t and q = b, respectively. Besides the OS scheme we also provide the option to use the DR scheme for all parameters and the wave function renormalization constants. In the DR scheme, all squarks receive loop-corrections to their masses. We will discuss this issue in Subsection 3.4.

Loop-corrected Higgs Boson Masses and Mixings
Since we use the mixed OS-DR renormalization scheme for the Higgs sector parameters together with the DR scheme for the Higgs fields, all Higgs bosons are mixed and receive loop corrections to their masses. For the evaluation of the loop-corrected Higgs boson masses and the Higgs mixing matrix, we use the numerical results obtained from NMSSMCALC [13]. In this code the two-loop corrected Higgs boson masses are obtained by determining the zeros of the determinant of the two-point functionΓ h (p 2 ) with where m h i are the tree-level masses andΣ h ij (p 2 ) is the renormalized self-energy of the h i → h j transition at p 2 . In NMSSMCALC, we have included in the renormalized Higgs self-energies the complete one-loop contributions with full momentum dependence [14,15] and the two-loop 14 Note that in principle all four masses of the stops and sbottoms can be renormalized on-shell simultaneously by adapting the input parameters appropriately. 15 In the Denner description ReΣb jj is replaced by ReΣb jj . 16 In the Denner description, Re Σq ii in Eq. contributions of O(α t α s ) [12] and of O(α 2 t ) [17] in the gaugeless limit at zero momentum. The loop-corrected masses of the Higgs bosons are then sorted by ascending masses 17 (3.100) We have improved the stability of the determination of the Higgs boson masses in NMSSMCALC by implementing two-point loop integrals with complex momentum. In the old version of NMSSMCALC, in order to take into account the contribution of the imaginary part of the complex momentum we expanded the renormalized Higgs self-energies around the real part of the complex momentum asΣ Note that this was done only for the one-loop correction with full momentum dependence. This expansion is not good when Re p 2 is close to threshold regions in which ∂Σ h ij (Re p 2 )/∂Re p 2 contains threshold singularities. To overcome this problem one can use complex masses for the loop particles or complex momenta. Using complex masses requires the decay widths of the particles. These have to be obtained in an iterative procedure which is very time consuming. We have decided to use the complex momenta and to keep the masses real. We have implemented the two-point loop integrals with complex momenta and therefore do not use any more the mentioned approximation. We have confirmed that the evaluation of the Higgs masses is stable in the singularities region and the differences between Higgs masses using the complex momentum and the expansion in Eq. (3.101), defined as (M expansion In processes with external Higgs bosons finite wave-function renormalization factors Z H have to be taken into account in order to ensure the on-shell properties of these Higgs bosons. The wave-function renormalization factor matrix performing the rotation to the OS states is given by [27] Here prime denotes the derivative with respect to p 2 . 17 We denote loop-corrected Higgs mass eigenstates by capital letters Hi (i = 1, ..., 5).

Loop-corrected Neutralino and Chargino Masses
Within the OS and DR schemes defined for the neutralino and chargino sector, the electroweakinos cannot all be renormalized OS, and there remain neutralinos and charginos that receive loop corrections to their masses. In the following we define our procedure to determine the loop-corrected masses for fermions in the general case where mixing contributions are present 18 . This procedure will be used for both the OS and the DR scheme. To give an intuitive definition of the loop-corrected fermion masses, we express the propagator matrix in terms of left-and right-handed Weyl spinors, ψ D = (ψ L , ψ R ) T . In this basis, the tree-level propagator matrix is given by where σ = (σ 1 , σ 2 , σ 3 ) denotes the three Pauli matrices. The mass matrix is given by The loop-corrected propagator matrix in the Weyl basis (ψ L , ψ R ) T is defined as The loop-corrected mass M is determined from the real pole p 2 = M 2 of the propagator matrix satisfying the equation The solution of Eq. (3.109) is given by which can be solved iteratively. When the fermion is OS, p 2 = m 2 , the above relation is nothing else but the OS condition obtained from Eq. (3.47). For the case of n Dirac spinors, the 1PI two-point function in the basis (ψ 1 The matrices a, b, c, d are 2 × 2 matrices in case of charginos and 5 × 5 matrices for neutralinos. With mχ i generically denoting the mass of an electroweakino with index i, the matrices are given by with i, j = 1, 2 for charginos and i, j = 1, ..., 5 for neutralinos.
The poles of the propagator matrix are the solutions of the equation This is equivalent to [30], In practice, we solve the equation numerically through iteration together with the diagonalization of the mass matrix Mχ = c( to obtain the complex poles. The loopcorrected masses are then obtained from the real parts of these complex poles. This procedure is applied for the calculation of the loop-corrected masses for neutralinos using the OS definition of the neutralino WFR constants. However, for the chargino sector the mass matrix Mχ± contains infrared (IR) divergences at arbitrary momentum. We have implemented the approximation used in Ref. [23] and calculate the loop-corrected chargino masses by using the formula (3.118)

Loop-corrected Squark Masses and Mixings
In our DR scheme, the counterterms contain only the UV divergent parts. For the renormalization of the squark fields we use a modified OS scheme. In the following, we will describe the details of this scheme. We first compute the DR squark WFR constants by taking the UV-divergent parts of the OS counterterms, defined in Eqs. (3.94) to (3.98), with tree-level mass values. Using these DR squark WFR constants, we then compute the loop-corrected squark masses Mq i that are obtained by solving iteratively the equations We have assumed here that the off-diagonal renormalized self-energies vanish and that the residues of the propagators are unity at the loop-corrected masses. This is equivalent to redefining the diagonal squark WFR constants at loop-corrected masses as and the off-diagonal WFR constants as In the above equations the renormalized self-energies are computed with the DR squark WFR constants. We keep also the imaginary part of the two-point loop integrals in the renormalized self-energies. Note that the WFR constants δZq ij (M 2 q k ) will enter the evaluation of the decay width. The diagonal WFR constants δZq ii (M 2 q k ) contain IR divergences evaluated at the loopcorrected masses. These IR divergences will cancel exactly with those arising from the virtual part and the real radiation contributions which are also evaluated at the loop-corrected masses. We have verified that this statement is true for both the EW and the QCD corrections.

Higher-Order Corrections to the Two-Body Decays of the Neutral Higgs Bosons
In this section, we present those two-body decay channels that we have improved by including the missing NLO EW and QCD corrections. These channels are the decays into OS SM fermion pairs, OS massive gauge bosons, into a pair of Higgs and gauge bosons, into chargino or neutralino pairs and into top or bottom squark pairs. We will not discuss decays into gluon pairs, photon pairs or Zγ which can be found in Ref. [13]. For these decays, NLO EW corrections are of two-loop order as the leading order (LO) decay widths are already loop-induced. The inclusion of the NLO EW corrections to Higgs-to-Higgs decays on the other hand have been presented in Ref. [11] and the dominant two-loop corrections of the O(α t α s ) have been provided in Ref. [12].
For our computation we have used several programs. The generation of the amplitudes was done by FeynArts [31,32] using a model file created by SARAH [33][34][35][36]. The output amplitudes were further processed using FeynCalc [37,38] for the simplification of the Dirac matrices and for the tensor reduction. The one-loop integrals were evaluated with the help of LoopTools [39].

Higgs Boson Decays into Fermion Pairs
In order to make use of the published results of higher-order corrections in the literature for CPeven and CP-odd Higgs bosons, it is convenient to write the interaction vertex of the complex NMSSM Higgs boson h i (i = 1, . . . , 5) and quarks as where the scalar and pseudoscalar coupling coefficients for the up-and down-type quarks at tree-level are given by where R ij (i, j = 1, 5) denotes the matrix elements of the mixing matrix rotating the tree-level Higgs gauge eigenstates to the mass eigenstates, see Eq. (2.8).
Following the prescription outlined in our publication [13], we improve the widths of the Higgs boson decays into quark pairs by including the missing SUSY-QCD and SUSY-EW corrections. We decompose the EW corrections into the known QED corrections arising from a virtual photon exchange and a real photon emission and the remaining unknown EW corrections from the genuine EW one-loop diagrams.
The one-loop SUSY-QCD corrections originate from loop diagrams with the exchange of a gluino, while the SUSY-EW corrections stem from loop diagrams with weak gauge bosons W, Z, fermions, Higgs bosons and their superpartners in the internal lines. They are both IR finite quantities. The computation of the Higgs boson decays into a bottom quark pair shows that the bottom quark mass counterterm contains terms proportional to t β . This contribution is large in the large-t β regime and universal. In many cases, this contribution is the leading part of the SUSY-QCD and SUSY-EW corrections and can be absorbed into an effective bottom quark Yukawa coupling. This can be done by using an effective Lagrangian formalism [40][41][42]. In Ref. [13], we have presented the effective bottom Yukawa couplings in the real and complex NMSSM. We do not repeat every detail here but only quote the relevant formulae. In Eq. (4.124) we have given the tree-level scalar and pseudoscalar coupling coefficients appearing in the Feynman rule for the CP-violating Higgs bosons h i to a bottom-quark pair. The Feynman rule for the effective coupling including the leading SUSY-QCD and SUSY-EW corrections [40][41][42][43][44][45][46][47][48][49] (denoted by a tilde to mark the inclusion of the corrections) is also decomposed into a scalar and a pseudoscalar part and reads [13] − The correction ∆ b including the leading SUSY-QCD and SUSY-EW corrections can be cast into the form with the one-loop corrections given by is the top-Yukawa coupling and C F = 4/3. The generic function I is defined as Note that the scale of α s in the SUSY-QCD corrections has been set equal to µ R = (mb 1 + mb 2 + |Mg|)/3, while in the SUSY-EW corrections it is µ R = (mt 1 + mt 2 + |µ|)/3. The strong coupling constant α s is evaluated with five active flavors.
The decay width of the CP-violating NMSSM Higgs bosons H i into qq, including the QCD, SUSY-QCD, QED, EW and SUSY-EW corrections, can then be cast into the form In the numerical analysis presented in Sec. 5 we will use the quantity Γ SEW(+QCD) for the decays into fermion pairs to denote the partial decay widths including the SUSY-EW and (for the quarks) SUSY-QCD corrections, i.e. exactly the partial decay width as defined in Eq. (4.134) with the loop-corrected Γ S and Γ P given in Eqs. (4.135) and (4.136), respectively. In contrast, we will denote by Γ tree the partial decay widths that only include the ∆ b corrections, i.e. Eq. (4.134) but with Γ S and Γ P given by the first lines in Eqs. (4.135) and (4.136), respectively.
Note that in δM rem,S/P , δ S/P sub and δM GZ,mix (which will be explained below) we use the tree-level Higgs couplings g h k qq to the quarks. We take the occasion to remind the reader that tree-level mass eigenstates are always denoted by h i and loop-corrected ones by H i . Unless stated otherwise, this means that we use tree-level couplings for external h i but with loop-corrected masses, and for particles inside loop diagrams we always use tree-level masses and tree-level couplings. We comment on the various terms appearing in Eqs. (4.134), (4.135) and (4.136) one by one. The one-loop QED corrections, denoted by ∆ QED have been known in the SM for a long time, cf. Refs. [50][51][52][53][54][55]. In the limit m q M H i they are given by Ref. [56] 19 , where Q q denotes the electric charge of the quark q. The one-loop SM QCD corrections are similar to the QED corrections, with the replacement of Q 2 q α by (4/3)α s (M 2 H i ). It is well-known that the SM QCD corrections are rather large. The large logarithmically enhanced part has been absorbed in Eq. (4.134) into the running MS quark mass m q (M 2 H i ) at the corresponding Higgs mass scale M H i to improve the convergence of the perturbative expansion. The QCD corrections can be taken over from the MSSM case by adapting the Higgs couplings [50][51][52][53][59][60][61][62][63][64][65][66][67][68]. After subtracting the enhanced part, the remaining QCD correction ∆ QCD reads (4.138) where N F = 5 active flavors are taken into account. In the CP-conserving case we also include top quark induced corrections ∆ S/P t by adding them to ∆ QCD . They can be taken over from the MSSM case and read [50][51][52][53][59][60][61][62][63][64][65][66][67], In the decay into a bottom quark pair, the large universal corrections proportional to O(α s t β , α b t β ) are resummed into the effective bottom Yukawa couplingsg S,P h i bb , given in Eqs. (4.127) and (4.128), while in the decay into top quarks we use the tree-level values of the effective couplings, i.e. with δM counter,S/P denoting the counterterm contributions. Since in the decay into bottom quarks we have resummed the dominant corrections into the effective couplings, in the remaining SUSY-QCD and SUSY-EW correction we have to subtract these corrections to avoid double counting by adding appropriate counterterms. This is taken care of by the last terms in Eqs. (4.135) and (4.136), respectively, to which we will come back below.
The term δM GZ,mix is the sum of the contributions from the mixing of the CP-odd component of the Higgs bosons with the neutral Goldstone boson G and with the Z boson, respectively. We use the tree-level masses for the Higgs bosons in the loops in order to ensure the proper cancellation of the UV-divergent pieces, but we use the loop-corrected Higgs boson masses for the external particles in the evaluation of the wave-function renormalization factors, of the amplitudes, and of the decay widths. It is well-known that the use of loop-corrected Higgs boson masses for the external particles violates the Slavnov-Taylor identity of the amplitude δM GZ,mix [11,27,69]. To restore the gauge symmetry one should use the tree-level masses for both the external and the internal Higgs bosons. This causes a mismatch with the phase-space factor (where we use the loop-corrected Higgs boson masses) and with the evaluation of the other amplitudes. We therefore use the loop-corrected masses for the external particles also in these contributions, which are then computed in the unitary gauge. Note that we apply the same method also for the other decays that contain the contribution δM GZ,mix . In Ref. [11], δM GZ,mix was computed only for the ith tree-level mass eigenstate of the decay H i → qq. This may cause instabilities in the case of large mixing between Higgs boson mass eigenstates at loop-level. We avoid this by multiplying it also with the WFR factor Z H .
The remaining SUSY-QCD and SUSY-EW corrections are computed in the Feynman diagrammatic approach. The corrections consist of the contributions from genuine one-loop diagrams and the counterterms. The counterterms are given by with the expressions for the scalar and pseudoscalar parts, δλ The counterterm δv for the VEV is given by (4.149) 20 Note, that the angle β in the sense of a mixing angle does not get renormalized.
The electric coupling e, the W and Z boson masses and tan β are renormalized according to the Higgs sector. The top and bottom quarks are renormalized OS using both the EMT and the Denner descriptions. The terms being related to left-handed and right-handed OS wave function renormalization constants for the quarks (δZ qL/R ) and anti-quarks (δZ qL/R ) (q = t, b) are 21 The DR wave function renormalization constants for the Higgs bosons are denoted by δZ h i h k . We have checked the UV-finiteness of the SUSY-QCD and SUSY-EW corrections to the decay amplitude.
As mentioned above, in the decay into bottom quarks we have to take care to avoid double counting after resumming the dominant part of the SUSY-QCD and SUSY-EW corrections into the effective bottom coupling. To subtract these contributions we add in the decays into a b-quark pair the following counterterms to Eq. (4.135) and Eq. (4.136), for the SUSY-QCD and SUSY-EW corrections, respectively. For the decays into a top-quark pair, these contributions are (4.153) In the decays into strange quarks we also include the one-loop SUSY-QCD corrections. They are obtained after substituting ∆ b as given in Eq. (4.129) with The decays into charm quarks are treated as the decays into top quarks, with the appropriate replacements.
The decays into lepton final states l = e, µ, τ do not receive QCD corrections. Their SUSY-EW corrected decay width is given by where x l = m 2 l /M 2 H i . The ∆ S,P QED are given by Eq. (4.137) after replacing (Q q , m q ) with (Q l , m l ). Furthermore, we resum the dominant SUSY-EW corrections into the effective couplingsg S,P h i ll .
They are obtained from the effective couplingsg h i bb in Eqs. (4.127) and (4.128) after replacing ∆ b with ∆ l , where ∆ l in the complex case is given by The contributions Γ S,P H i →ll are obtained from the ones given in Eqs. (4.135) and (4.136) after replacing q → l.

Higgs Boson Decays into W + W − and ZZ
We now address the higher-order corrections to the Higgs boson decays into gauge boson pairs. We consider corrections only for on-shell decays. Off-shell decays are still treated at tree-level as done in NMSSMCALC [13]. The one-loop corrected decay amplitude for the decay of a CPviolating NMSSM Higgs boson H i (i = 1, ..., 5) into a pair of massive gauge bosons V = Z, W ± , where µ (k 1 ) and ν (k 2 ) are the polarization vectors of the two gauge bosons with four-momenta k 1 and k 2 , respectively. Note that the GZ, mix contribution vanishes. The tree-level amplitudes for the two final-state pairs are And the tensor structure of the NLO corrections is given by 1L ε µνρσ k 1ρ k 2σ .
1L vanishes in the CP-conserving case. The decay width for the decay H i → V V , including the NLO corrections, is given by where R = 1/2 for V = Z and R = 1 for V = W . The improved tree-level decay width reads with The NLO partial width for the Z-boson pair final state contains only virtual contributions, For the W -boson pair final state it contains both virtual and real radiation contributions, The virtual part is given by The formulae for M 1L and M 1L are quite lengthy and we do not display them explicitly here. The counterterm contributions for V = W and Z, respectively, read For the decay H i → W + W − we have to include to the contribution from the radiation of a real photon in order to get an infrared-finite result. This contribution is given by where The formula is in agreement with the result given in Ref. [70]. Here, we have neglected the arguments of the Bremsstrahlung integrals I j 1 ,...,j 2 i 1 ,...,i 2 (M H i , M W , M W ) for the sake of readability. In terms of the photon momentum q, the W + momentum k 1 and the W − momentum k 2 , these integrals are defined as where k 0 is the four-momentum of H i and i l , j k = 0, 1, 2. The plus signs correspond to k 1 , k 2 , the minus sign to k 0 . Their analytic expressions are given in Ref. [22].
We have checked the UV finiteness of the NLO decay widths of both the H i → ZZ and the H i → W + W − decays. The IR divergence in the decay H i → W + W − , however, is more demanding. At strict one-loop level, i.e. one must use the tree-level Higgs boson mass for the external Higgs boson and the unity WFR factor matrix, the IR finiteness if fulfilled. However, the use of the loop-corrected Higgs masses and the WFR factors Z H ij breaks the IR finiteness, because different orders of perturbation theory are mixed in this case. At tree level there exists a relation between the coupling of the neutral Higgs boson with the charged Goldstone bosons and the coupling of the neutral Higgs boson with the W bosons. Defining the Lagrangian for the interaction between Higgs and Goldstone bosons by it is given by with the tree-level Higgs boson mass m h j . In order to obtain an IR-finite result while using the loop-corrected mass M H i , we chose to modify the coupling g h j G + G − as where the tree-level mass m 2 h j has been replaced by the loop-corrected mass M 2 H i of the decaying Higgs boson H i . We verified that the modification of this coupling ensures IR finiteness while not affecting UV finiteness. The same method has been used in Ref. [71]. While taking the loopcorrected mass for the external decaying Higgs boson ensures compatibility with the observation of a 125 GeV SM-like Higgs boson, this approach breaks gauge invariance, however. For more details on this issue, we refer to an upcoming publication [72].
For the non-SM-like neutral Higgs bosons, the tree-level coupling g h i V V is in general suppressed, in particular in case the Higgs boson with mass 125.09 GeV behaves very SM-like. In this case, the one-loop corrected decay width Γ 1L (H i → V V ) can be even larger than the tree-level improved one Γ tree (H i → V V ). This becomes a problem when the one-loop correction is negative, as then the one-loop corrected partial decay width becomes negative. In this case, we have to include the one-loop squared term, which is formally of higher order. For the decay H i → ZZ, we will include the one-loop squared contribution. 22 In particular, the decay width now is given by contains IR divergences so that we cannot treat it in the same way as in the decay H i → ZZ. Note, however, that the 1L decay width Γ 1L (H i → W W ) can always be divided into three parts that are separately UV and IR finite: the (s)fermion contribution arising from loops containing SM model fermions and their superpartners, the chargino/neutralino contribution from loops with internal charginos and neutralinos and the gauge/Higgs contribution from loops with gauge and Higgs particles. In many cases the dominant contribution is the (s)fermion part. That was also observed in the MSSM case [71,73]. We therefore add to the one-loop corrected decay with H i → W + W − the one-loop squared contribution from the (s)fermion part only. This part is IR finite as it solely involves fermions and sfermions but no photons. Both for the decays into ZZ and into W W we include the respective one-loop squared terms in case the one-loop contribution is larger than 80% of the tree-level decay width.

Higgs Boson Decays into a Z Boson and a Higgs Boson
The one-loop corrected amplitude for the decay of a heavy Higgs boson H i with four-momentum p into a light Higgs boson H j and a Z boson, with four-momenta k 1 and k 2 , respectively, H i (p) → H j (k 1 )Z(k 2 ), can be written as (4.180) 22 Note, however, that the thus obtained result has to be taken with caution. The complete two-loop calculation contributes further terms that might cause the complete two-loop result to differ considerably from the result obtained in the here applied approximation. Moreover, the inclusion of (part of) the two-loop corrections explicitly includes a dependence on the renormalization scheme chosen at one-loop order that would need to be cancelled by transforming the input parameters appropriately so as not to become inconsistent. We still use this approach in order to obtain physical, i.e. positive, partial decay widths and hence physical branching ratios. Since the partial decay width is suppressed here anyway, the effect of the difference between the approximation and the full two-loop result on the branching ratio is expected to be subleading. Still, the code NMSSMCALC will print out a warning to make the user aware of this issue. where The tree-level expression M h i h j Z consists of the genuine one-loop diagram contribution and the counterterm part given by Also here, the contribution from the one-loop diagrams with the transition h i → Z(G), M GZ,mix h i h j Z , is calculated in the unitary gauge. The improved tree-level decay width is given by and the NLO decay width by with the 2-particle phase-space factor Since the formulae for the 1-loop amplitudes are quite lengthy we do not display them explicitly here. Note that, as in the decay into massive gauge bosons, in Eq. (4.185) we also included, keeping in mind the caveat mentioned there, one-loop contributions squared as the one-loop corrections can be large and negative.

Higgs Boson Decays into Neutralinos and Charginos
The couplings of a neutral Higgs boson h i with the electroweakinos can be defined as whereχ stands generically for the neutralinos and charginos and g R h iχkχj = g L h iχkχj * . At tree level, the left-and right-handed coefficients for the Higgs-chargino couplings are given in terms of [13] g L h iχ where i = 1, . . . , 5, j, k = 1, 2, and for the Higgs-neutralino couplings we have with l, m = 1, . . . , 5. The decay width for the decay of a Higgs boson H i into a neutralino pair or a chargino pair including higher order corrections is given by where R = 1/2 for identical final states and R = 1 otherwise. The improved tree-level decay width reads where the 2-body phase space factor is 192) in terms of the improved tree-level amplitude The one-loop decay width for the decay into a neutralino pair is given by and for the decay into a chargino pair it is where the virtual contribution can be cast into the form with the left-and right-handed one-loop amplitudes containing genuine triangle, counterterm and 'GZ, mix' contributions, We do not display explicitly here the lengthy expressions for the triangle and 'GZ, mix' contributions. The explicit expressions of the counterterm amplitudes for the decays into neutralinos are and for the decays into charginos they read The right-handed counterterm amplitudes are equal to the complex conjugate of the corresponding left-handed parts after interchanging the indices of the charginos and neutralinos in the final state. The real photon contribution for the decays into a chargino pair is expressed in terms of the Bremsstrahlung integrals as where the arguments of the Bremsstrahlung integrals I j 1 ,...,jm i 1 ,...,in (M H i , Mχ+ j , Mχ− k ) have been neglected. Note that we use the loop-corrected masses for the external Higgs boson and the external charginos and neutralinos in the tree-level, virtual and real contributions. However, for particles inside loops we use the tree-level masses and tree-level couplings. This does not affect the UV-finiteness but can break the IR-finiteness in the decay into a pair of charginos. We overcome this problem by replacing the tree-level mass of the chargino in the loop diagrams with a photon by the corresponding loop-corrected chargino mass. Our treatment is different from Ref. [9] where the authors define an IR divergent counterterm to cancel the mismatch between the real and virtual contributions. Note finally, that in case the NLO decay width into neutralino final states becomes negative, the improved tree-level decay width is calculated instead in NMSSMCALCEW, including the Z H factor.

Higgs Boson Decays into Squark Pairs
The NLO corrections to the decay of a neutral Higgs boson into a squark-antisquark pair consist of the QCD and EW corrections. In the CP-conserving NMSSM, the NLO corrections to the decay of a CP-odd Higgs boson into a stop pair have been calculated and discussed in Ref. [6]. We extend this computation to the CP-violating case and include also the decay into a sbottomantisbottom pair in this paper. The NLO QCD corrections are positive and large. They can be larger than 100% as observed in Ref. [6] while the EW correction are negative and can be of up to −40%. In our calculation, we have implemented both the OS and the DR scheme. We have three options here. First, the seven parameters are renormalized in the OS scheme. Second, the parameters of the stop sector, m t , mQ 3 , mt R , A t , are renormalized in the OS scheme while the remaining parameters, m b , mb R , A b , are renormalized in the DR scheme. Third, all parameters are renormalized in the DR scheme. The loop-corrected decay width is decomposed into the improved tree-level, one-loop QCD and one-loop EW decay widths, (4.201) Denoting the color factor by N F , with N F = 3, the improved tree-level decay width is given by in terms of the improved tree-level amplitude The tree-level Higgs-squark-squark couplings are given by [13] with (4.207) The one-loop QCD and EW contributions to the decay width are given by the sum of the virtual and real contributions, respectively, For the virtual QCD contribution we have with the 2-body phase space factor R 2 defined in Eq. (4.192). The expression for the virtual EW contribution is different from the QCD one due to an extra contribution containing the transition h i → G, Z. Explicitly, we have The explicit expressions for the counterterm contributions are quite lengthy and given in Appendix A. We do not display, however, the more cumbersome amplitudes of the virtual QCD and EW contributions, M ∆,QCD The real photon radiation contribution in the EW corrections is expressed in terms of the Bremsstrahlung integrals as As usual, we have neglected the arguments of the Bremsstrahlung integrals I l (M 2 H i , M 2 q j , M 2 q k ) and I lm (M 2 H i , M 2 q j , M 2 q k ) (l, m = 1, 2 and j, k = 1, 2). The real gluon radiation contribution in the QCD corrections can be obtained from the EW real photon radiation contribution by replacing Q 2 q α with C F α 2 s , where C F = 4/3 for SU (3) C . We have checked the UV and IR finiteness of the EW and QCD corrections. We have compared numerically with the NLO EW and QCD corrections in the OS scheme for the decay A 2 →t 1t2 [6] using their description in the real NMSSM and found full agreement.

Numerical Results
To illustrate the importance of the higher-order corrections to the decays of the light and heavy neutral Higgs bosons and to test the stability of the NLO results in various regions of the parameter space we have performed a scan in the NMSSM parameter space. The parameter points are checked against compatibility with the experimental constraints from the Higgs data by using the programs HiggsBounds5.3.2 [74][75][76] and HiggsSignals2.2.3 [77]. These programs require as input 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 [13,78]. One of the neutral CP-even Higgs bosons is identified with the SM-like Higgs boson -it will be called h from now on -and its mass is required to lie in the range 123 GeV ≤ m h ≤ 127 GeV . Concerning the NMSSM sector, we follow the SUSY Les Houches Accord (SLHA) format [81] in which the soft SUSY breaking masses and trilinear couplings are understood as DR parameters at the scale This is also the renormalization scale that we use in the computation of the higher-order corrections. Note that we chose the charged Higgs boson mass as an OS input parameter. The computation of the O(α t α s + α 2 t ) corrections to the Higgs boson masses is done in the DR renormalization scheme of the top/stop sector. We have included the contribution of the gauge parameters g 1 , g 2 into the conversion from pole to DR top masses. In Table 1 we summarize the ranges applied in our parameter scan. In order to ensure perturbativity we apply the rough constraint The remaining mass parameters of the third generation sfermions that are not listed in the table are chosen as The mass parameters of the first and second generation sfermions are set to We have performed two scans. In the first (smaller) scan we took care to select only such  scenarios where the lightest CP-even Higgs boson H 1 is singlet-like and the second lightest CPeven Higgs boson is the SM-like Higgs boson. We refer to this scan as scan1 in the following. In the second (larger) scan, called scan2 in the following, we only retained scenarios where the SMlike Higgs boson is the lightest CP-even Higgs boson. Both scans allow for points that have a χ 2 computed by HiggsSignals-2.2.3 that is consistent with an SM χ 2 within 2σ. All the branching ratios shown in the following have been calculated by implementing the here presented higherorder corrections to the various decay widths in NMSSMCALC. In this way the new EW corrections are combined with the state-of-the-art higher-order QCD corrections already implemented in NMSSMCALC. Note, however, that the EW corrections are only taken into account if the respective decay is kinematically allowed. Otherwise, the corresponding decay width without the higherorder corrections discussed in this paper, which only apply for on-shell decays, are taken into account in the computation of the total decay width and branching ratios.

Decays into SM Fermion Pairs
In the old implementation in NMSSMCALC the tree-level couplings entering the various decay widths were improved by including loop effects in the Higgs mixing matrix elements. Thus, the tree-level rotation matrix R was replaced by the loop-corrected rotation matrix evaluated at zero external momentum both at one-loop and at two-loop order to ensure unitarity. 23 The implementation here differs by the fact that in the computation of Z H we include the momentum dependence at one-loop order and we do not apply the approximation of Ref. [14] to deduce Z H but proceed as described in Eqs. (3.102)-(3.104). In the following, we call the couplings where we apply R l as obtained from Eq. (5.218) with zero external momentum and by applying the approximation of Ref. [14] 'effective tree-level couplings' while those with Z H calculated according to Eqs. (3.102)-(3.104) including the momentum dependence at one-loop order are denominated 'improved couplings'.
The decays into SM fermion pairs in the old implementation in NMSSMCALC were calculated using the loop-corrected rotation matrix, R l , evaluated at zero external momenta and by including the ∆ b corrections 24 into the effective tree-level couplings, as specified in Ref. [13]. The thus obtained 'effective couplings' are given by 25 g eff,S/P h j qq =g S/P h j qq , . (5.219) Beyond the ∆ b approximation no further SUSY-EW nor SUSY-QCD corrections were included. To quantify the difference between the branching ratio computed in this paper and the old implementation in NMSSMCALC we introduce the relative change in the branching ratio for the decay H i → X j X k as , The branching ratio in the old implementation in NMSSMCALC is denoted by BR tree R l (H i → ff ) (although it also includes the ∆ b corrections where applicable). The quantity ∆ BR hence gives information on the importance of the improvement of the branching ratios by the Z H factor and the SEW(+SQCD) corrections. This quantity will also be used in the investigations of the decays into gauge boson pairs and into a pair of Z and Higgs bosons. The SM-like Higgs boson is given by the h u -like Higgs state 26 , and in our scans we found valid scenarios where this can be the lightest or the second lightest of the CP-even Higgs bosons. We first consider only the parameter points where the lightest CP-even Higgs boson H 1 is the singlet-like state, i.e. has a large h s component. These are points obtained in the above described scan1. Here and in the following we denote a Higgs boson H i to be dominantly x-like (x = h s , h u , h d , a s , a d ) if the corresponding matrix element squared |R l ix | 2 exceeds 80%. When H 1 is h s -like, the question which final state has the largest decay width strongly depends on the amount of admixture of h d and h u components to the singlet-state. In Fig. 1 we show, for all parameter points that pass our constraints, the scatter plot of the H 1 branching ratios into SM particles against its h d component represented by the element R l 1h d of the loop-corrected Higgs rotation matrix 27 . The mass of H 1 lies between 70 and 118 GeV for these points. As can be inferred from the plot, the dominant decays are those into bb, cc, ττ and gg. In most cases the branching ratio into a bottom-quark pair is dominant followed by the decay into τ τ . However, when the h d component of H 1 is very small the branching ratios into gg and cc become competitive and can even be larger than those for the decay into bb with values beyond 60% for the gg final state and of up to 35-39% for cc in some of the scenarios. In this case, i.e. for |R l 1h d | < ∼ 0.02, also the branching ratios into γγ and into the off-shell final state W + * W * − increase and can reach up to about 30% in the latter and about 2% in the former case. The branching ratio into the off-shell Z * Z * final state, which also increases then, is about one order of magnitude smaller than the one into W + * W − * . But already for |R l 1h d | > ∼ 0.02 the decay into bb takes over again and reaches branching ratio values of up to 90% followed by the branching ratio into τ τ with values of up to 10%.
In order to investigate the importance of the higher-order corrections we define for our new implementation the relative correction of the partial width for the decay H i → X j X k as with the higher-order decay widths for the decays H i → qq into quarks given in Eq. (4.134) and the higher-order decay widths for the decays H i → ll into leptons given in Eq. (4.155) and with the tree-level decay width Γ tree Z H including only the ∆ b corrections. The tree-level and higher-order decay widths are both evaluated with the new implementation of Z H . Note that the quantity δ gives information on the importance of the SEW(+SQCD) corrections in the decay width alone as the factor Z H cancels in the ratio. In Fig. 2 we show the scatter plot of the relative change of the branching ratios, ∆ BR (H 1 ff ), f = b, c, τ , for all the parameter points passing the constraints, against BR The color code in Fig. 2 as well as in Figs. 3-6 denotes the sizes of the relative corrections of the partial decay widths. The points where the absolute value of δ SEW(+SQCD) exceeds 10% are marked in blue, those with δ SEW(+SQCD) in the [5,10]% range in black and those with relative corrections less than 5% in red. For Figs. 5b and 6b we distinguish two regimes for the larger corrections, in blue where δ SEW+SQCD is in the [10,20]% range and in green where δ SEW+SQCD is in the [20,40]% range; in Fig. 4b we also add two other categories of points, in cyan where δ SEW+SQCD is in the [40,60]% range and in pink where δ SEW+SQCD is in the [60,80]  implementation ranges below about 30% with vertex corrections |δ| smaller than 5%. There are some very rare scenarios where ∆(H 1 bb) exceeds 50% and where at the same time the relative vertex corrections are between 5 and 10%. We investigated these cases and observed that there is an accidental cancellation either in the effective tree-level or in the improved couplings. These parameter points lead to similar results for the τ final states, i.e. |∆ BR (H 1 ττ )| > 50% and at the same time |δ| between 5 and 10%. The cancellation results in a suppression of the branching ratio, to less than 4% for the ττ final state and 10% at maximum for the bb final state. In most of the cases, the large ∆ BR (H 1 ff ) is due to the use of the wave-function renormalization factor Z H , however. There are also cases with a cancellation between the SUSY-EW/SUSY-QCD corrections and the wave-function renormalization factor Z H correction. This results in ∆ BR (H 1 cc) being less than 1%.
We have performed the same analysis for the heavier Higgs bosons, using the full set of points from our scan2. Figure 3 is the scatter plot of the ∆ BR (hff ), f = b, τ, c, t, against BR for a heavy a s -, a-and h d -like Higgs boson i.e. H i = H as , H a , H h d , respectively. As before, the SM-like Higgs boson h is always h u -like and decays dominantly into a bottom-quark pair with a branching ratio of about 60%, as expected, followed by the decay into a τ pair and the decay into a c-quark pair. As can be inferred from Fig. 3, the relative changes between the old and the new implementation are much smaller than for the singlet-like lightest Higgs boson and amount only to a few percent. The relative vertex corrections |δ| are below 10% for the b-quark pair final state and below 5% for the decays both into ττ and cc.
For the heavy Higgs bosons, the decay into a top quark pair can become kinematically possible. We start by discussing the decay pattern of the heavy singlet state H as , with a mass between 120 GeV and 1.7 TeV, into the b-quark and t-quark final states, presented in Fig. 4. 28 [%]

∆BR(H as bb)
[%] [%] For the bb final states the relative change in the branching ratios due to the new implementation is mostly between -20% and 20%. We also find points where the relative change is close to 100%, in particular for branching ratios close to 100%. Most points exhibit small relative vertex corrections (see, red points in Fig. 4a), so that the large changes of ∆ BR (H as bb) are due to the implementation of Z H . This is especially the case for large branching ratios close to 100%. There are a few points where the relative vertex corrections lie between 5 and 10% (black) and even above 10% (blue). This happens for the cases where the effective tree-level couplings H as ff are suppressed. The relative change ∆ BR can still be very small when the effects from the Z H factor and the vertex corrections cancel. The decay pattern for the tt channel, finally, is displayed in Fig. 4b. The branching ratio takes all values between almost 0 and 100%. The relative changes ∆ BR (H as tt) are mostly between -20% and 20% and close to 0% for large branching ratios above about 60%. We also observe large ∆ BR , in particular for branching ratios close to into ττ and cc. [%] [%] zero. This is mainly due to very suppressed effective tree-level H as tt couplings corresponding to the regions where the branching ratio BR(H as → bb) is enhanced. These regions correspond to large values of tan β close to the upper bound in our scan, or to smaller mass values of H as with not sufficient phase space to decay into an on-shell top-quark pair. In these regions the relative corrections |δ| are most of the time below 40%, and for cases where |δ| < 5% the large changes in ∆ BR are mostly due to the use of the wave-function renormalization factor Z H . For larger branching ratios the relative corrections |δ| are mostly below 10% (red and black points). Some rare scenarios display corrections above 40% and up to 80% (in cyan and in pink), again mostly in regions with lower branching ratios.
Similar observations can be made for the other heavier Higgs states H a (with a mass between 539 GeV and 2 TeV) and H h d (with a mass between 548 GeV and 2 TeV), with the notable exception that the relative changes ∆ BR (H a/h d X i X j ) are more reduced and never reach 100%. The relative changes ∆ BR (H a/h d bb) are most of the time positive and below 40% as seen in Figs. 5a and 6a. The decays into top-quark pairs can be dominant where the decays into bb are suppressed, and the relative changes ∆ BR between the old and new implementation are close to zero when BR(H a/h d → tt) → 100% as seen in Figs. 5b and 6b. For some rare scenarios the relative vertex corrections |δ| can reach 40%, depicted in green in the figures. Note that BR(H a/h d → bb) can reach 90%, corresponding to regions where the effective tree-level coupling H a/h d bb is strongly enhanced due to large values of tan β while at the same time the effective tree-level coupling H a/h d tt is strongly suppressed.

Decays into a Massive Gauge Boson Pair
In the CP-conserving case, the heavy Higgs boson that can decay into two on-shell massive gauge bosons is h d -like. The tree-level coupling of a Higgs boson H i to V V (V = W, Z) is proportional to [%] [%] because of sum rules. This leads to very suppressed tree-level partial decay widths Γ(H h d → V V ).
In order to compare the results obtained in this paper with the old implementation in NMSSMCALC using the tree-level coupling together with the loop-corrected rotation matrix R l , we show in Fig. 7a the relative change ∆ BR (H h d W W ) of the branching ratio into W W between the old and the new implementation including the NLO-EW vertex corrections as described in Sec. 4.2 and the improvement with the Z H factor, as a function of the loop-corrected branching ratio BR SEW Z H (H h d → W + W − ). The plotted points are those of our scan that pass the constraints we have applied. We display in Fig. 7b the same but for the decay into ZZ. The color and symbol code denotes the magnitude of the relative NLO electroweak vertex corrections alone,

∆BR(H hd ZZ)
[%] with δ defined in Eq. (5.221). As can be inferred from the plots both the Z H factor and the NLO electroweak corrections can be responsible for the large relative changes in the branching ratios. This is in particular reflected by the black points for which the vertex corrections are below 20% and at the same time the relative changes ∆ BR can reach up to 100%. In some cases there is a cancellation between the two contributions (the vertex corrections and the Z H factor) leading to relatively small relative changes in the branching ratios. These cases are the pentagon-marked full (red) points in Fig. 7a for the decay H h d → W + W − , which are mostly located in regions where |∆ BR (H h d W W )| < ∼ 25% while the relative vertex correction |δ| is at least 100%. In the case of the decay into a Z boson pair, however, the bulk of these pentagon-marked full (red) points, indicating again a relative vertex correction |δ| of at least 100%, induces large relative changes of ∆ BR (H h d ZZ) close to 100%. These points correspond to a region which is discussed in more detail in the next paragraph.
We also note that there are two regions concentrating many points for the decay into a Z boson pair, the region for which ∆ BR (H h d ZZ) 0% and the one for which ∆ BR (H h d ZZ) 100%. This is in contrast to the decay into a W boson pair which is mostly centered around ∆ BR (H h d W W ) 0% for vertex corrections |δ| < 80% and much more scattered for the points where 80 < |δ| < 100%, displayed with cyan-pentagon-marked points (for |δ| ≥ 100% the above described cancellation takes place in the decay H h d → W + W − ). This presence of the second region in the Z boson final state, for which ∆ BR is close to 100% can be explained by the occurrence of many parameter points having a very suppressed tree-level coupling H h d ZZ.
They also correspond to regions where the loop-corrected partial decay width Γ(H h d → ZZ) is higher, up to 1 GeV, while the decay width is at most 5 MeV for the region centered around ∆ BR (H h d ZZ) = 0. Note, that while the tree-level couplings H h d ZZ and H h d W W are the same, the loop-corrected decay widths differ by the fact that the decay into W W bosons receives real corrections and that in the one-loop squared contributions to the decay width H h d → W + W − we only include the (s)fermion contributions in contrast to the decay H h d → ZZ. 29

Decays into a Z Boson and a Higgs Boson
In the searches for heavy pseudoscalars, this decay can be an important search channel [82,83]. We are interested here in how large the branching ratio can be and how important are the newly included higher-order corrections, in the case for which on-shell decays are possible. With the obtained valid set of parameter points, we present in Fig. 8a scatter plots of the relative changes of the branching ratios between the old and new implementation for the decay of a heavy pseudoscalar-like Higgs boson H a into ZH 1 and in Fig. 8b for the decay into ZH 2 , against the respective loop-corrected branching ratios. The color and marker codes denote the relative sizes of the one-loop vertex corrections within specific ranges, identical to the ranges used in the previous sub-section for the decays into gauge boson pairs. The mass values of the individual involved CP-even Higgs bosons in the final state range in 123 GeV ≤ m H 1 ≤ 127 GeV and 463 GeV ≤ m H 2 ≤ 1.73 TeV, while the mass of the decaying Higgs boson ranges in 539 GeV ≤ m Ha ≤ 2.0 TeV for the on-shell decay into ZH 1 pairs and in 713 GeV ≤ m Ha ≤ 2.0 TeV for the on-shell decay into ZH 2 pairs. As can already be inferred from the mass ranges, the H 1 state is the SM-like Higgs boson h, while the H 2 state is the singlet-like scalar Higgs boson H hs .
We observe that the branching ratios into the ZH 1 final state remain very small, below 0.4%, while those of the decay into ZH 2 can reach 11%. This is due to the nature of the H 1 Higgs boson [%] [%] that is SM-like, with very suppressed tree-level H 1 H a Z couplings. The relative changes ∆ BR of the branching ratio for the decay H a → ZH 1 are mostly between 0 and -75% corresponding to relative vertex corrections |δ SEW | being at most 60% (black, blue, and pink points), while a few points corresponding to higher vertex corrections up to more than 100% (green, cyan, and red points) can reach ∆ BR = ±100%. These extreme points correspond to very small values for the branching ratios themselves which explains in turn the very large relative corrections |δ| that we observe. Note that in these decays we take into account the one-loop squared term as described in Eq. (4.185) which makes up for the main contribution to the very large relative corrections. As for the decay H a → ZH 2 the relative correction ∆ BR is most of the time between 0 and ±25%, corresponding to points where the relative vertex corrections |δ| are below 20%. A few points display larger ∆ BR values, and also larger relative vertex corrections |δ| that can reach 100% and even beyond, again for points that display very small branching ratios, below about 10 −4 %. Note that there are points for which the correction |∆ BR | is rather limited, below 25%, while the relative vertex correction can reach 40% (for one scenario even more than 100%). This can be explained by a sign compensation between the Z H factor and the vertex correction.
The corresponding results for the heavy singlet-like Higgs boson H hs decaying into ZA 1 and ZA 2 are shown in Figs. 9a and Figs. 9b, respectively. We see that the maximum achieved branching ratios for the decay H hs → ZA 1 are below 20% and for a large number of parameter points are tiny. In the case of the decay H hs → ZA 2 the branching ratio can reach around 15%. In most of the cases A 1 is singlet-like, corresponding to points where the branching ratio is small (below 10%), while A 2 is doublet-like. The mass values of the individual involved Higgs bosons in the final state range in 120 GeV ≤ m A 1 ≤ 1.50 TeV and 562 GeV ≤ m A 2 ≤ 1.63 TeV, while the mass of the decaying Higgs boson ranges in 464 GeV ≤ m H hs ≤ 1.75 TeV for the on-shell decay into ZA 1 pairs and in 696 GeV ≤ m H hs ≤ 1.75 TeV for the on-shell decay into ZH 2 pairs. The cases with larger branching ratios for the decay H hs → ZA 1 (larger than 10%) correspond mostly to the few A 1 pseudoscalar Higgs bosons with doublet-like admixture and mass values above 400 GeV. [%] [%] Figure 9: Scan2: Same as in Fig. 8 but for a heavy singlet-like Higgs boson H hs decaying into ZA1 and ZA2 pairs. In the right plot, however, red: |δ SEW | < 5%; black: 5 ≤ |δ SEW | < 10%.
The relative changes ∆ BR in the branching ratios are mostly positive, and can reach values of 100%. For some very rare scenarios we get ∆ BR (H hs ZA 1 ) close to -100%, while ∆ BR (H hs ZA 2 ) is not below -15%. The relative vertex corrections are moderate for the decay H hs → ZA 1 , mostly |δ| ≤ 20% (black points). This means that the large changes in ∆ BR (H hs ZA 1 ) are mostly due to the Z H factor. For very small branching fractions below 10 −4 % larger vertex corrections are possible, mainly because the denominator in the definition of δ is very small in these regions and can lead to sharp changes in δ. Note that the bulk of the changes between the old and the new implementation in these cases stems from the vertex corrections. In the case of the decay H hs → ZA 2 the relative vertex corrections are mostly small, with values |δ| < 5% (red triangle-marked points). For large relative changes ∆ BR (H hs ZA 2 ) the Z H factor is responsible then.

Decays into Charginos and Neutralinos
We start by investigating the loop corrections to the masses of the charginos and neutralinos using the three renormalization schemes OS1, OS2 and DR, imposed on the two gaugino masses M 1 and M 2 , as defined in Section 3.1.2. According to the SLHA format that we apply in our code, M 1 , M 2 are DR parameters given at the scale M SUSY = √ mQ 3 mt R . When we use the OS schemes we have to translate the DR input parameters to the OS values by applying the approximate transformation formulae      Table 2.
For scenario1, we present in Table 3 the tree-level and loop-corrected masses of the neutralinos and charginos in the three different renormalization schemes and for the Denner description. As expected the wino-like neutralino and the wino-like chargino which couple to the electroweak gauge bosons, get significant loop corrections in the DR scheme. The one-loop corrected masses themselves, however, barely differ in the three renormalization schemes so that the remaining theoretical error due to missing higher-order corrections is very small.
We vary the phases of the gaugino mass parameters M 1 and M 2 , ϕ M 1 and ϕ M 2 , in order to study their effect on the loop-corrected neutralino and chargino masses. Note that these complex phases have negligible impact on the Higgs sector [15,17]. We use the DR scheme for    In order to study the loop corrections on the decay widths, we computed the tree-level and loop-corrected decay widths, defined in Section 4.4, for the three different renormalization schemes OS1, OS2 and DR for the scenario1. Note that we use the loop-corrected masses for external Higgs bosons, charginos and neutralinos not only in the loop-corrected decay widths but also in the tree-level ones. For illustration, we present in Table 4   Γ SEW Z H (H i →χ jχk ), the relative corrections δ(H iχjχk ) as defined in Eq. (5.221), and the loopcorrected branching ratios BR l ≡ BR SEW Z H (H i →χ jχk ) using scenario1, for the OS1 and DR renormalization schemes. We found that Γ l is almost identical in the two OS schemes. The relative size of the loop corrections is below 10% in the OS scheme. The relative corrections in the DR scheme are always larger than in the OS schemes. For some channels with small decay widths, we see significant corrections in the DR scheme. For example in the decay channel H 4 →χ 0 2χ 0 2 , the relative loop correction is −0.13% in the OS scheme while it is −20% in the DR scheme. Based on our investigation, we conclude that it is better to use the OS scheme in the decays of the neutral Higgs bosons into electroweakinos. The largest uncertainty due to missing higher-order corrections that we estimate from the variation of the renormalization schemes is found to be 17% in the decay H 5 →χ 0 1χ 0 2 . We also studied the difference between the Denner and EMT descriptions and did not observe any significant difference. Defining the difference as with Γ x being the loop-corrected decay width for some Higgs decay into an electroweakino final state, we see that the differences are of per mille level for all investigated decays.

Decays into a Squark Pair
We start by discussing the top and bottom squark masses in the OS and DR renormalization schemes defined in Section 3.1.3. We follow the SLHA convention where the input parameters mQ 3 , mt R , mb R , A t , A b are DR parameters at the scale M SUSY . When we apply the OS scheme, these parameters must be translated into OS parameters by applying the approximate transformation formula (i = 1, 2) with X = mQ 3 , mt R , mb R , A t , A b and the finite part of their OS counterterms denoted as δX finOSi . We have used an iterative method to obtain a stable value of δX finOSi . Note that we include both the NLO QCD and the full NLO EW contribution in the conversion Eq. (5.228). 31 Using these OS parameters together with the OS top mass in the tree-level mass matrices, we obtain the top and bottom squark masses. When we apply the DR renormalization scheme, the top pole mass has to be translated to the DR top mass for which we follow the description in appendix C of Ref. [17]. The DR top and bottom masses together with the DR parameters mQ 3 , mt R , mb R , A t , A b are then used in the tree-level mass matrices to get the tree-level rotation matrices. They are subsequently used in the computation of the renormalized self-energies of the top and bottom squarks to obtain the loop-corrected squark masses as described in Section 3.4. In principle, we would expect that the loop-corrected masses computed in the DR scheme are closer to the OS masses in the OS description if one includes more higher order corrections. We consider here a parameter point (scenario2) given by the following soft SUSY breaking masses and trilinear couplings 32  tree  334  1166  1122  2297  1loop  334  1166  1125  2297   DR  tree  585  1216  1181  3000  1loop  341  1152 1152 2294  with the CP phases given by With the given DR parameters of the squark sector, we obtain their corresponding OS parameters mQ 3 = 1120 GeV , mt R = 402 GeV , mb R = 2997 GeV , A t = 1720 GeV , A b = −581 GeV . (5.232) The tree-level and loop-corrected masses of the stops and sbottoms in the OS and DR scheme are shown in Table 5. We see that for the DR scheme there are large changes between the tree-level and one-loop masses, in particular for the lightest stopt 1 . The loop-corrected masses, however, are then closer to each other in both schemes, as expected. The maximum difference if found for the light sbottomb 1 mass, where the OS and DR results differ by 27 GeV at one-loop order (compared to 59 GeV at tree level). The two-loop corrected neutral Higgs boson masses at O(α t α s + α 2 t ) together with their respective main component are displayed in Table 6. The SM-like Higgs boson mass is around 124 GeV while the remaining Higgs spectrum is quite heavy with masses above 1.6 TeV.
We now turn on the complex phase of A t . In the left plot of Fig. 11, we show the tree-level (black), NLO EW (green), NLO QCD (blue), and full, i.e. NLO QCD+EW, (red) corrections to the partial width of the decay H 4 →t 1t * 1 as function of the phase ϕ At in both the OS (full lines) and the DR (dashed lines) schemes while their corresponding branching ratios are depicted in the right plot. The decay H 4 →t 1t * 1 vanishes in the CP-conserving limit where H 4 , which is a-like in scenario2, is a CP-odd Higgs boson. (Note that CP-odd Higgs bosons at tree-level only couple to two different stops.) In the OS scheme, the relative EW corrections δ (see Eq. (5.221) for the definition) vary in the range (-6%,-10 %) and the QCD corrections in the range (−4%, −8%) depending on the phase ϕ At that is varied from zero to ±π/2. In the DR scheme, the relative EW corrections are of order −10% and the relative QCD corrections of 22% and depend slightly on the phase ϕ At . We define the relative differences between the OS and DR decay widths and branching ratios, as respectively. For ϕ At = −π/2, the relative difference ∆ Γ of the partial decay widths between the OS and DR schemes is then about 40% at tree-level while it reduces to 4% when both QCD and EW corrections are included, so that at one-loop level we clearly see a reduction of the theoretical error due to missing higher-order corrections For the relative error in the branching ratios we find values between 32% and 27% at tree level and between 0.3% and 3% at one-loop order including both the EW and QCD corrections while the phase ϕ At is varied from ±π/2 to zero.
For the decay H 4 (≡ H a ) →t 1t * 2 , we show in the upper panels of Fig. 12 the partial decay width (left) and branching ratio (right) at tree-level (black), NLO EW (green), NLO QCD (blue), and NLO EW+QCD (red) as a function of ϕ At , both for the OS (full) and DR (dashed) scheme . In the middle panels we show the relative NLO EW, NLO QCD and NLO EW+QCD corrections which are defined as respectively. The lower panels display the relative differences between the OS and DR decay widths and branching ratios, ∆ Γ and ∆ BR , as defined in Eq. (5.233) and Eq. (5.234), respectively.   The corrections vary slightly with the phase ϕ At . The EW corrections are negative in both schemes while the QCD corrections are positive and of the same order of magnitude. This shows the importance to include both types of corrections to make reliable predictions. Overall, the relative corrections δ in the DR scheme are larger than in the OS scheme. As can be inferred from the bottom left panel of Fig. 12, for ϕ At = 0, the relative difference in the partial width, ∆ Γ , between the OS and DR scheme at tree level is about 37% and decreases dramatically to less than 12% when both the EW and the QCD corrections are included. 33 Similar results are found for the branching ratios, presented on the right plots of Fig. 12, with smaller values of 21% at tree level and 7% at full one-loop order. Note that in the right hand side plots we treated the decays H a →t 1t * 2 and H a →t 2t * 1 at the same level of precision while all other decays are computed at the highest possible precision.
In the CP-invariant scenario the decay width of decay H a →t 1t * 2 is equal to the one of its charge conjugate decay H a →t 2t * 1 . For non-vanishing ϕ At , however, the CP asymmetry, defined as is non-zero. In Fig. 13 we show the CP asymmetry as a function of ϕ At . We see that the CP asymmetry appears already at tree-level, which results from the imaginary part of the WFR factor Z H and the imaginary part of the tree-level couplings g h iqjq * k . The relative change of the asymmetry due to loop corrections ranges between 18% and -9% in the OS scheme while in the  DR scheme it is about 8% when the phase ϕ At is varied from −π/4 to π/4.
We present in Fig. 14 the same plots for the decay widths and branching ratios as in Fig. 12 but for the decay of H 3 →t 1t * 2 which is the dominant decay channel of H 3 . In scenario2 H 3 is h d -like. For both the OS and the DR scheme the NLO EW corrections to the decay width are negative and the relative corrections δ Γ are around -3% in the OS and -8% in the DR scheme. The NLO QCD corrections on the other hand are positive and their relative size can reach 8% in the OS scheme and around 37% in the DR scheme. For ϕ At = 0, the difference between the decay widths in the OS and the DR scheme, ∆ Γ , is about 36% at tree-level and reduces to 12% at full one-loop level. The corresponding numbers for the branching ratios are similar.
We have observed that for this parameter point the EW and QCD corrections in the OS scheme are smaller in the DR scheme. In the OS scheme we have seen that this is due to a cancellation between the genuine triangle diagram contribution and the counterterm contribution while this does not happen in the DR scheme. Overall the inclusion of the QCD corrections in addition to the EW corrections reduces the difference between the OS and DR results.

Conclusions
In this paper, we have presented in the framework of the CP-violating NMSSM our calculation of the complete (SUSY-)EW to the on-shell Higgs boson decays into fermion pairs, massive gauge boson final states, gauge and Higgs boson final state pairs, electroweakino and stop and sbottom pairs. Where applicable we have included the SUSY-QCD corrections. We have implemented these new corrections into NMSSMCALC, a Fortran code for the computation of the Higgs mass spectrum up to presently two-loop order O(α t α s + α 2 t ) and the calculation of their branching ratios. The code already included in the branching ratios the state-of-the-art QCD corrections and the ∆ b corrections as well as decays into off-shell massive gauge boson pairs and decays with off-shell heavy quarks in the final state. The new code is called NMSSMCALCEW.
The consistent implementation of our newly computed corrections provides the presently highest level of precision in the calculation of the NMSSM Higgs boson decays. In contrast to the previous NMSSMCALC version, we have included the Z H factor with full momentum dependence in order to render the loop-corrected masses on-shell. For the decays into the electroweakinos and into squark pairs different renormalization schemes were implemented. The numerical analysis has demonstrated that the relative change in the branching ratios due to this new treatment and the newly implemented corrections is significant. The analysis of the decays into chargino/neutralino and squark pairs for different renormalization schemes has shown that the one-loop corrections reduce the theoretical uncertainty due to missing higher-order corrections. The new code NMSSMCALCEW can be obtained at the url: https://www.itp.kit.edu/∼maggie/NMSSMCALCEW/.

A Counterterm contribution to the decay of a neutral Higgs boson into a squark pair
In this appendix, we give explicit expressions of the counterterm couplings for the decays of a neutral Higgs boson into a squark pair. Note that we have used the redefined WFR constants for the squarks given in Eqs. (3.121) and (3.122). For the EW corrections, the counterterm entering Eq. (4.210) for the Higgs decay into a stop pair is given by Ut k2 F * Further information on the organization of the files for the code and their functionalities as well as modifications on the code (which are constantly updated) can be found at the webpage of NMSSMCALCEW. The code has been tested on a Linux machine.