Gauge-independent Renormalization of the N2HDM

The Next-to-Minimal 2-Higgs-Doublet Model (N2HDM) is an interesting benchmark model for a Higgs sector consisting of two complex doublet and one real singlet fields. Like the Next-to-Minimal Supersymmetric extension (NMSSM) it features light Higgs bosons that could have escaped discovery due to their singlet admixture. Thereby, the model allows for various different Higgs-to-Higgs decay modes. Contrary to the NMSSM, however, the model is not subject to supersymmetric relations restraining its allowed parameter space and its phenomenology. For the correct determination of the allowed parameter space, the correct interpretation of the LHC Higgs data and the possible distinction of beyond-the-Standard Model Higgs sectors higher order corrections to the Higgs boson observables are crucial. This requires not only their computation but also the development of a suitable renormalization scheme. In this paper we have worked out the renormalization of the complete N2HDM and provide a scheme for the gauge-independent renormalization of the mixing angles. We discuss the renormalization of the $\mathbb{Z}_2$ soft breaking parameter $m_{12}^2$ and the singlet vacuum expectation value $v_S$. Both enter the Higgs self-couplings relevant for Higgs-to-Higgs decays. We apply our renormalization scheme to different sample processes such as Higgs decays into $Z$ bosons and decays into a lighter Higgs pair. Our results show that the corrections may be sizeable and have to be taken into account for reliable predictions.


Introduction
Even after the discovery of the Higgs boson by the LHC experiments ATLAS [1] and CMS [2] there remain many open questions that cannot be solved within the Standard Model (SM). This calls for New Physics (NP) extensions, which feature predominantly extended Higgs sectors. The precise investigation of the Higgs sector has become an important tool in the search for NP, in particular since its direct manifestation through the discovery of new non-SM particles remains elusive. Among the beyond-the-SM (BSM) Higgs sectors those with singlet and doublet extensions are particularly attractive as they are at the same time rather simple and compatible with custodial symmetry. The 2-Higgs-doublet model (2HDM) [3][4][5] is interesting due to its relation to supersymmetry and has been extensively studied and considered as a possible benchmark model in experimental analyses. It features 5 physical Higgs bosons, 2 CP-even and 1 CP-odd neutral states and a charged Higgs pair. The next-to-minimal 2HDM (N2HDM) is obtained upon extension of the 2HDM by a real singlet field with a Z 2 parity symmetry. It contains in its symmetric phase a viable Dark Matter (DM) candidate. The N2HDM has been the subject of numerous investigations, both in its symmetric [6][7][8][9][10][11][12][13][14][15][16][17][18][19] and in its broken phase [20][21][22]. The Higgs sector of the latter consists after electroweak symmetry breaking (EWSB) of 3 neutral CP-even scalars, 1 pseudoscalar and a charged Higgs pair. With the Higgs mass eigenstates being superpositions of the singlet and doublet fields the N2HDM entails an interesting phenomenology, namely the possibility of a light Higgs boson, which is not in conflict with the experimental Higgs data in case of a sufficiently large singlet admixture so that its couplings to SM particles are suppressed. The enlarged Higgs sector together with the possibility of light Higgs states allows for cascade Higgs-to-Higgs decays that provide alternative production channels for the heavier Higgs bosons and also give access to the trilinear Higgs self-couplings. Their measurement provides important insights in the understanding of the Higgs mechanism [23][24][25].
Obviously, any NP extension has to comply with the relevant theoretical and experimental constraints. Thus, also the N2HDM has to provide at least one Higgs boson with a mass of 125 GeV compatible with the LHC data on the discovered Higgs resonance [26]. The additional Higgs bosons must not violate the LHC exclusion limits. The compatibility with the electroweak (EW) precision data has to be guaranteed as well as the compatibility with B-physics and lowenergy constraints. The symmetric N2HDM furthermore has to provide a DM candidate that complies with the DM observables. From the theoretical point of view, the N2HDM Higgs potential has to be bounded from below, its vacuum has to be the global minimum and perturbative unitarity has to be respected. In [21], part of our group investigated the N2HDM in great detail with respect to these constraints. The allowed parameter space was determined and the phenomenological implications were investigated. In the course of this work the model was implemented in HDECAY [27,28]. The generated code, N2HDECAY [29], computes the N2HDM Higgs decay widths and branching ratios including the state-of-the-art higher order QCD corrections and off-shell decays. The model was furthermore included in ScannerS [30,31] along with the theoretical conditions and the available experimental constraints, which then allowed to perform extensive parameter scans for the model. In [22], the work was extended and we compared the N2HDM to other NP extensions with the aim to work out observables that can be used to distinguish between various well-motivated BSM Higgs sectors by using collider data.
Since the discovered Higgs bosons behaves very SM-like [32][33][34][35], the search for NP in the Higgs sector requires on the theoretical side precise predictions for parameters and observables including higher-order (HO) corrections. In the framework of the 2HDM, some of the authors of this work provided an important basis for the computation of HO corrections in the 2HDM by working out a manifestly gauge-independent renormalization of the two 2HDM mixing angles α and β, which is also numerically stable and process independent [36]. These angles, which diagonalise the neutral CP-even and the neutral CP-odd or charged Higgs sectors, respectively, enter all Higgs couplings so that they are relevant for Higgs boson phenomenology. We completed the renormalization of the 2HDM Higgs sector in [37] by investigating Higgs-to-Higgs decays at EW next-to-leading order (NLO). Subsequent works [38][39][40] on the 2HDM renormalization applied different approaches and renormalization conditions, confirming our findings where they overlapped. 1 The renormalization of the N2HDM is more involved due to the additional mixing angles and the additional vacuum expectation value related to the singlet field in the broken phase. One of our authors worked on the renormalization of the SM extended by a real singlet field, cf. [42]. In this paper, we combine our expertise gained in the renormalization of the 2HDM and the singlet-extended SM and provide the complete renormalization of the N2HDM. The renormalization of the mixing angles α i (i = 1, 2, 3) of the neutral sector and the angle β of the CP-odd/charged sector is manifestly gauge independent as well as process independent.
Where not parametrically enhanced, it is furthermore numerically stable with respect to missing higher order corrections. We will demonstrate this in the numerical analysis where we explicitly compute the NLO EW corrections to sample Higgs decays. We also use the occasion and clarify in this paper the notion of the alternative tadpole approach with regard to the renormalization framework applied to achieve a manifestly gauge-independent renormalization of the mixing angles. With this paper we provide another important step in the program of precise predictions for BSM Higgs sector parameters and observables including higher order corrections, an indispensable requisite for the correct interpretation of the experimental results.
The paper is organised as follows: In Sec. 2 we introduce our model, set our notation and provide the relevant couplings. Starting with Sec. 3, we describe the renormalization of the model. In Sec. 4 we explain the way we treat the tadpoles in our renormalization procedure, before we give in Sec. 5 the renormalization conditions. Section 6 is dedicated to the computation of the one-loop EW sample decay widths. In Sec. 7 we present our numerical analysis before we conclude in Sec. 8. The paper is accompanied by an extensive appendix presenting the details of the computation of the pinched self-energies in the N2HDM.

