Precision calculations for $h \to WW/ZZ \to 4$ fermions in a Singlet Extension of the Standard Model with Prophecy4f

We consider an extension of the Standard Model by a real singlet scalar field with a $Z_2$-symmetric Lagrangian and spontaneous symmetry breaking with vacuum expectation value for the singlet. Considering the lighter of the two scalars of the theory to be the 125GeV Higgs particle, we parametrize the scalar sector by the mass of the heavy Higgs boson, a mixing angle $\alpha$, and a scalar Higgs self-coupling $\lambda_{12}$. Taking into account theoretical constraints from perturbativity and vacuum stability, we compute next-to-leading-order electroweak and QCD corrections to the decays $h\to WW/ZZ \to 4$ fermions of the light Higgs boson for some scenarios proposed in the literature. We formulate two renormalization schemes and investigate the conversion of the input parameters between the schemes, finding sizeable effects. Solving the renormalization-group equations for the MSbar parameters $\alpha$ and $\lambda_{12}$, we observe a significantly reduced scale and scheme dependence in the next-to-leading-order results. For some scenarios suggested in the literature, the total decay width for the process $h \to 4f$ is computed as a function of the mixing angle and compared to the width of a corresponding Standard Model Higgs boson, revealing deviations below 10%. Differential distributions do not show significant distortions by effects beyond the Standard Model. The calculations are implemented in the Monte Carlo generator Prophecy4f, which is ready for applications in data analyses in the framework of the singlet extension.


Introduction
or cosine of a mixing angle. The mass of the additional Higgs boson, the Higgs mixing angle, and one coupling factor of the scalar self-interactions parametrize the extended sector. In our phenomenological study, the lightest of the two Higgs bosons of the theory is considered to be the 125 GeV resonance observed at the LHC, but our theoretical approach is not restricted to this case.
Our goal is to perform next-to-leading order (NLO) computations within the SESM, including both EW and quantum chromodynamics (QCD) corrections. To this end, it is necessary to renormalize the theory. Recently, the renormalization of SM extensions has been subject of discussion, since very often there is no obvious formulation of on-shell (OS) renormalization conditions to define BSM parameters that are fully based on physical S-matrix elements. To define these parameters, it is customary to use renormalization conditions in the modified Minimal Subtraction (MS) scheme. While a consistent use of OS conditions based on physical S-matrix elements guarantees a gauge-independent parametrization of physical observables by renormalized input parameters, making use of both OS and MS conditions in the definition of a renormalization scheme can lead to gauge-dependent renormalization constants if the (gaugedependent) tadpoles are not treated properly. In Ref. [39], Fleischer and Jegerlehner proposed a renormalization scheme to avoid this problem in the SM, followed by other approaches such as, e.g., described in Ref. [40]. Recently, these strategies have been applied to Two-Higgs-Doublet Models (THDMs) in Refs. [41][42][43][44].
Different renormalization schemes for the SESM with a real singlet scalar were already considered in Refs. [44][45][46][47][48]. Among these proposals there is no convincing scheme that is fully based on OS conditions. Most schemes are based on ad hoc or on MS conditions, and many variants still suffer from gauge dependence issues. In this paper, we build on OS conditions as far as possible and take MS conditions for those parameters for which no distinguished OS conditions are available. 1 We formulate two renormalization schemes for the SESM, using two different ways to treat tadpole contributions, one of them based on the FJ variant [39], similar to a proposal made in Ref. [44]. We analyze, in both cases, the dependence of our NLO results on the renormalization scale µ r , which is due to MS definitions of the Higgs mixing angle and the Higgs self-coupling λ 12 . We study the conversion of input parameters between the two schemes and compare the results obtained in the two schemes to inspect the perturbative consistency of the chosen region for µ r . Parameter conversions between different renormalization schemes, and the corresponding scheme dependence of NLO results, were not yet discussed for the previously proposed schemes and their applications. In a situation where no distinguished renormalization scheme has yet emerged, discussions of renormalization scale and renormalization scheme dependences are very important in applications.
In this work, we compute decay observables for the decays of the 125 GeV Higgs boson of the SESM into four fermions via intermediate (off-shell) WW or ZZ states. These processes played a central role in the discovery of the Higgs boson, and are very important channels in Higgs couplings analyses. NLO computations for these processes were performed, including both EW and QCD corrections, in the SM with the Monte Carlo generator Prophecy4f [51][52][53] (and matched to a QED parton shower in Ref. [54]), as well as in presence of a fourth fermion generation [55] and in the THDM [44,56], but there are no corresponding results available yet in the SESM.
Using the Mathematica package FeynRules [57,58], we have implemented the Feynman rules for the SESM into a FeynArts [59] model file, which includes one-loop counterterms (CTs) for both renormalization schemes. The model file has been used to compute, with the Mathematica package FormCalc [60], the matrix elements for the decays of the light Higgs into four fermions (including EW and QCD NLO corrections) and to produce Fortran routines to extend the Monte Carlo program Prophecy4f. Finally, we have used the extended Prophecy4f version to compute partial widths and to generate differential distributions for a selection of benchmark scenarios, proposed in Refs. [3,37]. Other computations of NLO EW corrections relevant for phenomenology in SESMs can be found in Refs. [44][45][46][47][48].
The paper is organized as follows. In Sect. 2, we describe the Lagrangian of the SESM considered in this work and discuss the basic features of the model. We set up the renormalization procedure introducing renormalization transformations for fields and parameters in Sect. 3 and fixing the renormalization constants in Sect. 4. In Sect. 5, we briefly describe the capabilities of the Monte Carlo program Prophecy4f and discuss the computation of decay observables for the processes h → WW/ZZ → 4 fermions. In Sect. 6, we provide the input parameters and discuss the benchmark scenarios used to produce the numerical results shown in Sect. 7. In Sect. 8 we draw our conclusions. Additional numerical results are reported in the appendices.

Singlet Extension of the Standard Model
Singlet extensions of the SM add one or more scalar singlet fields to the SM Higgs sector. In general, the scalar fields can be complex, but we here consider the simplest variant of a SESM with one additional real field. The Lagrangian of the model can be easily obtained by modifying the scalar sector of the SM. As the extension does not modify strong interactions, we present only the contributions to the total SESM Lagrangian that are relevant for the EW sector, while the QCD part can be taken over from any standard reference, such as Ref. [61]. For convenience, we split the EW Lagrangian of the model as follows: (2.1) By gauge symmetry and renormalizability constraints, only small deviations from the SM are allowed: The additional scalar field, which has mass dimension 1 and is neutral under SU(2) W × U(1) Y gauge transformations, can be coupled to the SM fields only through gauge-invariant terms with mass dimension of at most 3, i.e. the singlet scalar can only couple to Φ † Φ, with Φ denoting the SM-like Higgs doublet. Moreover, the Higgs Lagrangian contains a kinetic term and self-interaction terms of the singlet field in addition to the terms that are already present in the SM. Since the singlet enters only the Higgs Lagrangian L Higgs , the gauge, fermionic, Yukawa, gauge-fixing, and ghost terms of Eq. (2.1) only need little adaption from the SM. For these contributions we make use of the formulation of Ref. [62]. In Sect. 2.1 we introduce the conventions used for the Higgs Lagrangian, then we briefly discuss the gauge (Sect. 2.2) and the fermion sectors (Sect. 2.3). Finally, we define our input parameter set in Sect. 2.4.

Mass spectrum
In complex SESMs, the singlet field is supposed to be responsible for the symmetry breaking in a hypothetical hidden sector, where interactions of the hidden sector are governed by an exotic U(1) or even higher-rank gauge symmetry [12,14,18,32,34]. Considering our SESM as a downgrade of such a more comprehensive theory which still carries salient features of the more complete theory, we enforce a Z 2 symmetry on the Higgs Lagrangian (under sign change of the singlet field). Moreover, we require a non-vanishing vev for the singlet. The model obtained in this way is very simple, but still phenomenologically interesting, as it involves mixing of the new singlet scalar with the Higgs field of the SM scalar sector, which is a generic feature of many SM extensions. The most general renormalizable, gauge-and Z 2 -invariant Higgs Lagrangian in presence of one real singlet and one doublet is given by [17,28] L Higgs = L Higgs,kin − V (Φ, σ), where Φ is the complex SM scalar doublet with hypercharge Y W,Φ = 1, and σ is the real singlet field. Splitting off the vevs v 2 and v 1 , we parametrize Φ and σ by where φ + , φ − = (φ + ) † , and χ denote the would-be Goldstone-boson fields for the W ± and Z bosons. The covariant derivative, which is relevant for the couplings of Φ to the EW gauge bosons, is given by In Eq. (2.4), g 2 , I a W (a = 1, 2, 3), and W a µ are, respectively, the gauge coupling, the generators, and the gauge fields of the weak isospin SU(2) W group; g 1 , Y W , and B µ are the gauge coupling, the generator, and gauge field of the weak hypercharge U(1) Y group. 2 For convenience, in Eq. (2.3), we introduce the subscripts 1 and 2 to label vevs and fields for the singlet and the doublet sectors, respectively. The role of v 2 is the same as in the SM, since only Φ couples to the EW gauge bosons. The vev of the singlet, v 1 , can be vanishing or non-vanishing, providing different phenomenology. As stated before, we consider the case v 1 = 0, which is phenomenologically more interesting. Two non-vanishing vevs can only arise for µ 2 2 , µ 2 1 > 0 and if the following vacuum stability conditions are fulfilled, λ 2 > 0, λ 1 > 0, λ 2 λ 1 − λ 2 12 > 0. (2.5) Expanding the potential V with the decomposition (2.3) and keeping terms containing the fields h 2 , h 1 up to second order leads to with the tadpole parameters and the non-diagonal mass matrix (2.8) In order to work with fields related to mass eigenstates, we diagonalize the mass matrix by a rotation about an angle α with −π/2 ≤ α ≤ π/2, where the shorthand notations s nα ≡ sin (nα) , c nα ≡ cos (nα) are used. Using the parameterization (2.15), the vacuum stability conditions (2.5) are automatically fulfilled, while µ 2 2 > 0 and µ 2 1 > 0 lead to a restriction on the parameter space of the theory by . (2.16)

