Precision calculations for h->WW/ZZ->4 fermions in the Two-Higgs-Doublet Model with PROPHECY4F

We have calculated the next-to-leading-order electroweak and QCD corrections to the decay processes h ->WW/ZZ ->4 fermions of the light CP-even Higgs boson h of various types of Two-Higgs-Doublet Models (Types I and II,"lepton-specific"and"flipped"models). The input parameters are defined in four different renormalization schemes, where parameters that are not directly accessible by experiments are defined in the MSbar scheme. Numerical results are presented for the corrections to partial decay widths for various benchmark scenarios previously motivated in the literature, where we investigate the dependence on the MSbar renormalization scale and on the choice of the renormalization scheme in detail. We find that it is crucial to be precise with these issues in parameter analyses, since parameter conversions between different schemes can involve sizeable or large corrections, especially in scenarios that are close to experimental exclusion limits or theoretical bounds. It even turns out that some renormalization schemes are not applicable in specific regions of parameter space. Our investigation of differential distributions shows that corrections beyond the Standard Model are mostly constant offsets induced by the mixing between the light and heavy CP-even Higgs bosons, so that differential analyses of h ->4f decay observables do not help to identify Two-Higgs-Doublet Models. Moreover, the decay widths do not significantly depend on the specific type of those models. The calculations are implemented in the public Monte Carlo generator PROPHECY4F and ready for application.


Introduction
The CERN Large Hadron Collider (LHC) was built to explore the validity of the Standard Model (SM) of particle physics at energy scales ranging from the electroweak scale ∼100 GeV up to energies of some TeV and to search for new phenomena and new particles in this energy domain. The discovery of a Higgs particle at LHC Run 1 in 2012 [1,2] was a first big achievement in this enterprise. Since first studies of the properties of this Higgs particle (spin, CP parity, couplings to the heaviest SM particles) show good agreement between measurements and SM predictions, the SM is in better shape than ever to describe all known particle phenomena up to very few exceptions. Assuming the absence of spectacular new-physics signals in LHC data, this means that any deviation from the SM hides in small and subtle effects. To extract those differences from data, both experimental analyses and theoretical predictions have to be performed with the highest possible precision. On the other hand, assuming that a new signal materializes at the 5σ level, the properties of the newly discovered particle have to be investigated with precision, in order to tell different models apart that can accommodate the new phenomenon.
Most of the promising candidates for models beyond the SM modify or extend the scalar sector of electroweak (EW) symmetry breaking, which introduces the Higgs boson in the SM. Lacking clear evidence of the realization of a specific model extension, it is well motivated to prepare experimentally testable predictions within generic SM extensions which are building blocks in larger models. Two-Higgs-Doublet Models (THDMs) [3,4], where a second Higgs doublet is added to the SM field content, provide an interesting class of such generic models. While the gauge structure of the SM is kept, THDMs contain five physical Higgs bosons in contrast to the one of the SM. Three out of the five are neutral and two are charged. In the CPconserving case, considered in this paper, one of the neutral Higgs bosons is CP-odd and two are CP-even. An important issue in identifying or constraining a THDM, thus, consists in telling the SM Higgs boson from an SM-like CP-even Higgs boson of THDMs. To this end, several phenomenological studies in THDMs have been carried out recently by various groups [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21].
In this paper we investigate the decay observables of the SM-like neutral, light, CP-even Higgs boson h decaying into four fermions, h → WW/ZZ → 4f , in the THDM, including nextto-leading-order (NLO) corrections of the EW and strong interactions. The fermions in the final state of these processes can either be quarks or leptons, and especially the latter can be resolved very well in the detector. These four-body decays were already crucial in the Higgsboson discovery, but also play a major role in precision studies of the Higgs boson, in particular, to determine the couplings to the EW gauge bosons W and Z. They also provide a window to physics beyond the SM [22][23][24][25][26][27][28][29], as, due to its high precision, small deviations from the SM can be measured, and differential distributions can be investigated and tested against the SM. The Monte Carlo program Prophecy4f [30][31][32] performs the calculation of the full NLO EW and QCD corrections for all h → WW/ZZ → 4f channels in the complex-mass scheme [33][34][35] to describe the intermediate W-and Z-boson resonances. It provides differential distributions as well as unweighted events for leptonic final states. The corrections to h → ZZ → 4 leptons were also calculated and matched to a QED parton shower in Ref. [36]. In the following, we describe results obtained with an updated version of Prophecy4f, extended to the computation in the THDM in such a way that the usage of the program and its applicability as event generator basically remain the same.
It is our goal to analyze the Higgs decay in the context of the most relevant THDM scenarios. To compute phenomenologically relevant results, we need to take into account current constraints which also restrict the large parameter space. The constraints come from direct LHC searches for heavy Higgs bosons and from theoretical aspects like vacuum stability, perturbative unitarity, or perturbativity of the couplings, which are required for a meaningful perturbative evaluation. The recent report of the LHC Higgs Cross Section Working Group [37] summarizes a selection of relevant benchmark scenarios proposed in other papers among which we study the most relevant. The results are compared with the SM prediction, and deviations are quantified. In addition to the usual investigation of residual scale uncertainties, we compare the results of different renormalization schemes recently presented and discussed in the literature [38][39][40][41]. Specifically, we employ the four different schemes described in Ref. [40] and vary the renormalization scale to investigate the perturbative stability of the predictions in the benchmark scenarios. Similar to what has already been found in the Minimal Supersymmetric SM [42], the different renormalization schemes may suffer from problems like gauge dependence, singularities in relations between parameters, or unnaturally large corrections. The comparison of the results obtained with different renormalization schemes allows us to determine regions where they behave well and yield reliable results. Electroweak corrections to other Higgs-boson decay channels in the THDM were investigated in Refs. [38,43]. Generic tools to calculate Higgs decay widths in the THDM, such as Hdecay [44] and THDMC [45] are currently restricted to QCD corrections (see, e.g., Refs. [46,47] for more details).
This paper is structured as follows: In Sec. 2 we briefly describe the program Prophecy4f, on which our implementation is based, and give some details on our NLO calculation of corrections to the h→4f decays in the THDM, including a survey of Feynman diagrams, the salient features of the calculation, and a short outline of the implementation into Prophecy4f. In Sec. 3, we describe the setup of our numerical analysis and the chosen THDM scenarios. The numerical results are presented and discussed in detail in Sec. 4. We conclude in Sec. 5 and provide some supplementary results in the appendix.
2 Predicting h → WW/ZZ → 4f in the THDM with the Monte Carlo program Prophecy4f

Preliminaries and functionality of Prophecy4f
The Monte Carlo program Prophecy4f [30][31][32] provides a "PROPer description of the Higgs dECaY into 4 Fermions" by calculating the decay observables of the process h → WW/ZZ → 4f at NLO EW+QCD accuracy in the SM. The original Prophecy4f code contains the matrix elements of all 19 possible 4f final states, which are listed in Tab. 1, in a generic way. It takes into account the full off-shell effects of the intermediate W and Z bosons and treats the W-and Z-boson resonances in the complex-mass scheme [33][34][35], which maintains gauge invariance and NLO precision both in resonant and non-resonant phase-space regions. For the evaluation of the one-loop integrals in the virtual corrections we have replaced the original internal integral library by the publicly available Fortran library Collier [48]. Ultraviolet (UV) divergences are treated in dimensional regularization, while the (soft and collinear) infrared (IR) divergences of the loop integrals and in the real photon or gluon emission are regularized by small photon, gluon, and external fermion masses. The final-state fermions are considered in the massless limit, i.e. small fermion masses are only kept as regulators in the singular logarithms. 1 However, in diagrams with a closed fermion loop the full mass dependence of those fermions is kept which allows to extend the calculation to include a heavy fourth fermion generation, as done in Ref. [50]. The cancellation of the IR divergences can be performed via phase-space slicing [51] or dipole subtraction [52][53][54].
The integration over the phase space is done using an adaptive multi-channel Monte Carlo integrator, where the integrand is evaluated at pseudo-random phase-space points, and the density of the points is adapted iteratively to the integrand to provide a better convergence. The Monte Carlo generator can also be used to generate samples of unweighted events for leptonic final states, which is particularly interesting for experimental analyses. Prophecy4f automatically provides distributions for leptonic and semi-leptonic final states. Distributions  (6) uūcc (1) e − e + µ − µ + (3) ν eνe dd (9) ddss (3) ν eνe µ − µ + (6) e − e + uū (6) uūss (4) e − e + dd (9) neutral current with interference e − e + e − e + (3) uūuū (2) ν eνe ν eνe (3) dddd (3) charged current ν e e + µ −ν µ (6) ν e e + dū (12) udsc (2) charged and neutral current ν e e + e −ν e (3) uddū (2) for fully hadronic final states are not predefined, since this should be done in the hadronic production environment. The h→4f decay width is the sum of all partial widths of the 19 independent final states listed in Tab. 1. All other final states differ only by generation indices and yield the same result, since the external fermion masses are neglected. One can weight these independent final states with their multiplicity (given in parentheses in Tab. 1) instead of computing partial widths for all existing final states. However, it is also of interest to separate the contributions from ZZ or WW intermediate states and WW/ZZ interferences in the partial width, as, e.g., described in Refs. [55,56], (2.1) The decomposition is trivial for 4f states to which only WW or ZZ intermediate states contribute; only one of the first two terms contributes in this case. Both WW and ZZ intermediate states can only contribute if all four final-state fermions are in the same generation (in the absence of quark mixing, which does not play a role in these processes). In such cases the WW and ZZ parts can be extracted by replacing the 4f state by ff F F and ff FF states with the same flavours as in the original ff f f state, but taking f and F from different generations. The interference term is then obtained by subtracting the WW and ZZ parts from the full squared matrix element. Exemplarily for the ν e e + e −ν e final state this reads