Model setup
The N2HDM is obtained from the CP-conserving (or real) 2HDM with a softly broken Z 2 symmetry upon extension by a real singlet field Φ S with a discrete symmetry, under which Φ S → −Φ S . The kinetic term of the two SU (2) L Higgs doublets Φ 1 and Φ 2 and the singlet field Φ S is given by in terms of the covariant derivative where τ a denote the Pauli matrices, W a µ and B µ the SU (2) L and U (1) Y gauge bosons, respectively, and g and g the corresponding gauge couplings. The scalar potential built from the two SU (2) L Higgs doublets and the scalar singlet can be written as 3) The first two lines correspond to the 2HDM part of the N2HDM, and the last line contains the contribution of the singlet field Φ S . The potential is based on two Z 2 symmetries, where the first one is the trivial generalization of the usual 2HDM Z 2 symmetry to the N2HDM, It is softly broken by the term involving m 2 12 . Its extension to the Yukawa sector ensures the absence of tree-level Flavour Changing Neutral Currents (FCNC). The second Z 2 symmetry on the other hand, under which is not explicitly broken. After EWSB the neutral components of the Higgs fields develop vacuum expectation values (VEVs), which are real in the CP-conserving case. Expanding the elementary field excitations around the doublet VEVs v 1 and v 2 and the singlet VEV v S , we may write where the field content of the model is parametrised in terms of the charged complex fields φ + i (i = 1, 2), the real neutral CP-even fields ρ 1 , ρ 2 , ρ 3 ≡ ρ S and the CP-odd fields η i . The minimisation conditions of the Higgs potential, where the brackets denote the vacuum state, require the terms linear in the Higgs fields, the tree-level Higgs tadpole parameters T i (i = 1, 2, 3), to vanish in the vacuum. Equation (2.7) leads to the three minimum conditions with At lowest order, the three tadpole conditions can be used to trade the mass terms m 2 11 , m 2 22 and m 2 S in favor of the other parameters of the potential. However, non-vanishing tadpole contributions are relevant at higher orders and must be included in the renormalization procedure, this being the reason why we shall retain them in our notation. The mass matrices of the Higgs fields in the gauge basis are obtained from the second derivatives with respect to these fields after replacing the doublet and singlet fields in the Higgs potential by the parametrisations (2.6). Due to charge and CP conservation the 7 × 7 mass matrix decomposes into three blocks. These are given by 2 × 2 matrices for the charged and the CP-odd fields, respectively, and a 3 × 3 matrix for the CP-even states. The former two are identical to the 2HDM case and read where we have kept explicitly the dependence on the tadpole parameters. They can be diagonalised as with the rotation matrix where we have introduced the abbreviations sin x ≡ s x and cos x ≡ c x . This yields the neutral CP-odd mass eigenstates, G 0 and A, and the charged mass eigenstates, G ± and H ± , respectively. The would-be Goldstone bosons G 0 and G ± are massless. Due to the additional real singlet field, the CP-even neutral sector differs from the 2HDM, now featuring a 3 × 3 mass matrix. In the basis (ρ 1 , ρ 2 , ρ 3 ) it can be cast into the form where t β stands for the ratio and v is defined as with v ≈ 246 GeV denoting the SM VEV. We have furthermore used Eqs. (2.8)-(2.10) to trade the mass parameters m 2 11 , m 2 22 and m 2 S for v, t β and v S . It is diagonalised by the rotation matrix R(α i ), which can be parametrised in terms of three mixing angles α 1 to α 3 as (2.20) Without loss of generality the angles can be chosen in the range (2.21) The mass eigenstates H 1 , H 2 and H 3 are obtained from the gauge basis (ρ 1 , ρ 2 , ρ 3 ) as 22) and the diagonal mass matrix D 2 ρ is given by We use the convention where the mass eigenstates are ordered by ascending mass as The full set of the N2HDM parameters is given by the parameters of the N2HDM potential Eq. (2.3), the VEVs and the free parameters of the SM. We hence have the following set of free parameters in the gauge basis of the N2HDM where y Ψ denotes the Yukawa couplings. For the renormalization of the model it is convenient to relate as many parameters as possible to physical parameters, like for example masses and the electric charge. This allows then to apply physical conditions in the renormalization of the respective parameters. Furthermore, the minimum conditions can be used to trade m 2 11 , m 2 22 and m 2 S for the tadpole parameters T 1,2,3 . Denoting by m Ψ the fermion masses, by m W and m Z the W and Z boson masses, respectively, and by e the electric charge, the 'physical' set of N2HDM parameters is given by We will specify in the following sections how these parameters get renormalized in our way of treating the tadpoles. Note also that later in our renormalization procedure we will express v S through a physical quantity that depends on it, given by a Higgs-to-Higgs decay width.
For the computation of the electroweak corrections to the Higgs decays we need the Higgs couplings, which we briefly summarise here. Since the singlet field ρ 3 does not couple directly to the SM particles, any change in the tree-level Higgs couplings with respect to the 2HDM is due to the mixing of the three neutral fields ρ i (i = 1, 2, 3). This means that any coupling not involving the CP-even neutral Higgs bosons remains unchanged compared to the 2HDM and can be found e.g. in [5]. Introducing the Feynman rules for the coupling of the Higgs fields H i to the massive gauge bosons V ≡ W, Z via  where g H SM V V denotes the SM Higgs coupling factor, we obtain the effective couplings The SM coupling in terms of the gauge boson masses m W and m Z , the SU (2) L gauge coupling g and the Weinberg angle θ W , is given by In Table 1 we list the effective couplings after replacing the R ij by their parametrisation in terms of the mixing angles.
In the Yukawa sector there exist four types of coupling structures after extending the Z 2 symmetry (2.4) to the Yukawa sector to avoid tree-level FCNCs. They are the same as in the 2HDM and summarised in Table 2. The CP-even H i Yukawa couplings can be derived from the N2HDM Yukawa Lagrangian The effective coupling factors κ H i f f in terms of the mixing matrix elements R ij and the mixing angle β are provided in Table 3. Replacing the R ij by their parametrisation in terms of the α i results in the effective coupling expressions given for type I and II in Table 4.
For the H i couplings to the Z boson and the pseudoscalar A or the Goldstone G 0 the Type I Feynman rules read where p A , p G 0 and p H i are the incoming four-momenta of the pseudoscalar, the Goldstone boson and the H i , respectively. The tilde over the coupling factor for the pseudoscalar indicates that it is not an effective coupling in the sense introduced above, as it is not normalized to a corresponding SM coupling, since there is no SM counterpart. The Feynman rules for the H i couplings to the charged pairs W ∓ and H ± or G ± read where p H ± and p G ± denote the four-momenta of H ± and G ± and again all momenta are taken as incoming. The coupling factorsκ H i V H are listed in Table 5.
The trilinear Higgs self-couplings needed for the Higgs decays into a pair of lighter Higgs bosons are quite lengthy. For their explicit form, we refer the reader to the appendix of Ref. [21]. Note finally, that by letting α 1 → α + π/2 and α 2,3 → 0, we obtain the limit of a 2HDM with an additional decoupled singlet. By the shift π/2 the usual 2HDM convention is matched, and α diagonalises the 2 × 2 mass matrix in the CP-even Higgs sector yielding the two CP-even mass eigenstates h and H, respectively, with m h ≤ m H by convention. Hence, (2.35)

Renormalization
The computation of the EW corrections to the Higgs decays involves ultraviolet (UV) divergences. Decays with external charged particles additionally induce infrared (IR) divergences. The UV divergences are canceled by the renormalization of the parameters and wave functions involved in the process. In the following we will present the renormalization of the N2HDM Higgs sector. For the purpose of this work we must deal with the renormalization of the electroweak and the Higgs sectors. With the main focus being on the renormalization of the N2HDM Higgs sector, in the sample decays presented in the numerical analysis we do not include processes that require the treatment of IR divergences or the renormalization of the fermion sector. Note also that we do not need to renormalize the gauge-fixing Lagrangian since we choose to write it already in terms of renormalized fields and parameters [43][44][45]. In the renormalization of the N2HDM Higgs sector we closely follow the procedure applied in the 2HDM renormalization of Refs. [36,37]. There, for the first time, a gauge-independent renormalization has been worked out for the 2HDM mixing angles by applying the treatment of the tadpoles of Ref. [46], which we call the alternative tadpole scheme, in combination with the pinch technique. The pinch technique allows to unambiguously extract the gauge-parameter independent parts of the decay amplitude and in particular of the angular counterterms. The N2HDM encounters four mixing angles instead of only two in the 2HDM. This leads to more complicated renormalization conditions compared to the 2HDM, as will be shown below. Additionally, the pinched self-energies needed in this renormalization program have to be worked out explicitly for the N2HDM. This has been done here for the first time. Since the formulae are quite lengthy, we defer them to App. A.2, which is part of App. A that is dedicated to the detailed presentation of the pinch technique in the N2HDM. We hope our results to be useful for further works on this subject in the future.
For the renormalization we replace the bare parameters p 0 , that are involved in the process and participate in the EW interactions, by the renormalized ones, p, and the corresponding counterterms δp, Denoting generically scalar and vector fields by Ψ, the fields are renormalized through their field renormalization constants Z Ψ as Note that in case the different field components mix Z Ψ is a matrix.
Gauge sector: The counterterms to be introduced in the gauge sector are independent of the Higgs sector under investigation. For convenience of the reader and to set our notation, we still repeat the necessary replacements here. The massive gauge boson masses and the electric charge are replaced by 2 (3.5) The gauge boson fields are renormalized by their field renormalization constants δZ, Fermion sector: Although not needed in the computation of our sample decay widths in the numerical analysis, for completeness we also include the renormalization of the fermion sector. The counterterms of the fermion masses m f are defined through And the bare left-and right-handed fermion fields are replaced by their corresponding renormalized fields according to Higgs sector: The renormalization is performed in the mass basis and the mass counterterms are defined through The field Φ stands generically for the N2HDM Higgs mass eigenstates, Φ ≡ H 1 , H 2 , H 3 , A, H ± . The replacement of the fields by the renormalized ones and their counterterms differs from the 2HDM case only by the fact that the wave function counterterm matrix in the CP-even neutral Higgs sector is now a 3 × 3 instead of a 2 × 2 matrix. Hence, (3.14) And for the mixing angles we make the replacements For the soft Z 2 -breaking mass parameter m 2 12 , finally, we replace m 2 12 → m 2 12 + δm 2 12 . (3.17) The tadpoles vanish at leading order, but the terms linear in the Higgs fields get loop contributions at higher orders. It must therefore be ensured that the correct vacuum is reproduced also at higher orders. As outlined in the following, there are two different approaches, depending on whether one chooses the tadpoles or the VEVs to be renormalized. The tadpole parameters T i (i = 1, 2, 3) and the VEVs v 1,2,S are correspondingly replaced by