Couplings
We insert the expressions (2.3) and (2.9) as well as Eq. (2.7) with t 2 = t 1 = 0 into the potential of Eq. ( with the coupling constants given by The couplings for the full Higgs Lagrangian (2.2) can be obtained from the Higgs Lagrangian of the SM following a simple procedure: Firstly, all the couplings coming from the SM Higgs potential are removed and replaced with the couplings reported in Eqs. (2.17) and (2.18); then the couplings of the light Higgs h (heavy Higgs H) to the other fields are obtained by rescaling the SM Higgs couplings by a factor of c α (s α ). Thanks to its simplicity, many tree-level computations in the SESM can be easily obtained from the SM results via rescaling coupling factors appropriately. However, when considering multi-Higgs processes or loop contributions, some care has to be taken.

Gauge sector
In the total Lagrangian (2.1) we use the gauge, gauge-fixing, and ghost Lagrangians of Ref. [62] (with corresponding generalized Higgs couplings in the latter). Since the singlet field σ does not couple to the gauge bosons, the gauge-boson mass terms are generated through the interactions of the EW gauge bosons with the vev v 2 of the Higgs doublet in the same way as in the SM. After a rotation about the EW mixing angle θ W (quantified here by c W ≡ cos θ W and s W ≡ sin θ W ) of the gauge fields into fields related to mass and charge eigenstates, the usual relation among the EW coupling constants, e = g 2 s W = g 1 c W , (2.19) ensures the presence of a massless field A µ , the photon, coupling to fermions as in pure QED via the electric unit charge e. To avoid confusion with the mixing angle α, we denote the electromagnetic coupling constant with α em = e 2 /(4π). Moreover, the vev v 2 is related to the W-and Z-boson masses M W and M Z and to the EW mixing angle as follows, (2.20)

Fermion sector
The form of the fermion and Yukawa Lagrangians, L Fermion and L Yukawa , of Eq. (2.1) is the same as in the SM (see Ref. [62]). In contrast to the SM case, the Higgs doublet contains a mixture of the fields h and H (according to Eqs. (2.3) and (2.9), h 2 = c α h + s α H). Consequently the Yukawa Lagrangian provides two copies of the SM Higgs couplings to fermions, one rescaled by a factor c α for the light field h, and another rescaled by s α for the heavy field H. The free parameters of the fermion sector are the elements of the Yukawa matrices, which are related to the fermion masses and to the elements of the Cabibbo−Kobayashi−Maskawa (CKM) matrix.
In the SESM the CKM matrix can be treated exactly as in the SM case (see e.g. Ref. [62]).

Input parameters
The Lagrangian of the gauge and Higgs sector, L Gauge + L Higgs , contains the seven free input parameters Building on the standard procedure in the SM, we derive the gauge couplings g 1 , g 2 , and one combination of parameters of the Higgs potential (viz. the vev v 2 ) via the relations ( The additional free parameters of the fermionic sector, contained in L Fermion + L Yukawa , all originate from the Yukawa coupling matrices as in the SM and are derived from the fermion 3 Actually, e is derived from the Fermi constant as described below. masses m f and the CKM matrix elements V ij as usual. The input parameter set for the theory used in this work is then given by: (2.22)

The counterterm Lagrangian
Moving beyond leading order (LO) by including NLO effects, the relations between the parameters in the Lagrangian and observables change and, in order to restore the physical meaning of the parameters, it is necessary to renormalize the theory. From now on, we denote the bare quantities introduced in the previous section with a subscript "0", in order to distinguish them from the renormalized quantities defined in this section. We split bare parameters into renormalized parts and corresponding renormalization constants additively and split bare fields multiplicatively into renormalized parts and field renormalization constants. After performing this renormalization transformation in the Lagrangian, all parts containing renormalization constants define the CT Lagrangian δL, which at NLO is linear in all renormalization constants. We use dimensional regularization to treat ultraviolet (UV) divergences. In the SESM, the renormalization of the QCD sector is the same as it is in the SM, so that we will not describe it. For the EW sector, in analogy to Eq. (2.1), the CT Lagrangian can be divided into several contributions, Among the various components of Eq. (3.1), we focus on the contributions deriving from the renormalization of L Higgs , i.e. δL Higgs,kin and δV . For the renormalization of the gauge and fermion sectors we proceed as described in Ref. [62]. Note that in Eq. (3.1) there is no counterpart of L Fix of Eq. (2.1), as the gauge-fixing term is introduced in the Lagrangian after the renormalization. Moreover, for our purpose it is not necessary to compute δL Ghost , since the ghost CTs enter only beyond the one-loop level.

Counterterms from the Higgs potential
In the previous section we have replaced the parameters appearing in the Lagrangian (2.2) by the input parameters M h , M H , λ 12 . Moreover, we introduced a field rotation about the mixing angle α in order to work with the fields h and H, which have diagonal propagators in lowest order. For the proper treatment of tadpoles at NLO, the bare tadpole terms t h,0 and t H,0 should be restored in the tree-level relations of the previous section, most notably in Eqs. (2.15), (2.17), and (2.18). We perform the following renormalization transformations on the free parameters of the scalar potential, Instead of using the renormalization constant δα for the angle α, it is often more handy to use the the following derived renormalization constants for trigonometric functions of α, which are related to δα as follows, Formally all renormalization constants count as O(α em ) corrections, and terms beyond O(α em ) will be dropped. To determine the CTs of all couplings, we have to express the renormalization constants of the original parameters in terms of the renormalization constants of the chosen independent input parameters. Defining renormalization transformations for the original parameters by the renormalization constants, in terms of the renormalization constants defined in Eq. (3.2), are given by where δv 2 is determined in the renormalization of the gauge sector below. In some places we kept the dependent parameters v 2 , v 1 , λ 2 , and λ 1 , in order to keep the result compact. Note that to derive Eq. (3.7) non-vanishing tadpole contributions in Eq. (2.15) had to be restored. However, at NLO the relations given in Eq. (2.15), which are valid for vanishing tadpoles, can be consistently used in all coefficients of renormalization constants in δL, independent of any detail of the renormalization of the tadpoles, since tadpole terms are of O(α em ) or even zero. The tadpole renormalization constants δt 2 and δt 1 are related to δt H and δt h by Applying the renormalization transformations (3.2) and (3.3) to the bare Higgs potential (2.2) and keeping terms linear in the renormalization constants leads to V + δV , where V denotes now the renormalized potential, which has the same analytic form as the bare potential, but contains renormalized quantities (including tadpole parameters) instead of the bare ones. The CT potential δV , which is of O(α em ), reads where the interaction terms can be obtained with the same procedure used to derive the terms linear and bilinear in the fields h, H. Note that the CTs to the scalar self-interaction terms contain explicit tadpole renormalization constants δt h,H after the renormalization transformation.
As, e.g., discussed in Ref. [43], the UV-divergent parts of the field renormalization constants appearing in the kinetic terms of Eq. (3.9) are not all independent. Indeed, a renormalization transformation for the fields Φ and σ would be sufficient to absorb the UV divergences of all Higgs field renormalization constants. In this sense, the field renormalization condition (3.11) We will use some of these expressions to compute UV-divergent contributions to specific renormalization constants. Moreover, these expressions can be used to check internal consistency after the application of the renormalization conditions, which will be discussed in Sect. 4.

Counterterms from the Higgs kinetic term and from the gauge and fermion sectors
Starting from the kinetic term L Higgs,kin of the Higgs Lagrangian (2.2), written in terms of bare parameters and fields, we derive the corresponding CT Lagrangian δL Higgs,kin . For this purpose, we apply the renormalization transformations given in Sect. 3.1, supplemented by the renormalization transformations that are relevant for the gauge sector. The bare W-and Zboson squared masses M 2 W,0 and M 2 Z,0 , and the bare electric charge e 0 are transformed according to Following the "complete on-shell renormalization" of the gauge sector [62], the gauge-boson fields W, Z, and A are renormalized by i.e. we apply a matrix-valued renormalization transformation to the photon−Z system. This has the advantage that no further wave-function or γ-Z mixing corrections need to be applied for external electroweak gauge bosons W, γ, Z if the corresponding field renormalization constants are fixed appropriately, as done in Section 4.
where, for the sake of brevity, we again do not spell out the interactions explicitly. It is important to note that the scalar−vector mixing terms appearing in Eq. (3.14), in contrast to what happens in the bare Lagrangian, are not canceled by the gauge-fixing terms, since the gauge-fixing contribution L Fix is introduced after renormalization, i.e. directly in terms of renormalized quantities. The complete set of renormalization transformations necessary to renormalize the gauge and fermion sectors can be found in Ref. [62]. For a better bookkeeping, renormalization constants for the sine and the cosine of the weak mixing angle are introduced according to Since M W = c W M Z is valid both for bare and renormalized quantities, the renormalization constants δs W and δc W are related to the W-and Z-mass renormalization constants by In the Yukawa Lagrangian, the Higgs field h 2 appearing in the doublet Φ is consistently rotated using Eq. (2.9) and renormalized by the transformations (3.2) and (3.3). In the transition from the SM to the SESM, the renormalization of the CKM matrix does not change; we refer to Refs. [62,65,66] for different formulations.
The fermion fields are renormalized using the simple transformations (in the absence of flavour mixing) where f = ν, l, u, d identify the fermion type, i = 1, 2, 3 is the generation index, and σ = L, R identify left-and right-handed fermion fields.

Renormalization conditions
To fix the renormalization constants introduced in the previous section we adopt, as far as possible, OS renormalization conditions. On-shell renormalization can be performed for all the parameters that are directly accessible by experiments, such as the masses M h and M H in the Higgs sector, but there is no obvious prescription for the mixing angle α and the coupling constant λ 12 . Our choice is to fix these parameters using MS renormalization conditions, so that the corresponding renormalization constants are proportional to the UV divergence where γ E is the Euler−Mascheroni constant and the loop integrals are regularized in D = 4 − 2ǫ dimensions. In this way, we have a mixed OS−MS renormalization scheme, and particular care has to be taken in the renormalization of the tadpole constants.
In the following, we describe two renormalization schemes used in our analysis, which differ in the treatment of the tadpoles. The two schemes have the same OS renormalization conditions in common, which are presented in Sect. 4.1. However, because of the different tadpole handling, there are differences between the two schemes in the renormalization of the parameters α and λ 12 , which are renormalized with MS conditions. In Sect. 4.2.1 we present the MS renormalization conditions for a scheme in which the renormalized tadpole constants t h and t H are set to zero. In this scheme, which is referred to as MS scheme in the following, gauge-dependent tadpole terms enter the relations among bare parameters. While such gauge-dependent terms drop out in the parametrization of observables in terms of input quantities if the latter are defined by OS conditions, gauge dependences can remain in the parametrization of observables if some input parameters are defined by MS conditions. Note that this does not automatically ruin the consistency of this scheme. This shortcoming simply demands that all calculations should be carried out in the same gauge; we always use the 't Hooft−Feynman gauge in this scheme.
MS conditions can be used without introducing gauge dependences in the parametrization of observables if bare (rather than renormalized) tadpoles are set to zero. Following the renormalization of the THDM described in Ref. [43], in Sect. 4.2.2 we describe the changes in the renormalization constants to obtain this gauge-independent scheme and refer to it as FJ scheme, since this variant was proposed by Fleischer and Jegerlehner in Ref. [39]. 4 Variants of this scheme, which are technically different, but phenomenologically equivalent, were also described in Refs. [41,42]. Note that we use the terminology "MS" in two different contexts: as a name for a renormalization scheme and for a type of renormalization condition. Both our "MS scheme" and our "FJ scheme" are based on a mixture of OS and MS renormalization conditions.  are demanded to be zero in our MS scheme. 5 At NLO, these contain contributions from the diagrams illustrated in Fig. 1  Higgs self-energies: The scalar sector of the SESM is characterized by the presence of two Higgs bosons, h and H, and loop corrections lead to a mixing between the two scalars. Therefore, the renormalized one-particle irreducible two-point function for two external scalar fields is not diagonal and can be split into a diagonal LO term plus a non-diagonal NLO contribution,

On-shell renormalization conditions
The functionsΣ ab are the renormalized self-energies (containing loop and counterterm contributions) with the fields a and b on the external legs. These can be cast in the form where the unrenormalized self-energies Σ ab , for a, b = h, H, contain loop contributions of the types shown in Fig. 2 The diagonal field renormalization constants are fixed by requiring that the (real parts of the) residues of the propagators at their respective poles are not changed by NLO corrections, i.e.
The renormalization constants δZ hh and δZ HH are then given by where Σ ′ (k 2 ) is the derivative of the unrenormalized self-energy with respect to the argument k 2 . Finally, to fix the mixing renormalization constants, we enforce the conditions that fields on their mass shells do not mix, i.e.
Using the expressions for the renormalized self-energies given in Eq. (4.5), the conditions (4.10) lead to In summary, the use of these on-shell conditions to fix δZ ij for i, j = h, H ensures that no wave function renormalization or Higgs mixing corrections for external Higgs states needs to be taken into account (these corrections are shifted to self-energy and vertex counterterms). Note also that this matrix field renormalization does not fix the mixing angle counterterm δα, although its UV divergences are connected to the ones in δZ ij via Eq. (3.11). In fact these relations will be used below to determine δα in the MS scheme, where δα only receives contributions from UV divergences.

Gauge-boson sector
The renormalization conditions for the gauge-boson sector of the SESM are identical to the ones used in the SM. The mass renormalization constants are fixed by imposing OS conditions on the W-and Z-boson masses M W and M Z , so that the renormalized squared masses correspond to the real parts of the locations of the propagator poles. 6 The field renormalization constants are fixed by requiring that the residues of OS propagators are not changed by NLO corrections; mixing renormalization constants are fixed in such a way that OS gauge bosons do not mix. The renormalization constants for the EW sector are then given by [62] are the transverse parts of the self-energies for generic gauge bosons V, V ′ . The renormalized electric charge e is defined as the electron−photon coupling in the Thomson limit of OS electrons with a photon at zero momentum transfer. The resulting charge renormalization constant reads [62]

Fermion sector
To fix the renormalization constants in the fermion sector, the procedure is identical to the SM case, described in detail in Ref. [62]. We require that the real parts of the locations of the poles of the fermion propagators correspond to the squares of the renormalized fermion masses, and that the residues of the fermion propagators do not receive loop corrections. Setting the CKM matrix to the unit matrix, the renormalization constants are given by where Σ f,L , Σ f,R , Σ f,S are, respectively, the left-handed, right-handed, and scalar parts of the fermion self-energy, as defined in Ref. [62]. The generalization to a non-trivial CKM matrix can be found in Refs. [62,65,66].

MS renormalization conditions
The mixing angle α and the coupling constant λ 12 still need to be fixed, but there is no obvious formulation of OS conditions, which are based on physical S-matrix elements, thereby avoiding any problems with gauge dependences. 7 In principle, these parameters could be extracted from the Higgs couplings to other particles and Higgs self-couplings, but we are far from having the precision required for such measurements. Also, requiring vanishing NLO contributions to a specific process could lead to artificially large contributions when computing other observables, as pointed out for other SM extensions [41,67].
Here, we present the MS renormalization conditions for α and λ 12 adopted in the two schemes considered in this work. Each scheme employs the same OS conditions for the other parameters and for the fields as described in Sect. 4.1. Imposing MS conditions (or conditions involving off-shell quantities) on mixing angles or couplings in spontaneously broken gauge theories is prone to introduce gauge dependences in the relations between physical observables and input parameters. Detailed discussions of this issue, which is intrinsically linked to the treatment of tadpole contributions, can, e.g., be found in Refs. [39][40][41][42][43][44]. In the following, we describe two different schemes, called "MS scheme" and "FJ scheme", which both renormalize α and λ 12 with MS conditions, but differ in the treatment of tadpole contributions. The former involves gauge dependences, while the latter does not.
In our MS scheme, the renormalization conditions for α and λ 12 are fixed using MS conditions, requiring UV finiteness for certain loop vertex functions, and demanding vanishing renormalized tadpoles. In spite of the issue of involving gauge dependences, this scheme is known to produce results that are rather stable with respect to variations of the renormalization scale. This was, e.g., observed in the THDM [41,43] and supersymmetric models [67]. In the second renormalization scheme, called here FJ scheme, we keep MS conditions for α and λ 12 , but change the tadpole treatment a la Fleischer and Jegerlehner [39] by setting bare tadpoles to zero consistently, which eliminates the gauge dependences. Technically, we follow the procedure described in Ref. [43] in detail for the THDM, i.e. we implement the FJ scheme by including appropriate finite terms in the renormalization constants obtained for α and λ 12 in 7 The use of physical S-matrix elements is crucial here to avoid gauge dependences. If instead renormalization conditions are imposed on off-shell Green functions or parts thereof (such as mixing self-energies, Green functions involving unphysical fields, etc.) at some momentum transfer, in general gauge dependences will result. Employing, for instance, the last two equations of Eq. (3.11) to derive δα from δZ hH and δZ Hh as given in Eq. (4.11) including UV-finite terms, leads to a gauge-dependent result, since Σ hH (M 2 h ) and Σ hH (M 2 H ) are not directly derived from S-matrix elements. the MS scheme. In applications to the THDM [41][42][43]56], it was observed that this scheme is prone to introduce large corrections that may also spoil the stability of predictions with respect to renormalization scale variations. To distinguish the two schemes we mark the renormalized parameters α and λ 12 and the corresponding renormalization constants with the superscripts MS and FJ if it is not clear from the context. The parameters in the two schemes are related by the coincidence of the respective bare parameters, which define the original Lagrangian, where we left implicit the dependence of the renormalization constants on the renormalized parameters. We will address the conversion between the two schemes in Sect. 7.1.1.

MS scheme (with vanishing renormalized tadpoles)
Mixing angle α: The renormalization constant for the mixing angle α can be determined from Higgs-boson self-energies using the relations given in Eq. (3.11). The last two relations yield Using the explicit expressions of Eq. (4.11) for the mixing renormalization constant, and recalling that MS renormalization constants contain only UV-divergent terms proportional to ∆ UV , the counterterm δα is given by Contributions induced by closed fermion loops are given by with c quark = 3, c lepton = 1, and the remaining bosonic contributions are where, to keep the notation compact, we introduced the functions F i , given by This counterterm has been computed in the 't Hooft−Feynman gauge and should be entirely used in this gauge.
h Figure 3: Diagram types contributing to the hhh vertex used for the renormalization of λ 12 .
Higgs self-coupling λ 12 : We fix the renormalization constant δλ MS 12 considering the loop corrections to the vertex function with three external light Higgs bosons, similar to the procedure pursued in Ref. [47]. Typical diagrams contributing to this vertex function are illustrated in Fig. 3. We require the one-loop renormalized vertex functionΓ hhh to be UV finite, This automatically renders all scalar three-and four-point vertex functions UV finite, and completes the set of renormalization conditions for the SESM. An explicit calculation of the counterterm yields for the contribution from closed fermion loops and for the bosonic contribution. Since λ 12 is a fundamental parameter of the original Lagrangian, the MS definition given above leads to a gauge-independent counterterm δλ MS 12 . We have checked our results on δλ MS 12 against a simpler derivation, which makes use of the fact that UV divergences in the CTs of dimensionless couplings are the same in the broken and unbroken phase of the theory. In the SESM, we can, thus, deduce δλ MS 12 in the unbroken phase where v 1 = v 2 = 0. In this phase, h 1 ≡ σ, and the coupling λ 12 only appears in the quartic couplings σσh 2 h 2 , σσχχ, and σσφ + φ − . At tree level, the σσφ + φ − vertex function is given by Figure 4: UV-divergent diagrams contributing to the σσφ + φ − vertex correction in the unbroken phase of the SESM, with V denoting the EW gauge bosons. For the first diagram, a crossed version exists as well.
and its MS CT reads Here we have used the fact that the σ field renormalization constant δZ σ = 0 in the unbroken phase, because σ appears only in quartic couplings, so that the σ self-energy is momentum independent. The MS field renormalization constant δZ MS Φ can be easily determined from the UV divergences in any of the Higgs-or Goldstone-boson self-energies, using Eq. (3.11). Only graphs with intermediate Goldstone-gauge-boson pairs or fermion-antifermion pairs contribute, yielding The UV-divergent diagrams contributing to the unrenormalized vertex function at one loop, Fig. 4. The corresponding divergences are easily calculated to where the order of terms follows the order of diagrams in Fig. 4. Demanding that the renormalized vertex functionΓ σσφ + φ − is UV finite, directly leads to (4.29) This is in agreement with the results (4.22) and (4.23) for δλ MS 12 given above, as can be checked by trading the couplings λ 1 and λ 2 for the chosen independent input parameters of the SESM with the help of Eq. (2.15).

FJ scheme (with vanishing bare tadpoles)
To obtain gauge-independent relations between observables and renormalized input parameters, we make use of the FJ scheme, proposed in Refs. [39,40] for the SM and applied to the THDM in Refs. [41][42][43]. In this scheme, the bare tadpoles t h,0 and t H,0 are set to zero, and any kind of reshuffling of tadpole terms is a mere question of taste, which does not change the results for observables. In principle, it is even possible to include explicit tadpole diagrams wherever they appear.
We have performed the FJ renormalization in two independent, but equivalent ways: Firstly, following the strategy proposed in Refs. [41,42], we have set the bare tadpole terms t h,0 and t H,0 to zero and omitted the introduction of tadpole renormalization constants δt h and δt H . Instead we have reintroduced the tadpole counterterms by shifting the Higgs fields according to where the additional finite terms depend on the tadpole contributions T h and T H . To compute these terms, we consider a variant in which the bare tadpole constants vanish, and the tadpole contributions are explicitly included in Green functions. Denoting the quantities defined in this scheme with a superscript "t", the same physical results are obtained using the counterterm where δα t can be obtained by setting to zero the tadpole renormalization constants δt h,H in δα MS and including tadpole diagrams in the related Green functions. For consistency, ∆α t contains also the finite terms coming from tadpole diagrams, otherwise the new terms could not be compensated by the tadpole contributions occurring elsewhere, leading to different renormalized amplitudes. The FJ renormalization scheme in the "t-variant" is obtained by reducing relation (4.32) to UV divergences only, since δα t,FJ has to be proportional to ∆ UV . Since δα MS is proportional to ∆ UV as well, we can write Taking the FJ version of Eq. (4.32) leads to where δα FJ is the counterterm we have to use in our counterterm Lagrangian to compute renormalized amplitudes in the FJ scheme. Finally, combining Eqs. (4.33) and (4.34), the relation between the α renormalization constant in the two schemes is given by The term ∆α t , according to Eq. (4.32), is the difference between δα t and δα MS , and is given by the tadpole contributions (that must be included in the "t-variant") to the self-energies used to define δα MS in Eq. (4.17), leading to where the superscript "t" in the self-energy Σ t,hH indicates that it is computed in the "tvariant", i.e. includes explicit tadpoles. Representing the unrenormalized tadpoles with black blobs, Eq. (4.36) leads to the expression with the factors for the hhH and hHH tree-level couplings given by 38) which are related to the couplings of Eq. (2.18) by c hhH = −eC hhH /2 and c hHH = −eC hHH /2. Therefore, in order to reproduce the result in the FJ scheme in the framework of our MS scheme, where we use vanishing renormalized tadpoles, we use the counterterm where the finite part of the last term is obtained by dropping the contributions proportional to ∆ UV from the expression (4.37). When computing a physical observable the use of this renormalization constant ensures a gauge-independent result.
5 Predictions for h → WW/ZZ → 4f in the SESM with the Monte Carlo program ÈÖÓÔ Ý