2)
With this procedure the contribution of all final states to the WW, ZZ partial widths, and the WW/ZZ interference contribution can be computed [55,56],

Details of the NLO calculation and implementation into Prophecy4f
We extend Prophecy4f to the calculation of the corresponding decays of the light CP-even Higgs boson h in THDMs. Specifically, we consider THDMs of Type I and II, as well as models of "lepton-specific" and "flipped" types. For the calculation of the decay h → WW/ZZ → 4f , we identify the decaying Higgs boson h with the discovered Higgs boson of mass 125 GeV. The calculation is similar to the one in the SM described in detail in Refs. [30,32]. As the particle content of the SM is extended in the THDM, all diagrams of the SM appear also in the THDM calculation. However, coupling factors of interactions involving scalar particles are modified in the THDM and have to be adapted. In addition, new diagrams including heavy Higgs bosons appear and need to be taken into account. In the following, we discuss the leading-order (LO) matrix elements and the EW and QCD NLO corrections.

Lowest order
At LO, the decay of the Higgs boson proceeds via a pair of (off-shell) gauge bosons V = W, Z which subsequently decay into fermions, as shown in Fig. 1. The diagrams involving a coupling of scalars to external fermions vanish and can be omitted due to the neglect of the external fermion masses. The only change in the THDM w.r.t. the SM is that the hV V coupling acquires an additional factor of sin (β − α), so that the LO matrix element becomes where α is the mixing angle between the two neutral CP-even Higgs bosons h and H, 2 and β is the mixing angle in both the neutral CP-odd as well as in the charged scalar sector which is related to the ratio of the vacuum expectation values of the two scalar doublets. We consistently follow the conventions of Ref. [40] for all quantities of the THDM.

Electroweak corrections
In the EW corrections, heavy Higgs bosons appear in loop diagrams, viz. in self-energy and vertex corrections. Generic diagrams are shown in Fig. 2. Four-and five-point diagrams do not contain heavy Higgs bosons, as this would require an hff coupling which is proportional to the mass m f of an external fermion f . The one-loop diagrams that do not include these heavy particles are in direct correspondence to the SM diagrams described in detail in Ref. [30]. However, the coupling factors of the internal (massive) fermions and the vector bosons to the light Higgs field need to be adapted in the THDM. The counterterm contribution can be split into two parts. The first one, M CT SM , is analogous to the counterterm contribution in the SM, although all renormalization constants appearing in this part are defined within the THDM using the renormalization conditions described in Ref. [40] and in general receive contributions from the exchange of heavy Higgs bosons. The second part is composed of the renormalization constants of the mixing angles α, β, entering via the factor sin (β − α) in M V V SM,LO , and the field renormalization constant of the mixing of the neutral CP-even fields. The full counterterm can be written as where we introduced the abbreviations s x ≡ sin x, c x ≡ cos x, t x ≡ tan x. Following Ref. [40], we employ four different renormalization schemes in order to define the mixing angles at NLO, i.e. to fix the renormalization constants δα, δβ: • MS(α) scheme: In this scheme α and β are independent parameters and fixed in the MS scheme. Tadpole parameters are renormalized in such a way that renormalized tadpole parameters vanish. As discussed in detail in Refs. [38,39], this treatment introduces gauge dependences in the relation between bare parameters and, thus, the relation between renormalized parameters and predicted observables inherit some gauge dependence. Since we work in the 't Hooft-Feynman gauge, all predictions (not only the ones presented in this work) should be made in the same gauge to obtain a meaningful confrontation of theory with data.
• MS(λ 3 ) scheme: This scheme coincides with the MS(α) scheme up to the point that α is traded for the (dimensionless) Higgs self-coupling parameter λ 3 as independent parameter. The coupling λ 3 is fixed by an MS condition, and α can be calculated from λ 3 and the other free parameters by tree-level relations. Renormalized tadpoles are again forced to vanish, however, in this scheme the relations between free parameters and predicted observables are gauge independent within the class of R ξ gauges at NLO, because λ 3 as basic coupling in the original Higgs potential does not introduce gauge dependences and the MS renormalization of β is known to be gauge independent in R ξ gauges at NLO [38,39].
• FJ(α) scheme: In this scheme, which is also described in Refs. [38,39] in slightly different technical realizations, again α and β are independent parameters, but gauge dependences are avoided by treating tadpole contributions differently, following a method proposed by Fleischer and Jegerlehner (FJ) [57] a long time ago already in the SM. 3 The basic idea is that bare tadpoles are defined to be zero, which preserves gauge independence in the relations between bare parameters of the theory. As a consequence, explicit tadpole loop contributions have to be taken into account in all loop calculations. This somewhat unpleasant feature can be avoided by introducing a new set of renormalization constants upon splitting off constant contributions from those fields that develop vacuum expectation values by field transformation in the functional integral (see Refs. [38,39]). Equivalently, the whole procedure of the MS(α) scheme, i.e. the full counterterm Lagrangian including tadpole counterterms, can be kept, but the renormalization constants δα, δβ, which contain only pure UV divergences in the MS(α) scheme, now receive some finite contributions from the different renormalization prescription of the tadpoles. This procedure is described in Ref. [40] in detail.
• FJ(λ 3 ) scheme: In this scheme β and λ 3 are independent parameters, as in the MS(λ 3 ) scheme, but tadpoles are treated following the gauge-independent FJ prescription.
More details on the different schemes and explicit results for the renormalization constants can be found in Ref. [40]. In all four schemes the parameters α, β, and the Higgs-quartic-coupling parameter λ 5 are defined directly in the MS scheme or are indirectly connected to MS parameters, i.e. all α, β, and λ 5 depend on a renormalization scale µ r . In Ref. [40] the µ r dependence of α, β, and λ 5 was taken into account by numerically solving the renormalization group equations in the four different renormalization schemes. Using these results on the running of α, β, and λ 5 , we will investigate the scale dependence of the NLO-corrected h→4f decay widths. In particular, we check whether the implicit µ r dependence of α and β, which already enters the LO amplitude, is compensated by the explicit µ r dependence contained in the loop corrections. The diagrams of the real emission can be obtained from the LO diagrams by adding photon radiation. The photon couplings in the THDM and in the SM are identical, i.e. the real emission matrix element M R,EW THDM of the THDM results from the matrix element M R,EW SM of the SM by multiplication with the coupling factor sin (β − α), The calculation of the SM amplitude M R,EW SM is described in detail in Ref. [30]. The IR-singular structure is not altered in the transition from the SM to the THDM, so that the subtraction and slicing procedures can be applied straightforwardly in the same way as it was done in the SM calculation [30]. Figure 3: Exemplary diagrams for the one-loop virtual QCD corrections. In the semi-leptonic case, only the first diagram type exists. Only the interference of the last four diagrams with the crossed LO diagram of Fig. 1(b) has a non-vanishing colour structure, demanding the quarks to be identical.

QCD corrections
As the THDM does not change the strongly interacting part of the theory, the computation of the QCD corrections is much simpler than for the EW corrections. Some diagrams of the virtual QCD corrections are shown in Fig. 3. In the diagrams (a), (b), (d), (e) the only coupling that changes w.r.t. the SM is the hV V coupling with an additional factor of s β−α . In the diagrams represented by Fig. 3(c), hqq couplings appear instead of the hV V , where q is any quark. The hqq couplings depend on the type of THDM and are given in Tab. 2. The QCD counterterm contribution is identical to the one appearing in the SM calculation up to the overall coupling factor s β−α in the matrix elements.
The diagrams for the real QCD corrections can be obtained from Fig. 1 by adding gluon emission off quarks and antiquarks. Similar to the EW case, the THDM does not affect the additional gluon emission, so that (2.11) The quark loop diagrams do not contain IR singularities, so that the singular structure of the one-loop matrix element matches the one of the corresponding SM amplitude multiplied by s β−α .
As the LO amplitude contains the same factor, the SM subtraction and slicing algorithms can be applied without modification.

Complex-mass scheme
To treat the vector-boson resonances in a proper way, we employ the complex-mass scheme which is explained in detail in Refs. [33][34][35]. This prescription consists of an analytic continuation of the masses of unstable particles into the complex plane which preserves gauge invariance as well as all algebraic relations between amplitudes or Green functions that do not involve complex conjugation (such as Ward and Slavnov Taylor identities). The complex mass µ V of V is directly connected to the real pole mass M V and the decay width Γ V , with V = W, Z. For our process at NLO, it is sufficient to treat only the W and Z boson in the complex-mass scheme even though the other scalar particles are not stable. We assume that in the THDM, the light Higgs boson has properties similar to the SM Higgs boson, i.e. its width is very small . Effects of this order can be neglected, as they are smaller than the contributions from NLO and have the same size as the uncertainties due to the separation of h production and decay. The other unstable Higgs bosons of the THDM enter only in loop diagrams, and the corrections from the complex masses are negligible as Γ S M S where Γ S and M S are the decay width and real pole mass of the considered Higgs boson. A fully consistent replacement of the real masses by its complex counterparts includes also a complex definition of the weak mixing angle θ W , The generalization of this prescription to the one-loop level leads to complex renormalization constants [34], for instance to the complex mass renormalization constants where Σ V V T (p 2 ) denotes the transverse parts of the V -boson self-energy with momentum transfer p and Σ V V T its derivative w.r.t. p 2 . As a consequence the renormalization constants of the weak mixing angle are The field renormalization constants of the vector bosons are given in Eq. (4.30) of Ref. [34]. In particular, they enter the calculation of the electric charge renormalization constant. The field renormalization constants of the stable fermions and the scalars are also affected by treating W and Z bosons in the complex-mass scheme as we do not take the real parts of the self-energies. Due to the appearing complex parameters, the self-energies and also the renormalization constants become complex. However, the field renormalization constants of internal fields drop out and those of external fields factorize from the LO, so that the imaginary parts drop out after squaring the matrix element at NLO.