Treatment of the tadpoles
The renormalization conditions fix the finite parts of the counterterms. Throughout this paper we will fix the renormalization constants for the masses and fields through on-shell (OS) conditions. Using an OS scheme provides an unambiguous interpretation of the bare parameters in the classical Lagrangian in terms of physically measurable quantities. In Ref. [36] it has been shown that the renormalization of the 2HDM mixing angles requires special care. Schemes used in the literature before, which are based on the definition of the counterterms through off-diagonal wave function renormalization constants and a naive treatment of the tadpoles, were shown to lead to gauge-dependent quantities. In order to cure this problem, in [36] for the first time a renormalization scheme has been worked out in which the angular counterterms are explicitly gauge independent. This guarantees the gauge independence of the decay amplitudes also in case the angular counterterms are not defined via a physical scheme as given e.g. by the renormalization through a physical process. The renormalization scheme developed in [36] is based on the combination of the alternative tadpole scheme with the pinch technique. The pinch technique allows for the extraction of the truly gauge-independent parts of the angular counterterms and requires the use of the alternative tadpole scheme.
As alluded to above, we treat the tadpoles in the alternative tadpole scheme in order to be able to define the angular (and also mass) counterterm in a gauge-independent way. While this procedure has been introduced in [36], we take here the occasion to explicitly pin down the differences between the standard and the alternative tadpole scheme. This, in particular, also reveals how these differences reflect in the renormalization of the singlet VEV.
The basic difference between the two schemes is the fact that in the alternative scheme as introduced by Fleischer and Jegerlehner in [46], also referred to by 'FJ' in the following, the VEVs are renormalized, while in the standard scheme the tadpole parameters are renormalized. We call the proper VEV the all-order Higgs vacuum expectation value Φ = v/ √ 2. It represents the true ground state of the theory and is connected to the particle masses and electroweak couplings. At tree level the proper VEV and the bare VEV coincide while at arbitrary loop orders the proper VEV corresponds to the renormalized VEV. In the alternative tadpole scheme the proper VEV coincides with the tree-level VEV and hence is gauge-parameter independent. In this scheme one renormalizes the VEV explicitly and its counterterm δv is fixed by ensuring the proper VEV to be v/ √ 2 = v tree / √ 2 to all orders. This renormalization condition yields δv = T loop /m 2 H , where T loop denotes the tadpole parameter at loop level. The condition generalises to multi-Higgs sectors, and we will show below in the example of the N2HDM, how the renormalization condition for the VEV counterterm is obtained. In practice, this scheme is equivalent to inserting tadpole graphs explicitly in the calculations. Since at loop level the proper VEV is given by the renormalized one, and in the FJ scheme coincides with the tree-level VEV, we have When a given v-dependent Lagrangian is used at higher orders these tree-level parameters {g, m W } tree still have to be renormalized, and they are then replaced by their corresponding renormalized parameters as It is important to note that ∆v is a mere label and not a VEV counterterm as such. This makes obvious that δv and ∆v are completely unrelated. In particular, they feature a totally different divergence structure. Figure 1 depicts the renormalization condition for the alternative tadpole scheme. In the standard scheme, on the other hand, the proper VEV is obtained from the minimisation of the gauge-dependent loop-corrected potential and hence is in principle gauge dependent. An equivalent condition to the FJ scheme requires the renormalized tadpole to vanish. Together with the requirement of the tree-level tadpole to be zero, this fixes the tadpole counterterm, which features here explicitly, as the tadpole is an input parameter in the standard scheme, cf.  For the singlet VEV v S a similar distinction, i.e. δv S versus ∆v S has to be made. When v S is related to measurable parameters the NLO VEV shift ∆v s denotes the corresponding combination of parameter counterterms, similarly to Eq. (4.21). In Ref. [47] it was shown that, in an R ξ gauge, a divergent part for ∆v S in the standard scheme is precluded at one loop if the scalar field obeys a rigid invariance. This is the case for typical singlet-extended Higgs sectors, e.g. the real singlet model [42], and thereby the N2HDM singlet scalar. In all these cases the singlet field is disconnected from the gauge sector and hence invariant under global gauge transformations. The conclusion of Ref. [47] relies on the use of the standard scheme, where the renormalized VEV coincides with the loop-corrected one as the renormalized tadpoles are set to zero 3 . However, this no longer applies if the VEVs are renormalized in the alternative tadpole scheme. In this case ∆v FJ S becomes indeed a UV-divergent quantity. We can prove it to cancel part of the UV poles that genuinely appear if one-loop amplitudes are computed in the FJ-scheme, when the corresponding tree-level amplitudes are directly sensitive to the singlet VEV v S . Salient examples are the Higgs-to-Higgs decays, which we discuss in detail in Section 6.