Features of ÈÖÓÔ Ý
The program Prophecy4f (Proper description of the Higgs decay into 4 f ermions) [51][52][53] is a Monte Carlo generator for the computation of any partial width for the decay of the Higgs boson into four light fermions at NLO, including both EW and QCD corrections. The generator can be used to produce differential distributions for any leptonic and semi-leptonic final state, as well as unweighted events for the leptonic final states. The first versions of Prophecy4f dealt with the decay of a SM Higgs boson and supported the presence of a fourth generation of massive fermions [55]. Recently, the program has been extended to allow for the same calculations in THDMs [43,56].
In the implementation, the final-state fermions are considered to be massless, but the physical mass values are kept in closed fermion loops which contribute to the virtual corrections. In the considered massless limit, the results are the same for final-state fermions of different generations (given that the same diagrams contribute), so that only the 19 independent final states reported in Table 1 need to be considered. In the table, these are classified by the intermediate gauge bosons appearing in the LO matrix element of the corresponding decay and by the number of lepton pairs in the final state. The W-and Z-boson resonances are treated in the complex-mass scheme [68][69][70], and the vector bosons are kept off-shell, so that the results have NLO accuracy both in resonant and non-resonant phase-space regions. The proper inclusion of e − e + uū (6) uūss (4) e − e + dd (9) neutral current with interference (12) udsc (2) charged and neutral current ν e e + e −ν e (3) uddū (2)  [71], which makes use of dimensional regularization to handle the UV divergences. Infrared (IR) divergences are regulated using small masses for the final-state fermions as well as for the emitted photon or gluon, and the divergences are canceled between virtual and real corrections using some slicing or the dipole-subtraction method [72][73][74].
The phase-space integral is performed by the adaptive algorithm implemented in the original Prophecy4f version, which evaluates the integrand at pseudo-random phase-space points, adapting iteratively the selection of channels in order to provide a better convergence.
Computing the widths for the decays of the Higgs boson into all the possible final states listed in Table 1 allows to get the total width for the inclusive decay of the Higgs boson into four fermions, Γ h→4f . The width Γ h→4f is the sum over the decay widths for the 19 independent final states, each of them weighted with the corresponding multiplicity given in Table 1.
In order to define a width for the decay of the Higgs boson into a pair of W or Z bosons, it is possible to separate contributions to Γ h→4f for which, in the LO matrix element, the intermediate vector bosons are two W or Z bosons. If both WW and ZZ are possible intermediate final states at LO, the WW and ZZ decay parts are defined by formally taking the two respective fermion−antifermion pairs of the W-or Z-boson decays from different generations. This procedure attributes all contributions to WW or ZZ channels except for terms that are interferences of WW-and ZZ-mediated contributions or corrections thereof. The sum of these interferences is denoted by Γ WW/ZZ−int (see also Refs. [75,76]), Note that interference contributions between ZZ channels with different fermion-number flow are included in Γ h→ZZ→4f . As a trivial example, consider the decay into ν e e + µ −ν µ , for which the LO process is entirely mediated by two W bosons, On the other hand, the leptonic final state ν e e + e −ν e contributes to all three parts of Eq. (5.1), Figure 5: Charged-and neutral-current LO diagrams contributing to the process h → 4f . The primed fermions f ′ and F ′ stand for the isospin partners of f and F , respectively. For F = f ′ in the final state, both contributions must be taken into account. For F = f the neutral-current diagram on the right-hand side appears in a second version with the f and F lines (or thef andF lines) interchanged.
Following this procedure for all four-fermion final states leads to the definition of Γ h→4f into WW-and ZZ-mediated parts and corresponding interference, 5.2 Details on the calculation of the process h → WW/ZZ → 4f at NLO

