Full next-to-leading-order calculations of Higgs boson decay rates in models with non-minimal scalar sectors

We present a complete set of decay rates of the Higgs boson with the mass of 125 GeV at the full next-to-leading order in a variety of extended Higgs models; i.e., a model with an additional real singlet scalar field, four types of two Higgs doublet models and the inert doublet model. All the one-loop contributions due to QCD and electroweak interactions as well as scalar interactions are taken into account, and the calculations are systematically performed. Branching ratios for all the decay modes are evaluated in these models, and patterns of deviations in each decay mode from the standard model predictions are comprehensively analyzed. We show how these models with extended Higgs sectors can be distinguished by using our calculation of the branching ratios and future precision measurements of the Higgs boson decays.


I. INTRODUCTION
After the discovery of a Higgs boson, the Standard Model (SM) has been completed in a sense that the existence of all the predicted particles was confirmed experimentally. In the SM, the minimal form with an isospin doublet scalar field is introduced as the Higgs sector. Although the discovered Higgs boson shows similar properties to that of the SM under the current experimental and theoretical uncertainties, the possibility that the Higgs sector takes a non-minimal form is not excluded at all, and its exploration is one of the central interests of current and future high-energy physics. If the Higgs sector is extended from the minimal form, it has a different structure which can be classified by the number of scalar fields, their representations, symmetries and strength of coupling constants. There should be strong connection between these properties of extended Higgs sectors and the physics behind the electroweak (EW) symmetry breaking. Furthermore, a non-minimal structure of Higgs sectors could solve the problems which cannot be explained in the SM such as neutrino masses, dark matter and baryon asymmetry of the Universe. Therefore, the Higgs sector is one of the most important probes of new physics beyond the SM.
The most important property of a non-minimal Higgs sector is the prediction of multiple scalar bosons. Thus, discovery of additional scalar bosons will be a clear evidence of extended Higgs sectors. At the LHC, direct searches for a new particle are being performed continuously. On the other hand, existence of such additional scalar fields generally affects the couplings of the SMlike Higgs boson to the particles in the SM by the effect of mixing and the quantum correction, yielding deviations from the predictions of the SM. Therefore, detecting such deviations by precision measurements is also a strong signature for models with extended Higgs sectors. Moreover, from the pattern of deviations in various Higgs boson couplings we can indirectly distinguish the shape of the Higgs sector and further determine new physics models [1].
Although the current data from these measurements are consistent with the SM predictions, the experimental and theoretical uncertainties are not small yet; e.g., about 20% for the Higgs boson couplings to weak bosons and typically 20-50% for the Yukawa couplings of the third generation at the 95% confidence level. Above experimental uncertainties can be much reduced at future colliders; e.g., at the High-Luminosity LHC (HL-LHC) [15,16], the International Linear Collider (ILC) [17][18][19], the Future Circular Collider (FCC) [20], the Circular Electron Positron Collider (CEPC) [21], the Compact LInear Collider (CLIC) [22] and so on. For example, at the ILC with the collision energy of 250 GeV and the luminosity of 2 ab −1 , some of the Higgs boson couplings are expected to be measured with O(1%) level or better [18].
In order to extract information on new physics from these precision measurements at future experiments, accurate calculations with higher-order corrections are required in models with various extended Higgs sectors. Radiative corrections to the SM-like Higgs boson vertices have been studied in various Higgs sectors such as, for example, a model with a real isospin singlet Higgs field (HSM) [23][24][25][26], two Higgs doublet models (THDMs) [27][28][29][30][31][32], the inert doublet model (IDM) [33,34], the Higgs triplet model [35,36] and the Georgi-Machacek model [37,38]. In order to see differences of the prediction among these models, it is quite important to calculate the renormalized Higgs boson vertices with a consistent and systematic way. Recently, we have published a numerical program H-COUP (version 1.0) [39] to compute a set of SM-like Higgs boson vertices at one-loop level in various extended Higgs models; i.e., the HSM, four types of THDMs and the IDM. Other numerical tools are also available to calculate Higgs boson decays with radiative corrections in models with extended Higgs sectors; e.g., Prophecy4f [40,41] and 2HDECAY [42,43].
In this paper, we present a complete set of the decay rates of the SM-like Higgs boson (h) including the h → W W * mode at the full next-to-leading order (NLO) in QCD and EW as well as scalar interactions in the HSM, four types of THDMs and the IDM 1 . Some important results have already been highlighted in our letter paper [44], in which the calculation of the partial decay rate of the h → W W * mode was not yet included. We then calculate the branching ratios of the SM-like Higgs boson at NLO in these models. One-loop calculations are consistently performed based on the on-shell renormalization scheme for EW parameters [45][46][47] and the modified minimal subtraction (MS) scheme for QCD corrections [48] in these models with the extended Higgs sectors.
We discuss the amount of the NLO corrections of the Higgs boson decay rates in each model with detailed descriptions of the computation. We show various correlations of the deviation in the branching ratios from the SM predictions under constraints of the perturbative unitarity [49], vacuum stability [50], conditions to avoid wrong vacua [51] and experimental constraints. Finally, we investigate the possibility to discriminate the extended Higgs sectors from the difference of the prediction among the models. This paper is organized as follows. In Sec. II, we introduce the HSM, the THDMs and the IDM. In Sec. III, we present analytic formulae for the decay rates of the Higgs boson at NLO. EW corrections in each decay mode are discussed in detail. In Sec. IV, we show numerical results of the total width, the branching ratios and correlations of the branching ratios. Conclusions are given in Sec. V. In Appendix, explicit formulae for the NLO calculations are presented.

