Renormalization schemes for the Two-Higgs-Doublet Model and applications to h ->WW/ZZ ->4fermions

We perform the renormalization of different types of Two-Higgs-Doublet Models for the calculation of observables at next-to-leading order. In detail, we suggest four different renormalization schemes based on on-shell renormalization conditions as far as possible and on MSbar prescriptions for the remaining field-mixing parameters where no distinguished on-shell condition exists and make contact to existing schemes in the literature. In particular, we treat the tadpole diagrams in different ways and discuss issues of gauge independence and perturbative stability in the considered schemes. The renormalization group equations for the MSbar parameters are solved in each scheme, so that a consistent renormalization scale variation can be performed. We have implemented all Feynman rules including counterterms and the renormalization conditions into a FeynArts model file, so that amplitudes and squared matrix elements can be generated automatically. As an application we compute the decay of the light, CP-even Higgs boson of the Two-Higgs-Doublet Model into four fermions at next-to-leading order. The comparison of different schemes and the investigation of the renormalization scale dependence allows us to test the perturbative consistency in each of the renormalization schemes, and to get a better estimate of the theoretical uncertainty that arises due to the truncation of the perturbation series.


Introduction
After the discovery of a Higgs boson at the Large Hadron Collider (LHC) [1,2] at CERN, the complete identification of this particle is ongoing. The properties of the discovered particle, such as its couplings, are determined experimentally in order to fully identify its nature. For the endeavour of the identification of this particle, input from the theory side is needed in form of precise predictions for the production and decay processes in the Standard Model (SM) as well as in its extensions that are to be tested. It is also crucial to provide reliable uncertainty estimates of the theoretical predictions. Underestimating this uncertainty might lead to wrong conclusions. In the SM, predictions and error estimates are well advanced, and in SM extensions they are consolidating as well (see, e.g., the reviews in Refs. [3][4][5][6][7][8]).
One of the simplest extensions of the SM is the Two-Higgs-Doublet Model (THDM) [9,10] where a second Higgs doublet is added to the SM field content. The underlying gauge group SU(3) C × SU(2) W × U(1) Y as well as the fermion content of the SM are kept. After spontaneous symmetry breaking, there are five physical Higgs bosons where three of them are neutral and two are charged. In the CP-conserving case, which we consider, one of the neutral Higgs bosons is CP-odd and two are CP-even with one of them being SM-like.
Even such a simple extension of the SM can help solving some questions that are unanswered in the SM. For example, CP-violation in the Higgs sector could provide solutions to the problem of baryogenesis [11][12][13][14], and inert THDMs contain a dark matter candidate [15,16]. An even larger motivation comes from the embedding of the THDM into more complex models, such as axion [17,18] or supersymmetric models [19]. Some of the latter are promising candidates for a fundamental theory, and supersymmetric Higgs sectors contain a THDM (in which the doublets have opposite hypercharges). Even though the THDM is unlikely to be the fundamental theory of nature, it provides a rich phenomenology, which can be used in the search for a non-minimal Higgs sector without being limited by constraints from a more fundamental theory.
In this sense, it is obvious that the THDM should be tested against data, and phenomenological studies have been performed recently, e.g., in Refs. [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. In order to provide precise predictions within this model, not only leading-order (LO), but also next-to-leading-order (NLO) contributions have to be taken into account. For the calculation of NLO contributions, a proper definition of a renormalization scheme is mandatory. There is no unique choice, and applying different renormalization schemes can help to estimate the theoretical uncertainty of the prediction that originates from the truncation of the perturbative series. The renormalization of the THDM has already been tackled in several publications: First, in Ref. [35], the fields and masses were renormalized in the on-shell scheme, however, the prescription given there does not cover all parameters. In Refs. [36,37], a minimal field renormalization was applied, and the field mixing conditions were used to fix some of the mixing angles. In view of an automation of NLO predictions within the THDM a tool was written by Degrande [38] where all finite rational terms and all divergent terms are computed using on-shell conditions or conditions within the "modified minimal subtraction scheme" (MS). Though automation is very helpful, often specific problems occur depending on the model, the process, or the renormalization scheme considered, and, it might be necessary to solve these "manually". Specifically, spontaneously broken gauge theories with extended scalar sectors pose issues with the renormalization of vacuum expectation values and the related "tadpoles", jeopardizing gauge independence and perturbative stability in predictions. Renormalization schemes employing a gauge-independent treatment of the tadpole terms were described recently in Refs. [39][40][41].
In this paper, we perform the renormalization of various types of THDMs (Type I, Type II, "lepton-specific", and "flipped"), describe four different renormalization schemes (for each type), and provide explicit results facilitating their application in NLO calculations. The comparison of results obtained in these renormalization schemes allows for checking their perturbative consistency, i.e. whether the expansion point for the perturbation series is chosen well and no unphysically large corrections are introduced. Knowing in which parts of the parameter space a renormalization scheme leads to a stable perturbative behaviour is important for the applicability of the scheme. In addition, we investigate the dependence on the renormalization scale µ r which is introduced by defining some parameters via MS conditions. In order to investigate the µ r dependence consistently, we solve the renormalization group equations (RGEs) and include the running effects. We also make contact to different renormalization schemes suggested in the literature, including the recent formulations [39,40] with gauge-independent treatments of tadpoles. 1 To facilitate NLO calculations in practice, we have implemented our renormalization schemes for the THDM into a FeynArts [44] model file, so that amplitudes and squared matrix elements can be generated straightforwardly. Finally, we apply the proposed renormalization schemes in the NLO calculation of the partial decay width of the lighter CP-even Higgs boson decaying into four fermions, h → WW/ZZ → 4f , a process class that is a cornerstone in the experimental determination of Higgs-boson couplings, but for which electroweak corrections in the THDM are not yet known in the literature. The impact of NLO corrections on Higgs couplings in the THDM was, for instance, investigated more globally in Refs. [45,46]. However, a full set of electroweak corrections to all Higgs-boson decay processes in the THDM does not yet exist in the literature, so that current predictions (see, e.g., Ref. [47]) for THDM Higgs analyses globally neglect electroweak higher-order effects. Our calculation, thus, contributes to overcome this shortcoming; electroweak corrections to some 1 → 2 particle decays of heavy Higgs bosons were presented in Ref. [39]. 2 The structure of the paper is as follows. We introduce the four considered types of THDMs and our conventions in Sect. 2, and the derivation of the counterterm Lagrangian is performed in Sect. 3. Afterwards we fix the renormalization constants with renormalization conditions (Sect. 4). The on-shell conditions, where the renormalized parameters correspond to measurable quantities, are described and applied in Sect. 4.1. The renormalization constants of parameters that do not directly correspond to physical quantities are fixed in the MS scheme, so that they contain only the UV divergences and no finite terms. We describe different renormalization schemes based on different definitions of the MS-renormalized parameters in Sect. 4.2. The RGEs of the MS-renormalized parameters are derived and numerically solved in Sect. 5. The implementation of the results into an automated matrix element generator is described in Sect. 6, and numerical results for the partial decay width h → 4f are presented in Sect. 7. Finally, we conclude in Sect. 8, and further details on the renormalization prescription as well as some counterterms are given in the appendix.

The Two-Higgs-Doublet Model
The Lagrangian of the THDM, L THDM , is composed of the following parts, L THDM = L Gauge + L Fermion + L Higgs + L Yukawa + L Fix + L Ghost . (2.1) The gauge, fermionic, gauge-fixing, and ghost parts can be obtained in a straightforward way from the SM ones, e.g., given in Ref. [48]. The Higgs Lagrangian and the Yukawa couplings to the fermions are discussed in the following and are mostly affected by the additional degrees of freedom of the THDM. A very elaborate and complete discussion of the THDM Higgs and Yukawa Lagrangians, including general and specific cases, can, e.g., be found in Refs. [49,50]. 1 In our work we do not consider the "tadpole-pinched" scheme suggested in Ref. [39]. Following the arguments of Refs. [42,43] we consider the "pinch technique" just as one of many physically equivalent choices to fix the gauge arbitrariness in off-shell quantities (related to the 't Hooft-Feynman gauge of the quantum fields in the background-field gauge) rather than singling out "its gauge-invariant part" in any sense. 2 Since we consider the decays of the light Higgs boson h via W-or Z-boson pairs, where at least one of the gauge bosons is off its mass shell, we have to consider the full 1 → 4 process with all off-shell and decay effects, rendering a comparison to results on H → WW/ZZ not meaningful.

The Higgs Lagrangian
The Higgs Lagrangian, L Higgs , contains the kinetic terms and a potential V , (2.2) with the complex scalar doublets Φ 1,2 of hypercharge Y W = 1, and the covariant derivative where I a W (a = 1, 2, 3) are the generators of the weak isospin. The SU(2) and U(1) gauge fields are denoted W a µ and B µ with the corresponding gauge couplings g 2 and g 1 , respectively. The sign in the g 2 term is negative in the conventions of Böhm, Hollik and Spiesberger (BHS) [48,51] and positive in the convention of Haber and Kane (HK) [19]. We implemented both sign conventions, but used the former one as default. In general, the potential involves all hermitian functions of the two doublets up to dimension four and can be parameterized in the most general case as follows [50,52], The parameters m 2 11 , m 2 22 , λ 1 , λ 2 , λ 3 , λ 4 are real, while the parameters m 2 12 , λ 5 , λ 6 , λ 7 are complex, yielding a total number of 14 real degrees of freedom for the potential. However, the component fields of the two Higgs doublets Φ 1 and Φ 2 do not correspond to mass eigenstates, and these doublets can be redefined using an SU(2) transformation without changing the physics, so that only 11 physical degrees of freedom remain [50]. For each Higgs doublet we demand that the fields develop a vacuum expectation value (vev) in the neutral component, It is non-trivial that such a stable minimum of the potential exists, restricting the allowed parameter space already strongly [53]. In general, the vevs are complex (with a significant relative phase). The Higgs doublets can be decomposed as follows, with the charged fields φ + 1 , φ + 2 , the neutral CP-even fields η 1 , η 2 , and the neutral CP-odd fields χ 1 , χ 2 .
Additional constraints: Since the 11-dimensional parameter space of the potential is too large for early experimental analyses, we restrict the model in our analysis by imposing two additional conditions, motivated by experimental results: • absence of flavour-changing neutral currents at tree level, • CP conservation in the Higgs sector (even though this holds only approximately).
The former requirement can be ensured by adding a discrete Z 2 symmetry Φ 1 → −Φ 1 (see Sect. 2.2). This condition implies that the parameters λ 6 and λ 7 vanish. Permitting operators of dimension two that violate this Z 2 symmetry softly, non-zero values of m 12 are still allowed [10,54]. Concerning the second condition, the potential is CP-conserving if and only if a basis of the Higgs doublets exists in which all parameters and the vevs are real [55]. For our description we assume that a transformation to such a basis has been done already (if the parameters or vevs were initially complex), so that we only have to deal with real parameters. This renders m 12 and λ 5 real. However, at higher orders in perturbation theory CP-breaking terms and complex phases in the Higgs sector are generated radiatively through loop contributions involving the quark mixing matrix. For our NLO analysis, this does not present a problem as they appear only beyond NLO in the specific processes we consider. In addition we assume that a basis of the doublets is chosen in which v 1 , v 2 > 0 (which is always possible as a redefinition Φ i → −Φ i changes the sign of the vacuum expectation value). The potential (2.5) has then the following form, Expanding the potential using the decomposition (2.7) and ordering terms with respect to powers of the fields, leads to the form with the tadpole terms proportional to the tadpole parameters t η 1 , t η 2 and linear in the fields. The mass terms contain the mass matrices M η , M χ , and M φ and are quadratic in the CPeven, CP-odd, and charged Higgs-boson fields, respectively. Terms cubic or quartic in the fields are suppressed in the notation here. Only the neutral CP-even scalar fields can develop nonvanishing tadpole terms, since they carry the quantum numbers of the vacuum. Further, in the mass terms, only particles with the same quantum numbers can mix, so that the three different types of scalars (neutral CP-even, neutral CP-odd, and charged) do not mix with one another. Of course, through the cubic and quartic terms, which are not shown here, these particles interact with each other. The tadpole parameters are where we introduced the abbreviations λ ij... = λ i + λ j + . . ., and the mass matrices are given by where general abbreviations for the trigonometric functions s x ≡ sin x, c x ≡ cos x, t x ≡ tan x are introduced. After the elimination of m 11 , m 22 using the above equations, the entries of the mass matrices contain also the tadpole parameters t h , t H . Using we obtain for the mass parameters of the CP-even Higgs bosons, the ones of the mass matrix of the CP-odd Higgs fields result in At tree level we demand vanishing tadpole terms corresponding to t H = t h = 0 and diagonal propagators, so that the mixing terms proportional to M 2 Hh , M 2 G 0 A 0 , M 2 GH + vanish at LO. This yields β c = β n = β and fixes the angle α as well. At NLO such a diagonalization is not possible, as the propagators receive also mixing contributions from the field renormalization and from (momentum-dependent) one-loop diagrams, so that there is no distinct condition to define the mixing angles. Therefore we keep the bare mass mixing parameters M 2 Hh , M 2 G 0 A 0 , M 2 GH + and the tadpole terms t H , t h in this section, and specify defining conditions for the bare parameters α, β n , β c and the tadpole terms later. With the above equations, m 12 , λ 1 , λ 2 , λ 3 , λ 4 can be traded for the masses of the physical bosons M H , M h , M A 0 , M H + , and the mixing angle α. The parameter λ 5 cannot be replaced by a mass or a mixing angle as it appears only in cubic and quartic Higgs couplings and acts like an additional coupling constant. Explicit relations can be obtained by inverting Eqs. (2.16), (2.17a), and (2.18a). The other parameters are related to the masses by The tree-level relations are easily obtained by setting β c = β n = β and t Parameters of the gauge sector: Mass terms of the gauge bosons arise through the interaction of the gauge bosons with the vevs, analogous to the SM. After a rotation into fields corresponding to mass eigenstates, one obtains relations similar to the SM ones: where the electric unit charge e is identified with the coupling constant of the photon field A µ in the covariant derivative. Inverting these relations and introducing the weak mixing angle θ W via cos θ W = g 2 / g 2 1 + g 2 2 , one can replace v and the gauge couplings g 1 and , v 1 , v 2 , g 1 , g 2 }, (2.23) in favour of the bare mass parameters including one parameter from scalar self-interactions which we take as λ 3 , Additionally, one has to keep in mind that we keep the mixing parameters α, β n , and β c generic, and they have to be fixed by additional conditions (which will be given later). One can use Eq. (2.19c) to trade λ 3 for α, in which case the mixing angle becomes a free parameter of the theory. Then, one obtains the parameter set

Yukawa couplings
The Higgs mechanism does not only give rise to the gauge-boson mass terms (which are determined by the vevs), but via Yukawa couplings, it introduces masses to chiral fermions. Since both Higgs doublets can couple to fermions, the general Yukawa couplings have the form with the mixing matrices ζ f,k , k = 1, 2, in generation space for the gauge-invariant interactions with Φ 1 and Φ 2 , respectively, and the generation indices i, j = 1, 2, 3. The left-handed SU(2) doublets of quarks and leptons are denoted Q ′L = u ′L , d ′L T and L ′L = ν ′L , l ′L T , while the right-handed up-type quark, down-type quark, and lepton singlets are u ′R , d ′R , and l ′R , respectively. The primes indicate that we deal with fields in the interaction basis here; fields without primes correspond to mass eigenstates. The fieldΦ k , k = 1, 2, is the charge-conjugated field of Φ k . Since, in the general THDM, there are two mass mixing matrices for each type f of fermions, flavour-changing neutral currents (FCNC) can occur at tree level, which, however, are experimentally known to be strongly suppressed. According to the Paschos-Glashow-Weinberg theorem [56,57], FCNC are absent at tree level if each type of fermion couples only to one of the Higgs doublets. This can be achieved by imposing an additional discrete Z 2 symmetry. It should be noted that the soft-Z 2 -breaking term proportional to m 12 in the Higgs potential does not introduce FCNC. The Yukawa Lagrangian reduces then to with n i being either 1 or 2. Depending on the exact form of the symmetry, one distinguishes four types of THDMs. In Type I models, all fermions couple to one Higgs doublet (conventionally Φ 2 , but this is equivalent to Φ 1 due to possible basis changes) which can be ensured by demanding a Φ 1 → −Φ 1 symmetry. In Type II models, down-type fermions couple to the other doublet, which can be enforced by the symmetry The other two possibilities are called "lepton-specific" (Type X) and "flipped" (Type Y) models. An overview over the couplings and symmetries of the different models is given in Tab. 1. For each of the fermion types, a redefinition of the fields can be performed in order to get diagonal mass matrices, analogously to the SM case. Similar to the SM, the coupling of fermions to the Z boson is flavour conserving, and a CKM matrix appears in the coupling to the charged gauge bosons. By specifying the model type, the Higgs-fermion interaction is determined, and one can write them, widely following the notation of Ref. [58], as  . Note that the sign of ξ f A 0 is defined relative to the coupling of the Goldstoneboson field G 0 and the relation β n = β c = β is used here.
where m l,i , m u,i , and m d,i are the lepton, the up-type, and the down-type quark masses, respectively, and V ij are the coefficients of the CKM matrix. Left-and right-handed fermion fields, f L and f R , are obtained from the corresponding Dirac spinor f by applying the chirality projectors The coupling coefficients ξ f H,h,A 0 are defined as the couplings relative to the canonical SM value of m f /v and are shown in Tab. 2. Note that we have used β n = β c = β in Eq. (2.28) and Tab. 2, which is most relevant in applications; the generalization to independent β n , β c , β is simple.

The counterterm Lagrangian
The next step in calculating higher-order corrections is the renormalization of the theory. In this section we focus on electroweak corrections of O(α em ). The QCD renormalization of the THDM is straightforward and completely analogous to the SM case, since all scalar degrees of freedom are colour singlets and do not interact strongly. In the formulation of the basic Lagrangian in the previous section, we dealt with bare parameters and fields. To distinguish those from renormalized quantities, in the following we indicate bare quantities by subscripts "0" consistently. We perform a multiplicative renormalization, i.e. we split bare quantities into renormalized parts and corresponding counterterms, use dimensional regularization, and allow for matrix-valued field renormalization constants in the case that there are several fields with the same quantum numbers. The counterterm Lagrangian δL containing the full dependence on the renormalization constants can be split into several parts analogous to Eq. (2.1), where the Higgs part of the Lagrangian δL Higgs is split up into the kinetic part δL Higgs,kin and the Higgs potential part δV Higgs . Since the gauge fixing is applied after renormalization, no gaugefixing counterterms occur, and since ghost fields occur only in loop diagrams, a renormalization of the ghost sector is not necessary at NLO for the calculation of S-matrix elements. Though, for analyzing Slavnov-Taylor or Ward identities a complete renormalization procedure would be advisable. Our renormalization procedure, thus, widely parallels the treatment described for the SM in [48]; an alternative variant that is based on the transformation of fields in the gauge eigenstate basis, as suggested for the SM in Ref. [51] and for the MSSM in Ref. [59], is described in Ref. [60].

Higgs potential
According to Eq. (2.8), the Higgs potential contains 8 independent parameters which have to be renormalized. In addition, there are two vevs and two gauge couplings completing the set of input parameters of (2.23). We have carried out different renormalization procedures and in this paper discuss the renormalization of the Lagrangian in the mass parameterization (a) with renormalization of the mixing angles, (b) alternatively, taking the mixing angles as dependent parameters and applying the field transformation after renormalization.
Additionally we have performed a renormalization of the basic parameters and a subsequent transformation to renormalization constants and parameters of the mass parameter set similar to Dabelstein in the MSSM [59] 4 . The latter has been used, after changing to our conventions for field and parameter renormalization, to check the counterterm Lagrangian. In the following, we give a detailed description of method (a), while method (b) is briefly described in App. A.

Renormalization of the mixing angles
In this section we show that the counterterms of mixing angles that are not used to replace another free parameter of the theory are redundant in the sense that they can be absorbed by field renormalization constants. We sketch the argument for generic scalar fields ϕ 1 , ϕ 2 , which are transformed into fields h 1 , h 2 corresponding to mass eigenstates by a rotation by the angle θ, where we added subscripts "0" to indicate bare quantities. The general argument can be applied to the neutral CP-even, the neutral CP-odd, and the charged Higgs fields of the THDM by replacing h 1 , h 2 by H, h, or G, A 0 or G ± , H ± , respectively, and by substituting the angle θ by α, β n , or β c . The fields corresponding to mass eigenstates are renormalized using matrix-valued renormalization constants, so that the renormalization transformation reads Applying this renormalization transformation to Eq. (3.2) leads to One can easily remove the dependence on the mixing angle with a redefinition of the mixing field renormalization constants by introducing Then, the Eq. (3.4) reads 4 Details about this method can be found in Ref. [60].
Obviously, the dependence on δθ can always be removed from the Lagrangian by a redefinition of the field renormalization constants. As a simple shift of the mixing field renormalization constant is performing the task, the renormalization of the mixing angle θ can be seen as an additional field renormalization (as it is done, e.g., in Ref. [58]). This argument is general and holds for any renormalization condition on θ. Without loss of generality one can even assume that such a redefinition has already been performed and set δθ = 0 from the beginning, as done in method (b) in App. A. Of course, the bookkeeping of counterterms depends on the way δθ is treated. This can be seen by considering the mass term of the potential. The general bare mass term can be written using the rotation matrix R ϕ (θ 0 ) as This expression can be expanded in terms of renormalized and counterterm contributions, yielding where we obtain off-diagonal terms from the counterterm δθ of the mixing angle and from the renormalization of the mass matrix. The latter contribution depends on the independent counterterms {δp} and is abbreviated by f θ ({δp}). At NLO, the mixing entry of the mass matrix reads and since the renormalization of the redundant mixing angle can be chosen freely, counterterm contributions in the Lagrangian can be shifted arbitrarily from mixing terms to mixing-angle counterterms. Note that f θ ({δp}) does not change by such redistributions, since it is fixed by the remaining renormalization constants.

Renormalization with a diagonal mass matrix -version (a)
In this prescription, we use the angle α as an independent parameter instead of λ 3 . To define α at NLO, we demand that the mass matrix of the CP-even Higgs bosons (in the Lagrangian), written in terms of bare fields, is diagonal at all orders (i.e. in terms of bare or renormalized parameters). Equation (2.19c) with then defines the parameter α. Note that (momentum-dependent) loop diagrams tend to destroy the diagonality of the matrix-valued two-point functions (inverse propagators) in the effective action as well. Below, the field renormalization will be chosen to compensate those loop effects at the mass shells of the propagating particles. It is relation (3.10) that distinguishes α from the case of a redundant mixing angle, such as θ in the previous section, which can be chosen freely or absorbed by field renormalization constants. The relation between α and λ 3 of Eq. (2.19c) is the same for bare and renormalized quantities and can be used to eliminate λ 3 from the theory. For each of the independent parameters of the mass parameter set (2.25) we apply the renormalization transformation so that the 12 parameter renormalization constants are corresponding to {p mass }. In this part we describe the commonly used tadpole renormalization (which is gauge dependent) and describe a gauge-independent scheme in Sect. 4.2.3. The higher-order corrections of the mixing angles β n and β c are irrelevant according to Sect. 3.1.1 and we can choose which defines the mixing terms uniquely and ensures that the angles β n , β c , and β do not have to be distinguished at any order. From these conditions, we can compute the mass mixing terms from Eq. (2.17c) and Eq. (2.18c) to The field renormalization is performed for each field corresponding to mass eigenstates, with the field renormalization constants δZ H , δZ Hh , δZ hH , δZ h , The interaction terms are derived in the same way, but they are very lengthy and not shown here.
The prescription for the field renormalization (3.15) is non-minimal in the sense that a renormalization of the doublets with two renormalization constants, actually would be sufficient to cancel the UV divergences. However, the prescription with matrix-valued renormalization constants allows us to renormalize each field on-shell. The UVdivergent parts of the renormalization constants in Eq. (3.15) cannot be independent, and relations between the UV-divergent parts of the two prescriptions exist. They can be obtained by applying the renormalization prescription (3.17) to the left-hand side and (3.15) to the righthand side of Eqs. (2.12), transforming thereafter the interaction states on the left-hand side to mass eigenstates and comparing both sides. This results in (3.18) We will use these relations to derive UV-divergent parts for specific renormalization constants in Sect. 4.2. In App. A we discuss a different choice of the mixing angles which is suited for the renormalization with λ 3 as an independent parameter.

The Higgs kinetic part
After expressing the Higgs kinetic term in terms of bare physical fields, mixing angles, and parameters, one can apply the renormalization transformations (3.11) and (3.15) to obtain the counterterm part of the kinetic Lagrangian which introduces scalar-vector mixing terms. The explicit terms are stated in App. B.1.

Fermionic and gauge parts
Since the THDM extension of the SM does not affect the gauge and the fermion parts of the Lagrangian, the renormalization of these parts is identical to the SM case. It is described in detail in Ref. [48] in BHS convention, which is included in the standard implementation of the FeynArts package [44]. Other renormalization prescriptions can, e.g., be found in Refs. [61,62]. Therefore, we here do not repeat the renormalization procedure of the CKM matrix, which does not change in the transition from the SM to the THDM, and spell out the renormalization of the fermionic parts only for the case where the CKM matrix is set to the unit matrix, i.e. V ij = δ ij . The transformation of the left-and right-handed fermions and of the gauge-boson fields are where the bare photon field is denoted A bare to distinguish it from the neutral CP-odd field A 0 .
Mixing between left-handed up-and down-type fermions does not occur, owing to charge conservation. Inserting this into the Lagrangian directly delivers the renormalized and the counterterm Lagrangians.

Yukawa part
The renormalization of the Yukawa sector is straightforward in Type I, II, lepton-specific, and flipped models and can be done by taking the Lagrangian of Eq. (2.28), replacing the vev v, and applying the renormalization transformations of Sect. 3.1.2, as well as a renormalization of the fermion masses, The corresponding counterterm couplings are stated in App. B.2.

Renormalization conditions
The renormalization constants are fixed using on-shell conditions for all parameters that are accessible by experiments. However, not all parameters of the THDM correspond to measurable quantities, so that we renormalize three parameters of the Higgs sector in the MS scheme, where the renormalization constants only contain the standard UV divergence in D = 4−2ǫ dimensions and with the Euler-Mascheroni constant γ E . In Sect. 4.2, four different options, resulting in four different renormalization schemes, are presented. An overview over the renormalization constants introduced in the previous section is shown in Tab. 3. In the following, we adapt the notation of Ref. [48], i.e. we use the same symbols for the renormalized and the corresponding unrenormalized Green function, self-energies, etc., but denoting the renormalized quantities with a caret.

Higgs sector
Tadpoles: We start with the (irreducible) renormalized one-point vertex functionŝ At NLO the renormalized tadpoleT consists of a counterterm contribution δt and an unrenormalized one-loop irreducible one-point vertex function T resulting from the diagrams shown in Fig. 1 (12):  two contributions cancel each other, which means that explicit tadpole diagrams can be omitted from the set of one-loop diagrams for any process. However, as a remnant of the tadpole diagrams the tadpole counterterms appear also in various coupling counterterms and need to be calculated. It should be noted that the condition on the tadpoles does not affect physical observables as long as physically equivalent renormalization conditions are imposed on the input parameters. This is, in particular, the case for on-shell renormalization, where input parameters are tied to measurable quantities. That means, changing the tadpole renormalization condition shifts contributions between Green functions and counterterms and merely changes the bookkeeping, but the dependence of predicted observables on renormalized input parameters remains the same. The situation changes if an MS renormalization condition is used, where the counterterm is not fixed by a measurable quantity, but by a divergence in a specific Green function, so that the gauge-dependent tadpole terms can affect the relation between renormalized input parameters and observables. The gauge-independent treatment of tadpole contributions is based on a different renormalization condition and discussed in Sect. 4.2.3.
Scalar self-energies: For scalars, the irreducible two-point functions with momentum transfer k areΓ where both fields a, b are incoming and a, b = H, h, A 0 , G 0 , H ± , G ± . The first term is the LO two-point vertex function, while the functionsΣ ab are the renormalized self-energies containing loop diagrams and counterterms. Generic diagrams contributing to the self-energies are shown in Fig. 2. Mixing occurs only between H and h, between A 0 and G 0 , and between H ± and G ± . For the neutral CP-even fields we obtain and for the CP-odd fieldŝ The charged sector involves the following self-energies, The mass mixing constants are given in Eqs. (3.10) and (3.14). On these two-point functions we now impose our renormalization conditions. First, we fix the renormalized mass parameters to the on-shell values, so that the zeros of the real parts of the one-particle-irreducible two-point functions are located at the squares of the physical masses: Using Eqs. (4.5), (4.6), and (4.7), fixes the mass renormalization constants to For the propagators of the fields, we demand that the residues of the particle poles are not changed by higher-order corrections. This determines the diagonal field renormalization constants by conditions on the one-particle-irreducible two-point functions, which implies where we introduced Σ ′ (k 2 ) as the derivative w.r.t. the argument k 2 . To fix the mixing renormalization constants, we enforce the condition that on-mass-shell fields do not mix, i.e.
After inserting the renormalized self-energies we obtain Since Goldstone-boson fields do not correspond to physical states, we do not render Green functions with external Goldstone bosons finite, so that we need not fix the constants δZ G 0 , δZ A 0 G 0 , δZ G + , δZ HG + ; we could even set them to zero consistently. The possible ZA 0 and W ± H ∓ mixings vanish for physical on-shell gauge bosons due to the Lorentz structure of the two-point function and the fact that polarization vectors ε µ are orthogonal to the corresponding momentum. Using the convention where all fields are incoming and k is the incoming momentum of the gauge-boson fields, the vector-scalar mixing self-energies obey The mixing self-energies on the other on-shell points k 2 = M 2 A 0 and k 2 = M 2 H + , respectively, are connected to the mixing of A 0 or H ± with the Goldstone-boson fields of the Z or the W boson and can be calculated from a BRST symmetry [63]. The BRST variation of the Green functions of one anti-ghost and a Higgs field implies Slavnov-Taylor identities. While the variation of the anti-ghost fields yields the gaugefixing term, the variation of the Higgs fields introduces ghost contributions which vanish for on-shell momentum resulting in 5 We have verified these identities analytically and numerically. Together with the renormalization condition of Eq. (4.12) we conclude that This set of renormalization conditions ensures that no on-shell two-point vertex function obtains any one-loop corrections, and the corresponding external self-energy diagrams do not have to be taken into account in any calculation.

Electroweak sector
The fixing of the renormalization constants of the electroweak sector is identical to the SM case. The mass renormalization constants are fixed in such a way that the squares of the masses correspond to the (real parts of the) locations of the poles of the gauge-boson propagators. The field renormalization constants are fixed by the conditions that residues of on-shell gauge-boson propagators do not obtain higher-order corrections, and that on-shell gauge bosons do not mix. For a better bookkeeping we also keep the dependent renormalization constants δc W and δs W in our calculation. This results in [48] The electric charge e is defined via the eeγ coupling in the Thomson limit of on-shell external electrons and zero momentum transfer to the photon, which yields in BHS convention [48] (4.23) 5 A particularly simple, alternative way to derive these identities is to exploit the gauge invariance of the effective action in the background-field gauge, as done in Ref. [43] for the SM. The respective Ward identities for the background fields differ from the Slavnov-Taylor identities only by off-shell terms, which vanish on the particle poles. Generalizing the derivation of Ref. [43] to the THDM and adapting the results to our conventions for self-energies, the desired Ward identities for the unrenormalized background fields read TĤs β−α − Tĥc β−α , (4.17)

Fermions
The renormalization conditions for the fermions are identical to the ones in the SM, described in detail in Ref. [48]. We demand that the (real parts of the locations of the) poles of the fermion propagators correspond to the squared fermion masses, and that on-shell fermion propagators do not obtain loop corrections. Assuming the CKM matrix equal to the unit matrix, the results for the renormalization constants simplify to where we have used the usual decomposition of the fermion self-energies into a left-handed, a right-handed, and a scalar part, Σ f,L i , Σ f,R i , and Σ f,S i , respectively. The expressions for a non-trivial CKM matrix can be found in Ref. [48].

MS renormalization conditions
In the four renormalization schemes we are going to present, the imposed on-shell conditions are identical, and the differences only occur in the choice of different MS conditions. The parameters α or λ 3 governing the mixing of the CP-even Higgs bosons, and the parameters β and λ 5 need to be fixed. A formulation of an on-shell condition for these parameters is not obvious. One could relate the parameters to some physical processes, such as Higgs-boson decays, and demand that these processes do not receive higher-order corrections. However, so far, no sign of further Higgs bosons has been observed, hence, there is no distinguished process, and such a prescription does not only require more calculational effort, but could introduce artificially large corrections to the corresponding parameters, which would spread to many other observables, as discussed in Refs. [39,64]. Therefore, we choose to renormalize these parameters within the MS scheme, though different variables (such as α or λ 3 ) can be chosen to parameterize the model. Imposing an MS condition on either of the parameters leads to differences in the calculation of observables. In addition, gauge-dependent definitions of MS-renormalized parameters spoil the gauge independence of the relations between input parameters and observables. However, gauge dependences might be even acceptable if the renormalization scheme yields stable results and a good convergence of the perturbation series. The price to pay is that subsequent calculations should be done in the same gauge or properly translated into another gauge. We will discuss different renormalization schemes based on different treatments of α or λ 3 parameterizing the CP-even Higgs-boson mixing, of the parameter β, and of the Higgs coupling constant λ 5 in the following. We begin with the so-called MS(α) scheme.

MS(α) scheme
In this scheme the independent parameter set is {p mass } of Eq. (2.25), so that the parameters β, α, and λ 5 are renormalized in MS. The corresponding counterterm Lagrangian was derived in Sect. 3.1.2.
The renormalization constant δβ: The renormalization constant δ tan β = δβ/c 2 β of the mixing angle β is related to the renormalization constants of the vevs by demanding the defining relation tan β = v 2 /v 1 for bare and renormalized quantities. In MS, δβ can be most easily calculated using the minimal field renormalization (3.17) with the following renormalization transformation of the vevs, v 1,0 = Z 1/2 Using the well-known relation [65] δv 1 /v 1 − δv 2 /v 2 = finite, (4.26) the general form of δβ in the MS scheme where UV indicates that we take only the UV-divergent parts, which are proportional to the standard divergence ∆ UV . The explicit calculation of the UV-divergent terms of δZ h , δZ H according to Eqs. (4.11) in 't Hooft-Feynman gauge reveals that only diagrams with closed fermion loops contribute to the counterterm, with the colour factors c quark = 3, c lepton = 1 and the coupling coefficients ξ f A 0 as defined in Tab. 2. In the class of R ξ gauges this result is gauge independent at one-loop order [39,40].
Higgs self-coupling: The Higgs self-coupling counterterm δλ 5 has to be fixed via a vertex correction. We define this renormalization constant in MS as well, as there is no distinguished process to fix it on-shell. Any 3-or 4-point vertex function with external Higgs bosons is suited to calculate the divergent terms. Since the HA 0 A 0 vertex correction involves fewest diagrams, it is our preferred choice. The condition iŝ Solving this equation for δλ 5 fixes this renormalization constant. The generic one-loop diagrams appearing in this vertex correction are shown in Fig. 3, the contribution of the diagrams involving closed fermion loops is The diagrams containing only bosons lead to Since λ 5 is a fundamental parameter of the Higgs potential, an MS definition leads to a gaugeindependent counterterm.

MS(λ 3 ) scheme
In this scheme, the independent parameter set is {p ′ mass } defined in Eq. (2.24). The renormalization of β and λ 5 is identical to the previous renormalization scheme and not stated again, but now the parameter λ 3 (instead of α) is an independent parameter being renormalized in MS. This has the advantage that this parameter is gauge independent, as it is a defining parameter Figure 3: Generic diagrams contributing to the HA 0 A 0 vertex correction used for the renormalization of λ 5 .
of the basic parameterization of the Higgs potential and thus is safe against potentially gaugedependent contributions appearing in relations between bare parameters. As stated above, the MS renormalization of the parameter β generally breaks gauge independence, but in R ξ gauges the gauge dependence cancels at one loop [39,40], so that this scheme yields gauge-independent results at NLO. We take the counterterm potential of Sect. 3.1.2, but treat δα as a dependent counterterm. As α is a pure mixing angle, we choose to apply the renormalization prescription of Sect. 3.1.2, where the mixing angle diagonalizes the potential to all orders. The relation between δα and the independent constants is given in Eq. (3.9) with δM 2 Hh = 0, where f α ({δp ′ mass }) can be obtained from Eq. (2.16c) by applying the renormalization transformation of Eq. (A.2) (which is identical to the renormalization transformation of Sect. 3.1.2, but renormalizing λ 3 instead of α). This yields The UV-divergent term of δα has been calculated in Eq. (4.31), and by renormalizing δλ 3 in MS scheme, it is clear that the dependent δα must now have a finite part in addition. We choose this finite term in such a way that the finite part in δλ 3 (which results from δλ 3 by setting ∆ UV to zero) vanishes and obtain where δλ 3 drops out as it has no finite part. The divergent part of δλ 3 can be calculated by solving Eq. (4.37) and using the knowledge about the divergent parts of δα from Eqs. (4.32) and (4.33). This results in (4.40) The fermionic contribution to δλ 3 is given by with the massive fermions f = e, . . . , t and the generation index i. For the bosonic contribution we obtain (4.42)

The FJ tadpole scheme
Since tadpole loop contribution T S are gauge dependent [66], the connection among bare parameters potentially becomes gauge dependent if δt S = −T S enters the relations between bare parameters, as it is the case if renormalized tadpole parameters t S are forced to vanish. Note that these gauge dependences systematically cancel if on-shell renormalization conditions are employed, i.e. if predictions for observables are parameterized by directly measurable input parameters. If some input parameters are renormalized in the MS scheme this cancellation of gauge dependences does not take place anymore in general, and the gauge dependence is manifest in relations between predicted observables and input parameters at NLO. In the MS(α) and the MS(λ 3 ) renormalization schemes, the bare definitions of α and β contain tadpole terms leading to a gauge dependence (although the MS(λ 3 ) scheme is gauge independent at NLO in R ξ gauges). Fleischer and Jegerlehner [67] proposed a renormalization scheme for the SM, referred to as the FJ scheme in the following, that preserves gauge independence for all bare parameters, including the masses and mixing angles 6 . In this scheme, the parameters are defined in such a way that tadpole terms do not enter the definition of any bare parameter so that all relations among bare parameters remain gauge independent. This can be achieved by demanding that bare tadpole terms vanish, t S,0 = 0, for all fields S with the quantum numbers of the vacuum. Since tadpole conditions have no effect on physical observables and change only the bookkeeping, such a procedure is possible. The disadvantage is that now tadpole diagrams have to be taken explicitly into account in all higher-order calculations. In particular, the one-particle reducible tadpole contributions destroy the simple relation between propagators and two-point functions. In the SM, the FJ scheme does not affect observables if all parameters are renormalized using onshell conditions-as usually done-except for the strong coupling constant α s , which is, however, directly related to the strong gauge coupling, a model defining parameter. A gauge-independent renormalization scheme for the THDM can be defined by applying the FJ prescription and imposing the MS condition on mixing angles [39,40]. The bare physical parameters defined in the FJ scheme differ by NLO tadpole contributions (including divergent and finite terms) from the gauge-dependent definition of the bare parameters {p mass } given in Eq. (2.25). Exceptions are e and the parameter λ 5 , which is a parameter of the basic potential and therefore gauge independent by construction. The renormalization of λ 5 in MS is identical to the one in the previous schemes.
It should be noted that in Refs. [39,40] m 2 12 is chosen as independent parameter in contrast to our choice of λ 5 . The latter, however, is closer to common practice used in the MSSM [69,70]. Moreover, in Refs. [39,40] tadpole counterterms are reintroduced by shifting the Higgs fields according to η i → η i + ∆v η,i , i = 1, 2, where the constants ∆v η,i can be chosen arbitrarily, since physical observables do not depend on this shift, which can be interpreted as an unobservable change of the integration variables in the path integral. In Refs. [39,40], this freedom of choice is exploited, and ∆v η,i are chosen in such a way that the fields η 1 , η 2 do not develop vevs at all orders. This affects the form of the counterterm Lagrangian and the definition of the renormalization constants with the consequence that the formulae given in Eq. (3.16) and Sect.4.1 cannot be applied.
We have implemented the FJ scheme following the strategy of Ref. [40] by performing the shifts η i → η i + ∆v η,i and in an alternative, simpler (but physically equivalent) way. In this simplified approach we keep the dependence of the Lagrangian in terms of gauge-dependent masses and couplings. In addition we keep the tadpole renormalization condition (4.3), so that the definitions of the renormalization constants of the on-shell parameters and the Z factors according to Sect. 4.1 remain valid (otherwise we needed to take into account actual tadpole diagrams everywhere). In this simplified approach the counterterms for α and β which reproduce the results in the FJ scheme result from the previously derived δα and δβ by adding appropriate finite terms, δα FJ = δα + finite terms, δβ FJ = δβ + finite terms, (4.43) which depend on the (finite parts of the) tadpole contributions T H and T h .
Before performing the full calculations, we outline the strategy of the derivation of those finite terms for β; for α everything works analogously. We start by exploiting the fact that the form of the tadpole renormalization cannot change physical results if all counterterms for independent parameters are determined by the same physical conditions. This means, as mentioned above, that we can simply define the bare tadpoles to vanish, but this forces us to include all explicit tadpole contributions to Green functions. We indicate quantities in this variant by a superscript "t" in the following. We get the same physical predictions in this "t-variant" if we use the counterterm instead of δβ, where δβ t is calculated in the same way as δβ, but with tadpole counterterms omitted and explicit tadpole diagrams (including divergent and finite parts) in the occurring Green functions taken into account. Note that the MS prescription to include only divergent terms, which is employed to define δβ, is not applied to the new tadpole contribution ∆β t (T H , T h ). Otherwise the new ∆β t terms could not be fully compensated by explicit tadpole contributions occuring elsewhere, so that there would be differences in the renormalized amplitudes. In fact, applying the MS prescription to ∆β t (T H , T h ) as well defines the FJ renormalization scheme, The quantity δβ t FJ is the gauge-independent counterterm for β introduced in Ref. [40] which is to be used in the t-variant, where all explicit tadpole diagrams are included in Green functions (or equivalently are redistributed by the ∆v shift as described Ref. [40]). We can translate the FJ renormalization prescription back to our renormalization scheme (with vanishing renormalized tadpoles) by the counterpart of Eq. (4.44), but now formulated in the FJ scheme, (4.46) i.e. δβ FJ is the counterterm for β to be used in our counterterm Lagrangian in order to calculate renormalized amplitudes in the FJ scheme. Combining the above formulas, we obtain the finite difference between δβ in the (gauge-dependent) MS scheme and δβ FJ in the (gaugeindependent) FJ scheme, (4.47) The renormalization constant δβ FJ : We begin our calculation of δβ| FJ with an alternative computation of δβ in the MS(α) scheme, because Eq. (4.26) cannot be applied in the FJ scheme.
To avoid the use of Eq. (4.26), we calculate the counterterm in the MS(α) scheme from the field renormalization constant by employing the last two equations of Eq. (3.18). This results in with δM 2 G 0 A 0 as given in Eq. (3.14). In the second step, δZ G 0 A 0 from Eq. (4.13) and have been used. Equation (4.49) results from demanding finiteness of the G 0 A 0 mixing selfenergy at zero-momentum transfer, k 2 = 0, but actually any other value of k 2 would be possible as well, since we only have to remove all UV-divergent terms in the mixing. The non-vanishing tadpole counterterms in the MS(α) scheme are δt S = −T S . In the transition to the t-variant, δβ gets modified by two kind of terms: First, there are no tadpole counterterms, i.e. the δM 2 G 0 A 0 term is absent, and second, there are explicit tadpole contributions to Σ G 0 A 0 . This implies where the superscript "t" indicates that one-particle-reducible tadpole diagrams are included in the self-energies. The subscript "T H , T h " means that only the those explicit tadpole contributions are taken into account here. Inserting δM 2 G 0 A 0 from Eq. (3.14) and evaluating the (momentum-independent) tadpole diagrams for the G 0 A 0 mixing, ∆β t evaluates to The counterterm δβ t FJ of the FJ scheme in the t-variant, thus, reads which is in agreement with Ref. [39,40]. This translates to our treatment of tadpoles as where again "finite" means that ∆ UV is set to zero in the tadpole contribution. Using this counterterm, it is possible to keep the form of the counterterm Lagrangian derived in Sect. 3 to obtain results in the gauge-independent FJ scheme, although the above counterterm Lagrangian employs a gauge-dependent (but very convenient) tadpole renormalization. Alternatively, δβ could be fixed by an analogous consideration of the ZA 0 mixing, leading to which is independent of k 2 and does not use relation (4.26). The transition to the FJ scheme then simply amounts to replacing the one-particle-irreducible self-energy Σ ZA 0 by Σ t,ZA 0 , which includes tadpole diagrams. The result δβ t FJ of this procedure is again given by Eq. (4.52), as it should be.
The renormalization constant δα FJ : We apply the same method to the renormalization constant δα FJ , starting from Eq. (4.31). The difference between δα and δα t is entirely given by the explicit tadpole diagrams that appear in the change from Σ Hh to Σ t,Hh in Eq. (4.31), which evaluates to with the coupling factors of the hHH and hhH vertices The counterterm δα t FJ of the FJ scheme in the t-variant, thus, reads which is again in agreement with Ref. [39,40]. This translates to our treatment of tadpoles as Concerning the use of δα FJ in our counterterm Lagrangian to obtain renormalized amplitudes in the gauge-independent FJ scheme, the same comments made above for δβ FJ apply.

The FJ(λ 3 ) scheme
In the MS(λ 3 ) scheme, the parameters λ 3 and λ 5 are defining parameters of the basic parameterization and gauge independent by construction. Therefore, the condition on δβ is the only renormalization condition potentially being gauge dependent. To provide a fully gaugeindependent renormalization scheme where λ 3 is an independent quantity, we apply the FJ scheme to the parameter β and keep the renormalization of λ 3 and λ 5 as in the MS(λ 3 ). We call the resulting scheme the FJ(λ 3 ) scheme. The renormalization of the parameters reads:

Conversion between different renormalization schemes
In the previous section, we have presented four different renormalization schemes, which treat the mixing parameters differently. When observables calculated in different renormalization schemes are compared, particular care has to be taken that the input parameters are consistently translated from one scheme to the other. The bare values of identical independent parameters are equal and independent of the renormalization scheme. Exemplarily, for a parameter p, the renormalized values p (1) and p (2) in two different renormalization schemes 1 and 2 are connected via the bare parameter p 0 , within the considered order. If p is a dependent parameter in one or both schemes, it must be calculated from the independent renormalized parameters and their counterterms from the relations between bare and renormalized quantities. For converting an input value from one scheme to another, one can solve for one renormalized quantity p (1) = p (2) + δp (2) (p (2) ) − δp (1) (p (1) ). The MS parameters are defined at the scale The motivation for this choice will become clear below. The remaining input parameters for the SM part are given in App. C. In both plots, we highlight the phenomenologically relevant region in the centre. The solid lines are the result obtained by solving the implicit equations (4.61) numerically, the dashed lines correspond to a linearized conversion. All curves show only minor conversion effects in the parameter values, i.e. the solution of the implicit equations agrees well with the approximate linearized conversion, affirming that the contributions of the higherorder functions ∆α t , ∆β t , and f α of Eqs. (4.64a)-(4.64c) are small, and perturbation theory is applicable. Since the values of the parameters change when going from one renormalization scheme to another, the alignment limit does not persist in these transformations, i.e. in this scenario the alignment limit sensitively depends on the definition of the parameters at NLO. For the schemes with λ 3 as input parameter, some singular behaviour in the parameter conversion can be observed in the phenomenologically disfavoured region where c β−α < ∼ − 0.3.
This artifact in the conversion appears when c 2α → 0 (see, e.g., Eq. (4.38)), indicating the breakdown of the MS(λ 3 ) and FJ(λ 3 ) schemes in such parameter regions. Already this casespecific study shows that stability issues of different renormalization schemes have to be carefully carried out for all interesting parameter regions and that the applicability of a specific scheme in general does not cover the full THDM parameter space, a fact that was also pointed out in Ref. [39] for the THDM and that is known from NLO calculations in the MSSM (see, e.g., Refs. [73,74]). Specifically, if the MS(λ 3 ) and FJ(λ 3 ) schemes are not applicable in a region that might be favoured by future data analyses, it would be desirable and straightforward to replace λ 3 by λ 1 or λ 2 as independent parameter, thereby defining analogous schemes like MS(λ 1 ), etc.. To address this issue properly, was our basic motivation to introduce and compare different renormaliztion schemes. We will continue this discussion in more detail in a forthcoming publication, where further THDM scenarios are considered.

The running of the MS parameters
Parameters renormalized in the MS scheme depend on an unphysical renormalization scale µ r . The one-loop β-function of a parameter p can be obtained from the UV-divergent parts of its counterterm δp, Since the renormalization constants are computed in a perturbative manner, the β-functions have a perturbative expansion in the coupling parameters. Note that the last equality in Eq. (5.1) holds in the FJ schemes for α and β only in the t-variant explained above, because the finite contributions ∆α t and ∆β t depend on the scale µ r . As discussed in the previous sections, the ratio of the vevs, tan β, the Higgs mixing parameter α or λ 3 , and the Higgs self-coupling λ 5 are renormalized in the MS scheme. For each renormalization scheme described in Sect. 4.2, one obtains a set of coupled RGEs involving the β-functions of the independent parameters. Therefore, the scale dependence varies when different schemes are applied. In the perturbative expansion of the β-function we consider only the one-loop term, being second order in the coupling constants, e.g., in the MS(α) scheme The dependence on the strong coupling constant vanishes at one-loop order as the parameters renormalized in MS appear only in couplings of particles that do not interact strongly. The coefficients A p , B p , C p of the respective renormalized parameter can be easily read from the divergent terms which have been derived in the previous section. We have checked them against the β-functions given for λ 3 and λ 5 in Ref. [50] and for β in Ref. [65] (supersymmetric contributions need to be omitted).
In general, RGEs, which are a set of coupled differential equations, cannot be solved analytically. Usually numerical techniques, such as a Runge-Kutta method, need to be employed to solve the RGEs and to compute the values of the parameters at a desired scale. Moreover, we emphasize that the renormalization-group flow of a running parameter depends on the renormalization scheme of the full set of independent parameters. That means the fact that we use on-shell quantities, such as all the Higgs-boson masses, to fix most of the scalar self-couplings has a significant impact on the running of our MS parameters. The renormalization-group flow in other schemes was, e.g., investigated in Refs. [50,[75][76][77][78].
The scale dependence of c β−α for µ = 100−900 GeV is plotted in Fig. 5, for the scenario defined in Eq. scheme has a large impact on the scale dependence. While the MS(α) scheme introduces only a mild running, the other schemes show a much stronger scale dependence, so that excluded and unphysical values of input parameters can be reached quickly. A similar observation has also been made in supersymmetric models for the parameter tan β [64]. Gauge-dependent MS schemes have a small scale dependence while replacing the parameters by gauge-independent ones like in the FJ schemes introduce additional terms in the β-functions which induce a stronger scale dependence. In Fig 5(b) one can also see that the curves for the MS(λ 3 ) and the FJ(λ 3 ) schemes terminate around 250 GeV. At this scale, the running of λ 3 yields unphysical values for which Eq. (2.20c) with the given Higgs masses becomes overconstrained, and no solution with |s 2α | ≤ 1 exists. This is unique to the λ 3 running as only there an implicit equation needs to be solved to obtain the input parameter α. For the other cases we prevent the angles from running out of their domain of definition by solving the running for the tangent function of the angles.

Implementation into a FeynArts Model File
The Mathematica package FeynRules (FR) [79] is a tool to generate Feynman rules from a given Lagrangian, providing the possibility to produce model files in various output formats which can be employed by automated amplitude generators. We have inserted the Lagrangian into FR in its internal notation to obtain the corresponding counterterm Lagrangian after the renormalization transformations. Before this insertion, we have computed and simplified the Higgs potential and the corresponding counterterm potential (3.16) with inhouse Mathematica routines. Using FR, the tree-level and the counterterm Feynman rules as well as the renormalization conditions in the MS(α) and MS(λ 3 ) schemes have been implemented into a model file for the amplitude generator FeynArts (FA) [44]. The renormalization conditions of the FJ(α) and the FJ(λ 3 ) have not been included in the model file, because using Eqs. (4.62), (4.63) it is straightforward to implement the corresponding finite terms of δα and δβ. With such a model file, NLO amplitudes for any process can be generated in an automated way. The FA NLO model file for the THDM, obtained with FR, has the following features: • Type I, II, flipped, or lepton-specific THDM; • all tree-level and counterterm Feynman rules; • renormalization conditions according to the MS(α) and MS(λ 3 ) schemes; • all renormalization constants are implemented additionally in MS as well, which allows for fast checks of UV-finiteness; • BHS and HK conventions; • CKM matrix set to the unit matrix (the generalization is straightforward).
This model file has been tested intensively, including checks of UV-finiteness for several processes, both numerically and analytically. This allows for the generation of amplitudes (and further processing with FormCalc [80]) for any process at the one-loop level, at any parameter point of the THDM. The model file can be obtained from the authors upon request.

Numerical results for h → WW/ZZ → 4f
In this section we present first results from the computation of the decay of the light, neutral CP-even Higgs boson of the THDM into four fermions at NLO. The computer program Prophecy4f [81-83] 8 provides a "PROPer description of the Higgs dECaY into 4 Fermions" and calculates observables for the decay process h→WW/ZZ→4f at NLO EW+QCD in the SM. We have extended this program to the calculation of the corresponding decay in the THDM in such a way that the usage of the program and its applicability as event generator basically remains the same. Owing to the fact that LO and real-emission amplitudes in the THDM receive only the multiplicative factor s β−α with respect to the SM, the bremsstrahlung corrections as well as the treatment of infrared singularities could be taken over from the SM calculation [81,83] via simple rescaling. The calculation in the THDM, the implementation in Prophecy4f, as well as results of the application will be described in detail in an upcoming publication. We just mention that we employ the complex-mass scheme [84] to describe the W/Z resonances, as already done in Refs. [81][82][83] for the SM. Note that the W/Z-boson masses as well as the weak mixing angle are consistently taken as complex quantities in the complex-mass scheme to guarantee gauge invariance of all amplitudes in resonant and non-resonant phase-space regions. Consequently all our renormalization constants of the THDM inherit imaginary parts from the complex input values, but the impact of these spurious imaginary parts is beyond NLO and negligible (as in the SM). Moreover, we mention that the modified version of Prophecy4f makes use of the public Collier library [85] for the calculation of the one-loop integrals. Apart from performing two independent loop calculations, we have verified our one-loop matrix elements by numerically comparing our results to the ones obtained in Ref. [40] for the related Wh/Zh production channels (including W/Z decays) using crossing symmetry.
In this paper, we present first results in order to demonstrate the use and the self-consistency of our renormalization schemes, employing again the scenario inspired by the first benchmark scenario of Ref. [72] where the additional Higgs bosons are not very heavy. The input values of the THDM parameters for a Type I THDM are given in Eqs. (4.65) and (4.66). Since c β−α is the only parameter of the THDM appearing at LO, our process is most sensitive to this parameter. We vary c β−α in the range [−0.2, +0.2] in scenario A for the computation of the partial decay width for h → WW/ZZ → 4f , Γ h→4f THDM , which is obtained by summing the partial widths of the h boson over all massless four-fermion final states 4f . The parameters of the SM part of the THDM are collected in App. C. Note that a non-trivial CKM matrix would not change our results, since quark mass effects of the first two generations as well as mixing with the third generation are completely negligible in the considered decays.
To perform scale variations we take two distinguished points named Aa and Ab with c β−α = ±0.1. For the central renormalization scale we use the average mass µ 0 defined in Eq. (4.67) of all scalar degrees of freedom. The scale µ of α s is kept fixed at µ = M Z which is the appropriate scale for the QCD corrections (which are dominated by the hadronic W/Z decays).

Scale variation of the width
The running of the MS-renormalized parameters α and β is induced by the Higgs-boson selfenergies (and some scalar vertex for λ 5 ), i.e. the relevant particles in the loops are all Higgs bosons, the W/Z bosons, and the top quark. If all Higgs-boson masses are near the electroweak scale, say ∼ 100−200 GeV, where the W/Z-boson and top-quark masses are located, then the scale M h turns out to be a reasonable scale, as expected. However, if some heavy Higgs-boson masses increase to some generic mass scale M S and the mixing angle β − α stays away from the alignment limit, there is no decoupling of heavy Higgs-boson effects, so that M S acts as generic UV cutoff scale appearing in logarithms log(M S /µ r ). The renormalization scale µ r has to go up with M S to avoid that the logarithm drives the correction unphysically large. The optimal choice of µ r , though, is somewhat empirical. A good choice of the central scale µ 0 should come close to the stability point (plateau in the µ r variation) in the major part of THDM parameter space. Our choice (4.67) of µ 0 effectively takes care of this and is eventually justified by the numerics.
To illustrate this and to estimate the theoretical uncertainties due to the residual scale dependence, we compute the total width while the scale µ r is varied from 100−900 GeV. Results with central scale M h are shown in App. D, proving that this would be not a good choice. The parameters α and β are defined in the MS(λ 3 ) scheme, and to compute results in other renormalization schemes their values are converted using Eqs. (4.64a)-(4.64c), which are solved numerically without linearization. Thereafter the scale is varied, the RGEs solved, and the width computed using the respective renormalization scheme. The results are shown in Fig. 6 at LO (dashed) and NLO EW (solid) for the benchmark points Aa and Ab. The QCD corrections are not part of the EW scale variation and therefore omitted in these results. The benchmark point Aa shows almost textbook-like behaviour with the LO computation exhibiting a strong scale dependence for all renormalization schemes, resulting in sizable differences between the curves. However, each of the NLO curves shows a wide extremum with a large plateau, reducing the scale dependence drastically, as it is expected for NLO calculations. The central scale µ r = (M h + M H + M A 0 + 2M H + )/5 lies perfectly in the middle of the plateau regions motivating this scale choice. In contrast, the naive scale choice µ 0 = M h is not within the plateau region, leads to large, unphysical corrections, and should not be chosen. The breakdown of the FJ(α) curve for small scales can be explained by the running which becomes unstable for these values (see Fig. 5(a)). For all renormalization schemes, the plateaus coincide and the agreement between the renormalization schemes is improved at NLO w.r.t. the LO results. This is expected, since results obtained with different renormalization schemes should be equal up to higher-order terms, after the input parameters are properly converted. The relative renormalization scheme dependence at the central scale,   renormalization scale, which attests a good absorption of further corrections into the NLO prediction.
The situation for the benchmark point Ab is more subtle. For negative values of c β−α the truncation of the schemes involving λ 3 at µ r = 250−300 GeV as well as the breakdown of the running of the FJ(α) scheme, which both were observed in the running in Fig. 5(b), are also manifest in the computation of the h→4f width. Therefore, the results vary much more, and the extrema with the plateau regions are not as distinct as for the benchmark point Aa. They are even missing for the truncated curves. Nevertheless, the situation improves at NLO. As for scenario Aa, the central scale choice of µ 0 is more appropriate in contrast than the choice of M h . For both benchmark points, the estimate of the theoretical uncertainties by varying the scale by a factor of two from the central value for an arbitrary renormalization scheme is generally not appropriate. A proper strategy would be to identify the renormalization schemes which yield reliable results, and to use only those to quantify the theoretical uncertainties from the scale variation. In addition, the renormalization scheme dependence of those schemes should be investigated. This procedure should be performed for different parameter regions (and corresponding benchmark points) separately, which is beyond the scope of this work.

c β−α dependence
The decay width for h → 4f in dependence of c β−α in scenario A is presented in Fig. 7 for all renormalization schemes with the input values α and β defined in the MS(λ 3 ) scheme. The LO (dashed) and the full NLO EW+QCD total widths (solid) are computed in the different renormalization schemes after the NLO input conversion (without linearization) and using the constant default scale µ 0 of Eq. (4.67). The SM values are illustrated in red. At tree level the widths show the suppression w.r.t. to the SM with the factor s 2 β−α originating from the HWW and HZZ couplings. The differences between the renormalization schemes are due to the conversion of the input. As the conversion induces NLO differences in the LO results, a pure LO computation is identical for all renormalization schemes as the conversion vanishes at this order and is represented by the LO curve of the MS(λ 3 ) scheme. The suppression w.r.t. the SM computation does not change at NLO, while the shape becomes slightly asymmetric, and the NLO results show a significantly better agreement between the renormalization schemes. Deviations of the THDM results from the SM expectations can be investigated when the SM Higgs-boson mass is identified with the mass M h of the light CP-even Higgs boson h of the THDM. The relative deviation of the full width from the SM is then which is shown in Fig. 8 at LO (dashed) and NLO (solid) in percent for parameters defined in the the MS(λ 3 ) scheme. The SM exceeds the THDM widths at LO and NLO. The LO shape which is just given by c 2 β−α shows minor distortions due to the parameter conversions. At NLO, the shape is slightly distorted by an asymmetry of the EW corrections, and a small offset of −0.5% is visible even in the alignment limit where the diagrams including heavy Higgs bosons still contribute. The NLO computations show larger negative deviations, and this could be used to improve current exclusion bounds or increase their significance. Nevertheless, in the whole scan region the deviation from the SM is within 6% and for phenomenologically most interesting region with |c β−α | < 0.1 even less than 2%, which is challenging for experiments to measure.

Conclusions
Confronting experimental results on Higgs precision observables with theory predictions within extensions of the SM, provides an important alternative to search for physics beyond the SM, in addition to the search for new particles. The THDM comprises an extended scalar sector with regard to the SM Higgs sector and allows for a comprehensive study of the impact of new scalar degrees of freedom without introducing new fundamental symmetries or other new theoretical structures.
In this article, we have considered the Type I, II, lepton-specific, and flipped versions of the THDM. We have introduced four different renormalization schemes which employ directly measurable parameters such as masses as far as possible and make use of fields that directly correspond to mass eigenstates. In all the schemes, the masses are defined via on-shell conditions, the electric charge is fixed via the Thomson limit, and the coupling λ 5 is defined with the MS prescription. The fields are also defined on-shell which is most convenient in applications. The renormalization schemes differ in the treatment of the coupling λ 3 and the mixing angles α and β: In the MS(α) scheme, α and β are renormalized using MS conditions. In the MS(λ 3 ) scheme instead λ 3 and β are MS-renormalized parameters. In addition to the conventional treatment of tadpole contributions, we have implemented an alternative prescription suggested by Fleischer and Jegerlehner where the mixing angles α and β obtain extra terms of tadpole contributions, rendering these schemes gauge independent to all orders. It should, however, be noted that the MS(λ 3 ) scheme is also gauge independent at NLO in the class of R ξ gauges. We have also discussed relations to renormalization procedures suggested in the literature for the THDM.
A comparison of these four different renormalization schemes allows for testing the perturbative consistency, and, for the parameter regions and renormalization schemes fulfilling this test, estimating the theoretical uncertainty due to the truncation of the perturbation series. To further investigate the latter, we have investigated the scale dependence, solving the corresponding RGEs. One important observation is that it is crucial to be very careful and specific about the definitions of the parameters applied, i.e. the declaration of the renormalization scheme of the parameters is vital if one aims at precision. This is already relevant in the formulation of benchmark scenarios, because a conversion to a different scheme might alter the physical properties of the scenario significantly. For example, the alignment limit may be reached with one specific set of parameters defined in a specific renormalization scheme, but converting these parameters consistently to parameters in a different renormalization scheme might shift the parameters away from the alignment limit.
The different renormalization schemes have been implemented into a FeynArts model file and are thus ready for applications. 9 As a first example, we have applied and tested the different schemes in the calculation of the decay width of a light CP-even Higgs boson decaying into four massless fermions. We discuss the dependence of the total h→4f decay width on the renormalization scale and advocate a scale that is significantly higher than the naive choice of µ r = M h , taking care of the different mass scales in the THDM Higgs sector. In addition, results for various values of cos(β − α), a parameter entering the prediction already at LO, are presented. The deviations of the SM are relatively small, in the phenomenologically interesting region they are about 2−6%-a challenge for future measurements.
The detailed description of the calculation of the decay width in the THDM and a survey of numerical results will be given in a forthcoming paper. This includes a deeper investigation in the renormalization scale dependence and the comparison of different renormalization schemes for more benchmark points as well as differential distributions.
to the mixing angles exist, we can write their behaviour in the renormalization transformation schematically as α 0 = α + 0, β c,0 = β + 0, β n,0 = β + 0. where G µ is the Fermi constant, α s the strong coupling constant at the Z pole, Γ Z and Γ W the total decay widths of the Z and W boson, respectively, and m e , . . . , m b the fermion masses. The W/Z masses are "on-shell masses", which are combined with the W/Z decay widths to complex pole masses; all Higgs-boson and fermion masses are (real) pole masses. The electromagnetic coupling is fixed in the G µ scheme, i.e. calculated from the muon decay constant according to since this choice is appropriate in the NLO calculation for h → 4f . In the G µ scheme, the charge renormalization constant δZ e of Eq. (4.23) receives an additional contribution ∆r, which quantifies the NLO corrections to muon decay (see, e.g., Ref. [86]). The correction ∆r was calculated in the THDM, for instance, in Ref. [87]. For the conversion of THDM parameters between the different renormalization schemes and the calculation of the MS parameter running choosing the G µ scheme plays only a minor role.
D Results for the h → 4f decay width with central renormalization scale M h . In contrast to Fig. 6(a), we observe big discrepancies between the results in the different renormalization schemes (with proper scheme conversion) at LO and NLO, with no tendency of improvement in the transition from LO to NLO. The large differences in the LO predictions at the central scale already signal huge scheme conversion effects due to unnaturally large corrections that cannot be made up by NLO effects. Figure 9(b) shows the respective results without parameter scheme conversion, so that the LO predictions coincide at the central scale and reflect the µ r dependence of s 2 β−α . Lacking the parameter conversion, no reduction of scheme dependence can be expected here. We rather include this figure to check whether and where the different schemes show some reduction of the µ r dependence in the transition from LO to NLO. Such stabilizations are observed at scales about 300−400 GeV, but not near M h = 125 GeV.
Choosing µ 0 = 1 5 (M h + M H + M A 0 + 2M H + ) = 361 GeV, the conversion effects and the NLO corrections, however, are nicely under perturbative control, as discussed in Sec. 7. Note that the results at µ r = 361 GeV neither in Fig. 9(b), nor in Fig. 9(a) correspond to µ 0 = 361 GeV in Fig. 6(a), since the input parameters α, β, and λ 5 are defined at different renormalizations scales µ 0 .  6). In (a) the input for β and α is defined in the MS(λ 3 ) scheme and converted to the other schemes at NLO (also for the LO curves). In (b) the input for β and α is taken without conversion between the schemes, so that c β−α (µ 0 ) = 0.1 in all schemes.