Leading order
At the Born level, in the massless limit for the final-state fermions, the decay h → 4f is mediated by a pair of (off-shell) gauge bosons, each of them decaying into two fermions. The contributions to the matrix element for the generic process are given by the Feynman diagrams reported in Fig. 5. Compared to the SM case, there are no additional diagrams, and the matrix element for the LO process is simply rescaled by a c α factor, so that LO predictions for the decay widths in the SESM can be easily obtained rescaling the SM results by a factor c 2 α . Depending on the fermions in the final state, the tree-level process either involves only the first diagram of Fig. 5 ("charged current"), only the second diagram ("neutral current"), both diagrams ("charged and neutral current", for f = F ′ ), or two diagrams of the second kind ("neutral current with interference", for f = F ).

Virtual corrections
Moving beyond the LO computation, loop and real-emission contributions must be taken into account. The one-loop virtual corrections to the decay process h → 4f in the SESM receive contributions from self-energy, vertex, box, and pentagon diagrams, as well as from counterterms in the self-energy and vertex corrections. These are very similar to the contributions arising in the SM case, which are described in detail in Refs. [51][52][53]. Indeed, the set of Feynman diagrams that contribute to the decay in the SESM is given by all the SM diagrams, supplemented by A survey of the generic diagrams contributing to the QCD matrix element of Eq. (5.6) is reported in Ref. [51].
EW loops: Comparing the SESM to the SM, the presence of the singlet has an impact on the EW corrections, giving rise to a higher number of loop diagrams and changing the analytic expressions of the SM-like contributions. For a list of the generic diagrams contributing to the EW matrix element, see Ref. [53]. Since, in these diagrams, no internal scalar lines appear in box and pentagon graphs, the heavy Higgs yields only additional self-energy and vertex diagrams of the type reported in Fig. 6. The computation of the EW loops can be performed with the standard machinery.

Real-emission corrections
At NLO, the final-state fermions can emit a photon or a gluon, so that it is necessary to include diagrams as depicted in Fig. 7. In the real-emission contributions the IR structure is the same as in the SM, and the extraction of the soft and collinear divergences appearing in the phase-space integration of the squared matrix elements of Eqs. (5.7) and (5.8) can be performed with the same methods used in the standard case [53]. Note that all LO and (real and virtual) QCD amplitudes are related to the corresponding SM counterparts by the factor c α , so that the relative QCD corrections to the partial widths (normalized to LO) are the same in the SESM and SM.

Complex-mass scheme
Within the complex-mass scheme [68][69][70], the renormalized masses of the W and Z bosons are replaced by the complex masses µ W and µ Z , defined via the real pole masses M W , M Z , and the decay widths Γ W , Γ Z , of the gauge bosons by The cosine of the weak mixing angle, which is defined by the ratio of the W-and Z-boson masses, is replaced by the complex quantity c W = µ W /µ Z . This relation ensures the gauge independence of NLO matrix elements in spite of the use of complex W and Z masses, whose imaginary parts result from a partial resummation of self-energy contributions. Even though the Higgs particles are unstable, we treat only the W and the Z bosons in the complex-mass scheme. Indeed, effects induced by a complex Higgs-boson mass are of the order Γ h /M h and, assuming that the light Higgs boson of the SESM has a small width (as it happens for the SM Higgs), these are negligible compared to the NLO contributions considered in this work. The heavy Higgs enters only loop diagrams, so that corrections from a complex mass are negligible, as long as Γ H ≪ M H . In the complex-mass scheme, the renormalization constants for the W-and Z-boson masses are complex to guarantee that µ 2 W and µ 2 Z correspond to the complex locations of the W and Zpropagator poles. This implies that the renormalization constant for the weak mixing angle becomes complex as well. The W -and Z-field renormalization constants are defined in the complex-mass scheme by self-energies (and not only by their real parts) which depend on the complex parameters, so that the field renormalization constants are complex. The electric charge renormalization constant, which depends on the (complex) field renormalization constants, becomes complex as well. Explicit definitions for the renormalization constants in the complex-mass scheme are given in Ref. [70] for the SM. In the SESM, the definitions of the additional renormalization constants δα and δλ 12 are not changed in the complex-mass scheme, but the constants are treated as complex quantities, since they are defined using two-and three-point loop functions which contain the complex W-and Z-boson masses and complex couplings.

G µ scheme
Adopting the so-called "G µ scheme", we use the Fermi constant G µ as input parameter and compute the electromagnetic coupling constant α em = e 2 /(4π) according to In this way, a large universal part of the O(α em ) corrections is absorbed into the LO prediction. More precisely, this choice absorbs the running of α em from zero-momentum transfer to the weak scale and the universal corrections to the ρ parameter into the lowest-order coupling α em /s 2 W . Following this procedure, to avoid double-counting, we have to subtract the EW corrections to muon decay from the explicit NLO contributions to the electric charge renormalization constant δZ e (see also Ref. [53]), where the renormalization constant δZ e is given by Eq. (4.13), and (∆r) 1-loop is the one-loop weak correction to the muon decay ∆r [62,77], but now calculated in the SESM, as, e.g., done in Ref. [78]. For consistency, both contributions are computed in the complex-mass scheme.
Using nevertheless the real value for α em defined in Eq. (5.10) is consistent at NLO.

