About the bosonic decays of heavy Higgs states in the (N)MSSM

The heavy, doublet-dominated Higgs bosons expected in Two-Higgs-Doublet-Model-like extensions of the Standard Model are obvious targets for searches at high-energy colliders, and considerable activity in this sense is currently employed in analyzing the results of the Large Hadron Collider. For sufficiently heavy states, the $SU(2)_{\mathrm{L}}$ symmetry drastically constrains the decays of these new scalars. In models that contain only additional Higgs doublets, fermionic decay channels are expected to dominate, although some bosonic modes may lead to cleaner signals. The situation is very different if the Higgs sector is further extended by singlet states and the singlet-dominated scalars are kinematically accessible, in which case the Higgs-to-Higgs widths may even become dominant. Yet, a quantitative interpretation of the experimental data in terms of a specific model also runs through the control of radiative corrections in this model. In this paper, we examine the bosonic decays of the heavy, doublet-dominated Higgs states in the MSSM and the NMSSM at full one-loop order, insisting on the impact of the $SU(2)_{\mathrm{L}}$ symmetry and the effects of infrared type, in particular, the emergence of sizable non-resonant contributions to the three-boson final states.


Introduction
and decay calculations, avoiding large spurious effects associated with an artificial dependence on gauge-fixing parameters and field renormalization, as commonly found in earlier estimates. In Ref. [35], we examined the Higgs phenomenology in the presence of very light singlet states, showing in particular that a tachyonic tree-level spectrum does not necessarily yield an unphysical point in parameter space, but may simply correspond to a problematic choice of renormalization scheme. With the current paper focusing on bosonic decays of heavy, doublet-like Higgs states, we conclude this overview of radiative corrections applying to Higgs decays into SM final states in the (N)MSSM. Similar projects were presented in Refs. [36][37][38][39][40][41][42][43][44][45][46][47].
Our implementation of the bosonic Higgs decays narrowly follows the steps outlined in Refs. [32][33][34][35] and we simply discuss the consequences of the 1L corrections to these channels for the Higgs phenomenology in two scenarios illustrative of the MSSM and the NMSSM. In particular, we attempt to quantify higher-order uncertainties, which may remain considerable for some channels at 1L due to a suppressed tree-level width. We first consider bosonic Higgs decays in the MSSM in Sect. 2, insisting on the implications of the approximate SU (2) L symmetry for both two-body and three-body decays. We then turn to the Z 3 -conserving NMSSM in Sect. 3, showing how channels with singlet scalars in the final state may come to dominate the total width. We briefly summarize the achievements of this paper in Sect. 4.

The bosonic Higgs decays in the MSSM
We first focus on the status of bosonic decays of heavy Higgs states in the MSSM. We denote the CP-even Higgs states as h = −s α h 0 d + c α h 0 u and H = c α h 0 d + s α h 0 u , the CP-odd one as A = s β a 0 d + c β a 0 u , and the charged ones as H ± = s β H ± d + c β H ± u , with transparent notations for the components of the doublets H d , H u , respectively taking a vacuum expectation value (v.e.v.) v c β , v s β , with v ≡ (2 √ 2 G F ) −1/2 and G F denoting the Fermi constant. The lightest state h is identified with the SM-like Higgs boson observed at the LHC. We will employ the notation M EW to represent any mass of electroweak proportions, i. e. M W,Z,h,t ∼ M EW . The masses of all heavy doublet-Higgs states scale like M H,A ∼ M H ± , which we will regard as falling above 0.5 TeV.