Alternative tadpole scheme for the N2HDM
In the following, we elaborate in detail the implications of the alternative tadpole scheme. We derive the necessary relations for the N2HDM, highlighting the differences with respect to the 2HDM case, derived in [36]. At tree level the minimum conditions of the N2HDM potential lead to the three relations Eqs.  These can be used to replace the parameters m 2 11 , m 2 22 and m 2 S by the VEVs v 1 , v 2 and v S . Note, however, that at arbitrary loop order, this may only be done after the proper VEVs are 3 Let us also notice that Ref. [47] distinguishes two (equivalent) parametrisations for the renormalization transformation of a generic scalar field VEV, where √ ZΦ is the field renormalization constant of the respective scalar field, whereas δv quantifies how the VEV itself is shifted differently by higher-order contributions with respect to the field. In our current conventions, δv → ∆v and δv → ∆v. The results of Ref. [47], together with [42], show that for a gauge-singlet scalar the quantity ∆vs in the standard scheme is UV finite at one loop order. taken into account in the Higgs potential. More precisely, at NLO the VEVs are modified in order to take into account the NLO effects, as v bare In the alternative tadpole scheme, δv 1 , δv 2 and δv S correspond to the proper doublet and singlet VEV counterterms in the gauge basis. In turn, v ren i are the proper VEVs, i.e. in the FJ scheme the renormalized VEVs (coinciding with the tree-level VEVs), and hence the VEVs that generate the necessary mass relations for the gauge bosons, fermions and the scalars. The VEVs are called the proper VEVs if the gauge-invariant relations presented in Fig. 1 (for the SM case) are fulfilled at all orders, which means that the VEVs represent the true vacuum state of the theory at all orders in perturbation theory. At NLO, we insert the relations Eq. (4.24) into the tadpole relations Eqs. (2.8)-(2.10). At NLO, the left-hand side of the equations is given by (4.25) We then get the NLO expressions for Eqs. (2.8)-(2.10), . By comparing with the squared mass matrix M 2 ρ of Eq. (2.17) we find analytically Rotation to the mass basis yields The latter identity is helpful in practice, as the calculation of the tadpole diagrams is usually performed in the mass basis, but the VEV shifts are introduced most conveniently in the gauge basis. Rewriting Eq. (4.30), the quantities δv H i can be interpreted as connected tadpole diagrams, containing the Higgs tadpole and its propagator at zero momentum transfer, (4.32) We want to emphasize again that in the alternative tadpole scheme Eq. (4.32) defines the counterterms of the vacuum expectation values. In contrast to the standard scheme, no tadpole counterterms are introduced. Tadpole graphs appear through the gauge-invariant condition in Fig. 1.
Once the leading-order VEVs are promoted to higher orders, namely by inserting Eq. (4.24) into a generic VEV-dependent Lagrangian L(v 1 , v 2 , v S ), the contribution of the VEV counterterms δv 1 , δv 2 and δv S , as given by Eq. (4.31), is equivalent to introducing explicit tadpole graphs in all loop amplitudes. Moreover, all tree-level relations between the VEVs and the weak sector parameters (masses, coupling constants) hold again. In particular, for the doublet VEVs this means with (v 2 By applying the renormalization conditions for the VEVs, the tree-level VEVs ensure the true ground state of the potential. Since they are not directly related to a physical observable, we express the FJ-renormalized doublet VEVs in terms of physical tree-level parameters, here m W , g and the mixing angle β. In higher order calculations, these parameters are then renormalized by choosing physical renormalization conditions. 4 To better illustrate the implications of the alternative tadpole scheme, we consider the scalar-vector-vector vertex between the physical H 1 and a W boson pair. We first define the Feynman rules, needed in the following, by The coupling constants for the triple vertex in terms of the mixing angles and the VEVs v 1 and v 2 are and for the quartic vertices When expressing the couplings in terms of the VEVs, care has to be taken to differentiate between the angle β in the sense of a mixing angle and β in the sense of the ratio of the VEVs.
Only the latter is to be replaced by the VEVs that are to be renormalized. The same distinction must be applied for the α i . Note that in all couplings but the trilinear and quartic Higgs self-couplings the angles α i have the roles of mixing angles. Only in the Higgs self-couplings, the α i partly appear in the sense of the ratio of N2HDM potential parameters. Bearing these considerations in mind, we see that the quartic couplings do not receive any δv i , whereas g H 1 W W contains β as ratio of the VEVs. Instead, the angles α 1 and α 2 are mixing angles here. At NLO, we therefore have to make the replacement The subscript 'trunc' means that all Lorentz structure of the vector bosons as well as the Lorentz structure of the coupling has been suppressed here for simplicity. The second term in the second line generates, through the VEV counterterms δv i , the tadpole diagrams contributing to the scalar-vector-vector vertex. On the other hand, as the VEVs in this expression have already been to fix the (FJ-renormalized) VEVs v tree i in terms of the tree-level weak sector parameters and the angle β. 5 At loop level the EW parameters and mixing angles that enter the coupling, here g, m W , β, α 1 and α 2 have to be renormalized, i.e. we replace them by their renormalized values plus the corresponding counterterms, cf. Eq. (3.1). We then get for the vertex of Eq. (4.39) The exact form of these counterterms 6 depends on the renormalization conditions, which will be given in the next section.
Our derivation also shows the difference with respect to the 2HDM, namely the last two terms in Eq. (4.39) do not arise in the 2HDM. They are due to the additional singlet-doublet mixing and have no counterpart in a pure 2HDM structure (cf. Eq. (A.61) of [36]).
As a final remark, let us summarise the key differences with respect to the standard tadpole scheme. In the latter case, VEV counterterms of the form of Eq. (4.32) are strictly speaking not introduced. Instead, one introduces renormalized tadpoles and tadpole counterterms, fulfilling the same condition as in Fig. 1 -that is, T ren In doing so, the VEVs correspond to the ground state of the loop-corrected scalar potential, and the corresponding VEV relations to weak sector parameters hold order-by-order. Due to the fact that in the standard tadpole scheme one considers the VEVs from the one-loop corrected potential (in contrast to the alternative scheme, where one considers the tree-level VEVs), VEV diagrams in the self-energies and vertices explicitly vanish and thus need not be taken into account, at the expense of defining mass counterterms which become manifestly gauge dependent.
In practice, the rigorous introduction of the VEV counterterms in the alternative tadpole scheme yields the following rules for its application in the renormalization of a generic process within the N2HDM: 1. Include explicit tadpole contributions in all self-energies used to define the (off-diagonal) wave function renormalization constants 7 and wherever the self-energies appear in the counterterms, such that now Σ tad (p 2 ) contains the additional tadpole contributions, cf. 2. Include explicit tadpole contributions in the virtual vertex corrections, if the tadpole insertions are connected to an existing coupling. This is applicable e.g. to all triple Higgs self-interactions as well as to the Higgs couplings to gauge bosons.
In the alternative tadpole scheme not only the angular counterterms but also the mass counterterms become gauge independent. This has been shown for the electroweak sector in [48]. All counterterms of the electroweak sector have exactly the same structure as in the standard scheme. Only the self-energies Σ have to be replaced by the self-energies Σ tad containing the tadpole contributions. Note however, that there are no tadpole contributions to the transverse photon-Z self-energy Σ T γZ nor to the transverse photon self-energy Σ T γγ so that Having introduced the tadpole scheme, we now list explicitly the counterterms needed in the computation of the electroweak corrections. In particular, we illustrate the renormalization of the N2HDM Higgs sector.

Renormalization conditions
With the previous section we are now able to specify the counterterms needed in the renormalization of the N2HDM. Those of the EW and Yukawa sector correspond to the ones of the SM, while differences obviously arise in the Higgs sector itself. For completeness, all counterterms of the model will be listed, although not all of them will be necessary to study the sample processes discussed in Section 7.

Counterterms of the gauge sector
The gauge bosons are renormalized through OS conditions implying the mass counterterms where T denotes the transverse part of the self-energy including the tadpole contributions. The wave function renormalization constants that guarantee the correct OS properties are given by Note that in Eqs. (5.43) and (5.44) they are the same in the standard and in the alternative tadpole scheme introduced above. The reason is that the tadpoles are independent of the external momentum so that the derivatives of the self-energies do not change. Furthermore, Σ T γZ is identical in both schemes, as alluded to above. For better readability we therefore drop the superscript 'tad' here and wherever possible. For the same reasons the counterterm for the electric charge is invariant with respect to the choice of the tadpole scheme. The electric charge is renormalized to be the full electron-positron photon coupling for OS external particles in the Thomson limit. This implies that all corrections to this vertex vanish OS and for zero momentum transfer. The counterterm for the electric charge in terms of the transverse photon-photon and photon-Z self-energies reads [49] δZ α(0) The sign in the second term of Eq. (5.45) differs from the one in [49] because we have adopted different sign conventions in the covariant derivative of Eq. (2.2). In our computation we will use the fine structure constant at the Z boson mass α(m 2 Z ) as input. This way the results are independent of large logarithms due to light fermions f = t. The counterterm δZ e is therefore modified as [49] δZ where the transverse part of the photon self-energy Σ T γγ in Eq. (5.47) includes only the light fermion contributions. The calculation of the EW one-loop corrected Higgs decay widths also requires the renormalization of the weak coupling g, which can be related to e and the gauge boson masses as Its counterterm can therefore be expressed in terms of the electric charge and gauge boson mass counterterms through

Counterterms of the fermion sector
Defining the following structure for the fermion self-energies the fermion mass counterterms applying OS conditions are given by The fermion wave function renormalization constants are determined from

Higgs field and mass counterterms
The OS conditions for the physical Higgs bosons yield the mass counterterms (i = 1, 2, 3) Having absorbed the tadpoles into the self-energies, no tadpole counterterms appear explicitly in the mass counterterms any more, in contrast to the corresponding expressions in the standard tadpole scheme. The OS conditions for the Higgs bosons yield the following wave function renormalization counterterm 3 × 3 matrix for the CP-even neutral N2HDM scalars, And in the CP-odd and charged sector we have the 2 × 2 matrices

Angular counterterms
As in the 2HDM, we renormalize the mixing angles based on the definition of the counterterms through off-diagonal wave function renormalization constants and combine this with the alternative tadpole approach together with the application of the pinch technique in order to arrive at an unambiguous gauge-independent definition of the mixing angle counterterms. Let us note that a process-dependent renormalization of the mixing angles would also lead to a gauge-independent renormalization, as shown in [36] for the 2HDM case. In the N2HDM the situation becomes more involved as four different processes need to be identified to fix all mixing angle counterterms δα i and δβ. Moreover, the construction of such a process-dependent scheme is complicated by the fact that the different Higgs decay modes typically rely on more than one mixing angle, implying that the different angular counterterms appear as linear combinations in each individual vertex counterterm. It is therefore imperative to choose a set of processes where the angular counterterm dependences enter as a linearly independent combination, such that they can be fixed unambiguously through linear combinations of the different decay widths. Moreover, all these processes have to be phenomenologically accessible. The process-dependent renormalization of the N2HDM mixing angles is hence rather unpractical from a physical point of view, and we will therefore not consider it any further.
While the expression for the counterterm in the charged and CP-odd sector, δβ, in terms of the off-diagonal wave function renormalization constants does not change with respect to the 2HDM, this is not the case for the mixing angle counterterms δα i in the CP-even sector.
We therefore present their derivation here. It is based on the idea of making the counterterms δα i (and also δβ) appear in the inverse propagator matrix and thereby in the wave function renormalization constants in a way that is consistent with the internal relations of the N2HDM. 8 This can be achieved by performing the renormalization in the physical basis (H 1 , H 2 , H 3 ), but temporarily switching to the gauge basis (ρ 1 , ρ 2 , ρ 3 ), and back again. For the CP-even sector of the N2HDM this means, (5.59) The field renormalization matrix in the mass basis can be parametrised as By identifying the off-diagonal elements with the off-diagonal wave function renormalization constants δZ H i H j (i = j), the three neutral CP-even angular counterterms are obtained as while the auxiliary counterterms δC ij do not play a role in the remainder of the discussion.
The definition of the counterterm δβ can be taken over from the 2HDM. It is derived analogously to the δα i , but from the charged and CP-odd Higgs sectors. In this case, there are altogether four off-diagonal wave function constants, while only three free parameters to be fixed. For details, we refer to Ref. [36]. There we proposed two different possible counterterm choices for β, one based on the charged and the other on the CP-odd sector. Also here we will apply these two possible choices, given by and . While the use of the alternative tadpole scheme ensures that the angular counterterms can be expressed in a gauge-independent way, at this stage they still contain a dependence on the gauge-fixing parameter. We therefore combine the virtues of the alternative tadpole scheme with the pinch technique [51][52][53][54][55][56][57][58]. The pinch technique allows us to extract the truly gauge-independent parts of the angular counterterms.

Gauge-independent pinch technique-based angular counterterm schemes
By the application of the pinch technique it possible to define pinched self-energies Σ which are truly gauge independent. They are built up by the tadpole self-energies evaluated in the Feynman gauge and extra pinched components Σ add , i.e.
where ξ V stands for the gauge fixing parameters ξ Z , ξ W and ξ γ of the R ξ gauge. By Σ add we dub the additional (explicitly ξ V -independent) self-energy contributions obtained via the pinch technique. It is important to notice that, in order to apply the pinch technique, it is necessary to explicitly include all tadpole topologies, i.e. to use the alternative tadpole scheme. In App. A we present the basic idea of the pinch technique (see also Refs. [51][52][53][54][55][56][57][58] for a detailed exposition). We exemplarily show, for the CP-even sector, how to proceed in the derivation of the pinched self-energy. Additionally, we give useful formulae on the gauge dependences of the scalar self-energies and for the application of the pinch technique in the N2HDM.
On-shell tadpole-pinched scheme The self-energy Σ add in Eq. (5.64) is explicitly independent of the gauge fixing parameter ξ V . By replacing the wave function renormalization constants in the counterterms Eqs.
And for the two chosen renormalization prescriptions of δβ we get (5.67) With this procedure we have now obtained angular counterterms that are explicitly gauge independent.
The additional contribution Σ add Hh has been given for the MSSM in [59], and the ones for the 2HDM in [36,40,60]. We have derived the contributions necessary in the N2HDM, given here for the first time (i, j = 1, 2, 3), where B 0 is the scalar two-point function [61,62], while the shorthand notation O

(x)
H i H j (x = 1, ..., 4) stands for different coupling combinations in the Higgs-gauge sector, We note that in the N2HDM the following sum rules hold, H i H j becomes the Kronecker delta δ H i H j and hence, for i = j, the additional pinched contributions in Eq. (5.68) become UV-finite by themselves as well.
In the general N2HDM case instead, Σ add H 3 H 2 contain UV-divergent poles, which nevertheless cancel as they enter the mixing angle counterterms Eq. (5.65) via the additive structure Σ add p tadpole-pinched scheme Along the same lines followed for the 2HDM in Ref. [36], we now generalise the p tadpole-pinched scheme to the N2HDM Higgs sector. Again, we replace the scalar self-energies within the mixing angle counterterms with the corresponding pinched self-energies, Σ, Eq. (5.64), which we evaluate this time at the average of the particle momenta squared [63], H j ), (G ± , H ± ) and (G 0 , A), respectively. In this way the additional selfenergies Σ add vanish, and the pinched self-energies are given by the tadpole self-energies Σ tad computed in the Feynman gauge, i.e.
The angular counterterms δα i in Eq. (5.61) then read For the counterterm δβ we get or alternatively (5.78)