Implementation into ÈÖÓÔ Ý
To take advantage of the capabilities of the original Prophecy4f version, we have modified the code in order to include the expressions for the SESM matrix elements described in Sect. 5.2. For the LO contributions, the QCD corrections, and the photonic real-emission contributions to the decay process h → 4f , this can be easily achieved by rescaling the SM Higgs couplings to vector bosons by the appropriate prefactor c α , according to Eqs. (5.5), (5.6), (5.7), and (5.8).
For the EW virtual corrections, we computed the matrix elements in two independent ways. In the first computation of the NLO matrix elements contributing to the decay h → 4f , we constructed a model file for the SESM, including all the one-loop counterterm vertices and the definitions for the renormalization constants, for the amplitude generator FeynArts [59]. To produce the model file we used the Mathematica package FeynRules [57,58]. FeynRules allowed us to get the Feynman rules from the SESM Lagrangian (including the vertices from the counterterm Lagrangian described in Sect. 3) and to generate the FeynArts model file in an automated way. Afterwards, we have added the definitions of the renormalization constants to the FeynArts model file as they are reported in Sect. 4, using the FormCalc format [60], both for the MS and the FJ renormalization schemes. The model file can be used to generate, to compute, and to simplify one-loop matrix elements with the packages FeynArts and FormCalc for (in principle) any process within the SESM. The model file has been tested by checking UV finiteness for many processes, both analytically and numerically, devoting special attention to the multi-scalar vertex functions, which involve the renormalization constants δα and δλ 12 . We adapted the model file to the demands of Prophecy4f, using complex masses for the gauge bosons and keeping the full mass dependence in the closed fermion loops, and used it to generate the Fortran routines for the computation of the virtual matrix elements contributing to the decay h → 4f . Finally, we have incorporated the Fortran code in Prophecy4f.
In the second calculation, we generated the amplitudes using a tree-level FeynArts 1 [79] model file, and we inserted the counterterms by hand and processed them further with in-house Mathematica routines. The results from both calculations are UV-and IR-finite and in good mutual numerical agreement.

Input parameters and benchmark scenarios
In this section we fix the input parameters used to derive our numerical results. In Sect. 6.1 we present the SM input parameter set and in Sect. 6.2 we discuss how the parameter space of the theory is constrained by the requirements of vacuum stability and perturbativity of the couplings. In Sect. 6.3 we define the benchmark scenarios used for the numerical evaluations.

SM parameters
We identify the light Higgs boson h with the known Higgs particle and set M h = 125.1 GeV, (6.1) in agreement with the mass value measured by ATLAS and CMS [80]. The numerical values of the other parameters are fixed according to the recommendations of the LHC Higgs Cross Section Working Group (HXSWG) [3], mostly based on Ref. [81]. The Fermi and the strong coupling constants are We simply take α s at the scale of the Z-boson mass, i.e. we do not change the QCD renormalization scale in the scale variations discussed below, because it merely leads to changes at next-to-next-to-leading order, which are part of the residual theoretical uncertainty from missing higher orders. The OS gauge-bosons masses and widths and the fermion masses are In the numerical analysis, we use the W-and Z-boson masses obtained from Eq. (6.4). The decay widths Γ W and Γ Z are calculated from the given experimental input, taking into account O(α em ) corrections and using real masses. We do not use the pole widths of Eq. (6.4), but we compute Γ V at NLO, in order to ensure that the effective W/Z branching ratios add up to one in the sum over all decay channels. In this step, we neglect effects due to the presence of the singlet; in principle, it contributes to the NLO corrections to the W and Z widths, but for the small α value we consider, the effect is negligible. As in the original SM version of Prophecy4f the full dependence on the fermion masses given in Eq. (6.3) is kept in corrections induced by closed fermion loops, while external (light) fermions are treated in the massless limit. Since quark mixing to the third generation as well as the differences in the (internal) light-quark masses are negligible, the CKM matrix drops out in the calculation of the inclusive (flavour-summed) width Γ h→4f . We, thus, set the CKM matrix to the unit matrix in the following.

Constraints on BSM parameters
Even without taking into account the data collected from the experiments, the parameter space of the SESM is limited by theoretical constraints [33,35,37,45,46,82]. Before choosing the input values for the free parameters of the theory, it is worth recalling these constraints and how the free parameters can fulfill such conditions.
Perturbativity of the couplings: The scalar couplings of the SESM must not exceed a certain value, so that the perturbative approach used in the calculations remains valid. Thus, we require that the contributions from the coupling constants λ i to the coefficients of the quartic coupling terms in the Higgs potential, respect some limit O(|λ i |/π) 1. The following choice is made in order to replicate the results of Ref. [37], where a similar analysis was performed using a different input parameter set and different conventions. The conditions from there translate into the bounds These values are meant to be rough estimates that are used to show where perturbativity problems can arise, rather than sharp boundaries on the allowed values.
Vacuum stability: As discussed in Sect. 2.1, vacuum stability at LO is guaranteed by the conditions given in Eq. (2.16).
In Fig. 8 we show the effects of the requirements of perturbativity (6.6) and vacuum stability (2.16) on the input parameter space for different heavy Higgs masses in the range M H = 200−800 GeV. The allowed region is restricted to the white area, and it is possible to see how the condition on λ 1 presented in Eq. (6.6) plays an important role both for negative and positive λ 12 values. The vacuum stability condition (2.16) has an impact only on negative λ 12 values, ruling out a large part of space that is not excluded by perturbativity requirements. Green and yellow areas indicate regions where one can expect perturbativity problems, but should not be intended as sharp-cut regions. For the M H values considered here (200, 400, 600 and 800 GeV), the perturbativity of λ 12 does not affect the parameter space, but it becomes relevant for higher M H values. The perturbativity constraint on the coupling λ 2 is irrelevant in the considered regions.

Benchmark scenarios
For the numerical analysis we consider some of the benchmark scenarios proposed in Ref. [3], which were originally suggested in Ref. [37], adapting the input values to our needs. In Ref. [37], different values for the mass M H are considered (both lighter and heavier than M h ), and for each mass the mixing angle α is fixed to the maximal allowed value. Moreover, for mass values M H ≥ 2M h (i.e. when the H → hh decay is kinematically allowed), two values are proposed for tan β ≡ v 2 /v 1 , corresponding to the maximal and the minimal branching ratios for the H → hh decay. Among these possibilities, we only consider scenarios in which M H > M h (since the other possibility is phenomenologically disfavoured) and vary the heavy Higgs mass in the interval 200−800 GeV with 200 GeV steps. When two different tan β values are proposed, we consider the average. Since tan β enters in our calculation only at NLO and the two proposed values are always quite close, we expect negligible differences due to this choice.
We have to convert the numerical values of the input parameters given in Ref. [37] to our conventions. The SESM Higgs Lagrangian used here, as given in Eq. (2.2), is equivalent to the one given in Ref. [37] using the following substitutions, where the label "ref" indicates the parameters used in Ref. [37], in which the numerical input is given in terms of M H , α, and tan β ref .  taken over directly. The scalar coupling λ 12 , which is a free parameter in our conventions, can be obtained from tan β ref using the relations We convert the benchmark points using Eqs. (6.7) and (6.8), rounding the λ 12 values to two decimal digits. The input values, in our convention, for the scenarios considered in our analysis are reported in Table 2  In the following we will make use of these scenarios in each of the renormalization schemes proposed in this paper.

Numerical analysis
In the following, we present the numerical results relevant for the decay h → 4f of the light Higgs boson of the SESM. Starting from benchmark scenario BHM200 + , we show the effects of the conversion of the input variables between the two renormalization schemes presented in Sect. 4. Then we investigate the scale dependence of the parameters α and λ 12 , which are defined by MS renormalization conditions, by solving numerically the corresponding renormalization group equations (RGEs). Afterwards, we present the results for the decay width Γ h→4f computed at different renormalization scales and show the deviations from the SM results as a function of the mixing angle. The same analysis is presented for benchmark scenario BHM600, while results for the scenarios BHM200 − and BHM400 are reported, respectively, in Appendices A and B. Finally, we show some differential distributions, comparing the results in the SM with the ones in the benchmark scenarios of Table 2. 7.1 BHM200 +

Scheme conversion
When computing a physical observable at NLO accuracy, starting from a set of input parameters, it is crucial to realize that the input values correspond to a specific renormalization scheme adopted in the calculation. This becomes even more important when comparing NLO results for the same observable obtained using different renormalization schemes. In different schemes, the same numerical values for the input parameters represent different physical scenarios and, in order to have a sensible comparison of predictions for an observable in a given scenario, a proper conversion of the input parameters between the schemes is required.
In general, defining N renormalized parameters p i in two renormalization schemes, denoted, respectively, by p (1) i and p (2) i , the relation between them is given by the solution of the following system of equations, where the connection between the parameters in the two schemes is given by the bare parameters p i,0 , which are independent of the renormalization scheme. In our particular case, converting the input values from the MS to the FJ scheme is quite simple, since, apart from the mixing angle α, all the other input parameters of the SESM have the same definition in the two schemes. Ignoring effects beyond NLO, the input parameters p i = α are defined by identical renormalization conditions in the two schemes, i.e.
To solve the equation and find the relation between α MS and α FJ , we adopt two strategies. In the first approach we linearize Eq. (7.3) and obtain Since our computations are performed at NLO, the O(α 2 em ) term in Eq. (7.4) can be neglected. An analogous procedure can be applied to determine α MS when α FJ is given as input. Using this method, converting an input value for the mixing angle from one scheme to the other and repeating the procedure to go back to the initial scheme, the final numerical result for α will change by contributions that are formally beyond NLO.
In the second approach, we solve Eq. (7.3) numerically, in order to keep the contributions of O(α 2 em ), which can become relevant for large counterterms or small tree-level values. Using this method, converting α to the other scheme and back, does not change the value of α. In the following results we use, as much as possible, the second method, i.e. we include the O(α 2 em ) terms.
In Sect. 4.2 we have derived the counterterm δα in the two schemes, which differs by finite contributions, where the tadpoles T h , T H , and the coupling factors C hhH , C hHH in the last term depend on α FJ . Using these expressions, it is straightforward to get the conversion from the FJ to the MS scheme, while the conversion from MS to FJ requires a numerical solution of Eq. (7.6) for α FJ , which appears also in the ∆α t term. In Fig. 9, we show the results for the conversion of the sine of the mixing angle, s α , between the two schemes, both for the full solution of Eq. using the linearized solution (7.4). 8 The curves on the right sides inside the plots are obtained fixing the mass M H of the heavy Higgs boson and the coupling λ 12 according to their values in the scenario BHM200 + , reported in Table 2. On the left sides, similar curves show the conversion effects for negative s α values. For consistency, we adjust the sign of λ 12 so that sgn(s α ) = sgn(λ 12 ) (for the input s α ) and Eq. (2.12) is not violated. The renormalization scale is fixed to the mass of the light Higgs boson, µ r = M h ; the motivation for this choice will become clear in Sect. 7.1.3. The dark-gray shaded areas in the plots mark the values of s α for which the perturbativity constraint (6.6) on λ 1 is violated; from the last line of Eq. (2.15) it is easily seen that λ 1 necessarily violates its perturbativity bound for s α → 0, since we keep M h , M H , v 2 , λ 12 fixed. The light-gray shaded areas denote regions where the sign of s α is flipped by the conversion and becomes inconsistent with the sign of the considered λ 12 . The conversion effects in the perturbative regions are small: The red line is, in general, very close to the dashed diagonal line, which corresponds to the absence of any conversion effect (i.e. s MS α = s FJ α ), and the linearized solution reproduces the full conversion very well. Large effects (and deviations between full and linearized solutions) are only observed when approaching the non-perturbative regime, corresponding to small values of the mixing angle. In both plots of Fig. 9, a slight asymmetry can be observed between positive and negative s α values, due to the different NLO contributions obtained by changing the sign of the input values for s α and λ 12 .