Two-body decays
In the MSSM, the squared-mass splitting among members of the heavy Higgs doublet (H, A, H ± ) is of order M 2 EW ∼ M 2 Z,W , implying a quasi-degeneracy at large M H ± . This is a consequence of the restored SU (2) L symmetry at energies sufficiently far above the electroweak-breaking scale. Correspondingly, the accessible decay channels of these doublet scalars into pairs of SM bosons are restricted to • decays into pairs of gauge bosons: H, A → gg, γγ, Zγ, ZZ, W + W − ; H ± → ZW ± , γW ± ; • decays into a pair of SM-like Higgs bosons: H → hh; • decays into a SM-like Higgs and an electroweak gauge boson: A → Zh; H ± → W ± h; where we have neglected a possible CP-violating mixing, which would combine the decay modes of H and A. All these channels have in common to violate the SU (2) L symmetry, resulting in a suppression of the associated widths by the ratio M EW /M H ± at large M H ± . Correspondingly, all tree-level amplitudes involving a massless gauge boson vanish. In the case of H → hh, the trilinear Higgs coupling is directly proportional to the electroweak v.e.v. v. For other channels, the tree-level amplitudes are mediated by the SU (2) L -breaking mixing angle in the CP-even sector α − β + π/2 = O M 2 EW /M 2 H ± . Given that the suppression of these amplitudes in the large-mass limit is protected by a symmetry, it is also necessarily observed by loop corrections. For this reason, heavy-Higgs decays in the MSSM are dominated by the fermionic channels, which we studied in depth in Ref. [32], while SU (2) L -conserving bosonic decay modes, such as H → ZA, are kinematically inaccessible and contribute to final states of higher multiplicity, e. g. H → Zbb, after interfering with other diagrams. We compute radiative corrections to the decay amplitudes at full 1L order, taking careby the means of a strict perturbative expansion-not to introduce SU (2) L -breaking artifacts in the calculation [33]. We will not narrowly investigate the pure radiative decays here: our implementation has been described in Ref. [35]-see also Ref. [31]-and includes known higherorder QCD corrections to the diphoton [48] and digluon [49,50] widths-see also the summary in Ref. [51]. Instead, we focus on final states involving weak gauge or lighter Higgs bosons. A few channels are displayed for illustration in Fig. 1. The SUSY spectrum is decoupled with masses at 10 TeV while t β = 10 (see Appx. B for a summary of the input parameters). At the tree level, the considered decay widths read: As expected, all these decay widths scale like M 2 EW M H ± in the limit of large doublet masses. The born-level amplitude for the H → hh decay does not vanish in the limit α → β − π 2 and is merely t β -suppressed. On the contrary, H → ZZ and A → Zh become purely radiative in this limit. For these latter channels, the alignment of the CP-even sector in the limit M H ± M EW is crucial to recover the correct behavior.
In the left column of Fig. 1, the tree-level prediction (short-dashed blue lines) for the channels of Eqs. (1) is shown against M H ± and compared to our estimates at 1L (red or green solid curves). All curves fall at large M H ± , as expected for SU (2) L -violating channels. We observe that the magnitude of 1L contributions is comparable to that of the born-level amplitudes, indicating that the perturbative series converges slowly. In these conditions, the inclusion of a |1L| 2 term (long-dashed curves) seems numerically justified. Nevertheless, such an operation is in general dangerous [33], because the |1L| 2 term is a partial 2L contribution which generally contains loose UV-logarithms or uncontrolled symmetry-violating contributions, thus spoiling a quantitative interpretation. Yet, in the case of channels that are purely radiative in the alignment limit, the 1L amplitude is purged of its undesirable features when explicitly applying the limit α → β − π 2 : the corresponding piece typically dominates at the numerical level, thus endowing the |1L| 2 term with some degree of legitimacy. In any case, comparing the strict 1L widths with their counterparts including a |1L| 2 term provides us with a first assessment of the higher-order uncertainty, which reaches as much as 100% for the considered channels.
Another measurement of the large uncertainties at stake at 1L order can be derived from the variation of input in the quark sector: quark masses indeed enter the calculation only in 1L diagrams-both in triangles and self-energy corrections on the external legs-so that a choice of scheme input appears as a higher-order concern. Thus, we compare an on-shell input (green   colors) and an MS input with QCD running at the scale of the heavy Higgs masses (red/orange shades). Once again, the variation reaches a 100% magnitude, due to the substantial weakening of Yukawa interactions at high energy. Specializing on the example of H → hh, the contribution of the top quark to the 1L amplitude is of order g (1L),t implying a quartic dependence on the quark mass input, hence a dramatic dependence on QCD corrections of higher order. In practice, the choice of MS masses at the high scale is expected to properly resum the associated logarithms of UV type, and we correspondingly choose it as our prediction per default. Nevertheless, our discussion should emphasize the impossibility to claim accurate results in the prediction of the bosonic Higgs decays in a 1L analysis: all that is achieved is setting the order of magnitude of the suppressed two-body widths-which is anyway constrained by the symmetries and dimensional analysis.
Despite these bleak conclusions as to the possibility of quantitatively exploiting the bosonic two-body decays of a heavy doublet-Higgs state in the MSSM-they are suppressed and imprecise at 1L-it is crucial to properly estimate their magnitude. On the right-hand side of Fig. 1, still for the channels H → hh, H → ZZ and A → Zh, we display the tree-level (short dashed blue line) and 1L (solid red line) estimates-identical to the corresponding curves on the left plots-and compare them to the predictions of FeynHiggs-2.18.1 [52][53][54][55][56][57][58][59] (long dashed purple curves). FeynHiggs employs a partial 1L description of the considered widths, where the treelevel amplitude is rotated with an effective mixing matrix diagonalizing the Higgs-mass system at 1L (and leading 2L) order [60]. In the case of H → hh, this approach does not particularly improve the prediction with respect to a tree-level description (as it is not closer to the full 1L result) but, in view of the sizable uncertainties at stake, it does not notably worsen it either, with, in particular, the correct decoupling behavior at high mass.
However, if we turn to H → ZZ, we observe a problematic growth of the width predicted by FeynHiggs for M H ± M EW . This response is the direct consequence of the inclusion of a partial 1L order containing SU (2) L -violating pieces that are not controlled by the electroweak v.e.v.: non-decoupling mixing effects actually need to be carefully balanced with the vertex piecewhich obviously fails if the latter is not included. We already described this issue in Ref. [33] (see e. g. Fig. 10 in this reference). At this point, one may wonder why the same cause does not result in similarly devastating consequences for H → hh: the reason simply rests with the fact that, in this latter case, the SU (2) L -breaking parameter is already contained within any of the (tree-level) triple-Higgs couplings, i. e. g hhh and g Hhh both scale like M EW ; on the contrary, for H → ZZ, the hZZ coupling is unsuppressed, so that a non-decoupling behavior emerges from a call to this quantity via a (non-decoupling) effective mixing. A similar issue appears in A → Zh, with the difference that the non-decoupling mixing is introduced at the level of the SM-like state, through a call to the SU (2) L -conserving AZH coupling.
Finally, we warn the reader against computing off-shell decays H → Z * A * , A → Z * H * , etc., as an integration over off-shell momenta of the final state bosons, after the fashion of Eq. (117) of Ref. [51]. Such expressions are meant to estimate the contribution to four-fermion final states, but the latter involve several destructively interfering diagrams in the case of heavy doublet-Higgs states-with e. g. Z-boson emission from a fermion instead of a Higgs line. In our approach, corresponding effects are described in terms of three-body decay widths, e. g. A → ff Z, ff h, including all tree-level contributions and resumming Sudakov double logarithms (which appear in phase-space integrals) [32]: this results in a less dramatic behavior at high mass, consistent with the SU (2) L symmetry (contrarily to the off-shell approach). In the case of decays above threshold, such as A → Zh, use of the integration formula is possible, since the on-shell contribution dominates the integral; on the other hand, it also brings little, far from threshold, as compared to a simple on-shell calculation.

