Renormalization of mixing angles

We discuss the renormalization of mixing angles for theories with extended scalar sectors. Motivated by shortcomings of existing schemes for mixing angles, we review existing renormalization schemes and introduce new ones based on on-shell conditions or symmetry requirements such as rigid or background-field gauge invariance. Considering in particular the renormalization of the mixing angles in the Two-Higgs-Doublet Model and the Higgs-Singlet Extension of the Standard Model, we compare electroweak corrections within these models for a selection of renormalization schemes. As specific examples, we present next-to-leading-order results on the four-fermion decays of heavy and light CP-even Higgs bosons, $H_1/H_2\to WW/Z Z\to4f$, and on electroweak Higgs-boson production processes, i.e. Higgs-strahlung and vector-boson fusion. We find that our new proposals for on-shell and symmetry-based renormalization conditions are well-behaved for the considered benchmark scenarios in both models.


Introduction
After the discovery of a Higgs boson at the Large Hadron Collider (LHC) [1,2], the investigation of the Higgs sector is still of prime importance for particle physics. Theories with extended Higgs sectors typically contain additional scalar multiplets leading to physical scalar states that mix. Simple examples of such extensions are the Two-Higgs-Doublet Model (THDM) [3,4] and the Higgs-Singlet Extension of the Standard Model (HSESM) [5][6][7]. For a precise study of such theories, next-to-leading-order (NLO) QCD and electroweak (EW) corrections have to be taken into account. This requires a renormalization of these models and thus the renormalization of mixing angles or, more generally, of mixing matrices.
The need for renormalization of mixing matrices already appears in the Standard Model (SM) where the quark-mixing matrix has to be renormalized. While this is phenomenologically unimportant owing to the smallness of the down-type quark masses, the problem has nevertheless found quite some interest in the literature, and the corresponding theoretical developments have also influenced the work on the renormalization of mixing matrices in scalar sectors, which is the subject of this paper. A first renormalization condition for the quark-mixing matrix based on on-shell field-renormalization constants of the quark fields was proposed in Refs. [8,9]. This prescription is simple, symmetric in the fields that mix, and smoothly connected to the limit of degenerate quark masses. Later it was discovered [10] that the straightforward use of the renormalization condition of Refs. [8,9] gives rise to gauge-parameter-dependent counterterms for the quark-mixing matrix and thus to a gauge-parameter-dependent parametrization of S-matrix elements in terms of renormalized parameters. In the sequel, various proposals were made for a gauge-parameter-independent renormalization of the quark-mixing matrix [10][11][12][13][14][15]. Typically, these are cumbersome to apply, their generalization beyond one-loop order remains unclear, and/or they potentially lead to singularities in the S-matrix elements for degenerate quark masses. The last problem occurs, in particular, in the modified minimal subtraction (MS) scheme. A gaugeindependent, symmetric, physical renormalization condition was proposed in Ref. [16]. It was also suggested to define the quark-mixing matrix counterterm from the quark-field renormalization constants calculated in the 't Hooft-Feynman gauge [17]. Generalizing this idea, it was argued in Ref. [12] that any renormalization scheme for the quark-mixing matrix may be viewed as a gauge-invariant scheme by definition, in the sense that S-matrix elements remain invariant if the gauge used in the calculation of the loop corrections and all other renormalization constants is changed, while keeping the defining gauge for the renormalization constants of the quark-mixing matrix fixed.
The need for suitable renormalization schemes for mixing angles becomes more important in extensions of the SM. Specific examples are models with additional Higgs bosons, additional vector bosons, or additional fermions. In particular, for the renormalization of mixing angles in the scalar sector a variety of schemes were used in the literature. Specifically, the renormalization of the mixing angle β in the Minimal Supersymmetric Standard Model (MSSM) was discussed in Refs. [18][19][20][21][22].
The renormalization of the mixing angles α and β in the THDM was considered in Refs. [23][24][25][26][27][28][29]. Kanemura and collaborators [23,25,30] used the vanishing of the renormalized non-diagonal "on-shell" scalar 2-point functions to fix the mixing-angle renormalization. 1 Despite the choice of "on-shell" momenta these conditions do not derive from S-matrix elements and are gauge dependent. In Ref. [24], the mixing angle α was fixed by the condition that the mixing self-energy of the CP-even Higgs bosons vanishes "on-shell", while β was renormalized requiring that the ratio of vacuum expectation values (vevs), v 2 /v 1 , is expressed in terms of the "true" vacua following the treatment of Refs. [18,19] in the MSSM. The authors of Ref. [26] employ the renormalization conditions of Ref. [23] within the FJ Tadpole Scheme by Fleischer and Jegerlehner [31] and define gauge-independent counterterms based on the "pinch-technique" prescription [32]. 2 In the series of papers [27,28] the gauge dependence of the MS definition of mixing angles with respect to different tadpole counterterm schemes was investigated, and predictions were compared against before-mentioned schemes based on mixing energies. Finally, in Ref. [29] new (gauge-independent) MS schemes were introduced, replacing a mixing-angle definition by the MS renormalization of a coupling parameter of the Higgs potential. scalar mixing energy or MS renormalization. While in Ref. [28] an MS scheme was compared with schemes based on on-shell self-energies, in Refs. [37] different MS schemes were studied.
The purpose of this paper is a discussion of renormalization prescriptions and schemes for mixing angles in general scalar sectors of gauge theories and a comparison of different schemes in concrete phenomenological applications in the THDM and the HSESM. We critically review existing renormalization prescriptions and introduce new ones that exhibit several desirable properties [20]. In particular, we introduce genuine on-shell renormalization conditions for mixing angles based on combinations of suitable S-matrix elements and put renormalization schemes based on symmetry requirements such as rigid invariance or background-field gauge invariance on a general footing.
This paper is organized as follows. In Section 2 we specify some useful definitions and conventions. In Section 3 we review and introduce renormalization schemes for mixing angles in scalar sectors. We begin with a discussion of existing MS renormalization schemes, followed by sections where we construct new on-shell renormalization schemes and renormalization schemes based on symmetries. In Section 4 we provide a numerical discussion of renormalization schemes in applications to Higgs decays into 4 fermions as well as Higgs production at the LHC via Higgs-strahlung or weak vector-boson fusion. After the conclusion in Section 5, we give translation rules of our conventions to other formulations in the literature in App. A. Further appendices provide explicit analytical results for quantities used in the various renormalization schemes, including scalar self-energies in the background-field method, the tadpole contributions to the scalar self-energies, vertex corrections for on-shell schemes, and a discussion of background-field Ward identities in different tadpole counterterm schemes.

Renormalization transformation for mixing fields
As specific examples for the renormalization of mixing angles we consider the mixing of scalar fields in the THDM and the HSESM. Both theories involve two physical CP-even scalar bosons. Let the corresponding fields in the symmetric basis be η 1 and η 2 and the fields in the physical mass-eigenstate basis H 1 and H 2 . The fields are related via a rotation where we use the shorthand notations c α = cos α and s α = sin α. In the general case of more than two mixing fields H i , the matrix R(α) depends on a set of mixing angles α = {α i }.
Performing the renormalization in the physical basis in the complete on-shell scheme [9,38], the renormalization transformations for the mixing angle and the field renormalization constants of the scalar fields read where bare quantities carry an index B, and δα and δZ H ij represent the renormalization constants for the mixing angle α and the scalar fields corresponding to mass eigenstates. In the complete on-shell scheme of Refs. [9,38], the non-diagonal field-renormalization constants are fixed as where M H i and M H j denote the masses of the corresponding scalar bosons H i and H j , and Σ ij their mixing energy. In this paper we consistently use the convention that the self and mixing energies include explicit and implicit tadpole contributions, i.e. they are defined as the higher-order contributions to the inverse propagators which at one-loop order can be depicted as The first contribution is the bare loop one-particle-irreducible (1PI) energy and the second term the corresponding 2-point tadpole counterterm. The third and forth terms are the explicit tadpole loop and tadpole counterterm contributions, respectively. Note that no counterterms other than the ones from the tadpoles (such as mass or field renormalization constants) are included in Σ ij . In the renormalization of the SM Higgs sector described in Refs. [9,34], the vev v is renormalized in such a way that the Higgs field has vanishing vev, and consequently, all explicit tadpole contributions vanish, i.e. the third and forth term in (2.6) add up to zero. However, implicit tadpole counterterms may remain in the second term, but their explicit form depends on the considered fields and on the tadpole counterterm scheme in use. If we use self-energies without any tadpole contributions, we will indicate this with an extra index 1PI for one-particle irreducible. If not stated otherwise, the formulas are given in the Fleischer and Jegerlehner tadpole counterterm scheme as used in Ref. [27] for which the construction of the implicit tadpole contributions to 2-point functions is given in App. C. This convention for the tadpole contributions in Σ ij is in agreement with the one of Refs. [31,39].
In the presence of mixing, the well-known problem of degenerate states in time-independent perturbation theory of quantum mechanics also appears in quantum field theory and becomes first apparent in the one-loop renormalization. The condition (2.5) shifts corrections from the mixing of states that are induced by external (non-diagonal) self-energy insertions to vertex counterterms. For degenerate masses, M H i → M H j , the constants δZ H ij become singular, and thus also the S-matrix elements with external H i , H j fields unless this singularity is cancelled by some other contribution. As the loop diagrams are regular in this limit, the cancellation should come from another counterterm. This is where the mixing angles come into play, and one expects that an appropriate renormalization of the mixing angles can make the S-matrix well-behaved as we will sketch in the following. 3 Combining (2.1) for bare fields with (2.4) and using yields to one-loop accuracy Thus, if the mixing matrix R in the Lagrangian results only from the rotation between symmetric and mass eigenstates as in (2.1), the counterterm to the mixing matrix appears only in the combination in S-matrix elements. For the specific case of a single mixing angle this gives rise to the two combinations −δα (2.10) Choosing the mixing-angle counterterm appropriately allows one to cancel the singularity for degenerate masses in S-matrix elements arising from the denominators in (2.5).
If the mixing angle α is promoted to a physical parameter upon eliminating some bare parameter of the Lagrangian in its favour, the counterterm δα appears also independently of δZ H 12/21 in S-matrix elements. This is the case in the THDM and the HSESM if a quartic coupling parameter λ i of the Higgs potential is eliminated in favour of the mixing angle α in the sector of CP-even Higgs bosons. However, it turns out that in this case the mixing-angle counterterm appears in addition to (2.10) only in the combination (M 2 H 1 − M 2 H 2 )δα. This can easily be verified by considering the coupling parameter λ i as a function of Higgs masses and mixing angles (see, for instance, Eq. (3.10) in Ref. [27] or Eq. (2.20) in Ref. [29]). This statement can be generalized to mixing angles in more general theories. Note, however, that the situation is different for mixing angles defined via vevs, such as β in the THDM or the MSSM.

