Decays of the neutral Higgs bosons into SM fermions and gauge bosons in the $\mathcal{CP}$-violating NMSSM

The Next-to-Minimal Supersymmetric Standard Model (NMSSM) offers a rich framework embedding physics beyond the Standard Model as well as consistent interpretations of the results about the Higgs signal detected at the LHC. We investigate the decays of neutral Higgs states into Standard Model (SM) fermions and gauge bosons. We perform full one-loop calculations of the decay widths and include leading higher-order QCD corrections. We first discuss the technical aspects of our approach, before confronting our predictions to those of existing public tools, performing a numerical analysis and discussing the remaining theoretical uncertainties. In particular, we find that the decay widths of doublet-dominated heavy Higgs bosons into electroweak gauge bosons are dominated by the radiative corrections, so that the tree-level approximations that are often employed in phenomenological analyses fail. Finally, we focus on the phenomenological properties of a mostly singlet-like state with a mass below the one at $125\,$GeV, a scenario that appears commonly within the NMSSM. In fact, the possible existence of a singlet-dominated state in the mass range around or just below $100\,$GeV would have interesting phenomenological implications. Such a scenario could provide an interpretation for both the $2.3\sigma$ local excess observed at LEP in the $e^+e^-\to Z(H\to b\bar{b})$ searches at $\sim 98\,$GeV and for the local excess in the diphoton searches recently reported by CMS in this mass range, while at the same time it would reduce the"Little Hierarchy"problem.

The Next-to-Minimal Supersymmetric Standard Model (NMSSM) offers a rich framework embedding physics beyond the Standard Model as well as consistent interpretations of the results about the Higgs signal detected at the LHC. We investigate the decays of neutral Higgs states into Standard Model (SM) fermions and gauge bosons. We perform full one-loop calculations of the decay widths and include leading higher-order QCD corrections. We first discuss the technical aspects of our approach, before confronting our predictions to those of existing public tools, performing a numerical analysis and discussing the remaining theoretical uncertainties. In particular, we find that the decay widths of doublet-dominated heavy Higgs bosons into electroweak gauge bosons are dominated by the radiative corrections, so that the tree-level approximations that are often employed in phenomenological analyses fail. Finally, we focus on the phenomenological properties of a mostly singlet-like state with a mass below the one at 125 GeV, a scenario that appears commonly within the NMSSM. In fact, the possible existence of a singlet-dominated state in the mass range around or just below 100 GeV would have interesting phenomenological implications. Such a scenario could provide an interpretation for both the 2.3 σ local excess observed at LEP in the e + e − → Z(H → bb) searches at ∼ 98 GeV and for the local excess in the diphoton searches recently reported by CMS in this mass range, while at the same time it would reduce the "Little Hierarchy" problem.