Renormalization of m 2 12
The soft Z 2 breaking parameter m 2 12 enters the Higgs self-couplings. For the computation of higher-order corrections to Higgs-to-Higgs decays it therefore has to be renormalized as well. We may consider two different renormalization schemes.
Modified Minimal Substraction Scheme: One possibility is to use a modified MS scheme, cf. [37], where the counterterm δm 2 12 is chosen such that it cancels all residual terms of the amplitude that are proportional to where γ E denotes the Euler-Mascheroni constant. These terms obviously contain the remaining UV divergences given as poles in together with additional finite constants that appear universally in all loop integrals. The renormalization of δm 2 12 in this scheme is thereby given by The right-hand side of the equation symbolically denotes all terms proportional to ∆ that are necessary to cancel the ∆ dependence of the remainder of the amplitude.
Process-dependent renormalization: Alternatively, one could resort to a process-dependent scheme, in which case the divergent parts of δm 2 12 , along with additional finite remainders, are related to a physical on-shell Higgs-to-Higgs decay. While this method provides a physical definition for the counterterm, it relies on having at least one kinematically accessible on-shell Higgs-to-Higgs decay. For a generic Higgs-to-Higgs decay process H i → H j H k , where the final state pair H j H k can also be a pair of pseudoscalars, if kinematically allowed, the counterterm δm 2 12 is then fixed by imposing as renormalization condition Note that δm 2 12 is gauge independent in either of the proposed schemes, and also independently on how the tadpole topologies are treated. The key reason is that m 2 12 is indeed a genuine parameter of the original N2HDM Higgs potential before EWSB, and hence unlinked to the VEV, this being the source for the potential gauge-parameter dependences that arise at higher orders in certain schemes. In this paper we will apply the MS renormalization scheme.

One-Loop EW Corrected Decay Widths
Having elaborated in detail the renormalization scheme for the N2HDM, we compute the NLO EW corrections to a selected set of decay widths, in order to illustrate their impact. The chosen decays widths are All processes require the renormalization of the mixing angles. The Higgs-to-Higgs decays demand in addition the renormalization of m 2 12 . And the Higgs decays into CP-even pairs, Eq. (6.84), additionally involve the renormalization of v S . The chosen processes are structurally different and involve the various mixing angles in different more or less complicated combinations, allowing us to study the impact of our renormalization scheme in different situations, and enabling us to study the renormalization of the Higgs potential parameter m 2 12 as well as of the singlet VEV v S . Note finally that all these decays only involve electrically neutral particles, so that we do not encounter any IR divergences in the EW corrections.