Higgs-Singlet Extension of the Standard Model
In the HSESM we associate η 2 in Eq. (2.1) with the scalar component of the Higgs doublet where v 2 is the corresponding vev and φ + , χ the would-be Goldstone-boson fields. The field η 1 is the field excitation of the (canonically normalized) Higgs singlet field which acquires the vev v 1 . The Higgs potential of the considered variant of the HSESM is given by [5][6][7]40] which possesses a Z 2 symmetry under σ → −σ. Translation rules of these conventions to the ones used in Refs. [28,37] are given in App. A. The masses M W and M Z of the weak gauge bosons are given by where g 2 and g 1 are the SU(2) and U(1) gauge couplings, respectively.

Two-Higgs-Doublet Model
In the THDM, η 1 and η 2 are the neutral CP-even scalar components of the two Higgsdoublet fields which has a Z 2 symmetry w.r.t. Φ 1 → −Φ 1 and Φ 2 → Φ 2 that is only softly broken by the m 2 12 term. The two scalar doublets in the THDM contain, in addition to the scalar fields η i , also pseudoscalar χ i and charged scalar φ + i fields. These are transformed to the physical basis as follows, (2.17) Here G ± and G 0 are the charged and neutral would-be Goldstone-boson fields, and H ± and A 0 the physical charged and pseudoscalar Higgs fields, respectively. The mixing angle β is related to the vevs v i of the two scalar doublets via t β ≡ tan β = v 2 /v 1 , and we use the abbreviations c β = cos β and s β = sin β. The masses M W and M Z of the weak gauge bosons are given by The renormalization transformations for the mixing angle β and the pseudoscalar fields in the complete on-shell scheme can be written as where δβ and δZ ... represent the renormalization constants for the mixing angle β and the physical pseudoscalar A 0 and the would-be-Goldstone field G 0 . Similar equations can be written for the charged scalar fields.

Renormalization schemes for mixing matrices
In this section we review existing renormalization schemes for mixing angles and propose and discuss new ones. In Ref. [20] desirable properties for the renormalization of mixing matrices were formulated: • The mixing-angle renormalization should be gauge independent, i.e. renormalized Smatrix elements should be gauge-independent functions of the renormalized mixing angles.
• The mixing-angle renormalization should be symmetric with respect to the mixing degrees of freedom. Moreover, the renormalized mixing angle should be independent of a specific physical process.
• The mixing-angle renormalization should not spoil the numerical stability of the perturbative expansion; in particular, the running of parameters and radiative corrections to physical observables should be accessible via perturbation theory.
We add a further condition, which could be viewed as a refinement of the third condition of Ref. [20]: • In the limit of degenerate masses of the mixing particles or in the limit of extreme mixing angles, no singularities should be introduced in physical observables, i.e. Smatrix elements should behave smoothly in these limits. Furthermore, there should be no "dead corners" in the parameter space of the model where a renormalized input parameter nominally goes to infinity. 5 Focusing on these requirements, as far as possible, we propose new renormalization schemes, specifically for, but not limited to the THDM and the HSESM, and compare NLO predictions for some processes obtained with some old and the new schemes. We write the relations between renormalization constants and self-energies without taking real parts leading to renormalization constants with imaginary parts. This is appropriate for the complex-mass scheme [42,43]. In the usual on-shell scheme, the real part of the self-energies should be taken in all renormalization conditions.

Renormalization of mixing angles in MS schemes
A straightforward, universal renormalization scheme, which does not distinguish a specific mass scale in the case of the renormalization of a mixing angle, is provided by MS renormalization, where the renormalization constants contain only ultraviolet(UV)-divergent parts along with some universal finite constants, i.e. the combination in dimensional regularization, where D is the space-time dimension and γ E the Euler-Mascheroni constant. MS renormalization can be straightforwardly applied to arbitrary mixing matrices. The MS renormalization of mixing angles is, by construction, symmetric in the fields that mix and does not depend on a specific observable. Since mixing angles are defined in the physical basis, their MS renormalization depends on the precise treatment of tadpoles [27,31]. If tadpoles are treated in a conventional way, i.e. if renormalized tadpoles are set to zero in the course of parameter renormalization and tadpole counterterms partially absorbed into bare masses (see, e.g., the renormalization of the SM in Refs. [9,39]), the counterterms for the mixing angles become gauge dependent. The tadpole scheme based on Ref. [9] where the bare masses are defined as the coefficients of the terms quadratic in the physical fields in the Lagrangian is referred to as PRTS (Parameter-Renormalized Tadpole Scheme) in the following. 6 We denote the MS scheme applied to mixing angles based on the PRTS as MS(PRTS) in the following. If tadpoles are treated in the FJ Tadpole Scheme (FJTS) [26,27,31], i.e. removed by a suitable field redefinition, the resulting MS renormalization scheme applied to mixing angles, denoted in the following as MS(FJTS), is by construction gauge independent. The MS(PRTS) and MS(FJTS) schemes have been worked out for the THDM and the HSESM in different variants in Refs. [26][27][28][29] and [28,37], respectively. The "Tadpole scheme" for the renormalization of tan β in the Minimal Supersymmetric Standard Model studied in Ref. [20] is equivalent to the MS(FJTS) scheme for the mixing angle β. It was found that this scheme leads to an unacceptably large scheme uncertainty in the renormalization of tan β. In Ref. [28] the MS(FJTS) scheme has been applied to the renormalization of mixing angles α and β in the THDM and the HSESM in an NLO study of heavy and light Higgs production in Higgs-strahlung, pp → H 1,2 µ − µ + + X, and vector-boson fusion, pp → H 1,2 jj + X. The results of the MS(FJTS) scheme turned out to be unstable and suffering from large scale uncertainties in many scenarios, while results in schemes based on on-shell self-and mixing energies remained well-behaved. NLO results obtained with the MS(PRTS) and MS(FJTS) schemes were compared in quite some detail for the light Higgs-boson decays H 2 → WW/ZZ → 4 fermions in the THDM [29,41] and the HSESM [37]. While the results of the two schemes are perturbatively stable in the HSESM, with nice plateaus showing up in the variation of the renormalization scale at NLO, the results on H 2 → 4f in the MS(FJTS) scheme lead to serious perturbative instabilities in THDM scenarios that are away from the alignment limit, where | cos(β − α)| is not small.
All MS renormalization schemes for mixing angles give rise to large corrections in the limit of degenerate masses. This enhancement results from terms of the form (2.10) in Smatrix elements. In MS schemes, δα cancels only the UV-divergent parts, but the remaining UV-finite terms in (2.10) resulting from the field (or wave-function) renormalization become singular for M H i → M H j . The size of these terms, which are also present in the MS(PRTS) scheme, is enhanced by additional tadpole contributions in the MS(FJTS) scheme. While for tan β similar enhancements due to additional tadpole contributions take place, the limit M H i → M H j does not introduce singularities connected with the renormalization of β, so that no corresponding singular contributions to S-matrix elements can result. In this sense, the situation for true mixing angles (such as α) is more involved.
Instead of imposing the MS condition on the mixing angles, MS renormalization can be applied directly to parameters of the Higgs potential. This idea was, e.g., pursued in the "λ 3 " schemes of Refs. [29,41] as an alternative to the MS renormalization of the mixing angle α in the THDM; specifically, α was replaced as input parameter by the coupling λ 3 , which was MS renormalized. While such a renormalization condition is gauge independent and does not lead to singularities for degenerate masses, problematic regions ("dead corners") in the parameter space show up where the relation between the distinguished coupling and the mixing angle cannot be inverted. In the "λ 3 " schemes of Refs. [29,41], for instance, a singularity in the parametrization of observables by λ 3 occurs in scenarios in which cos(2α) → 0. To circumvent this problem in the THDM, one would have to patch the parameter space by switching from λ 3 to another scalar coupling as renormalized parameter.
In summary, MS renormalization schemes for mixing angles have some desirable properties (simplicity, symmetry, process independence), but suffer, in general, from problems with perturbative stability in certain parameter regions, such as for mass degeneracy of the mixing fields. Moreover, care has to be taken in view of gauge dependence. On the other hand, it should be mentioned that MS renormalization offers a simple way to estimate perturbative stability by varying the renormalization scale in predictions and checking for a stabilization of results in the transition from leading order (LO) to NLO.

