Probing extended Higgs sectors by the synergy between direct searches at the LHC and precision tests at future lepton colliders

We discuss a possibility that the parameter space of the two Higgs doublet model is significantly narrowed down by considering the synergy between direct searches for additional Higgs bosons at the LHC and its luminosity upgraded operation and precision measurements of the Higgs boson properties at future electron-positron colliders such as the International Linear Collider. We show that, in the case where the coupling constants of the discovered Higgs boson are slightly different from the predicted values in the standard model, most of the parameter space is explored by the direct searches of extra Higgs bosons, in particular for the decays of the extra Higgs bosons into the discovered Higgs boson, and also by the theoretical arguments such as perturbative unitarity and vacuum stability. This can be done because there appears an upper limit on the mass of the extra Higgs bosons as long as the deviation exists in the Higgs boson coupling. We also show that in the alignment limit where all the Higgs boson couplings take the standard model like values most of the parameter space cannot be excluded because most of the Higgs to Higgs decays are suppressed and also there is no upper limit on the masses from the theoretical arguments.


I. INTRODUCTION
The current observations at the LHC experiments indicate that properties of the discovered Higgs boson with the mass of 125 GeV coincide with those predicted in the standard model (SM) [1,2]. This, however, does not mean that the Higgs sector in the SM, which plays an essential role in the electroweak (EW) symmetry breaking, is verified. While the minimal Higgs sector is composed of one Higgs doublet field in the SM, there is no principle to determine the structure of the Higgs sector. In fact, it is possible to consider a variety of non-minimal Higgs sectors. Extended Higgs sectors are often introduced in new physics models which can explain observed phenomena beyond the SM, such as neutrino oscillations, dark matter and baryon asymmetry of the Universe. In addition, they also appear in some of the new paradigms motivated from a theoretical problem in the SM; e.g., the hierarchy problem. Therefore, new physics beyond the SM can be revealed by thoroughly testing the Higgs sector. In addition to the direct searches, extended Higgs sectors can be explored by measuring various properties of the discovered Higgs boson such as cross sections, the width, branching ratios and coupling constants. If deviations from the SM are observed, we can extract upper limits on the mass scale of the second Higgs boson by taking into account theoretical consistencies. Furthermore, by looking at the pattern of the deviation we can extract the structure of the Higgs sector; e.g., the representation of the weak isospin, the number of Higgs fields, and symmetries. To this end, precision measurements of the Higgs boson couplings are most important. Although the current accuracy of the measurements is not enough, typically order 10 (20) percent level for the Higgs boson coupling to weak bosons (third generation fermions) [1,2], it is expected to be improved at the HL-LHC [18] and further significantly at future lepton colliders; e.g., the International Linear Collider (ILC) [19][20][21][22], the Future Circular Collider (FCC-ee) [23] and the Circular Electron Positron Collider (CEPC) [24].
It goes without saying that accurate calculations of the Higgs boson couplings are inevitable in order to compare theory predictions with the future precision measurements. It has been well known that QCD corrections to Higgs boson couplings with quarks or gluons can be quite large. For example, QCD corrections to the decay rate of the Higgs boson into gluons at the next-to-leading order (NLO) is about 70% level [25][26][27]. Thus, QCD corrections must be included for calculations, by which we can discuss the deviation from the SM prediction. On the other hand, EW corrections are typically much smaller than QCD ones, but they have a sensitivity to the structure of the Higgs sector, particularly non-decoupling nature of extra scalar fields. So far, EW corrections to Higgs boson couplings and/or decays have been investigated in models with extended Higgs sectors such as those with extra singlets [28][29][30][31][32][33], doublets [29,[32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] and triplets [50][51][52][53][54]. Therefore, calculations with both QCD and EW corrections are quite important for the precision measurements in near future, and several numerical tools have been available; e.g., H-COUP [55,56], 2HDECAY [57] and Prophecy4f [58].
In this paper, we investigate the impact of the combined study of direct searches for new particles at hadron colliders and precision measurements of Higgs boson couplings at future lepton colliders. We perform such study including higher-order QCD corrections. We consider two Higgs doublet models (THDMs) as a representative extended Higgs model.
The observed Higgs boson couplings are consistent with those in the SM under current experimental and theoretical uncertainties [1,2], so that this fact gives a strong motivation to investigate the alignment scenario where the Higgs boson couplings are nearly or exactly SM like. In the near alignment region, the decays of the extra Higgs bosons into the discovered Higgs boson such as A → Zh and H → hh can be dominant, and at the same time the discovered Higgs boson couplings can deviate from the SM predictions. These decay modes of extra Higgs bosons can be well tested at the HL-LHC [93], by which we can set a lower limit on the masses of extra Higgs bosons. In addition, we can impose an upper limit on the masses [87,94,95] when deviations of the Higgs boson couplings are found at future lepton colliders.
We show that by utilizing the synergy between the direct search for additional Higgs bosons and the precision measurement of the Higgs boson couplings a large portion of the parameter space can be explored in the near alignment region. We also show that in the alignment limit; i.e., all the Higgs boson couplings are exactly same as the SM values, plenty of the parameter space still remains even if the mass of the additional Higgs bosons are around the EW scale. This is because most of the Higgs to Higgs decays are prohibited and also there is no upper limit on masses of additional Higgs bosons. This paper is organized as follows. In Sec. II, we define the THDMs and give the Higgs potential, the kinetic terms and the Yukawa interactions. Theoretical constraints from perturbative unitarity and vacuum stability are also discussed. Constraints from flavor physics and previous colliders are summarized. Sec. III is devoted to the discussion for decays of the Higgs bosons. We first give the analytic expressions of the decay rates with higherorder QCD corrections and then numerically show total widths and branching ratios of the Higgs bosons. In Sec. IV, we show the excluded region of the parameter space from the direct searches at the LHC Run-II experiments. In Sec. V, we discuss how the parameter space is widely explored by combining direct searches at the HL-LHC and precision measurements of the Higgs boson couplings at future lepton colliders. Conclusions are given in Sec. VI. In Appendix, we present the analytic expressions for the perturbative unitarity and the vacuum stability conditions (Appendix A) and the decay rates of the Higgs bosons at the leading order (LO) (Appendix B).

II. MODEL
We discuss the THDM, whose Higgs sector is composed of two isospin doublet scalar fields Φ 1 and Φ 2 . In order to avoid flavor changing neutral currents (FCNCs) at tree level, we impose the Z 2 symmetry [96] (Φ 1 → +Φ 1 , Φ 2 → −Φ 2 ) which can be softly-broken by a dimensionful parameter in the Higgs potential. The most general Higgs potential under the Z 2 symmetry is given by where m 2 3 is a soft-breaking parameter of the Z 2 symmetry. Throughout this paper, we assume CP-conservation in the Higgs sector, so that the m 2 3 and λ 5 parameters are taken to be real. It is convenient to define the Higgs basis [97][98][99] as with the rotation angle β being determined by tan . We introduced short-hand notation for trigonometric functions s θ ≡ sin θ and c θ ≡ cos θ. In where h can be identified as the discovered Higgs boson with the mass of 125 GeV.
The squared masses of the physical Higgs bosons are expressed as follows: Type-X (lepton specific) where M 2 ≡ m 2 3 /(s β c β ) and M 2 ij are the elements of the squared mass matrix in the basis of (h 1 , h 2 ) given by with λ 345 ≡ λ 3 + λ 4 + λ 5 . The mixing angle β − α can also be expressed by these matrix elements as We can choose the following six variables as the free parameters: where we define 0 < β < π/2 and 0 < β − α < π such that tan β > 0 and 0 < s β−α ≤ 1.
The kinetic terms for the Higgs doublets are written in the Higgs basis as The covariant derivative is defined by D µ = ∂ µ −igI a W a µ −ig Y B µ , with the SU(2) L generator I a (a = 1-3) and the hypercharge Y , from which electric charge Q is derived by Q = I 3 + Y .
In the expression of D µ , W a µ (g) and B µ (g ) denote the SU(2) L and U(1) Y gauge bosons (coupling), respectively. The W ± bosons and the neutral gauge bosons are then identified as Under the Z 2 symmetry, the Yukawa interaction terms are expressed as whereΦ u = iσ 2 Φ * u with σ 2 being the second Pauli matrix, and Φ u,d,e denote Φ 1 or Φ 2 . We here do not explicitly show the flavor indices. Using the Higgs basis Eq. (2), they can be rewritten as the so-called Type-I, Type-II, Type-X and Type-Y as shown in Table I. For the later convenience, we introduce the scaling factors κ φ X which are defined by the ratio of the Higgs boson couplings at tree level: where h SM is the Higgs boson in the SM. From the above Lagrangians, the scaling factors can be extracted as follows: where V represents W and Z, and I 3 f = 1/2 (−1/2) for f = u (d, e). For the loop induced couplings φγγ, φZγ and φgg, we define κ φ XY ≡ Γ(φ → XY )/Γ(φ → XY ) SM with Γ(φ → XY ) being the decay rate of φ → XY and XY = γγ, Zγ and gg. For the charged Higgs bosons H ± , their Yukawa couplings are expressed as where P L (P R ) is the projection operator for left-(right-) handed fermions and V ud is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element. Here, we also give the scalar trilinear couplings λ φφ φ defined by the coefficient of the corresponding Lagrangian term, which are relevant to the decay rates discussed in Sec. III: where Eq. (23) is followed from the CP-invariance.
Let us discuss the important limits of the parameters in the THDM, the decoupling limit and the alignment limit. First, the decoupling limit is realized by taking M → ∞, by which all the masses of the additional Higgs bosons become infinity, and only h remains at the EW scale 1 . In this limit, new physics effects on low energy observables disappear due to the decoupling theorem [99,103]. Second, on the other hand, the alignment limit can be defined by taking s β−α → 1, in which the h 1 state in Φ coincides with the mass eigenstate h, and κ h V = κ h f = 1 is satisfied at tree level. We note that this limit is automatically realized in the decoupling limit. Here, the important thing is that if κ h V = 1 and/or κ h f = 1 are found at future collider experiments, we cannot take the decoupling limit. This provides us a new no-loose theorem [28,87,95], where we can extract the upper bound on the mass scale of the second Higgs boson 2 . Quantitatively, such a bound is given by imposing the constraints from perturbative unitarity and vacuum stability as we will discuss them below. Notice here that the inverse of the above statement does not hold in general, namely, alignment without decoupling can be considered. Such scenario is well motivated by; e.g., the successful EW baryogenesis [62][63][64][65][66][67][68][69].
In this case, m 2 H and m 2 h are determined only by M 2 22 and M 2 11 , respectively. 2 An original no-loose theorem was discussed for the SM in Ref. [104].
As mentioned above, we take into account the perturbative unitarity and the vacuum stability bounds. For the unitarity bound, we impose |a i | ≤ 1/2, where a i are independent eigenvalues of the s-wave amplitude matrix for two-body to two-body scattering processes in the high-energy limit [94,[105][106][107]. The analytic expressions for a i are given in Appendix A.
In this limit, due to the equivalence theorem [108], only the contact scalar interaction terms contribute to the s-wave amplitude, which can be written in terms of the scalar quartic couplings. Thus, the unitarity bound gives constraints on the masses of additional Higgs bosons and the mixing angle through the relations given in Eqs. (5)-(12), see also Eqs. (A8)-(A11). On the other hand, the vacuum stability is the requirement that the Higgs potential is bounded from below in any direction with large field values. The sufficient and necessary conditions are given in Refs. [109][110][111][112]. We give a comment on the true vacuum condition.
The Higgs potential can have several extrema besides the EW true vacuum. In such a case, we need to ensure that the true vacuum is the deepest vacuum than all the other ones. In Ref. [113], it has been shown that most of the parameter regions with M 2 < 0 are excluded by the true vacuum condition. Thus, throughout the paper we simply assume M 2 to be a positive value in order to satisfy the true vacuum condition.
Before closing this section, we briefly mention constraints from various flavor observables which particularly sensitive to the mass of the charged Higgs bosons. The comprehensive studies of these constraints in the Z 2 symmetric THDMs have been carried out in Refs. [88,114]. In Type-II, the B → X s γ process gives the lower bound of m H ± 800 GeV at 95% confidence level (CL) almost independently of the value of tan β for tan β 2 [115]. On the other hand in Type-I, the severe constraint on m H ± is given particularly for smaller tan β; e.g., m H ± 450 GeV for tan β = 1 [116]. However, above tan β 2, the bound becomes weaker than the lower bound from the direct search at LEP; i.e., m H ± 80 GeV [76]. Because the lepton Yukawa couplings are irrelevant to the B → X s γ process, similar bounds given in Type-I and Type-II can be obtained in Type-X and Type-Y, respectively.
In Type-II, B → τ ν and B s → µμ processes give an upper limit on tan β; e.g., tan β 20 for m H ± = 800 GeV [88]. In Type-X, constraint by τ → µνν becomes important for large tan β [102,117,118]. In the small tan β region, the neutral meson mixing processes B 0 −B 0 give a stronger bounds for m H ± compared to the bound from B → X s γ, and these exclude the wide region in all the types of THDMs.

III. DECAYS OF THE HIGGS BOSONS
In this section, we give the analytic expressions for the decay rates of the Higgs bosons including higher-order corrections in QCD. In addition, some numerical results for the decays of the Higgs bosons are shown.

A. Running parameters
We give the expressions for the running strong coupling α s (µ) and the running quark masses m q (µ) at the scale µ in the MS scheme. In order to compute these variables, we need the coefficients of the β function for α s (µ) and those of the anomalous dimension for m q (µ).
Their formulae at the three-loop level are given by [119,120] with N f being the number of active flavors and with ζ(n) indicating the Riemann zeta function. The running strong coupling α s at the scale µ is expressed as where µ = ln(µ 2 /Λ 2 QCD ), and Λ QCD is the asymptotic scale parameter [121]. The running quark mass at the scale of the pole mass m q is given by [122][123][124][125] The running quark mass at the scale µ is expressed as where the function c(x) is given by [126,127] c(x) = xγ 0 1 + (γ 1 −β 1γ0 )x + 1 2 If the renormalization group (RG) evolution crosses the flavor threshold, we need to take into account the matching condition of the running strong coupling [128] and the running quark mass [129,130] α where ζ g and ζ m are the matching coefficients. We note that the matching coefficients are unity up to NLO, and we use ζ g = ζ m = 1 in the following. For example, the running charm quark mass at m h can be evaluated as

B. QCD corrections to the neutral Higgs decays
In the following, we describe how to include QCD corrections for processes of the neutral Higgs bosons φ (= h, H, A) in our calculations. For the decay rates of h, we adopt the formulae of incorporating those QCD corrections in H-COUP v2 [56].
The decay rate into a pair of light quarks (q = t) including next-to-next-to-leading order (NNLO) QCD corrections in the MS scheme is given by [120,[131][132][133][134] where with the color factor C F = 4/3. The last term ∆ φ t-loop indicates top-quark loop contributions, which calculated in the case with m t m φ and µ = m φ as In the LO decay rate Γ 0 , mass parameters arising from Yukawa couplings are replaced by the running masses m q (µ). Thereby, large logarithmic corrections induced by the light quark masses are resummed [135].
For the top pair, the QCD correction factor ∆ φ t depends on the CP property of the Higgs boson. We obtain the decay rate at the NLO in the on-shell scheme as where [61,136] where Li 2 is the dilog function. In the chiral limit β t → 1, we obtain Contributions of the top quark mass in the NLO QCD corrections are significant near the threshold region. On the other hands, dominant contributions in m φ m t can be the logarithmic contribution, ln(m 2 t /m 2 φ ), which appears in the QCD corrections in the MS scheme. In order to take into account both of the effects, we use interpolation for the corrections to φ → tt as discussed in Ref. [137].
For the decays into an off-shell gauge boson φ → V V * and φ → φV * (V = W, Z), the QCD correction can enter in the V * → qq part. This effect can be included by [138] where 13 The fermion loop contribution to the decay rate of φ → γγ receives QCD corrections.
At the NLO, the QCD correction can be implemented by the following replacement of the quark loop function I φ F (τ q ) in the MS scheme [27,139] where I φ F (τ q ) is defined in Appendix B, and the factor C φ is determined by the scale µ and the mass ratio τ q ≡ m 2 φ /(4m 2 q ). In our computation, we adopt the analytic expression of C φ given in Ref. [140], in which C φ is written in terms of the polylog functions, up to the Li 4 function. It has been known that the factor C φ becomes the simple form in the large top mass limit, τ t → 0, as [27,139,141] On the other hand, in the large Higgs mass limit or equivalently the massless fermion limit, the factor C φ is common to the case for the CP-even and CP-odd Higgs boson [27]: For H/A → Zγ decays, we calculate them at the LO.
For the φ → gg decays, we take into account the decay rate corrected up to NNLO expressed as, For the NLO QCD corrections to the φ → gg decays, there are contributions from virtual gluon loops and those from real emissions of a gluon (φ → ggg) and a gluon splitting into (49) can be decomposed as [27], The first and second terms respectively denote the contribution from virtual gluon loops and that from real gluon emissions in the large top-quark mass limit. These are expressed The last term ∆E φ vanishes in the large top-quark mass limit, which can be decomposed into the following three parts: Similar to the φ → γγ decays, we adopt the analytic expression for the virtual correction ∆E virt φ given in Ref. [140]. Those for the real emissions ∆E ggg φ and ∆E gqq φ are given in Ref. [27], which are expressed in the form with a double integral with respect to phase space variables. According to Ref. [27], the factor ∆E φ is dominantly determined by the contribution from the virtual gluon loop ∆E virt φ , so that in our computation we neglect the contributions from ∆E ggg φ and ∆E gqq φ . From Eq. (53), E φ is given to be about 18 at µ = m φ and N f = 5, and it gives sizable correction to the decay rate; e.g., ∼ 70% for m φ = 100 GeV. For NNLO contributions; i.e., E (2) φ , we incorporate those in the limit with m t m φ and setting as µ = m φ , which are expressed as [142,143]

C. QCD corrections to the charged Higgs decays
The QCD corrections to charged Higgs decays into light quarks are presented in the MS scheme. The expression can be written in the same way with the neutral Higgs boson decays as where the ∆ H ± q is given by Eq. (35) but without the last term ∆ φ t-loop . For the the decays into quarks including the top quark, we apply the QCD correction in the on-shell scheme.
It is given in [61,144] where where . A function B qq is given in Ref. [61]. In these expressions quark pole masses are used. Similar to φ → tt, we incorporate the corrections with interpolation to consider the effect of the top quark mass and the logarithmic corrections due to light down-type quark masses.
For the off-shell decays into a neutral Higgs boson and a W boson, H ± → φW * , the QCD correction can be applied as similar to φ → φ V * . It can be written as where the QCD correction factor is given in Eq. (45). For loop induced decay processes of the charged Higgs bosons, H ± → W ± V (V = Z, γ), which have been studied in Refs. [145][146][147][148][149][150], we calculate them at the LO.

D. Total decay widths and decay branching ratios
We here discuss total widths and branching ratios for the neutral Higgs bosons and the charged Higgs bosons in four types of the THDMs in order for later discussion about direct searches of heavy Higgs bosons. We describe the behavior of the total widths and the branching ratios in cases with the alignment limit, s β−α = 1 and without taking the alignment limit, s β−α = 0.995. In the numerical computations, we use the beta version of H-COUP v3 [151], where the QCD corrections presented in previous subsections are included. In the QCD correction functions C φ and E virt φ , polylog functions appear. We use CHAPLIN [152] for the numerical evaluation of such polylog functions. We have confirmed that our numerical results for the total widths and the branching ratios are consistent with We here show the case that masses of the additional Higgs bosons as well as M are GeV or 800 GeV, tan β is scanned in the following range, 0.5 < tan β < 50.
We note that, without depending on tan β, results with m Φ = 200 GeV for Type-II and Type-X are already excluded by the constraint from the flavor physics (also, for Type-I and Type-Y in lower tan β regions, tan β 2) [115,116]. Nevertheless, we show them in order to compare results among four types of the THDM. For the SM parameters, we use the  for N f = 6, 5, 4, and 3, respectively. The input value of the CKM matrix elements and the total width for the weak gauge bosons as well as the top quark are taken as [154], The former is relevant for the charged Higgs decays into quarks, H ± → tb, H ± → ts and H ± → cb. The latter is used in computation of the Higgs boson decays into off-shell particles.
Before we show numerical behaviors of the total widths and the branching ratios, we mention the loop induced decays of the charged Higgs bosons. The branching ratio of H ± → W ± Z can be enhanced when the mass difference between H ± and A is taken to some extent [146,150]. Whereas, in the following numerical results, where the additional Higgs bosons are degenerate, the branching ratio of H ± → W ± Z is at most O(10 −4 ) in the present parameter choices. Furthermore, the branching ratio of H ± → W ± γ is smaller than that of The following numerical results for the total widths and the branching ratios are similar to those given in Ref. [87], where the systematic studies have been done. Nevertheless, we here show them because there are some developments from the previous study. Main difference from Ref. [87] is that we compute the decay processes including higher-order QCD corrections. Also, we incorporate the above mentioned decay processes for the charged Higgs bosons, H ± → W ± Z and H ± → W ± γ, in the evaluation of the total width.
In Fig. 1 The branching ratios including QCD corrections are discussed in the above paragraphs.

IV. DIRECT SEARCHES AT THE LHC
In this section, we present current constraints on the parameter space in the THDMs from direct searches for heavy Higgs bosons with the LHC Run-II data.
Let us briefly summarize the procedure how we obtain the constraints on the parameters in the THDMs from model-independent analyses for heavy Higgs boson searches at the LHC. First, we compute production cross sections of heavy neutral Higgs bosons, φ = H and A, in the THDMs for the gluon-fusion process (pp → φ) and for the bottom-quark associated (or bottom-quark annihilate) process (pp → φ(bb)) at the NNLO in QCD by using Sushi-1.7.0 [162,163]. For the charged Higgs boson production pp → tH ± , we use the values given at the NLO QCD by the Higgs cross section working group (HXSWG) [164], based on Refs. [165][166][167][168]. Second, we calculate decay branching ratios of the Higgs bosons in the THDMs, including higher-order QCD corrections, as described in Sec. III. Finally, we compute the production cross sections times the branching ratios for each parameter point for each search channel at the LHC listed in Table II,  Although we use the ATLAS data, listed in Table II, the similar limits have been reported

Constrained quantity
Applicable mass region Reference Fig. 7 Fig. 9(a)

A. Production cross sections for the additional Higgs bosons
Before we discuss current constraints on the parameter space from direct searches, we present production rates for the heavy Higgs bosons at the 13 TeV LHC.  For the gluon-fusion process, shown in the left two columns in Fig. 6, the Higgs bosons are produced via quark loops. Therefore, the difference of the Yukawa sector between Type-I and Type-II in Eq. (19) leads to significantly different dependence on the model parameters.
In Type-I, where the top-quark loop is entirely dominant, the larger tan β is, the smaller the cross section is for a fixed mass. One can also see the threshold enhancement of the top-quark loop at m Φ ∼ 2m t . In Type-II, the top-quark loop is dominant for small tan β, small for small tan β. In the large tan β region, on the other hand, the cross sections for a fixed mass tend to be larger as s β−α deviates from the alignment limit. The production via the bottom-quark associated process, shown in the right two columns in Fig. 6, is entirely subdominant in Type-I, while that becomes dominant for large tan β in Type-II.
In Fig. 7, similar to Fig. 6, but for c β−α > 0, we show the production rates.  (left two columns) and the bottom-quark associated process (right two columns). Different from the CP-even Higgs bosons, the production rates only depend on tan β because of the Yukawa structure in Eq. (19). The global parameter dependence of the cross sections via the gluon fusion is similar to that for H with s β−α = 1, but the production rate for A is slightly larger than that for H at each point on the m Φ -tan β plane. The parameter dependence of the cross sections via the bottom-quark annihilation is as same as for H with s β−α = 1.
In Fig. 9, at the LHC charged Higgs bosons H ± are mainly produced in association with a top quark via gb → tH ± for m H ± > m t , whose cross sections are shown. Similar to the productions for A, the cross section only depends on tan β. For a fixed mass, in Type-I, the larger tan β is, the smaller the production rate is. In Type-II, on the other hand, up to tan β ∼ 7, the larger tan β is, the smaller the production rate is, similar to the Type-I case.
However, for tan β 7, the production rate becomes larger for larger tan β due to tan β enhancement of the bottom-Yukawa coupling.
We Regarding to the CP-odd Higgs boson A; • For large tan β, exclusion regions only appear in the Type-II and the Type-Y THDMs, in which the production via the bottom-quark loop as well as the bottom-quark associated production becomes dominant.
• The A → ττ channel is significant only for m A < 2m t or for large tan β in Type-II.
We note that, although the branching ratio of the A → ττ decay is even dominant for large tan β in Type-X, the production rate is too small to be constrained.    case. The region of the exclusion from the A → Zh channel becomes larger from s β−α = 0.995 to 0.98, since the decay rate for A → Zh is proportional to c 2 β−α .
Regarding to the CP-even heavier Higgs boson H; • The production rate via the gluon fusion for the heavier CP-even Higgs boson H is smaller than that for the A production, as mentioned above. Moreover, in the non-alignment case, the fermionic branching ratios of H for low tan β is smaller than those for A due to the decays into a pair of the weak gauge bosons, which are forbidden for A. Therefore, the constraints are slightly weaker than the A case, and we do not present the exclusions explicitly for the H → ττ , H(bb) → ττ , H(bb) → bb and H → tt channels.
• Regarding to the charged Higgs boson H ± ; • For the near alignment scenario, in the low tan β region (tan β 5), the H ± → tb decay is dominant for all the types, therefore the exclusions of the low-mass and lowtan β region from the H ± → tb channel are almost same for all the panels.
• In the large tan β region, the constraint from the H ± → τ ν channel can be significant only in Type-II. Although the branching ratio of the H ± → τ ν is even dominant for large tan β in Type-X, the constraint is insignificant due to the small production rate.
• We note that, as mentioned in Sec. II, in Type-II and Type-Y there is an independent constrain from flavor observables on the mass of charged Higgs bosons, m H ± 800 GeV. Figure 11 shows the same as in Fig. 10, but for the c β−α > 0 case. The global picture of the exclusion regions is same as for the c β−α < 0 case. A remarkable difference is that the constraints for H in the non-alignment case are much weaker for around tan β ∼ 7 − 10 due to the strong suppression of the production rates. Although σ(A → Zh) does not depend on the sign of c β−α , the exclusion regions for c β−α > 0 in Type-II and Y are smaller than those for c β−α < 0. This is because the analysis includes the h → bb decay, whose branching ratio has a singular behavior for c β−α > 0; see Figs. 4 and 5. Before closing this section, we briefly discuss the signal strength for the discovered Higgs boson measured at the LHC Run-II experiment, which provides independent constraints on the parameter space from those given by the direct searches discussed in this section.
Measurements of the signal strength set constraints on the Higgs boson couplings; i.e., the κ values defined in Sec. II, which can be translated into those on s β−α and tan β. In Table III, we summarize the 95% CL allowed range of tan β in the THDMs with fixed values of s β−α .
The κ values are extracted from Ref.
[1], which are presented in Table IV as a reference.
We see that except for the Type-I THDM it gives severe constraints on tan β, because κ b and/or κ τ can significantly differ from unity in the Type-II, Type-X and Type-Y THDMs even for the approximate alignment case.

V. COMBINED RESULTS OF DIRECT SEARCHES AT THE HL-LHC AND PRE-CISION TESTS AT THE ILC
Now, let us turn to investigate how the current parameter space in the THDMs discussed in the previous section can be explored further in future experiments, especially by direct searches for heavy Higgs bosons at the HL-LHC as well as by precision measurements of the Higgs boson couplings at the ILC. We note that complementarity for direct searches for heavy Higgs bosons between at the LHC and the ILC500 was discussed for the THDMs in Ref. [86].
In order to obtain the sensitivity projection to the HL-LHC with 3000 fb −1 of integrated luminosity, we rescale the current expected sensitivity by 3000/36 ∼ 9.1. We also perform a further rescaling of the sensitivity from √ s = 13 TeV to √ s = 14 TeV by taking into account the ratio of the signal cross sections, σ(m Φ ) 14TeV /σ(m Φ ) 13TeV . Here, we assume    [20]. For the ILC500, the expected accuracies are based on the results of the ILC250 combining the simulations at √ s = 350 GeV with 200 fb −1 and those at √ s = 500 GeV with 4000 fb −1 [20].
that signal and background increase by the same amount from 13 TeV to 14 TeV, which can be conservative particularly for the high-mass region. Detailed projection with systematic uncertainties for the φ → ττ channel was performed in the report for the HL-LHC [93], where one can see the higher sensitivity for m Φ 1200 GeV.
In addition, from precision measurements of the 125 GeV Higgs boson couplings, we can further constrain the parameter space in the THDMs. In Table IV  these uncertainties can be reduced significantly at those future collider experiments; e.g., κ Z is expected to be measured with a few percent at the HL-LHC and less than 1% at the ILC.
As we explained in Sec. II, if a nonzero deviation in a Higgs boson coupling is confirmed, an upper limit on the mass of the additional Higgs bosons can be given because the decoupling   limit is no longer realized. In the following discussion, we numerically derive the upper limit on the common mass of the additional Higgs bosons m Φ by imposing the bounds from perturbative unitarity and vacuum stability, which are discussed in Sec. II. We will see that the upper limit appears for the non-alignment case s β−α = 1, depending on the value of tan β.
In Fig. 12, we show regions on the m Φ -tan β plane expected to be excluded at 95% CL   in the Type-I, Type-II, Type-X and Type-Y THDMs (from the left For these constraints, we scan the value of M 2 with M 2 > 0, so that the black shaded region indicates that there is no value of M 2 which simultaneously satisfies the unitarity and the vacuum stability bounds. In the above sense, the black region can be regarded as a conservative excluded region. Interestingly, it is seen that a non-zero deviation for the 125 GeV Higgs couplings from the SM prediction sets an upper limit of the heavy Higgs masses. For s β−α = 0.995, the alignment limit is included by the 2σ error, so that the dashed curve does not appear. Details of the behavior of the upper limit from precision measurements on m Φ , shown in For c β−α < 0, the third condition of the vacuum stability bound given in Eq. (A7) sets an upper limit on M which is slightly smaller than m Φ almost without depending on the value of tan β; e.g., M 680, 730 and 780 GeV being excluded for x = −0.1 and m Φ = 800, 900 and 1000 GeV, respectively, where x ≡ π/2 − (β − α). The important point here is that the required value of m 2 Φ − M 2 (> 0) gets larger for a larger value of m Φ . Whereas, the unitarity bound excludes a larger difference between M 2 and m 2 Φ , which makes magnitudes of the λ parameters larger, as seen in Eqs. (A8)-(A11). Therefore, for a fixed value of s β−α and tan β we can find a critical value of m 2 Φ , above which the solution of the value of M 2 to satisfy the both unitarity and vacuum stability bounds vanishes. Such an upper limit on m Φ becomes stronger when the value of tan β differs from unity because the λ 1 or λ 2 parameter becomes significant so that the unitarity bound sets more severe constraint on We here emphasize that the entire parameter space we consider is explored by combining the constraints from the direct searches at the HL-LHC and from the precision measurements of the 125 GeV Higgs boson couplings at the ILC. Figure 13 shows the same as in Fig. 12, but for the c β−α > 0 case. Because of the singular behaviors of the production cross section for H and of the branching ratios for h around tan β ∼ 7 − 10, shown in Figs. 7 and 5, a narrow parameter region in the Type-II and the Type-Y models remains without any constraints from the direct searches even for low m Φ . Similar to Fig. 12, there appears an upper limit on m Φ by the constraints of unitarity and vacuum stability in Fig. 13. A remarkable difference, however, arises from the vacuum stability bound as compared with the case for c β−α < 0. In this case with a low tan β region, the condition λ 2 > 0 sets an upper limit on M 2 for a fixed value of m 2 Φ with M 2 m 2 Φ . This upper limit on M 2 gets milder when tan β becomes larger. When tan β exceeds a certain value, the upper limit on M 2 is almost fixed to be m 2 Φ due to the condition λ 1 > 0 instead of λ 2 > 0. Such a non-trivial tan β dependence on the vacuum stability bound provides two peaks of the upper limit on m Φ as seen in Fig. 13. As a result, some small parameter regions remain uncovered by both the HL-LHC and the ILC250.
We here give a comment on the case, where the degeneracy between the common mass of the additional Higgs bosons m Φ and M is relaxed. In the above analysis, we have set M = m Φ in the analysis of the exclusion region by the direct searches for simplicity. As we have mentioned in Sec. III, the decay width for H → hh depends on the value of M , and the exclusion region for H might change if we consider the case of M = m Φ . We note, however, that most of the parameter regions excluded by H → hh are also excluded by the A → Zh decay mode, which does not depend on the value of M . Therefore, our main conclusion does not change even if we relax the degeneracy among m Φ and M .
To summarize, the entire parameter space in the THDMs can be explored by the synergy between the direct searches at the HL-LHC and the precision measurements of the 125 GeV Higgs boson couplings at the ILC. In other words, if we observed any deviations for the Higgs boson couplings at the ILC, we would be able to find additional Higgs bosons at the HL-LHC, or reject a certain type of new physics models. In order to quantify the above statement, we have also checked the 5σ discovery sensitivity by naive rescaling. We find that the discovery regions are certainly smaller than the 95% CL excluded region shown in Figs. 12 and 13. Consequently, for c β−α < 0, we find that most of the parameter space is covered by the direct searches at the HL-LHC and the precision tests at the ILC250. For c β−α > 0, on the other hand, some parameter regions appear, which requires more data and/or more precision to be explored.

VI. CONCLUSIONS
We have discussed the possibility that a wide region of the parameter space in the four types of the THDMs can be explored by the combination of the direct searches for the additional Higgs bosons at the LHC and precision measurements of the discovered Higgs boson couplings at future lepton colliders. The direct searches give lower limits on the masses of the additional Higgs bosons, while the precision measurements set upper limits by using the perturbative unitarity and the vacuum stability bounds. Thus, these two searches play an complementary role to explore the parameter space. We first have shown that the parameter region excluded by the direct search at the LHC Run-II, and then shown that the exclusion expected by using the synergy between the direct searches at the HL-LHC and the precision tests assuming the accuracy expected for the measurements of the Higgs boson couplings at the ILC with a collision energy of 250 GeV. It has been found that in the  The unitarity bound is defined by |a i | ≤ 1/2 as we discuss in Sec. II, where the independent eigenvalues of the s-wave amplitude matrix are given by [94,[105][106][107] a ± 4 = 1 16π (λ 3 + 2λ 4 ± 3λ 5 ), (A4) a ± 5 = 1 16π (λ 3 ± λ 4 ), (A5) a ± 6 = 1 16π (λ 3 ± λ 5 ).
As we see the above expression, the unitarity and the vacuum stability bounds constrain the value of the λ parameters. Thus, it would be convenient to express these parameters in terms of the physical parameters as follows: where m Φ = m H = m A = m H ± , and x ≡ π/2 − (β − α) such that x = 0 corresponds to the alignment limit s β−α = 1.
Appendix B: Decay rates at the leading order We present the analytic expressions of the decay rates of the Higgs boson at the LO. In order to specify the LO formula, the subscript 0 is put in the decay rate, Γ 0 .
The decay rates into a pair of on-shell weak bosons (V = W, Z) are given by where c V = 1 (2) for V = W (Z). When one of the weak bosons is off-shell, we obtain where F (x) = −|1 − x 2 | 47 2 x 2 − 13 2 + 1 x 2 + 3(1 − 6x 2 + 4x 4 )| log x| The loop induced decay rates are given by where