Three-body decays
Order counting sets a three-body width at tree level on the same footing as 1L contributions to a two-body decay. When enough phase space is available, it is thus impossible to perform a complete 1L analysis without studying three-body widths-even after discarding possible contributions with an on-shell internal line, which we count as two body. In the case of fermionic decays, we have shown in Ref. [32] how IR effects (Sudakov logarithms) compensate between three-body and two-body channels, resulting in sizable effects when assessing the full width or branching ratios. Here we focus on final states that strictly involve bosonic particles. For the MSSM, the relevant channels are • two gauge bosons and one Higgs boson: • one gauge boson and two Higgs bosons: A → hhZ; H ± → hhW ± ; • three Higgs bosons: H → hhh; where we have neglected CP-violating mixing in the neutral sector (indeed absent at tree level) as well as photon radiation-which we include in the two-body widths so that these are IR-safe with respect to QED.
Contrarily to the two-body decays, many of these channels are SU (2) L -conserving-it is indeed possible to build a doublet out of e. g. three doublets, or one doublet and two triplets. Therefore, no suppression associated to the gauge group arises at the amplitude level, and these channels only entail a three-body phase-space suppression, as well as a t β -suppression: the corresponding widths then scale linearly with M H ± . For example, in the limit M H ± M EW , the H → hhh width is dominated by the contribution of the quartic Higgs coupling: Due to this behavior, the three-body widths involving bosonic final states overtake their two-body counterparts as soon as M H ± reaches a few TeV. The decays of the neutral heavy-doublet states are shown for illustration in Fig. 2. For each state, CP-even or odd, three channels are SU (2) L -conserving and the associated widths grow linearly with the Higgs mass, while one-channel is SU (2) L -violating, hence suppressed-CPviolating channels vanish. One can easily identify the status of each channel with respect to the SU (2) L symmetry by checking the existence of a corresponding quartic coupling in the tree-level lagrangian; in such a game, gauge bosons should be identified with their associated Goldstone bosons. Furthermore, the behavior at large M H ± is mostly determined by these quartic couplings and one deduces the following relations in this limit: The theoretical uncertainty in all these channels can be expected to be sizable (of order 100%), since the widths are evaluated at leading order-consistently with the 1L order for two-body decays-and the Higgs potential is known to receive sizable radiative corrections from e. g. Yukawa interactions. Therefore, as for the two-body decays, an evaluation at the considered order only fixes the general magnitude of the three-body widths, without allowing for precision tests.
The decays of the charged Higgs bosons are constrained by the SU (2) L symmetry to follow a pattern similar to the one of the neutral states, hence these channels exhibit limited novelty. For completeness, we collect a few plots in Appx. A.