Implementation and checks
The implementation of our calculation is performed in two independent ways: In the first method, we use the FeynArts [59] model generated as described in Ref. [40] and adapt it to the specific demands of the Prophecy4f calculation, so that masses of fermions belonging to closed loops are treated with the full mass dependence and the complex vector-boson masses are implemented. The amplitudes are generated and processed using FeynArts [59] and Form-Calc [60,61] and implemented into the Prophecy4f code. Additionally, the coupling factors in the Prophecy4f code of the HV V and Hff couplings are adapted to the THDM so that the LO, real photonic corrections, and the QCD corrections can be obtained by simple rescaling. In a second, independent calculation the amplitudes are generated via a FeynArts 1 [62] model file, and the counterterms are inserted by hand. These two implementations allow us to check the results and ensure their correctness. 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 Refs. [39,41] for the related Wh/Zh production channels (including W/Z decays) using crossing symmetry.

Input parameters and scenarios for the THDM
For the SM-like input parameters we take the values recommended by the LHC Higgs Cross Section Working Group [37] which essentially follow the Particle Data Group [63]: We employ the G µ scheme where the electromagnetic coupling is derived from the Fermi constant G µ , which absorbs the running of the electromagnetic coupling α em from the Thomson limit to the electroweak scale and accounts for universal corrections to the ρ-parameter. The CKM matrix is consistently taken as the unit matrix, since all quark-mixing effects drop out in flavour sums, since we work with massless light quarks without mixing to the third quark generation.
Prophecy4f performs its calculation in the complex-mass scheme and automatically converts the experimentally measured on-shell gauge boson masses M OS V to pole masses M pole V of the propagators according to From these measured input values, the program recalculates the widths of the vector bosons in O(α em ) in the SM using real mass parameters everywhere. This recalculation ensures that the branching ratios of the vector bosons are correctly normalized and add up to one for the SM. In the THDM, the heavy Higgs bosons enter the width in the mass counterterms, however, as we are close to the alignment limit (c β−α → 0) the effects are negligible. For an easier reproducibility of our results, we keep the SM values, which are also compatible with the measured W/Z widths. The final-state fermions are considered massless. The fermion masses are only inserted as regulator masses in soft and/or collinear divergent terms and in mass terms of closed fermion loops.
The strong coupling constant α s appears only in the QCD corrections (see Sec. 2.2.3), and we take its value at the Z-boson mass. As central renormalization scale µ 0 we use the average mass of all scalar degrees of freedom, This scale choice might seem surprising at first glance, since the light Higgs-boson mass is the centre-of-mass energy of our process. However, the loop diagrams including heavy scalar particles S = H, A 0 , H ± with mass M S introduce potentially large terms of ln (M 2 S /µ 2 ) in the amplitudes as long as the mixing angle β − α stays away from the alignment limit. Therefore we adapt the choice of the scale to the arithmetic mean of the Higgs-boson masses, and the scale variations performed in Sec. 4 confirm this choice. The input values of the additional parameters of the THDM depend on the investigated scenario and are given in the following. The scale-dependent input parameters c β−α , t β , λ 5 are defined at the central scale µ 0 by default.
To provide collinear-safe differential observables, a photon recombination is performed in the real corrections. This procedure invokes the addition of the photon momentum to the one of the fermion in the histogram if the invariant mass of a photon and a charged fermion is smaller than 5 GeV. When this is possible with more than one fermion, the photon is added to the fermion that yields the smallest invariant mass. We apply the photon recombination in all our calculations; further details about its impact are discussed in Ref. [30].
The recent report [37] of the LHC Higgs Cross Section Working Group summarizes a selection of benchmark scenarios proposed in other papers. We study the most relevant for our process. In particular, the scenarios proposed as BP1A in Ref. [13] are relevant for our work. For these benchmark scenarios experimental constraints from direct LHC searches, shown in Fig. 4, as well as theoretical constraints from vacuum stability and perturbative unitarity, illustrated in Fig. 5, are taken into account. Additionally, we employ perturbativity constraints to improve this scenario. Large coupling factors can lead to a breakdown of perturbation theory, so that we demand sufficiently small coupling factors. To this end, we compute the size of each coupling factor λ(S 1 S 2 S 3 S 4 ) of all the four-point Higgs-boson vertices at tree level, where S i = h, H, A, G, H ± , G ± for i = 1, . . . , 4, and use the largest coupling factor, λ/(4π) = max |λ(S 1 S 2 S 3 S 4 )|/(4π), as a measure. We use Mathematica and our FeynArts model files exploiting the hybrid basis (c.f. Ref. [13]). The parameters Z 4 , Z 5 , and Z 7 of the hybrid basis are related to our input parameters via with v 2 = 1/( √ 2G µ ). As the masses and mixing angles appear in the couplings, the perturbativity criterion gradually constrains the parameter space. Since the convergence of the perturbation series becomes worse with increasing coupling factors a clear discrimination of perturbative and nonperturbative parameter points is impossible. However, values of λ/(4π) larger than 1 indicate that higher-order corrections do not systematically become smaller and perturbativity is not given anymore which rules out such parameter points. Values between 0.5 and 1 usually still yield large higher-order corrections and need to be taken with care. The result of the perturbativity analysis is given in Fig. 6 for M H = 300 GeV (left) and M H = 600 GeV (right). Excluded areas are shown in gray, while blue (0.5 < λ/(4π) < 1) and yellow (0.3 < λ/(4π) < 0.5) indicate different sizes of the coupling factors. Parameter points where all couplings are smaller than 0.3 do not appear. The excluded trench at tan β = 1 is a singularity of the hybrid basis used in Ref. [13], since t 2β in Eq. (3.7) and, hence, the coupling factors diverge at this point. Overlaying these results with the previous experimental and theoretical constraints shows a significant reduction of the allowed parameter region. Nevertheless, we need to modify the scenarios proposed by Ref. [13] only slightly to obtain the low-and high-mass scenarios as well as the benchmark plane scenario described later. We do not consider heavy Higgs masses in the TeV range, because the allowed parameter space is dramatically reduced in this region, which can be seen in Fig. 7.   Only parameters close to the alignment limit and with t β ≈ 2 and t β ≈ 0.5 remain allowed for |c β−α | ∼ 0.1.

Low-mass scenario:
The low-mass scenario, which we already introduced in Ref. [40], consists of a heavy neutral CP-even Higgs boson H of mass M H = 300 GeV. The input values are based on a benchmark scenario of Ref. [13] and consist of a THDM of Type I with Specifically, scenario A contains a scan in c β−α , as this is the only parameter of the THDM appearing at LO. The range of the scan is limited by constraints from experiments and perturbative unitarity. These constraints indicate that values of |c β−α | exceeding 0.1 are phenomenologically disfavoured [13]. However, we perform our analysis with less stringent bounds to get a more complete picture. We take two distinguished points of the scan region named Aa and Ab with c β−α = ±0.1 to perform scale variations: High-mass scenario: The high-mass scenario is again based on a Type I THDM, however, with heavier Higgs bosons, 3. Different THDM types: In this scenario, we compare different types of THDMs. Yukawa couplings appear in our process only in closed fermion loops in Higgs-boson twopoint functions and in hV V and hgg vertex corrections, so that the top-quark contribution is dominant. The couplings to up-type quarks is identical in all types of THDM, so that we expect negligible effects from changing the type. The comparison is performed for scenarios Aa and B1a.

Benchmark plane:
For this scenario, we analyze a large area of the M H − tan β plane: The fixed parameters are based on the Type I non-alignment scenario of Ref. [13], and are given in the hybrid basis (c.f. Ref. [13]) by

Baryogenesis:
The BP3 B scenario of Ref. [37] was initially proposed in Ref. [64]. With a second Higgs doublet, a first-order electroweak phase transition is possible, which could explain the baryon asymmetry in the universe. The main signature for this model is the decay of a pseudoscalar Higgs boson. Nevertheless, the non-alignment of the benchmark points BP3 B renders these scenarios also interesting for our study. The parameterization in the original form uses m 12 for the Higgs self-coupling parameter from which we compute λ 5 using (3.14) The input parameters are 6. Fermiophobic heavy Higgs: By choosing a Type I THDM as well as a vanishing mixing angle α, the heavy Higgs boson H decouples from the fermions. Such a scenario was proposed in Ref. [10] with a direct detection of the heavy Higgs bosons as the leading signature. However, the alignment limit (c β−α = 0) cannot be reached in this model as this would require large values of tan β which are ruled out by stability constraints. This gives rise to possibly sizable effects on the light Higgs-boson decay. Different tan β values can be chosen, and with larger values the alignment limit is approached: We transform the input parameter m 12 to our convention using Eq. (3.14) and use the same fixed λ 5 for all tan β, and