Introduction
The signal that was discovered in the Higgs searches at ATLAS and CMS at a mass of ∼ 125 GeV [1][2][3] is, within the current theoretical and experimental uncertainties, compatible with the properties of the Higgs boson predicted within Standard-Model (SM) of particle physics. No conclusive signs of physics beyond the SM have been reported so far. However, the measurements of Higgs signal strenghts for the various channels leave considerable room for Beyond Standard Model (BSM) interpretations. Consequently, the investigation of the precise properties of the discovered Higgs boson will be one of the prime goals at the LHC and beyond. While the mass of the observed particle is already known with excellent accuracy [4,5], significant improvements of the information about the couplings of the observed state are expected from the upcoming runs of the LHC [3,[6][7][8][9] and even more so from the high-precision measurements at a future e + e − collider [10][11][12][13].
Motivated by the "Hierarchy Problem", Supersymmetry (SUSY)-inspired extensions of the SM play a prominent role in the investigations of possible new physics. As such, the Minimal Supersymmetric Standard Model (MSSM) [14,15] or its singlet extension, the Next-to-MSSM (NMSSM) [16,17], have been the object of many studies in the last decades. Despite this attention, these models are not yet prepared for an era of precision tests as the uncertainties at the level of the Higgs-mass calculation [18][19][20] are about one order of magnitude larger than the experimental uncertainty. At the level of the decays, the theoretical uncertainty arising from unknown higher-order corrections has been estimated for the case of the Higgs boson of the SM (where the Higgs mass is treated as a free input parameter) in Refs. [21,22] and updated in Ref. [23]: depending on the channel and the Higgs mass, it typically falls in the range of ∼ 0.5-5%. To our knowledge, no similar analysis has been performed in SUSY-inspired models, but one can expect the uncertainties from missing higher-order corrections to be larger in general-with many nuances depending on the characteristics of the Higgs state and the considered point in parameter space: we provide some discussion of this issue at the end of this paper. In addition, parametric uncertainties that are induced by the experimental errors of the input parameters should be taken into account as well. For the case of the SM decays those parametric uncertainties have been discussed in the references above. In the SUSY case the parametric uncertainties induced by the (known) SM input parameters can be determined in the same way as for the SM, while the dependence on unknown SUSY parameters can be utilised for setting constraints on those parameters. While still competitive today, the level of accuracy of the theoretical predictions of Higgs-boson decays in SUSY models should soon become outclassed by the achieved experimental precision on the decays of the observed Higgs signal. Without comparable accuracy of the theoretical predictions, the impact of the exploitation of the precision data will be diminished-either in terms of further constraining the parameter space or of interpreting deviations from the SM results. Further efforts towards improving the theoretical accuracy are therefore necessary in order to enable a thorough investigation of the phenomenology of these models. Besides the decays of the SM-like state at 125 GeV of a SUSY model-where the goal is clearly to reach an accuracy that is comparable to the case of the SM-it is also of interest to obtain reliable and accurate predictions for the decays of the other Higgs bosons in the spectrum. The decays of the non-SM-like Higgs bosons can be affected by large higher-order corrections as a consequence of either large enhancement factors or a suppression of the lowest-order contribution. Confronting accurate predictions with the available search limits yields important constraints on the parameter space.
In this paper we present an evaluation of the decays of the neutral Higgs bosons of the Z 3conserving NMSSM into SM particles. The extension of the MSSM by a gauge-singlet superfield was originally motivated by the 'µ problem' [24], but also leads to a richer phenomenology in the Higgs sector (see e. g. the introduction of Ref. [25] for a recent summary of related activities). Several public tools provide an implementation of Higgs decays in the NMSSM: HDECAY [26][27][28], focussing on the SM and MSSM, was the object of various extensions to the CP-conserving or -violating NMSSM for NMSSMTools [29][30][31][32] and NMSSMCALC [33,34]; SOFTSUSY [35,36] recently released its own set of routines [37], which generally are confined to the leading order or leading QCD corrections; a one-loop evaluation of two-body decays in the DR MS scheme for generic models [38] has been recently presented for SPHENO [39][40][41][42], which employs SARAH [43][44][45][46]. Moreover, the non-public code SloopS [47] has been extended to the NMSSM [48,49] and applied to the calculation of Higgs decays [49,50].
The current work focussing on NMSSM Higgs decays is part of the effort for developing a version of FeynHiggs [18,[51][52][53][54][55][56][57] dedicated to the NMSSM [25,58]. The general methodology relies on a Feynman-diagrammatic calculation of radiative corrections, which employs FeynArts [59,60], FormCalc [61] and LoopTools [61]. The implementation of the renormalization scheme within the NMSSM [25] has been done in such a way that the result in the MSSM limit of the NMSSM exactly coincides with the MSSM result obtained from FeynHiggs without any further adjustments of parameters (cases where the NMSSM result is more complete than the current implementation of the MSSM result will be discussed below). Concerning the Higgs decays, our routines, in their current status, contain an evaluation of all the two-body decays of neutral and charged Higgs bosons. In the present paper, we wish to focus on decays of the neutral Higgs bosons into SM final states, where we have obtained results including higher-order contributions as detailed below as well as further refinements. The channels of the type Higgs-to-Higgs and Higgs-to-SUSY are currently only implemented at leading order. 1 For the evaluation of the decays of the neutral Higgs bosons of the NMSSM into SM final states we have followed the same general approach as for the implementation of the MSSM Higgs decays in FeynHiggs, which have been described in Refs. [63,67]. At the level of the external Higgs fields, mixing effects are consistently taken into account as explained in Refs. [25,68]. Full one-loop contributions are considered in the fermionic decay channels, supplemented by QCD higher-order corrections. For the bosonic decay modes generated at the radiative order, leading QCD corrections are taken into account. For decays into massive electroweak gauge bosons, the implementation in the MSSM is such that FeynHiggs first extracts the loop-corrected width that Prophecy4f [69][70][71] calculates in the SM for a Higgs boson at a given mass, and then rescales this result by the squared coupling of the MSSM Higgs boson to W W and ZZ, normalized to the SM value. We go beyond this approach and include full one-loop on-shell results for these decay widths-however, the on-shell kinematical factor of the tree-level contribution is replaced by its off-shell counterpart, leading to a tree-level estimate below threshold. 2 More generally, the refinements that we implement, e. g. the inclusion of higher-order corrections, often surpass the assumptions made in public codes dedicated to the NMSSM. However, we stress that we strictly confine ourselves to the free-particle approximation for the final states. Accordingly, dedicated analyses would be needed for a proper treatment of decays close to threshold, since finite-width effects need to be taken into account in this region, and effects of the interactions among the final states can be very large [72][73][74][75][76].
The case where the Higgs spectrum contains a singlet-dominated Higgs state with a mass below the one of the signal detected at 125 GeV is a particularly interesting scenario that commonly emerges in the NMSSM. Among the appealing features of this class of scenarios it should be mentioned that the somewhat high mass (in an MSSM context) of the SM-like Higgs state observed at the LHC could be understood more naturally as the result of the mixing with a lighter singlet state-see e. g. Ref. [77] for a list of references and an estimate of the possible uplift in mass. Furthermore, the existence of a mostly singlet-like state in the range of 100 GeV has been suggested [78,79] as a possible explanation of the 2.3 σ local excess observed at LEP in the e + e − → Z(H → bb) searches [80]. More recently, the CMS collaboration has reported a local excess in its diphoton Higgs searches for a mass in the vicinity of ∼ 96 GeV [81]. This excess reaches the (local) 2.9 σ level in the Run II data and has already received attention from the particle-physics community [82][83][84][85]. A similar excess (2.0 σ) was already present in the Run I data in the same mass range. In an NMSSM context, the possibility of large diphoton signals is well-known [86][87][88]. Here we show that it is in fact possible to describe simultaneously the LEP and the CMS excesses within the NMSSM. However, it should be kept in mind that the excesses that were observed at LEP and CMS at 100 GeV could of course just be statistical fluctuations of the background and that their possible explanation in terms of an NMSSM Higgs state remains somewhat speculative. In particular, the diphoton excess observed at CMS would of course require confirmation from the ATLAS data as well [89].
In the following section, we discuss the technical aspects of our calculation, describing the conventions, assumptions and higher-order corrections that we address. Sect. 3 illustrates the workings of our decay routines in several scenarios of the NMSSM, and we perform comparisons with existing public tools. We also investigate the NMSSM scenario with a mostly singlet-like state with mass close to 100 GeV. Furthermore we discuss the possible size of the remaining theoretical uncertainties from unknown higher-order corrections. The conclusions are summarized in Sect. 4.

Higgs decays to SM particles in the CP-violating NMSSM
In this section, we describe the technical aspects of our calculation of the Higgs decays. Our notation and the renormalization scheme that we employ for the Z 3 -conserving NMSSM in the general case of complex parameters were presented in Sect. 2 of Ref. [25], and we refer the reader to this article for further details.

Decay amplitudes for a physical (on-shell) Higgs state -Generalities
On-shell external Higgs leg In this article, we consider the decays of a physical Higgs state, i. e. an eigenstate of the inverse propagator matrix for the Higgs fields evaluated at the corresponding pole eigenvalue. The connection between such a physical state and the tree-level Higgs fields entering the Feynman diagrams is non-trivial in general since the higher-order contributions induce mixing among the Higgs states and between the Higgs states and the gauge bosons (as well as the associated Goldstone bosons). The LSZ reduction fully determines the (non-unitary) transition matrix Z mix between the loop-corrected mass eigenstates and the lowest-order states. Then, the amplitude describing the decay of the physical state h phys i (we shall omit the superscript 'phys' later on), into e. g. a fermion pair ff , relates to the amplitudes in terms of the tree-level states h 0 j according to (see below for the mixing with gauge bosons and Goldstone bosons): Here, we characterize the physical Higgs states according to the procedure outlined in Ref. [25] (see also Refs. [53,63,68]): • the Higgs self-energies include full one-loop and leading O(α t α s , α 2 t ) two-loop corrections (with two-loop effects obtained in the MSSM approximation via the public code FeynHiggs 3 ); • the pole masses correspond to the zeroes of the determinant of the inverse-propagator matrix; • the (5 × 5) matrix Z mix is obtained in terms of the solutions of the eigenvector equation for the effective mass matrix evaluated at the poles, and satisfying the appropriate normalization conditions (see Sect. 2.6 of Ref. [25]).
In correcting the external Higgs legs by the full matrix Z mix -instead of employing a simple diagrammatic expansion-we resum contributions to the transition amplitudes that are formally of higher loop order. This resummation is convenient for taking into account numerically relevant leading higher-order contributions. It can in fact be crucial for the frequent case where radiative corrections mix states that are almost mass-degenerate in order to properly describe the resonance-type effects that are induced by the mixing. On the other hand, care needs to be taken to avoid the occurrence of non-decoupling terms when Higgs states are well-separated in mass, since higher-order effects can spoil the order-by-order cancellations with vertex corrections. We stress that all public tools, with the exception of FeynHiggs, neglect the full effect of the transition to the physical Higgs states encoded within Z mix , and instead employ the unitary approximation U 0 neglecting external momenta (which is in accordance with leading-order or QCDimproved leading-order predictions). We refer the reader to Refs. [25,53,68] for the details of the definition of U 0 or U m (another unitary approximation) as well as a discussion of their impact at the level of Higgs decay widths.
Higgs-electroweak mixing For the mass determination, we do not take into account contributions arising from the mixing of the Higgs fields with the neutral Goldstone or Z bosons since these corrections enter at the sub-dominant two-loop level (contributions of this kind can also be compensated by appropriate field-renormalization conditions [92]). We note that, in the CP-conserving case, only external CP-odd Higgs components are affected by such a mixing. Yet, at the level of the decay amplitudes, the Higgs mixing with the Goldstone and Z bosons enters already at the one-loop order (even if the corresponding self-energies are cancelled by an appropriate field-renormalization condition, this procedure would still provide a contribution to the h i ff counterterm). Therefore, for a complete one-loop result of the decay amplitudes it is in general necessary to incorporate Higgs-Goldstone and Higgs-Z self-energy transition diagrams [62,63,66]. In the following, we evaluate such contributions to the decay amplitudes in the usual diagrammatic fashion (as prescribed by the LSZ reduction) with the help of the FeynArts modelfile for the CP-violating NMSSM [25].
The corresponding one-loop amplitudes (including the associated counterterms) will be symbolically denoted as A 1L G/Z . These amplitudes can be written in terms of the self-energies Σ h i G/Z with Higgs and Goldstone/Z bosons in the external legs. In turn, these self-energies are connected by a Slavnov-Taylor identity (see e. g. App. A of Ref. [93]): 4 where the T h i correspond to the tadpole terms of the Higgs potential, and (U n ) ij are the elements of the transition matrix between the gauge-and tree-level mass-eigenstate bases of the Higgs bosonsthe notation is introduced in Sect. 2.1 of Ref. [25]. Similar relations in the MSSM are also provided in Eqs. (127) of Ref. [63]. We checked this identity at the numerical level.
Inclusion of one-loop contributions The wave function normalization factors contained in Z mix together with the described treatment of the mixing with the Goldstone and Z bosons ensure the correct on-shell properties of the external Higgs leg in the decay amplitude, so that no further diagrams correcting this external leg are needed. Moreover, the SM fermions and gauge bosons are also treated as on-shell particles in our renormalization scheme. Beyond the transition to the loopcorrected states incorporated by Z mix , we thus compute the decay amplitudes at the one-loop order as the sum of the tree-level contribution A tree (possibly equal to zero), the Higgs-electroweak oneloop mixing A 1L G/Z and the (renormalized) one-loop vertex corrections A 1L vert (including counterterm contributions)-we note that each of these pieces of the full amplitude is separately ultravioletfinite. In the example of the ff decay, the amplitudes with a tree-level external Higgs field h 0 j -on the right-hand side of Eq. (2.1)-thus symbolically read: All the pieces on the right-hand side of this equation are computed with the help of FeynArts [59,60], FormCalc [61] and LoopTools [61], according to the prescriptions that are encoded in the modelfile for the CP-violating NMSSM. However, we use a specific treatment for some of the contributions, such as QED and QCD one-loop corrections to Higgs decays into final state particles that are electrically and/or color charged, or include certain higher-order corrections. We describe these channel-specific modifications in the following subsections.
Goldstone-boson couplings The cubic Higgs-Goldstone-boson vertices can be expressed as The doublet vacuum expectation value (vev), v = M W s w / √ 2 π α, is expressed in terms of the gauge-boson masses M W and M Z s w = 1 − M 2 W /M 2 Z , as well as the electromagnetic coupling α. The symbol m 2 h j , (j = 1, . . . , 5), represents the tree-level mass squared of the neutral Higgs state h 0 j , and m 2 H ± the mass squared of the charged Higgs state. The use of the tree-level couplings of Eq. (2.4) together with a physical (loop-corrected) external Higgs leg h i = j Z mix ij h 0 j is potentially problematic regarding the gauge properties of the matrix elements. The structure of the gauge theory and its renormalization indeed guarantee that the gauge identities are observed at the order of the calculation (one loop). However, the evaluation of Feynman amplitudes is not protected against a violation of the gauge identities at the (incomplete) two-loop order. We detected such gauge-violating effects of two-loop order at several points in our calculation of the neutral-Higgs decays, e. g.: • the Ward identity in h i → γγ is not satisfied (see also Ref. [87]); • infrared (IR) divergences of the virtual corrections in h i → W + W − do not cancel their counterparts in the bremsstrahlung process h i → W + W − γ (see also Ref. [94]); • computing h i → ff in an R ξ gauge entails non-vanishing dependence of the amplitudes on the electroweak gauge-fixing parameters ξ Z and ξ W .
As these gauge-breaking effects could intervene with sizable and uncontrolled numerical impact, it is desirable to add two-loop order terms restoring the gauge identities at the level of the matrix elements. Technically, there are different possible procedures to achieve this: one would amount to replace the kinematic Higgs masses that appear in Higgs-gauge-boson couplings by tree-level Higgs masses; we prefer the alternative procedure consisting in changing the Higgs-Goldstone-boson couplings of Eq. (2.4): for the Higgs mass associated to the external Higgs leg the loop-corrected Higgs mass M h i is used instead of the tree-level one. This is actually the form of the Higgs-Goldstone-boson couplings that would be expected in an effective field theory of the physical Higgs boson h i . Using the definition of Z mix ij as an eigenvector of the loop-corrected mass matrix for the eigenvalue M 2 h i -see Sect. 2.6 of Ref. [25]-one can verify that the effective Higgs-Goldstoneboson vertices employing the physical Higgs mass differ from their tree-level counterparts by a term of one-loop order (proportional to the Higgs self-energies) so that the alteration of the one-loop amplitudes is indeed of two-loop order. Employing this shift of the Higgs-Goldstone couplings cures the gauge-related issues that we mentioned earlier.
Another issue with gauge invariance appears in connection with the amplitudes A 1L G/Z . The Goldstone and Z-boson propagators generate denominators with pole M 2 Z (or ξ Z M 2 Z in an R ξ gauge): in virtue of the Slavnov-Taylor identity of Eq. (2.2a) these terms should cancel one another in the total amplitude at the one-loop order-we refer the reader to Sect. 4.3 of Ref. [63] for a detailed discussion. However, the term ( h i (the loop-corrected Higgs mass), the cancellation is spoilt by a term of two-loop order. In order to address this problem, we re-define A 1L G/Z by adding a two-loop term: where Γ tree Gff represents the tree-level vertex of the neutral Goldstone boson with the fermion f (in the particular example of a Higgs decay into ff ). Then, it is straightforward to check thatÃ 1L G/Z is gauge-invariant. The transformation of Eq. (2.5) can also be interpreted as a two-loop shift redefining Σ h i Z , so that it satisfies a generalized Slavnov-Taylor identity of the form of Eq. (2.2a), but applying to a physical (loop-corrected) Higgs field, with the term Numerical input in the one-loop corrections As usual, the numerical values of the input parameters need to reflect the adopted renormalization scheme, and the input parameters corresponding to different schemes differ from each other by shifts of the appropriate loop order (at the loop level there exists some freedom to use a numerical value of an input parameter that differs from the tree-level value by a one-loop shift, since the difference induced in this way is of higher order). Concerning the input values of the relevant light quark masses, we follow in our evaluation the choice of FeynHiggs and employ MS quark masses with three-loop QCD corrections evaluated at the scale of the mass of the decaying Higgs, m MS q (M h i ), in the loop functions and the definition of the Yukawa couplings. In addition, the input value for the pole top mass is converted to m MS t (m t ) using up to two-loop QCD and one-loop top Yukawa/electroweak corrections (corresponding to the higher-order corrections included in the Higgs-boson mass calculation). Furthermore, the tan β-enhanced contributions are always included in the defining relation between the bottom Yukawa coupling and the bottom mass (and similarly for all other down-type quarks). Concerning the Higgs vev appearing in the relation between the Yukawa couplings and the fermion masses, we parametrize it in terms of α(M Z ). Finally, the strong coupling constant employed in SUSY-QCD diagrams is set to the scale of the supersymmetric particles entering the loop. We will comment on deviations from these settings if needed.

Higgs decays into SM fermions
Our calculation of the Higgs decay amplitudes into SM fermions closely follows the procedure outlined in the previous subsection. However, we include the QCD and QED corrections separately, making use of analytical formulae that are well-documented in the literature [95,96]. We also employ an effective description of the Higgs-bb interactions in order to resum potentially large effects for large values of tan β. Below, we comment on these two issues and discuss further the derivation of the decay widths for this class of channels.
Tree-level amplitude At the tree level, the decay h 0 j → ff is determined by the Yukawa coupling Y f and the decomposition of the tree-level state h 0 j in terms of the Higgs-doublet components: The δ-s are Kronecker symbols selecting the appropriate Higgs matrix element for the fermionic final state, u k = u, c, t, d k = d, s, b or e k = e, µ, τ . We have written the amplitude in the Diracfermion convention, separating the scalar piece g S h j f f (first two terms between curly brackets in the first line) from the pseudoscalar one g P h j f f (last two terms). The fermion and antifermion spinors are denoted asū f (p f ) and v f (pf ), respectively.
Case of the bb final state: tan β-enhanced corrections In the case of a decay to bb (and analogously for down-type quarks of first and second generation, but with smaller numerical impact), the loop contributions that receive a tan β enhancement may have a sizable impact, thus justifying an effective description of the Higgs-bb vertex that provides a resummation of large contributions [33,63,[97][98][99][100][101][102]. We denote the neutral components of H 1 and H 2 from Eq. (2.2) of Ref. [25] by H 0 d and H 0 u , respectively. The large tan β-enhanced effects arise from contributions to the (H 0 u ) * b P L b operator-P L,R are the left-and right-handed projectors in the Dirac description of the b spinors-and can be parametrized in the following fashion: Here, ∆ b is a coefficient that is determined via the calculation of the relevant (tan β-enhanced) one-loop diagrams to the Higgs-bb vertex, involving gluino-sbottom, chargino-stop and neutralinosbottom loops. 5 The symbol µ eff represents the effective µ term that is generated when the singlet field acquires a vev. The specific form of the operator, (S H 0 u ) * b P L b, is designed so as to preserve the Z 3 symmetry, and it can be shown that this operator is the one that gives rise to leading contributions to the tan β-enhanced effects. We evaluate ∆ b at a scale corresponding to the arithmetic mean of the masses of the contributing SUSY particles: this choice is consistent with the definition of ∆ b employed for the Higgs-mass calculation. From the parametrization of Eq. (2.7), one can derive the non-trivial relation between the 'genuine' Yukawa coupling Y b and the effective bottom mass m b : . Then, the effective couplings of the neutral Higgs fields to bb read: This can be used to substitute A tree h 0 j → bb in Eq. (2.3) by: where this expression resums the effect of tan β-enhanced corrections to the h 0 j bb vertex. However, if one now adds the one-loop amplitude A 1L vert , the one-loop effects associated with the tan β-enhanced contributions would be included twice. To avoid this double counting, the terms that are linear in ∆ b in Eq. (2.8) need to be subtracted. Employing the 'subtraction' couplings we define the following 'tree-level' amplitude for the Higgs decays into bottom quarks: QCD and QED corrections The inclusion of QCD and QED corrections requires a proper treatment of IR effects in the decay amplitudes. The IR-divergent parts of the virtual contributions by gluons or photons in A 1L vert are cancelled by their counterparts in processes with radiated photons or gluons. We employ directly the QCD and QED correction factors that are well-known analytically (see below) and therefore omit the Feynman diagrams involving a photon or gluon propagator when computing with FeynArts and FormCalc the one-loop corrections to the h 0 j ff vertex and to the fermion-mass and wave-function counterterms. The QCD-and QED-correction factors applying to the fermionic decays of a CP-even Higgs state were derived in Ref. [95]. The CP-odd case was addressed later in Ref. [96]. In the CP-violating case, it is useful to observe that the h j ff scalar and pseudoscalar operators do not interfere, so that the CP-even and CP-odd correction factors can be applied directly at the level of the amplitudes-although they were obtained at the level of the squared amplitudes: Here, Q f is the electric charge of the fermion f , C 2 (f ) is equal to 4/3 for quarks and equal to 0 for leptons, M h i corresponds to the kinematic (pole) mass in the Higgs decay under consideration and the functions ∆ S,P are explicated in e. g. Sect. 4 of Ref. [105].
4 . As noticed already in Ref. [95], the leading logarithm in the QCD-correction factor can be absorbed by the introduction of a running MS fermion mass in the definition of the Yukawa coupling Y f . Therefore, it is motivated to factorize m MS f (M h i ), with higher orders included in the definition of the QCD beta function.
The QCD (and QED) correction factors generally induce a sizable shift of the tree-level width of as much as ∼ 50%. While these effects were formally derived at the one-loop order, we apply them over the full amplitudes (without the QCD and QED corrections), i. e. we include the oneloop vertex amplitude without QCD/QED corrections A 1L wo. QCD/QED vert and A 1L G/Z in the definitions of the couplings g S,P h j f f that are employed in Eq. (2.12)-we will use the notation g S,P 1L h j f f below. The adopted factorization corresponds to a particular choice of the higher-order contributions beyond the ones that have been explicitly calculated.
Decay width Putting together the various pieces discussed before, we can express the decay amplitude at the one-loop order as: Summing over spinor and color degrees of freedom, the decay width is then obtained as: (2.14) At the considered order, we could dismiss the one-loop squared terms in A h i → ff 2 . However, in order to tackle the case where the contributions from irreducible one-loop diagrams are numerically larger than the tree-level amplitude, we keep the corresponding squared terms in the expression above (it should be noted that the QCD and QED corrections have been stripped off from the one-loop amplitude that gets squared). The approach of incorporating the squared terms should give a reliable result in a situation where the tree-level result is significantly suppressed, since the other missing contribution at this order consisting of the tree-level amplitude times the two-loop amplitude would be suppressed due to the small tree-level result. In such a case, however, the higher-order uncertainties are expected to be comparatively larger than in the case where one-loop effects are subdominant to the tree level.
The kinematic masses of the fermions are easily identified in the leptonic case. For decays into top quarks the 'pole' mass m t is used, while for all other decays into quarks we employ the MS masses evaluated at the scale of the Higgs mass m MS q (M h i ). We note that these kinematic masses have little impact on the decay widths, as long as the Higgs state is much heavier. In the NMSSM, however, singlet-like Higgs states can be very light, in which case the choice of an MS mass is problematic. Yet, in this case the Higgs state is typically near threshold so that the free-parton approximation in the final state is not expected to be reliable. Our current code is not properly equipped to address decays directly at threshold independently of the issue of running kinematic masses. Improved descriptions of the hadronic decays of Higgs states close to the bb threshold or in the chiral limit have been presented in e. g. Refs. [74-76, 106-108].

Decays into SM gauge bosons
Now we consider Higgs decays into the gauge bosons of the SM. Almost each of these channels requires a specific processing in order to include higher-order corrections consistently or to deal with off-shell effects.
Decays into electroweak gauge bosons Higgs decays into on-shell W -s and Z-s can be easily included at the one-loop order in comparable fashion to the fermionic decays. However, the notion of W W or ZZ final states usually includes contributions from off-shell gauge bosons as well, encompassing a wide range of four-fermion final states. Such off-shell effects mostly impact the decays of Higgs bosons with a mass below the W W or ZZ thresholds. Instead of a full processing of the off-shell decays at one-loop order, we pursue two distinct evaluations of the decay widths in these channels.
Our first approach is that already employed in FeynHiggs for the corresponding decays in the MSSM. It consists in exploiting the precise one-loop results of Prophecy4f for the SM-Higgs decays into four fermions [69][70][71]. For an (N)MSSM Higgs boson h i , the SM decay width is thus evaluated at the mass M h i and then rescaled by the squared ratio of the tree-level couplings to gauge bosons for h i and an SM Higgs boson H SM (V = W, Z): represents the decay width of the physical Higgs state h i in the NMSSM, denotes the decay width of an SM-Higgs boson with the mass M h i . The matrix elements R ij reflect the connection between the tree-level Higgs states and the physical states. This role is similar to Z mix . However, decoupling in the SM limit of the model yields the additional condition that the ratio in Eq. (2.15a) reduces to 1 in this limit for the SM-like Higgs boson of the NMSSM. For this reason, FeynHiggs employs the matrix U m (or U 0 ) as a unitary approximation of Z mix -see Sect. 2.6 of Ref. [25]. An alternative choice consists in using X ij ≡ Z mix ij k |Z mix ik | 2 . However, the difference of the widths when employing U 0 , U m , Z mix or X ≡ (X ij ) corresponds to effects of higher order, which should be regarded as part of the higher-order uncertainty. The rescaling of the one-loop SM width should only be applied for the SM-like Higgs of the NMSSM, where this implementation of the h i → V V widths is expected to provide an approximation that is relatively close to a full one-loop result incorporating all NMSSM contributions. However, for the other Higgs states of the NMSSM one-loop contributions beyond the SM may well be dominant. Actually, the farther the quantity [R ij · (U n ) j2 ] [R ij · (U n ) j1 ] departs from tan β, the more inaccurate the prediction based on SM-like radiative corrections becomes.
Our second approach consists in a one-loop calculation of the Higgs decay widths into on-shell gauge bosons (see Ref. [94] for the MSSM case), including tree-level off-shell effects. This evaluation is meant to address the case of heavy Higgs bosons at the full one-loop order. The restriction to on-shell kinematics is justified above the threshold for electroweak gauge-boson production (offshell effects at the one-loop level could be included via a numerical integration over the squared momenta of the gauge bosons in the final state-see Refs. [109,110] for a discussion in the MSSM). Our implementation largely follows the lines described in Sect. 2.1, with the noteworthy feature that contributions from Higgs-electroweak mixing A 1L G/Z vanish. In the case of the W + W − final state, the QED IR-divergences are regularized with a photon mass and cancel with bremsstrahlung corrections: soft and hard bremsstrahlung are included according to Refs. [111,112] (see also [94]). We stress that the exact cancellation of the IR-divergences is only achieved through the replacement of the h i G + G − coupling by the expression in terms of the kinematical Higgs mass, as discussed in Sect. 2.1. This fact had already been observed by Ref. [94]. In order to extend the validity of the calculation below threshold, we process the Born-order term separately, applying an off-shell kinematic integration over the squared external momentum of the gauge bosons-see e. g. Eq. (37) in Ref. [113]. Thus, this evaluation is performed at tree level below threshold and at full one-loop order (for the on-shell case) above threshold. The vanishing on-shell kinematical factor multiplying the contributions of one-loop order ensures the continuity of the prediction at threshold. Finally, we include the one-loop squared term in the calculation. Indeed, as we will discuss later on, the treelevel contribution vanishes for a decoupling doublet, meaning that the Higgs decays to W W/ZZ can be dominated by one-loop effects. To this end, the infrared divergences of two-loop order are regularized in an ad-hoc fashion-which appears compulsory as long as the two-loop order is incomplete-making use of the one-loop real radiation and estimating the logarithmic term in the imaginary part of the one-loop amplitude.
Radiative decays into gauge bosons Higgs decays into photon pairs, gluon pairs or γZ appear at the one-loop level-i. e. A tree = 0 for all these channels. We compute the one-loop order using the FeynArts modelfile, although the results are well-known analytically in the literature-see e. g. Ref. [87] or Sect. III of Ref. [50] ( [113] for the MSSM). The electromagnetic coupling in these channels is set to the value α(0) corresponding to the Thomson limit. The use of tree-level Higgs-Goldstone couplings together with loop-corrected kinematic Higgs masses M h i in our calculation would induce an effective violation of Ward identities by two-loop order terms in the amplitude: as explained in Sect. 2.1, we choose to restore the proper gauge structure by re-defining the Higgs-Goldstone couplings in terms of the kinematic Higgs mass M h i . Since our calculation is restricted to the leading-here, one-loop-order, the transition of the amplitude from tree-level to physical Higgs states is performed via U m or X instead of Z mix in order to ensure the appropriate behavior in the decoupling limit.
Leading QCD corrections to the diphoton Higgs decays have received substantial attention in the literature. A frequently used approximation for this channel consists in multiplying the amplitudes driven by quark and squark loops by the factors respectively-see e. g. Ref. [114]. However, these simple factors are only valid in the limit of heavy quarks and squarks (compared to the mass of the decaying Higgs boson). More general analytical expressions can be found in e. g. Ref. [115]. In our calculation, we apply the correction factors 1 + C S (τ q ) α s (M h i )/π and 1 + C P (τ q ) α s (M h i )/π to the contributions of the quark q to the CP-even and the CP-odd h i γγ operators, respectively, and 1 + C(τQ) α s (M h i )/π to the contributions of the squarkQ (to the CP-even operator). Here, τ X denotes the ratio 4 m 2 The coefficients C S,P and C are extracted from Ref. [116] and Ref. [117]. In order to obtain a consistent inclusion of the O(α s ) corrections, the quark and squark masses m X entering the one-loop amplitudes or the correction factors are chosen as defined in Eq. (5) of Ref. [116] and in Eq. (12) of Ref. [117] (rather than MS running masses).
The QCD corrections to the digluon decays include virtual corrections but also gluon and lightquark radiation. They are thus technically defined at the level of the squared amplitudes. In the limit of heavy quarks and squarks, the corrections are known beyond NLO-see the discussion in Ref. [113] for a list of references. The full dependence in mass was derived at NLO in Refs. [116,117], for both quark and squark loops. In our implementation, we follow the prescriptions of Eqs. (51), (63) and (67) of Ref. [113] in the limit of light radiated quarks and heavy particles in the loop. For consistency, the masses of the particles in the one-loop amplitude are taken as pole masses. Effects beyond this approximation can be sizable, as evidenced by Fig. 20 of Ref. [116] and Fig. 12 of Ref. [117]. As the CP-even and CP-odd Higgs-gg operators do not interfere, it is straightforward to include both correction factors in the CP-violating case. Finally, we note that parts of the leading QCD corrections to h i → gg are induced by the real radiation of quarkantiquark pairs. In the case of the heavier quark flavors (top, bottom and possibly charm), the channels are experimentally well-distinguishable from gluonic decays. Therefore, the partial widths related to these corrections could be attached to the Higgs decays into quarks instead [118]. The resolution of this ambiguity would involve a dedicated experimental analysis of the kinematics of the gluon radiation in h i → gqq (collinear or back-to-back emission). In the following section, we choose to present our results for the h i → gg decay in the three-radiated-flavor approach, while the contributions from the heavier quark flavors are distributed among the h i → cc, bb, and tt widths (provided the Higgs state is above threshold). These contributions to the fermionic Higgs decays are of O(α 2 s ). The QCD corrections to the quark loops of an SM-Higgs decay into γZ have been studied in Refs. [119][120][121], but we do not consider them here.

Numerical Analysis
In this section, we present our results for the decay widths of the neutral Higgs bosons into SM particles in several scenarios and compare them with the predictions of existing codes. While a detailed estimate of the uncertainty associated to missing higher-order corrections goes beyond the scope of our analysis, we will provide some discussion at the end of this section, based on our observations and comparisons.
Throughout this section, the top-quark pole mass is chosen as m t = 173.2 GeV. Moreover, all DR parameters are defined at the scale m t , and all stop parameters are treated as on-shell parameters. Concerning the Higgs phenomenology, we test the scenarios presented in this section with the full set of experimental constraints and signals implemented in the public tools HiggsBounds-4.3.1 (and 5.1.1beta) [122][123][124][125][126][127] and HiggsSignals-1.3.1 (and 2.1.0beta) [127][128][129]. We refer the reader to the corresponding publications for a detailed list of experimental references. The input parameters employed in our scenarios are summarized in Tab. 1.

Comparison with FeynHiggs in the MSSM-limit
The MSSM limit of the NMSSM is obtained at vanishingly small values of λ and κ: the singlet superfield then decouples from the MSSM sector but the µ eff term remains relevant as long as κ ∼ λ. It is then possible to compare our results for the Higgs decays to the corresponding predictions of FeynHiggs-2.13.0. The settings of FeynHiggs are thus adjusted in order to match the level of higher-order contributions and renormalization conditions of our NMSSM mass calculation: the corresponding FeynHiggs input flags read FHSetFlags[4,0,0,3,0,2,0,0,1,1]. We will denote the MSSM(-like) Higgs bosons as h and H, for the CP-even states, and A for the CP-odd one.
First, we consider the Higgs decays into SM fermions. We turn to a region of the parameter space of the CP-conserving NMSSM characterized by the input provided in the column ' Fig. 1' of Tab. 1, where we vary the masses and trilinear couplings associated with the squarks of the third generation.  #1 Sect   In this setup, the lightest Higgs state is SM-like, with a mass in the range [124, 126.5] GeV, while the heavy doublet states both have masses of about 997 GeV in this scenario. The full range under study is found to be in agreement with constraints in the Higgs sector, as implemented in HiggsBounds and HiggsSignals.
In Fig. 1, we show the variations of the Higgs decay widths into bb and tt for the doublet states (when kinematically allowed). The solid lines correspond to our predictions in several approximations: at the 'tree level' with Yukawa couplings defined in terms of running MS quark masses at the scale of the physical Higgs mass (black); including QCD and QED corrections (but without SQCD contributions) as well as the transition to the physical Higgs state via Z mix (blue); replacing also the Higgs couplings to the down-type quarks by their effective form as expressed in Eq. (2.9) (green); at full one loop, with higher-order improvements as described above (red). The green dotted line is similar to the solid green line, up to the replacement of Z mix by its unitary approximation U m . The purple diamonds are obtained with FeynHiggs. As the Higgs masses vary little over the range of the scan, the modification of the decay widths is essentially driven by radiative effects: this explains the relatively flat behavior of the tree-level results-at least in the case of the heavy states; for the light state, the mass and decay widths vary somewhat more. The same remark applies to the QCD/QED-corrected widths with Z mix , although the H → tt decay width displays a more pronounced variation due to the h 0 -H 0 mixing at the loop level. The one- where U m has been used instead of Z mix ; the red curve corresponds to our prediction at the full one-loop level (with higher-order improvements). The contribution from h i → g(g * → qq) is included within the results in the solid red curves. The purple diamonds mark the corresponding evaluation by FeynHiggs. We also consider the widths without the contribution from gg * , as well as using CP-even QCD/QED-correction factors in the case of the A decays (in order to match FeynHiggs) as a dashed red line (which is not visible in the case H → tt).
loop corrections to the bb decay width of the SM-like state shift this quantity upwards by ∼ 20%. The bulk of this effect, beyond the QCD corrections included in the running b quark mass, is already contained within the QCD/QED-corrected width. For the heavy-doublet states, however, QCD and QED corrections (beyond the effect encoded within the running Yukawa coupling) do not lead to a significant improvement of the prediction of the decay widths, as the 'tree-level' result often appears closer to the full one-loop widths than the blue curve-for the bb final state, this deviation is only partially explained by the radiative corrections that can be resummed within the effective Higgs couplings to the bottom quark (solid green curve).
In the case of the tt final state, the solid green curve and the blue curve are essentially identical as the effective couplings to the down-type quarks only play a secondary part. The difference between the solid and dotted green lines originates from the treatment regarding the transition matrix employed for the description of the external Higgs leg-Z mix or the approximation U m : the consequences for the predicted width reach O(5%). We note that, with the exception of FeynHiggs, public tools usually neglect the effects associated with the momentum dependence of the Higgs self-energies (only properly encoded within Z mix ). On the other hand, the estimate using U m is somewhat closer to the full one-loop result (using Z mix ), which means that, as long as one restricts to an 'improved tree-level approximation', the choice of a unitary transition matrix might provide slightly more reliable results than the same level of approximation with Z mix . As we mentioned earlier, the ∼ 10% difference between the red and green curves in the case of the heavy doublets hints at sizable electroweak effects of one-loop order. This can be understood in terms of large electroweak Sudakov logarithms for the heavy Higgs bosons. The impact of the sfermion spectrum on the decay widths into SM fermions consists in a suppression in the presence of light stops and sbottoms. This effect is of order 10% over the considered range of squark masses. Concerning the comparison with FeynHiggs, we observe a very good agreement of the predicted one-loop widths for the CP-even states: this is expected since we essentially apply the same processing of the parameters. However, a small discrepancy in the h → bb width is noticeable: it is related to our inclusion of the contribution from h → g(g * → bb) to the decay width. Subtracting this contribution (dashed red curve), we recover the FeynHiggs prediction. This effect is negligible for the other Higgs states. 6 Furthermore, for the CP-odd state, we find a deviation of a few percent in the tt channel: this difference is due to FeynHiggs employing CP-even QCD/QED-correction factors instead of the CP-odd ones. The agreement is restored if we adopt the same approximation (dashed red curve). No such discrepancy appears for the bb final state, as the CP-even and CPodd QCD/QED-correction factors converge in the limit of very light fermions (as compared to the mass of the Higgs state). Furthermore, the processing of the effective Higgs couplings to downtype quarks differs between FeynHiggs and our implementation, leading to small numerical effects (below 1%) for light SUSY spectra: we evaluate ∆ b at the scale defined by the arithmetic mean of the SUSY masses involved, while FeynHiggs employs a geometric mean ('modified Y eff b '). Another (numerically minor) difference with FeynHiggs is the slightly different treatment of the Goldstone-Higgs couplings regarding the restoration of gauge invariance of the result-see Sect. 2.1.
Then, we consider the Higgs decays into electroweak gauge bosons in the MSSM-limit. To this end, we perform a scan over the charged-Higgs mass in the range [150,2000] GeV; the rest of the parameters is set as in Fig. 1, but the stop and sbottom soft masses and trilinear couplings are frozen to 1.5 TeV and 2.3 TeV, respectively. Correspondingly, the lightest CP-even Higgs is SMlike with a mass of 124-125 GeV as soon as m H ± 350 GeV. The mass of the heavy doublet states varies over the scan and is comparable to m H ± . For clarity, the masses are shown in Fig. 2. Except for m H ± 300 GeV, where the SM-like Higgs is not in the desired experimental window, this scenario is consistent with experimental bounds implemented in HiggsBounds and HiggsSignals. In Fig. 3, we show our results for the decay widths of the neutral Higgs states into W W and ZZ. The black curves correspond to the prediction obtained by rescaling the SM width of Prophecy4f by the tree-level Higgs-W W/ZZ couplings (relative to the SM). The rotation to the loop-corrected Higgs state is then described in the unitary approximation U m . This is the approach employed by FeynHiggs (purple diamonds), leading to a good agreement. In the case of the SM-like state (upper plots), the decays are off-shell and both the (superposed) blue and red curves correspond to the tree-level results with off-shell kinematics but with Z mix applied to the external Higgs leg. Thus, the difference between the red and black lines must be interpreted as the magnitude of the SM radiative corrections implemented in Prophecy4f (which are not included in our result displayed by the red curve) and amounts to somewhat less than 10%. For the heavy CP-even Higgs (plots in the middle), the red and blue curves are differentiated by the inclusion of onshell one-loop contributions to the widths (red curves). While the results essentially agree with the 'Prophecy4f approach' at low mass (where the tree-level Higgs-W W/ZZ coupling remains sizable due to a substantial mixing of the heavy CP-even Higgs with the SM-like state), large deviations are observed for m H ± 500 GeV. Indeed, for this state the tree-level Higgs coupling to electroweak gauge bosons is very small (as a consequence of the decoupling limit), so that the decay width is largely dominated by one-loop effects arising both from the contributions to the external Higgs leg (via Z mix ) and in the vertex corrections. The dip of the 'Prophecy4f prediction' (black line) at m H ± 1.1 TeV is due to an exactly vanishing Higgs-gauge coupling (tree level, unitarily rotated) at this point in parameter space. This does not happen when the couplings are transformed to the physical Higgs state via Z mix , due to the imaginary part of this transition matrix. On the other hand, the 'tree-level · Z mix ' approach (blue line) does not offer an accurate estimate of the Higgs → W W/ZZ decay widths either, due to sizable vertex corrections. The 'sudden drop' of the H decay widths for low values of m H ± (for all the curves) is associated to the off-shell regime (the Higgs state has a mass below threshold). Finally, the CP-odd Higgs decays into W W and ZZ (lower plots) are generated at the radiative level. In this case, all the approximations based on tree-level predictions ('Prophecy4f', 'Z mix ', FeynHiggs) vanish, and therefore only our one-loop on-shell description (red curve) is displayed in the plots. The peculiar shape of the predicted decay widths-in particular the peak at m H ± 1 TeV-is associated to various thresholds-in particular the two-wino threshold at 1 TeV. We stress that one-loop effects dominate the decays of the heavy doublet Higgs states into W W/ZZ, which means that our results, albeit formally of next-to-leading order, come with a large uncertainty due to QCD (two-loop) corrections: in particular, we observed that the use of pole instead of MS quark masses within the loop could shift the widths by ∼ 50%.
We choose to discuss the diphoton and digluon decay widths in the MSSM-limit in a scenario where tan β scans the range [1,40]-m H ± is set to 1 TeV again; the details of the input are available in Tab. 1. The mass of the SM-like Higgs state is of order 100 GeV at tan β = 1, but settles in the interval [123.5, 126.5] when tan β 7. Correspondingly, the low-tan β limit is disfavored by HiggsSignals in this scenario. The heavy doublet states have a mass of about 996-997 GeV. HiggsBounds excludes the large-tan β ∼ 40 endpoints, due to constraints on heavy-Higgs searches in the τ τ channel [130].
We display the Higgs decay widths into gg and γγ in Fig. 4. The transition to the physical Higgs states is performed using various approximations: U 0 (black curves), U m (blue curves) and X (red curves). These three descriptions agree rather well. The predictions for the digluon final state employ the QCD corrections for three radiated quark flavors-as radiated cc, bb or tt are regarded as contributions to the fermionic Higgs decays. Nevertheless, in order to compare with FeynHiggs, we also show the results in the five-radiated-flavor approach (solid green curves) and cutting off universal QCD corrections beyond NLO (dashed green curves). The resulting deviation from the predictions of FeynHiggs (purple diamonds) is resolved when replacing the pole quark masses by MS masses in the one-loop amplitude (dotted green curves). The appropriate choice for  the considered QCD-correction factor is that of pole masses in the loop, however. The difference between the widths depicted by the black/blue/red lines and by the solid green lines is of order 20%: this enhancement is due to the larger number of radiated flavors. Universal QCD corrections beyond NLO (difference between the solid and the dashed green curves) represent almost 15% of the width. In the case of the diphoton widths, the results of FeynHiggs (purple diamonds) should be compared to our predictions employing U m (blue curves): a deviation of somewhat less than 10% is noticeable for the heavy states. We checked that this discrepancy can be interpreted in terms of the heavy-quark/squark approximation that FeynHiggs employs for the NLO QCD corrections as well as the use of MS running masses (instead of the running masses defined in Eq. (5) of Ref. [116]): simplifying our processing of the widths to this approximation (dashed blue curves) yields a very good agreement with the results of FeynHiggs.
Finally, we consider a CP-violating scenario in Fig. 5: the parameters are set as indicated in the column ' Fig. 5' of Tab. 1. We perform a scan over φ At ∈ [−π, π]. The doublet Higgs states have masses of about ∼ 125 GeV, 493 GeV and 494 GeV. This scenario is phenomenologically consistent with the limits on the Higgs sector as implemented in HiggsBounds and HiggsSignals. Ideally, limits from the measured electric dipole moments (EDMs) should be considered as well: in particular, Barr-Zee contributions involving squark loops are sensitive to variations of φ At . On the other hand, such effects are relatively suppressed given the high mass of the stops (∼ 1.5 TeV). In any case, we mainly consider this scenario for the sake of comparison. Fig. 5 displays the Higgs decay widths into bb and W W . For the bb final state, we consider the tree-level widths (corresponding to an MS running Yukawa coupling; black lines), incorporate QCD/QED corrections and transform to the physical Higgs state using Z mix (blue), subtract and redistribute the SUSY corrections to the bottom-quark mass in the definition of the Higgsbottom couplings (green) and finally evaluate the widths at full one-loop order (red, including two-loop pieces as described in Sect. 2). The predictions of FeynHiggs (purple diamonds) are in good agreement with our full one-loop results. While the inclusion of the corrections associated to QED/QCD effects and the definition of effective Higgs couplings to bb notably improve the tree-level width of the SM-like state, as compared to the one-loop results-from a discrepancy of ∼ 20% to less than ∼ 4%-the performance of these 'leading' corrections is less convincing in the case of the heavy doublet states: the deviations with respect to the full one-loop results remain of order 5-10%.
Turning to the W W final state, we expectedly recover the results of FeynHiggs in the approximation with the rescaled SM widths of Prophecy4f (black lines). The impact of the one-loop corrections (red curves) on the width of the mostly CP-even heavy doublet state h 2 are rather mild, which we should put in perspective with the fact that this state is comparatively light. However, for the mostly CP-odd state h 3 , the tree-level approximations sizably underestimate the one-loop widths.
To summarize, in this comparison with FeynHiggs in the MSSM limit of the NMSSM, we were able to quantitatively recover the widths predicted by FeynHiggs and interpret the origin of the differences with our results. In the case of the Higgs decays into electroweak gauge bosons, our one-loop approach goes beyond the current approach of FeynHiggs and shows that the tree-level approximation for the SUSY contributions-even though rescaled from the Prophecy4f SM widthsleads to significant deviations for the heavy doublet states. In the case of the digluon decay, universal QCD corrections beyond NLO as well as the number of radiated quark flavors have a sizable impact on the widths. Finally, we observed that accounting for the mass dependence in the QCD corrections to the diphoton widths has a mild effect on the decay of the heavy states. It is planned to include all the refinements that go beyond the current status of FeynHiggs and that we have employed here into the predictions of the MSSM Higgs decays of a future update of FeynHiggs.

Comparison with NMSSMCALC
We now depart from the MSSM-limit of the NMSSM. We first investigate the Higgs decays in scenarios that are similar to those that we considered in the MSSM-limit. We will compare our estimates for the Higgs decay widths to the predictions of NMSSMCALC-2.1 [33,34]. In order to minimize the impact of the Higgs-mass calculation, we bypass the mass-evaluation routines of NMSSMCALC-we refer the reader to Refs. [20,25] for a comparison of the Higgs-mass predictions-and directly inject our Higgs spectrum (two-loop masses and U m mixing matrix) in the SLHA [131,132] file serving as an interface between the mass-evaluation and the decayevaluation routines. We proceed similarly with the squark masses and mixing angles. The subroutine of NMSSMCALC that evaluates the Higgs decays is based on a generalization of HDECAY_6.1 [26,27]. The corresponding widths include leading NLO effects (e. g. QCD/QED corrections, effective Higgs couplings to the bottom quark, etc.) but do not represent a complete one-loop order evaluation. Furthermore, internal parameters, renormalization scales or RGE runnings are not strictly identical to our choices. For instance, the decay widths predicted by NMSSMCALC are systematically normalized to G F , while we employ other parametrizations in terms of M W , M Z and e. g. α(M Z ): the corresponding tree-level contributions numerically differ by a few percents. In the case of NMSSMCALC, such effects are of one-loop electroweak order, hence beyond the considered approximation. In our calculation, the one-loop electroweak order is consistently implemented with respect to our renormalization scheme. Other small numerical differences appear at the level of e. g. running quark masses. Therefore, the numerical comparison between the two sets of results is expected to show a certain level of deviations.
We set κ = 0.4, λ = 0.3 and A κ = −100 GeV. In Fig. 6, we consider the Higgs decay widths into bb and tt in a scenario where the squark masses of third generation vary between 0.7 and 2 TeV while the trilinear couplings are in the interval [1.4, 4] TeV (see Tab. 1 for details). Correspondingly, the lightest CP-even Higgs state is SM-like with a mass of ∼ 120-126 GeV (the lowest range in m Q is in tension with the measured Higgs data); the second-lightest CP-even and the lightest CP-odd Higgs states are singlet-like, with masses of ∼ 643 GeV and ∼ 319 GeV respectively; the heavy doublet states have masses of ∼ 999 GeV. This scenario satisfies the tests of HiggsBounds and HiggsSignals in most of the mQ range. The decay widths are shown at the tree level (with MS Yukawa couplings; black curves), in the approximation of QCD/QED corrections and including the transition by Z mix (blue curves), including the leading SUSY corrections to the relation between the bottom-quark mass and the bottom Yukawa coupling, and furthermore substituting Z mix by U m (green curves), and finally for our full one-loop calculation (red curves). These various approximations perform in the same fashion as what we observed in the MSSM-limit: while the inclusion of QCD/QED-and Yukawa-driven leading corrections improve the tree-level-based predictions for the width of the SM-like state h 1 , the improvement is less obvious for the other Higgs states, pointing at sizable electroweak effects. The purple diamonds represent the predictions of NMSSMCALC/HDECAY: they should correspond to the approximation of our results shown in green. These results qualitatively agree at the level of a few percent (as expected, given e. g. the differing parametrizations at the tree level). These 'improved tree-level' predictions of the decay widths typically remain 5-10% away from our full one-loop implementation.
Then, we study the Higgs decays to electroweak gauge bosons as a function of the charged-Higgs mass-the squark masses and trilinear couplings of third generation are frozen to 1.5 TeV and 3 TeV respectively. The lightest CP-even Higgs with mass ∼ 122-127 GeV is SM-like. The mostly singletlike CP-even and CP-odd states have masses of order 650 GeV and 300 GeV, respectively. With growing m H ± , the masses of the heavy doublet states increase in proportion and successively cross the singlet masses, leading to strong mixing regimes. For clarity, the Higgs masses are plotted in Fig. 7. The decay widths into W W and ZZ are displayed in Fig. 8   where the SM widths of Prophecy4f are rescaled (black curves), in the tree-level approximation using Z mix (blue curves), in our on-shell one-loop description (with off-shell tree-level contributions; red curves). The predictions of NMSSMCALC are in agreement with our result obtained by a rescaling of the Prophecy4f widths. Again, we find that this approximation is not appropriate for heavy doublet Higgs states (h 3 for m H ± 650 GeV), in which case one-loop contributions dominate the decays into W W and ZZ. Differences may reach a factor 100 in certain mass ranges/scenarios. On the other hand, the various predictions for the width of the CP-even singlet at ∼ 650 GeV (h 3 for low m H ± , then h 2 ) are all consistent with one another: at least in the scenario under consideration, the tree-level contribution remains dominant for this state. At m H ± ∼ 400 GeV, we observe a suppression of the decay widths of the heavy doublet state h 2 : once again, it is associated with a vanishing h 2 -W W/ZZ coupling. One-loop effects do not lead to large deviations for this state. Furthermore, the 'drop' at low m H ± in the h 2 → ZZ channel can be traced back to the crossing of the kinematical threshold for on-shell Z bosons (the width below threshold is suppressed). Finally, the decays of the CP-odd Higgs are generated at the one-loop level: consequently, they are nontrivial only in our one-loop approach.
We investigate the Higgs decay widths into gluon and photon pairs as a function of tan β. The results are displayed in Fig. 9. The Higgs masses are of order ∼ 125 GeV for the SM-like state (except for tan β 3, in which case the corresponding mass is too low to satisfy the observed Higgs properties), ∼ 1 TeV for the heavy-doublet states, ∼ 320 GeV and ∼ 645 GeV for the CP-odd and CP-even mostly singlet-like states, respectively. The points with tan β 35 appear to be in tension with LHC searches for heavy Higgs bosons in the τ + τ − decay channel. Also in this case, we consider several descriptions of the transition to the physical Higgs state (black, blue and red lines): all the predictions agree to a good accuracy. In the case of the digluon decays, we also show the widths obtained with five radiated quarks (green curves): this result is essentially in agreement with the prediction of NMSSMCALC (purple diamonds), using the same approach. For the diphoton decays, NMSSMCALC also applies full NLO QCD corrections. Again, we observe a good agreement with our results. We checked that the remaining discrepancies-at the percent level-between the predictions from NMSSMCALC for the digluon and diphoton widths and ours are largely accounted for by a relative normalization factor G F M 2 W s 2 w √ 2 (π α) = [1 + O(α)] and minor deviations in the running of α s and the quark masses.
Finally, we consider a CP-violating scenario where we set λ = 0.2, |κ| = 0.6, tan β = 25 GeV and A κ = −750 GeV. We scan over φ κ ∈ [−π, π], which is phenomenologically realistic in the sense that EDMs in principle allow for large variations of this phase [133,134]; however, scenarios with large CP-violating mixing of the Higgs states-as the one we consider-tend to be constrained. The Higgs spectrum consists of an SM-like state with a mass of ∼ 124.5 GeV, a triplet of CP-even/CP-odd doublet and CP-even singlet states near ∼ 1 TeV with large and fluctuating mixing depending on φ κ , as well as a mostly CP-odd singlet at ∼ 1.2 TeV. The masses are depicted   Figure 9: Comparison with NMSSMCALC of the decay widths of the Higgs states into gg and γγ. The input parameters are provided in Tab. 1. The black, blue and red curves correspond to a transition to the physical Higgs state via U 0 , U m and X respectively. In the case of the gg final state, the black/blue/red lines are obtained for three radiated quark flavors in the QCD corrections. The green lines are obtained for five radiated quark flavors, which agrees with the approach of NMSSMCALC (purple diamonds).
in Fig. 10. The Higgs decays into bb and W W are plotted in Fig. 11. For the bb final state we find a relatively good agreement between the predictions of NMSSMCALC (purple diamonds) and our predictions employing the same approximation (shown in green: transition via U m , QCD/QED corrections and effective SUSY-corrected Higgs-bb couplings). A sizable discrepancy with our full one-loop result appears for |φ κ | 2.6 in the case of h 2,3 : this difference originates in large Higgsmixing effects encoded within Z mix that are not captured by the U m approximation (see Sect. 3.3 and Sect. 3.4 of Ref. [25] for a detailed discussion). In the case of the W W channel, deviations can be large when the tree-level approximation fails to capture the leading contribution to the decay width (as e. g. for h 2,3,4 ): the NMSSMCALC predictions are comparable to our rescaled widths from Prophecy4f.
We have thus observed a qualitative agreement of our results with the decay widths predicted by NMSSMCALC, when employing the same approximations as this code. Our calculation goes beyond the approach of NMSSMCALC/HDECAY at the level of the Higgs decays into SM fermions or into W W/ZZ, since we consider the full one-loop order. Sizable differences may thus appear, e. g. in the case of the decays into electroweak gauge bosons. Concerning the decays into gluon or photon pairs, the two codes are considering the decay widths at the same order, and the results are comparable.

A singlet dominated state at 100 GeV and possible explanation of slight excesses in the CMS and LEP data
The possible presence of a singlet-dominated state with mass in the ballpark of ∼ 100 GeV is a long-standing phenomenological trademark of the NMSSM-see e. g. Refs. [78,79]. One motivation for such a scenario is for instance the 2.3 σ local excess observed in Higgs searches at LEP in the e + e − → Z(H → bb) channel [80], which would be consistent with a scalar mass of ∼ 98 GeV (but with a rather coarse mass resolution). It would correspond to a signal strength with respect to the SM at the level of ∼ 10%. A natural candidate to explain this excess consists in a mostly singlet-like Higgs with a doublet component of about 10% (mixing squared). Interestingly, recent LHC Run II results [81] for CMS Higgs searches in the diphoton final state show a local excess of ∼ 3 σ in the vicinity of ∼ 96 GeV, while a similar upward fluctuation of 2 σ had been observed in the CMS Run I data at a comparable mass. A hypothetical signal of this kind would amount to 60% ± 20% of that of an SM Higgs boson at the same mass. In the NMSSM, relatively large Higgs branching fractions into γγ are possible due to the three-state mixing, in particular when the effective Higgs coupling to bb becomes small-see e. g. Refs. [86,87]. However, it then appears more difficult to interpret the LEP excess simultaneously. 7 Below, we consider the decays of a light 7 For instance, the authors of Ref. [135] came to a negative answer when considering the Run I 'excess' together with Dark Matter constraints in a specific region of the parameter space. singlet-like Higgs in this regime employing our calculation, to show that it is indeed possible to describe both 'excesses' simultaneously (without exploiting all possibilities within the NMSSM to desribe these effects).
We first consider a region of the NMSSM parameter space characterized by the input displayed in the column ' Fig. 13' of Tab. 1: we scan over A κ with tan β = 12. In this regime, the lightest CP-even state is singlet-like (except at the upper boundary in A κ ). Its mass varies within the range [50,125] GeV. The second-lightest Higgs has SM-like properties with a mass of ∼ 125 GeV. The variations of these masses as functions of A κ are shown in Fig. 12. We do not discuss the CP-odd singlet at ∼ 700 GeV and the heavy doublet states at ∼ 1.4 TeV in this context. The experimental constraints from Higgs searches, as summarized in HiggsBounds and HiggsSignals, are satisfied over the whole interval. In Fig. 13, we display the branching ratios of the lightest CP-even Higgs state into bb, ZZ, gg and γγ-the total width is calculated as the sum of the decay widths in the bb, τ + τ − , cc, W W , ZZ, gg and γγ channels. The tree-level results (including a tree-level evaluation of the full width, e. g. excluding the gg or γγ channels) are shown in black while the full one-loop results appear in red. Various approximations of the one-loop branching ratios are shown in blue and green: they are in good agreement with the full one-loop prediction for all the considered channels. We observe sizable deviations for the tree-level and one-loop predictions of the branching fraction into bb: this is due to the existence of a region with suppressed h 1 -bb coupling arising as a consequence of the mixing between the three neutral CP-even Higgs states. This suppressed region occurs below the displayed range of M h 1 at the tree level, but is shifted to M h 1 ∼ 90-100 GeV at the one-loop order. At its minimum for M h 1 95 GeV, the h 1 → bb decay is in fact dominated by its O(α 2 s ) contribution from h 1 → g(g * → bb). Concerning the ZZ final state, the difference between tree-level and one-loop branching fractions is largely driven by the magnitude of the full width: the large h 1 → gg contributions are absent at the tree level, and the large h 1 → bb contributions differ sizably. The branching fractions into gg and γγ exist only at the loop level. They show a peak at M h 1 ∼ 90-100 GeV, which mostly originates from the suppression of the h 1 → bb decay in this range. In particular, the gg final state contributes about 40% to the total width. Finally, we display the quantities ξ b and ξ γ , defined as follows: . (3.1b) These definitions of ξ b,γ give loose estimates of the signals that h 1 would generate in the LEP searches for e + e − → Z(H → bb) and the LHC searches for pp → H → γγ, normalized to the SM crosssections. The decay widths of an SM-Higgs boson are evaluated with our code, taking the SM-limit of the NMSSM. In our example, the signal in the bb searches would thus reach a few percent of an SM signal for M h 1 ∼ 100-110 GeV, but would be very suppressed in the range M h 1 ∼ 90-100 GeV, due to the small branching fraction into bb. On the other hand, the signal in the diphoton channel would amount to ∼ 20% of the corresponding SM cross-section for M h 1 ∼ 90-100 GeV. The Higgs properties for a specific point in this scenario are provided in the column ' Fig. 13' of Tab. 2. The production cross-sections, approximated by Γ ZZ for LEP and Γ gg for the LHC (in gluon-gluon fusion), amount to a comparable fraction-a few percent-of the corresponding SM-Higgs cross-sections. The branching ratio into bb is considerably reduced, as compared to an SM-Higgs state at the same mass, whereas the diphoton branching ratio is enhanced by up to a factor ∼ 10. This is caused by the mixing between the three neutral Higgs states, which in this parameter region gives rise to a very small H 0 d component of the mostly singlet-like state implying suppressed couplings to downtype quarks (and leptons). Thus, the mechanism leading to a large diphoton signal in this scenario     1.2 · 10 −7 8.8 · 10 −4 8.4 · 10 −6 7.9 · 10 −4 1.2 · 10 −6 7.9 · 10 −4 Γ ZZ 1.4 · 10 −8 1.1 · 10 −4 1.0 · 10 −7 9.7 · 10 −5 1.5 · 10 −7 9.9 · 10 −5 Γ gg 2.6 · 10 −6 2.5 · 10 −4 2.2 · 10 −5 2.2 · 10 −4 2.8 · 10 −5 2.1 · 10 −4 Γ γγ 9.0 · 10 −8 9.0 · 10 −6 6.3 · 10 −7 7.6 · 10 −6 7.2 · 10 −7 7.  is purely driven by the decay properties of the light singlet state. This scenario would therefore give rise to a γγ signal at the LHC without a corresponding bb excess at LEP. In fact, the reverse is also possible: an NMSSM-Higgs singlet could cause an excess in the e + e − → Z(H → bb) channel without generating a significant signal in pp → H → γγ. A relatively large bb decay rate can naturally occur for a Higgs state in the 100 GeV mass range, in which case the diphoton branching fraction remains small.
We now turn to another scenario in the low-tan β regime with large λ: the chosen input is provided in column ' Fig. 15' of Tab. 1. We vary µ eff in a narrow interval. The masses of the two lightest Higgs states are shown in Fig. 14. In the lower range of values for µ eff , the mixing between the light singlet at ∼ 101 GeV and the SM-like state at ∼ 121 GeV almost vanishes. These points are in tension with the measured properties of the SM-like state, as tested with HiggsSignals, because of the relatively low mass of the SM-like state. However, with growing µ eff , the mixing between the two light CP-even states increases, eventually pushing the singlet mass down to ∼ 90 GeV and the mass of the SM-like state up to ∼ 128 GeV. Consistency with the experimental results obtained on the observed state at 125 GeV is achieved for a mass of the SM-like state that is compatible with the LHC discovery within experimental and theoretical uncertainties. The CP-odd singlet has a mass of ∼ 150 GeV, while the heavy doublet states are at ∼ 1 TeV in this scenario. The decay properties of h 1 are documented in Fig. 15, as well as in the column ' Fig. 15' of Tab. 2 (for a specific point). The branching ratio into bb changes very significantly between the treelevel and the one-loop approach: again, the point of vanishing   only moderately enhanced with respect to the SM branching fraction due to an H 0 u -dominated doublet composition of h 1 , while BR[h 1 → bb] remains dominant, albeit slightly suppressed. This scenario would thus simultaneously address the LEP and the CMS excesses in a phenomenologically consistent manner.
Finally, we consider a CP-violating scenario, still in the large-λ, low-tan β regime. We scan over φ κ ∈ − π 8 , π 8 . The lightest Higgs state is dominantly CP-odd, with a mass of ∼ 100 GeV in the CP-conserving limit. The second-lightest Higgs is SM-like, with a mass of ∼ 120 GeV for φ κ = 0. With increasing |φ κ |, these states mix and the masses draw apart, reaching ∼ 60 GeV and ∼ 150 GeV at |φ κ | ∼ π 8 , which can be seen in Fig. 16. Correspondingly, appropriate Higgs properties, as tested with HiggsBounds and HiggsSignals, are obtained for |φ κ | 0.1-0.2. The CP-even singlet and the heavy doublet states have masses of the order of 210 GeV and 1180 GeV, respectively. In Fig. 17, we show the branching ratios of the lightest Higgs state as a function of its mass. Contrarily to the previous cases, BR[h 1 → bb] is nearly constant over the whole range; the tree-level and one-loop results agree reasonably well with each other. The branching ratios into gauge bosons show an abrupt decrease near M h 1 ∼ 102 GeV: this corresponds to the CP-conserving limit of our scenario-in which case h 1 is a pure CP-odd state. The estimated signals in the e + e − → Z(h 1 → bb) and pp → h 1 → γγ channels reach ∼ 20-25% of their SM counterparts at M h 1 ∼ 95 GeV. As can be seen in Tab. 2, the decays of h 1 for this point (|φ κ | ∼ 0.14) approximately stay in SM-like proportions. The bb and γγ signals would thus remain comparable in magnitude. In particular, a diphoton signal as large as 60% of the one for an SM-Higgs boson could not be accommodated in this configuration, as the LEP limits on h 1 → bb indirectly constrain the diphoton rate. Yet, it is remarkable that a mostly CP-odd Higgs state would trigger sizable signals at both LEP and the LHC.
To summarize this discussion, a light CP-even-or dominantly CP-even and even dominantly CPodd in the CP-violating case-and mostly singlet-like Higgs state in the vicinity of 100 GeV could have interesting consequences for the phenomenology of the SM-like state. Sizable signals in the e + e − → Z(h 1 → bb) and/or pp → h 1 → γγ channels are possible, independently or simultaneously. They could thus explain the corresponding 'excesses' reported by LEP and CMS, respectively. Further experimental effort is required to investigate this scenario. 8

Discussion concerning the remaining theoretical uncertainties
Below, we provide a summary of the main sources of theoretical uncertainties from unknown higherorder corrections applying to our calculation of the NMSSM Higgs decays. We do not discuss here the parametric theoretical uncertainties arising from the experimental errors of the input parameters. For the experimentally known SM-type parameters the induced uncertainties can be determined in the same way as for the SM case (see e. g. Ref. [21]). The dependence on the unknown SUSY parameters, on the other hand, is usually not treated as a theoretical uncertainty but rather exploited for setting indirect constraints on those parameters.
Higgs decays into quarks (h i → qq, q = c, b, t) In our evaluation, these decays have been implemented at full one-loop order, i. e. at QCD, electroweak and SUSY next-to-leading order (NLO). In addition, leading QCD logarithmic effects have been resummed within the parametrization of the Yukawa couplings in terms of a running quark mass at the scale of the Higgs mass. The Higgs propagator-type corrections determining the mass of the considered Higgs particle as well as the wave function normalization at the external Higgs leg of the process contain full one-loop and dominant two-loop contributions.    For an estimate of the remaining theoretical uncertainties, several higher-order effects should be taken into account: • First, we should assess the magnitude of the missing QCD NNLO (two-loop) effects. We stress that there should be no large logarithms associated to these corrections, since these are already resummed through the choice of running parameters and the renormalization scale. For the remaining QCD pieces, we can directly consider the situation in the SM. In the case of the light quarks, the QCD contributions of higher order have been evaluated and amount to ∼ 4% at m H = 120 GeV (see e. g. Ref. [136]). For the top quark, the uncertainty due to missing QCD NNLO effects was estimated to 5% [21].
• Concerning the electroweak corrections, Fig. 1 suggests that the one-loop contribution is small-at the percent level-for an SM-like Higgs, which is consistent with earlier estimates in the SM [21]. For the heavy Higgs states, Fig. 1 indicates a larger impact of such effects-at the level of ∼ 10% in the considered scenario. Assuming that the electroweak NNLO corrections are comparable to the squared one-loop effects, our estimate for pure electroweak higher orders in decays of heavy Higgs states reaches the percent level. In fact, for multi-TeV Higgs bosons, the electroweak Sudakov logarithms may require a resummation. Furthermore, mixed electroweak-QCD contributions are expected to be larger than the pure electroweak NNLO corrections, adding a few more percent to the uncertainty budget. For light Higgs states, the electroweak effects are much smaller since the Sudakov logarithms remain of comparatively modest size.
• Finally, the variations with the squark masses in Fig. 1 for the heavy doublet states show that the one-loop SUSY effects could amount to 5-10% for a sub-TeV stop/sbottom spectrum. In such a case, the two-loop SUSY and the mixed QCD/electroweak-SUSY corrections may reach the percent level. On the other hand, for very heavy squark spectra, we expect to recover an effective singlet-extended Two-Higgs-Doublet model (an effective SM if the heavy doublet and singlet states also decouple) at low energy. However, all the parameters of this low-energy effective field theory implicitly depend on the SUSY radiative effects, since unsuppressed logarithms of SUSY origin generate terms of dimension ≤ 4-e. g. in the Higgs potential or the Higgs couplings to SM fermions. On the other hand, the explicit dependence of the Higgs decay widths on SUSY higher-order corrections is suppressed for a large SUSY scale. In this case, the uncertainty from SUSY corrections reduces to a parametric effect, that of the matching between the NMSSM and the low-energy lagrangian-e. g. in the SM-limit, the uncertainty on the mass prediction for the SM-like Higgs continues to depend on SUSY logarithms and would indirectly impact the uncertainty on the decay widths.
Considering all these higher-order effects together, we conclude that the decay widths of the SMlike Higgs should be relatively well controlled (up to ∼ 5%), while those of a heavy Higgs state could receive sizable higher-order contributions, possibly adding up to the level of ∼ 10%.
Higgs decays into leptons Here, QCD corrections appear only at two-loop order in the Higgs propagator-type corrections as well as in the counterterms of the electroweak parameters and only from three-loop order onwards in the genuine vertex corrections. Thus, the theory uncertainty is expected to be substantially smaller than in the case of quark final states. For an SM-like Higgs, associated uncertainties were estimated to be below the percent level [23]. For heavy Higgs states, however, electroweak one-loop corrections are enhanced by Sudakov logarithms and reach the ∼ 10% level for Higgs masses of the order of 1 TeV, so that the two-loop effects could amount to a few percent. In addition, light staus may generate a sizable contribution of SUSY origin, where the unknown corrections are of two-loop electroweak order.
Higgs decays into W W/ZZ The complexity of these channels is illustrated by our presentation of two separate estimates, expected to perform differently in various regimes.
• In the SM, the uncertainty of Prophecy4f in the evaluation of these channels was assessed at the sub-percent level below 500 GeV, but up to ∼ 15% at 1 TeV [21]. For an SM-like Higgs, our Fig. 3 shows that the one-loop electroweak corrections are somewhat below 10%, making plausible a sub-percent uncertainty on the results employing Prophecy4f. On the other hand, the assumption that the decay widths for an NMSSM Higgs boson can be obtained through a simple rescaling of the result for the width in the SM by tree-level couplings, is in itself a source of uncertainties. We expect this approximation to be accurate only in the limit of a decoupling SM-like composition of the NMSSM Higgs boson. If these SM-like characteristics are altered through radiative corrections of SUSY origins or NMSSM-Higgs mixing effects-both of which may still reach the level of several percent in a phenomenologically realistic setup-the uncertainty on the rescaling procedure for the decay widths should be of corresponding magnitude.
• In the case of heavier states, Fig. 3 and Fig. 8 indicate that the previous procedure is unreliable in the mass range 500 GeV. In particular, for heavy doublets in the decoupling limit, radiative corrections dominate over the-then vanishing-tree-level amplitude, shifting the widths by orders of magnitude. In such a case, our one-loop calculation captures only the leading order and one can expect sizable contributions at the two-loop level: as we already mentioned in discussing Fig. 3, shifting the quark masses between pole and MS values-two legitimate choices at the one-loop order that differ in the treatment of QCD two-loop contributionsresults in modifications of the widths of order ∼ 50%. On the other hand, one expects the decays of a decoupling heavy doublet into electroweak gauge bosons to remain a subdominant channel, so that a less accurate prediction may be tolerable. It should be noted, however, that the magnitude of the corresponding widths is sizably enhanced by the effects of one-loop order, which may be of interest regarding their phenomenological impact.
Radiative decays into gauge bosons As these channels appear at the one-loop order, our (QCDcorrected) results represent (only) an improved leading-order evaluation. Yet the situation is contrasted: • In the SM, the uncertainty on a Higgs decay into γγ was estimated at the level of 1% in Ref. [21]: however, the corresponding calculation includes both QCD NLO and electroweak NLO corrections. In our case, only QCD NLO corrections (with full mass dependence) are taken into account. The comparison with NMSSMCALC in Fig. 9 provides us with a lower bound on the magnitude of electroweak NLO and QCD NNLO effects: both evaluations are at the same order but differ by a few percent. The uncertainty on the SUSY contribution should be considered separately, as light charginos or sfermions could have a sizable impact.
In any case, we expect the accuracy of our calculation to perform at the level of 4% (the typical size of the deviations in Fig. 9).
• In the case of the Higgs decays into gluons, for the SM prediction-including QCD corrections with full mass dependence and electroweak two-loop effects-an uncertainty of 3% from QCD effects and 1% from electroweak effects was estimated in Ref. [21]. In our case, the QCD corrections are only included in the heavy-loop approximation, and NLO electroweak contributions have not been considered. Consequently, the uncertainty budget should settle above the corresponding estimate for the SM quoted above. In the case of heavy Higgs bosons, the squark spectrum could have a significant impact on the QCD two-loop corrections, as exemplified in Fig. 5 of Ref. [117].
• For h i → γZ, QCD corrections are not available so far, so that the uncertainty should be above the ∼ 5% estimated in the SM [21].
Additional sources of uncertainty from higher orders For an uncertainty estimate, the following effects apply to essentially all channels and should be considered as well: • The mixing in the Higgs sector plays a central role in the determination of the decay widths. Following the treatment in FeynHiggs, we have considered Z mix in all our one-loop evaluations, as prescribed by the LSZ reduction. Most public codes consider a unitary approximation in the limit of the effective scalar potential (U 0 , in our notation). The analysis of Ref. [25] and our discussion on Fig. 1-employing U m , a more reliable unitary approximation than U 0indicate that the different choices of mixing matrices may affect the Higgs decays by a few percent (and far more in contrived cases). However, even the use of Z mix is of course subject to uncertainties from unknown higher-order corrections. While the Higgs propagator-type corrections determining the mass of the considered Higgs boson and the wave function normalization contain corrections up to the two-loop order, the corresponding prediction for the mass of the SM-like Higgs still has an uncertainty at the level of about 2%, depending on the SUSY spectrum.
• In this paper, we confined ourselves to the evaluation of the Higgs decay widths into SM particles and did not consider the branching ratios. For the latter an implementation at the full one-loop order of many other two-body decays, relevant in particular for the heavy Higgs states, would be desirable, which goes beyond the scope of the present analysis. Furthermore, in order to consider the Higgs branching ratios at the one-loop order, we would have to consider three-body widths at the tree level, for instance h i → bbZ, since these are formally of the same magnitude as the one-loop effects for two-body decays. In addition, these three-body decays-typically real radiation of electroweak and Higgs bosons-exhibit Sudakov logarithms that would require resummation in the limit of heavy Higgs states.
• At decay thresholds, the approximation of free particles in the final state is not sufficient, and a more accurate treatment would require the evaluation of final-state interactions. Several cases have been discussed in e. g. Refs. [75,76,137].
In this discussion we did not attempt to provide a quantitative estimate of the remaining theoretical uncertainties from unknown higher-order corrections, as such an estimate would in any case sensitively depend on the considered region in parameter space. Instead, we have pointed out the various sources of higher-order uncertainties remaining at the level of our state-of-the-art evaluation of the Higgs decays into SM particles in the NMSSM. For a decoupling SM-like Higgs boson one would ideally expect that the level of accuracy of the predictions approaches the one achieved in the SM. However, even in this limit, missing NNLO pieces-that are known for the SM, but not for the NMSSM-give rise to a somewhat larger theoretical uncertainty in the NMSSM. Furthermore, uncertainties of parametric nature (for instance from the theoretical prediction of the Higgs-boson mass) need to be taken into account as well. For heavy Higgs states, the impact of electroweak Sudakov logarithms and SUSY corrections add to the theoretical uncertainty to an extent that is strongly dependent on the details of the spectrum and the characteristics of the Higgs state. For a decoupling doublet at ∼ 1 TeV, an uncertainty of ∼ 5-15% may be used as a guideline for the fermionic and radiative decays, while the uncertainty may be as large as ∼ 50% in h i → W W/ZZ.

Conclusions
In this paper, we have presented our evaluation of neutral Higgs decay widths into SM final states in the (CP-conserving or CP-violating) NMSSM. Full one-loop corrections have been included for all the considered channels, as well as higher-order QCD corrections to the decays that are generated at the radiative level. The inclusion of one-loop contributions to the decays into SM fermions or electroweak gauge bosons goes beyond the usual approximation amounting to a QCD/QEDcorrected tree level. In addition, QCD corrections to the digluon and diphoton decay widths have been carefully processed, including the mass dependence in the γγ case and corrections beyond NLO in the gg case. In its current form, this state-of-the-art implementation of the neutral Higgs decays into SM particles is available as a Mathematica package, but should also be integrated into the FeynHiggs code in the near future.
In order to illustrate this calculation of the Higgs decay widths, we have presented our results in several regimes of the parameter space of the NMSSM. In the MSSM limit, we were able to recover the predictions of FeynHiggs and trace the origins of deviations from our new results. In particular, we emphasized the relevance of one-loop contributions in the decays of heavy doublet states into electroweak gauge bosons, for which the usual estimates based on the tree-level Higgsgauge couplings are not appropriate. Minor effects in the treatment of QCD and QED corrections have also been noted. Beyond the MSSM limit, we have compared our decay widths to the output of NMSSMCALC. We observed a qualitative agreement wherever this could be expected. We also gave an account of the various sources of theoretical uncertainties from higher-order corrections and discussed the achieved accuracy of our predictions.
As a phenomenological application, we investigated in particular the case of a mostly singlet-like state with mass in the vicinity of 100 GeV. The decays of such a state can be notably affected by suppressed couplings to down-or up-type quarks which can occur in certain parameter regions as a consequence of the mixing between the different Higgs states. In particular, an additional Higgs boson h i of this kind could manifest itself via signatures in the channels e + e − → Z(h i → bb) and/or pp → h i → γγ. The presence of such a light Higgs boson could thus explain the slight deviations from the SM predictions reported by LEP and CMS in those channels.