Branching ratios
Having estimated all Higgs decays into SM two-body final states at 1L and three-body final states at tree level, it is possible to consistently evaluate the full widths and branching ratios at this order in scenarios where exotic particles (SUSY) are kinematically inaccessible.
The total widths for the heavy neutral states at t β = 10 are shown in the first row of Fig. 3. As expected from the SU (2) L symmetry, the predicted widths are comparable for both states, either at the tree level (short-dashed lines) or at the loop level (long-dashed lines). The results of FeynHiggs (solid curves) do not satisfy this property at large M H ± , due to the issues mentioned earlier in the assessment of bosonic decay widths, mainly. In addition, the apparent agreement of the expectations of FeynHiggs for the CP-even state (in green) with our full 1L result is in fact coincidental. Given that FeynHiggs does not consider three-body decays, its predictions should be compared to the pseudo-widths, where only 1L two-body widths are summed (not IR-safe with respect to weak corrections; plotted in dashes of intermediate size): the mismatch is of order 15% and can be explained, in part by the problematic H → ZZ, W W widths, but also by the contamination of the fermionic widths, on the side of FeynHiggs, by symmetry-violating terms of higher order-we refer the reader to the more detailed discussion of Ref. [33].
The branching ratios at full 1L order are shown in the last row of Fig. 3 and can be compared to the (QCD-corrected) tree-level version (in the middle). The main difference is related to the emergence of the three-body fermionic decays (reaching O(10%) at large M H ± ; dashed orange lines), which complement the dominant two-body fermionic channels (in solid red). As explained in Ref. [32], these modes essentially involve the radiation of electroweak gauge bosons from the initial Higgs or fermionic final states and are dominated by double logarithms of Sudakov type. Bosonic branching ratios (blue and green curves) remain below percent level-and even permil, for M H ± 2 TeV. The three-body bosonic decays, in particular those modes involving at least one SM-like Higgs boson in the final state (dashed cyan lines), overtake the bosonic two-body decays for M H ± 4 TeV. The linear growth of the associated widths with the Higgs mass ensures that such channels are not further suppressed at large M H ± , contrarily to the SU (2) L -violating two-body modes.
Bottom: Branching ratios at 1L: three-body fermionic (dashed orange) and bosonic (involving one SM-like Higgs; in dashed cyan) modes, as well as the two-body gluonic rate (dashed magenta) are added to the picture.
The picture described above remains comparatively stable under variation of t β , with effects associating to the top Yukawa coupling at low t β replacing those of the bottom Yukawa at high t β . Therefore, as a tribute to the SU (2) L -suppression, bosonic decays of the heavy Higgs states always represent sub-percent-level channels for masses above TeV in the MSSM. Of course, if SUSY final states become accessible, some of them are SU (2) L -conserving and may compete with the SM-fermionic channels.

The bosonic Higgs decays in the NMSSM
The presence of singlet degrees of freedom in the NMSSM may noticeably affect the phenomenology of heavy doublet-Higgs states, in particular if these new scalars are kinematically accessible. Indeed, the heavy Higgs bosons may then decay into bosonic final states without breaking the SU (2) L symmetry, hence opening channels that are competitive with respect to the fermionic ones. Below, we denote the singlet-dominated CP-even and odd scalars as h S and a S respectively.

