Heavy to light Higgs boson decays at NLO in the Singlet Extension of the Standard Model

We study the decay of a heavy Higgs boson into a light Higgs pair at one loop in the singlet extension of the Standard Model. To this purpose, we construct several renormalization schemes for the extended Higgs sector of the model. We apply these schemes to calculate the heavy-to-light Higgs decay width at next-to-leading order electroweak accuracy, and demonstrate that certain prescriptions lead to gauge-dependent results. We comprehensively examine how the NLO predictions depend on the relevant singlet model parameters, with emphasis on the trademark behavior of the quantum effects, and how these change under different renormalization schemes and a variable renormalization scale. Once all present constraints on the model are included, we find mild NLO corrections, typically of few percent, and with small theoretical uncertainties.


Introduction
With the discovery of the Higgs boson [1,2], the LHC has finally reached the very frontiers of the Standard Model (SM). Dedicated analyses based on Run I data have so far shown excellent agreement between the observed 125 GeV bosonic resonance and the scalar particle originally postulated by Higgs [3,4], Englert and Brout [5], and Guralnik, Hagen, and Kibble [6]. Notwithstanding, it is an ongoing task to decipher whether such a state corresponds indeed to the SM agent of electroweak (EW) symmetry breaking [7][8][9], or if alternatively the LHC has unveiled just one Higgs-like state among many others, a composite state, or the overlap of multiple resonances, just to mention few possibilities. Moreover, the state-of-the-art precision in the Higgs coupling extraction lies within the 10 − 20% level [10], right at the ballpark of the deviations predicted by popular new physics models. The overall picture strengthens the belief that perhaps the Higgs discovery is in fact our first glimpse at a more fundamental UV complete structure.
The arguably most simple, renormalizable extension of the Higgs sector, is constructed by expanding the SM Lagrangian with one additional spinless real electroweak singlet [11][12][13]. This adds up one extra Higgs companion to the physical spectrum of the model, providing an excellent framework to guide collider searches for exotic scalars, either via direct production or through off-shell effects. Moreover, the coupling between the doublet and the singlet fields mixes the two neutral Higgs states, leading to rescaled Higgs couplings to the SM particles. A chief prediction for collider phenomenology are the universally suppressed cross sections and partial amplitudes in all Higgs production modes and decay channels.
Another paramount signature of a second Higgs is the possibility of the heavy-to-light Higgs decay mode H → hh [12,14]. The process is governed by the Higgs self-coupling λ Hhh , and as such it constitutes a direct probe of the scalar potential of the model. If the new Higgs boson is lighter than the SM Higgs, the novel H → hh decay distorts all Higgs branching ratios, and thereby its signal strengths, from the SM expectations. Alternatively, if the extra scalar is identified with a heavier Higgs companion and m H > 2m h , H → hh can significantly contribute to the heavy Higgs width and lineshape and modify its decay pattern.
In this paper we present a thorough study of the heavy-to-light Higgs decay mode H → hh in the singlet extension of the SM, including the complete set of radiative corrections at one loop. Aside from being relevant on its own, we use this process as a physics case to construct and compare different renormalization schemes for the extended Higgs sector of the model. Using the Sloops [15][16][17][18] general non-linear gauge fixing setup, we illustrate how certain prescriptions still exhibit gauge dependence for physical quantities. Our task carries to completion the steps initiated in our previous publication [19], and sharpen up all theoretical tools necessary to completely characterize the singlet model phenomenology at next-to-leading order (NLO) accuracy.
The remainder of the paper is organized as follows. In Section 2 we provide a brief description of the model setup and constraints. In Section 3 we discuss in full detail our renormalization setup. We devote Section 4 to characterize the general aspects of heavy-to-light Higgs decays at leading order and at one loop, while in Section 5 we present a detailed phenomenological analysis. We summarize and conclude in Section 6. Additional analytical details are provided in the Appendix.

Model setup at the classical level
We construct the singlet extension of the SM by adding one colorless, real scalar field, which transforms as a singlet under the SU (2) L ⊗ U (1) Y gauge charges [11][12][13][19][20][21]. Such a simple renormalizable extension can be viewed as a simplified model approach to the low-energy Higgs sector of a more fundamental UV completion, for instance the decoupling limit of the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [22], some realizations of GUTs [23], models with additional gauge sectors [24] or hidden valleys [25]. The implications of this model were addressed for the first time in Refs. [11][12][13], and it has been the object of dedicated investigation for the past two decades, displaying a very attractive phenomenology, especially in the context of collider physics, see e.g. . It has also been subject to many dedicated searches by the LHC experiments, cf. e.g. [59][60][61][62][63][64] for recent studies.

Classical Lagrangian
The singlet scalar extension of the SM (denoted as xSM ) is defined by the Lagrangian L xSM = L gauge + L fermions + L Yukawa + L scalar + L GF + L ghost (1) where the gauge boson and fermionic kinetic parts L gauge,fermions , as well as the Yukawa interaction L Yukawa , are given by the respective SM contributions. The gauge-fixing and ghost Lagrangians L GF,ghost will be defined in more detail below. The scalar sector is given by where D µ is the covariant derivative and V(Φ, S) the scalar potential The latter corresponds to the most general, SU L (2) ⊗ U (1) Y -invariant, renormalizable Lagrangian involving the Higgs doublet Φ and the singlet S fields, and compatible with an additional discrete Z 2 symmetry, that precludes other terms with odd (e.g. cubic) field powers in the potential. By assuming all parameters in Eq. (3) to be real, we disregard additional sources of CP violation.

Mass spectrum
The doublet and singlet fields are expanded as where v ≡ √ 2 Φ = ( √ 2 G F ) −1/2 246 GeV and v s ≡ √ 2 S stand for their respective vacuum expectation values (vevs). The fields G ± , G 0 denote the charged and neutral Goldstone bosons. Since the singlet transforms trivially under the electroweak gauge group, only the doublet vev v takes part in electroweak symmetry breaking, which proceeds exactly as in the SM case.
The linear terms in the fields φ h and φ s from Eq. (3) lead to the tadpole relations by which we can express the minimization condition of the Higgs potential (3) as T φ h ,φs = 0.
In turn, the quadratic terms in the Higgs fields may be arranged into a 2 × 2 squared mass matrix M 2 hs . In the gauge basis (φ h , φ s ) these take the form where the tadpole component T φ h ,φs and the squared mass matrix M 2 φ h ,φs are defined by Requiring this matrix to be positively defined leads to the stability conditions Once we impose the tadpoles T φ h,s to vanish, Eq. (8) can be readily transformed into the (tree-level) Higgs mass basis through the rotation U (α) · M 2 φ h ,φs · U −1 (α) = M 2 hH = diag(m 2 h , m 2 H ), the physical Higgs masses reading We identify the corresponding mass-eigenstates as a light [h] and a heavy [H] CP-even Higgs boson, which are related to the gauge eigenstates through where the rotation angle α is defined in the range −π/2 ≤ α ≤ π/2 by sin 2α = λ 3 vv s Likewise the tadpoles in the gauge basis [T φ h , T φs ] may be rotated into the mass basis [T h ,T H ] through U (α): The above equations are of service to recast the quartic couplings in the Higgs potential (3) The mixed mass term m 2 hH denotes the (symmetric) off-diagonal element of the squared mass matrix M 2 hH , defined in the mass-eigenstate basis. While at tree-level we have T h,H = 0 and m 2 hH = 0, keeping these dependencies explicit in Eqs. (14)-(16) will be useful to link the Lagrangian parameter counterterms to the corresponding mass counterterms, as we discuss in Section 3. Similarly, it is practical to rephrase the bilinear mass terms µ and µ s in Eq.(3) as

Input parameters
The Higgs sector of the model is determined at tree-level by i) the doublet vev, bilinear mass term and quartic self-coupling; ii) their counterparts for the singlet field; and iii) the portal coupling λ 3 between both. The singlet vev is traded as customary via the parameter tan β ≡ vs v . a With the help of the above relations (5)-(18), we can conveniently recast them in terms of the following five independent parameters: Two of them are readily fixed in terms of experimental data: the doublet vev is linked to the Fermi constant through v 2 = ( √ 2 G F ) −1 , while one of the Higgs masses is given by the LHC value 125.09 GeV [65]. Overall, we are left with three quantities which define the relevant directions in the new physics parameter space. This is also helpful to identify which physical parameters can be used to fix the renormalization conditions.

