Next-to-leading order corrections to decays of the heavier CP-even Higgs boson in the two Higgs doublet model

We investigate the impact of electroweak (EW), scalar and QCD corrections to the full set of decay branching ratios of an additional CP-even Higgs boson ($H$) in the two Higgs doublet model with a softly broken $Z_2$ symmetry. We employ the improved gauge independent on-shell scheme in the renormalized vertices. We particularly focus on the scenario near the alignment limit in which couplings of the discovered 125 GeV Higgs boson ($h$) coincide with those of the standard model while the $Hhh$ coupling vanishes at tree level. The renormalized decay rate for $H\to hh$ can significantly be changed from the prediction at tree level due to non-decoupling loop effects of additional Higgs bosons, even in the near alignment case. We find that the radiative corrections to the branching ratio of $H\to hh$ can be a few ten percent level in the case with the masses of additional Higgs bosons being degenerate under the constraints of perturbative unitarity, vacuum stability and the EW precision data. Further sizable corrections can be obtained for the case with a mass difference among the additional Higgs bosons.


I. INTRODUCTION
Current data at LHC show that properties of the discovered Higgs boson (h) [1,2] are consistent with those of the standard model (SM) Higgs boson within experimental uncertainties [3,4]. In addition, there is so far no report which clearly indicates signatures of new particles from collider experiments. On the other hand, new physics beyond the SM must exist, because there are phenomena which cannot be explained in the SM such as neutrino oscillations, dark matter and baryon asymmetry of the Universe. In such new physics models, the Higgs sector is often extended from the minimal one assumed in the SM, and its structure depends on new physics scenarios.
Therefore, clarifying the structure of the Higgs sector is important to determine the direction of new physics.
One of the useful ways to determine the structure of extended Higgs sectors is to measure deviations in various observables of h from SM predictions. We can extract information on the structure such as the number of additional Higgs fields and their representations from the pattern of the deviations [5]. Furthermore, the mass scale of extra Higgs bosons can be deduced from the size of the deviations [6][7][8][9][10][11][12][13][14][15]. Because the precise measurements of cross sections, decay branching ratios (BRs) and the width of the Higgs boson will be performed at future collider experiments such as the high luminosity LHC (HL-LHC) [16], International Linear Collider (ILC) [17][18][19][20], Circular Electron Positron Collider (CEPC) [21] and Future Circular Collider (FCC-ee) [22], precise predictions of these observables including radiative corrections are necessary to compare these measurements.
In this case, decays of additional Higgs bosons into a lighter Higgs boson, "Higgs to Higgs decays", are suppressed due to the small mixing with h, which can be regarded as a golden channel for the direct searches, see e.g., H → hh and A → Zh in Refs. [48][49][50][51][52][53][54][55][56]. It has been shown in Ref. [56] that by considering the synergy between direct searches for additional Higgs bosons at the HL-LHC and precise measurements of the h couplings at the ILC, wide regions of the parameter space on extended Higgs sectors can be explored.
Such a synergy analysis discussed in Ref. [56] has been performed at leading order (LO) in EW interactions. However, tree level analyses might not be sufficient because of the following reasons.
First, they have analyzed the excluded and testable parameter regions inputting the tree level hV V coupling (V = W, Z) as constraints from indirect searches. If we take into account loop corrections, the hV V coupling is modified. Second, decay rates of the Higgs to Higgs decay processes can be significantly changed at one-loop level from the prediction at tree level, because some of one-loop diagrams do not vanish at the alignment limit and can be sizable due to the non-decoupling effects of additional Higgs bosons. Therefore, in order to see the synergy in a more realistic way, radiative corrections to both the discovered Higgs boson couplings and the decay BRs of additional Higgs bosons should be taken into account. Radiative corrections to decays of additional Higgs bosons have been studied in Refs. [13,[57][58][59][60][61] for a heavy CP-even Higgs boson, in Refs. [14,62] for a CP-odd Higgs boson and in Refs. [12,[63][64][65][66] for singly charged Higgs bosons.
In this paper, we investigate the impact of next-leading-order (NLO) corrections in EW and scalar interactions to the decay BRs of an additional CP-even Higgs boson (H) in the two Higgs doublet model (THDM) with a soft broken Z 2 symmetry [67,68] as a simple but important example.
We compute decay rates of H → hh, H → ff , H → V V , H → γγ, H → Zγ and H → gg at one-loop level based on the improved on-shell renormalization scheme without gauge dependences, while QCD corrections are also included in decay rates of H → qq, H → γγ, H → Zγ and H → gg.
We study the relation between model parameters; e.g., such as masses of additional Higgs bosons and mixing parameters, and the effects of radiative corrections for the H → hh decay. 1 Moreover, we clarify how the one-loop corrections can significantly change the tree level predictions in the scenario with nearly alignment, and the correlation between the deviation of the h → W W * decay from the SM prediction and the one-loop corrected decay rate of H → hh under theoretical and experimental constraints. We also study the correlation between the BR of the H → hh decay and the deviation of the one-loop corrected hhh vertex of the THDM [6, 8,9] from that of the SM, which is important for testing the EW baryogenesis scenario [69][70][71][72][73]. Finally, we give results of decay BRs of other decay processes in several scenarios. This paper is organized as follows. In Sec. II, we define the THDMs and we mention the current situation of experimental constraints. In Sec. III, we define our renormalization scheme, and give renormalized vertices Hhh, Hff (f = t, b, τ, c) and HV V (V = W, Z). In Sec. IV, we give formulae of one-loop corrected decay rates of H → hh, H → ff and H → V V . Numerical evaluations of these decay rates are shown in Sec. V. Discussions and conclusions are respectively given in Sec. VI and VII. In Appendix, we present explicit analytic formulae of 1PI diagrams for the Hhh vertex.

II. MODEL
The Higgs sector is composed of two isospin doublet fields Φ 1 and Φ 2 with hypercharge Y = 1/2.
In order to avoid flavor changing neutral currents mediated by Higgs bosons at tree level, we impose a discrete Z 2 symmetry to the model, where Φ 1 and Φ 2 are transformed to +Φ 1 and −Φ 2 [67,68], respectively.
The Higgs potential is given under the Z 2 symmetry by where the m 2 3 term softly breaks the Z 2 symmetry. In general, m 2 3 and λ 5 are complex, but we assume these parameters to be real. The Higgs fields are parameterized as where v 1 and v 2 are the vacuum expectation values (VEVs) with v = v 2 1 + v 2 2 ≃ 246 GeV, and their ratio is parameterized as tan β = v 2 /v 1 . The mass eigenstates of the Higgs bosons are defined as follows, where with s θ ≡ sin θ and c θ ≡ cos θ. In Eq. (3), H ± , A, H and h are the physical mass eigenstates, while G ± and G 0 are the Nambu-Goldstone bosons. After imposing the tadpole conditions, the masses of the charged and CP-odd Higgs bosons are expressed as where M 2 is defined as M 2 = m 2 3 /(s β c β ). The masses of the CP-even Higgs bosons and the mixing angle β − α are given by where with M 2 ij being elements of the mass matrix in the basis of R(β)(h 1 , h 2 ) T . In the following, we identify h as the discovered Higgs boson with a mass of 125 GeV. The 8 parameters in Eq. (1) can be expressed by the following 8 parameters, where we define tan β > 0 and 0 ≤ s β−α ≤ 1, while c β−α can be either positive and negative.
For the later discussion, we here give expressions of the triple scalar couplings as where they are defined as The kinetic term of the Higgs fields is where D µ is the covariant derivative given by The gauge-gauge-scalar interaction terms in the mass eigenstates are extracted as where g Z is defined as g Z ≡ g/ cos θ W , and κ φ V (φ = h, H) is the scaling factor obtained as The most general Yukawa interaction is given under the Z 2 symmetry as where Φ i , Φ j and Φ k are either Φ 1 or Φ 2 . These labels of the Higgs doublet fields are determined by fixing the Z 2 charge of each field as summarized in Tab. I [83,84]. The Yukawa interaction terms are given in terms of the mass eigenstates of the Higgs bosons as where I f = 1/2 (−1/2) for f = u (d, e), and V ud is the Cabibbo-Kobayashi-Maskawa matrix element. Scaling factors in Eq. (20) are expressed as, where ζ f are given in Tab. I.
In the limit of s β−α → 1, all couplings of h become the same as those of the SM at tree level.
We call this limit as the alignment limit. As mentioned in Sec. I, we are interested in the nearly alignment case, so that we introduce a parameter x defined by where the limit x → 0 corresponds to the alignment limit. In the nearly alignment case; i.e., x ≃ 0, the scaling factors and the triple scalar couplings can be expressed as with We mention constraints from current experimental data on the THDM.
• Electroweak precision data New physics effects on the EW oblique parameters are expressed by the S, T and U parameters [85,86]. It has been known that the T parameter represents the violation of the custodial SU (2) V symmetry [87][88][89][90]. One of the ways to restore the SU (2) V symmetry in the potential is to take m H ± = m A , by which quadratic-power like dependences of Higgs boson masses on the T parameter disappear. 2 By using the global fit of EW parameters [93], new physics effects on the S and T parameters under U = 0 are constrained by with the correlation factor of +0.91 [93].

• Signal strength of the Higgs boson h
Measurements of the signal strengths of h constrain mixing parameters α and β as seen in Eqs. (18) and (21). According to Refs. [4,94], the parameter regions with −0.25 c β−α 0.01 for tan β ≃ 2 are allowed in the Type-I THDM. In the Type-II, X and Y THDMs, the constraint is given by; e.g., |c β−α | 0.1 for tan β ≃ 2, and stronger bounds are taken for larger tan β.
• Direct searches for the additional Higgs bosons at LHC So far, there has been no report for the discovery of the additional Higgs bosons at LHC, so that lower limits on their masses have been taken. In the THDMs, constraints from the direct searches have been studied by using LHC Run-II data in; e.g., Refs. [54,56,95]. In the nearly alignment region, the search for the pp → A → Zh process typically excludes the largest region of the parameter space in the four types of THDMs. For instance, m A 2m t (900 GeV) are excluded in the case with s β−α = 0.99 and tan β ≃ 5 (2) [56]. For larger values of 1−s β−α , the constraint by this mode becomes stronger. In Type-II and Type-Y, the search for the pp → bbA → bbZh process also provides severe constraints on the parameter space particularly for the case with larger tan β values. For instance, in the case with s β−α = 0.99 and tan β = 10, m A 500 GeV have been excluded [56]. In addition to the A → Zh mode, the H → hh mode can also exclude wide regions of the parameter space. This, however, strongly depends on the value of M 2 , and can be changed by the radiative corrections which will be discussed later.

• Flavor experiments
The mass of the charged Higgs bosons is restricted by the B → X s γ decay. In Type-II and Type-Y, m H ± 800 GeV is excluded for tan β > 2 [96,97]. On the other hand, in Type-I and Type-X, m H ± 400 GeV is excluded for tan β = 1 [96,97]. Moreover, in the Type-II case, data of the B s → µµ decay exclude regions with large tan β [98-100]; e.g., tan β 25 for m H ± = 1 TeV.

III. RENORMALIZED VERTICES
In this section, we discuss the renormalized vertices for H based on the improved on-shell scheme [9], where gauge dependences in counter terms of mixing angles are removed by adding pinch terms. This treatment is not applied to wave function renormalizations but counter terms from shifts of parameters in the Lagrangian.
In the following discussion, we decompose each renormalized form factorΓ HXX (X = h, W, Z, t, b, τ, c) asΓ where with Γ Tree HXX , δΓ HXX , Γ 1PI HXX and Γ Tad HXX being contributions from tree level diagrams, counter terms, 1PI diagrams and tadpole diagrams inserted to tree level vertices, respectively.

A. Renormalization conditions
We discuss renormalization conditions in order to determine the counter terms for the renormalized HXX vertices. The renormalized two-point functions for neutral scalar bosons are expressed where δZ i , δC ij , δα f and δβ f come from the field shifts defined as We here distinguish δα f (δβ f ) from δα (δβ) which appears from the parameter shift; i.e., α → α+δα By imposing the following on-shell conditions for the above renormalized two-point functions; the counter terms δm 2 i and δZ i are determined as By imposing the on-shell conditions for the mixing two-point functions aŝ we obtain where we take δC hH = δC Hh ≡ δC h and δC AG 0 = δC G 0 A ≡ δC A .
As mentioned above, the counter terms δα and δβ should include the pinch term; where explicit formulae of Π PT Hh (p 2 ) and Π PT AG (p 2 ) are given in Ref. [9]. The counter term δM 2 cannot be determined by using the on-shell conditions discussed above.
Thus, we apply the MS scheme to the renormalization of the hhh vertex, and then δM 2 is determined as where ∆ Div = 1/ǫ − γ E + log 4π with γ E being the Euler's constant. In the second line, Div(T 1PI h(H) ) indicates the divergent part of tadpole type diagrams whose explicit formulae are given in Ref. [8].
In this prescription, renormalized quantities including δM 2 such as the decay rate of H → hh have a dependence on the renormalization scale µ.
Finally, the counter terms of the EW parameters δv, δm 2 V , δZ V , δm f and δZ f V are given in Ref. [8].

B. Renormalized Hhh vertex
The contributions from the tree level diagram and from counter terms are given by where λ hhh , λ Hhh and λ HHh are given in Eqs. (13)- (15), and δλ Hhh is given by with G α and G β being Explicit formulae of contributions of 1PI diagrams Γ 1PI Hhh are given in Appendix. As mentioned in the previous subsection, the Hhh coupling depends on the renormalization scale µ.

C. Renormalized Hff vertex
The renormalized Hff (f = t, b, c, τ ) vertex is expressed in terms of the following form factors, where p µ 1 and p µ 2 are the incoming momenta of external particles f andf , respectively, and q µ (= p µ 1 + p µ 2 ) is the outgoing momentum of H. For the case with on-shell fermions; i.e., p 2 1 = p 2 2 = m 2 f , the following relations hold: The contributions from the tree level diagram and from counter terms are given by where the index i runs over i = {P, V 1 , V 2 , A 1 , A 2 , T, P T }, and The mixing factor ζ f and κ φ f (φ = h, H) are given in Sec. II.

D. Renormalized HV V vertex
The renormalized HV µ V ν (V = W, Z) vertex is composed of three types of form factors ex- where p µ 1 and p ν 2 are incoming momenta of the weak bosons, and q µ is the outgoing momentum of H.
The contributions from the tree level diagram and from counter terms are given by Since we drop the one-loop squared term |M 1-loop | 2 which corresponds to next-to-NLO (NNLO) where λ Hhh and Γ 1-loop Hhh are given in Eq. (14) and (31), respectively. The LO contribution Γ LO [H → hh] is given by In Eq. (63), ∆r represents EW radiative corrections to the VEV v, which should be added to the decay rate, because we choose α em , m Z and G F as the EW input parameters in the renormalization calculation. The analytic expression for ∆r is given by [101] The same procedure is also applied to the other decay rates given below.
In the nearly alignment region; i.e. x ≃ 0, the decay rate given in Eq. (63) can be expanded in withΓ In Eq. (66), λ Hhh,1 and λ Hhh,2 are given in Eq. (27). The bosonic loop and fermionic loop contri-butions to ∆ EW 0 are respectively expressed as with N f c being the color factor and The functions A[X], B 0 [q 2 ; X, Y ] and C 0 [X, Y, Z](≡ C 0 [p 2 1 , p 2 2 , q 2 ; X, Y, Z]) represent Passarino-Veltman functions [102]. In Eq. (69), terms in the second line are from Γ 1PI Hhh , those in the third line are from δC h and δα, and those in the fourth line are from δβ. In particular, for ∆ EW 0,B diagrams shown in Fig. 1 give the dominant contribution. For ∆ EW 0,F , light fermion loops can be neglected because there is the overall factor m 2 f , so that the top quark loop gives the dominant effect, which is proportional to ζ t = cot β in four types of Yukawa interaction. Therefore, the major quantum effects in the nearly alignment scenario do not depend on the types of Yukawa interaction.
Let us here consider the x 2 term in Eq. (66), in which the first term in the parentheses is the contribution from the tree level diagram, while the others are those from the one-loop diagrams which are suppressed by the loop factor (1/16π 2 ). If we regard the effect of the loop suppression factor as the small expansion parameter x, the decay rate can be approximately rewritten as The renormalization scale µ appears in the one-loop corrected decay rate of H → hh as where The tadpole terms T µ-part φ (φ = h, H) give the dominant contribution to the µ dependence, which are proportional to cot 2β, so that the magnitude of ∆ µ grows as tan β becomes large. We also mention that ∆ µ is roughly proportional to m 4 Φ , which comes from the contributions from the tadpole diagram of h. In the following, we fix as µ = m H in the numerical calculations.

B. Decay rate of H → ff
The decay rate of the process H → ff with NLO corrections in EW and scalar interactions and QCD-NNLO corrections can be written as where Γ LO [H → ff ] is the decay rate at tree level expressed as For the hadronic decays; i.e., f = q, the QED correction is given in the MS scheme as [106] ∆ QED where we fix the renormalization scale µ as µ = m H in numerical calculations in this paper.
For decays into a quark pair, we implement NNLO-QCD corrections in the MS scheme according to the formulae summarized in Ref. [56].

C. Decay rate of H → V V
We give the formulae of decay rates into a pair of on-shell gauge bosons; i.e., H → V V , with NLO corrections in EW and scalar interactions. We express the decay rates as with c V = 1 (2) for W (Z). In Eq. (81), ∆ EW HV V indicates EW loop contributions expressed as In the above expression, the second term is the contribution from wave function renormalizations of external vector bosons, which are non-zero in our on-shell renormalization scheme. The tree level contribution of the HV V vertex Γ Tree,1 HV V is given in Eq. (59), and the explicit formula ofΠ V V is given in Eq. (56) of Ref. [8]. The term ∆ brem in Eq. (81) indicates the contribution from the real photon bremsstrahlung which is needed in order to remove infrared (IR) divergence from virtual photon loop diagrams. The explicit formula of the contribution is given by [107] ∆ brem = α em π 2 − 1/r with r = m 2 H /(4m 2 W ) and ρ ± = √ r ± √ r − 1. In Eq. (85)  contributions can be expanded as with ǫ ≡ m 2 h /m 2 H and C ≡ 2 − π/ √ 3 ≃ 0.186. From this approximate formula, it is seen that ∆ EW 0,B is enhanced by m 4 H for the case with M ≪ m H due to the non-decoupling effect. On the other hand, for M ≃ m H such an enhancement is highly suppressed by the factor of (1 − M 2 /m 2 H ) 2 and thus ∆ EW 0,B is roughly proportional to (m 2 H − M 2 ) 2 /v 4 . It can also be seen that Eq. (86) has the factor cot 2β, so that the magnitude of the NLO corrections grows as tan β increases. We note that the NLO contributions come from the cross term of the amplitude from the tree level and one-loop contributions, so that the sign of the NLO contributions changes depending on the sign of c β−α . Namely, if c β−α is positive (negative), the NLO contributions increase (decrease) the decay rate. Features of the loop corrections as those described here can be concretely confirmed by the following Figs. 2 and 3.     Next, we investigate the correlation between the BR of H → hh and the BR of h → W W * → W ff ′ which is expected to be measured with 2% accuracy at the ILC with the collision energy being 250 GeV [18]. In order to parametrize the deviation in the BR of h → W W * from the SM prediction, we introduce ∆µ W W defined as GeV, respectively. In some of the panels, several colored regions do not appear because of no allowed region by the constraints. Results with ∆µ W W < 0 (∆µ W W > 0) correspond to those with c β−α > 0 (c β−α < 0), because the partial decay width of h → bb is enhanced (suppressed). In the case with heavier H, predictions are well determined to be narrower regions, because the allowed range of M 2 by the theoretical constraints is shrunk. It is seen that for ∆m = 0 and ∆µ W W < 0 In the case with non-zero mass difference among additional Higgs bosons; i.e., ∆m = 0, the behavior of loop corrections can drastically be different from that in case with ∆m = 0. As ∆m increases, the difference from the tree level prediction becomes more significant than that in the degenerate mass case. For results of m H = 300 GeV, predictions including the NLO corrections can be about 30 % larger than tree level predictions if ∆m is larger than 300 GeV. If ∆m is non-zero, the effect of NLO corrections can increase the BR in the both cases with c β−α > 0 and c β−α < 0.
However, due to theoretical constraints, allowed regions become smaller as m H increases.
If ∆µ W W is larger than 2%, it can be observed as a deviation from the SM prediction by the precision measurements at the ILC with the center of mass energy √ s to be 250 GeV [18]. However, even if the deviation of the h → W W * decay is too small to be observed by the ILC, it might be possible to explore H via the H → hh process at the HL-LHC. If H is lighter than 2m t and ∆m is non-zero, the H → hh decay mode can be dominant in large parameter regions.
It is known that similar non-decoupling effects also appear in the loop corrected hhh coupling [6, 8,9]. The physics of the hhh coupling is strongly related with the EW baryogenesis, because the strong first order phase transition can lead to a large deviation in the hhh coupling from the SM prediction at zero temperature [69][70][71][72][73]. The hhh coupling can be extracted from the measurements of the double-Higgs production at hadron, lepton and photon colliders as discussed in Ref. [108].
The measurement accuracy of the hhh coupling is expected to be about 27% at the ILC with √ s = 500 GeV [18]. In Fig. 5, we show the correlation between ∆µ W W and the deviation in the renormalized hhh vertex in the Type-I THDM from that in the SM, in the case with m H = 300 GeV and tan β = 3. We calculate the renormalized hhh vertex using H-COUP [23,24], excepting parameter regions causing BR(H → hh)< 0. The color definitions of the regions are the same as specified in Fig. 4. In order to examine the correlation with the BR of H → hh, we also place the panel which is shown in Fig. 4. It can be seen that the deviation of the hhh coupling is almost determined by the magnitude of ∆m. The larger ∆m causes the larger deviation in the hhh coupling, since it is caused by a larger non-decoupling effect. Namely, the structure of the non-decoupling effects is the same as those of the H → hh decay. Such parameter regions are common with regions where the hhh coupling shifts from the SM predictions significantly so that the H → hh search at the HL-LHC might also be used to test the EW baryogenesis scenario multi-directionally.  at LHC in Refs. [50,[109][110][111][112][113][114][115][116][117][118].

VI. DISCUSSIONS
We discuss the direct search for the additional Higgs bosons at future collider experiments. At the HL-LHC, H or A is mainly produced by the gluon fusion process and the associated production with bb. 3 The parameter region expected to be explored via these single productions has been studied in Ref. [56], where the analysis has been done at LO in the EW interaction. It goes without saying that the search for the additional Higgs bosons can also be done at the ILC energy upgrade, where the collision energy √ s can be extended to be up to 1 TeV [17]. They can mainly be produced in pairs as e + e − → HA and e + e − → H + H − up to 500 GeV for the degenerate mass case. As we have shown in this paper, the BRs of H can significantly be changed by the NLO corrections in EW and scalar interactions. Thus, it is quite important to include such effects in the exploration of the additional Higgs bosons at the HL-LHC and the ILC. We will upgrade the numerical program H-COUP [23,24] such that the decay BRs for the additional Higgs bosons are calculated including EW, scalar and QCD corrections based on this paper (H), Ref. [120] (A) and Ref. [66] (H ± ). We will then be able to discuss the synergy between the direct search and the precise measurement of h, which will be left for future works.
Finally, we would like to comment that a portion of the parameter space shown in our numerical results is excluded by the current experimental data from the additional Higgs boson searches at LHC and the measurement of the signal strength of h given in Sec. II. For example, the region ∆µ W W 0 are excluded by taking into account the constraints on α and β in the THDMs from the signal strength data [94]. It is, however, seen that there is a large discrepancy between the region excluded by the observed data and that by the expectation of the MonteCarlo analysis.
Thus, the observed exclusion can be drastically changed by accumulating more data. Therefore, we have investigated wider regions of the parameter space than the allowed ones by the current experimental data in the numerical calculations.

VII. CONCLUSIONS
We have computed the decay rates of the additional CP- and where definitions of combinations of C -functions are given as [4] CMS collaboration, Combined Higgs boson production and decay measurements with up to 137 fb −1 of proton-proton collision data at √ s = 13 TeV, .