Two-body decays
When h S and/or a S take a mass below M H ± , the following two-body decay channels may add up to the MSSM modes: • decays into two Higgs bosons: • decays into a Higgs and a gauge boson: Again, we have assumed CP-conservation for a cleaner classification, but the inclusion of CPviolating mixings does not fundamentally change the situation. Among these new channels, several are again SU (2) L -violating, in substance all those involving only singlet scalars in the final state. However, in contrast to the MSSM, some other decay modes-involving one singlet and one doublet or weak gauge boson (i. e. its Goldstone component)-are SU (2) L -conserving and need not be small in the limit M H ± M EW provided the 'portal' between singlet and MSSM sectors, i. e. the coupling λ, is sufficiently intense.
We consider a scenario with a comparatively light CP-odd singlet-dominated Higgs boson, setting its mass to 300 GeV, while M H ± ! = 3 TeV. The SUSY sector (as well as the CP-even singlet) remains decoupled, taking masses in the 10 TeV range. All input parameters are summarized in Appx. B. The singlet a S primarily decays into fermion pairs, so that the two-body decays of the heavy doublet states that we study below can be understood as the resonant contribution to four-fermion decay channels. In Fig. 4, we vary λ = κ between low values-where NMSSM effects are negligible and the singlet decouples-and 0.2 with an active connection between singlet and MSSM sectors. Large radiative corrections to the mass of the singlet pseudoscalar develop for sizable values of λ = κ: this is the simple consequence of the quadratic dependence of scalar masses on the UV spectrum below the SUSY-breaking scale. Technically, this causes complications with the requirement of M a S ! = 300 GeV in our original renormalization scheme (DR with the renormalization scale M t ), as the pseudoscalar mass may be vastly different at the tree level and at 1L. In fact, a 'tachyonic tree-level syndrome'-see Ref. [35]-emerges as soon as λ 0.04, making evaluations at the 1L order impossible in this scheme. For this reason, it is convenient to directly work in the renormalization scheme where the counterterm for A κ is fixed by an on-shell condition on the mass of a S [35], after the tree-level input is correspondingly translated: the bare A κ is kept unchanged in the scheme translation.  We then focus on the two SU (2) L -conserving channels H → Za S and A → ha S . At the tree level, one may express the widths as In the SU (2) L limit with M EW M H ± , one can directly work out Obviously, g Aha S does not vanish and may reach sizable values for κ, λ = O(1). This approximation is less immediate for g HZa S : with X R,I h i d,u denoting the components of the Higgs field h i with respect to the scalar and pseudoscalar doublet gauge eigenstates; the gauge interaction couples only doublet, not singlet states. Thus, the H-Z-a S coupling necessarily emerges from singlet-doublet mixing and is SU (2) Lviolating. Yet, the prefactor M −2 Z in Eq. (4a) shows that the relevant object to consider in the SU (2) L limit is not g HZa S but g HZa S M Z . Then, it is convenient to use the connection between gauge and Goldstone couplings, providing the relation |g HZa S | M Z = 2 |g HG 0 a S | (M 2 H − M 2 a S ), and to apply the SU (2) L limit to the Goldstone coupling, i. e.
At this point, we see that the H → Za S transition actually amounts to a Higgs-to-Higgs decay in the SU (2) L limit, the Z boson serving as proxy for its Goldstone counterpart. In addition, both channels of Eq. (4) are related by the SU (2) L symmetry, so that the widths are approximately equal. Radiative corrections are expected to abide by such a symmetry-protected relation.
The decay widths are displayed in the first row of Fig. 4 at the tree level (blue short-dashed line) and at 1L for different definitions of the Yukawa couplings (solid green and red curves). The width corresponding to the sum of the MSSM-like decay modes is also shown (long-dashed gray line) and remains comparatively stable with varying κ = λ. Contrarily to the MSSM case, Higgsto-Higgs decays are found to dominate over the fermionic modes, provided λ = κ is sufficiently large. The reason for this behavior is the dependence of g Aha S ,HG 0 a S on the scale κ µ eff , which can compete with Y f M H ± in the fermionic case (with Y f denoting the Yukawa coupling of the fermion species f ). Now, focusing on the bosonic decay widths, we observe that these are largely determined by the tree-level amplitudes; 1L corrections only amount to O(10%). This again contrasts with the MSSM-like bosonic channels. A fast convergence of the perturbative series is expected and the inclusion of |1L| 2 terms appears unjustified (and unpredictive) in such a context. The choice of input for the fermion couplings also appears largely unimportant, with an impact of O(1%).
In the second row of Fig. 4, we test the scheme dependence in Γ[A → ha S ]. In the plot on the left, we study the impact of the treatment of A κ , considering the variations between a calculation in the scheme with on-shell M a S (black solid line), the original DR scheme (dotted red) and several intermediate choices. The 'collapse' of the predictions for non-OS choices is related to the tree-level mass m a S becoming tachyonic; this issue should not distract our attention too much: perturbative calculations are expected to fail for an ill-defined spectrum, and the OS choice is a priori more predictive. More interestingly, in the regime where these schemes are well-defined, their predictions for the decay width do not differ from that obtained with the OS choice by more than a few percent-and this magnitude is reached only in the limit of validity of the non-OS schemes. This confirms the stability of the predicted width under varying definitions of the renormalized A κ . Yet, given the form of the tree-level couplings-see the paragraph after Eq. (4)-we expect the main uncertainty to originate in the definition of κ (or µ eff ). Thus, we also vary the scheme associated with κ in the plot on the righthand side: keeping the bare κ unchanged, we consider variations of the renormalized parameter by ∆κ ≈ 3 (16 π 2 ) κ (λ 2 + κ 2 ) ln M 2 SUSY M 2 t , as suggested by the anomalous dimension at 1L. The widths at 1L (dashed green curves) are affected at permil level, at most, against percent for the born level predictions (dotted red lines). This is another sign of the convergence of the perturbative series and we can conclude as to a higher-order uncertainty by at most a few percent on the SU (2) L -conserving channels. The approximate SU (2) L relation between the two channels H → Za S and A → ha S is roughly satisfied in Fig. 4. Violating effects amount to at most ∼ 5% at the tree level and ∼ 10% at 1L (in this specific scenario). This apparently increasing magnitude of the SU (2) L -breaking at 1L should not be over-interpreted: indeed, at this order, three-body decays should be counted at the same level as 1L corrections to the two-body decays. In the next subsection, we will find that the larger three-body decay width of A somewhat compensates the larger suppression of Γ[A → ha S ] as compared to Γ[H → Za S ].
In Fig. 5, we consider the SU (2) L -violating decay H → a S a S (a channel that does not exist in the THDM) in the same scenario as before. As expected, this width is considerably suppressed as compared to those of Fig. 4, though growing with λ = κ. We actually recover features that are largely similar to the MSSM-like Higgs-to-Higgs decays, with a slow convergence of the perturbative series and a sizable dependence on the definition of the SM Yukawa couplings, all hinting at an uncertainty of the order of 100% at 1L. On the other hand, it is always possible to tailor scenarios where a large mixing between H and h S ensures a sizable decay of their admixtures into a S a S -exploiting the SU (2) L -conserving decay h S → a S a S ; see e. g. Ref. [61]. As the approximate degeneracy between singlet and doublet-dominated states is accidental, such setups are generally fine-tuned.
Attempts at a comparison with the public tools NMSSMCALC [62] and NMSSMCALCEW [45] did not lead to satisfactory results. First, we stress that a preliminary translation of the input is needed, since the renormalization scheme in these codes is different from our choice. Then, electroweak corrections do not seem to produce noticeable IR-like effects in the predictions of NMSSMCALCEW, neither in fermionic nor bosonic two-body modes, which contrasts with our calculation. On the other hand, contrarily to FeynHiggs, we did not observe any obvious violation of the SU (2) L symmetry in the results of NMSSMCALC and NMSSMCALCEW. Tree-level predictions of bosonic two-body decays are in rough agreement with ours. A deeper understanding of the two implementations of decay widths would thus be needed for a meaningful comparison of the results, which goes beyond our purpose in the current paper.