The NLO EW corrected decay H i → ZZ
The LO decay width for the decay of a CP-even Higgs boson H i into a pair of Z bosons, is given by and depends on the mixing angles through the coupling factors The generic diagrams describing the virtual corrections contributing to the NLO decay width together with the counterterm diagram introduced to cancel the UV divergences are displayed in Fig. 4. With the decay width involving only neutral particles there are neither IR divergences nor real corrections. The corrections to the external legs in Fig. 4 (c), (f) and (g) vanish due to the OS renormalization of H i and Z, respectively, and the mixing contributions (d) and (e) are zero because of the Ward identity satisfied by the OS Z boson. The one-particle irreducible (1PI) diagrams contributing to the vertex corrections originate from the triangle diagrams with scalars, fermions, massive gauge bosons and ghost particles in the loops, depicted in the first three rows of Fig. 5, and from the diagrams involving four-particle vertices, as given by the last four diagrams of Fig. 5.
To work out the vertex counterterms, the relations The H i ZZ vertex counterterm in terms of the different parameter counterterms and wave function renormalization constants is obtained from the corresponding counterterm Lagrangian with the various counterterms given in Section 3 and the δR ij defined in Eq. (6.89). Since we apply the alternative tadpole scheme, tadpole contributions to the H i ZZ vertex have to be taken into account explicitly in the computation of the decay width. They are shown in Fig. 6. The formulae for the vertex corrections and counterterms in terms of the scalar one-, two-and three-point functions are quite lengthy so that we do not display them explicitly here. Figure 5: Generic diagrams contributing to the vertex corrections in Hi → ZZ with fermions F , scalar bosons S, gauge bosons V and ghost particles U in the loops. Figure 6: Tadpole contributions to the vertex diagrams to be included in the decay Hi → ZZ in the alternative tadpole scheme. (a)

The decay H i → AA at NLO EW
The LO decay width of the CP-even H i decay into a pair of CP-odd scalars, It is governed by the trilinear coupling where M 2 ≡ m 2 12 /(s β c β ). The EW one-loop corrections consist of the virtual corrections and the counterterm contributions ensuring the UV-finiteness of the decay amplitude. Again we do not have to deal with IR divergences nor real corrections. The virtual corrections, consisting of the corrections to the external legs and the pure vertex corrections, are shown in Fig. 7. The corrections to the external legs in Fig. 7 (c), (d) and (e) are zero because of the OS renormalization of the external fields, while diagrams (f) and (g) vanish due to a Slavnov-Taylor identity [64]. The 1PI diagrams of the vertex corrections are depicted in Fig. 8. They are given by the triangle diagrams with fermions, scalars and gauge bosons in the loops and by the diagrams containing  Figure 9: Tadpole contributions to the vertex diagrams to be included in the decay Hi → AA in the alternative tadpole scheme.
four-particle vertices. The counterterm contributions consist of the genuine vertex counterterm δg vertex H i AA and the counterterm insertions on the external legs δg field H i AA , with the δR ij given in Eq. (6.89). Working in the alternative tadpole scheme, we additionally have to take into account the vertices dressed with the tadpoles, displayed in Fig. 9.
The one-loop correction to the decay is obtained from the interference of the The NLO corrections factorise from the LO amplitude so that the loop-corrected partial width can be cast into the form Again we refrain from giving the explicit expressions for the various contributions to Γ NLO as they are quite lengthy.

Electroweak one-loop corrections to H j → H i H i
The LO decay width for the decay of a neutral CP-even Higgs boson into two identical CP-even scalars is given by (i, j = 1, 2, 3 ) with the trilinear Higgs coupling where ijk denotes the totally antisymmetric tensor in three dimensions with 123 = 1. At variance with the processes discussed so far, Higgs-to-Higgs decays in the CP-even sector are directly sensitive to the singlet VEV v S at tree level. As discussed in section 4.1, this explicit dependence must be handled with care when the NLO calculations are performed in the alternative tadpole scheme. Here, a non-vanishing UV-divergent singlet VEV shift ∆v S cancels a subset of the UV poles in the NLO Higgs-to-Higgs decay amplitude which genuinely arise in this scheme. To fix ∆v S we proceed along the same lines as for the doublet VEV. First, we identify the singlet VEV input value in this scheme with the (would-be) experimental input, to be extracted eventually through the measurement of an observable Higgs-to-Higgs decay width Γ H i →H j H j . When promoted to higher orders, the tree-level relation v tree in such a way that the (would-be) experimental value v exp S is properly written in terms of the renormalized (physical) width from which it would be extracted. Notice that the quantity ∆v S is simply a shorthand for the combination of counterterm contributions contained in Γ ct H i →H j H j -the same role that ∆v plays in Eq. (4.21) for the doublet VEV case. For our sample processes H 3 → H 2 H 2 and H 2 → H 1 H 1 discussed in the numerical analysis we assume the v S input values to be extracted from the decay H 3 → H 1 H 1 . 9 The choice of this process is of course not unique. Therefore, given that the finite parts included in ∆v S are to some degree arbitrary, we could formally resort to MS-like conditions to fix ∆v S by retaining only the UV-divergent parts contained in Γ ct H i →H j H j . In this case the v S input values could not be extracted directly from the experimental data. The relation to the to be measured v exp S would be given by a scheme-dependent finite shift. In the process-dependent framework ∆v S can be fixed through the requirement Γ NLO Factorising the NLO decay width as The choice of the process relies on the experimental feasibility of measuring it and on its dependence on vS itself. For some scenarios the parameter configurations can be such that the decay is not measurable or the dependence on ∆vS is almost vanishing, cf. also the discussion in [65] on the renormalization of the NMSSM where similar issues arise. and isolating the v S -dependent part of the corresponding self-interaction Lagrangian, The diagrams contributing to the virtual corrections of our process H j → H i H i are shown in Fig. 10. The 1PI diagrams contributing to the vertex corrections are depicted in Fig. 11 and the tadpole diagrams are shown in Fig. 12. They have to be included in the alternative tadpole scheme. The counterterm is given by the genuine vertex counterterm and the counterterm insertions on the external legs, 108) Figure 11: Generic diagrams contributing to the vertex corrections in Hj → HiHi. and and implemented by hand. The amplitudes have been analytically processed via FormCalc [70]. The dimensionally regularised loop form factors have been evaluated in the 't Hooft-Veltman scheme [71,72] and written in terms of standard loop integrals. These have been further reduced through Passarino-Veltman decomposition and evaluated with the help of LoopTools [70].
In the following we give the input parameters for the numerical evaluation. As explained in section 5 we use the fine structure constant α at the Z boson mass scale, given by [73] α(m 2 Z ) = as recommended by [74]. We consider the CKM matrix to be unity. This approximation has negligible impact on our results. The SM-like Higgs mass value, denoted by m h , has been set to [26] m h = 125.09 GeV .
Note that, depending on the parameter set, in the N2HDM any of the three neutral CP-even Higgs bosons can be the SM-like Higgs boson.
In the subsequently presented analysis we only used N2HDM parameter sets compatible experimental and theoretical constraints. These data sets have been generated with the tool ScannerS [30,31]. 10 The applied theoretical constraints require that the vacuum state found by ScannerS is the global minimum, that the N2HDM potential is bounded from below and that tree-level unitarity holds. On the experimental side, compatibility with the EW precision constraints is guaranteed by requiring the oblique parameters S, T and U to be compatible with the SM fit [77] at 2σ, including the full correlations. The constraints from B physics observables [78][79][80][81][82] and the measurement of R b [79,83] have been taken into account, as well as the most recent bound of m H ± > ∼ 580 GeV for the type II and flipped (N)2HDM [82]. For the compatibility with the LHC Higgs data we require one of the scalar states, denoted by h 125 , to have a mass of 125.09 GeV and to match the observed LHC signal rates. Furthermore, the remaining Higgs bosons have to be consistent with the exclusion bounds from the collider searches at Tevatron, LEP and LHC. For further details on these checks and the scan, we refer to [21,22].
Note that in all scenarios presented in the following we stick to the N2HDM type I, with the type II scenarios leading to the same overall results. The only difference between the models comes from the fermion loops. The Yukawa couplings are, in all Yukawa types, well-behaved functions of the α i and β because extreme values of β are already disallowed by all the constraints imposed on the model. Therefore, this is sufficient for our analysis to illustrate the effects of the EW corrections, without aiming at a full phenomenological analysis of N2HDM Higgs decays.