Physical (on-shell) renormalization conditions for mixing angles
The renormalization of mixing angles can be directly fixed from observables or S-matrix elements that depend on these mixing angles at LO. Such on-shell renormalization conditions are evidently gauge independent. 7 Fixing mixing angles by specific processes, however, has a number of potential disadvantages: • By construction, on-shell renormalization conditions are process dependent and often destroy the symmetries between the particles that mix.
• On-shell conditions are only directly applicable to S-matrix elements that do not involve charged particles; otherwise the counterterms would become infrared (IR) singular. 8 The problem of IR singularities can, in principle, be avoided by imposing the renormalization condition on a full physical observable, e.g. by demanding that a partial decay width does not receive any correction, but this procedure shifts processspecific real-radiation effects into the renormalization constant.
• Typically, observables and S-matrix elements depend not only on a mixing angle, but on other parameters as well. Upon defining the mixing-angle renormalization from such quantities, one thus absorbs corrections to the considered observable or S-matrix element into the mixing-angle counterterm that are related to other parameters of the model. This can be a source for unnaturally large corrections. In Ref. [26] it was, e.g., demonstrated that on-shell renormalization conditions based on specific observables lead to numerically unstable results in the THDM.
The situation can be improved by considering combinations of physical observables or S-matrix elements that depend exclusively on a specific mixing angle and on no other parameters, so that renormalization contributions of other parameters or normalization effects systematically drop out. For the quark-mixing matrix such a renormalization scheme was proposed in Ref. [16].
In order to fix the renormalization of the Higgs mixing angle α introduced in (2.1), we consider a set of processes involving the fields H i that have a simple dependence on the mixing angle α in LO. If the dependence on α in the considered observables only results from the transformation (2.1), this is typically the case. 7 Note that we do not consider renormalization conditions as "on shell" that are based on mixing energies, Green functions, or formfactors as well as "matrix elements" involving unphysical degrees of freedom at some "on-shell" configurations of momenta. This includes, in particular, mixing energies of scalar bosons, of would-be Goldstone bosons with Higgs bosons or of would-be Goldstone bosons with gauge bosons. In the literature, schemes of this kind are often called "on shell" as well. 8 Considering the analytic structure of one-loop 3-point functions, it can be seen that the corresponding IR singularities in decay S-matrix elements can only be cancelled by those of the field-renormalization constants of the external charged particles if one of the three external particles is massless and neutral, as for the photon in the electron-positron-photon vertex in QED. In other words, potential IR singularities in the S-matrix element used to fix a renormalization constant only cancel in this very specific case. In all other cases, IR singularities would enter the parameter renormalization constants.

Higgs-Singlet Extension of the Standard Model
As a first specific example, we choose the decay of scalar bosons H 1 and H 2 into pairs of Z bosons in the HSESM. 9 The LO vertices read where δ H i ZZ = δ H i ZZ (M 2 H i ) are the unrenormalized relative one-loop corrections to the respective decays, and the counterterms have been written explicitly. In particular, the two scalar fields are renormalized according to (2.4). Inserting (3.5) into (3.4) and expanding to NLO yields for the counterterm of the mixing angle α: Using this counterterm, the renormalized NLO matrix elements become Since all differences have been incorporated in the renormalization of α, the relative corrections to the two different Higgs-boson decays become equal. The renormalization condition (3.4) has several desirable properties: • It is gauge independent, because it is based on physical S-matrix elements; • it is symmetric with respect to the scalar fields H 1 and H 2 ; • it is numerically stable for degenerate masses M H 1 ∼ M H 2 ; • it has smooth limits for extreme mixing angles, i.e. for c α → 0 or s α → 0.
which is finite for degenerate masses, M H 1 → M H 2 , and moreover all momentum-independent contributions to the mixing energy, such as tadpole contributions, cancel therein. These statements hold for all other on-shell schemes for α discussed in this section. Still, the condition (3.4) has some disadvantages: • It can be directly applied to decay processes involving only electrically neutral external particles. For charged external particles, the renormalization constant δα becomes IR singular.
• Depending on the masses of the external particles, the form factors have to be evaluated at phase-space points in the unphysical region. 10 This is, for instance, the case in (3.6) if H 1 or H 2 is identified with the observed Higgs state of mass 125 GeV.
All these drawbacks can be lifted upon introducing extra neutral fields with a simple coupling structure that allows us to fix the renormalization of the mixing angles while recovering the original theory upon sending the extra couplings to zero.
As an example for this procedure, we consider the HSESM and add an additional fermion singlet field ψ with the Lagrangian (3.10) 10 At least for decays at the one-loop level this does not constitute an obstacle. In this case the relevant 3-point functions can be analytically continued to the unphysical region as discussed in Ref. [44].
In terms of scalar fields in the basis of mass eigenstates, this becomes Considering the limit of a vanishing Yukawa coupling, y ψ → 0, we recover the original HSESM with an additional massless fermion ψ (m ψ = y ψ v 1 → 0), which completely decouples from all other particles, i.e. we effectively recover the original theory. We require that the ratio of matrix elements for the decays of the two scalar Higgs bosons into a ψψ pair of singlets is equal to its leading-order value −c α /s α in the limit of vanishing coupling y ψ : where the proportionality factor, which is not spelled out, only contains the ratio of spinor chains (with different kinematics), but no other model parameters. Since the ratio is based on matrix elements for the decay of massive scalars into massless neutral fermions these are in the physical region, and no IR singularities occur. Moreover, since all NLO vertex corrections in M H i →ψψ tend to zero at least quadratically in y ψ and thus drop out in the limit y ψ → 0 11 , we obtain for the counterterm Thus, owing to the simple structure of the model, all vertex corrections drop out, and the mixing-angle counterterm is fixed by a gauge-independent combination of field-renormalization constants only. Note also that the spinor chains suppressed in (3.12) do not enter the final result for δα, since they cancel in the ratios

Two-Higgs-Doublet Model
In order to formulate on-shell renormalization conditions for the THDM, we add two righthanded fermion singlets to the Lagrangian: ν 1R transforming under the extra Z 2 symmetry as ν 1R → −ν 1R and ν 2R transforming as ν 2R → ν 2R , so that ν iR can only receive a Yukawa coupling to Φ i . The additional Lagrangian reads where y ν i are new Yukawa couplings that are considered in the limit y ν i → 0. The fields L iL = (ν i , l i ) T L are left-handed lepton doublets of the SM, say the electron-neutrino and muon-neutrino doublets, Φ i are the two Higgs-doublet fields, and σ 2 the second Pauli matrix. Upon inserting the representations of the doublet fields this leads to
where we suppressed terms involving charged scalar Higgs and would-be Goldstone-boson fields. For non-zero couplings y ν i , the neutrinos ν 1 and ν 2 correspond to massive Dirac fermions without generation mixing owing to the conserved lepton number. In the limit y ν i → 0, the neutrinos become massless, and the right-handed parts decouple.
Renormalization of α The renormalization of the angle α can be fixed upon requiring that the ratio of the matrix elements for the decays H i → ν 1ν1 is the same at LO and NLO, This leads to the renormalization constant where δ H i ν 1ν1 represent the unrenormalized relative one-loop corrections to the decays H i → ν 1ν1 . Alternatively, if the ratio of the matrix elements for the decays H i → ν 2ν2 is used to fix the renormalization of α, this leads to the renormalization constant The explicit form of the vertex correction factors δ H 1 ν 1ν1 , etc., is given in App. D. At NLO, these loop corrections respect the chiral structures of the respective underlying LO couplings; beyond NLO this might not be the case anymore, so that one would have to write the renormalization conditions in terms of the form factors that correspond to the LO couplings.
Renormalization of β The renormalization of the angle β can be fixed by demanding that the ratio of the matrix elements for the decays H 1 → ν iνi and A 0 → ν iνi is the same at LO and NLO for one of the two neutrinos ν i . Specifying ν i to ν 1 , means where again the suppressed proportionality factor is given by a ratio of spinor chains, but does not contain further model parameters. This results in where δ A 0 ν 1ν1 denotes the unrenormalized relative one-loop corrections to the matrix element for the decay A 0 → ν 1ν1 . Upon inserting the renormalization constant δα from (3.17) this becomes which involves only vertex corrections for neutrino ν 1 . Note that we would get the same relation if we had fixed δβ from the ratio of A 0 → ν 1ν1 and H 2 → ν 1ν1 matrix elements by virtue of (3.16). Alternatively, δα and δβ can be fixed from analogous matrix elements involving only neutrino ν 2 . In this case, for δβ we demand which can be further processed with δα from (3.18) to yield The conditions (3.21) and (3.24) become singular for c β → 0 or s β → 0, respectively. Since in the phenomenological applications of the THDM, c β and s β are always non-vanishing, this does not lead to a singularity but can cause artificial enhancements. The renormalization scheme based on the conditions (3.17) and (3.21) involving only ν 1 is called OS1 in the following, the one based on (3.18) and (3.24) involving only ν 2 is called OS2.
The conditions (3.19) and (3.22) do not directly apply to β, but to a combination of α and β. A condition that fixes β directly can be obtained by using the decays into both neutrino singlets and requiring Note that the spinor chains all cancel within the multiple ratio. This condition leads to the counterterm which is regular for s β → 0 of c β → 0, but potentially singular for s α → 0 or c α → 0.
A condition leading to a counterterm that is regular in all these limits can be constructed upon using linear combinations of observable quantities from different processes.
To this end, we consider on-shell form factors instead of complete matrix elements. For the scalar decays into (nearly) massless fermions we write the matrix elements as whereū ν and v ν are the spinors of the final-state fermions and antifermions, and [. . . ] H i /A 0 indicates the decay kinematics of the spinor chain. The functions F H i →ν jνj and F A 0 →ν jνj denote the formfactors for the decays into neutrinos ν j of type j = 1, 2. The LO formfactors follow directly from the Lagrangian (3.15): As renormalization condition we require that the following LO relation holds also at higher orders: This fixes the counterterm for the mixing angle β to This result is non-singular in all limits s α → 0, c α → 0, s β → 0, or c β → 0.
The above on-shell renormalization conditions for mixing angles in the HSESM and THDM depend on the introduction of specific auxiliary fields. While this method can in principle be generalized to more general theories, we are not able to provide a simple specific recipe for the on-shell renormalization of mixing angles or mixing matrices in general. As far as we can see, this has to be investigated anew for each theory, but the shown examples can certainly serve as guidelines.