Three-body decays
Numerous bosonic three-body final states open up in addition to the THDM ones, when a singletdominated (pseudo)scalar is kinematically accessible to the decays of heavy Higgs bosons. It is again necessary to calculate these widths at tree level for a meaningful 1L analysis of the total widths and branching ratios of Higgs bosons. There exist fermionic three-body decay channels involving h S and a S as well: these are typically dominated by the emission of a fermion pair via an off-shell h, G 0 or G ± in the two-body topology of the previous subsection-see Ref. [32]. We do not study them in detail below.
We consider the three-body decay widths into bosonic final states in the scenario of Fig. 4. As before, we exclude the contributions from resonant internal lines, which are counted in the two-body widths, and keep only off-shell effects at this level. We do not discuss the modes of THDM type here, since they continue to behave in a fashion comparable to that observed in the MSSM. Instead, we focus on final states involving a pseudoscalar a S . The results are displayed in Fig. 6. Some channels acquire a sizable width at large λ = κ, which may actually exceed the full width of THDM type. In fact, as could be anticipated, the magnitude of the bosonic three-body widths is comparable (with opposite sign) to that of the 1L effects in the dominant (SU (2) L -conserving) bosonic two-body widths, amounting to 5-10% for this specific scenario. Given that the 1L corrections lead to a larger suppression of the two-body width in the pseudoscalar case, which we mentioned earlier, the three-body widths are also more important for A, as compared to H.
This balance between virtual corrections and boson radiation underlines the importance of the IR effects in the three-body channels. Indeed, contrarily to the THDM three-body modes encountered in Sect. 2.2, the widths involving a S are not characterized by the quartic couplings, but by the 'dressing' of the leading two-body channels H → Za S , A → ha S and H ± → W ± a S with radiated light bosons. The latter evidently include electroweak gauge bosons and explain the large contributions to e. g. H, A → W + W − a S , ZZa S . From the SU (2) L perspective, one of the gauge bosons then replaces the associated Goldstone boson while the other behaves vector-like. The width H → Zha S could be understood in similar terms-at least for moderate values of λ = κ. However, for the larger range of λ = κ, this channel shares some characteristics with A → a S hh, which obviously violates the SU (2) L symmetry, but emerges as an important channel. The effect at stake here is the radiation of h from the a S line in the soft/collinear regime. Kinematical integration in the off-shell regime close to the a S pole indeed generates a scaling ∝ m −2 a S of the decay width, which compensates the suppression by the SU (2) L -breaking trilinear Higgs coupling a S -a S -h: in the aftermath, Γ[A → a S hh] ∝ M v 2 /m 2 a S , where M ∼ µ eff , M H ± is one of the large masses of the problem. Therefore, the importance of such SU (2) L -violating channels emerges as a consequence from the accidental fact that m 2 a S ∼ v 2 . Finally, we stress that, from the IR nature of these effects, they largely compensate between virtual and real corrections, so that SU (2) L remains a pertinent approximate symmetry for the inclusive width.
The precision on these (off-resonance) three-body decay widths calculated at the tree level cannot be very competitive. Given the dominance of IR effects, we can expect the absolute uncertainty to match the residual one in virtual corrections to the two-body decays, i. e. amount to percent level of the two-body widths, or 20-50% of the three-body widths. Yet, the general order of magnitude should be under control.