Results for H 2/3 → ZZ
In this section we investigate the relative size of the NLO EW corrections as well as the impact of the different renormalization schemes for the mixing angles on the decay H i → ZZ. We base our numerical analysis upon a set of representative N2HDM scenarios of phenomenological interest. To this aim we select among the generated parameter points compatible with the theoretical and experimental constraints scenarios that either have a large or a small LO branching fraction (BR) into ZZ. Discarding the SM-like decay of the H 1 fixed to be the 125 GeV Higgs boson, we select hence four scenarios, two for H 2 and H 3 , respectively, which we denote by 'BRH2/3high' and 'BRH2/3low' for high and low branching ratio scenarios. The corresponding input parameters are listed in Table 6. Note that, if not stated otherwise, the mixing angles are understood to be the angles defined in the OS tadpole pinched scheme (pOS) with δβ defined via the charged sector, denoted by the superscript 'c'. 11 The suppressed branching fractions in the BRlow scenarios are due to a small tree-level coupling to ZZ of the decaying Higgs boson. The branching fractions given in this table have been obtained with the Fortran code N2HDECAY. 12 We insured to consider purely OS decays into massive gauge bosons in N2HDECAY, as we do not include any gauge boson off-shell effects in the NLO computation.
In Table 7 we present for all four benchmark scenarios the results for the LO and the NLO width as well as the relative corrections ∆Γ. They are given for four different renormalization schemes. These consist of the p and the pOS tadpole pinched schemes, which employ two different renormalization scales, and for these additionally the two possibilities to renormalize β, either via the charged sector (denoted by 'c') or the CP-odd sector (denoted by 'o'). The relative corrections are defined as When computing the NLO EW corrected decay width Γ NLO in a different renormalization scheme b than the one of the input parameters p, scheme a, these parameters first have to be 11 While the scheme choice is not relevant for the LO width alone, it becomes important when the NLO EW corrections are included. The renormalization of the parameters then fixes the scheme of the input parameters at LO. 12   BRH3ZZlow 4.2 5.0 4.8 9.3 Table 7: Higgs decay widths (in GeV) at LO and NLO EW accuracy as well as the relative corrections for the N2HDM benchmarks presented in Table 6 and four different renormalization schemes.
converted to the scheme that is applied. We perform this conversion for the mixing angles α and β through (p = α, β) where δp denotes the counterterm in either scheme a or scheme b. With the thus obtained input parameters in scheme b we compute the quantity ∆Γ NLO and the LO width Γ LO , to which we normalize the relative correction. 13 The relative corrections for the scenarios with relatively large branching ratios turn out to be of moderate size with values between 13.2 and 15.6%. The variation due to different renormalization schemes is at most 1.9%, indicating a relatively small theoretical error due to missing higher order corrections. For the low branching ratio scenarios on the other hand, the differences between the results for the various renormalization schemes are substantial. This points towards a large theoretical error due to missing higher order corrections. A reliable prediction in these cases would require the inclusion of corrections beyond one-loop order. This is to be expected as the tree-level widths are very small in these scenarios so that the one-loop correction effectively becomes the leading contribution to the width. When changing from the charged to the CP-odd based renormalization of β, the change in the relative corrections is rather mild for most of the scenarios. This is because the two different scales, m H ± or m A , involved in these two renormalization schemes of β are close in our scenarios.

Results for H 2/3 → AA
Here we study the decay into a pair of pseudoscalars and again concentrate on the decays of the heavier Higgs bosons H 2 and H 3 and choose scenarios where H 1 is the 125 GeV Higgs boson 14 and with low and high branching ratios for H 2/3 → AA, respectively. The corresponding benchmark scenarios are called 'BRH2/3AAhigh' and 'BRH2/3AAlow', with the input values summarised in Tab. 8 together with the LO total widths and branching ratios computed with N2HDECAY. The input mixing angles are given in the pOS scheme and the β renormalization is based on use as input parameters mW , mZ and α, while in N2HDECAY all decay widths are expressed in terms of the Fermi constant GF as input value. Including in our LO results the SM correction ∆r SM [84][85][86], which relates mW to GF , would bring the derived Fermi constant numerically very close to the PDG value GF = 1.166 · 10 −5 GeV −1 used in N2HDECAY. 14 We do not consider H1 decays into AA. They would require mA to be below about 65 GeV and care would have to be taken to keep the decay H1 → AA small enough to still be compatible with the LHC Higgs data.
6.  Table 9: Higgs decay widths (in GeV) at LO and NLO EW accuracy as well as the relative corrections for the N2HDM benchmarks presented in Table 8 and four different renormalization schemes. The renormalization scale of m 2 12 is set to µR = 2mA.
In Table 9 we display for all four benchmark scenarios the LO and NLO widths as well as the relative corrections ∆Γ. They are given for the four different renormalization schemes, p c/o , pOS c/o . As can be inferred from the table, for the BRhigh scenarios we obtain moderate corrections of O(10)%, i.e. of the same order as for H 2/3 → ZZ. The associated theoretical uncertainties are very mild, as indicated by the the rather small influence of the renormalization schemes of the mixing angles, which lead to a change of at most 2.8%, when considering all four schemes. The H 2/3 → AA decays in the BRlow scenarios, on the contrary, are dominated by the loop effects. Here, the small trilinear Higgs coupling suppresses the tree-level width. At one loop, however, the Higgs decay is also sensitive to the additional trilinear Higgs couplings, some of which being very large as a result of the heavy Higgs masses and the large m 2 12 scale -yet in agreement with the unitarity constraints. This results in very large NLO effects, and also induces the strong dependence on the renormalization scheme and the renormalization scale. This reflects the fact that the H 2/3 → AA decays in these benchmarks are effectively loop-induced and higher order corrections beyond the one-loop level need to be considered to make reliable predictions. These sizable higher-order effects are particularly apparent in the BRH2AAlow scenario, where some of the renormalization schemes even lead to negative and hence unphysical NLO widths. Note, furthermore, that the change when switching from the charged to the CP-odd based renormalization schemes for β is now larger when compared to the results in Table 7, due to the now wider separation between the scales given by the charged and the CP-odd Higgs mass as compared to the scenarios studied for the decays into a Z boson pair. Finally, we consider the decay of a heavy neutral CP-even Higgs boson into a pair of lighter CPeven Higgs bosons. We evaluate the NLO EW corrections for a number of illustrative scenarios, given in Table 10. The scenarios have been chosen such that their Higgs mass spectra allow simultaneously for the OS H 3 → H 2 H 2 and H 2 → H 1 H 1 decays. Furthermore, the chosen large m 2 12 parameter insures these heavy Higgs mass scenarios to be in agreement with the unitarity and vacuum stability constraints. All scenarios feature Higgs-to-Higgs decay branching ratios that are of moderate size. Only HHHIV features a H 2 branching ratio into H 1 H 1 that is dominating. All input mixing angles are assumed to be given in the pOS scheme, with charged sector-based renormalization for the angle β, and m 2 12 is assumed to be defined at the renormalization scale given by the total final state mass, µ R = 2m H i . The LO total widths and branching ratios in this table have been obtained from N2HDECAY.
In Table 11 we summarise the relative NLO corrections for the various decays. Note, that the decay process H 3 → H 1 H 1 appears only at LO because we use it for the renormalization of v S , as explained in detail in Section 6. The sizeable m 2 12 and heavy Higgs mass values imply large Higgs self-couplings and thereby enhanced contributions from the virtual Higgs exchanges. On the other hand, these enhancements are partly damped by the inverse Higgs mass powers in the Higgs-mediated loops. The balance between these dynamical features governing the Higgs-mediated loops, and how they interplay with the remaining gauge boson and fermionmediated one-loop contributions, determines the overall size of the NLO EW effects. For most of the decays, the relative NLO corrections are moderate and reach at most 21%. Accordingly, they show a mild renormalization scheme and scale dependence with changes in the predicted   Table 11: Higgs decay width predictions (in GeV) at LO and NLO EW accuracy as well as the relative corrections for the N2HDM benchmarks presented in Table 10 and four different renormalization schemes.
NLO widths typically at the percent level or below. Some decays, however, exhibit a stronger renormalization scheme and scale dependence. This implies a larger theoretical uncertainty and can be explained by the mass hierarchies and couplings governing these cases, which lead to loop-dominated decays.