Rigid symmetry and wave-function renormalization for physical states
The renormalization of mixing matrices can be related to the wave-function renormalization of external fields upon using rigid symmetry, i.e. the symmetry under global gauge transformations, of the Lagrangian. Again, we specifically consider the THDM and the HSESM, where the CP-even scalar fields in the symmetric and mass-eigenstate bases are related via (2.1).
A theory with a spontaneously broken gauge symmetry can be renormalized using the renormalization constants for fields and dimensionless parameters from the symmetric phase [45][46][47][48][49][50][51]. In particular, the counterterms for the dimensionless parameters and fields can be directly taken from the symmetric formulation. Once the counterterms for the mass parameters are adjusted appropriately, all UV divergences cancel while the external lines in S-matrix elements still require additional finite wave-function renormalization. For the SM such a renormalization scheme was used in Ref. [52]. In the case considered here, the relevant renormalization transformations read where the components η 1 and η 2 of η belong to different multiplets of the gauge group. Consistency of Eqs. (2.1) and (3.33) requires that at least the divergent parts of the renormalization constants in both schemes are related via This implies Thus, we find, in particular, a relation between the renormalization constant of the mixing angle δα and the non-diagonal field renormalization constants of the scalar Higgs-field pair. While the relations (3.37)-(3.40) hold for the UV-divergent parts, as for instance discussed in Ref. [29], not all of them can be required simultaneously for the finite parts if the field renormalization is fixed in the complete on-shell scheme.
On the other hand, (3.40) can be used to fix the renormalization of the mixing angle α in terms of on-shell field renormalization constants of the scalar Higgs fields, i.e. we can define where we expressed the non-diagonal field renormalization constants by the non-diagonal mixing energy upon using the on-shell renormalization conditions (2.5). This renormalization condition has been introduced by Kanemura et al. in Ref. [23] and used in Ref. [26] both in the tadpole scheme of Ref. [23] and the FJTS. As discussed in Section 2, the counterterm to the mixing angle appears, besides in the regular expression (M 2 H 1 − M 2 H 2 )δα, only in the combinations (2.10) in S-matrix elements. Thus, when using the renormalization condition (3.41), only the combination remains. According to (3.9), this is finite for degenerate masses, M H 1 → M H 2 , and moreover all momentum-independent contributions to the mixing energy, such as tadpole contributions, cancel therein. Renormalizing mixing angles through appropriately chosen field renormalization constants has several advantages: • The symmetry between the different states is respected; • if the limit of vanishing mixing is protected by a symmetry implying Σ H 12 → 0 for α → 0, this is not violated by the renormalization condition, i.e. δα → 0 for α → 0; • there is a smooth limit for degenerate masses, i.e. the renormalization is numerically stable and does not lead to enhanced corrections; • there is no problem with IR singularities since the mixing energies are free of such contributions.
An apparent drawback is the gauge dependence of the field renormalization constants. However, we can choose a specific gauge to calculate the counterterms for the mixing angles and fix these counterterms in this gauge. 12 Then, we can vary the gauge as usual (but keeping δα fixed in the original gauge), and S-matrix elements are gauge independent for fixed δα. It remains to pick a suitable gauge to fix the mixing-angle counterterm. In order not to introduce artificially large parameters, the 't Hooft-Feynman gauge (ξ = 1) can be chosen. This can be done in the conventional formalism or, preferably, in the backgroundfield method (BFM) where rigid symmetry holds for the effective action by construction (see Section 3.3.3). Once the mixing-angle renormalization is fixed in this way, relations between observables can be calculated as usual. The gauge independence of S-matrix elements ensures that no singularities appear if the calculation is done in a different gauge. This procedure relies on the fact that unrenormalized S-matrix elements with appropriate LSZ factors are gauge independent as functions of the bare parameters. This requires that all relations between bare parameters are gauge independent which is the case in the FJTS scheme, but not in the tadpole schemes of Refs. [9,39].
For the sake of clarity, let us assume that the mixing angle counterterm δα(ξ 0 ) has been fixed from δZ H ij (ξ 0 ) in a specific gauge, with fixed gauge parameter ξ 0 . If the S-matrix elements are calculated in a different gauge (with gauge parameter ξ = ξ 0 ), but with the fixed counterterm δα(ξ 0 ), the cancellations of potential singularities for degenerate masses occurring in (3.42) and (3.9) are not obvious anymore. The contributions of such terms can be studied by considering the gauge dependence of (3.41). The gauge-dependent parts of the one-loop mixing energy obey the Nielsen identity [53] ij are one-loop Green functions involving the operators for the BRST transformations of the fields H i . This implies, i.e. the gauge-dependent part is regular for M H i → M H j and free of momentum-independent contributions. Using (2.5) and (3.35) for higher-dimensional matrices, the renormalization condition (3.41) can be straightforwardly generalized to models with more physical scalar fields as which fixes the finite parts in δα by definition. Moreover, (3.35) implies In the light of the discussion of this section, the original proposal for the renormalization of the quark-mixing matrix [8] in the SM turns out to be viable. Strictly speaking, this requires the use of the gauge-independent FJTS scheme. However, since there are no tadpole contributions to the quark mixing energies and thus to the non-diagonal fermion field renormalization constants in the SM, the renormalization condition of Ref. [8] for the quarkmixing matrix is not affected by the tadpole scheme. Thus, after fixing the counterterms for the quark-mixing matrix in the 't Hooft-Feynman gauge, (gauge-independent) observables can be calculated as usual. Moreover, for the relevant one-loop fermion self-energies the results of the BFM are equivalent to those of the conventional formalism.

Rigid symmetry and wave-function renormalization for unphysical states
The renormalization condition (3.41) can also be used for the mixing with the wouldbe Goldstone bosons as for instance in the pseudoscalar and charged-scalar sector of the THDM. Using the fact that would-be Goldstone bosons are massless upon disregarding the gauge-fixing term that is not renormalized in linear gauges, we obtain for instance for mixing of the pseudoscalar with the would-be Goldstone boson, where M 2 is an arbitrary scale at which the G 0 A 0 mixing energy Σ G 0 A 0 is required to vanish. Choosing, for instance, M = 0 leads to (3.48) and [in analogy to (3.41)] In the BFM (see Section 3.3.3) this simplifies owing to the Ward identity 13 This expression formally coincides with the renormalization of β proposed in Ref. [23] within the Landau gauge of the conventional formalism. We stress that the simple form (3.50) for the counterterm δβ holds only in special gauges, while in general (3.49) has to be used. Moreover, (3.50) requires the PRTS scheme of Ref. [9], while in the FJTS scheme extra tadpole contributions appear [see Eq. (3.82) below].