II. MODELS WITH NON-MINIMAL HIGGS SECTORS
In this section, we define the HSM, the THDM and the IDM in order. Before moving on to the discussion on each extended Higgs sector, let us briefly explain constraints on a parameter space, as their basic notion are common to models with the extended Higgs sectors.
First of all, the size of dimensionless parameters in the potential can be constrained by imposing the perturbative unitarity bound which has originally been introduced in Ref. [49] to obtain the upper limit on the Higgs boson mass in the SM. Using the equivalence theorem [52], this bound requires that the magnitude of partial wave amplitudes for the elastic scatterings of 2 body to 2 body scalar boson processes, including the Nambu-Goldstone (NG) bosons, does not exceed a certain value. Each eigenvalue of the s-wave amplitude denoted as a i 0 is required to satisfy: where ξ = 1 [49] or 1/2 [53]. We here take ξ = 1/2. We note that each a i 0 only depends on the scalar quartic couplings as only scalar contact interactions contribute to the scattering process at the high-energy limit.
Next, the vacuum stability bound provides an independent constraint on scalar quartic couplings. It requires that the Higgs potential is bounded from below in any direction with large field values. This condition is schematically written by where V (4) represents quartic terms of the Higgs potential. Although in the SM this condition is trivially satisfied by taking the scalar quartic coupling to be positive, in models with non-minimal Higgs sectors it is given by a set of inequalities in terms of scalar quartic couplings [50].
Furthermore, in extended Higgs sectors, wrong local vacua can generally appear in addition to the true vacuum giving the correct value of the Fermi constant G F . Thus, we have to avoid parameter regions which realize the depth of such wrong vacua to be deeper than that of the true one 2 . The condition to avoid the wrong vacua can be written by combinations of dimensionful and 2 There is still a possibility of a meta-stable electroweak vacuum even if such situations are realized. By taking into account such a possibility a constraint from wrong vacua is more moderate. For example, see Ref. [54].
dimensionless parameters in the potential [51], so that it can provide an independent constraint from the above two constraints.
Apart from these theoretical constraints, we need to take into account bounds from experimental data. At the LEP/SLC experiments, various EW observables have been precisely measured such as the masses and widths of the weak gauge bosons. These precise measurements can be used to constrain new physics effects which can indirectly be entered into the self-energy diagrams for weak gauge bosons. Such indirect effect, so called oblique corrections, is conveniently parameterized by the S, T and U parameters introduced by Peskin and Takeuchi [55,56], which are expressed in terms of two point functions of the weak bosons. From the global fit of EW parameters [57], new physics effects on the S and T parameters under U = 0 are constrained by with the correlation factor of +0. 91  Flavor experiments also provide important constraints on a parameter space of extended Higgs models, particularly models with a multi-doublet structure. We will discuss these constraints in more detail in Sec. II B about THDMs. As mentioned in Introduction, additional scalars have been directly searched at the LHC [2][3][4][5][6][7][8][9][10][11], and some of the constraints are interpreted in THDMs. Moreover, the Higgs coupling measurements also constrain the mixing parameters in THDMs [12][13][14] by using the so-called κ framework [59]. We note that the κ framework is constructed by the leading order (LO) relation, and hence the interpretation of such constraints at higher-order level might not be straightforward.
The application of these constraints to each extended Higgs sector will be discussed in the following subsections.

A. Higgs singlet model
The Higgs sector of the HSM is composed of an isospin doublet scalar field Φ with the hypercharge Y = 1/2 and a real singlet field S with Y = 0. These scalar fields are parameterized where v is the vacuum expectation value (VEV) of the doublet filed which is related to the Fermi while v S is the VEV of the singlet field. Because the singlet field does not contribute to the EW symmetry breaking, the component fields G ± and G 0 in the doublet field correspond to the NG bosons which are absorbed into the weak bosons.
The most general Higgs potential is written as where all the parameters are real. By the reparameterization of the Higgs potential, we can take any value of v S without changing physical results [60]. Hence, we take v S = 0 throughout the paper.
In the HSM, there are two physical neutral Higgs bosons. Their mass eigenstates are defined as where α is the mixing angle, and we define the domain of α by −π/2 ≤ α ≤ π/2. Throughout the paper, we use the shorthand notation for the trigonometric function as s θ ≡ sin θ and c θ ≡ cos θ.
We identify h as the discovered Higgs boson with a mass of 125 GeV. After solving the tadpole conditions, the squared masses of these Higgs bosons are expressed as where M 2 ij (i, j = 1, 2) are the squared mass matrix elements in the basis of (s, φ). Each element is given by with M 2 = 2m 2 S . The seven parameters in the potential are then expressed by the following five input parameters and the two parameters m h and v that are fixed by experiments.
Let us briefly discuss the other relevant parts of the Lagrangian in the HSM. The kinetic term is given by where D µ is the covariant derivative for the Higgs doublet. Because the singlet field S does not have the gauge interaction, the additional Higgs boson H couples to weak bosons only through the φ component of H. Thus, the gauge-gauge-scalar type interactions are given as where g is the weak gauge coupling and g Z = g/ cos θ W with θ W being the weak mixing angle. The Yukawa interactions are written by the same form as those in the SM: whereΦ = iσ 2 Φ * , and we do not show the flavor indices here. In the above equation, Q L , L L , u R , d R and e R are respectively the left-handed quark doublet, lepton doublet, right-handed up-type quark singlet, down-type quark singlet and charged lepton singlet. The singlet field does not couple to fermions, so that the interaction terms for h and H are extracted as As it is seen in Eqs. (13) and (15), the SM-like Higgs boson h couplings are universally suppressed by the factor of c α as compared to the corresponding SM values.
As we already mentioned at the beginning of this section, the parameters in the potential can be constrained by the unitarity, the vacuum stability and the condition to avoid wrong vacua. For the unitarity bound, there are four independent eigenvalues given in Refs. [25,61]. In this paper, we use the expression for the eigenvalues given in Ref. [25], where the same notation of the potential parameters as that in this paper is applied. The necessary and sufficient condition to satisfy the vacuum stability is given by [62] For the condition to avoid these wrong vacua is found in Ref. [60,63,64]. We use the expression given in Ref. [25]. In the HSM, one-loop corrected two point functions for weak bosons are found in Ref. [65]. Imposing the bound from the S and T parameters, we can obtain the upper limit on m H depending on the value of c α . Constraints on the mass of the additional Higgs boson and the mixing angle from the LHC data have been studied in Refs. [66][67][68][69].