Running of s α and λ 12
Since we have defined the parameters α and λ 12 by MS renormalization conditions, they depend on an unphysical renormalization scale µ r . The dependence on this scale is governed by the RGEs ∂ ∂ ln µ 2 r α µ 2 r = β α µ 2 r , ∂ ∂ ln µ 2 r λ 12 µ 2 r = β λ 12 µ 2 r , (7.7) 8 If Eq. (7.3) is used for the conversion, the corresponding curves in the two plots are related by a simple reflection about the diagonal s MS α ≡ s FJ α (apart from the different truncation of the curves in the non-perturbative region). The reflection symmetry is not there in the linearized version (7.4), but broken by effects beyond NLO. Note also that both versions coincide on the r.h.s., because δα MS does not contain finite contributions (UV divergent terms are canceled analytically). where the β α and β λ 12 functions can be extracted from the expressions of the counterterms δα and δλ 12 , taking the coefficients of the UV divergence ∆ UV . These functions are different for the two considered renormalization schemes: For the MS scheme, the β functions can be obtained considering the following derivatives with respect to the UV divergence, where the counterterm δα MS is given in Eqs. (4.18) and (4.19), and δλ MS 12 in Eqs. (4.22) and (4.23). For the FJ renormalization scheme, also the UV contributions due to the tadpoles must be taken into account, leading to the β functions where β λ 12 is not changed due to the fact that δλ 12 is the same in the two schemes. The RGEs are coupled differential equations for which, in general, an analytical solution is not possible. We solve the equations numerically, using a Runge−Kutta algorithm, obtaining the scale dependence for the sine of the mixing angle, s α , and the coupling λ 12 , as shown in Fig. 10. These results are obtained fixing the parameters M H , α, and λ 12 according to the values reported in Table 2 for benchmark scenario BHM200 + at the scale µ 0 = M h , and changing the scale µ r in the range 50−350 GeV, using the β functions for the two schemes. Since the purpose here is the assessment of the scale dependence of the MS parameters, no conversion between schemes is applied on the input values. In the two schemes, the scale dependence of s α shows a completely different behaviour, while the running of λ 12 displays the same trend in the two schemes. As we will discuss below, the scale dependence of the mixing angle has a big impact on the scale variation of the decay width Γ h→4f .

Scale dependence of the inclusive decay width Γ h→4f
In order to assign a sensible value to the renormalization scale µ r , we study the impact of the scale choice on the results for Γ h→4f , the inclusive decay width of the light Higgs into four fermions. Taking Fig. 10). Comparing LO and NLO results, it is evident how the inclusion of loop contributions drastically reduces the scale dependence, as well as the scheme dependence: The solid lines are not only much flatter than the dashed lines, proving a reduced scale dependence, but solid curves of different colours are much closer to each other than the corresponding dashed curves. This is expected, since results obtained within different renormalization schemes should be equal up to higher-order contributions, provided that the input parameters are properly converted. The NLO predictions, in particular, perfectly coincide at the central scale µ r = µ 0 = M h . Quantifying the scale dependence by the change in Γ h→4f obtained by varying the scale µ r by factors of two up and down (µ r = 2µ 0 and µ r = µ 0 /2), we observe a reduction from ∼ 3−4% at LO to < ∼ 0.5% at NLO. The scheme dependence, defined by the relative difference between the results in the MS and FJ schemes at the central scale, on the other hand, reduces from ∼ 1% at LO to < ∼ 0.1% at NLO.
As discussed in Refs. [43,56] for NLO predictions of Γ h→4f in the THDM, when the computation involves multiple mass scales in the loops, it is not clear a priori that the central scale µ 0 = M h is an appropriate choice. In principle, this applies also to the SESM, where the heavy Higgs boson appears in loop diagrams. However, it is evident from Fig. 11 how the scale dependence is minimized for values around the light Higgs mass. On the other hand, using the alternative scale µ 0 = (M h + M H )/2 analogous to the scale choice advocated in Ref. [43] for the THDM, would not make a big difference for the scenario BHM200 + , due to the relatively small value for the heavy Higgs mass, fixed at M H = 200 GeV. As discussed later in the  Figure 12: Dependence of the decay width Γ h→4f at LO and NLO EW+QCD with respect to the variation of s α . Free parameters other than s α are fixed according to benchmark scenario BHM200 + (see Table 2). On the left (right) the input parameters are defined in the MS (FJ) scheme and converted to the FJ (MS) scheme (both for LO and NLO predictions). The vertical dashed line signals that for s α 0.06 perturbativity problems might arise.
section, repeating the same analysis on the scale variation with higher M H values, we will see how the scale µ 0 = M h is better than the alternative scale given by the arithmetic mean of the Higgs-boson masses. Note that finding optimal scales, to some extent, is empirical, i.e. the scale dependence should be investigated whenever qualitatively new scenarios are considered.

Mixing-angle dependence of the inclusive decay width Γ h→4f
Among the free parameters of the SESM, the mixing angle α plays the central role in the computation of the decay width Γ h→4f . Its value affects already the LO result, while the heavy Higgs mass M H and the coupling λ 12 enter only the NLO decay amplitudes. For this reason, we compute the decay width varying s α in the range 0.01−0.3, keeping M H and λ 12 fixed according to the values for benchmark scenario BHM200 + . 9 We consider both MS and FJ as input schemes and compute the decay width in the two schemes, using the renormalization scale µ r = M h . The results for Γ h→4f are reported in Fig. 12, where we also show the SM value of the decay width identifying h with the SM Higgs boson. Dashed and solid lines denote, respectively, the LO and the NLO results, where the latter include both EW and QCD corrections. The dashed vertical line indicates the minimal s α value for which the perturbativity conditions of Eq. (6.6) are satisfied. Differences between the LO results in the two schemes (within the same plot) are due to the scheme conversion, which is done at NLO. Comparing the LO results, it is possible to observe the suppression with respect to the SM given by the c 2 α factor, coming from the square of Eq. (5.5). The proportionality to c 2 α = 1 − s 2 α is exact if no conversion of the input is done, i.e. on the left (right) side for the LO MS (FJ) curve. The NLO contributions modify the s α dependence, so that for s α ∼ 0.08 the NLO decay width is the same as in the SM, and in general the difference between SESM and SM are smaller for the NLO decay width. Note also the reduction of the scheme dependence at NLO, visible by the fact that MS and FJ curves practically lie on top of each other.  Figure 13: s α dependence of the relative EW+QCD NLO corrections to the decay width Γ h→4f . Free parameters other than s α are fixed according to benchmark scenario BHM200 + (see Table 2). The vertical dashed line signals that for s α 0.06 perturbativity problems might arise.  Figure 14: Deviation of the LO and NLO decay widths Γ h→4f in the SESM from the SM values, as a function of s α . Free parameters other than s α are fixed according to benchmark scenario BHM200 + (see Table 2). The vertical dashed line signals that for s α 0.06 perturbativity problems might arise.
In Fig. 13 we show the relative corrections δ NLO to the inclusive h → 4f decay width, defined by where the NLO result includes both EW (δ EW ) and QCD corrections (δ QCD ). The relative QCD corrections do not depend on the mixing angle and are equal to the SM case, providing an offset of about 5%. The relative corrections in the SESM are, for the scenario considered here, bigger than the relative corrections in the SM case, somewhat compensating the LO suppression factor c 2 α mentioned above. For s α = 0.29, as defined for the scenario BHM200 + , the relative corrections are 9.6% in the FJ scheme and 8.6% in the MS scheme. Figure 13 illustrates that in the two schemes contributions (related to the tadpole terms) are shared differently between LO and NLO parts. Recall that we have seen in Fig. 12 how the NLO decay widths are in good agreement in the two renormalization schemes, independent of the s α value.
In Fig. 14 we compare the SESM result with the SM prediction, defining the relative deviation by Again, it is possible to see how the negative deviations induced by the c 2 α factor are somewhat reduced at NLO due to the positive loop contributions, which in the SESM are bigger than in the SM. For s α = 0.29, corresponding to the scenario BHM200 + , the deviations from the SM  Table 3: Partial widths for scenario BHM200 + in the FJ renormalization scheme. The integration errors are given in parentheses.
are ∆ SM ∼ 7−8%. The difference between the results shown in the two panels, which correspond to s α values in the two different renormalization schemes, is due to the fact that, changing the input scheme, the same numerical values correspond to different physical scenarios.