Background-field gauge invariance
The method of the previous sections is not applicable to the renormalization of parameters that are not directly related to the mixing of fields. An example is the renormalization of the singlet sector of the HSESM. Besides the mass of the second Higgs scalar, a second parameter has to be selected. One can choose one of the quartic couplings λ 1 or λ 3 of the Higgs-singlet field, its vev v 1 , or the quantity tan β = v 2 /v 1 , in analogy to the THDM. In such cases, the background-field gauge invariance that appears when quantizing in the BFM can be exploited to fix renormalization constants. The BFM (see, e.g., Refs. [54][55][56]) for the EW SM was introduced in Ref. [34], and its application to the THDM and the HSESM was described in Ref. [28]. Following these references we denote background fields with carets. All relations discussed in Sections 3.3.1 and 3.3.2 hold in the BFM for the corresponding quantities of background fields as well.
In the conventional formalism, the gauge-fixing term breaks rigid invariance. 14 In the BFM, rigid gauge invariance is maintained for the background fields upon choosing an appropriate gauge-fixing term for the quantum fields [28,34]. Rigid invariance for the background field gives rise to Ward identities (c.f. App. E) and restrictions on the renormalization constants for the background fields. In particular, in the BFM the relations (3.37)-(3.40) can all be maintained including the finite parts. The relations resulting from rigid gauge invariance in the BFM were presented for the SM in the renormalization scheme 13 For the validity of (E.27) it is crucial that the self-energy is the quantity appearing in the inverse of the propagator, i.e. that it includes tadpole contributions and corresponding counterterms.
14 This does not affect the discussion in Sections 3.3.1 and 3.3.2, since it relies only on the fact that a symmetric field renormalization is possible in the spontaneously broken phase.
of Ref. [9] in Ref. [34]. In the FJTS scheme, as defined in Refs. [27,31], additional tadpole contributions ∝ ∆v appear in the Ward identities and thus also in these relations. For instance, using the conventions of Ref. [28] the last line of Eq. (46) of Ref. [34] for the SM in the FJTS scheme becomes where the extra term, the shift ∆v in the vev, is related to the tadpole by δtĤ denotes the tadpole counterterm and TĤ the unrenormalized tadpole, i.e. the renormalized Higgs one-point vertex function is given by ΓĤ = (TĤ + δtĤ ) = 0. The counterterms δZ e and δc 2 w are fixed according to (48) and (46) of Ref. [34] and can be shown to be independent of the tadpole counterterm scheme. In the FJTS scheme the W-boson mass counterterm gets additional implicit tadpole contributions ∝ ∆v as described in Ref. [27] which is equivalent to including explicit tadpoles in the self-energy that determines the counterterm, i.e.
As a consequence, in the FJTS scheme the W-boson mass counterterm is gauge independent, and the gauge dependence of δZη matches the one of ∆v. In the renormalization scheme of Ref. [9] the tadpole contribution in (3.53) is absent, and the mass counterterm is gauge dependent (since the definition of the bare mass involves the shifted vev, see the discussion in Ref. [28]). Using as well as the renormalization transformations of the SM parameters, the relation (3.51) implies Thus, the vev is renormalized as the corresponding scalar doublet field apart from the explicit shift ∆v introduced to ensure a vanishing tadpole at one-loop order.
In the rest of this section we use the FJTS scheme.
Higgs Singlet Extension of the Standard Model In the HSESM within the BFM, (3.51) and (3.56) still hold with the components of the SM doublet replaced by those of the doublet of this model, where (3.58) In addition, we obtain a corresponding relation between the renormalization of the component η 1 of the singlet field σ and the corresponding vev v 1 , The shifts of the vevs can be expressed in terms of the tadpole counterterms as Fixing δZη 2 from (3.57), using the on-shell field renormalization constants (2.5) to determine δZη 1 from (3.39) including the finite parts, where δZĤ ij are the elements of the field renormalization matrix for the scalars as defined in (2.4), we can use (3.59) to fix the renormalization of the singlet vev. This results in Defining this translates to where Alternative definitions of the counterterms can be obtained upon using (3.37) and (3.38) instead of (3.39) to fix (δZη 1 − δZη 2 ). In particular, upon using an appropriate linear combination of (3.37), (3.38) and (3.39), one finds Using this expression avoids the potential singularity for c α → 0 or s α → 0 in (3.64) and leads to Two-Higgs-Doublet Model Considering now the THDM and using again the BFM and the corresponding background-field gauge invariance, we obtain instead of (3.51): where the shifts of the vevs can be expressed via the tadpole counterterms as The terms involving shifts in the vevs in (3.74) can be expressed by the tadpoles as (3.75) 15 Note that in the first preprint version of this paper the scheme BFMS in the HSESM was based on (3.64) and (3.65). For the scenarios considered in Section 4 the differences between the two choices in the numerical results are marginal.
As for the HSESM, alternative counterterms can be obtained using (3.37) and (3.38) instead of (3.39) to fix (δZη 1 − δZη 2 ). If s α or c α become small, the counterterm defined by (3.74) becomes artificially large, which is actually the case in the THDM scenarios B1 and B2 considered in Section 4 below. This can be avoided by using (3.66) in Eq. (3.72), resulting in The renormalization scheme based on (3.76) and (3.75) for β and (3.41) within the BFM for α in the THDM is denoted as BFMS in the following. 16 We note that the renormalization of α in the "on-shell tadpole-pinched scheme" of Ref. [26] is equivalent to the one based on (3.41) within the BFM.
There are a number of further possibilities to fix the counterterm δβ using rigid invariance in the BFM. The discussion of Sections 3.3.1 and 3.3.2 holds also for the mixing between the pseudoscalar or charged scalar fields. This gives rise to the relations where we can use the field renormalization constants of either  can be fixed upon choosing a specific gauge, such as the 't Hooft-Feynman gauge (within the BFM) and a specific equation.
For the gauge dependence of the renormalization schemes based on the BFM the same remarks as in Section 3.3.1 apply.

Summary of renormalization schemes
The various renormalization schemes used for the HSESM in this paper are summarized in Table 1. The MS and FJ schemes discussed in Ref. [37] are identical to the MS(PRTS) and MS(FJTS) schemes of this paper, up to the point that λ [37] 12 = λ 3 /2 is used as MSrenormalized parameter instead of λ 1 . The MS scheme of Ref. [28] is identical to the MS(FJTS) scheme of this paper.
The renormalization schemes used for the THDM in this paper are summarized in Table 2. The MS(λ 3 ) and FJ(λ 3 ) schemes of Ref. [29], where the scalar coupling λ 3 replaces the angle α as MS-renormalized parameter of the MS(α) and FJ(α) schemes, are not considered in this paper. The MS scheme of Ref. [27] coincides with the MS(FJTS) scheme of this paper up to the fact that λ 5 is used as MS-renormalized parameter instead of M 2 sb = M 2 A 0 + 4M 2 W s 2 w λ 5 /e 2 .

Matching procedure and running couplings
For a comparison of predictions based on different renormalization schemes a conversion of renormalized parameters is necessary. The matching between different schemes is based on the fact that the bare parameters defining the model are renormalization-scheme independent. Following Ref. [29], we describe two variants that can be used to convert renormalized parameters at NLO. Denoting a set of input parameters generically as {p i } defined in two different renormalization schemes "(1)" and "(2)", the two different renormalization schemes are connected via j } by {p (1) j } in the last term, so that which is valid up to terms beyond NLO. The results of the two versions agree in NLO accuracy. The advantage of the full conversion is mainly the exact invertibility, i.e. converting from scheme (1) into (2) and back into (1), reproduces the parameters {p (1) j } exactly, while the linearized version reproduces those parameters only in NLO accuracy.
In the subsequent section, we make use of the full conversion to translate benchmark scenarios between the various renormalization schemes. For the THDM, this means that up to three parameters are converted at a time; for the HSESM, only up to two parameters are concerned. For each benchmark scenario we perform the conversion at a chosen central scale µ = µ 0 , where µ 0 is chosen according to the process and model. In the extension of Prophecy4f [57,58] to the HSESM [37] and THDM [29,41], which is used to calculate Higgs decay widths into four fermions at NLO, (3.85) is solved numerically in double precision by minimizing the χ 2 of the error when solving for {p (2) j } by iteration. 17 In Recola2, quasi-Newton methods (based on the L-BFGS algorithm) are used to find the full solution compatible in double precision. This procedure is used to prepare run cards for the Monte Carlo program Hawk 2.0 [59], with the parameters already converted from a specific input scheme and evolved to the scale, which are then used to evaluate Higgs-production cross sections at NLO.
In order discuss scale uncertainties, we perform the scale variation by including effects of running MS-renormalized parameters. The running is given by the following coupled system of differential equations where β p i denotes the β-function of p i which is non-zero if p i is defined in an MS scheme. For the integration we use standard Runge-Kutta techniques.

Higher-order ambiguities in the conversion
The "full" parameter conversion described in the previous section suffers from ambiguities that are connected to higher-order contributions beyond NLO accuracy. Firstly, the whole renormalization programs addressed in this paper are worked out to NLO, i.e. many relations between renormalizations constants are worked out to linear order only. Beyond NLO, many missing terms should be completed. Secondly, taking the matching equation (3.85) at NLO (or any fixed order) in the "full" conversion, it should be noted that not all UV divergences cancel exactly, because the parameters in the coefficients in front of the UVdivergent contributions in δp (i) are not the same (but defined in different schemes). That means that some UV scale has to be fixed in the evaluation of δp (i) (which is typically set to the renormalization scale), the effect of which, however, is beyond NLO. We, thus, cannot claim that the full conversion is more precise than the linearized version; both are equivalent at NLO. As already mentioned, the full conversion has the advantage of being exactly invertible. Another benefit of supporting both versions is the fact that the comparison of the respective results gives the typical size of effects beyond NLO, which often helps to identify scenarios which are perturbatively unstable.
Looking beyond NLO, it should be mentioned that again specific care has to be taken with respect to the inclusion of tadpole contributions in the matching equation (3.85). For instance, in the PRTS scheme tadpole contributions enter the relations between bare parameters, while in the FJTS scheme they do not. Thus, if the matching of the schemes is done in the original parametrization of the theory (i.e. at the level of µ 2 and λ parameters in the scalar potential), in general there will be tadpole contributions in δp where ∆T i δt (PRTS) j is some function of PRTS tadpole counterterms δt PRTS j . Besides the treatment of tadpoles other sources of ambiguities exist. Since the conversion is performed at fixed loop order the precise choice of the independent parameters matters in higher orders. For instance, the matching of bare parameters for the mixing angles can be performed directly on the mixing angles α B = α + δα, β B = β + δβ, (3.90) or it can be performed on derived parameters such as tan α B = tan α + δ tan α, tan β B = tan β + δ tan β, (3.91) which will result in a conversion equivalent only at NLO. In App. F the conversion tables for the input parameters of Section 4 are given. The conversion is performed in two variants which exemplify the before-mentioned ambiguities and their impact on final results.