Gauge-Fixing Lagrangian
Gauge invariance will play an important role when discussing the renormalization of the singlet Higgs sector, in particular in defining a gauge-independent mixed mass counterterm, as we discuss in detail in Sections 3.3.6 and 5.2. Such non-linear gauges have proven useful to check the gauge independence of higher order calculations within the SM [66][67][68], and its supersymmetric extensions [15][16][17][18][69][70][71][72][73]. The gauge-fixing Lagrangian can be written in general as where the functions F depend non-linearly on the Higgs and gauge fields, In the above equations e is the electromagnetic coupling constant, g the SU (2) L coupling constant and θ W denotes the weak mixing angle. The quantities {α,β · · ·˜ 2 } correspond to the generalized gauge-fixing parameters. Setting these parameters to zero leads to the standard linear R ξ gauge fixing, with ξ i = 1 defining the familiar 't Hooft-Feynman gauge. In our renormalization setup we will take these gaugefixing terms already as renormalized quantities -in such a way that no additional counterterms δL GF are introduced for this part of the Lagrangian.
In turn, the ghost Lagrangian L ghost is derived by requiring the complete Lagrangian at the quantum level to be invariant under BRST transformations. This means δ BRST L xSM = 0 and hence δ BRST L GF = −δ BRST L ghost . This follows from the fact that by construction the gauge, fermionic, Yukawa and scalar a Note the different conventions for tan β employed in the literature. The definition we adopt herewith is the inverse of that from Refs. [19][20][21].
components of L xSM are invariant under BRST transformations, as the latter are equivalent to gauge transformations. We consider both L GF and L ghost to be written in terms of renormalized quantities. The BRST transformations specific to the singlet extension are given by where c θ , s θ are shorthand notations for cos θ, sin θ, while c Z , c ± , c A stand for the Faddeev-Popov ghost fields associated to the Z 0 , W ± and the photon field respectively. Within this particular gauge fixing, we set in practice ξ W,Z,A = 1. This is convenient since the gauge boson propagators take a very simple form, while we still keep the possibility to check the gauge independence of the final result, at the expense of adding new gauge parameter-dependent vertices to the model [66]. The gauge independence of a calculation can be examined numerically by varying the non-linear gauge parametersα,β · · ·˜ 2 . Technically, we perform our implementation of the singlet model with non-linear gauge fixing using LanHEP [74,75] and Sloops [15][16][17][18].

Interactions
The key theoretical structure in the model is the doublet-singlet portal coupling L xSM ⊃ λ 3 (Φ † Φ) S 2 , which is responsible for the Higgs mass eigenstates to be admixtures of the doublet φ h and the singlet φ s neutral components. One main consequence of this mixing is the universal depletion of all Higgs boson couplings to the SM particles as This global rescaling is ultimately due to the fact that, owing to electroweak gauge invariance, only the doublet component can couple to the fermions (via ordinary Yukawa interactions) and the gauge bosons (via the gauge covariant derivative). The limits of sin α = 0 (resp. cos α = 0) correspond to the decoupling scenarios for the light (resp. heavy) Higgs bosons, in which all couplings of the additional scalar to SM fields identically vanish. The Higgs self-interactions do not obey such a plain rescaling pattern. Instead, they depend non-trivially on the cross-talk between the singlet and the doublet fields. Analytic expressions for their Feynman rules can be found in the Appendix.

Constraints
As discussed above, the singlet model specified by the Lagrangian (1) is subject to numerous constraints, which have been explicitly discussed in [19][20][21]. Although our primary focus in this paper is the structure of the higher order corrections in this model irrespectively of the parameter constraints, we briefly remind the reader which ranges are still feasible from the theoretical and experimental sides -and include all of them in our phenomenological analysis of Section 5. We consider
• perturbativity of the self-couplings in the scalar potential, i.e. |λ i (µ run ) | ≤ 4 π, where the couplings are evolved through standard renormalization group equations [79] and evaluated at a reference high-energy scale µ run ∼ 4 × 10 10 GeV c for the high-mass and at µ run = v = 246 GeV for the low-mass scenario (see [20,21] for a more detailed discussion). The high-mass (resp. low-mass) regions correspond to m H > 2m h , with m h = 125.09 GeV (resp. m H > 2m h , with m H = 125.09 GeV); • vacuum stability (cf. Eqn. (9)) up to the same high-energy scale.
It is interesting to observe the interplay of these different constraints on the overall parameter space. We here only summarize the main features e -a dedicated discussion can be found in [21].
• in the high mass region the leading constraints stem from i) direct searches (for m H 300 GeV) f ; ii) the difference between the experimental W-mass measurement and its theoretical prediction [19] (in the intermediate range M H ∈ [300 GeV; 800 GeV]); and iii) perturbativity of the self-couplings in the scalar potential (for m H ≥ 800 GeV). All these features are summarized in Figure 1 and Table 1; b Electroweak gauge bosons are replaced by Goldstone scalars according to the Goldstone equivalence theorem [78]. c This scale has been chosen such that the model still guarantees vacuum stability at a scale slightly larger than the SM breakdown scale for which λ1 ≤ 0 in the SM limit (sin α = 0). Requiring validity up to higher scales leads to stronger constraints, cf. the discussion in [20]. d A detailed discussion of the determination of limits from the Higgs signal strength can be found in [21]. For mH ≤ 152 GeV, we test a two-scalar hypothesis versus the LHC data, leading to an mH -dependence for the respective χ 2 . Results in [89], however, are derived under an SM assumption. For mH ≥ 152 GeV, the χ 2 is independent of the second resonance mass and in this range we therefore adopt the improved combined experimental limit. e See also [90]. f Note that the most recent experimental searches published in 2015 have not been included. These potentially influence the allowed regions for mH 300 GeV. Indeed, preliminary studies show that results from [62] especially modify constraints for mH ≤ 250 GeV [91].   [19]; ii) electroweak precision observables (EWPOs) tested via the oblique parameters S, T and U (orange, dashed); iii) perturbativity, of the RG-evolved coupling λ 1 (blue, dotted), evaluated for an exemplary choice tan β = 10, iv) perturbative unitarity (grey, dash-dotted), v) direct LHC searches (green, dashed), and vi) Higgs signal strength measurement (magenta, dash-dotted). For Higgs masses m H ∈ [300 GeV; 800 GeV] the W boson mass measurement yields the strongest constraint [19]. The present plot corresponds to an update of figure 8 from [21], where the latest constraints from the combined signal strength [89] have been taken into account.
• in the low mass region where m H ∼ 125 GeV, the parameter space is extremely constrained, especially from demanding agreement with the LHC Higgs signal strength measurement and the LEP constraints. In Table 2 we summarize these constraints. Note that in this regime the SM limit corresponds to | sin α| = 1. Tables 1 and 2 show the current constraints for the maximal (minimal) allowed values of sin α and tan β, following the analysis presented in [21]. Note that the minimal tan β values shown here were taken at a fixed value of sin α, so results from a generic scan might slightly differ. All the constraints mentioned above have been taken into account when considering viable parameter space regions of the model for our numerical analysis in Section 5. Also the results from the combined ATLAS and CMS signal strength fit have been included when applicable. We expect the results from the most recent LHC searches to influence the global picture in the mass region m H 350 − 400 GeV, while for higher values the W boson mass still poses the strongest constraint on the mixing angle.