Decay widths for individual four-fermion final states
In Table 3 we compile our results on the widths for the decays h → 4f into the various different final states for benchmark scenario BHM200 + , computed in the FJ scheme. The contributions Γ h→WW→4f , Γ h→ZZ→4f , and Γ WW/ZZ−int are calculated according to Eq. (5.4), and the total decay width Γ h→4f using Eq. (5.1). In Table 3, we also show the relative EW and QCD corrections, δ EW and δ QCD , and in the last two columns the deviation ∆ SM from the SM both at LO and NLO. For all the quantities, we report the integration uncertainty in parentheses.
To determine the errors of the decay widths to WW, ZZ, WW/ZZ interference and of the total width Γ h→4f , we apply the standard error propagation to Eqs. (5.1) and (5.4), making use of the integration uncertainties for each single final state.
The main contribution to the total h → 4f decay width originates from the charge-current final states, while the neutral-current processes have a smaller impact, and the WW/ZZ interference gives a very small negative contribution. The EW corrections to the WW contributions are about 5%, and lead to a similar value for the inclusive decay h → 4f . The EW corrections to neutral-current final states range from 0 to 5%, depending on the flavour of the final-state fermions. The QCD corrections are mostly due to the corrections to the decays W/Z → qq, and amount to α s /π for each quark pair in the final state. Exceptions are the final states uūuū and dddd, where interference contributions from two different topologies of the ZZ channel occur, and the QCD corrections to these final states are only about 4%. The SM deviation, at LO, comes from the c α rescaling factor of the hV V coupling with respect to the SM coupling, and is equal to c 2 α − 1 = −0.0841 in the considered scenario. As already observed in the previous section, at NLO the deviations from the SM are about 1.5% smaller.
We have computed the same quantities as in Table 3 in the MS scheme, observing somewhat smaller values for the EW corrections, since in the two schemes contributions are shared differently between the LO and the NLO (as observed also in Fig. 13). Moreover, using the same numerical input in the MS input scheme leads to NLO deviations from the SM about 1% higher, since the same numerical input corresponds to a slightly different physical scenario. In total, the MS results follow the same qualitative pattern as in the FJ scheme, and are not reported here.

Scheme conversion
In Fig. 15 we show the conversion of the input value for s α from the MS to the FJ scheme (left panel) and vice versa (right panel). The values of M H and |λ 12 | are given in the scenario BHM600 as stated in Table 2. For consistency, for negative s α input values we consider a negative Higgs self-coupling λ 12 . The renormalization scale µ r is again chosen to be equal to the light Higgs mass M h . For the considered values of the input parameters, the perturbativity constraint on λ 1 of Eq. (6.6) is violated for |s α | 0.04. The dark-gray shaded area corresponds to the region where the perturbativity condition breaks down. The light-gray areas denote regions where the vacuum stability condition (2.16) is violated and where the sign of s α becomes different from the sign of λ 12 by effect of the conversion. As discussed in Sect. 7.1.1, the conversion from the FJ to the MS scheme can be computed straightforwardly from the matching condition (7.3), while to compute the inverse conversion the matching condition has to be solved numerically. Alternatively, the linear approximation (7.4) may be used. In the left panel of Fig. 15 we use the linear approximation (shown in blue) in the non-perturbative region, where the numerical inversion does not provide a solution. The red lines correspond to the results obtained from Asymmetries of the plots are due to the different sign used for λ 12 , so that sgn(s α ) = sgn(λ 12 ).

Running of s α and λ 12
The solution of the RGEs for the MS parameters s α and λ 12 is shown in Fig. 16 for scenario BHM600. The solutions are again obtained numerically by a Runge−Kutta algorithm. The values for the BSM parameters s α , M H , and λ 12 are fixed according to scenario BHM600 at the scale µ 0 = M h and evolved to the renormalization scale µ r in the range µ r = 50−350 GeV. We report the results both for the MS and the FJ schemes, where we can observe a different behaviour in the running of s α : The sine of the mixing angle increases with the scale in the FJ scheme, while it slowly decreases in the MS scheme. The running of the coupling λ 12 , in the considered range, is almost identical in the two schemes. Compared to the scenario BHM200 + , the running of s α in the MS scheme is strongly reduced.

Scale dependence of the inclusive decay width Γ h→4f
In Fig. 17 we present the results for the inclusive decay width Γ h→4f for scenario BHM600 at scales µ r in the range 50−350 GeV. The results obtained using the MS (magenta) and the FJ (blue) input schemes are reported, respectively, on the left and the right panels of the figure.  Figure 18: Dependence of the decay width Γ h→4f at LO and NLO EW+QCD with respect to the variation of s α . Free parameters other than s α are fixed according to benchmark scenario BHM600 (see Table 2). On the left (right) the input parameters are defined in the MS (FJ) scheme and converted to the FJ (MS) scheme (both for LO and NLO predictions). The vertical dashed line signals that for s α 0.04 perturbativity problems might arise.

Mixing-angle dependence of the inclusive decay width Γ h→4f
In Figs. 18 to 20 we present, respectively, the decay width Γ h→4f , the relative corrections to the decay width, and the deviations with respect to the SM result as a function of the parameter s α using our default scale choice µ r = µ 0 = M h . The heavy Higgs mass and the coupling λ 12 are fixed according to scenario BHM600, as defined in  Figure 19: s α dependence of the relative EW+QCD NLO corrections to the decay width Γ h→4f . Free parameters other than s α are fixed according to benchmark scenario BHM600 (see Table 2). The vertical dashed line signals that for s α 0.04 perturbativity problems might arise.  Table 2). The vertical dashed line signals that for s α 0.04 perturbativity problems might arise. s α value for which the perturbativity constraints of Eq. (6.6) are fulfilled. The results reported in the left (right) panel are obtained with input parameters defined in MS (FJ) scheme with a proper parameter conversion at NLO applied if the input is not directly given in the scheme used in the calculation.
In the LO results of Fig. 18, the c 2 α dependence of the decay width in the SESM is manifest, where the proportionality to c 2 α is exact if no scheme conversion is involved. For small mixing angles, where c α → 1, the deviations from the SM are due to the conversion effects. As observed in the other scenarios, the inclusion of NLO corrections slightly compensates the c 2 α suppression for small s α values. Moreover, we observe that the two schemes are in much better agreement after the inclusion of NLO corrections, i.e. the inclusion of NLO corrections reduces the renormalization scheme dependence drastically.
The relative NLO corrections, as defined in Eq. (7.10), are shown in Fig. 19 including both EW and QCD contributions (the latter are independent of s α and about 5%). In the considered s α range the relative corrections amount to 7−9% (slightly depending on the input scheme used) and the difference with the SM relative corrections |δ SESM − δ SM | does not exceed the 1% level in the considered region.
In Fig. 20 we illustrate the deviations from the SM, defined by Eq. (7.11). It is possible to see how the NLO results converge nicely in the two schemes. Note that the large schemedependence due to missing higher-orders for small s α values occurs only in the non-perturbative

Decay widths for individual four-fermion final states
In Table 4 we compile the results for the h → 4f decay widths for each final state listed in Table 1, together with the widths for the decays into two vector bosons and the inclusive decay width Γ h→4f . For each decay channel, we report the NLO result, the EW and QCD relative corrections, and the deviations from the SM result. The latter are reported both for LO and NLO. The qualitative picture is basically the same as for the scenarios BHM200 + , with slightly different values for the EW corrections, which are between −2% and 3%, i.e. a bit smaller than in the other cases with smaller M H . The deviation from the SM is about −5%, in the FJ scheme, at LO and NLO. Using the MS input scheme, the deviations ∆ NLO SM are only slightly different, around −4.5%, and display a similar pattern.

Differential distributions
Footprints of potential BSM physics often can be found by looking into the shapes of differential distributions, even if integrated results do not deviate from the SM predictions significantly. In the following, we discuss some of the distributions produced with Prophecy4f for the SESM, reporting results both for charged-and neutral-current final states. The generator provides invariant-mass and angular distributions for leptonic and semi-leptonic final states, while dis-   Figure 21: Invariant-mass and angular distributions for the neutral-current decay into leptons h → µ − µ + e − e + , in the FJ renormalization scheme for the various SESM scenarios.
tributions for fully hadronic final states are not interesting, since they are not experimentally accessible. A detailed survey of distributions in the SM can be found in Refs. [52][53][54] for fully leptonic final states and in Ref. [51] for the semi-leptonic case. There, the treatment of the final-state radiation and the photon recombination for nearly collinear fermion−photon pairs are also discussed in detail. In the SESM, relative corrections induced by final-state radiation are identical to the SM case and thus not discussed in greater detail here.

Leptonic final states
We consider the fully leptonic final states µ − µ + e − e + and ν µ µ + e − ν e , which involve, respectively, intermediate (off-shell) WW and ZZ states. The distributions discussed in the following have been computed in the FJ renormalization scheme, the results obtained in the MS scheme show the same features and are not reported here. All the distributions are generated using the renormalization scale µ r = M h .
For the neutral-current final state µ − µ + e − e + , the left panel of Fig. 21 shows the NLO distribution for the invariant mass of the muon pair around the Z-boson resonance, both for the SM and for the considered SESM scenarios. In general, the invariant mass of a fermionic pair is defined by where k a and k b are the four-momenta of the fermion f a and the anti-fermionf b and, in case of photon recombination, the photon momentum is added to the momentum of the fermion or anti-fermion. The distributions obtained for the SESM display the Z-boson resonance peak and the "radiative tail" observed in the standard case: The real-emission contributions lead to positive corrections for invariant masses below the Z-boson peak. This is explained by the fact that outgoing leptons lose momentum when radiating a photon, so that these events are shifted towards lower masses. This can be seen in the lower panel of the plot, where the relative NLO corrections δ NLO are shown bin-by-bin. The shapes of the distributions in the SESM are the   Figure 22: Invariant-mass and angular distributions for the charged-current decay into leptons h → ν µ µ + e − ν e , in the FJ renormalization scheme for the various SESM scenarios.
same as in the SM, and the only difference is given by an offset, which depends on the SESM scenario. The difference between the NLO distribution in the SM and in the SESM equals, for each considered scenario, the quantity ∆ NLO SM obtained for the corresponding integrated result (see Tables 3 and 4 for BHM200 + and BHM600). We observe this pattern in all the generated distributions. The right panel of Fig. 21 shows, for the same final state, the differential decay width with respect to the angle φ, which is defined as the angle between the decay planes of the two intermediate Z bosons in the Higgs rest frame, with the sign convention (7.14) The distribution resembles a cos(2φ) oscillation with some constant offset and can be used to determine the parity of the Higgs boson and to set bounds on BSM couplings to EW gauge bosons of the decaying scalar (see Refs. [83][84][85][86][87][88][89][90]). As observed for the invariant-mass distributions, the shape of the distributions in the SESM scenarios is the same as in the SM, and the NLO relative corrections differ just by a constant offset which is equivalent to ∆ NLO SM observed for the integrated widths. Figure 22 shows differential distributions for the charged-current final state ν µ µ + e − ν e for the SM and the SESM. The invariant-mass distribution for M νµµ , shown in the left panel of the figure, is not experimentally accessible, but is interesting from the theoretical point of view. The distribution shows the W-boson resonance and the radiative tail already described for the µ − µ + e − e + final state, with an enhancement below the W-boson peak driven by the real photon emission. The SESM does not induce any additional distortion on top of the SM shape. The corrections and the differences between the SM and the SESM scenarios match always the values   cos φ µe,T = k µ,T · k e,T |k µ,T | |k e,T | , sgn(sin φ T ) = sgn{e z · (k µ,T × k e,T )}, (7.15) where k i,T are the projections of the lepton momenta onto the plane orthogonal to the unit vector e z , which denotes the beam direction of the Higgs production process. The distributions, in the SESM, have the same shape as in the SM, and the relative NLO corrections, reported in the lower panel, are the same as in the SM up to constant offsets. As observed in the other cases, the differences between the SM result and the ones in the various SESM scenarios can be quantified by the corresponding ∆ NLO SM obtained for the integrated decay width.