Higgs-boson decays H 1,2 → WW/ZZ → 4 fermions
The Monte Carlo program Prophecy4f [57, 58] provides a "PROPer description of the Higgs dECaY into 4 Fermions" and calculates observables for the decay process Higgs → WW/ZZ → 4 fermions at NLO EW+QCD. Its first version [57,58] was designed for the SM and slightly generalized to include a possible fourth fermion generation in Ref. [60]. In Refs. [29,41] and [37], Prophecy4f was extended to the corresponding decays of the light CP-even Higgs bosons of the THDM and the HSESM, respectively, keeping the functionality and applicability of the program basically the same. The THDM and the HSESM were renormalized using MS schemes for the mixing angles α and β. In the following, we present first results from a further extension of Prophecy4f 18 which covers the decays of the heavy CP-even scalar bosons as well and which additionally supports the on-shell and symmetry-inspired renormalization schemes described in this paper. All Higgs-boson decays via intermediate on-or off-shell EW gauge bosons W/Z into all light fermions (all other than top quarks) are supported, but potential decays of a heavy Higgs boson into a pair of light Higgs bosons such as H 1 → H 2 H 2 , or to tt pairs are considered as separate processes and not included in the calculation.

Higgs-Singlet Extension of the Standard Model
In Ref. [37], the HSESM was renormalized in two schemes, called MS and FJ there, which are identical with the MS(PRTS) and MS(FJTS) schemes of this paper up to the point that different MS-renormalized quartic scalar couplings λ i are used to parametrize the Higgs sector. In Ref. [37], the coupling λ 12 = λ 3 /2, which mixes the doublet and singlet scalars was chosen, while we choose the quartic coupling λ 1 of the singlet sector in this paper instead. For the application to the decays H 1,2 → 4f , this choice makes only a marginal difference, which is even beyond NLO, since the renormalization of those quartic Higgs couplings does not enter in this case. Only a minor effect from the different running of the couplings in some one-loop corrections remains. In Ref. [37], five different HSESM scenarios were considered, which are still compatible with current LHC results. These scenarios, which are called BHM200 ± , BHM400, BHM600, and BHM800, identify the light Higgs boson H 2 with the discovered state with a mass of ∼ 125 GeV and contain a heavy Higgs boson H 1 of mass M H 1 = 200 GeV, . . . , 800 GeV, as   suggested by the names of the scenarios. The HSESM parameters of scenarios BHM200 ± , BHM400, and BHM600 are summarized in Table 3. In the following, we take over these scenarios and refer to Ref. [37] for the precise values of the SM-like input parameters. Table 4 shows a comparison of the LO and NLO results for the decay width of the light Higgs boson, H 2 → 4f , obtained in the different renormalization schemes for each scenario. Figure 1 illustrates the scale dependence of Γ H 2 →4f on the l.h.s. for scenario BHM200 + ; the results for the other scenarios look qualitatively similar. The input parameters of the scenarios are defined in the OS scheme and consistently converted into the other schemes, as described in   respectively. One crucial observation in Table 4 is that all renormalization schemes deliver central NLO results differing only below the permille level, while the LO results deviate significantly by up to 2%, i.e. the renormalization-scheme dependence reduces drastically in the transition from LO to NLO. Another important observation concerns the dependence on the renormalization scale. For the MS(PRTS) and MS(FJTS) schemes, in which the MS definitions of the angle α lead to a running α(µ), the residual scale dependence of the decay widths is reduced from 5% at LO to 0.7% at NLO. This observation was already made in Ref. [37], since the MS and FJ schemes from there almost coincide with the MS(PRTS) and MS(FJTS) schemes used here. The residual NLO scale uncertainty of 0.7%, as estimated from rescaling µ 0 by factors of 1/2 and 2, is, thus, of the order of magnitude typically expected from a well-behaved EW NLO calculation. Clearly, the higher-order scale dependence of the OS and BFMS schemes cannot be used as an estimate for the theoretical uncertainty, while the difference of these schemes reflects part of the renormalization-scheme dependence. Table 5 shows a comparison of the LO and NLO results for the decay width of the heavy Higgs boson, H 1 → 4f , obtained in the different renormalization schemes for each scenario. Figure 1 illustrates the scale dependence of Γ H 1 →4f on the r.h.s. for scenario BHM200 + , while again the results for the other scenarios look qualitatively similar. Note that the phenomenology of the decays H 2 → 4f and H 1 → 4f is rather different. Firstly, the chosen values for the mass M H 1 of the heavy Higgs boson are larger than 2M W and 2M Z , so that the decays can proceed via two resonant W or Z bosons, while the decays H 2 → 4f involve at least one off-shell W/Z boson. This effect and the larger phase space for H 1 leads to a large enhancement of the decay width for H 1 → 4f . This enhancement is somewhat damped by the second difference of the two decay types. While the LO decay width of H 2 → 4f involves an explicit factor of c 2 α w.r.t. the SM case, the LO decay width of H 1 → 4f contains the complementary factor s 2 α . Since s α ∼ 0.2−0.3 in the considered scenarios, Γ H 1 →4f is reduced by a factor of ∼ 0.04−0.1 w.r.t. to the SM decay width with the same (hypothetical) Higgs-boson mass M H 1 . The renormalization-scale dependence in the MS(PRTS) and MS(FJTS) schemes is reduced by a factor 5−10 at NLO as compared to LO. A similar reduction is observed for the differences between the renormalization schemes.

Two-Higgs-Doublet Model
In Refs. [29,41], the THDM was renormalized in four schemes, called MS(α), FJ(α), MS(λ 3 ), and FJ(λ 3 ) there, where the first two are identical with the MS(PRTS) and MS(FJTS) schemes of this paper, respectively. In the other two schemes, MS(λ 3 ) and FJ(λ 3 ), the angle α is replaced by the scalar self-coupling λ 3 as MS-renormalized input parameter; these two schemes are not considered in the following.
In Refs. [29,41], following suggestions in the literature, various different THDM scenarios were considered, which are still compatible with current LHC results and identify the light CP-even Higgs boson H 2 with the discovered state with a mass of 125 GeV. In the following, we pick four of those scenarios, covering typical cases with light or heavy Higgs bosons in addition to the known of mass 125 GeV. Table 6 summarizes the corresponding THDM input parameters, while the remaining SM-like parameters are taken over from Refs. [29,41].    Table 7 shows a comparison of the LO and NLO results for the H 2 → 4f decay width obtained in the different renormalization schemes for each scenario. Figure 2 illustrates the scale dependence of Γ H 2 →4f in greater detail.
Among the two schemes with MS-renormalized mixing angles α and β, the MS(PRTS) scheme delivers NLO results with a sufficiently well reduced renormalization-scale dependence of 1−2% in scenarios A1, A2, and B1, while this is the case for the MS(FJTS) scheme only in scenarios A1 and B1. For scenario B2, neither the MS(PRTS), nor the MS(FJTS) scheme produces results with a visible stability region in µ. A possible reason is the relatively small difference M H 1 − M H 2 in the scalar masses in scenario B2 which can lead to enhanced corrections owing to the mixing renormalization constants Eq. (2.5)    parameter conversion (c.f. Table 15 in App. F) into the schemes, which involves very large corrections and shows significant differences between linearized and full conversions. It should be mentioned, however, that the versions of these MS-like renormalization schemes in which α is replaced as input parameter by the coupling parameter λ 3 (the schemes MS(λ 3 ) and FJ(λ 3 ) of Refs. [29,41]) behave somewhat better and at least produce some narrow plateau regions in the vicinity of the scale µ 0 (not shown in this paper).
The OS schemes and the BFMS variant, on the other hand, deliver NLO results which agree within very few per mille, which is also the size of the residual scale dependence resulting from the running of λ 5 in the loop corrections. An exception is scenario B2, which shows a residual scale dependence of up to 1%. Now we turn to the discussion of the decays of the heavy CP-even Higgs boson H 1 . Table 8 and Fig. 3     unexpectedly, the renormalization schemes based on MS-renormalized mixing angles deliver problematic results. At least for the scenarios A1, A2, and B1, the NLO results of the MS(PRTS) scheme still show extrema in the µ dependence that may be interpreted as plateaus of stability, but those regions are rather small and do not cover a range in µ from µ 0 /2 to 2µ 0 . Though the extrema are located near µ 0 , and the NLO results for Γ H 1 →4f in this region are compatible with the results obtained in the other stable schemes. The MS(FJTS) scheme fails to produce stability regions in general, as a result of the very strong running of β − α. For scenario B2, both the MS(PRTS) and the MS(FJTS) scheme fail badly, as expected already from the bad behaviour observed for the decay of the light Higgs boson H 2 above. The OS and BFMS schemes, however, mostly deliver NLO results that are in nice mutual agreement, with the only exception of the OS1 scheme in scenario B1, where the running of λ 5 in the loop corrections introduces a scale uncertainty of the order of 20−30%. Note that also the NLO correction in this scheme is extremely large.