Branching ratios
With SUSY decays kinematically inaccessible in the considered scenario, we have completed the calculation of all the relevant widths for a 1L evaluation of the total widths and branching ratios of the heavy doublet-Higgs states. The corresponding results for the scenario of Fig. 4 are displayed in Fig. 7.
Focusing on the total width first, we stress the relevance of the three-body channels for a consistent evaluation: the difference between the full widths and the pseudo-widths, considering only two-body channels at 1L, is indeed comparable to the magnitude of the 1L effects, implying an error of full 1L magnitude when the three-body modes are overlooked. In addition, we observe that the magnitude of the SU (2) L -violating effects, measured as the difference between the widths for the scalar and pseudoscalar states (see plot in the top right-hand corner of Fig. 7), remains comparable at the tree level and at 1L, but about doubles when considering the pseudo-widths. This feature is the natural consequence of SU (2) L -violating effects of IR type developing in virtual corrections and Higgs radiation, but compensating between the two.
The general magnitude of the widths considerably varies with the choice of λ = κ, as NMSSMspecific effects range from negligible to dominant. This is most visible at the level of the branching ratios where the THDM fermionic modes (red solid curves) decrease from a significance of about 100% at low λ to percent level at λ ≈ 0.2. While the absolute magnitude of the corresponding widths (about 5 GeV in this scenario) remains constant, the competition of the bosonic Higgs decays involving a S (blue solid line) steadily increases with λ, resulting in a completely different phenomenology. The general picture qualitatively changes little between the tree level and the 1L order. Only the emergence of the three-body decays diminishes the significance of the two-body modes, substituting them with more complicated final states for a weight of ∼ 5-10%. Both fermionic and bosonic three-body modes involving a S become important at large λ and may compete in their own right with the THDM modes. As in the MSSM case, however, bosonic final states of THDM type-i. e. involving gauge bosons and the SM-like Higgs boson-remain subdominant in the whole considered range of parameters.
As a concluding note, let us stress that the relative weight of the three-body modes in the decays of heavy Higgs bosons and their mostly IR origin would make it necessary to closely scrutinize the experimental cuts in order to quantitatively exploit a hypothetical measurement of a new Higgs boson and identify the relative importance of the various channels. 1L order