Semi-leptonic final states
In the following, we present differential distributions obtained for the semi-leptonic final states dde − e + and ν e e + dū, computed in the FJ renormalization scheme. In the MS renormalization scheme, the distributions show similar features and are not reported here. The renormalization scale used to generate the distributions is again µ r = M h . In Fig. 23 we depict differential distributions for the neutral-current final state dde − e + , including both EW and QCD corrections. The left panel shows the differential width with respect to the invariant mass of the quark pair around the Z-boson peak. Below the peak, positive NLO corrections are driven by photon radiation from the final-state quarks. Compared to the leptonic case, the radiative tail is less pronounced due to the smaller charge factor of the quarks. Since all the gluons are recombined with the quark pair, gluon radiation does not contribute to the tail [51]. The presence of the singlet does not induce any shape distortion with respect to the SM distribution, and the difference between each SESM scenario and the SM equals the difference obtained for the corresponding integrated result.   The right panel of Fig. 23 shows the angular distribution in the cosine of the angle φ between the two Z-boson decay planes. For events without gluon radiation, the final-state quarks are identified with two jets. When gluon radiation occurs, the two QCD partons with the smallest invariant mass are recombined into a single jet, so that we always obtain events with two outgoing jets. Since the jets cannot be distinguished, any observable must be invariant under the exchange of the two jets, and the cosine of φ can be reconstructed only up to a sign. Thus, in the figure, it is defined by [51] | cos φ| = k jet 1 + k jet 2 × k 1 k jet 1 × k jet 2 | k jet 1 + k jet 2 × k 1 | |k jet 1 × k jet 2 | . (7.16) Note that, in the binning of the distribution, cos φ is used instead of φ, so that the result looks different from the leptonic case reported in Fig. 21. The difference between the SESM and the SM is given by an offset, which depends on the considered scenario and is equal to the difference obtained for the corresponding integrated results. In Fig. 24 invariant-mass and angular distributions for the charged-current final state ν e e + dū are illustrated, including both EW and QCD corrections. In the left panel, the invariant-mass distribution displays the W-boson peak, and the radiative tail induced by photon radiation can be observed. As observed in the leptonic case, there is no shape difference between the SM and the considered SESM distributions, and the difference in the normalization equals the relative difference obtained for the integrated decay widths. The angular distribution, reported in the right panel of Fig. 24, shows the cosine of the angle between the electron and the hadronically decaying W boson, in the Higgs rest frame. The NLO corrections slightly deform the distribution, but there is no difference in the distortion induced by the SESM and the SM. The difference between the two equals the difference encountered for the integrated results.
In general, we observe that the presence of the singlet does not change the shape of the SM distributions. Consequently, the study of differential distributions is not helpful to discriminate the SESM from the SM. The difference between the models is given by a different normalization in the distributions and equals the relative difference observed for the integrated decay widths. Similar results were obtained for the THDM in Ref. [56].

Conclusions
To explore the nature of EW symmetry breaking at the LHC, precise theory predictions are needed. This is valid not only within the SM, but also for its extensions, since BSM physics might show small deviations from SM predictions, below the 10% level. Precision is also required in case of new discoveries, as different BSM theories lead to comparable effects, and a very high accuracy would be necessary to tell the right theoretical framework underlying the newlyobserved states.
In this paper, we have considered a SESM characterized by a Z 2 -invariant Lagrangian and non-vanishing vevs both for the SU(2) W doublet and the real singlet scalar field. The model provides two Higgs bosons, which couple to the SM fields with the same couplings of the SM Higgs, rescaled by the sine or cosine of a mixing angle α. We parametrize the extended Higgs sector by the mass M H of the heavy Higgs boson, the mixing angle α, and the scalar selfcoupling λ 12 connecting the scalar singlet and doublet. The renormalization of the theory has been performed adopting two schemes, in which the renormalized parameters α and λ 12 are defined by MS conditions, since these parameters are not directly experimentally accessible. In both schemes, all the quantities other than α and λ 12 are renormalized using OS conditions. In the first scheme the renormalized tadpoles are set to zero, as it is customary in OS renormalization schemes. This scheme introduces gauge dependences in the parametrization of observables by input parameters if MS parameters are involved, and thus should be used only within a fixed gauge. In the second scheme, following a prescription suggested by Fleischer and Jegerlehner, gauge dependences are avoided by setting unrenormalized (bare) tadpoles to zero, so that all relations between bare parameters of the theory are gauge independent.
Identifying the lighter Higgs boson h with the observed Higgs boson of mass 125 GeV, we have computed NLO EW and QCD corrections to the decays h → WW/ZZ → 4 fermions. Using the Mathematica package FeynRules, we have implemented the SESM into a FeynArts model file including the expressions of the renormalization constants, so that the model file can be used to perform NLO computations within the two considered renormalization schemes. Employing FeynArts and FormCalc, the model file has been used to produce Fortran code for the numerical computation of the matrix elements for the decays h → WW/ZZ → 4 fermions. The Fortran routines have been embedded in Prophecy4f to extend the capabilities of the Monte Carlo generator, which allows now for the computation of observables relevant for the Higgs decays to four fermions in the SESM at NLO.
The class of decay processes h → WW/ZZ → 4 fermions played a central role in the discovery of the Higgs boson and is important for the accurate characterization of the Higgs particle. We have analyzed the decays for some SESM benchmark scenarios proposed in the literature. For each scenario, we have computed the total decay width Γ h→4f for the decay of the light Higgs boson into four fermions and studied the dependence of the results on the renormalization scale µ r , solving the RGEs for the MS parameters α and λ 12 . We observe that the inclusion of NLO corrections drastically reduces the scale dependence and, consequently, the related theoretical uncertainty. Changing the scale up and down by factors of two, reduces the scale uncertainty of Γ h→4f typically from < ∼ 3−4% at LO to only < ∼ 0.5% at NLO. All these analyses have been performed using both renormalization schemes, properly converting the numerical input values between the two schemes in order to ensure a consistent comparison between the predictions obtained for specific scenarios. To this end, we have investigated the conversion of the mixing angle α between renormalization schemes and found sizeable effects which become large when approaching non-perturbative regimes. The inclusion of NLO corrections improves the agreement between the results computed in the two schemes, i.e. the renormalization scheme dependence is reduced at NLO. In the considered scenarios, the scheme dependence at the central scale typically reduces from ∼ 1% at LO to < ∼ 0.1% at NLO. Note that the inclusion of conversion effects is essential in a consistent comparison of NLO predictions obtained in the two schemes.
Comparing the NLO decay widths Γ h→4f in the SESM with the corresponding quantities in the SM, we find deviations from the SM that reach about −7% in the scenarios BHM200 ± with a heavy Higgs boson of mass 200 GeV. The NLO corrections are typically about 5−10%, but only 1−2% of those are due to effects beyond the SM. For higher values of the heavy mass M H , the (absolute values of the) deviations from the SM become smaller. Differential distributions have been produced for the SESM scenarios and compared to the SM case, and no distortions are observed on top of the SM shapes. The only observed difference is given by a constant offset, implying that differential distributions are not helpful to observe traces of the SESM. Both the FeynArts model file and the new version of Prophecy4f are ready for further applications and can be obtained from the authors upon request.

B Results for the scenario BHM400
In this appendix we show results for the scenario BHM400, defined by the input values given in Table 2. The calculations and the results are very similar to the ones described in Sects. 7.1 and 7.2 for the scenarios BHM200 + and BHM600, respectively, thus we do not discuss the details again.
The conversions of the mixing angle from the MS to the FJ scheme, and vice versa, are computed for M H and |λ 12 | values corresponding to the scenario BHM400 and shown in Fig. 30. The sign of λ 12 is fixed according to the sign of s α on the x-axis, as indicated in the figure. Red and blue lines correspond, respectively, to the complete and linearized solutions of Eq. (7.6). For values of s α in the dark-gray area, the perturbativity constraints (6.6) are violated. The light-gray areas denote where the vacuum stability condition (2.16) is violated, or where the sign of s α is flipped by the conversion (and becomes inconsistent with the sign of λ 12 ). The conversion effects are in general small and become larger when approaching the non-perturbative region.
The running of the parameters defined by MS renormalization conditions is shown, for the scenario BHM400, in Fig. 31. The scale dependence of the mixing angle is more accentuate in the FJ scheme, while the running of λ 12 is very similar in the two renormalization schemes.
In Fig. 32, we show the renormalization scale dependence of the h → 4f decay width. The left (right) panel is obtained using MS (FJ) input parameters and converting to the FJ (MS) scheme at the scale µ 0 = M h . Dashed lines correspond to LO results, solid lines include NLO EW+QCD contributions. Both scale and scheme dependence are reduced by the inclusion of NLO corrections. The scale dependence of the width reduces from ∼ 0.5% at LO to ∼ 0.1% at NLO in the MS scheme, and from ∼ 2% at LO to ∼ 0.2% at NLO in the FJ scheme, while the scheme dependence at the central scale is ∼ 0.2% at LO and < ∼ 0.1% at NLO. Figures 33, 34, and 35 show, respectively, the absolute values, the relative NLO corrections, and the deviations from the SM for the decay width Γ h→4f in the SESM as functions of s α . The parameters M H and λ 12 are fixed according to scenario BHM400. The plots on the left (right) are obtained using the MS (FJ) input scheme and converting s α to the FJ (MS) scheme. As usual, dashed and solid lines represent, respectively, LO and NLO EW+QCD results. Where relevant, the SM result is reported in green. The dashed vertical line marks the minimal s α value for which the perturbativity constraints (6.6) are fulfilled. The inclusion of NLO corrections improves the agreement between the results computed in the two schemes. For s α = 0.26, corresponding to the scenario BHM400, the decay width deviates 6−7% from the SM value, slightly depending on the input scheme used.   Figure 32: As in Fig. 11, but for benchmark scenario BHM400.  Figure 33: As in Fig. 12, but for benchmark scenario BHM400.