Higgs-boson production processes at the LHC
As a second application of our renormalization schemes, we consider the EW corrections to the production of BSM Higgs bosons in association with a vector boson, also known as Higgs-strahlung, and in association with two jets, usually referred to as vector-boson fusion (VBF). For the numerical integration we have used a modified version of the Monte Carlo program Hawk 2.0 [59]. Hawk 2.0 is a Monte Carlo integrator for Higgs-strahlung and VBF in the SM, including the full fixed-order NLO QCD and EW corrections [61][62][63]. The modification of Hawk 2.0 concerns an interface with the Recola2 library which has been introduced in Ref. [27]. Recola2, the successor of Recola [64], is a highly efficient one-loop amplitude provider for the SM and BSM theories which is publicly available [65]. All necessary ingredients for the integration of the before-mentioned processes for BSM theories can be automatically generated by Recola2 which, in turn, upgrades the original Hawk 2.0 to a general integrator for BSM Higgs production in Higgs-strahlung and VBF as long as no external charged Higgs boson is considered. A new release of Hawk 2.0 will be made available shortly.

Cut setup and parameters
For the analysis of Higgs-strahlung we focus on the final state with one charged muon and a muon-neutrino, pp → Hµ + ν µ + X at 13 TeV. The muon is not recombined with collinear photons, and is assumed to be perfectly isolated, treated as bare muon as described in Ref. [63]. We use similar cuts to the ones given in Ref. [66], i.e. we demand the muon • has transverse momentum p T,µ + > 20 GeV, • be central with rapidity y µ + < 2.4, and require a missing transverse momentum of p miss T > 25 GeV. For the Higgs-boson production in VBF we require two hards jets emitted from partons i which fulfil • pseudo-rapidity |η i | < 5.   The jet definition is performed using the anti-k T algorithm [67] with jet size D = 0.4. The jets j i , i = 1, 2, are required to fulfil typical VBF cuts (see e.g. Ref. [68]): • transverse momentum p T,j i > 20 GeV, • rapidity |y j i | < 5, • rapidity difference |y j 1 − y j 2 | > 3, • opposite hemispheres y j 1 y j 2 < 0, For Higgs-strahlung we investigate all the benchmark scenarios given in Table 3 and Table 6 in the HSESM and THDM, respectively. For VBF we only consider the points A1 and A2 in Table 6 in the THDM, since in general the results look pretty similar to those for Higgsstrahlung.

Higgs-strahlung in the HSESM
In Tables 9 and 10 we show the results for pp → Hµ + ν µ +X for H = H 2 and H 1 , respectively. For the light-Higgs-boson production in Table 9 the LO scale dependence is small, ranging from ∼ 1.5% for BHM400 and BHM600 to ∼ 3% for BHM200 ± . The scale uncertainty gets significantly reduced at NLO by a factor of 4-5. The NLO scale uncertainty for all on-shell schemes and the MS(PRTS) are below 0.05%. Overall, the observed EW corrections are small, and central results are consistent within all schemes, with a scheme dependence at the permille level. A detailed scale dependence of BHM200 + and BHM200 − is shown in Fig. 4.
For heavy-Higgs-boson production H 1 in Table 10 the picture is qualitatively the same with the difference being that the size of the corrections and the scale uncertainties are amplified with respect to light-Higgs-boson production. However, phenomenologically, the scenario corresponds to an entirely different situation, since the heavy Higgs boson is only coupling weakly to the vector bosons and loop contributions significantly contribute to the production, resulting in large EW corrections. In Fig. 5 we depict the scale dependence for BHM400 and BHM600. While for BHM400 the MS schemes show a significant reduction in the scale uncertainty, for BHM600 the absolute scale variation does not change from LO to NLO, though a reduction of the relative scale uncertainty results from the increase of the integrated cross section. It is notable that even though the corrections reach up to 100% all schemes agree at the central scale within less than about 2%. We conclude that also for the heavy-Higgs-boson production all schemes are mutually consistent with no sign of artificially enhanced corrections.

Higgs-strahlung in the THDM
The results for Higgs production in Higgs-strahlung in the THDM are shown in Tables 11  and 12 for the light (H 2 ) and heavy (H 1 ) Higgs-boson production, respectively, with the OS12 as input scheme. Figures 6 and 7 show the scale dependence for pp → H 1,2 µ + ν µ + X  Fig. 4, but for H 1 production in the benchmark scenarios BHM400 (left) and BHM600 (right).   Table 9, but for the cross section σ [fb] of heavy Higgs-boson production in Higgs-strahlung, pp → H 1 µ + ν µ + X, in the HSESM.
for different renormalization schemes.
As illustrated in Fig. 6, the on-shell schemes yield stable results for light-Higgs-boson production. Table 11 shows that neglecting B2 for the moment and comparing to Higgsstrahlung in the HSESM, the reduction of the scale uncertainty is either less strong (e.g. for A1, A2, B1 in MS(PRTS), and A1 in MS(FJTS)) or not observed (A2, B1 in MS(FJTS)). Nonetheless, the results for the benchmark scenarios A1, A2, and B1 visibly agree in all   As can be seen in Table 12, for the heavy-Higgs-boson production in general the scale uncertainty within the traditional [µ 0 /2, 2µ 0 ] window remains large within the MS schemes. Nonetheless, extrema or at least regions of smaller scale dependence show up for the benchmark scenarios A1 and A2 (c.f. Fig. 7) which can be viewed as narrow plateaus.
Focusing on scenario A1 the scale uncertainty is still large at NLO in the MS schemes, while the results in the on-shell schemes are well consistent and their spread decreases at NLO. Yet uncertainties of 4% are visible for OS1 which are not unexpected due to the genuinely large corrections for heavy-Higgs-boson production in all on-shell schemes. The large corrections are due to the proximity of the scenario A1 to the alignment limit and also +4.9% Table 12: As in Table 11, but for the cross section σ [fb] of heavy Higgs-boson production in Higgs-strahlung, pp → H 1 µ + ν µ + X, in the THDM.    cause large conversion effects between the on-shell or BFM and MS schemes (see Table 15 in App. F). For scenario A2, which is further away from the alignment limit, the picture is slightly better.
For scenario B1 the conversion effects between the schemes MS(PRTS), OS2, OS12, and BFMS are small (see Table 15). This is reflected in a LO prediction of similar size and a surprisingly good agreement at NLO between those schemes, even though the corrections exceed 100% in heavy-Higgs-boson production. The MS(FJTS) scheme fails to give any reasonable result in heavy-Higgs-boson production due to the large conversion effects pushing into the alignment limit c αβ ≈ 0. The large conversion effects observed for the OS1 scheme can be traced back to the renormalization of β, and replacing (3.21) with e.g. (3.24), but keeping the renormalization of α fixed, results in small conversion effects. While, we could not completely disentangle the reason for the large effects in (3.21), this is a least partly caused by the enhancement proportional to t β in (3.21) which is not present in the OS2 (3.24), OS12 (3.30), and BFMS schemes. While the OS2 scheme might still be affected by similar problems in certain parameter regions, the scheme OS12 is free of such artificial enhancements and therefore preferable. The seemingly fair agreement between the OS1 and the other on-shell and BFM schemes at NLO is just accidental, since for the Higgs decay (Table 8) no good agreement is found. We conclude that the schemes OS2, OS12, BFMS, and MS(PRTS) give consistent predictions even though the K factor is about 2.
For scenario B2, conversion effects between MS and on-shell or BFM schemes (see Table 15) are sizeable, while within the on-shell and BFM schemes these are small. For the on-shell and BFM schemes, even for heavy-Higgs-boson production the corrections are well behaved, ranging between −11% and +1%. The MS schemes suffer from very large scale uncertainties and large corrections at the central scale.

Higgs production via vector-boson fusion in the THDM
For Higgs production via VBF we provide only some exemplary results. We show the scale dependence for light-Higgs-boson production in A1 and A2 in Fig. 8, and summarize the usual scale variation in Table 13. Since the behaviour of these results resembles closely the one for Higgs-strahlung, apart from the magnitude of the cross sections, we did not investigate any other benchmark scenarios. We plan to provide more detailed phenomenological results on VBF elsewhere.

Conclusions
Models with extended Higgs sectors are of prime importance for investigating the mechanism of electroweak symmetry breaking. Precision investigations of these models require the inclusion of higher-order corrections at least at NLO and thus renormalization. In particular, prescriptions for the renormalization of mixing angles in the scalar sector are needed.
In this paper we have discussed a variety of renormalization prescriptions for scalar mixing angles with particular regard to symmetry, gauge independence, and numerical stability. In detail, we have considered and compared three types of renormalization schemes for mixing angles in the Higgs sector: • MS renormalization conditions for mixing angles are easy to implement. They depend, however, on the treatment of tadpoles and require care in view of gauge dependence, but have the benefit that the size of missing higher-order corrections can be investigated by renormalization scale variation.
• We have formulated on-shell renormalization conditions for the mixing angles in the Two-Higgs-Doublet Model and the Higgs-Singlet Extension of the Standard Model based on combinations of physical observables. More precisely, ratios of matrix elements or formfactors depending on the desired mixing angles only, deliver appropriate renormalization conditions. To obtain such ratios, it is often useful to introduce spurious particles with infinitesimal couplings, which do not change the physical theory.
• Rigid gauge invariance and/or the background-field method allow to introduce renormalization conditions for mixing angles in general theories.
We numerically studied and compared various renormalization conditions for mixing angles in the Two-Higgs-Doublet Model and the Higgs-Singlet Extension of the Standard Model for Higgs decays into four fermions and for Higgs-production in association with a vector boson or in vector-boson fusion. While renormalization schemes based on MS subtraction tend to become unstable in delicate scenarios, in particular for heavy Higgsboson production in the Two-Higgs-Doublet Model, the proposed on-shell schemes and the schemes motivated by the background field method behave decently. They do not allow for a realistic estimate of missing higher-order corrections via scale variation, but instead the renormalization-scheme dependence can be investigated by comparing results obtained with different schemes after a consistent conversion of input parameters between the schemes. In general, one should avoid renormalization schemes that introduce artificial enhancement factors, e.g. for degenerate masses or small mixing angles.
Based on our study, we propose to use on-shell or symmetry-based schemes for the central predictions, as these turn out to be more robust. The reliability of the results can be checked by using two or more different on-shell schemes, where already the consistent parameter conversion between different schemes provides uncertainty estimates. Finally, schemes based on MS renormalization can be used, with due care, to study scale uncertainties.
The schemes introduced in this paper are based on genuine on-shell conditions or simple symmetry principles. Thus, their generalization to higher orders is well-defined (similar to the MS scheme), even if the explicit results for the counterterms might become more involved owing to reducible contributions. This is in contrast to some of the schemes proposed in the literature that rely on the form of the one-loop expressions.
The proposed schemes are not restricted to the considered models, but can be generalized to other extensions of the Standard Model, involving additional scalars, vector bosons, or fermions.