Numerical results
In this section we present our numerical results for the decay h → 4f of the light CP-even Higgs boson h in the THDM for the different scenarios described in the previous section, beginning with the low-mass scenario. There, we investigate at first the conversion of the renormalized input parameters between different renormalization schemes and the running of the couplings. Afterwards we discuss the scale dependence of the h→4f width and show the dependence on c β−α . First results of this study have already been published in Ref. [40]. Finally, we study the partial widths and differential distributions in order to identify deviations from the SM expectations. The same procedure is performed for the high-mass scenario (split into two regions with positive or negative c β−α ), while we do not perform such a detailed analysis for the other scenarios.

Conversion of the input parameters
The numerical values of an input parameter defined via different renormalization conditions in two renormalization schemes do in general not coincide in one and the same physics scenario, but have to be properly converted from one scheme to the other. This means that the mixing angles α and β have to be converted in the transitions between the four renormalization schemes described in Sec. 2.2.2, as already discussed in Ref. [40]. For a generic parameter p, the renormalized values p (1) and p (2) in two different renormalization schemes 1 and 2 are connected via where p 0 is the corresponding bare parameter and δp (1,2) are the NLO renormalization constants of O(α em ). In case of more parameters, this is a set of coupled equations. For the conversion of p (2) to p (1) , we can either use the linearized solution where δp (1) (p (1) ) is approximated by δp (1) (p (2) ), or solve the set of implicit equations (4.1) numerically. The full solution of the implicit set of equations (4.1) has the advantage that converting parameters from one scheme to another and back is an identity, while this is only approximately the case in the linearized approach. The error of the linearized approximation is beyond our desired NLO accuracy as long as the perturbation series behaves well and higherorder terms are small. The comparison of the results obtained with the two methods allows for a consistency check of the computation. For the conversion of α and β, we employ the MS(α) scheme as one of the two schemes, so that we only have to deal with one set of finite counterterm contributions at a time.
In scenario A, we extend the range of c β−α values to −0.4 to 0.4, so that we get a larger picture even though the regions with large |c β−α | are ruled out by phenomenology. The results are shown in Fig. 8 with a conversion from (l.h.s.) and to the MS(α) scheme (r.h.s.), while the MS(λ 3 ) (green), FJ(α) (pink), and FJ(λ 3 ) (turquoise) schemes are employed as the second scheme. All other conversions can be seen as a combination of the presented ones. On the left-hand side, the gray dashed lines are the result obtained using the linearized approximation 4 . In both plots, we highlight the phenomenologically relevant region in the centre.
All curves show only small changes in the parameter values, and the numerical solution agrees well with the linearized conversion, affirming that the finite higher-order contributions of the counterterms are small, and perturbation theory is applicable. However, we would like to mention that a parameter set in the alignment limit (c β−α → 0) is not preserved in the conversion to other renormalization schemes. The alignment limit, thus, depends on the choice of the renormalization scheme. Outside the phenomenologically relevant region, the curves for the transformation involving the schemes with λ 3 as an independent parameter have a small region where large effects occur. This is an artifact as these schemes become singular at c 2α = 0.
In the vicinity of the singularity (corresponding here to c β−α ≈ −0.32), the MS renormalization of λ 3 introduces large finite contributions to the conversion equation resulting in a breakdown of perturbation theory 5 . This occurs in the MS(λ 3 ) and the FJ(λ 3 ) renormalization schemes which limits their use. For scenario A, the phenomenologically relevant region is, however, not affected by this artifact.

The running of c β−α
Parameters renormalized in MS depend on a renormalization scale µ r , where the dependence is governed by the renormalization group equations (RGEs) of the THDM (see, e.g., Refs. [65][66][67][68][69]). For each renormalization scheme we solve the RGEs using a classical Runge-Kutta algorithm. We isolate the effects of the running from the conversion by considering each renormalization scheme separately, but do not convert the input values in this investigation. The scale dependence of c β−α from µ r = 100 GeV to 900 GeV is plotted in Fig. 9  ing the parameters by gauge-independent ones introduces additional terms in the β-functions, which arrange for a stronger scale dependence of such schemes. We find similar results comparing the gauge-dependent MS schemes with the gauge-independent FJ schemes in Fig. 9. It is remarkable that the sign of the slope differs for the different renormalization schemes. This is another consequence of the additional terms in the β-functions and shows that the choice of the renormalization scheme has large effects. As some curves hit the c β−α = 0 axis and therefore run into the alignment limit, we explicitly see that this limit depends both on the renormalization scheme and on the renormalization scale in a given scheme. In Fig 9(b) one can also see that the curves for the MS(λ 3 ) and the FJ(λ 3 ) scheme terminate around 250 GeV. At this scale, the running of λ 3 yields unphysical values for the parameters of the theory. This is unique to the λ 3 running as only there an equation needs to be solved in the relation to c β−α . For the other cases we prevent the angles from running out of their domain of definition by solving the running for the tangent of the angles.