Conclusions
In this paper we worked out the renormalization of the N2HDM, which is an interesting benchmark model for studying extended Higgs sectors involving Higgs-to-Higgs decays. For the mixing angles, we provided a renormalization scheme that is manifestly gauge independent by applying the alternative tadpole scheme combined with the pinch technique. We explained in great detail the notion of the alternative tadpole scheme in our renormalization framework, and for the first time provided the formulae for the pinched self-energies in the N2HDM. Apart from the addi-tional mixing angles as compared to the 2HDM, in the N2HDM we encounter a singlet VEV that needs to be renormalized as well. We elaborated in detail the implications of the alternative tadpole scheme for the renormalization of the singlet VEV that we renormalize through a physical quantity, given by a Higgs-to-Higgs decay width. The soft Z 2 breaking parameter m 2 12 , which, like v S , enters the Higgs self-couplings and hence features in Higgs-to-Higgs decays, is renormalized in the MS scheme. We studied the impact of our renormalization scheme by computing the EW one-loop corrections to various Higgs decay widths, including the Higgs decays into a massive Z-boson pair and into lighter Higgs pairs.
The computation of the EW corrections to our different sample decay widths has shown that the corrections can be sizeable and have to be taken into account in order to make reliable predictions for the Higgs observables. For a broad range of phenomenologically representative scenarios we find a rather weak renormalization scale and renormalization scheme dependence, indicative of a rather small theoretical error due to missing higher order corrections. In some cases the EW corrections can be sizeable, in particular if the corresponding LO decay widths are suppressed, so that the NLO-corrected width effectively becomes the leading order decay width. Higher order corrections beyond NLO would then be necessary in order to reduce the theoretical error.
With this paper, we have provided an important contribution to the renormalization of extended Higgs sectors involving singlet fields. This is crucial input for the computation of the EW corrections to the Higgs bosons of such models and therefore indispensable for the correct prediction and interpretation of Higgs observables at the LHC.
where Φ stands for an arbitrary neutral or charged scalar particle and m Φ i,j = 0 in case Φ i,j is a Goldstone boson. We introduce the one-loop integrals where A 0 , B 0 and C 0 denote the usual scalar one-, two-and three-point integrals and C 2 denotes the coefficient integral of the tensor integral C µ , which can be expressed solely through A 0 and B 0 integrals, cf. Refs. [49,61]. The index V denotes a vector boson V ∈ {W ± , Z, γ}.
In what follows, we extract the gauge dependences of all self-energies via the definition iΣ tad (p 2 ) = iΣ tad (p 2 ) where iΣ tad (p 2 ) is the fully gauge-dependent modified self-energy with tadpole contributions included, cf. Fig. 3, iΣ(p 2 ) g.d. represents the truly gauge-dependent part of the self-energy and iΣ(p 2 ) ξ V =1 denotes the evaluation of the self-energy in the 't-Hooft Feynman gauge. The inclusion of tadpole contributions for the analysis of the self-energies with respect to gauge dependence is necessary for a consistent application of the pinch technique [52]. While the extraction of the gauge dependence via Eq. (A.7) is not unique, we show in the following by applying the pinch technique that iΣ(p 2 ) g.d. is considered to be the truly gauge-dependent part of the self-energies, since it is precisely these terms which are cancelled by the pinch contributions.
A.1.1 Gauge dependence of the CP-even scalar self-energies First, we consider the gauge dependence of the CP-even scalar self-energies, i.e. the self-energies of all possible combinations of H i and H j (i, j = 1, 2, 3). All Feynman diagrams contributing gauge-dependent terms are shown in Fig. 13. The evaluation of Eq. (A.7) for the CP-even scalars of the N2HDM sector yields H i H j β ZZ (p 2 ) + β ZξZ (p 2 ) H i H j α W H i H j β W W (p 2 ) + β W ξW (p 2 ) , Figure 13: All Feynman diagrams contributing to the gauge dependence of the CP-even self-energies Σ tad H i H j (p 2 ). For the tadpole diagram, a sum over intermediate Higgs states H k (k = 1, 2, 3) is assumed. Note that the ghost and vector boson contributions in the tadpole diagrams precisely cancel against each other, so that these are not shown.
where the combinations O (1) H i H j and O (4) H i H j have been defined in Eq. (5.71). We note that when evaluating these combinations in the 2HDM limit, i.e. by applying Eq. (2.35), where O (4) H i H j reduces to the Kronecker delta δ H i H j , the result in Eq. (A.8) coincides with the results presented in Refs. [40,60] for the 2HDM as well as with the result presented in Ref. [59] for the MSSM, since the structure of the gauge-dependence of the CP-even scalar self-energies does not differ between the MSSM and the 2HDM.
A.1.2 Gauge dependence of the charged scalar and vector self-energies Next, we consider the charged sector. Due to the mixing of the charged particles of the N2HDM, we have to consider not only all possible self-energy combinations of the scalar particles H ± and G ± , but additionally their mixing with the charged vector bosons W ± . In the SM, where only one Higgs boson exists, it was shown that the Higgs contributions to the gauge dependence of the charged sector form a gauge-dependent subset which is cancelled by a corresponding subset of pinch contributions [56]. In the N2HDM we follow the same approach, i.e. we focus only on gauge-dependent contributions stemming from the enriched scalar sector of the N2HDM, which form a subset with respect to gauge dependence as well.
We first consider the gauge dependence of the self-energies of all combinations of W ± and G ± . The relevant contributions from the Higgs sector are given by the Feynman diagrams in Fig. 14 for all possible self-energies. Note that since we consider only the subset where the scalars of the N2HDM appear in the loops, only terms containing the gauge-fixing parameter Figure 14: All Feynman diagrams contributing to the gauge-dependence of the charged self-energies Σ tad where φ ± i,j ∈ {W, G ± , H ± }. A sum over intermediate Higgs states H k is assumed wherever they appear. Overlapping dashed and twiggled lines denote a scalar or a gauge boson, respectively, depending on the chosen particles. Note that we only consider contributions to the extended scalar sector of the N2HDM. Depending on the particles φ ± i,j chosen, some of the diagrams shown may not exist in the N2HDM. ξ W contribute to these self-energies. They explicitly read 16 iΣ tad W W,µν (p 2 ) = iΣ tad W W,µν (p 2 ) iΣ tad W G ± ,µ (p 2 ) = iΣ tad W G ± ,µ (p 2 ) ξ V =1 (A.10) (A.11) Next, the gauge dependence of the self-energies of all combinations of H ± and G ± or W ± is given by the relevant contributions from the Higgs sector as given by the Feynman diagrams in Fig. 14 as well. In the case of the self-energy for two H ± particles, additional dependences on 16 Note that in the case of the self-energy Σ tad G ± G ± we subtracted an additional term of f G ± G ± (p 2 )αW with respect to the diagrams shown in Fig. 14. This term stems from other gauge-dependent subsets of the gaugedependence of the self-energy, which we do not present explicitly here. This is in line with [56], where these additional terms are simply dropped since they cancel elsewhere. ξ Z and ξ γ appear even when focusing on the extended scalar sector of the N2HDM only, while The main idea of the pinch technique (cf. Refs. [51][52][53][54][55][56][57][58] for a detailed exposition) is to isolate the gauge dependences of an arbitrary toy scattering process, which features the to-be pinched self-energies in a unique way. This is achieved by applying the elementary Ward identities / k P L/R = S −1 1 (p + k)P L/R − P R/L S −1 2 (p) + m 1 P L/R − m 2 P R/L P L/R / k = P L/R S −1 1 (p + k) − S −1 2 (p)P R/L + m 1 P L/R − m 2 P R/L , (A. 17) where k denotes the loop momentum, m 1 and m 2 the masses of the external fermions of the considered toy process and S(p) the fermion propagator It turns out that the gauge dependences are all similar in structure, i.e. they are always selfenergy-like, independently of their origin within the scattering process. The isolation of all pinch contributions from the toy scattering process then allows for a manifestly gauge-independent definition of pinched self-energies. Since these self-energies are considered to be independent from the toy process chosen, cf. [52], the pinched self-energies are unique.
The application of these formulae to fermion-fermion-Higgs couplings enables the projection of the pinch contributions onto the desired CP-even Higgs couplings to the fermions, where "..." contains pinch contributions for other than the CP-even Higgs self-energies. Consequently, we can neglect them for the CP-even self-energies.
We consider the two contributions depicted in the Feynman diagrams in Fig. 17. The momenta are as defined in the diagrams. With the definitions given above, the sum of both diagrams Figure 16: All generic Feynman diagrams contributing to the CP-even pinched self-energies. between the external fermions and the scalar or vector particles of interest. In the case of the