B Translation of self-energies to the background-field method
In Section 3.3 we have derived renormalization conditions based on the background-field method which requires the evaluation of mixing and self-energies in this framework. Since the expressions differ from those in the conventional formalism, we provide the differences relevant for the renormalization of the mixing angles in the HSESM and the THDM using the following notation Some corresponding results in the SM can be found in Ref. [69]. The one-loop tadpoles (one-point functions) do not differ between the BFM and the conventional formalism, although the individual tadpole contributions are not the same diagram by diagram. As a consequence also all tadpole contributions to the self-energies are the same in the BFM and the conventional formalism, and the differences reported in the following only result from 1PI contributions.
In order to use the mixing-angle renormalization based on the BFM, only the parameter counterterms have to be calculated in the BFM while the loop diagrams and the wave-function renormalization constants can be calculated in the conventional formalism (owing to gauge independence of the sum of a bare amplitude and its corresponding wavefunction counterterms). As far as parameter counterterms are determined from on-shell field renormalization constants, the latter have to be calculated in the BFM as well. To calculate the counterterms in the BFM formalism it is sufficient to combine the results for the differences (B.1) given in this appendix with the conventional self-energies.
In order to distinguish the electromagnetic coupling α em from the mixing angle α, we label it with the index "em".

B.1 Higgs sector of the HSESM
In the HSESM the differences for the self-energies relevant for the renormalization of the mixing angle α read for H i,j ∈ {H 1 , H 2 } For the pseudo-scalar and charged sector we have

B.2 Higgs sector of THDM
In the THDM, the differences read for H i,j ∈ {H 1 , H 2 } For the pseudo-scalar fields we find and for the charged ones we have C Tadpole contributions to scalar mixing and self-energies In general, the vertex functions contain explicit and implicit tadpole contributions. The explicit tadpole contributions, i.e. tadpole loop diagrams, result from the expansion of the effective action about a point that does not correspond to the stationary point. Implicit tadpole contributions or tadpole counterterms originate from tadpole terms in the Lagrangian or shifts in the fields. Since the treatment of tadpoles involves some freedom, and as we employ two MS schemes differing in the treatment of tadpole counterterms, we briefly summarize their properties and provide, for completeness, the tadpole-counterterm expressions for mixing and self-energies that are needed in those renormalization schemes. In the PRTS, the expansion is performed about the stationary point at the one-loop level, and the bare (squared) masses of the physical fields are defined as the coefficients of the terms quadratic in the corresponding fields. Moreover, we require vanishing mixing between the physical (mass-eigenstate) scalar fields, i.e. there are no implicit tadpoles in the mixing energies of these "physical fields". In the THDM we define the mixing angle β from the ratio of the true vevs, tan β = v 2 /v 1 . In this scheme, the tadpole terms result exclusively from tadpole counterterms in the Lagrangian and are absent, by definition, in two-point functions that involve only physical fields.
In the FJTS, all bare parameters are defined in terms of the bare parameters of the symmetric Lagrangian. The bare scalar fields are shifted, H B,i → H B,i + ∆v i , so that the shifted fields describe excitations about the stationary point. The tadpole terms result exclusively from the shifts ∆v i of the scalar fields. Tadpole terms appear in almost all vertex functions, apart from the one-point functions.
In the following we provide the implicit tadpole counterterms for the scalar self-energies and mixing energies in the THDM and HSESM both in the PRTS and the FJTS, using the following convention for the tadpole counterterms, corresponding to the second terms on the r.h.s. in (2.6). While in the PRTS the sum of the third and fourth terms in (2.6) is zero owing δtĤ = −TĤ , in the FJTS this condition must not necessary be fulfilled but is convenient.

E Background-field Ward identities
In the BFM, the invariance of the effective action under background gauge transformations gives rise to simple Ward identities (WIs) for the vertex functions of the background fields [34]. These WIs depend on the treatment of tadpoles. 19 In general, the generating functional of vertex functions Γ is related to the generating functional T c of connected Green functions via the Legendre transformation This relation implies, in particular, that the presence of tadpoles in the propagators is directly connected to the presence of tadpoles in the 2-point vertex functions. Accordingly, we define the self-energies Σ as the higher-order contributions to the 2-point vertex functions (full inverse propagators) including all relevant tadpole contributions (c.f. Eq. (2.6)). The 1PI contribution to the self-energies Σ 1PI includes only the first term on the r.h.s. of (2.6).

E.1 A Standard Model example
We illustrate the influence of the tadpole treatment on the BFM WIs using an example for the SM. Invariance of the effective action Γ under background gauge transformations related to the parameterθ Z yields in general where {X} represents the set of all background fields of the SM and we suppressed some terms that are irrelevant in the following. Here,v +Ĥ(x) denotes the decomposition of the Higgs-boson field into a constant vevv and the field excitationĤ(x), which will be made more precise for the different renormalization schemes.
The WIs are obtained by taking derivatives of (E.4) with respect to some background fields and evaluating at specific field valuesX(x) =X. Differentiating for instance with respect to the would-be Goldstone-boson fieldχ yields X − e 2c w s w v +Ĥ δ 2 Γ δχ(x)δχ(y) X + e 2c w s w δ(x − y) δΓ δĤ(x) X + e 2c w s wχ i.e. the connected one-point functions, the explicit tadpoles, vanish and consequently the vertex functions are one-particle irreducible (1PI) as indicated in (E.7). Moreover, no implicit tadpoles occur since we expand about the bare vacuum. Note, however, that the so-defined vertex functions do not correspond to a generating functional of connected Green functions for vanishing sources and cannot be used to calculate the S-matrix in a straight-forward way unless additional tadpole contributions are included properly.
In the case of spontaneous symmetry breaking, the correct generating functional of connected Green functions is related to the vertex functional at the stationary pointX s , defined by Splitting the self-energies into irreducible parts and implicit tadpole terms results in In the last step we used the fact that the Higgs tadpole is cancelled by the corresponding counterterm (implicit tadpole), i.e. δtĤ = −TĤ . The last line of Eq. (E.12) coincides with (E.7) and the WI (30) in Ref. [34]. In the FJTS, two different approaches can be used. In the first approach, which is our default, the fields are expanded about the stationary point, but the full vev consists of its tree-level part v B and the field shift ∆v, Splitting the self-energies into irreducible parts and implicit tadpole terms results in = p 2 ΣẐχ 1PI (p 2 ) − iM Z Σχχ 1PI (p 2 ) + i e 2c w s w TĤ . (E.17)

E.2 A THDM example
We give a second example 21 that is relevant for the renormalization of the mixing angle β in the THDM. Starting from the analogue of the WI (E.4) in the THDM and differentiating with respect to the physical pseudoscalar field A 0 yields when expanding about the stationary point: δÂ 0 (x)δÂ 0 (y) X s + . . . .

F Parameter conversion tables
In this appendix we give results for the parameter conversion from the on-shell input schemes OS and OS12 in the HSESM and THDM, respectively, to the other renormalization schemes in two different variants, as used for the results presented in Section 4. The two variants are full conversions, as explained in Section 3.5, are equivalent at NLO, but differ in higher orders due to the tadpole counterterm scheme and the precise form of the input parameters. For the details see the caption of the tables in the following.
F.1 HSESM scenarios of Table 3 In Table 14 the results for the full conversion of parameters from the OS as input scheme to other schemes for the considered benchmark scenarios of Table 3 are shown. No large conversion effects are observed, and the two different conversion variants mutually agree.
F.2 THDM scenarios of Table 6 In Table 15 the corresponding conversion in the THDM from the OS12 as input scheme to other renormalization schemes is shown. As compared to the HSESM, the conversion effects for the scenarios in Table 6 are more pronounced. Generically the conversion effects are 5% for the on-shell schemes and for the BFMS scheme, with the only exception being the OS1 scheme in scenario B1 with conversion effects of the order of ∼ 40%. The conversion effects to the MS schemes are typically larger than those to the on-shell schemes, and the conversion to MS(FJTS) becomes perturbatively unstable (and thus fails completely) for scenarios B1 and B2 while for MS(PRTS) only B2 is unstable.
Comparing the two conversions, the differences in the scenarios A1 and A2 to MS(PRTS) (MS(FJTS)) amount to 13%(6%) and 1%(1%) for c αβ and t β , respectively. The different conversions to the on-shell schemes and the BFMS scheme agree on the level of 4%, with scenario B1 being again the only exception where we find a difference of 33% in the conversion of t β .