Setup
The renormalization program we present here sticks close to the general strategy followed in multidoublet Higgs extensions such as the MSSM [17,93] Table 1: Table II from [21], with adjusted conventions for tan β, and updated constraints on the maximally allowed mixing angle from the combined Higgs signal strength fit [89]. It presents allowed ranges for sin α and tan β in the high mass region for fixed Higgs masses m. The allowed interval of sin α was determined fixing (tan β) −1 = 0.15. The 95% C.L. limits on sin α from the Higgs signal rates are derived from one-dimensional fits and taken at ∆χ 2 = 4. The lower limit on sin α always stems from vacuum stability, and the upper limit on tan β always from perturbativity of λ 2 , evaluated at sin α = 0.1. The source of the most stringent upper limit on sin α is named in the third column. We fixed m h = 125.1 GeV and the stability and perturbativity were tested at the reference scale µrun ∼ 4 × 10 10 GeV.
the required counterterms by introducing multiplicative renormalization constants to the weak coupling constant, fields, tadpoles, masses and vevs. These are then fixed by as many renormalization conditions as independent parameters are present in the model [95]. We adopt on-shell conditions to renormalize the electroweak gauge parameters [18,[96][97][98][99] and the diagonal terms of the Higgs boson mass matrices. Using an on-shell scheme, as customary in this context, provides an unambiguous interpretation of the bare parameters in the classical Lagrangian in terms of physically measurable quantities. We also recall that field renormalization constants are not needed if we only require the observables derived from S-matrices to be finite, but not each of the Greens' functions individually. They are nonetheless convenient from the technical viewpoint, as they account for loop corrections to the external legs and less Feynman diagrams have to be explicitly included. We proceed as customary by splitting the bare Lagrangian of the model (1) into the renormalized and the counterterm pieces as L 0 ({X 0 }) → L ({X}) + δL ({δX}). Accordingly, we rewrite each of the bare parameters X 0 as a renormalized part X and its counterterm δX. For the purpose of this work we only need to deal with the scalar and gauge sectors L 0 scalar,gauge , as the other sectors do not feature for the remainder of our discussion. We also recall that the gauge-fixing Lagrangian L GF does not contribute to δL , since we choose to write it already in terms of renormalized fields and parameters [17,66]. The physical parameters of the gauge sector are the electromagnetic coupling constant e and the gauge boson masses m W , m Z that we split as [99], Also the bare parameters appearing in L 0 scalar in the gauge-eigenstate basis are decomposed as  Table 2: Table III from [21], with adjusted definition for tan β and updated constraints on the minimally allowed mixing angle from the combined Higgs signal strength fit [89]. It presents limits on sin α and tan β in the low mass scenario for various light Higgs masses m h . The limits on sin α have been determined at tan β = 1. The lower limit on sin α stemming from exclusion limits from LEP or LHC Higgs searches is obtained using HiggsBounds [85][86][87]92] and given in the second column. If the lower limit on sin α obtained from the test against the Higgs signal rates using HiggsSignals [88] results in stricter limits, we display them in the third column. The upper limit on tan β in the fourth column stems from perturbative unitarity for the complete decoupling case (| sin α| = 1). In the fifth column we give the tan β value for which Γ H→hh = 0 is obtained, given the maximal mixing angle allowed by the Higgs exclusion limits (second column). At this tan β value, the | sin α| limit obtained from the Higgs signal rates (third column) is abrogated.
A similar splitting is introduced for the Higgs tadpoles T 0 , which feature explicitly for calculations beyond the leading order. Equivalent expressions can be written trading some of the above bare parameters for more physical ones through the relations given by Eqs. (14)- (18).
In our setup we choose not to renormalize the mixing angle α. Instead, we promote the relation between the Higgs eigenstates in the gauge (φ h , φ s ) and mass basis (h, H) to be valid to all orders, Doing so, the bare and the physical mixing angle coincide and we need no additional mixing angle counterterm.
In turn, field renormalization constants for the physical Higgs states are introduced by shifting the bare Higgs fields in the mass basis as Finally, we introduce the (matrix-valued) Higgs mass counterterm via where the generic index φ applies to the squared mass matrix in both the gauge and the mass basis. Their respective matrix counterterms are linked through where the mixed mass counterterms are symmetric δm 2 Hh = δm 2 hH . Thus, aside from the purely SM inputs, the renormalization of the scalar sector in the singlet model is completely specified by four renormalization constants for the neutral Higgs fields, the respective singlet and doublet vev counterterms, and five additional counterterms linked to the parameters in the Higgs potential (3). In the mass eigenstate basis, these can be traded by: Defining a renormalization scheme is then tantamount to identifying a set of independent conditions by which to link the above quantities to physical inputs. The renormalization conditions by which we fix these counterterms will rely on the one-point and two-point Greens' functions of physical fields. Depending on the scheme we choose for Higgs field renormalization, not all of the above field renormalization constants will be independent from each other. A complete renormalization scheme fixes all the counterterms which are necessary to absorb the UVdivergent contributions from loop-level amplitudes, such that one obtains UV finite predictions for physical observables. Another important property of a renormalization scheme is gauge independence. More precisely, maintaining gauge independence when defining a scheme allows to write physical predictions as a function of the input parameters in a way that does not vary when the gauge-fixing is changed. Only in this case one can unambiguously relate physical observables to Lagrangian parameters. In this work we examine different strategies to extend the conventional SM renormalization to the singlet model case, and discuss in detail whether these comply with gauge independence. The Sloops non-linear gauge-fixing setup (cf. Eq. (20)) turns out to be instrumental in this task.

Gauge sector
We begin by introducing the on-shell definition of the electroweak mixing angle sin 2 . This relation fixes the weak mixing angle counterterm (s 2 The weak gauge boson masses are renormalized in the standard on-shell scheme [18,[96][97][98][99], i.e. by requiring the real part of the transverse renormalized weak gauge boson self-energies to vanish at the respective pole masses. The condition where δZ V stands for the weak gauge boson field renormalization V → Z 1/2 All renormalized self-energies are denoted hereafter by a hat. The transverse part of the gauge boson self-energies follows from the vacuum polarization tensor, The explicit form of the weak gauge boson two-point functions in the singlet model can be found in Ref. [19].
To renormalize the electromagnetic coupling constant, we require the electric charge to be equal to the full eeγ vertex in the Thompson limit. With the help of the QED Ward identities, this condition is given in terms of the photon and mixed Z − γ two-point functions To avoid large logarithms from light fermion masses, we rephrase as customary the photon vacuum polarization as where the superindex indicates that only the light fermion contributions (all leptons and quarks, except the top) are included in the photon self-energy, while the QED-induced shift to the fine structure constant, is known to very good accuracy [100,101]. The improved electric charge counterterm in the Thompson limit is thus given by in such a way that the value of the renormalized electric charge at zero momentum transfer e(0) = 4π α em (0) can be extracted from the measured fine-structure constant in this limit: α em (0) = 1/137.035999074(44) [102].
On the other hand, the very precise measurement of the muon lifetime provides a link between the weak gauge boson masses, the fine structure constant and the Fermi constant. This allows for different input choices to fix the electroweak sector. In our numerical analysis we shall use two alternative parametrizations: • The α em -parametrization, in which we select α em (0) and m W,Z as input parameters; • The G F -parametrization, in which we instead replace the W-boson mass by the Fermi constant G F = 1.1663787(6)×10 −5 GeV −2 [102], the latter being fixed by the muon lifetime via [98,[103][104][105].
These two parametrizations are related via the conventional parameter ∆r [98,[103][104][105][106][107] as where m W,Z and s W are renormalized in the on-shell scheme. For a detailed analysis of ∆r in the singlet model cf. Ref [19]. Since ∆r vanishes at leading-order, both parametrizations are trivially linked at tree-level as , while they depart from each other at higher perturbative orders. We will explicitly quantify these departures further down in Section 5.

Tadpole renormalization
For the tadpole renormalization we proceed as customary [18,[96][97][98][99] and imposê This is equivalent to requiring that v and v s are the physical vacuum expectation values of the doublet and the singlet fields respectively, so that they define the (renormalized) minimum of the Higgs potential.
In practice, this implies that no Higgs one-point insertions feature explicitly in our calculation.

Doublet vev renormalization
The vev v of the scalar doublet Φ is fixed as in the SM through its relation to the electroweak on-shell parameters where all needed counterterms are defined in (37), (42) and (35).

Singlet vev renormalization
The general renormalization transformation of a generic scalar field vev [108] can be particularized to the singlet vev case as where we have introduced for convenience the singlet field renormalization in the gauge basis . The additional counterterm δv s characterizes to what extent the singlet vev renormalizes differently from the singlet field φ s itself. In Ref. [108] it was shown that, in an R ξ gauge, a divergent part for δv s is forbidden if the scalar field obeys a rigid invariance (see also Ref. [109] and references therein). This is precisely the case in the singlet model, since the singlet field is unlinked from the gauge sector and hence invariant under global gauge transformations. In addition to that, the singlet field renormalization constant δZ S is also UV-finite. This can be easily shown by computing δZ S in the unbroken phase where Φ = S = 0. Such a scenario is analogous to a plain λφ 4 -theory, in which the (singlet) scalar field is coupled to a second scalar (doublet) field only through the gauge-singlet quartic In this case, all one-loop contributions to the singlet two-point function are momentum-independent, implying that δZ S does not get a UV pole (cf. e.g. [110]). g We thus conclude that the singlet vev counterterm δv s = δv s + δZ S /2 defined by Eqs. (29) and (46) gets at most a finite contribution at this order.
We finally note that, given the condition of vanishing tadpoles (44), the singlet vev v s corresponds to the physical minimum of the Higgs potential in the broken phase in the singlet field direction, viz. S = v s / √ 2, at a given order in perturbation theory. Since v s does not contribute to the electroweak symmetry breaking, it cannot be fixed in terms of SM observables. Instead, we must promote it to an independent input parameter (which should eventually be determined from a future measurement of e.g. the Hhh coupling). So doing, any finite shift δv s can be subsumed into the physical definition of v s itself at one-loop. Therefore, in our renormalization setup we can simply fix δv s in the MS scheme, δv MS s = 0, so that no singlet vev counterterm features in one-loop calculations.

Higgs masses renormalization
The (matrix-valued) Higgs mass counterterm in the gauge basis yields where we have already fixed δv s = 0, as justified above. This result can be linked as customary to the mass basis through Eq. (33).
To renormalize the physical Higgs masses we impose on-shell conditions on the renormalized diagonal Higgs self-energies, whereby we obtain The explicit form of the field renormalization constants δZ φ in different schemes is discussed below in Section 3.3.5.
In theories where the gauge eigenstates mix, the renormalization of the non-diagonal or mixing terms must be addressed with care (cf. [17,18,93] for an analogue discussion in the context of the squark sector in the MSSM). As we have seen in Section 2, a bare angle α 0 ≡ α rotates the scalar fields from the gauge basis to the mass basis through Eq. (11). While such diagonal form is valid at leading order, radiative corrections will in general misalign the (tree-level) mass eigenstates. This is reflected in the off-diagonal terms of the loop-corrected propagators, g We have numerically verified that δZS is UV-finite at one loop in all of the different renormalization schemes.
traded by the non-diagonal Higgs two-point function One possibility is to absorb these additional quantum effects into the renormalization of the mixing angle. This is equivalent to diagonalizing the loop-corrected mass matrix further through an additional rotation U (δα), where δα plays the role of a mixing angle counterterm, such that α 0 → α + δα. Alternatively, in our approach we take the mixing matrix U (α) as written in terms of the physical mixing angle and hence valid to all orders. The two alternative approaches are related through The residual mixing induced by the off-diagonal terms in the mass matrix is instead removed by the non-diagonal field renormalization constants, which we present below.

Higgs field renormalization: diagonal parts
Taking the Higgs boson masses m h,H as experimental inputs, we fix the diagonal field renormalization constants via the on-shell conditions where ReΣ φ (p 2 ) was defined in Eq. (48), while the familiar shorthand notation f (p 2 ) ≡ df (p 2 )/dp 2 denotes the derivative with the respect to the momentum squared. This leads to which set the Higgs propagator residues to unity in the limit p 2 → m 2 φ (φ = h, H).

Higgs field renormalization: non-diagonal parts
Fixing the non-diagonal field renormalization is a crucial step in setting up a gauge-invariant scheme, in which the renormalized one-loop amplitudes are independent of the gauge-fixing parameters, as discussed above. We first construct a set of schemes in analogy to the more familiar approaches in the literature. As we will show, these lead in general to gauge-dependent predictions for physical observables. To circumvent this problem, we introduce an additional (dubbed improved ) scheme, which is defined merely in terms of two-point functions and gives numerically stable results throughout the entire parameter space. Similar discussions are addressed e.g. when defining renormalization schemes for the parameter tan β in the MSSM [17,111].
• Minimal field: As a first setup to fix the non-diagonal Higgs field renormalization δZ hH we resort to a minimal field renormalization. We attach one single renormalization factor per field in the gauge basis, where we have expanded them to first order.
This procedure is in straight analogy to the conventional renormalization of the Higgs sector in multidoublet extensions such as the MSSM [93] and the Two-Higgs-Doublet Model [94]. Assuming symmetric off-diagonal components, and using the rotation matrix U (α) in Eq. (11), we can write the physical Higgs wave function renormalization constants in terms of the gauge basis ones δZ Φ,S as with the shorthand notation {s α , c α , t α } = {sin α, cos α, tan α}. The scheme is dubbed minimal as the non-diagonal field renormalization δZ hH is not independent. Instead, it is linked to the diagonal parts δZ h,H , which we have already fixed via on-shell conditions (53). Additionally, since at one-loop we have δZ MS S = 0 (cf. Section 3.3.3), we can further simplify the relations above to get Finally, for the mixed mass counterterm, which enters explicitly in the Hhh vertex counterterm, we demand the off-diagonal renormalized Higgs self-energy to vanish at an arbitrary renormalization scale, From Eq. (57) we see that in this scheme all vertices with external Higgs legs receive a finite wave-function renormalization correction, which absorbs the residual loop-induced h − H mixing for an external on-shell Higgs state. These finite wave-function factors are given in general by [93] where O(α 2 ew ) denote the contributions beyond one-loop accuracy. Since the diagonal field renormalization has been fixed via on-shell conditions Eq. (53), the diagonal finite factors at one loop yieldẐ h,H = 1 and hence we do not include them explicitly.
• On-shell: We define a second prescription in close analogy to squark renormalization [112][113][114] h . This time we attach one field renormalization constant δZ h , δZ H per Higgs field directly in the masseigenstate basis (31), in which case the off-diagonal field renormalization constants δZ hH and δZ Hh are not directly related to the diagonal terms. The diagonal parts δZ h,H are again given by the on-shell relations of Eq.
h While this work was being finalized, we learned of the work [115], which presents a study of the quantum corrections to the Higgs couplings to fermions and gauge bosons in a similar singlet model setup. The renormalization scheme for the extended Higgs sector used by these authors is equivalent to the on-shell scheme we discuss here, and which, as we analyse in the following, is not gauge-independent.
Using Eq. (50) leads to Therefore, to fully fix the non-diagonal renormalization constants one must provide a proper definition of the mixed mass counterterm δm 2 hH . One possibility, as inspired from [112][113][114], is to impose δZ hH = δZ Hh , which fixes δm 2 hH accordingly as The above condition removes the loop-induced H − h mixing when either of the two Higgs bosons are on shell, so that the physical states propagate independently and do not oscillate.
The customary on-shell scheme, as well as the minimal field scheme discussed above, show indisputable benefits, e.g. the fact that all counterterms are given in terms of two-point functions and related to physically measurable quantities. However, both of them lead to renormalized one-loop amplitudes which, albeit UV finite, may still have a left-over dependence on the parameters of the gauge-fixing Lagrangian (20). This is a well known fact for on-shell fermion [116][117][118][119][120] and sfermion mixing in supersymmetric theories [18,121,122]. Exploiting the non-linear gauge fixing of Eq. (20), we explicitly verify this drawback to appear in the singlet model case as well, and illustrate it numerically in Section 5.2. In this discussion, it is worthwhile recalling that gauge dependencies may well persist in general in all non-physical building blocks which are involved in the renormalization of any gauge theory (e.g. field renormalization constants). The key test for a given renormalization scheme is thus whether it leads to gauge-independent predictions for physical observables. In the minimal field and the on-shell schemes, renormalized one-loop amplitudes are proven to contain left-over gauge-dependent contributions. These can be traced back to the mixed mass counterterm δm 2 hH , which also enters the non-diagonal field renormalization constants δZ hH,Hh . The former is fixed in these schemes through Eqs. (57) and (62) respectively, and ultimately follows from the h − H mixing self-energy. Using the non-linear gauge of (20), we find where B 0 (p 2 , m 2 , m 2 ) is the two-point Passarino-Veltman scalar integral [123] and theδ i terms are a short-hand notation for the non-linear gauge parameters in Eq. (20). The first line of Eq. (63) is identical to the result of the self-energy computation in the 't Hooft-Feynman gauge. The second and third lines correspond to the genuine gauge-dependent contributions in the non-linear gauge (we recall that for practical calculations we always set ξ A,W,Z = 1). The latter enter the mixed mass counterterm definition through Eqs. (57) or (62) and are responsible for the uncancelled dependencies on the gaugefixing parameters in the renormalized H → hh one-loop amplitude, which we pin down numerically in Section 5.2.
One first roadway to construct a gauge-independent definition of δm 2 hH alternative to Eq. (62) would be to exploit the pole structure of a process-specific one-loop amplitude (e.g. a Higgs decay) in the limit m 2 h → m 2 H , as suggested by Ref. [18]. Such a limit corresponds though to a vanishing quartic coupling λ 3 (16) and hence to a vanishing mixing angle α (12). Therefore, in this no-mixing situation, δm 2 hH cannot be defined through the mixed self-energy Σ hH , because it is identically zero. A second possibility would be to link the problematic mixed mass counterterm to a physical observable directly -viz. using a per se gaugeindependent quantity such as a decay rate or scattering cross section [17,111]. The price one would pay would be a process-dependent scheme definition, and sometimes one would have to resort to quantities out of current experimental reach. A third option is retaining only the UV-divergent part of such a quantity via an MS prescription, which we examine next. Besides this possibility, we also propose an additional prescription leading to a gauge invariant scheme, which furthermore does not render artificially enhanced contributions in any part of the parameter space.
• Mixed MS/on-shell: In this case we trade δm 2 hH by one of the Higgs self-coupling counterterms δλ i from Eq. (3), and fix it using MS conditions. For convenience we choose λ 2 and compute the divergent part of the one-loop correction to the singlet field four-point coupling λ 2 S 4 . So doing we find where ∆ stands for the UV divergent part in dimensional regularization This result is manifestly gauge independent, as it should, and agrees with the beta function for the singlet quartic coupling λ 2 given in Ref. [21]. The corresponding gauge-invariant counterterms for λ 1,3 can now be traded by δm 2 h,H , δv, δT h , δT H and δλ MS 2 using the relations from Eqs. (14)- (16), We are thus left with Finally, we use the on-shell relations (60)-(61) to obtain the non-diagonal field renormalization constants, which are now fixed in terms of Eq. (68).
Since all of the renormalization constants within δm 2 hH are either related to physical observables and/or correspond to prefactors of gauge invariant operators (e.g. λ 1,2,3 ) the mixed mass counterterm δm 2 hH is by construction gauge-invariant -and leads in turn to gauge-independent renormalized one-loop amplitudes, as we prove numerically in Section 5. This observation, together with the analytical structure of the mixed self-energy Σ hH (p 2 ) from Eq. (63), reflects that the renormalization conditions chosen for δm 2 hH (and linked to them, for δZ hH,Hh ) are the ultimate origin of the uncancelled gauge dependence found in the minimal field and the on-shell schemes.
In spite of leading to gauge-independent results, this mixed MS/on-shell scheme tends to produce overestimated radiative corrections in the phenomenologically interesting regions (s α → 0, c α → 0), as manifest from the analytic dependencies of the counterterms (66)- (67), which are proportional to inverse powers of small trigonometric factors. We therefore refrain from using this scheme explicitly in our phenomenological analysis, and instead propose an improved gauge-independent setup right below.
• Improved on-shell A second alternative to sidestep the gauge-dependent δm 2 hH definition in the default on-shell scheme is to isolate the gauge invariant part of the mixed self energy of Eq. (63). In so doing, we can use it to define the problematic mixed mass counterterm through a gauge-independent improved self-energy [124]. This is actually possible if the mixed scalar self-energy (63) is computed in the linear 't Hooft-Feynman gauge and evaluated at the average geometrical mass p 2 * = (m 2 h + m 2 H )/2. As shown in Ref. [122] with the help of the so-called pinch technique, [124][125][126], the mixed scalar selfenergy (63) obtained in this way coincides with the gauge-invariant part of the pinched result. While the results proven in Ref. [122] are applied to the squark and Higgs sectors of the MSSM, the proof does not rely on Supersymmetry and hence can be exported to the more general case of a system of two gauge eigenstates which mix in the mass basis. In addition, self-energies computed using the pinch technique are independent of the gauge-fixing scheme [124]. With these arguments in mind, we thus retain only the first line in Eq. (63) and define the mixed mass counterterm through which must be therefore gauge-independent (as we again confirm numerically in Section 5.2). Finally, the non-diagonal field renormalization are once more fixed using OS conditions (59) and fully determined in terms of δm 2 hH . In Table 3 we provide a summarized overview of the different renormalization schemes discussed in this section. Notice that they differ from each other in the renormalization conditions used to fix the non-diagonal Higgs field renormalization δZ hH,Hh constants and the mixed mass counterterm δm 2 hH . Table 3: Overview of the scheme-dependent counterterms in the different renormalization setups considered in this paper.

Heavy-to-light Higgs decay width 4.1 Leading-order contribution
When kinematically accessible, the heavy-to-light Higgs decay mode H → hh proceeds at leading order (LO) via the tree-level contact interaction λ Hhh with partial width [12,14] Γ LO H→hh = λ 2 where Notice that, owing to the structure of the scalar self-coupling, the decay width is not symmetric under a sign flip of the mixing angle s α → −s α . As such, this decay mode constitutes a genuine new physics contribution to the total heavy Higgs width -aside from the global rescaling of its decay modes into SM particles. The opening of this novel channel is thus capable to alter the Higgs boson lineshape significantly, as well as its decay pattern. More specifically, the branching fractions of the heavy Higgs boson of mass m H to SM fields φ are modified as where Γ SM H (m H ) stands for the total width of a SM-like Higgs boson with mass m H . For the lighter Higgs boson with mass m h , the branching fractions are exactly as for a SM-like Higgs with that mass. Notice that for λ Hhh = 0, all partial decay widths are universally rescaled in terms of the Higgs mixing angle α, leading to the same branching ratios that a Higgs boson of that mass would experience in the SM. The characteristic O(m −1 H ) phase-space suppression is compensated by the trilinear Higgs coupling strength λ Hhh , which depends quadratically on both the light and the heavy Higgs masses. On the other hand, there are cot β-enhanced contributions which can invigorate these Higgs self-interactions for t β < 1, and push the H → hh rates even higher. We illustrate these effects in the right panel of Figure 2, in which we show the leading-order heavy-to-light Higgs decay width Γ LO H→hh in the s α − t β plane for a heavy Higgs boson with mass m H = 300 GeV. We can identify three different configurations in which the H → hh mode exactly vanishes [21]: i) the light Higgs decoupling limit, s α = 0; ii) the heavy Higgs decoupling limit, |s α | = 1; and iii) the line t β = −t α . In cases i (resp. ii), all couplings of the heavy (resp. light) Higgs boson eigenstate are identically zero.

Electroweak one-loop corrections
Since all external particles involved in this process are colorless and electrically neutral, the next-to-leading order (NLO) corrections are given by purely weak one-loop effects. These O(α ew ) corrections stem from the interference of the LO amplitude and different subsets of one-loop graphs. On the one hand we have the genuine one-particle irreducible (1PI) vertex corrections. These include triangle and bubble-like three-point topologies which involve the exchange of virtual heavy fermions, weak gauge bosons and Higgs bosons, as generically illustrated in Figure 3. The neutral Goldstone bosons and the SU (2) L Faddeev-Popov ghost contributions appear explicitly in the 't Hooft-Feynman gauge. In addition to the genuine 1PI topologies, the one-loop corrections involve as well the Hhh vertex counterterm, which relies on a combination of Higgs and gauge boson two-point functions, as discussed beforehand in Section 3. This contribution cancels the UV-divergent poles of the 1PI amplitude and allows us to write the complete one-loop amplitude in terms of physical (renormalized) parameters. Lastly, we must include the finite wave-function corrections to the external Higgs boson legs (58) in the minimal field scheme -while for the on-shell schemes these are identically zero. Combining all these pieces we may express the NLO heavy-to-light Higgs decay width as Hhh + 2 Re λ Hhh δΓ Hhh + δΓ WF Hhh + δλ Hhh .
By δΓ Hhh we denote the one-loop contribution from the 1PI three-point vertex graphs. The wave-function corrections yield where we have introduced the finite field renormalization constants Eq. (58) and expanded them to first order in α ew . Finally, δλ Hhh stands for the counterterm of the trilinear scalar coupling. The latter is constructed from the tree-level expression (71), expanding all the bare quantities as customary as X 0 → X + δX. Doing so we find δT H + c Hhh where the coefficients c i are quoted separately in the Appendix.
The relative size of the quantum effects is quantified through the ratio where all quantities are given in the α em -parametrization. The pure one-loop corrections ∆Γ 1-loop include all terms stemming from the LO-NLO interference.

Phenomenology
Hereafter we describe the phenomenology of heavy-to-light Higgs decays at NLO EW accuracy. We begin in Section 5.2 by completing the discussion on the gauge dependence issues that were pointed out qualitatively in Section 3. We here revisit them on quantitative grounds and justify the choice of the improved on-shell scheme as our default setup for the remainder of the analysis. Furthermore, we perform a dedicated numerical comparison of different schemes and show that these theoretical shortcomings have arguably a negligible impact in practice.
We continue in Sections 5.3 and 5.4 with a detailed presentation of our phenomenological analysis. In line with Ref. [21], we separately consider two regions of interest, where heavy-to-light Higgs decays are kinematically accessible.
• High-mass region: in which the lighter eigenstate is identified with the discovered SM-like Higgs of (fixed) mass m h , while the heavier mass-eigenstate corresponds to an additional heavy Higgs companion with a variable mass m H , such that m H > 2 m h .
• Low-mass region: where one instead identifies the heavier mass eigenstate with the SM-like Higgs of (fixed) mass m H , while h represents now a light Higgs companion and m H > 2 m h .
Specific scenarios with maximal H → hh branching fractions in agreement with all of the model constraints are analysed separately in Section 5.5.

Computational Setup
In the remainder of our numerical analysis, we fix the SM Higgs boson mass to the best-fit value based on the combined data samples of the ATLAS and CMS experiments m h = 125.09 GeV [65]. Whenever GeV. This is in fact equivalent to defining tan β in the G F -parametrization. To perform our calculation in the α em -parametrization, we must translate it accordingly through Eq. (43) tan β where Plugging the above relation along with Eq. (43) into the expression for the decay width (73), which relates the α em and G F parametrizations up to NLO EW accuracy through Feynman rules for the singlet model rely on two independent implementations. For one of them we use LanHEP [74,75] and Sloops [15][16][17][18] and include a non-linear gauge fixing Lagrangian (20). For the second one we generate UFO [128] and FeynArts [127] files using FeynRules [129], while the counterterms are derived analytically and implemented by hand. Both implementations are in perfect agreement.
The one-loop decay amplitude is generated with FeynArts and analytically processed via Form-Calc [127]. The loop form factors are handled with dimensional regularization in the 't Hooft-Veltman scheme, and written in terms of standard loop integrals. These are further reduced via Passarino-Veltman decomposition and evaluated with the help of LoopTools [130].    Table 5: Dependence on the gauge-fixing parameters of the mixed mass counterterm δm 2 hH (in GeV 2 ) within the different renormalization schemes introduced in Section 3.3.6. The model parameters are fixed as in Eq. (81). For the (scale dependent) minimal field scheme we set the renormalization scale at µ 2 R = (m 2 h + m 2 H )/2.

Scheme choice and gauge invariance
Gauge-fixing parameters may appear explicitly at intermediate stages in the calculation of S-matrix elements in gauge theories. Taken separately, counterterms and unrenormalized loop amplitudes may in general depend on the gauge-fixing parameters and are eventually also UV-divergent. We only expect these UV divergent contributions to cancel once all the different building blocks are combined together into predictions for physical observables. Nonetheless, depending on which renormalization conditions are chosen for a certain input parameter X, one may obtain loop amplitudes which, albeit finite, still depend on the gauge-fixing. These situations reflect that, for such a renormalization scheme, the definition for X is gauge-dependent. In the following we check the different renormalization schemes introduced in Section 3.3.6 in the light of gauge independence. We compute the one-loop correction to the heavy-to-light Higgs decay width which gives a leading-order width Γ LO H→hh = 0.137 GeV. In Table 4 Table 4. This breakdown can be ultimately traced back to the renormalization condition that determines the mixed mass counterterm δm 2 hH . Its definitions in the minimal field scheme (57) and the OS scheme (62) are not gauge-independent, and lead to a {nlgs}-dependent decay width. Instead, we find no residual {nlgs}-dependencies in the mixed MS/OS and the improved OS schemes, in which δm 2 hH is fixed via the gauge-independent definitions of Eq. (68) and (69) respectively. We make these observations patent in Table 5, where we display the numerical value of the mixing counterterm δm 2 hH corresponding to the four renormalization schemes under analysis. Since the counterterm is not UV finite, we split it into a finite and singular part as (with ∆ as defined in (65)) Neither the coefficient of the UV pole δm 2 hH ∞ nor the finite remainder δm 2 hH fin depend on the gaugefixing parameters when we fix δm 2 hH either in the mixed MS/OS or the improved OS conditions. Instead, both terms are shifted when we switch from the linear {nlgs} = 0 to the non-linear gauge-fixing choice {nlgs} = 10, when the calculation is performed using the minimal field or the OS schemes. In view of the fact that δm 2 hH (along with the mixed field renormalization δZ hH , cf. Table 3) are the only different ingredients between these four schemes, they are ultimately responsible for the finite {nlgs}-dependent remainders in δΓ 1-loop in the latter two schemes -and linked to them, of the uncancelled UV poles. We emphasize as well that these concomitant UV divergences vanish for {nlgs} = 0 and hence do not feature in the customary 't Hooft-Feynman gauge, where the results in all schemes are UV finite. Finally, let us also notice that, given the relation between the mixed mass counterterm and the mixing angle via Eq. (51), a gauge-independent δm 2 hH supports a more physical interpretation of the mixing angle, viz. as value that could be extracted from e.g. a deviation in the LHC Higgs signal strengths or, alternatively, an excess which points to the direct production of the heavy scalar. j For practical purposes, therefore, the proven robustness of the improved OS scheme (giving in all cases UV-finite, {nlgs}-independent, and numerically stable renormalized one-loop amplitudes) justifies its use as default scheme choice in our numerical analysis hereafter. Moreover, the excellent agreement between the δΓ 1-loop results for the different schemes in the linear 't Hooft-Feynman gauge -as explicitly shown further down -give convincing arguments that also the schemes where the mixed mass counterterm is gauge dependent render reliable results -at least as long as the linear 't Hooft-Feynman gauge is used and, in the case of the minimal field scheme the renormalization scale is chosen in the ballpark of the i Using double precision we expect an agreement of 14 to 15 digits. Given the variation ranges ∆ = 0 . . . 10 7 and {nlgs} = 0 . . . 10, we deem the test as satisfactory if 6 to 8 common digits are achieved.
j Similar lines of argument are used in the context of the renormalization of the tan β parameter in the MSSM, cf. e.g. Table 2 in Ref. [17].
relevant physical scales. This is again in line with analogue situations such as e.g. the squark sector of the MSSM [18,121,122].  Table 6: Heavy-to-light Higgs decay width Γ H→hh at LO and NLO EW accuracy, for representative parameter choices and different renormalization schemes, in the high-mass region. The total decay widths are obtained in the αem-parametrization, as defined in Eqs. (76), while the relative one-loop EW effects are quantified in both the αem-parametrization and the G F -parametrization, cf. Eq. (80). For the (scale-dependent) minimal field scheme, the renormalization scale is fixed to the geometrical average mass µ 2 R = p 2 * = (m 2 h + m 2 H )/2. The input value for tan β is linked to the singlet vev through vs = (

High-mass region
In Table 6 we evaluate Γ LO H→hh and Γ NLO H→hh for representative parameter choices and different renormalization schemes. The relative one-loop EW corrections are given in both the α em -parametrization and the G F -parametrization introduced in Section 4. Our results show decay rates that strongly vary with the relevant parameters of the model. The heavy-to-light Higgs decay width significantly depends on the decaying Higgs mass m H , changing by two orders of magnitude when sweeping the range m H = 300 . . . 700 GeV. For heavy Higgs masses close to the di-Higgs threshold, the partial Higgs widths lie in the ballpark of O(0.01 − 0.1) GeV. These results depend as well on the mixing angle, and change by roughly one order of magnitude from small (viz. sin α 0.1) to moderate mixing angles (viz. sin α 0.3). For larger heavy Higgs masses, the Γ NLO H→hh values may rise up to the few GeV level. The mild numerical discrepancies between the different schemes are indicative of small theoretical uncertainties in the 't Hooft-Feynman gauge. For a more general gauge-fixing choice, though, the minimal field and on-shell schemes are no longer reliable, in view of their proven gauge-dependent nature. It is also worth noticing that the radiative corrections in the G F -parametrization (δ G F ) are generically smaller than in the α em -parametrization. The reason is twofold: i) part of the NLO EW corrections in the latter case (δ α ) are contained in the ∆r parameter, and hence already embedded into the G F -scheme LO calculation (cf. Eq. (79)). Consequently, the quantum effects encoded by ∆r do not belong to δ G F ; ii) for phenomenologically relevant scenarios, ∆r is dominated by purely SM effects, for which ∆ SM > 0 [98,[103][104][105][106][107], and thereby δ α > δ G F , given the relation between both (80). The analysis is complemented in Figure 4 with a thorough survey of the parameter space dependencies. The NLO results are calculated in the improved OS scheme. The shaded regions are ruled out by different theoretical and experimental constraints on the model: i) the ranges m H > 840 GeV (left panel) and tan β < 1.27 (right panel) are excluded by perturbativity ii) | sin α| > 0.31 (central panel, green shading) is incompatible with electroweak constraints from the m W measurement; finally, the central range | sin α| < 0.06 (central panel, orange shading) is incompatible with vacuum stability. Fig. 4 can be readily traced back to the LO dynamics which governs the decay process. The two key players, as alluded to above, are the trilinear Higgs self-coupling λ Hhh and the characteristic 1 → 2 kinematics. The former is responsible for the quadratic growth Γ H→hh ∼ O(m 2 H ) (cf. Eq. (84)), which overcomes the phase space suppression at m 2 H m 2 h , and explains the power-like increase as a function of m H (left panel in Figure 4). The NLO-corrected result with respect to the mixing angle mimics the LO result, with the expected nodes in the decoupling limits | sin α| = 0 or 1 as well as for tan β = − tan α (cf. the central panel of Fig. 4).

Most features observed in
Unlike the stark changes observed for the decay width, the relative one-loop EW corrections are much more stable, positive, and of the order of few percent. Differences between the α em -parametrization and the G F -parametrization, as well as between the different renormalization schemes, are mild and remain typically below the percent level. The slight kink in δ α for m H 350 GeV (left panel, Fig. 4) reflects the top-quark threshold. The finite correction δ α ∼ 3% for sin α → 0 (cf. the lower subpanel of Figure 4, center) follows from the fact that both Γ LO H→hh and Γ NLO H→hh tend to zero in this limit, while its ratio remains roughly constant. The unphysical large effects at sin α −0.9 are due to the LO node in the limit tan α → − tan β, for which λ Hhh = 0. The pronounced NLO slope at low tan β is ultimately due to the exchange of virtual Higgs bosons, and constitutes a telltale imprint of the singlet model dynamics at the one-loop level. While the fermion and the gauge boson-mediated contributions are all controlled by (globally rescaled) gauge couplings, the size of the Higgs-mediated loops is governed by the Higgs self-couplings. These are strongly enhanced for tan β 1, specially the Higgs boson two-point graphs, which depend on them quadratically. For low enough tan β values, e.g. typically tan β 0.3 and for m H 300 GeV, the relative yield δ α exceeds ∼ 50%, indicating that the process becomes effectively loop-induced. Such sizable loop effects are nonetheless hampered in practice, owing to the unitarity and perturbativity bounds which severely constrain the phenomenologically viable low-tan β range. The limit tan β 1 corresponds in fact to the onset of a strongly-coupled regime, in which at least one of the scalar self-couplings becomes non-perturbative, cf. also the discussion in Sec. 2.6.
Complementary vistas to the H → hh landscape are displayed in Figure 5. Here we show the relative one-loop effects δ α (76) as density contours in the sin α − tan β plane. The yellow contour separates the allowed and excluded regions in the parameter space. Only the horizontal fringes enclosed by the contour are compatible with all constraints on the model. The white voids stand for values of δ α 100% and correspond to regions where δ α is no longer a meaningful measure of the relative quantum effects, while it instead indicates that the decay process becomes loop-induced. We find this situation: i) along the strip tan α − tan β, due to the suppressed tree-level couplings; and ii) for tan β < 1, due to the cot β-enhanced Higgs-mediated loops.
The impact of heavy-to-light Higgs decays on the decay pattern of the heavy Higgs state is portrayed in Figure 6. The branching ratios for the leading decay channels are represented as a function of the heavy Higgs mass. The mixing angle and tan β values are fixed in each panel such that they maximize the H → hh branching ratio for a given heavy Higgs mass [90], as explicitly indicated in the figure. In this plot, we show the partial decay widths to SM fields by rescaling the SM predictions [44], while for  the H → hh we use the LO result k . As well known, bosonic modes dominate the Higgs boson decays at high masses [7][8][9]. We find a rather featureless profile, with WW being the leading mode and with roughly no changes over the whole mass range. Only the decays into top-quark pairs are also competitive, and attain up to BR ∼ O(10)%. The remaining fermionic channels, as well as the loop-induced γγ, γZ and gg modes, stagnate at the O(0.1)% level or below and are not shown. In the lower subpanels we show the relative one-loop EW correction to the heavy-to-light Higgs decay width for the same parameter variation.
Finally, in Figure 7 we analyse how Γ NLO H→hh varies with the renormalization scale in the (scaledependent) minimal field scheme, as introduced in Eq. (57). We compare the minimal field to the (scale-independent) improved OS scheme, which we show as reference value. For tan β = 5, Γ NLO H→hh flattens not far from the geometrical average mass scale µ 2 R p 2 * = (m 2 h + m 2 H )/2. Precisely around this value, both the minimal field and the improved OS predictions tend to converge, suggesting that µ 2 R = p 2 * in Eq. (57) is indeed a convenient scale choice for the former. Moreover, the very stable NLO predictions around this scale, added up to the mild changes with the different renormalization schemes shown in Table 6, indicate a small theoretical uncertainty.
Much steeper scale variations arise instead in the tan β < 1 region. Here, the Γ NLO H→hh predictions become unstable, especially for heavy Higgs masses. Such instability may be once more traced back to the Higgs-mediated scalar two-point graphs: these become overly large owing to the enhanced Higgs selfcouplings, and artificially dominate the scale dependence in these regions. Such a stark scale dependence is simply the reflect of the poor perturbative behavior of the model in the vicinity of a strongly-coupled regime λ ∼ O(4π), which obviously translates into a huge theoretical uncertainty.
k A global study including all Higgs decays in the singlet model to state-of-the-art accuracy lies beyond the scope of the present study and will be discussed in a forthcoming publication [131].

Low-mass region
Assuming now m H = 125.09 GeV and a free light Higgs mass m h , direct LEP and LHC mass bounds, and most remarkably the measured LHC Higgs signal strengths, narrow the viable sin α region down to a slim fringe | sin α| 1. Constraints become particularly tight in the region of interest m H > 2m h , given the limited tolerable room for deviations in the total SM-like Higgs width when additional decay modes feature. State-of-the-art LHC constraints on the total Higgs width place an upper limit of Γ h ≤ 22 MeV [132,133]. For definiteness, we hereafter adopt the fiducial choice sin α = 0.998 [21].
In Figure 8 we examine the parameter space dependence of Γ H→hh in this scenario. Complementarily, in Table 7 we provide precise predictions for specific parameter space points, while comparing again the different renormalization schemes. In Figure 9 we analyse the m h − tan β interplay by showing the total NLO amplitude Γ NLO H→hh (73) and the relative NLO correction δ α (76) in the form of two-dimensional density maps. The obtained Γ NLO H→hh values span two orders of magnitude, ranging from O(10 −3 ) down to O(10 −5 ) GeV as we navigate throughout the different parameter space regions. This sharp variation is again connected to the behavior of the leading-order coupling λ Hhh : whilst it tends to zero in the limit sin α → 1, it can yet contribute if the cot β-enhanced terms are large enough. Either way, let us once more recall that a significant patch of the low-tan β range is in practice precluded by the different constraints on the model (see e.g. the top panels of Fig. 8 and Fig. 9). In particular, the shaded regions at small tan β and low Higgs masses are incompatible with the LHC Higgs signal strength measurements. Another salient feature is the steep rise of the quantum corrections at low tan β (see the top panels of Fig. 8): these are positive, tend to increase with the light Higgs mass, and may surmount the O(50%) level. Aside from the discussed Higgs-mediated loop enhancements, additional mechanisms reinforce this behavior in this case: i) the suppressed tree-level decay amplitude, due to the lesser phase space available, the closer we move to the kinematical threshold; ii) the vicinity of the di-Higgs loop threshold, which invigorates the light Higgs-mediated loops even further. For tan β > 1 the corrections are instead moderate and negative, becoming even more so for very light m h values. The latter effect may be traced back to the fermionic (viz. the top-mediated ) three-point loops, which are in this case the dominant source of quantum corrections and present a trademark logarithmic dependence ∼ log(m 2 t /m 2 h ).  The scale dependence of the NLO results in this low-mass region is analysed in Figure 10. The Γ NLO H→hh predictions obtained in the improved OS and the minimal field schemes schemes agree very well in the ballpark of the geometric average mass µ 2 R = p 2 * = (m 2 h + m 2 H )/2, and the latter barely varies with the scale. The very stable slope even in the tan β < 1 region, which is in contrast to the strong scale dependence in the high-mass region, is explained by the much lower scales µ 2 R p 2 * involved in this case, for which the finite Higgs-mediated contributions to the Higgs field two-point functions are much smaller.
In Figure 11 we recast the above analysis in terms of the heavy Higgs branching ratios. We track down their behavior as a function tan β for exemplary light Higgs masses and fiducial mixing sin α = 0.998. From values of tan β 1 onwards, the obtained decay pattern approaches that of a purely SM-like Higgs boson. The dominant mode is bb, while the di-Higgs final state is hampered due to the tiny tree-level  Table 7: Heavy-to-light Higgs decay width Γ H→hh at LO and NLO EW accuracy for representative parameter choices and renormalization schemes in the low-mass region. The total decay widths are obtained in the αem-parametrization, as defined in Eqs. (76), while the relative one-loop EW effects are quantified in both the α-parametrization and the G F -parametrization, cf. Eq. (80). For the (scale-dependent) minimal field scheme, the renormalization scale is fixed to the geometrical average mass µ 2 R = p 2 * = (m 2 h + m 2 H )/2. The input value for tan β is linked to the singlet vev through vs = ( √ coupling λ Hhh ∼ c α . In this case, the H → hh mode carries not more than a few percent of the total budget -on equal footing with the loop-induced decay H → gg. If we instead move towards lower tan β values, the cot β-enhanced terms overcome in part the sin α → 1 suppression and promote H → hh again to a chief role.

Maximal branching ratios
So far we have discussed the general behavior of the NLO EW corrections to the heavy-to-light Higgs decay width along the relevant parameter space directions sin α, tan β, m h/H . Before closing, we focus on the series of benchmarks with maximal tree-level heavy-to-light Higgs branching ratio proposed in [90]. These are defined as a function of the heavy Higgs through the parameter choices quoted in Tables 8, for the high and the low mass regions respectively. In these regimes, the decays of the heavy Higgs state provide a particularly interesting phenomenological ground for studying finite width effects and lineshape modifications in the production of a heavy scalar resonance, cf. e.g. [135][136][137][138][139][140] l . Numerical predictions for Γ LO H→hh and Γ NLO H→hh , together with the relative one-loop correction in the two parametrizations δ αem and δ GF (80) are provided in Tables 9 and 10, for the high and low mass regions respectively. Complementarily, we list down the corresponding branching fractions for the additional decay modes (barring those channels below 0.01%). Let us recall that for the latter we use the rescaled partial widths from [44], while for H → hh we quote the LO result in the α em -parametrization, in line with Figures 6 and 11 Table 8: Maximal branching ratios for the heavy-to-light Higgs decay mode H → hh in the high-mass (left) and low-mass regions (right) as proposed in Refs. [90,91]; the results quoted here are obtained in the setup of the mentioned references, evaluating Γ(H → hh) at LO in the G F -parametrization. Note that the maximal branching ratios are determined for a maximal mixing, to ensure a large production rate. In this case, the lower limit of tan β is mainly determined by the requirement of perturbativity for λ 2 , cf. the extensive discussion in [20]. The same strategy was followed for the low-mass region, where again for fixed sin α values the minimal value of tan β is determined. Here, the lower limit on tan β stems from the signal strength fit.  Table 9: Heavy-to-light Higgs decay width Γ H→hh at LO and NLO EW accuracy for the maximal branching fraction scenarios in the high-mass region given in Table 8. The relative one-loop EW effects are quantified in both the αem-parametrization and the G Fparametrization, as defined in Eqs. (76)- (80). Renormalization is performed in the improved on-shell scheme. In the right columns we document the branching ratios (in %) for the leading Higgs decay channels and the total Higgs width. Like in Figures 6 and 11, all partial decay widths to SM fields are evaluated by rescaling the SM predictions [44], while for H → hh we use the LO result evaluated in the αem-parametrization. All partial widths are given in GeV.
l See also e.g. Ref. [53] in the context of Higgs pair production.  Table 10: As in Table 9 for the low-mass region. Notice that in this case all partial widths are given in MeV.

Summary
Heavy-to-light Higgs decays H → hh are of undisputed relevance in the phenomenological characterization of extended Higgs sectors. When kinematically accessible, these may contribute to, and in some scenarios even dominate, the heavy Higgs lineshape, while at the same time they significantly modify its decay pattern with respect to the SM picture. On the other hand, both the tree-level and the leading one-loop contributions to this process are governed by the scalar self-interactions, which makes this decay a unique handle on the architecture of the scalar potential.
While electroweak corrections to the Higgs self-couplings have been the subject of dedicated analyses in the 2HDM [141][142][143], the NMSSM [144] or the Inert Doublet model [145], a corresponding study for the singlet extension was lacking. Extending upon previous work [19], we have presented herewith a detailed analysis of the heavy-to-light Higgs decays at NLO electroweak accuracy. To renormalize the singlet-extended Higgs sector we have proposed four renormalization schemes: i) a minimal field setup; ii) a traditional on-shell prescription; iii) a mixed MS/on-shell scheme; and iv) an improved on-shell scheme. Using the general non-linear gauge-fixing of Sloops, we have discussed in detail the gauge independence of the different renormalization setups. We have found that, while the minimal field and on-shell approaches still lead to a residual dependence on the non-linear gauge fixing parameters, the mixed MS/OS and improved OS schemes render gauge independent one-loop predictions for physical observables. Furthermore, the improved OS scheme is numerically stable in all regions of the parameter space. We therefore advocate for the use of this scheme to investigate the phenomenology of singlet extensions of the SM at higher orders.
We have applied the above schemes to compute the corresponding heavy-to-light Higgs decay widths Γ H→hh including one-loop electroweak corrections. We have performed a comprehensive phenomenological analysis, with a separate study of two possible realizations of the model: a high-mass and a low-mass region, in which the additional scalar field corresponds to a heavy (resp. a light) companion of the SM-like Higgs boson.
The phenomenological implications of our study can be summarized as follows: • The heavy-to-light Higgs decay width at LO is governed by two competing mechanisms: i) the Higgs self-coupling strength λ Hhh ; ii) the one-to-two body decay kinematics. We pinpoint a strongly varying width with the relevant model parameters. Overall, Γ H→hh may attain up to O(1) GeV for heavy Higgs masses above m H ∼ 500 GeV.
• Aside from the tan β < 1 region, the relative one-loop effects are mild and show tempered variations over the parameter space. In the high-mass region, electroweak corrections are positive, loosely variable, and stagnate in the ballpark of few percent. In the low-mass region, mainly for tan β > 1 and small m h values, these may be pulled down to δ α ∼ −10%.
• For certain parameter choices, the H → hh decay becomes effectively loop-induced: i) along the tree-level nodes where the LO contribution vanishes; ii) at low tan β, where the cot β-enhancements lead to increased scalar self-couplings, and thereby to large Higgs-mediated loop graphs. In practice, though, these sizable quantum effects are precluded once the constraints on the model are taken into account.
• Let alone extreme parameter space corners, the NLO predictions are robust under changes of renormalization schemes and renormalization scale choices, as indicative of a small theoretical uncertainty.
Having constructed a complete renormalization scheme for the Higgs sector, the path ahead is clear for further analyses on the topic. On the one hand, it will be interesting to further explore the role of the quantum effects on the Higgs self-couplings themselves, and whether these may have relevant implications e.g. for collider searches or in electroweak baryogenesis. On the other hand, the complete renormalization of the Higgs sector paves the way towards characterizing the singlet model phenomenology at one-loop electroweak accuracy, including all Higgs production modes and decay channels, and exploiting the rich possibilities of off-shell effects. Work in this direction is underway [131].