Scale variation of the width
Owing to the appearance of heavy Higgs bosons in the loop diagrams multiple scales occur in the calculation of the NLO EW corrections in the THDM. Therefore, a naive choice of the central renormalization scale of µ 0 = M h might not be appropriate. To choose and to justify our central scale of Eq. (3.4), and to estimate the theoretical uncertainties, we compute the total width according to Eq. (2.1) while the scale is varied by roughly a factor of two around µ 0 . As a definition of the input parameters in each of the four renormalization schemes represents a physical scenario on its own, we have four input prescriptions (MS(λ 3 ), MS(α), FJ(α), FJ(λ 3 )), and for each of them we compute the result in all renormalization schemes. After converting the input to the desired renormalization scheme, we evolve the MS parameters from the scale µ 0 to µ r by solving the RGEs, and finally compute the h→4f width. The results are shown in Figs. 10 and 11 at LO (dashed) and NLO EW (solid) for the scenarios Aa and Ab for each of the input prescriptions. The QCD corrections are not part of the EW scale variation and therefore omitted in these results.
The benchmark scenario Aa shows almost textbook-like behaviour, and the results are similar for all input prescriptions so that we discuss all of them simultaneously. First of all, the LO (d) Figure 11: As in Fig. 10, but for scenario Ab.
schemes is improved. This is expected since results obtained with different renormalization schemes should be equal up to higher-order terms, if the input parameters are properly converted. The relative renormalization scheme dependence at the central scale, expresses the dependence of the result on the renormalization scheme. It can be computed for a specific input prescription from the difference of the smallest and largest width of the four renormalization schemes, Γ h→4f min (µ 0 ) and Γ h→4f max (µ 0 ), normalized to their average. In Tab. 3, ∆ RS is given at LO and NLO for each of the input variants and confirms the reduction of the scheme dependence in the NLO calculation. In addition, the choice of the MS(α) scheme as an input scheme leads to the smallest dependence on the renormalization schemes in scenario Aa. This fits well to the observation perceived when the running was analyzed that the MS(α) scheme shows the smallest dependence on the renormalization scale, attesting a good absorption of further corrections into the NLO prediction. Table 3: The variation ∆ RS of the h→4f width in scenarios Aa and Ab at the central scale µ 0 using different renormalization schemes (with NLO parameter conversions). The columns correspond to the schemes in which the input parameters are defined. The technical uncertainty in brackets is calculated by exploiting the integration errors for the central values corresponding to the maximal and the minimal width.
The situation for benchmark scenario 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. 9(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 previous scenario, and even missing for some of the truncated curves. Nevertheless, the situation improves at NLO, and the relative renormalization scheme dependence reduces, as shown in Tab. 3. Also the central scale choice of µ 0 seems to be appropriate in contrast to a naive choice of M h .
For both benchmark scenarios, 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 that 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. However, this procedure must be performed for each benchmark scenario separately, which is beyond the scope of this work for a larger list of benchmark scenarios.

c β−α dependence
The dependence of the h→4f width on c β−α is one of the central results of our analysis, as the decay observables of the Higgs boson into four fermions in the THDM are most sensitive to this THDM parameter. The h→4f width in dependence of c β−α in scenario A is shown in Figs. 12(a)-(d) for the four different input prescriptions. The LO width with NLO conversion (dashed) and the full NLO EW+QCD total width (solid) are computed in the different renormalization schemes after the NLO input conversion, using the constant default scale µ 0 of Eq. (3.4). The SM values with a SM Higgs-boson mass of M h are illustrated in red. The results are similar for all input prescriptions so that we discuss them simultaneously. At LO they show the suppression w.r.t. to the SM with the factor s 2 β−α . The differences at LO between the renormalization schemes are due to the conversion of the input. A pure LO computation is identical for all renormalization schemes, since there is no conversion in a pure LO prediction. This pure LO prediction is represented in each plot by the LO curve for which no conversion to another scheme is performed. The suppression w.r.t. to 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. This is also confirmed by the relative renormalization scheme dependence shown in Fig. 13.
The relative corrections to the h→4f width, defined by which is shown in Fig. 15 at LO (dashed) and NLO (solid) for parameters defined in the MS(λ 3 ) scheme (other input definitions yield similar results). The SM exceeds the THDM widths at LO and NLO. At LO the shape of c 2 β−α can be seen, with modifications due to input conversion. This shape is slightly distorted at NLO by the 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, which 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 the parameter region with |c β−α | < 0.1 even less than 2%, which will be challenging for experiments to measure.
We also investigate the origin of the relative EW corrections. To this end, in Fig. 16 we plot different contributions to the full correction in the MS(λ 3 ) renormalization scheme (the breakup in the other schemes is qualitatively similar). The contribution called "SM-like+THDMvirt" comprises all diagrams that have a SM correspondence as well as the real corrections of Eq. (2.10), and the diagrams of Fig. 2 with at least one heavy Higgs boson in the loop and the contributions of heavy Higgs bosons to the counterterms. The part called "Hh-mix" is defined by the contributions of the Hh-mixing field renormalization constant δZ Hh and the renormalization constant δα in Eq. (2.9). Note that this splitting is neither gauge-independent, nor UV-finite. 7 However, the major contribution to the Hh-mix part is furnished by Higgs-boson loops, which do not depend on the gauge, so that some qualitative conclusions may be drawn. The SM-like+THDM-virt diagrams shows a small off-set from the SM result and only a small c β−α dependence, showing that the modification of the coupling factors in the THDM is small, but grows when the alignment limit is left. In the alignment limit, the off-set originates only from the heavy Higgs boson contribution, since only these diagrams introduce differences w.r.t. the SM. For values of |c β−α | > 0.05 the major deviation from the SM and the shape of the EW corrections are mainly due to the Hh Higgs mixing. These terms factorize from the LO contribution and thus lead to a uniform (i.e. phase-space independent) and universal (i.e. final-state independent) correction factor to the LO prediction.

Partial widths for individual four-fermion states
The partial h→4f widths (as defined in Sec. 2) at NLO are shown in the MS(λ 3 ) scheme for benchmark scenarios Aa and Ab in Tabs. 4 and 5, respectively. For other schemes, the numbers differ slightly, but show the same qualitative pattern, so that we do not show them here. In the tables, we do not only state the full NLO QCD+EW partial widths, but also the relative EW and QCD corrections. The qualitative picture is similar for the two benchmark scenarios. The WW contribution originating from charged-current final states yields the largest contribution, while the ZZ contribution is minor and the interference term yields a small negative contribution. The EW corrections to the WW-mediated final states are uniformly about 2−3%, which determine the EW corrections to the partial h→4f width. The EW corrections to the neutral-current final states strongly depend on the fermion flavour and range between ±3%. The QCD corrections are essentially the strong corrections to W/Z → qq and therefore amount to α s /π for each pair of quarks in the final state. The uūuū and dddd final states, where interference contributions between two different ZZ channels exist, are somewhat exceptional with QCD corrections of only about 4%. The deviations ∆ SM from the SM expectation are shown at NLO and LO in the last two columns. The LO deviation is due to the suppression factor s β−α of the hV V coupling w.r.t. the SM and therefore identical with s 2 β−α − 1 = −c 2 β−α = −10 −2 for all final     states, since c β−α = ±0.1. It should be noted that the indicated errors are integration errors, and the presented LO results are thus also a consistency check. At NLO the deviation is slightly larger, though still within only 1.3% (2%) for the Aa (Ab) benchmark scenario. The deviations from the SM are quite uniform, i.e. insensitive to the final state, so that they are described by the partial h→4f width well within a few per mille.

Differential distributions
The program Prophecy4f provides invariant-mass and angular differential distributions for the h→4f decays. The differential decay widths may serve as a window to observe BSM effects as the shape of distributions might be distorted significantly by new coupling structures. This might occur even if the partial widths do not change significantly, and therefore, the differential distributions of leptonic and semi-leptonic final states (see Sec. 2) are important observables.
In the following, we study them for both charged-and neutral-current processes, e.g. the fully leptonic final states e − e + µ − µ + (nc), ν e e + µ −ν µ (cc), and the semi-leptonic final states e − e + qq (nc), ν e e + dū (cc). Most likely, differential distributions for fully hadronic final states are not experimentally accessible. A detailed discussion of the SM distributions at NLO, including issues of final-state radiation (such as photon recombination), can be found for the fully leptonic final states in Refs. [30,31,36] and for semi-leptonic final states in Ref. [32]. In our study we emphasize the differences between the SM and the THDM results, while the features of photonic (and gluonic) corrections in the THDM and the SM are identical. The distributions discussed in the following are calculated using the MS(λ 3 ) renormalization scheme; the other renormalization schemes yield similar results.

Leptonic final states:
We begin with the leptonic final state e − e + µ − µ + which is a decay mediated by Z bosons. The invariant mass M faf b of a fermion-anti-fermion pair is defined by with the momentum k a of the fermion f a and k b of the anti-fermionf b , where the photon momentum is already added to the fermion momentum in case of recombination. The NLO invariant-mass distributions of the muon pair are displayed in the first panel of Fig. 17(a) for the SM and the THDM in scenarios Aa and Ab and show the Z-boson resonance. The relative corrections normalized to the LO are illustrated in the second panel and exhibit the well-known effects of final-state radiation near the Z resonance: Photons radiated off a final-state lepton lower the invariant mass of the lepton pair and lead to positive corrections-the "radiative tail"for invariant masses below the Z-boson peak and negative corrections above. These corrections would contain a logarithm of the form α ln(m µ /M h ) from collinear photon emission off muons if no photon recombination was applied. However, photon recombination mitigates this large effect by shifting events back to larger invariant masses for collinear emission and leads to the necessary level of inclusiveness required by the Kinoshita-Lee-Nauenberg theorem [70,71] to remove the collinear singularity. In case of photon recombination, the µ − µ + and e − e + invariantmass distributions are, thus, identical. Yet, non-collinear photons, which are not recombined, still lead to a sizable net effect which is observed in the relative corrections. The shapes of the invariant-mass distributions in SM and THDM are practically identical, i.e. the impact of new Lorentz structures in the NLO THDM diagrams is negligible. The relative difference between the SM and the THDM distributions is just given by the difference observed already in the integrated h→e − e + µ − µ + decay width, i.e. −1.15% for scenario Aa and −1.81% for Ab. We recall that those differences were traced back to the impact of Hh-mixing effects in Sec. 4.1.4, which are independent from the decay kinematics, and thus conclude that those mixing effects are the dominant higher-order effects visible in the distributions as well.   Figure 18: As in Fig. 17, but for the leptonic charged-current decay h → ν µ µ + e −ν e .
The differential decay width with respect to the angle φ between the µ − µ + and e − e + decay planes is defined by with the sign convention where k 1 , k 2 , and k 3 are the momenta of the muon, the anti-muon and the electron, respectively. The corresponding distribution is shown in the upper panel of Fig. 17(b). One observes a cos (2φ) pattern in the shape of the distribution, which can be used to set bounds on nonstandard couplings of the Higgs boson to the EW gauge bosons (see Refs. [22][23][24][25][26][27][28][29]). Note that the oscillation pattern in the distribution of a pseudo-scalar Higgs boson would have a different sign. We again observe that the SM shape is not distorted by THDM effects and that the difference between SM and THDM prediction just resembles the difference in the integrated widths.
We have also considered the invariant-mass and angular distributions of the e − e + e − e + final state (not shown), for which interference terms between different ZZ channels appear. There, the assignment of the lepton pairs to intermediate Z bosons is not unique; usually the electron and positron with an invariant mass closest to the Z-boson mass is combined to a pair. Again we find that the relative difference ∆ SM between THDM and the SM is practically constant over the phase space and given by its values for the integrated width.
For the W-boson-mediated ν e e + µ −ν µ final state, the respective distributions are shown in Fig. 18. The invariant-mass distribution of M νµµ is not experimentally accessible, but shown for theoretical interest. The plot shows the W resonance around M νµµ ≈ M W with the radiative corrections caused by photon radiation as discussed above. As already observed in the neutralcurrent final state, there is no significant shape distortion in the THDM w.r.t. the SM prediction. As the neutrinos cannot be detected, neither the Higgs nor the W boson can be fully reconstructed. However, projecting all lepton momenta into a fixed plane mimics the experimental situation at the LHC in the centre-of-mass frame of the Higgs boson in the plane transverse to the proton beams, where the sum of the neutrino momenta is measurable as missing momentum. We, thus, analyze the transverse angle φ µe,T between the two charged leptons [30], defined by cos φ µe,T = k µ,T · k e,T |k µ,T ||k e,T | , sgn(sin φ T ) = sgn{e z · (k µ,T × k e,T )}, (4.9) where k i,T are the parts of the full lepton momenta k i orthogonal to the fixed unit vector e z representing the beam direction of the Higgs-boson production process. The distribution in φ µe,T is shown in the first panel of Fig. 18(b). As expected, no shape distortion is seen w.r.t. the SM prediction, only the constant relative deviation which is identical to the deviation in the partial width of −1.04% (Aa) and −1.87% (Ab). Other fully leptonic final states show similar patterns, so that their distributions are not separately shown here.
Semi-leptonic final states: The invariant-mass distribution of the hadronic system of the neutral-current final state dde − e + is displayed in Fig. 19 (l.h.s), together with the corresponding NLO EW+QCD corrections. In case of gluon radiation, the invariant mass is built from the whole hadronic qqg system to obtain an IR-safe observable. The distribution and the corrections show similar characteristics to the ones of the corresponding leptonic final state: Photon radiation leads to a radiative tail, but SM and THDM distributions do not show any visible shape difference. Note that the effect of the photon radiation is less pronounced compared to the leptonic final state, as the quark charge factors are smaller than for leptons. There is no   Figure 20: As for Fig. 19, but for the charged-current semi-leptonic decay h → ν e e + dū.
radiative tail from gluon radiation, because all gluons are recombined with the quark pair, so that only a flat QCD correction remains [32]. In order to analyze angular distributions, we identify the quarks and antiquarks with jets for events without gluon radiation. In case of gluon radiation, we always combine the two QCD partons with the smallest invariant mass to a single jet, so that we again obtain an event with two jets. As the jets cannot be distinguished, any observable must be invariant under the permutation of the two jets. For this reason, the angle φ between the two Z-boson decay planes can only be reconstructed up to the sign of cos φ, so that we define [32] | cos φ| = The corresponding distribution, which is depicted on the r.h.s. of Fig. 19, looks rather different from the leptonic case, since | cos φ| instead of φ is used in the binning. Again, the major finding is the fact that the shape of the distribution does not change in the transition from the SM to the THDM. Only the flat offsets of −1.12% (Aa) and −1.76% (Ab) already encountered in the partial width are visible. The invariant-mass distribution of the hadronic system of the semi-leptonic W-boson-mediated final state ν e e + dū is pictured in Fig. 20 (l.h.s) and shows the same characteristics as the one of the neutral-current final state considered above: a moderate radiative tail from photon radiation, flat QCD corrections (not explicitly shown), and no shape difference between SM and THDM predictions. The distribution in the angle between the electron and the hadronically decaying W boson, φ eW , in the rest frame of the Higgs boson is shown in Fig. 20 (r.h.s). The electron is predominantly produced in the direction opposite to the W boson, and the EW corrections slightly distort the shape of the distribution. The difference between SM and THDM is described well by the deviation observed for the partial width of −1.05% for the Aa and −1.85% for the Ab scenario.
To summarize, the effects of the THDM on the shapes of distributions are negligible, only offsets in the normalization are observed. Thus, distributions for the Higgs decay into four massless fermions are not helpful in the search for deviations from the SM induced by effects of the THDM.

High-mass scenario B1
The high-mass scenario is divided into two branches which are valid for positive or negative c β−α and have different tan β. In this section we cover scenario B1 with positive values of c β−α , scenario B2 with negative values is discussed in the subsequent section. The perturbativity measure increases with rising M H , as can be seen in Fig. 6, restricting the range in c β−α and potentially affecting the stability of the results in the high-mass scenario in a negative way. Instead of relaxing the situation by moving closer to the alignment limit, where no problems with too large corrections are expected owing to decoupling, we delibarately keep this parameter point in order to investigate the robustness of the different renormalization schemes by checking the scale uncertainty in the various schemes and by studying the scheme dependence. Appendix A.1 supplements the discussion of this section by results with c β−α = 0.05, which are closer to the alignment limit and show better perturbative stability. The following discussion of the numerical results is structured in the same way as for the previous scenario, beginning with the conversion of the input parameters between different renormalization schemes.

Conversion of the input parameters
We compute the conversion between the input values in different renormalization schemes for c β−α = −0.1 to 0.3 and use MS(α) either as input or as target scheme. Using input parameters defined in the FJ(α) scheme leads to particularly large changes in the c β−α , indicating that the NLO terms are large and that the perturbative expansion converges poorly, but also for the FJ(λ 3 ) scheme sizable shifts are observed. Owing to these large effects, the linearization of the conversion equations suffers from large uncertainties, and a proper numerical solution is desirable and shown Fig. 21(a). Actually, the two conversions of Fig. 21 should be inverse to each other, and we perform this consistency check in (a) by plotting the curves of (b) mirrored at the diagonal with orange dotted lines. The conversions of the two parameters (α, β) from one scheme to another in fact is invertible if the implicit equations are solved. In Fig. 21(b), this invertibility is not fully respected, since we consider only a projection of the conversion to the c β−α line, suppressing the changes in β in the plot, i.e. we always take the input values from Eq. (3.11) in the start scenario of the conversion.

The running of c β−α
The running of c β−α in scenario B1a is investigated in Fig. 22 analogously to the low-mass scenario for each renormalization scheme independently, without any conversion. The scale dependence of c β−α (µ r ) with c β−α (µ 0 ) = 0.1 is computed from µ r = 300 GeV to 1500 GeV using a Runge-Kutta method. We indicated regions where perturbativity is not valid with dotted lines using a slightly different perturbativity measure than in Figs. 6 and 7. Perturbativity is considered to be intact unless the largest of the quartic coupling parameters λ k of the Higgs potential with k = 1, . . . , 5 becomes too large, λ max k /(4π) > 1. Compared to the previously used perturbativity measure, we found that with this measure a slightly larger part of the parameter space fulfills the perturbativity criterion.
In comparison to the low-mass scenario scenario ( Fig. 9(a)) the scale dependence in the FJ(λ 3 ) scheme increases. For scales above the central scale we obtain large values of c β−α for which predictions become unreliable. But also the FJ(α) scheme shows a remarkable behaviour as the alignment limit is approached for low as well as for high scales.

Scale variation of the width
We now turn to the calculation of the h→4f width where we perform a scale variation similarly to the previous section in order to investigate the perturbative stability of the results and the validity of the central scale choice for scenario B1a. The scale is varied from µ r = 300 GeV to 1000 GeV, and the results are shown in Fig. 23 with one plot for each input prescription. First, the input values are converted to the target scheme, and afterwards the scale is varied. In regions where perturbativity is not valid (λ max k /(4π) > 1) the NLO result is plotted with dotted lines. The results do not show such a clear and well behaved picture as for the low-mass scenario: • The MS(λ 3 ), Fig. 23(a), and the MS(α) input prescriptions, Fig. 23(b), yield similar results. In both cases, these schemes as target schemes show very good agreement, an extremum and a distinct plateau region in which the central scale fits perfectly. They only begin to deviate when perturbativity breaks down. The other renormalization schemes do not behave as nicely: The results of the FJ(α) scheme has a significant offset and drops dramatically for lower scales, until perturbativity breaks down at about 400 GeV. The FJ(λ 3 ) scheme suffers from the strong running and diverges at high scales as expected, while it shows relatively good (but not stable) agreement with the other schemes for lower scales.
• For input values defined in the FJ(α) scheme (Fig. 23(c)), the conversion transports the large NLO corrections to all other schemes, so that perturbativity is not given at all, and all curves disagree. Together with the behaviour of the FJ(α) scheme in the other plots, we conclude that the perturbative predictions using the FJ(α) scheme are not trustworthy for this benchmark scenario.
• The FJ(λ 3 ) input prescription (Fig. 23(d)) seems to yield the best agreement between the schemes, however, the conversion to other renormalization schemes results in particularly small values for c β−α and therefore corresponds in the other renormalization schemes to a scenario closer to the alignment limit.  and displayed using the colour code of Fig. 9. The breakdown of perturbativity (λ max k /(4π) > 1) is indicated by changing the NLO curve to dotted lines. reduction of the scale dependence, the development of plateau regions for all schemes, and an overlap of the results from the different schemes.
In the computation of the relative renormalization scheme dependence at the central scale ∆ RS only reliable renormalization schemes should be used. Therefore only the widths computed in the MS(α) and MS(λ 3 ) schemes enter this calculation for which the result is shown in Tab. 6. Omitting the FJ schemes, our estimate of the scheme dependence is, thus, just the difference of the two MS schemes, which is very small both at LO and NLO. Nevertheless a tendency towards a reduction of the scheme dependence is seen in the transition from LO to NLO.  Table 6: The variation ∆ RS of the h→4f width in scenario B1a at the central scale µ 0 using the reliable renormalization schemes MS(λ 3 ) and MS(α) (with NLO parameter conversions). The columns correspond to the schemes in which the input parameters are defined. Using parameters defined in the FJ(α) and FJ(λ 3 ) schemes, the results are unreliable and a computation of ∆ RS is not meaningful. The zeroes in brackets show that the integration errors are negligible.
values defined in one of the FJ schemes, this analysis cannot be performed as the results are unreliable.

c β−α dependence
The h→4f width of the Higgs decay as a function of positive c β−α is shown for all combinations of input prescriptions and renormalization schemes in Fig. 24 at the scale µ 0 . The results from all schemes agree very well in the alignment limits, where c β−α → 0. For |c β−α | > 0.05, differences in the results obtained with different schemes after conversion from a common parameter input scheme start to become significant. The patterns observed in the investigation of the scale dependence recur. The dashed lines represent the LO result with an NLO input conversion, while at pure LO, i.e. without conversion, all renormalization schemes deliver identical results. The respective curve is the one where no conversion is necessary. The well-known s 2 β−α pattern can be observed at LO while the conversion into the FJ(λ 3 ) scheme introduces large corrections leading to a breakdown (see Fig. 21(a)), so that this scheme is only applicable for very low values of c β−α . The NLO results away from the alignment limit are more complicated: • The MS(α) and the MS(λ 3 ) input prescriptions (Figs. 24(a),(b)) have similar characteristics which is due to the small shifts of the parameters in the conversion. The width in the MS(α) and the MS(λ 3 ) renormalization schemes agree very well, and the agreement improves from LO to NLO, as desired. The FJ(α) scheme (as target scheme) shows differences which can be explained by large higher-order terms shifting the input values towards the alignment limit (see Fig. 21(b)). Owing to the large corrections at NLO, sizeable corrections beyond NLO are expected in the FJ schemes as well, i.e. for reliable predictions for the B1 scenario in the FJ schemes the inclusion of leading corrections beyond NLO should be calculated and taken into account.
• Using input values defined in the FJ(α) scheme, Fig. 24(c), expresses this problem more clearly. The large corrections spread to the other renormalization schemes and affect perturbativity in a negative way. But also within the FJ(α) scheme the corrections are large and differ from the s 2 β−α shape seen for other input variants. This confirms the conclusion of the previous section that predictions obtained using the FJ(α) are not reliable for this scenario.
• The good agreement of the renormalization schemes in the FJ(λ 3 ) input prescription ( Fig. 24(d)) is based on the shift of the input values towards the alignment limit in the conversion. This shrinks the range effectively to 0 < c β−α < ∼ 0.05 for the other target schemes after the conversion and pushes the results together.
For the input defined in the MS(λ 3 ) scheme, the relative corrections separated in EW, QCD, and EW+QCD are shown in Fig. 25. Taking the input in the MS(α) scheme instead, the results look similar (not shown). The QCD corrections are similar for all schemes, because only the closed quark-loop diagrams in the hV V vertex corrections do not factorize from the SM LO amplitude with the coupling factor sin(β − α), but the impact of those diagrams is small. In contrast to the low-mass scenario, the EW corrections decrease with increasing c β−α , so that the deviations from the SM shown in Fig. 26 are larger than in the low-mass case, although they remain below 2% for c β−α < 0.1. The relative renormalization scheme dependence can only be applied using the MS(α) and MS(λ 3 ) schemes, as the results obtained using the two schemes involving FJ prescriptions are only reliable for very small c β−α . From Figs. 23(a),(b), one can see that the differences between the MS(α) and MS(λ 3 ) schemes decrease from LO to NLO, and the scheme dependence is reduced.

Partial widths for individual four-fermion states
The partial NLO widths, the relative corrections δ EW/QCD , and the deviations from the SM, ∆ LO/NLO SM , are shown in Tab. 7 for benchmark scenario B1a using the MS(λ 3 ) renormalization scheme. The MS(α) scheme yields similar results which differ only at the permille level, whereas the other schemes are not reliable at this benchmark scenario. The widths are slightly smaller than in the low-mass scenario (see Sec. 4.1.5) in spite of the identical values of c β−α . The negative deviation from the SM rises to almost 2%, however, no final state accounts for distinctively large THDM effects that could be exploited in experiments.
In addition, the differential distributions of scenario B1a, as defined in Sec. 2, do not change the shape w.r.t. to the SM significantly. They are shown in App. A.2, together with the SM ones and the ones from scenario B2b. As observed in the low-mass scenario, for each four-fermion final state the difference between the h→4f widths in the THDM and the SM resembles a constant shift in all distributions as well.

High-mass scenario B2
To complete the discussion of the high-mass scenario, we turn to negative values of c β−α for which the parameter space is strongly reduced by perturbativity, stability, and unitarity constraints, leaving only a small branch around tan β = 1.5 and leading to scenario B2. Being in the vicinity of excluded parameter sets potentially affects the conversion, the scale dependence, and the reliability of the full results. Hence, similar to scenario B1, scenario B2 is well suited to actually address possible problems in that respect. We discuss the results in the same manner as scenario B1 for positive c β−α above and present results for the less delicate case with c β−α = −0.05 in App. A.1.

Conversion of the input parameters
The conversion of the input parameter c β−α between different renormalization schemes is shown in Fig. 27 for an enlarged range with MS(α) either as input or target scheme. The conversion into the MS(α) scheme shows several ominous features ( Fig. 27(b)). First of all, divergences for the MS(λ 3 ) and FJ(λ 3 ) schemes occur at c β−α ≈ −0.19. We have seen such a divergence already in the low-mass scenario outside of the relevant region (c.f. Sec. 4.1.1), which is caused by the singularity at c 2α = 0. Since the ratio of the vacuum expectation values is lower in this scenario, the divergence moves towards the alignment limit and closer to the relevant region. It affects the conversion for c β−α < −0.15 in the MS(λ 3 ) and FJ(λ 3 ) schemes, so that such values should be taken with care. If experimental observations favour this region of parameter space, it becomes necessary to redefine the renormalization scheme and choose a different Higgs self-coupling parameter (e.g. λ 1 ) or a combination (e.g. λ 1 + λ 2 ) as independent parameter renormalized in MS.
The singularity then appears in other parameter regions and allows for predictions with c β−α < ∼ − 0.15. Not only the schemes involving λ 3 are problematic, but also the conversion from the FJ(α) scheme, as large shifts indicate problems with the perturbative expansion, analogous to scenario B1. Figure 27(a) shows the results for the "inverse" conversion from the MS(α) scheme, together with the inverse of (b) obtained graphically by mirroring the curves at the diagonal (dashed orange). Note that the linearized approximation for the conversion would involve large uncertainties here. Although the comparison of these curves projects the reduction of the conversion to one dimension (spanned by c β−α ), which is thus not exact, it gives a quick overview over the convergence of the numerical solution. As expected, the singularity in the relation between λ 3 and α reduces the domain of definition for the conversion in Fig. 27(a)

The running of c β−α
The running of c β−α (µ r ) with c β−α (µ 0 ) = −0.1 is computed from µ r = 300 GeV to 1500 GeV with a Runge-Kutta method for scenario B2b. The result is shown in Fig. 28 for each renormalization scheme without any conversion. The curves look similar to the ones of the low-mass scenario pictured in Fig. 9(b) (although the range of µ r is different) and we observe the same effects: the truncation of the MS(λ 3 ) and FJ(λ 3 ) schemes, but also the strong scale dependence of the FJ(α) scheme and the good stability of the MS(α) scheme.

Scale variation of the width
For the h→4f width we again perform a scale variation in order to estimate theoretical uncertainties and to motivate the central scale choice. The method is as described in Sec. 4.1.3, and the results are shown in Fig. 29. The FJ(λ 3 ) renormalization scheme is not included as target scheme here, since it is not possible to convert input values to it for c β−α = −0.1 (see Fig. 27(a)), however, it can be used when the input values are defined in it. The observations correspond in the most cases to the ones of scenario B1: • The first two plots using parameters defined in the MS(λ 3 ) (Fig. 29(a)) and MS(α) (Fig. 29(b)) schemes show, as in the previous scenarios, similar characteristics. The result obtained with the MS(α) renormalization scheme shows almost no scale dependence, and its value agrees with the extremum of the the MS(λ 3 ) renormalization scheme which lies at the central scale. However, through the truncation of the running a broad plateau region cannot be observed for the latter scheme with input in MS(α). The width in the FJ(α) scheme is consistent with the results of the MS(λ 3 ) and MS(α) schemes at the central scale, but shows an offset at the plateau and decreases for scales below µ 0 , as expected from the running of c β−α .
• The results using the FJ(α) input prescription (Fig. 29(c)) are not conclusive, since large corrections from the conversion spread to all other schemes.
• The scale variation of the FJ(λ 3 ) input prescription (Fig. 29(d)) corresponds again to an aligned scenario in the other renormalization schemes. Closer to the alignment, the renormalization scheme dependence decreases, which can also be seen from the separate scale variation of a more aligned scenario with c β−α = −0.05 given in App. A.1.
Generically, we obtain a somewhat better improvement compared to benchmark scenario B1a, which probably originates from smaller perturbativity measures (see Fig. 6). The central scale of Eq.   Fig. 29.

c β−α dependence
The dependence of the h→4f width on c β−α is shown for the different input prescriptions in the four panels of Fig. 30, for all renormalization schemes. Close to the alignment limit, −0.05 < ∼ c β−α < 0, the results from different renormalization schemes agree nicely. Away from this limit, however, the results deviate significantly, demanding some discussion: • The curves obtained using the MS(λ 3 ) and the FJ(λ 3 ) input prescriptions, Fig. 30(a) and Fig. 30(d), show the largest deviation from the s 2 β−α dependence of the LO width because of the large corrections inherited from the parameter conversions to the other schemes, which were observed in Fig. 27(b). Defining the input in the FJ(λ 3 ) scheme, the NLO width even slightly increases with smaller c β−α values, since the corresponding c β−α value in the MS(α) scheme stays close to the alignment limit. Owing to the large conversion effects, especially the predictions in the FJ(λ 3 ) scheme involve large uncertainties, which could only be reduced by systematically including the leading effects beyond NLO.
• Using input values defined in the MS(α) scheme yields the smooth curves of Fig. 30(b) which have the expected s 2 β−α shape. The relative renormalization scheme dependence reduces from LO to NLO, while the breakdown of the MS(λ 3 ) and FJ(λ 3 ) schemes is manifest, since values of c β−α smaller than ∼ −0.1 or ∼ −0.05 in the MS(α) scheme cannot be converted into the MS(λ 3 ) or FJ(λ 3 ) schemes, respectively (cf. Fig. 27(a)).
• The FJ(α) input prescription shows largest deviations from the SM as large NLO contributions spread to the other schemes through the conversion, shifting the values away from the alignment limit and increasing the deviations from the SM prediction.
However, all results show a significantly better agreement between the renormalization schemes at NLO than for LO for all regions including the problematic ones, suggesting that the perturbative expansion works for this scenario in the vicinity of our central scale µ 0 in spite of partially large NLO terms. As the MS(λ 3 ) scheme has a limited region of applicability, we show in Fig. 31 the relative corrections using the MS(α) scheme, which is reliable over the whole scan range. The QCD corrections are similar to the SM and renormalization scheme independent, while the EW corrections show the breakdown of the MS(λ 3 ) and FJ(λ 3 ) schemes. The difference between the FJ(α) and the MS(α) schemes is slightly larger than in the low-mass case, however, the sizes of the corrections are almost equal. This results in similar deviations from the SM as can be seen by comparing Fig. 32 with Fig. 15 (it should, however, be noted that a different input scheme has been used), so that it is difficult to distinguish these scenarios using the Higgs decay into four fermions.

Partial widths for individual four-fermion states
We give the partial widths in Tab. 9 for scenario B2b in the MS(α) scheme, as this scheme provides reliable results for c β−α = −0.1. All the partial widths are similar to the ones of the low-mass scenario Ab (Tab. 5) in size (note, however, different input schemes have been used). This observation applies to the EW and QCD corrections and to the differences to the SM predictions as well. Again, there is no final state particularly sensitive to the THDM contributions. The differential distributions analogous to Sec. 4.1.6 are shown together with

Different THDM types
In this section, we compare the h→4f decay widths of the Type I, Type II, lepton specific, and flipped THDMs for the two scenarios Aa and B1a using the MS(λ 3 ) renormalization scheme. We do not expect large differences in the results, because the considered THDM versions differ only in the Yukawa couplings of Higgs bosons to the leptons and to down-type quarks, which are not enhanced by large fermion masses. The by far largest contributions involving Yukawa couplings, however, result from diagrams with top-quark-Higgs couplings, which are identical in all four THDM versions for the h→4f processes with massless external fermions. The results are shown in Tab. 10 with the numerical errors in parentheses, confirming our expectation: The differences originating from the different THDM types in the NLO corrections are below a permille and not even significant over the integration error (although we employ large statistics with 190 million phase-space points). The h→4f decay observables are, thus, rather insensitive to the different types of THDMs, so that our predictions are universally valid for all types.

Benchmark plane
For the benchmark plane scenario defined in Sec. 3 we analyze only the relative deviation ∆ SM of the h→4f width with respect to the SM in the MS(λ 3 ) scheme. At LO, this is −0.01 as c β−α = 0.1 is kept constant. The NLO corrections differentiate this picture as they are dependent on both the heavy-Higgs-boson masses and tan β. We show the result for a wide range in the (M H , tan β) plane in Fig. 33 where the colour of the parameter point indicates the   deviation ∆ NLO SM and gray areas are excluded by perturbativity constraints (λ max k /(4π) > 1). We interpolate between the computed parameter points to obtain a smooth result, however, the original grid can be seen at the border between the computed area and the area excluded by non-perturbativity. The major deviation is between 0 and −5% and grows in magnitude with increasing tan β. For very large values of this parameter and close to the perturbative exclusion, values up to −8% occur. Very interesting is also the region with a small tan β, as very small enhancements with respect to the SM can be found around M H = 300 GeV (displayed in green). However, this region has a strong mass dependence because for large masses the negative corrections become −5%. We note that this effect is also visible using the MS(α) scheme and therefore not an artifact of the singularity of the MS(λ 3 ) scheme.

Baryogenesis
In this section we discuss the results for the benchmark sets BP3 [37,64], as defined in Sec. 3, which were proposed as a possible solution to the problem of baryogenesis. The results shown in Tab. 11 are computed in the MS(λ 3 ) scheme without considering the other schemes. In spite of the large distance to the alignment limit, the small heavy-Higgs-boson masses render both scenarios perturbatively stable with perturbativity measures of about 0.4. Already at tree level we observe a large negative deviation from the SM caused by the large values for c β−α suppressing the hV V coupling. These effects are enhanced at NLO for which we observe an increase of the negative deviation by 3 percentage points. This should be used in experiments measuring the Higgs decay into four fermions to put stronger bounds on these scenarios.

Fermiophobic heavy Higgs
The results for the fermiophobic heavy Higgs scenario [10,37], as defined in Sec. 3, are shown in Tab. 12 for the MS(λ 3 ) scheme. The three scenarios have a perturbativity measure of λ max k = 0.6 and differ by their value of tan β. Note that all those scenarios are close to the alignment limit, so that the SM width is almost reached at LO. The NLO corrections increase the deviation from the SM by about 1% for all scenarios.

Conclusions
We have investigated the decay processes h → WW/ZZ → 4f in the THDM, where we identify the light neutral CP-even Higgs boson h with the discovered Higgs boson of mass M h = 125 GeV. This signature contributed to the discovery of the Higgs boson and is important in the experimental investigation of the properties of the Higgs boson, such as the measurement of its couplings to other particles. The corresponding decay observables allow for precision tests of the SM and, thus, contribute to the search for any deviations from SM predictions. The calculation of strong and electroweak corrections in specific SM extensions, such as the one presented in this paper in the THDM, is an important theory input to successful data analyses. In our phenomenological discussion of numerical results, we have considered several THDM benchmark scenarios proposed by the LHC Higgs Cross Section Working Group. For the investigated scenarios, we generally observe that the THDM predictions for the h→4f width are bounded from above by the SM prediction and that the deviations from the SM typically increase at NLO, which might be used to improve exclusion limits in the THDM parameter space. The individual partial widths show similar deviations from the SM for all final states, but the shapes of differential distributions are not distorted by THDM contributions, so that the latter are not helpful to identify traces of the THDM. Moreover, we find that the h→4f widths do not discriminate between different types of THDMs (Types I and II, lepton-specific and flipped).
We employ different renormalization schemes to define the THDM (i.e. the precise physical meaning of its input parameters) at NLO. Specifically, we apply four different schemes which have in common that we use as many as possible input quantities that are directly accessible by experiment, such as the (on-shell) masses of all five Higgs bosons of the THDM. For the remaining three free parameters, which are Higgs mixing angles and Higgs self-couplings, we adopt MS prescriptions in four different variants. In detail, the MS(α) scheme defines the two mixing angles α and β of the CP-even and CP-odd Higgs bosons, respectively, as well as the quartic Higgs self-coupling parameter λ 5 in the MS scheme, and FJ(α) is a modified variant of this scheme with a different treatment of tadpole contributions in such a way that no gauge dependence between input parameters and predicted observables is introduced. Similarly, we define the two schemes MS(λ 3 ) and FJ(λ 3 ) in which we replace the angle α by another self-coupling parameter λ 3 as input. For a consistent comparison of results obtained in the different renormalization schemes, the MS-renormalized parameters have to be properly converted between the schemes. Depending on the scenario, we observe sizeable conversion effects on those parameters which can grow very large in scenarios close to the experimental exclusion limits or in parameter regions where perturbative stability deteriorates. These corrections, in particular, imply that the so-called alignment limit, in which one of the CP-even Higgs bosons of the THDM is SMlike, corresponds to different Higgs mixing angles in different renormalization schemes (even to different angles of a given renormalization scheme if the renormalization scale is changed). This shows that a proper definition of parameters at NLO is mandatory in future predictions and parameter fits in the THDM when precision is at stake.
While we observe a reduction of both the renormalization scheme and renormalization scale dependence of the h→4f width in the transition from LO to NLO as long as all Higgs-boson masses are moderate and the distance to the alignment limit is not too large, some renormalization schemes prove unreliable, i.e. prone to large corrections beyond NLO, for scenarios with heavy Higgs bosons or away from the alignment limit. Generically, the comparison of the different schemes reveals that the gauge-dependent MS(α) scheme shows a minimal scale dependence which reflects good perturbative stability. The MS(λ 3 ) scheme deviates only slightly from the former, yields reliable results and in addition is gauge independent at one loop in R ξ gauges. However, a singular region in the THDM parameter space exists in which the scheme is not defined. If this region is experimentally favoured, it is necessary to redefine the scheme by replacing λ 3 by another scalar self-interaction λ k =3 , so that the singularity is avoided. The gauge-independent FJ schemes partially suffer from large corrections and can only be applied for parameter points with sufficiently small coupling factors. Since the different schemes do not yield reliable results for all scenarios, self-consistency checks should be performed for every scenario when higher-order corrections are computed. In cases where NLO fails to be predictive, NLO calculations should be stabilized upon including the leading (renormalization-scheme-specific) corrections beyond NLO, a task that is, however, beyond the scope of this paper.
In more detail, in the low-mass scenario with a heavy CP-even Higgs boson of 300 GeV we obtain textbook-like results for the scale dependence, i.e. an improvement of the scale uncertainty and a reduction of differences between all four renormalization schemes at NLO, which indicates that perturbation theory works well. The deviation of the h→4f width from the SM prediction are, depending on the parameter set, between 0% and −6%, to which the NLO corrections contribute about 1−2%. In high-mass scenarios with heavy Higgs bosons with masses of about 600 GeV, the coupling factors are larger, resulting in less predictive results and larger differences between the renormalization schemes. Close to the alignment limit, the results of all four schemes are self-consistent and nicely agree, but away from it differences occur. While the MS(α) and the MS(λ 3 ) (in its domain of definition) schemes still yield trustworthy results, the FJ(α) and the FJ(λ 3 ) schemes suffer from large corrections, and their results should be taken with care. The deviations from the SM are similar to the low-mass case, and the NLO corrections similarly contribute 1−2% to the deviations. The other investigated scenarios support the described picture as they yield similar results.
The calculated NLO corrections to all h → WW/ZZ → 4f decays are integrated in a new version of the Monte Carlo program Prophecy4f, extending its applicability to the THDM. The new code can be obtained from the authors upon request and will be available from the public webpage 8 soon. (d) Figure 35: As in Fig. 29, but for scenario B2 with c β−α = −0.05.