Conclusions
In this paper, we have shown how the approximate SU (2) L symmetry in the (experimentallyfavored) limit M H ± M Z severely constrains the decays of the heavy doublet-dominated Higgs bosons of the (N)MSSM. In particular, most available MSSM-like bosonic channels involving two-body final states are purely radiative or SU (2) L -violating, meaning that an evaluation at 1L is only able to estimate the leading order effects, with a very large uncertainty-due to e. g. QCD corrections to fermion loops-persisting in the absolute determination of the widths and branching ratios. The fermionic decay modes thus emerge as the leading contributors to the widths of the heavy-doublet scalars of the MSSM. The bosonic three-body decays are not necessarily SU (2) L -violating and prove competitive with the two-body ones as soon as M H ± reaches a few TeV. Finally, comparisons with FeynHiggs justify the necessity to modernize this tool in order to study the phenomenology of heavy Higgs bosons.
The situation can drastically change in the presence of comparatively light singlet-dominated states, provided no singlet-doublet decoupling is enforced. Then, two-body decay modes involving Higgs-to-Higgs transitions exist that do not break the electroweak symmetry and may even dominate the decays of the heavy mostly-doublet states. Since no symmetry suppresses the tree-level contributions, such transitions are better controlled by the perturbative expansion (in constrast to their MSSM-like counterparts), with uncertainties at 1L order falling to the percent level. Of course, SU (2) L -violating modes with a large remaining uncertainty at the full 1L order are also present; the exact precision will depend on the details of the scenario. We stressed the importance of IR effects in the 1L corrections (reaching O(10%) in the considered example), underlining also that these (virtual) corrections do not necessarily preserve the SU (2) L symmetry independently from the corresponding tree-level three-body decays (radiation). The latter contribute at the same order as 1L two-body decays and emerge as an essential ingredient for a reliable determination of the widths and branching ratios at the complete 1L order.
With the completion of the evaluation of Higgs decays into SM final states at 1L order, it is possible for us to consistently predict the total widths and branching ratios of Higgs bosons in scenarios where the SUSY spectrum decouples. In the converse case, further competition can be expected from heavy Higgs decays into e. g. electroweakinos, sleptons or even squarks. As such, our results can also be generalized to other models based on a THDM structure.

A. Bosonic decays of a heavy charged Higgs
As dictated by the SU (2) L symmetry, the decay channels of the charged Higgs in the high-mass regime show little novelty as compared to the case of the neutral heavy-doublet states. From the perspective of electromagnetism, a W boson is always present among the bosonic final states (as this is the only lighter fundamental charged boson).
In the MSSM, the only kinematically accessible bosonic two-body decay at the tree level is H ± → hW ± . The channels H ± → ZW ± and H ± → γW ± are purely radiative. In fact, the decay H ± → hW ± is SU (2) L -violating, so that it would also vanish at the tree level in the exact alignment limit. Thus, similarly to the channel A → Zh, the associated width scales like Γ[H ± → hW ± ] ∝ M 2 EW M H ± , with the M 2 EW M 2 H ± suppression emerging from the mixing in the CP-even sector. As we noted in Sect. 2, a considerable uncertainty persists at the 1L order on such tree-level suppressed channels: this is illustrated by the dispersion of the corresponding predictions in the left plot of Fig. 8, where the choice of scheme for the Yukawa couplings or the inclusion of a 1L 2 term result in variations of the order of 100%. The non-resonant three-body decays in the limit M H ± M EW are determined by the quartic Higgs couplings and do not necessarily break the SU (2) L symmetry. Consequently, they are competitive with the two-body bosonic widths as soon as masses in the few TeV range are considered. The relevant widths are shown in the plot on the right-hand side of Fig. 8.
In the NMSSM, the accessibility of comparatively light singlet scalar states modifies the situation with respect to the SU (2) L symmetry, as unsuppressed channels correspondingly emerge. The latter can completely dominate the charged Higgs width if the singlet-doublet interactions are sufficiently intense. In this case, the tree-level contribution typically determines the widths and a percent-level precision can be expected at the 1L order. The width of the leading twobody channel H ± → W ± a S in the scenario of Sect. 3 is shown on the left-hand side of Fig. 9, with 1L corrections of order 10%. Once again, these 1L contributions (for M H ± M EW,a S ) involve sizable effects of IR-type, so that (non-resonant) three-body widths of comparable magnitude (∼ 10% of the two-body width) balance them at the level of the full width. These three-body channels involving a S in the scenario of Sect. 3 are shown on the right-hand side of Fig. 9. In this case, at large λ = κ, they may compete in magnitude with the MSSM-like channels.

B. Input parameters
The input parameters for the scenarios in Sect. 2 and Sect. 3 are summarized in Tab. 1. The SUSY scale m SUSY sets the bilinear SUSY-breaking parameters in the sfermion and electroweakino sectors; the gluino-mass parameter is given by mg. The tree-level value of A κ is fixed according to the on-shell renormalization scheme for the pseudoscalar singlet a s described in Ref. [35]; at λ = κ = 0 it is set to −3 GeV such that M as ∼ m as ∼ 300 GeV. At that point, it also agrees with the input value of A κ in the DR scheme at the renormalization scale µ ren = m t = 172.3 GeV.