B. Two Higgs doublet model
The Higgs sector of the THDM is composed of two isospin doublet scalar fields Φ 1 and Φ 2 with Y = 1/2. These doublets are parameterized as where v 1 and v 2 are the VEVs of two doublets with v = v 2 1 + v 2 2 , and their ratio is expressed by Having two doublet fields with the same quantum charges causes dangerous flavor changing neutral currents (FCNCs) at tree level, because both doublets couple to each type of fermions. In order to avoid such FCNCs, we impose a discrete Z 2 symmetry, where two doublets transform as Φ 1 → +Φ 1 and Φ 2 → −Φ 2 . One can introduce the soft breaking term of the Z 2 symmetry in the potential, retaining the good property of the flavor sector. In the following, we discuss the THDM with the softly-broken Z 2 symmetry and the CP-conservation.
The most general Higgs potential is given by where the m 2 3 term softly breaks the Z 2 symmetry. The m 2 3 and λ 5 parameters are taken to be real as we consider the CP-conserving case. The scalar mass eigenstates can then be defined as follows: where H ± and A are the charged and CP-odd Higgs bosons, respectively, while H and h are the CP-even Higgs bosons. Similar to the HSM case, we identify h as the discovered Higgs boson with a mass of 125 GeV. We define the domain of α and β to be −π/2 ≤ α ≤ 0 and 0 ≤ β ≤ π/2, respectively, so that the viable range for β − α is expressed as 0 ≤ β − α ≤ π.
After solving two tadpole conditions for h 1 and h 2 , squared masses of the charged and CP-odd Higgs bosons are given by where M 2 = m 2 3 /(s β c β ) which describes the soft-breaking scale of the Z 2 symmetry. For the two CP-even Higgs bosons, the squared mass matrix elements M 2 ij in the Higgs basis [70] are given by with λ 345 ≡ λ 3 + λ 4 + λ 5 . The squared masses of the two CP-even Higgs bosons and the mixing angle β − α are expressed in terms of the matrix elements M 2 ij as The eight parameters in the potential are then expressed by the following six input parameters and the two parameters m h and v that are fixed by experiments. In addition to these parameters, there is a degree of freedom of the sign of c β−α .
Let us discuss the other relevant parts of the Lagrangian. The kinetic term for the Higgs doublets are written as In the mass eigenbasis of the Higgs bosons, the gauge-gauge-scalar type interaction terms are extracted as The most general Yukawa interactions under the Z 2 symmetry are written as where the subscripts i, j and k are 1 or 2. These indices are fixed when we determine the Z 2 charge for right-handed fermions. As in Table I, there are four independent types of Yukawa interactions depending on the assignment of the Z 2 charge [71][72][73]. The interaction terms for the physical Higgs bosons are then extracted as Type-X (lepton specific) with I f = 1/2 (−1/2) for f = u (d, e) and V ud is the Cabibbo-Kobayashi-Maskawa matrix element.
Similar to the HSM, parameters in the THDMs can be constrained by both the theoretical and experimental constraints. For the unitarity bound, there are 12 independent eigenvalues of the s-wave amplitude matrix [74][75][76][77]. We use the expression for the eigenvalues given in Ref. [31].
The vacuum stability bound is sufficiently and necessarily satisfied by imposing the following conditions [50,[78][79][80] In addition, the wrong vacua can be avoided by taking M 2 ≥ 0 [54]. We thus only take the positive value of M 2 in the following discussion. The expressions of the two point functions for the weak bosons in the THDM are found in Refs. [81][82][83][84][85]. Constraints on the parameters in THDMs from the LHC data have been discussed in Refs. [68,69,[86][87][88][89][90].
Differently from the HSM, constraints from flavor experiments are important to be taken into account in the THDM. These bounds particularly provide the lower limit on the mass of the charged Higgs boson m H ± depending on the type of Yukawa interaction and tan β. For example, from the B s → X s γ data, m H ± hs to be greater than about 600 GeV at 95% confidence level in the Type-II and Type-Y THDMs with tan β 2, while O(100) GeV of m H ± is allowed in the Type-I and Type-X THDMs with tan β 2 [91]. Constraints on m H ± and tan β from various flavor observables are also shown in Ref. [92] in four types of the THDMs.

C. Inert doublet model
The contents of the scalar sector in the IDM are the same as those in the THDM. In the THDM, the Z 2 symmetry can be softly-broken by introducing the m 2 3 term in the potential, while in the IDM it is assumed to be unbroken even after the EW symmetry breaking. Thus, the potential is obtained from Eq. (18) with m 2 3 = 0. In addition, the second Higgs doublet Φ 2 is supposed not to develop the nonzero VEV, otherwise the Z 2 symmetry is spontaneously broken.
In the IDM, the component scalar fields of Φ 1 and Φ 2 do not mix with each other due to the unbroken Z 2 symmetry. Therefore, we can identify these scalar bosons (w ± 2 , z 2 , h 2 , h 1 ) defined in Eq. (17) with the mass eigenstates (H ± , A, H, h). The lightest neutral inert scalar boson can be a candidate of dark matter as it cannot decay into a pair of SM particles.
The mass formulae for the scalar bosons are changed from those of the THDMs, not just because of the absence of the m 2 3 term, but also the absence of the tadpole condition for h 2 . These are given as follows: where M 2 = m 2 2 . We then choose the following five parameters to be free input parameters of the IDM and the two fixed parameters m h and v.
The same conditions for the perturbative unitarity and vacuum stability in the THDM can be applied to the IDM, because these bounds are given in terms of the scalar quartic couplings. The condition to guarantee the inert vacuum with ( Φ 0 1 , Φ 0 2 = (v/ √ 2, 0)) is given by [93], Since the tadpole condition makes m 2 1 negative, and the vacuum stability condition constraints λ 1 and λ 2 to be positive, the condition given in Eq. (36) is satisfied by taking M 2 > 0. We refer to this condition as the one to avoid wrong vacua, according to the other two models discussed above.
For the constraints of the S and T parameters, we can use the same expressions as those in the THDM with s β−α = 1.
In the IDM, constraints on the masses of the additional Higgs bosons from collider experiments are relatively weak since the additional scalars do not couple to SM fermions. The constraints from the LEP and the LHC have been studied in Refs. [94,95] and Refs. [96,97], respectively. Dark matter constraints from relic density and direct detection also limit the parameter space; see, e.g., Refs. [97,98] for details.

III. DECAY RATES OF THE SM-LIKE HIGGS BOSON AT ONE LOOP
In this section, we discuss the decay rates of the SM-like Higgs boson; i.e., h → ff , h → ZZ * → Zff and h → W W * → W ff at NLO in EW and QCD. The loop induced decay rates h → γγ, h → Zγ and h → gg are also discussed at NLO in QCD.
We outline our one-loop calculations. For the computation of the EW corrections, we adopt the modified on-shell renormalization scheme which has been defined in Ref. [47], while for the QCD corrections we apply the MS scheme. In the on-shell scheme, all the counterterms appearing in the decay rates of the h → ff and h → V V * → V ff modes are determined in terms of the one particle irreducible (1PI) diagrams for one-and two-point functions of Higgs bosons, gauge bosons and fermions by imposing a set of the renormalization conditions. Adding these counterterms, one can obtain the ultra-violet (UV) finite one-loop corrected vertices.
The on-shell scheme is a physically quite natural renormalization scheme, and it is suitable to apply to EW observables as they include well defined scales such as the weak boson masses.
However, it has been known that the on-shell scheme introduces gauge dependent counterterms, particularly in some mixing parameters [99]. In extended Higgs sectors, a mixing between Higgs bosons can generally appear. We thus apply the so-called pinch technique to remove the gauge dependence in the renormalized vertex functions to our computations [23,32,47].
Apart from the UV divergences, there appear infrared (IR) divergences when we calculate virtual photon loop contributions. Such IR divergences can exactly be cancelled by adding contributions of real photon emissions, where the finite QED corrections are common to those in the SM. The analytic expressions of QED corrections are known for h → ff [100][101][102] and h → Zff [46]. Thus, we simply switch off the photon-loop contributions, and use these analytic formulae in our computation as the QED correction part. For h → W ff , on the other hand, we cannot separate EW corrections into QED and weak corrections. Therefore, by using the phase-space slicing method [103], we numerically evaluate both the virtual and real corrections to h → W ff , see Appendix D for details. For QCD corrections, we use the similar technique to remove IR divergences coming from virtual gluon loop contributions.
In our renormalization calculation, we choose the fine structure constant α em , the Fermi constant G F and the Z boson mass m Z as the input parameters for the EW parameters. In this scheme, all the other EW parameters such as v, m W and s W are outputs. Using the on-shell definition of the weak mixing angle; i.e., s 2 and the modified relation among the EW parameters: we can calculate the renormalized squared W boson mass as In Eq. (37), ∆r is calculated by [104] whereΠ W W is the renormalized W boson two-point function and the second term corresponds to vertex corrections and box diagram contributions to the muon decay rate. Including these three EW parameters, we choose the following parameters as the SM inputs: where ∆α em is the shift of the fine structure constant α em from zero energy to m Z . We also input the parameters in the potential given in Eqs. (11), (25) and (35) in the HSM, THDMs and IDM, respectively. We note that the parameters c α and s β−α in the HSM and THDMs, respectively, do not physically mean the mixing angles for the CP-even Higgs bosons after the renormalization.

A. Renormalized vertices
Important ingredients for calculations of decay rates of the Higgs boson are renormalized Higgs boson vertices. In our computations, the hff , hV µ V ν (V = W, Z) and hV µ V ν (VV = γγ, Zγ, gg) vertices are relevant, where the hV µ V ν vertices are one-loop induced. Each of these vertices can be decomposed into several form factors depending on their Lorentz structure as written below. The H-COUP program (ver. 1.0) [39] provides numerical values of these renormalized form factors in the extended Higgs sectors. 3 We fully use H-COUP in our numerical evaluation of the form factors.
The renormalized hff vertices can be decomposed into the following form factors: where p µ 1 (p µ 2 ) is the incoming four-momentum of the fermion (anti-fermion), and q µ (= p µ 1 + p µ 2 ) is the outgoing four-momentum of the Higgs boson. For the case with on-shell fermions; i.e., p 2 1 = p 2 2 = m 2 f , the following relations hold: These relations are used for the calculation of the Higgs boson decay into fermions discussed in Sec. III B.
Next, the renormalized hV µ V ν vertices are defined in terms of three renormalized form factors: where p µ 1 and p µ 2 are incoming four-momenta of the weak bosons, and q µ is the outgoing fourmomentum of the Higgs boson. Similarly, we can define the loop induced vertices aŝ For the on-shell photon and gluon with a four-momentum p µ i , the Ward identity holds, i.e, p iµΓ µν hVV = 0. This gives the following relation This relation can be applied to the computation of the loop induced decays of the Higgs boson.
For an off-shell photon appearing in the h → Zγ * → Zff decay mode at NLO, Eq. (45) cannot be used, so thatΓ 1 hVV andΓ 2 hVV separately appear, as it will be discussed in Sec. III C. We note that the form factorΓ 3 hVV is non-zero only when the Higgs boson is a CP-mixed state. In this paper, we consider the case with CP-conservation in the Higgs sector, so that this form factor becomes zero.
All the renormalized form factors for the hXX vertices defined above are further decomposed into the tree level and one-loop level parts as follows: HSM THDMs IDM Scaling factors for the Higgs boson couplings to fermions (κ f ) and weak bosons (κ V ) in the extended Higgs models at tree level. The ζ f factor in the THDMs is given in Table I.
The tree level contribution to each form factor, denoted as Γ i,tree hXX , is given as where the scaling factors κ f and κ V are given in Table II The first and second terms of the right-hand side are the contribution from 1PI diagrams and counterterms, respectively. As we mentioned at the beginning of this section, counterterms are determined by a set of on-shell conditions by which these are written in terms of 1PI diagrams for one-and two-point functions with a fixed value of squared momenta. Analytic expressions for the contributions from 1PI diagrams and counterterms to these renormalized Higgs boson vertices are presented in Ref. [24] for the HSM, Refs. [31,106] for the THDMs, and Refs. [34,106] for the IDM.
For the computation of the partial decay rates of h → V V * → V ff , we also need to calculate the one-loop corrected V ff vertices and box diagrams in addition to the above Higgs boson vertices.
In the massless limit for the external fermions, the renormalized V ff vertices are the same as those in the SM, while the contribution from the box diagrams is simply given by the SM expression multiplied by the scaling factor κ V . For the completeness, we present the analytic expressions for one-loop corrections to the V ff vertices in Appendix B and those for the box corrections in At NLO, the partial decay rate of the h → ff (f = t) process can be written in terms of the EW correction part ∆ f EW and the QCD correction part ∆ f QCD as where Γ 0 is the decay rate at LO expressed as with N f c being the color factor, i.e., N f c = 3 (1) for f to be quarks (leptons). The expression for the tree level form factor Γ S,tree hf f is given in Eq. (47). The EW corrections ∆ f EW can be further decomposed into weak corrections ∆ f weak and QED corrections ∆ f QED as: Here, the weak correction means contributions from W , Z and scalar boson loops, namely all the loop contributions except for the photon and gluon loops. We also use this terminology in the later discussion. The expression for ∆ f weak is given in terms of the form factors defined in Eqs. (41) and (46) as where ∆r is given in Eq. (39).
The QED (QCD) corrections are obtained by taking into account contributions from the virtual photon (gluon) loops and those from the real photon (gluon) emissions. We can obtain simple expressions for these corrections by neglecting the term proportional to m 2 f /m 2 h . 4 For the leptonic decays, f = , the QED correction in the on-shell scheme is given by [100][101][102] For the hadronic decays, f = q, the QED and QCD corrections in the MS scheme [107] with the renormalization scale µ are given by where C F = 4/3. In the numerical evaluation, we set µ = m h , and replace the quark mass in the tree level form factor Γ S,tree hf f in Eq. (50) by the running massm q (µ = m h ). From Eqs. (53) and (54), we can see that there is no additional Higgs boson mass dependence in their expression, so that these corrections do not provide deviations in the Higgs boson decay rate from the SM prediction at NLO.

C. h → ZZ * → Zff
We calculate the partial decay rate of the Higgs boson into a pair of weak bosons at NLO in this and next subsections. Because the mass of the discovered Higgs boson is about 125 GeV, at least one of the weak bosons must be off-shell. We thus calculate the process with 3-body final states, i.e., h → V V * → V ff . Throughout this paper, we neglect the masses of external fermions in the h → V V * → V ff processes. In Fig. 1, all the diagrams contributing to the process are shown.
Similar to the h → ff mode, the decay rate for the h → ZZ * → Zff mode at NLO is expressed and the EW correction can separately be expressed by the weak corrections and the QED corrections: The LO contribution to the decay rate of h → Zff is expressed by where s is the Mandelstam variable defined by (p µ f + p μ f ) 2 , and another variable u defined by (p µ Z + p μ f ) 2 is already integrated out in the squared tree level amplitude |M Z 0 | 2 expressed as The tree level form factor Γ 1,tree hZZ is given in Eq. (47), and that for the Zff vertex v f and a f is given in Eq. (B3).
The QED (∆ Z QED ) and QCD (∆ Z QCD ) corrections only enter the Zff vertex correction depicted in the diagram (c) in Fig. 1. Their expressions are common to the SM given as [46] Although the diagrams (d) and (e) can also receive QED and QCD corrections, they vanish in the massless limit for the external fermions.
The weak corrections ∆ Z weak are given by whereλ The similar expression in the SM is found in Ref. [46]. The first and second lines correspond to the contribution from the diagram (a) in Fig. 1, whereΓ 1,2 hZγ are the renormalized form factors for the loop induced hZγ vertex. The analytic expressions forΓ 1,2 hZγ are presented in Appendix A. The third line corresponds to the contribution from the diagrams (b) and (c). The diagram (c) contains the V ff vertex corrections, so that we need to prepare the calculation of the renormalized V ff vertex denoted asΓ V f f which will be implemented in the H-COUP ver. 2.0 [105]. In the massless limit of the external fermions, this correction becomes the same as the SM prediction. Details of the calculation ofΓ V f f are given in Appendix B. In the fourth line, the T Z hf f and B Z terms represent the contribution from the hff vertex corrections and the box diagrams shown as the diagrams (d) and (e) in Fig. 1, respectively. Both T Z hf f and B Z depend on the Mandelstam variable u in loop functions, which has to be integrated out. The integration range of u is given by The explicit formulae for T Z hf f and B Z are given in Appendix C. Although the hff vertex corrections can be calculated by using H-COUP ver. 1.0 [39], we present the explicit analytic formulae for the contribution from the diagram (d) in Eq. (C1) in the massless limit of the external fermions.
In this limit, both contributions from (d) and (e) become the SM predictions times the scaling factor κ V . We note that we need to add the contribution from the wave function renormalization of the external Z boson, i.e.,Π ZZ , because in our on-shell scheme, the counterterm for the wave function renormalization (normally denoted as δZ Z ) is not fixed by the condition which requires a vanishing derivative of the renormalized Z boson two-point function, but it is determined by the other conditions, see Ref. [47].
We compute the partial decay rate of h → W W * → W ff mode at NLO. Feynman diagrams are shown in Fig. 1. The decay rate is expressed as where Γ 0 , ∆ W EW and ∆ W QCD are the contributions from LO, EW corrections and QCD corrections, respectively. The expression for the LO contribution is obtained from Eqs. (57) and (58)   In addition to the h → ff and h → V V * → V ff decays, the Higgs boson can also decay into γγ, Zγ and gg. The LO contributions to the decay rates arise from one-loop diagrams, and they can be expressed in terms of the renormalized hVV (VV = γγ, Zγ, gg) vertices defined in Sec. III A as where Eq. (45) is used. The analytic expressions forΓ 1,2 hVV are given in Appendix A. Let us discuss QCD corrections to these loop induced decay rates at NLO. For h → gg, there are two sources of the QCD corrections: virtual gluon exchanges in the quark loop diagrams and real gluon emissions. In the MS scheme, the QCD corrected decay rate is given as [109]     In order to extract the new physics effects of the EW corrections to the partial decay rates, we where ∆ X EW | NP (∆ X EW | SM ) denotes the prediction of ∆ X EW in the models with the extended Higgs sectors (SM). Our results for ∆ X EW | SM are summarized in Table III. It is important to mention here that the dominant contribution to ∆ X EW comes from the first term of Eqs. (52), (60) and (64) where c ϕ = 2(1) for additional charged (neutral) scalar loop contributions. 5 The above equation indicates that scalar loop effects become non-decoupling when m ϕ is mostly determined by v, or equivalently the case with M 2 /v 2 1. In such a non-decoupling case, the right-hand side of Eq. (69) is nearly proportional to m 2 ϕ . Of course, there must be an upper limit on m ϕ from the unitarity bound, under which the magnitude of the non-decoupling effect can be typically a few percent level, as we will see in the plots below. 5 In the THDMs, the charged Higgs boson and top quark loop contribution can also be important for ∆ b EW . The analytic expression for the top quark mass dependence due to charged Higgs boson loops is found in Ref. [31]. EW turns out to be important to understand the behavior of the deviation in branching ratios from the SM prediction, which will be discussed in the next section. In the following discussion, we impose the bounds from the perturbative unitarity, the vacuum stability, the conditions to avoid wrong vacua (taking M 2 ≥ 0) and the S, T parameters discussed in Sec. II.
The flavor constraints discussed in Sec. II B are also important to be taken into account particularly in the THDMs, but we do not impose them here in order to study and compare the behavior of ∆ X EW among the extended Higgs sectors. In the later discussion given in Sec. IV C, we discuss the branching ratios imposing the flavor constraints as well.
In Fig. 2 On the other hand, in the case with s β−α = 0.99 shown in the lower panels, the decoupling limit cannot be taken, so that there appears the upper limit on the mass of the extra Higgs boson from the theoretical constraints depending on the value of tan β and the sign of c β−α . At m Φ ∼ 2m t , the threshold effects of tt appear, which push ∆ b EW into the positive direction. This peak comes from the top quark loop contribution to the Z-A mixing diagram which appears in the counterterm of the β parameter. More detailed discussions have been given in Ref. [30]. We can also see the dip above the tt threshold for tan β = 1.5. The origin of this dip can be explained in the same way as in Fig. 2 and tan β = 1.5 for all the types of the THDMs. It is seen that the direction of the peak at around m Φ = 2m t is determined to be positive (negative) if ζ f = cot β (− tan β), see Table I. The behavior of ∆ c EW and ∆ τ EW is also classified by the factor of ζ f , e.g. that of ∆ τ EW in the Type-II THDM is common to the Type-X THDM. Concerning ∆ b EW , the behavior is different from e.g. ∆ τ EW in the Type-II THDM even though these two depend on the same factor of ζ f . This can be understood by the fact that the top mass dependence enters in ∆ b EW when charged Higgs bosons run in the loop, while for the other ∆ f EW a dependence on fermion masses in loops is negligibly small. Detailed discussions have been given in Ref. [30] for the behavior of the radiative correction to the Yukawa couplings in the THDMs.
In Fig. 6, we show the value of ∆ Z EW in the Type-I THDM with a fixed value of tan β. The GeV for s β−α = 1 with tan β = 1.5, 3 and 5, respectively, is the same as that shown in the upper panel of Fig. 3, because it is determined by the unitarity bound. It is seen that for s β−α = 1, the possible values of ∆ Z EW with larger tan β are included in those with smaller tan β. This is simply because the unitarity bound more strongly constrains the possible non-decoupling effect for

IV. NUMERICAL RESULTS FOR THE HIGGS BOSON DECAY RATES
In this section, we numerically show predictions of the total width and the branching ratios of the Higgs boson at NLO in the HSM, the THDMs and the IDM. After we show these quantities, we demonstrate if these extended Higgs models can be distinguished by the difference of the pattern of deviations in the branching ratios from those of the SM predictions. Similar to the analysis in Sec. III F, we take into account the constraints from the unitarity, the vacuum stability, the conditions to avoid wrong vacua and the S, T parameters. Except for Sec. IV C, we dare not to impose the flavor constraints in order to just see the predictions of deviations in total width and branching ratios. For the THDMs, we introduce the common mass of the additional Higgs bosons

A. Total widths
We first discuss the total width of the Higgs boson h. In Fig. 8, we show the deviation in the total width from the SM prediction in the HSM. We scan the parameters c α , m H and M 2 within 0.95 < c α < 1, 300 ≤ m H ≤ 5000 GeV and 0 ≤ M 2 ≤ m 2 H , respectively. The dependences on c α and m H are then displayed in the left and right panels, respectively. At tree level, the deviation in the width is determined by s 2 α , and it almost corresponds to the upper edge of the distribution in the left panel. The loop effects typically reduce the width by at most about 2% level. In the right panel, it is seen that the magnitude of allowed deviations becomes smaller for larger mass regions, because the large mixing is excluded by the theoretical bounds. We note that the information of the width is important to identify the HSM, because the branching ratios of the Higgs boson are almost the same as those of the SM due to nearly universal suppression of the partial decay rates in the HSM.
In Fig. 9, the deviation in the total width is shown as a function of tan β in four types of the THDMs with s β−α = 0.99 and c β−α < 0 (> 0) in the left (right) panel. We scan the values of M 2 and m Φ . In the left plot, it is seen that except for the Type-I THDM, the width becomes larger as tan β increases, because some of the partial widths have a tan β enhancement, e.g., the h → bb (h → ττ ) mode in the Type-II and Type-Y (Type-II and Type-X) THDMs. In the Type-I THDM on the contrary, the total width approaches to the SM prediction, more precisely s 2 β−α Γ SM , at the large tan β region. We note that the curves are truncated at around tan β = 11 (the same thing also happens in the Type-II and Type-Y THDMs), because of the theoretical constraints. In the case with c β−α > 0 (the right panel), the situation is drastically different from the case with c β−α < 0. The total width has the minimal value at tan β ∼ 7 in the Type-II, Type-X and Type-Y THDMs, due to the cancellation between the s β−α term and the c β−α term in κ f , see Table II. This behavior is remarkably observed in the Type-II and Type-Y THDMs, because the h → bb mode, which is the biggest partial width of h in the SM, follows the behavior explained above. We can also see that at tan β 14, the deviation in the total width becomes zero, as we have κ 2 f 1 for all the types of Yukawa interaction. In the Type-I THDM, the width approaches to the SM value at a large value of tan β as also seen in the case with c β−α < 0. The typical amount of the loop corrections to the total width is a few percent level, which is shown by a width of each curve.
In Fig. 10, we show the m Φ dependence on the deviation in the total width in four types of the THDMs. Here, we scan the values of M 2 and s β−α , while we fix tan β to be 1.5. For c β−α < 0 (left panels), the deviation is distributed in the positive (negative) direction in the Type-II and Type-Y (Type-I and Type-X) THDMs, while for c β−α > 0 (right panels), the situation is opposite. This can be understood by focusing on the deviation in the decay rate of h → bb which is expressed by As we are considering s β−α ∼ 1, the 2ζ b s β−α c β−α term dominantly determines the behavior. We see that the allowed magnitude of the deviation is shrunk at around m Φ = 700 (450) GeV for c β−α < 0 (c β−α > 0), because in the region above m Φ = 700 (450) GeV the unitarity and/or vacuum stability bounds constrain s β−α to be closer to 1. As it is expected from the decoupling theorem, for larger m Φ the magnitude of the deviation is  Fig. 11, we show the total width in the IDM as a function of the charged scalar boson mass m H ± with m A = m H ± . We here take two cases; i.e., (i) m H is fixed to 63 GeV and (ii) m H = m H ± . The case (i) is motivated by the dark matter physics [97,98,110,111],

Finally in
where H can be a dark matter candidate. In this case, the value of M 2 is taken such that the HHh coupling normalized to v becomes around 10 −3 to avoid constraints from dark matter direct detection experiments. In the IDM, the total width does not change from the SM value at tree level, so that any deviation is purely due to loop effects. We can see that in the case (i), the total width monotonically decreases and the deviation is larger as m H ± is getting larger. The black curve is truncated at around m H ± = 700 GeV, because of the unitarity constraint. In the case (ii), the magnitude of the deviation becomes larger up to m H ± 600 GeV, while it becomes smaller above m H ± 600 GeV. The maximal deviation is given at M 2 = 0 for m H ± < 600 GeV, while the unitarity constrains the minimal value of M 2 above m H ± 600 GeV and possible deviations become smaller.

B. Branching ratios
We move on to the discussion of the branching ratios of the Higgs boson h at NLO. For reference, in Table IV we give our results for the branching ratios in the SM.
In the HSM and the IDM, the branching ratios for h are almost the same as those in the SM predictions, because the partial decay rates are universally suppressed by the radiative corrections and the mixing, where the latter does not happen in the IDM. Thus, in the following discussion, we concentrate on the THDMs. In the THDMs, branching ratios can be modified from those in the SM by both tree-level mixing effects parameterized by the scaling factor κ X and loop effects. When s β−α = 1 is taken, branching ratios can significantly be modified from the SM predictions due to the tree level mixing effects, and the pattern of deviations strongly depends on the type of Yukawa interactions. In this case, we may be able to determine the type from the pattern of deviations.
On the other hand, loop contributions to deviations in the branching ratios are relatively smaller than the tree level mixing effects, so that it would be relatively difficult to extract the loop effects.
When s β−α = 1, however, the pure loop effect can be extracted, because the tree level mixing vanishes. Therefore, in the following discussion, we first show the predictions of branching ratios in the THDMs with s β−α = 1 to see how the mixing effects modify them. We then show those with s β−α = 1 in order to extract the size of loop effects.
In Fig. 12 c β−α > 0, because of the fact that some κ f factors become zero, e.g. κ b in the Type-II and Type-Y THDMs, and it makes the value of the total width to be minimal as we saw it in the right panel of Fig. 9. Loop effects due to additional Higgs bosons appear as a width of each curve.
In order to see the deviation in the ratio of the branching ratio from the SM prediction, we introduce the following quantity for the h → XX mode Using the formulae of the partial decay rates at NLO discussed in Sec. III, ∆µ XX can be written in terms of the EW corrections ∆ X EW defined in Eq. (68) in the alignment limit as where BR 0 denotes the tree level branching ratio in the SM. We note that the second term of the right hand side is dominantly determined by ∆ b EW because the branching ratio of h → bb typically has the largest value among all the decay modes. This expression is helpful to understand the behavior of some plots which will be shown below.
In Fig. 13, we show ∆µ f f (f = b, c, τ ) as a function of m Φ in four types of the THDMs with s β−α = 1 and tan β = 1.5. Here, we fix the value of λv 2 = m 2 Φ − M 2 to be 0 (solid curves) and (200 GeV) 2 (dashed curves). We see the decoupling behavior in the large mass region, and observe the peak at around m Φ = 2m t depending on the type of Yukawa interaction and the type of fermions, where the direction of the peak is the same as that for the plots of ∆ f EW shown in Fig. 5. Notice that the peak appearing at m Φ > 2m t in Fig. 5 does not appear in this plot, as we here fix the value of λv 2 .  The m Φ dependence of ∆µ V V (V = W, Z) is shown in the THDMs with s β−α = 1 and tan β = 1.5 (Fig. 15) and 3 (Fig. 16), where the values of ∆µ W W and ∆µ ZZ are almost the same of each other in this case. As in Fig. 13, the value of λv 2 is fixed to 0 (solid curves) and (200 GeV) 2 FIG. 14. Same as Fig. 13, but for tan β = 3.
(dashed curves). The left and right panels show the results in the Type-I and Type-II THDMs, respectively, while the results in the Type-X (Type-Y) THDM are almost the same as those in the Type-I (Type-II) THDM. In all the panels, the value of ∆µ V V approaches to zero in the large m Φ region, because of the decoupling property of the additional Higgs bosons. In the left (right) panel, a peak appears at around m Φ = 2m t , because the EW correction to the partial width of h → bb mode has a peak in the Type-I and Type-X (Type-II and Type-Y) THDMs, see Figs. 3 and 4. We can also see that in the Type-I THDM with tan β = 1.5 the value of ∆µ V V with λv 2 = (200 GeV) 2 is almost the same as that with λv 2 = 0, because the change of ∆ b EW due to taking different values of λv 2 is accidentally cancelled by that of ∆ V EW . For tan β = 3, the value of ∆µ V V with λv 2 = (200 GeV) 2 is slightly smaller than that with λv 2 = 0, because the change of ∆ b EW becomes smaller than that for tan β = 1.5, while ∆

C. Correlations
In the discussions so far, we have seen the deviation in the total width and those in branching ratios in each extended Higgs sector. We now see correlations among the deviations in branching ratios in all the extended Higgs sectors discussed in this paper in order to clarify how we can  distinguish extended Higgs sectors from the precise measurements of the branching ratios.
The branching ratios of the Higgs boson will be measured as accurately as possible at future collider experiments. In particular at the ILC, we can measure the cross section of e + e − → Zh without depending on the decay of h by using the recoil method [112,113]. This makes the measurements of the branching ratios possible without depending on the cross section. At the ILC with the collision energy of 250 GeV and the integrated luminosity of 2 ab −1 (ILC250), the branching ratios are expected to be measured as shown in Table V. We thus consider the situation where the branching ratios are measured to some extent at ILC250, namely we impose the further constraint on the value of ∆µ XX with a given central value and an error in addition to the theoretical constraints which are imposed in the discussion above.
In order to take into account the constraints from B s → X s γ [91,92], we consider the case with larger masses of extra Higgs bosons, i.e, m Φ ≥ 600 GeV in the THDMs as well as that for m Φ ≥ 300 GeV. As discussed in Sec. II, the lower bound on m H ± is about 600 GeV in the Type-II and Type-Y THDMs.
In Fig. 17, we show the correlation between ∆µ τ τ and ∆µ bb in the HSM, four types of the THDMs and the IDM under the additional constraint from ∆µ W W = 0 ± 2% (left) and ∆µ W W = 0 ± 4% (right). The errors of 2% and 4% are taken to consider about 1σ and 2σ level region at ILC250, respectively [18]. Parameters of each model are scanned as it is described in the caption.
We see that the predictions in the HSM and the IDM are given almost at the origin of this plane, because in these models the partial decay rates are almost universally suppressed as we already the h → W W * , also strongly vary at the same time, and then such configurations are constrained by the bound on ∆µ W W . From this figure, we find that if ∆µ τ τ is found to be a several percent, we can distinguish the models considered in this paper. In the following discussion, we focus on the case with ∆µ W W to be constrained at the 2σ level, i.e., allowing 4% uncertainty.
Let us also consider the case where the central value of ∆µ W W is found to be nonzero, and ∆µ W W = 0 is excluded at the 2σ level. In Fig. 18, we show such situations with ∆µ W W = 5.0±4.0% (left) and ∆µ W W = −5.0 ± 4.0% (right). In this setup, predictions of ∆µ W W in the HSM and the IDM are almost zero, so that these models are excluded, while four types of THDMs can explain such a deviation. If the value of ∆µ W W is given to be 5.0 ± 4.0% (left), then four types of Yukawa interactions are well separated of one another, so that we can determine the type by measuring ∆µ τ τ and ∆µ bb in addition to ∆µ W W . We note that the positive value of ∆µ W W essentially comes from the reduction of the other decay rates, especially the h → bb mode, because the partial width of h → W W * reduces by the tree level scaling factor κ V ≤ 1 and the one-loop effect as seen from This is because the negative value of ∆µ W W can be explained by either decreasing the partial width of the h → W W * mode or increasing the other partial widths. Therefore, the branching ratio of h → bb can be either larger or smaller than the SM prediction, as seen in the right panel.
This makes discrimination among the four types of Yukawa interactions difficult as compared to the case with positive ∆µ W W .
In Fig. 19, we show the correlation between ∆µ τ τ and ∆µ cc under the constraint on ∆µ W W with the 4% uncertainty. The central value of ∆µ W W is supposed to be 0 in the upper panel, and to be +5% and −5% in the lower left and right panels, respectively. As compared to Fig. 17 (right), the allowed points on the upper panel are widely distributed in the ∆µ τ τ and ∆µ cc plane, because the branching ratio of h → cc typically has a smaller portion of the total width than that of h → bb.
If we only look at the correlation shown in the upper panel, it seems difficult to distinguish the models unless ∆µ τ τ is given to several percent level. However, we would like to emphasize that by looking also at the corresponding plot shown in Fig. 17 (right), we can distinguish the models. For example, if ∆µ τ τ and ∆µ cc are measured to be small negative and positive respectively, i.e., the second quadrant in this plane, both the Type-II and Type-X THDMs can explain such situation, but these two models may not be distinguished from each other. However, by looking at the correlation between ∆µ τ τ and ∆µ bb with a negative value of ∆µ τ τ , the Type-II and Type-X THDMs have allowed points in the different directions compared to each other. Therefore, we can distinguish all the four THDMs by using the combination of these correlations even if the central value of ∆µ W W is measured to be close to zero.
Similar to Fig. 18, we show the case for nonzero ∆µ W W at the 2σ level in the lower two plots in Fig. 19. We see that for the case with the central value of ∆µ W W to be +5%, four types of THDMs are clearly separated from one another, while the case with the central value of ∆µ W W to be −5% two models are overlapping at some regions of this plane. However, again the correlation between ∆µ τ τ and ∆µ bb helps for further discrimination of the models.
In the above analyses, we constrained the value of ∆µ W W . We now constrain ∆µ bb instead of ∆µ W W . In Fig. 20 Finally, we show the correlation between ∆µ τ τ and ∆µ cc in Fig. 21 using the same setup as in Fig. 20. The shape of the upper panel looks similar to that seen in the upper panel of Fig. 19, while the allowed points in Fig. 21 are distributed in smaller regions than those shown in Fig. 19. This is simply because the foreseen uncertainty of the measurements of ∆µ bb is smaller as compared to that of ∆µ W W . Interestingly, for both the cases of ∆µ bb = 2.5 ± 2% (left) and ∆µ bb = −2.5 ± 2% (right) four types of the THDMs are well separated, and it becomes clearer when m Φ is taken to be greater than 600 GeV.
In this subsection, we have discussed various correlations between the deviations in branching ratios from the SM predictions at NLO. First, if we observe a percent level deviation in one of the decay modes of h, then the HSM and IDM could be ruled out as the branching ratios are almost the same as the SM predictions in these models. Second, the discrimination of four types of the THDMs strongly depends on the situation. If we observe a positive deviation in the branching ratio for the h → W W * mode, the discrimination is possible by looking at the correlation between ∆µ τ τ and ∆µ bb or ∆µ τ τ and ∆µ cc . On the contrary, if we measure a negative deviation in the branching ratio of the h → W W * mode, then the discrimination of four types becomes more complicated as two of four models can overlap with each other in the correlation of ∆µ τ τ and ∆µ bb or ∆µ τ τ and ∆µ cc . However, using three observables ∆µ τ τ , ∆µ bb and ∆µ cc with the results from the direct searches of additional scalar bosons and from flavor experiments, we may be able to separate four types of the THDMs.

V. CONCLUSIONS
We have discussed the total width and the branching ratios of the 125 GeV Higgs boson h at NLO in EW and QCD in the HSM, four types of the THDMs and the IDM. These quantities can be measured at collider experiments as precisely as possible under a given machine performance.
Thus, accurate calculations for the total width and branching ratios are quite important to compare them with future precision data, e.g. at the HL-LHC and the ILC. For the one-loop computation, we systematically applied the on-shell renormalization scheme for each model, in which we adopted the H-COUP program to evaluate numerical values of the renormalized Higgs boson vertices. The analytic expressions for the decay rates of h → ff , h → ZZ * → Zff and h → W W * → W ff are presented at NLO, among which the h → W W * → W ff mode is newly computed in this paper.
We also provided the decay rates of the loop induced processes; i.e., h → γγ, h → Zγ and h → gg with QCD corrections at NLO.
We have shown that in the HSM and the IDM, all the partial decay rates are almost universally suppressed by both the tree level mixing (for the HSM) and the one-loop effect, so that the branching ratios remain almost the same values as those in the SM. Thus, if deviations in the branching ratio from the SM prediction (denoted as ∆µ XX for the decay h → XX) are found, then we may be able to exclude the HSM and IDM. On the contrary, when we observe the deviation in the total width but not in the branching ratios, then it could be a smoking gun signature to identify these two models. We also have found that if we observe a positive deviation in the branching ratio of the h → W W * mode, four types of the THDMs can be well separated from one another from the correlation between ∆µ τ τ and ∆µ bb or ∆µ τ τ and ∆µ cc . If we observe a negative deviation in the branching ratio of the h → W W * mode, some of the THDMs can overlap in the ∆µ τ τ and ∆µ bb plane, but this can be disentangled by further looking at another correlation, such as ∆µ τ τ and ∆µ cc . Using the synergy between the direct and indirect searches, the parameter space of extended Higgs models can be effectively narrow down.
Finally, we would like to mention that the H-COUP version 2.0, where all the NLO computations for the decay rates presented in this paper are implemented, will be publicly available in near future [105].
The analytic expressions for each contribution to the hZγ vertex are given bŷ where B i , C i and C ij are the Passarino-Veltman's functions [114]. In this paper, we follow the convention of these functions given in Ref. [31]. Here, we use the shorthand notation for the C functions defined by C i,ij (m) ≡ C i,ij (p 2 1 , p 2 2 , q 2 ; m, m, m) and C 1223 ≡ C 12 + C 23 . The form factors for the renormalized hγγ vertex is obtained from the expressions of the hZγ vertex by the replacement of (g Z , c 2 W , s 2 W , v F ) → (e, 1, −1, Q F ). Finally, the hgg vertex is induced only from the quark loop. Thus, it is expressed bŷ where a and b represent the color index.

Appendix B: Renormalized V ff vertices
We give analytic formulae of the renormalized V ff (V = Z, W ) vertices, which appear in the decay rates of h → V V * → V ff at NLO given in Eq. (60) and (64). In the limit of massless external fermions, expressions of these vertices are common to those in the SM.
The renormalized V ff (V = Z, W ) vertices can be decomposed in the massless limit for external fermions asΓ where p µ 1 (p µ 2 ) is the incoming four-momentum of the fermion (anti-fermion), and q µ is the outgoing four-momentum of the gauge boson. The gauge coupling g V is g Z and g/ √ 2 for Z and W , respectively. Similar to Eqs. (46) and (48), we can further decompose these vertices into the tree level part and 1-loop part: The tree level contribution is given by analytic expressions for T V hf f are given as follows: T where t V = m 2 h + m 2 V − s − u. For the calculation of box diagrams, we define the Passarino-Veltman's D functions: where N 1 = k 2 − m 2 1 , N 2 = (k + p 1 ) 2 − m 2 2 , N 3 = (k + p 1 + p 2 ) 2 − m 2 3 , N 4 = (k + p 1 + p 2 + p 3 ) 2 − m 2 4 . We note that in our calculation up to the second rank tensors D µν appear, and these functions are UV finite. The first and second rank tensor functions are decomposed into the following scalar coefficients: In the following, we shortly express the D functions by D 0,i,ij (p 2 1 , p 2 2 , p 2 3 , (p 1 + p 2 + p 3 ) 2 , (p 1 + p 2 ) 2 , (p 2 + p 3 ) 2 ; a, b, c, d) ≡ D 0,i,ij (p 2 1 , p 2 2 , p 2 3 , (p 1 + p 2 + p 3 ) 2 , (p 1 + p 2 ) 2 , (p 2 + p 3 ) 2 ; m a , m b , m c , m d ). The contribution from box diagrams B V can be expressed as where Φ 4 is the four body phase space function. The first integral denoted as S is performed up to the cutoff of the photon energy ∆E, while the second integral is done from ∆E to the maximal value of the photon energy. The ∆E dependence in the NLO decay rate, of course, disappears after summing up the soft and hard photon parts.
The soft-photon part is calculated using the eikonal approximation by which the amplitude can be expressed by the product of the Born amplitude and the soft-photon factor. Then, we can separately perform the integration with respect to the 3-body phase space and the photon phase space. Therefore, the soft-photon part is expressed as where