Type-II 2HDM under the Precision Measurements at the $Z$-pole and a Higgs Factory

Future precision measurements of the Standard Model (SM) parameters at the proposed $Z$-factories and Higgs factories may have significant impacts on new physics beyond the Standard Model in the electroweak sector. We illustrate this by focusing on the Type-II two Higgs doublet model (Type-II 2HDM). The contributions from the heavy Higgs bosons at the tree-level and at the one-loop level are included in a full model parameter space. We perform a multiple variable global fit and study the extent to which the parameters of non-alignment and non-degenerate masses can be probed by the precision measurements. We find that the allowed parameter ranges are tightly constrained by the future Higgs precision measurements, especially for small and large values of $\tan\beta$. Indirect limits on the masses of heavy Higgs can be obtained, which can be complementary to the direct searches of the heavy Higgs bosons at hadron colliders. We also find that the expected accuracies at the $Z$-pole and at a Higgs factory are quite complementary in constraining mass splittings of heavy Higgs bosons. The typical results are $|\cos(\beta-\alpha)|<0.008, |\Delta m_\Phi |<200\ {\rm GeV}$, and $\tan\beta \sim 0.2 - 5$. The reaches from CEPC, FCC-ee and ILC are also compared, for both Higgs and $Z$-pole precision measurements.


Introduction
With the milestone discovery of the Higgs boson (h) at the CERN Large Hadron Collider (LHC) [1,2], particle physics has entered a new era. All the indications from the current measurements seem to confirm the validity of the Standard Model (SM) up to the electroweak (EW) scale of a few hundred GeV, and the observed Higgs boson is SM-like. Yet, there are compelling arguments, both from theoretical and observational points of view, in favor of the existence of new physics beyond the Standard Model (BSM) [3]. As such, searching for new Higgs bosons would be of high priority since they are present in many extensions of BSM theories. One of the most straightforward, but well-motivated extensions is the two Higgs doublet model (2HDM) [4], in which there are five massive spin-zero states in the spectrum (h, H, A, H ± ) after the electroweak symmetry breaking (EWSB). Extensive searches for BSM Higgs bosons have been actively carried out, especially in the LHC experiments [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Unfortunately, no signal observation has been reported thus far. This would imply either the non-SM Higgs bosons are much heavier and essentially decoupled from the SM, or their interactions are accidentally aligned with the SM configuration [19,20]. In either situation, it would be challenging to observe those states in experiments.
Complementary to the direct searches, precision measurements of SM parameters, in particular, the Higgs boson properties could lead to relevant insights into new physics. There have been proposals to build a Higgs factory in the pursuit of precision Higgs measurements, including the Circular Electron Positron Collider (CEPC) in China [21,22], the electron-positron stage of the Future Circular Collider (FCC-ee) at CERN (previously known as TLEP [23][24][25]), and the International Linear Collider (ILC) in Japan [26]. With about 10 6 Higgs bosons produced at the Higgs factory, one would expect to reach sub-percentage precision determination of the Higgs properties, and thus to be sensitive to new physics associated with the Higgs boson. As an integrated part of the program, one would like to return to the Z-pole. With about 10 10 − 10 12 Z bosons, the achievable precisions on the SM parameters could be improved by a factor of 20 − 200 over the Large Electron Positron (LEP) Collider results [27]. Such a high precision would hopefully shed light on new physics associated with the electroweak sector.
In this paper, we set out to examine the impacts from the precision measurements of the SM parameters at the proposed Z-factories and Higgs factories on the extended Higgs sector. There is a plethora of articles in the literature to study the effects of the heavy Higgs states on the SM observables [4]. We illustrate this by focusing on the Type-II 2HDM 1 . In our analyses, we include the tree-level corrections to the SM-like Higgs couplings and one-loop level contributions from the heavy Higgs bosons. A global fit is performed in the full modelparameter space. In particular, we study the extent to which the parametric deviations from the alignment and degenerate mass limits can be probed by the precision measurements. We find that the expected accuracies at the Z-pole and at a Higgs factory are quite complementary in constraining mass splittings of heavy Higgs bosons. The reach in the heavy Higgs masses and couplings can be complementary to the direct searches of the heavy Higgs bosons at the LHC.
The rest of the paper is organized as follows. In Section 2, we summarize the anticipated accuracies on determining the EW observables at the Z-pole and Higgs factories. Those expectations serve as the inputs for the following studies for BSM Higgs sector. We then present the Type-II 2HDM and the one-loop corrections, as well as the existing constraints to the model parameters in Section 3. Section 4 shows our main results from the global fit, for the cases of mass degeneracy and non-degeneracy of heavy Higgs bosons. We summarize our results and draw conclusions in Section 5.

The EW and Higgs Precision Measurements at Future Lepton Colliders
The EW precision measurements are not only important in understanding the SM physics, but also can impose strong constraints on new physics models [30,31]. The benchmark scenarios of several proposed future e + e − machines and the projected precisions on Z-pole and Higgs measurements are summarized below. These expected results serve as the inputs for the later studies in constraining the BSM Higgs sector.

The electroweak precision measurements
The current best precision measurements for Z-pole physics came mostly from the LEP-I, and partially from the Tevatron and the LHC [32,33]. These measurements could be significantly improved by a Z-pole run at future lepton colliders with a much larger data sample [21,[23][24][25]34]. For example, the parameter sin 2 θ ef f can be improved by more than one order of magnitude at the future e + e − collider; the Z-mass precision can be measured four times better in CEPC. Precisions of other observables, including m W , m t , m h , A b,c,l F B , R b , etc., can be improved as well, depending on different machine parameter choices. Given the complexity of a full Z-pole precision fit, we study the implications of Z-pole precision measurements on the 2HDM adopting the Peskin-Takeuchi oblique parameters S, T and U [35].
The anticipated precisions on the measurements of α s , ∆α [32,[36][37][38][39] for various benchmark scenarios of future Z-factories with the indicated Z data samples. The corresponding constrained S, T and U ranges and the error correlation matrices are listed in Tab. 2. The results listed as "current" are obtained directly from the Gfitter results which use the current Z-pole precision measurements [32,33], with reference values of the SM Higgs boson mass of m h ,ref = 125 GeV and m t ,ref = 172.5 GeV [33]. The predictions for future colliders are obtained by using the Gfitter package [32] with corresponding precisions for different machines, using the best-fit SM point with the current precision measurements as the central value. For the Z-pole observables with estimated precisions not yet available at future colliders, the current precisions Current (1.7 × 10 7 Z's) CEPC (10 10 Z's) FCC-ee (7 × 10 11 Z's) ILC (10 9 Z's)  Table 2. Estimated S, T , and U ranges and correlation matrices ρ ij from Z-pole precision measurements of the current results, mostly from LEP-I [27], and at future lepton colliders CEPC [21], FCC-ee [23] and ILC [34]. Gfitter package [32] is used in obtaining those constraints.
are used instead. As seen from the at 1σ level. FCC-ee would further improve the accuracy. In our analyses as detailed in a later section, the 95% C.L. S, T and U contours are adopted to constrain the 2HDM parameter spaces, using the χ 2 -fit with error-correlation matrices .

Higgs precision measurements
At a future e + e − collider of the Higgs factory with the center-of-mass energy of 240−250 GeV, the dominant channel to measure the Higgs boson properties is the Higgsstrahlung process of Due to the clean experimental environment and well-determined kinematics at the lepton colliders, both the inclusive cross section σ(hZ) independent of the Higgs decays, and the exclusive ones of different Higgs decays in terms of σ(hZ)×BR, can be measured to remarkable precisions. The invisible decay width of the Higgs boson can also be very well constrained. In addition, the cross sections of W W, ZZ fusion processes for the Higgs boson production grow with the center-of-mass energy logarithmically. While their rates are still rather small and are not very useful at 240−250 GeV, at higher energies in particular for a linear collider, such fusion processes become significantly more important and can provide crucial complementary information. For √ s > 500 GeV, tth production can also be used as well. To set up the baseline of our study, we hereby list the running scenarios of various machines in terms of their center-of-mass energies and the corresponding integrated luminosities, as well as the estimated precisions of relevant Higgs boson measurements that are used in our global analyses in Tab. 3. The anticipated accuracies for CEPC and FCC-ee are comparable for most channels, except for h → γγ. There are several factors that contribute to the difference for this channel, which include the superior resolution of the CMS-like electromagnetic calorimeter that was used in FCC-ee analyses, and the absence of background from beamstrahlung photons [23]. In our global fit to the Higgs boson measurements, we only include the rate information for the Higgsstrahlung Zh and the W W fusion process. Some other measurements, such as the angular distributions, the diboson process e + e − → W W , can provide important information in addition to the rate measurements alone [41][42][43].  Table 3. Estimated statistical precisions for Higgs boson measurements obtained at the proposed CEPC program with 5 ab −1 integrated luminosity [21], FCC-ee program with 5 ab −1 integrated luminosity [23], and ILC with various center-of-mass energies [40].
3 Type-II Two Higgs Doublet Model
The 2HDM Lagrangian for the Higgs sector can be written as with the Higgs potential of by assuming CP-conserving and a soft Z 2 symmetry breaking term m 2 12 .
After EWSB, one of the four neutral components and two of the four charged components are eaten by the SM gauge bosons Z, W ± , providing their masses. The remaining physical mass eigenstates are two CP-even neutral Higgs bosons h and H, with m h < m H , one CP-odd neutral Higgs boson A, as well as a pair of charged ones H ± . Instead of the eight parameters appearing in the Higgs potential m 2 11 , m 2 22 , m 2 12 , λ 1,2,3,4,5 , a more convenient choice of the parameters is v, tan β, α, m h , m H , m A , m H ± , m 2 12 , where α is the rotation angle diagonalizing the CP-even Higgs mass matrix 2 .
The Type-II 2HDM is characterized by the choice of the Yukawa couplings to the SM fermions and is given in the form of After EWSB, the effective Lagrangian for the light CP-even Higgs couplings to the SM particles can be parameterized as (3.6) for i indicates individual Higgs coupling. Their values at the tree level are Our sign convention is β ∈ (0, π 2 ), β − α ∈ [0, π], so that sin(β − α) ≥ 0. The CP-even Higgs couplings to the SM gauge bosons are g hV V ∝ sin(β − α), and g HV V ∝ cos(β − α). The current measurements of the Higgs boson properties from the LHC are consistent with the SM Higgs boson interpretation. There are two well-known limits in 2HDM that would lead to a SM-like Higgs sector. The first situation is the alignment limit [19,45] of cos(β − α) = 0, in which the light CP-even Higgs boson couplings are identical to the SM ones, regardless of the other scalar masses, potentially leading to rich BSM physics. For sin(β − α) = 0, the opposite situation occurs with the heavy H being identified as the SM Higgs boson. While it is still a viable option for the heavy Higgs boson being the observed 125 GeV SM-like Higgs boson [46,47], the allowed parameter space is being squeezed with the tight direct and indirect experimental constraints. Therefore, in our analyses below, we identify the light CP-even Higgs h as the SM-like Higgs with m h fixed to be 125 GeV. The other well-known case is the "decoupling limit", in which the heavy mass scales are all large m A,H,H ± 2m Z [48], so that they decouple from the low energy spectrum. For masses of heavy Higgs bosons much larger than λ i v 2 , cos(β − α) ∼ O(m 2 Z /m 2 A ) under perturbativity and unitarity requirement. Therefore, the light CP-even Higgs boson h is again SM-like. Although it is easier and natural to achieve the decoupling limit by sending all the other mass scales to be heavy, there would be little BSM observable effects given the nearly inaccessible heavy mass scales. We will thus primarily focus on the alignment limit.
Note that while κ g , κ γ and κ Zγ are zero at the tree-level for both the SM and 2HDM, they are generated at the loop-level. In the SM, κ g , κ γ and κ Zγ all receive contributions from fermions (mostly top quark) running in the loop, while κ γ and κ Zγ receive contribution from W -loop in addition [49]. In 2HDM, the corresponding hf f and hW W couplings that enter the loop corrections need to be modified to the corresponding 2HDM values. Expressions for the dependence of κ g , κ γ and κ Zγ on κ V and κ f can be found in Ref. [50]. There are, in addition, loop corrections to κ g , κ γ and κ Zγ from extra Higgs bosons in 2HDM.
It is of particular importance to include a discussion for the triple couplings among Higgs bosons themselves. At the alignment limit, which is the parameter that enters the Higgs self-couplings and relevant for the loop corrections to the SM-like Higgs boson couplings. This parameter could be used interchangeably with m 2 12 as we will do for convenience. For the rest of our analysis, we fix v = 246 GeV and m h = 125 GeV. The remaining free parameters are tan β , cos(β − α) , m H , m A , m H ± and λ . (3.10) Note that while these six parameters are independent of each other, their allowed ranges under perturbativity, unitarity, and stability consideration are correlated. For simplicity with important consequences, one often starts from the degenerate case where all heavy Higgs boson masses are set the same. We will explore both the degenerate and non-degenerate cases specified as Given the current LHC Higgs boson measurements [51][52][53][54], deviations of the Higgs boson couplings from the decoupling and alignment limits are still allowed at about 10% level. All the tree-level deviations from the SM Higgs boson couplings are parametrized by only two parameters: tan β and cos(β −α). Once additional loop corrections are included, dependences on the heavy Higgs boson masses as well as λv 2 also enter. In our analyses below, we study the combined contributions to the couplings of the SM-like Higgs boson with both tree-level and loop corrections.
Before concluding this section, a special remark is in order. The model parameters introduced in this section and henceforth are all at the electroweak scale, identified as on-shell parameters to directly compare with experimental measurements. We do not consider the running effects due to other new physics at a higher scale such as in Supersymmetry or Grand Unified theories. This would become relevant if one asks whether the alignment behavior could be a natural result due to some symmetry or other principles [20]. In such scenarios, the alignment may take place at a higher scale but could be modified at the electroweak scale. Our results here, on the other hand, could be viewed as the acceptable deviations from the exact alignment conditions in a more fundamental theory.

Loop corrections to the SM-like Higgs couplings
We define the normalized SM-like Higgs boson couplings including loop effects as where κ tree ≡ g 2HDM tree /g SM tree . g 2HDM loop (Φ) and g 2HDM loop (SM) are the 2HDM Higgs boson couplings including loop corrections with heavy Higgs bosons or with SM particles only, respectively.
To the leading order in 1-loop corrections, Eq. (3.13) simplifies to (3.14) In the alignment limit of κ tree = 1, the term in the bracket is exactly zero, and κ 2HDM 1−loop | alignment = 1 + ∆κ 2HDM 1−loop . In our calculations, we adopt the on-shell renormalization scheme [55]. The conventions for the renormalization constants and the renormalization conditions are mostly following Refs. [55,56]. All related counter terms, renormalization constants and renormalization conditions are implemented according to the on-shell scheme and incorporated into model files of FeynArts [57] 3 . One-loop corrections are generated using FeynArts and FormCalc [63] including all possible one-loop diagrams. FeynCalc [64,65] is also used to simplify the analytical expressions. LoopTool [66] is used to evaluate the numerical value of all the loopinduced amplitude. The numerical results have been cross-checked with another numerical program H-COUP [67] in some cases.
For the couplings of the SM-like Higgs boson to a pair of gauge bosons and fermions, the general renormalized hf f and hV V vertices take the following formŝ where q µ , p µ 1 , and p µ 2 are the momenta of the Higgs boson and two other particles, respectively, and q 2 is the typical momentum transfer of the order m 2 h . κ i for each vertex is given byΓ S hf f andΓ 1 hV V for hf f and hV V , which includes both the tree-level and one-loop corrections: (3.17)

Loop corrections to Z-pole precision observables
The 2HDM contributions to the Peskin-Takeuchi oblique parameters [35] are given by [68] where we explicitly split these expressions into terms independent of or dependent on the alignment parameter of cos(β − α). The expression for various B and F -functions can be found in Ref. [68]. The mass splittings among heavy Higgs bosons of (m H , m A , m H ± ) violate the SU(2) custodial symmetry and thus will lead to contributions to the T and U parameters. In Fig. 1, we show the contributions to ∆S (left panel) and ∆T (right panel) in 2HDM varying ∆m A ≡ m A − m H and ∆m C ≡ m H ± − m H between ± 300 GeV, for cos(β − α) = 0. While the contribution to ∆S is typically small |∆S| 0.03, the contribution to ∆T quickly increases when m H ± is non-degenerate with either m A or m H . Therefore, an improved determination of ∆T from Z-pole precision measurement would severely constrain the mass splitting between the charged Higgs and its neutral partners. Furthermore, non-alignment case also breaks the symmetric pattern between ∆m A and ∆m C for ∆T contribution, preferring a slightly negative value of mass splittings.

Theoretical constraints and current experimental bounds
Heavy Higgs loop corrections would involve the Higgs boson masses and self-couplings, which are constrained by various theoretical considerations and experimental measurements, such as vacuum stability, perturbativity and unitarity, as well as electroweak precision measurements, flavor physics constraints, and LHC direct searches. We briefly summarize below the theoretical considerations and experimental constraints.

• Vacuum stability
In order to have a stable vacuum, the following conditions on the quartic couplings need to be satisfied [69]: (3.21) • Perturbativity and unitarity We adopt a general perturbativity condition of |λ i | ≤ 4π and the tree-level unitarity of the scattering matrix in the 2HDM scalar sector [70].
In Fig. 2, we show the constraints in the λv 2 -tan β plane once all the theoretical considerations are taken into account. For the upper panels, we work under the assumption with degenerate heavy Higgs boson masses m H ± = m H = m A ≡ m Φ . The left panel is for m Φ = 800 GeV and the right one is for m Φ = 2000 GeV, with cos(β − α) =0.005 (red curves), 0 (alignment limit, blue curves), and −0.005 (green curves). Regions enclosed by the curves are theoretically preferred. For a lower mass m Φ = 800 GeV, the constraints vary very little with the values of cos(β − α). The largest range on λv 2 ≡ m 2 Φ − m 2 12 / sin β cos β occurs at tan β = 1 [28]: which gives −0.29 < λ = −λ 4 = −λ 5 < 5.95 and 0 < λ 3 < 6.21. For a large value of m Φ = 2000 GeV, a slight shift of cos(β − α) leads to notable change in constraints on λv 2 , as shown by the red and green curves in the top right panel of Fig. 2.
The theoretically preferred region also depends on the individual heavy Higgs boson masses, as well as the deviation from the degenerate condition. In the lower panels of Fig. 2 The strongest bounds at large tan β come from A/H → τ + τ − mode, which excludes m A/H ∼ 300 − 500 GeV for tan β ∼ 10, and about 1500 GeV for tan β ∼ 50. The strongest bounds at small tan β 1 come from A/H → tt mode. The latest ATLAS search on such channel utilized the lineshape of tt invariant mass distribution, which exhibits a peak-dip structure due to the interference between the signal and the SM tt background [72,73]. A strong 95% C.L. bound of m A/H around 600 GeV can be reached for tan β = 1 for degenerate mass of m A = m H under the alignment limit. The direct searches for heavy charged Higgs bosons have been conducted with the H ± → (τ ν , tb) channels [16][17][18], and the bounds are relatively weak given the rather small leading production cross section for bg → tH ± , the large SM backgrounds for the dominant H ± → tb channel and the relatively small branching fraction of H ± → τ ν [74].  the alignment limit and mass-degenerate assumption. The strongest constraints for the large tan β region come from the A/H → τ + τ − searches: m A/H could be excluded to about 1000 GeV for tan β ∼ 10, and even larger masses for larger tan β. H ± → tb offers better exclusion at low tan β, which excludes m H ± to about 600 GeV for tan β ∼ 1. Possible A/H → tt mode might help to extend the exclusion reach to about 2000 GeV for tan β ∼ 1 [73,76]. At 100 TeV pp collider with 3 ab −1 luminosity, A/H → τ + τ − could extend the reach at large tan β to about 2000 GeV at tan β ∼ 10 and about 3 TeV for tan β ∼ 50. The coverage at low tan β could also be extended to about m H ± ∼ 1500 GeV via H ± → tb and m A ∼ 2500 GeV via A/H → tt for tan β ∼ 1 [75].
Since the branching fractions of the conventional search channels could be highly suppressed once other exotic decay channels of the non-SM Higgs boson to light Higgs bosons and/or SM gauge bosons open up [77][78][79], it is important to note that the current exclusion limits could be relaxed. Current LHC limits on m A,H via searches of exotic decay modes A/H → HZ/AZ are up to about 700 − 800 GeV, depending on the spectrum of non-SM Higgs bosons [12,13]. m A,H could be excluded to about 1500 GeV at HL-LHC and about 3000 GeV at 100 TeV pp collider [80].
While the exotic Higgs decay channel of A → h(→ bb, τ + τ − )Z is absent in the alignment limit, this channel could be used to constrain cos(β − α) and tan β when the deviation from the alignment limit is allowed. The projected A → hZ search results in the cos(β − α)-tan β plane of LHC 13 TeV with an integrated luminosity of 36 fb −1 (cyan) [11] and future HL-LHC 14 TeV with an integrated luminosity of 3 ab −1 (green) [81] for m A = 800 GeV (left panel) and m A = 2000 GeV (right panel) are shown in Fig. 3 with the colored survival regions. For the case of m A = 800 GeV, a narrow band within | cos(β − α)| 0.1 or | cos(β − α)| 0.02 is still allowed by the current LHC or the future HL-LHC data, as expected. Another branch from cos(β − α) = 0 to cos(β − α) = 1.0 with tan β decreasing from 5 − 10 to ∼ 0.1 is also allowed, which corresponds to the region with a suppressed BR(h → bb). The constraint for the m A = 2000 GeV case is far less stringent for the LHC 13 TeV case. Only the lower left region is excluded, in which both the production cross section σ(gg → A) and decay branching fraction of BR(A → hZ) × BR(h → bb) are enhanced. For the HL-LHC case, the tan β 1 regions are largely excluded, leaving the narrow band with | cos(β − α)| 0.1 or a branch stretching from cos(β − α) = 0 to cos(β − α) = 1.0 with tan β decreasing from ∼ 1 to ∼ 0.1 allowed by the future HL-LHC data. This is complementary to the SM-like Higgs boson signal strength measurements, which constrain the range of cos(β − α) to be less than about 0.1 around tan β ∼ 1 and even narrower regions for small and large tan β for Type-II 2HDM [28] with the current LHC measurements, except for a small wrong-sign Yukawa coupling region at tan β 2.
Flavor physics consideration usually constrains the charged Higgs mass to be larger than about 600 GeV for the Type-II 2HDM [74]. However, given the uncertainties involved in those flavor measurements, and that they are in general less stringent than the direct collider limits, we thus will not pursue the flavor bounds further.

Study Strategy and Results
In an earlier work [28], constraints from the tree-level effects on cos(β − α) and tan β, as well as from loop contributions in the degenerate mass case m H = m A = m H ± = m Φ under the alignment limit are analyzed. In this work, we extend the studies to more general cases of the non-degenerate masses and non-alignment, as well as including both the tree-level and one-loop contributions. We also incorporate the Z-pole precision results to show the complementarity between the Higgs and Z-pole precision measurements.

Global fit framework
To transfer the anticipated accuracy on the experimental measurements to the constraints on the model parameters, we perform a global fit by constructing the χ 2 with the profile likelihood method (4.1) Here, µ BSM i = (σ × BR) BSM /(σ × BR) SM for various Higgs search channels. We note that the correlations among different σ × BR are usually not provided, and are thus assumed to be zero in the fits. µ BSM i is predicted in each specific model, depending on model parameters. In our analyses, for the future colliders, µ obs i are set to be the SM value µ obs i = 1, assuming no deviations from the SM observables. The corresponding σ µ i are the estimated errors for each process, as already shown in Tab. 3 for the CEPC, FCC-ee and ILC. For the ILC with three different center-of-mass energies, we sum the contributions from each individual channel.
We fit directly to the signal strength µ i , instead of the effective couplings κ i . The latter are usually presented in most experimental papers. While using the κ-framework is easy to map to specific models, unlike µ i , various κ i are not independent experimental observables. Ultimately, fitting to either µ i or κ i should give the same results, if the correlations between κ i are properly included. Those correlation matrices, however, are typically not provided from experiments. Therefore, fitting to κ i only, assuming no correlations, usually leads to more relaxed constraints. For a comparison of µ-fit versus κ-fit results, see Ref. [28].
For Z-pole precision measurements, we fit into the oblique parameters S, T and U , including the correlations between those oblique parameters, as given in Tab. 2. We define the χ 2 as with X i = (∆S , ∆T , ∆U ) 2HDM being the 2HDM predicted values, andX i = (∆S , ∆T , ∆U ) being the current best-fit central value for current measurements, and 0 for future measurements. The σ ij are the error matrix, σ 2 ij ≡ σ i ρ ij σ j with σ i and correlation matrix ρ ij given in Tab. 2.

Case with degenerate heavy Higgs boson masses
We first consider the simple case of degenerate heavy Higgs boson masses m H = m A = m H ± ≡ m Φ such that the Z-pole precision are automatically satisfied. As shown in Ref. [28], in the Type-II 2HDM, the current LHC Higgs precision has already constrained cos(β −α) to be less than about 0.1. To explore the impact from the anticipated precision Higgs measurements at the CEPC, we perform a two-parameter global fit including the loop contributions. In Fig. 4, we show the 95% C.L. allowed region in the two-parameter cos(β − α)-tan β plane from the individual couplings by the colored curves: blue (κ b ), orange (κ c ), purple (κ τ ), green (κ Z ), cyan (κ g ), for a benchmark point of m Φ = 800 GeV, √ λv 2 = 300 GeV. κ γ does not have a notable effect therefore not shown. For large values of tan β, regions below the colored curves are allowed, while for small values of tan β, regions above the colored curves are allowed. The central red region is the global fit result with the best-fit point indicated by the black star. The two solid horizontal black lines represent the upper and lower limit for parameter tan β from theoretical constraints, as shown in Fig. 2

earlier. The region enclosed by the dashed black lines shows the tree-level only result for comparison.
For the Type-II 2HDM, the cos(β − α) region gets smaller for larger and smaller values of tan β. At large tan β, κ b and κ τ provide the strongest constraint since they are enhanced by a universal tan β factor. For small values of tan β, κ g (or effectively, κ t ) rules out large values of cos(β − α), followed by κ c for negative cos(β − α). Combining all the channels, the 95% C.L. region for the global fit leads to 0.2 ≤ tan β ≤ 30, −0.01 ≤ cos(β − α) ≤ 0.008, for √ λv 2 = 300 GeV is used here. The constraints from individual couplings are given with the color codes: blue (κ b ), orange (κ c ), purple (κ τ ), green (κ Z ), cyan (κ g ). The region enclosed by the dashed black lines shows the tree-level twoparameter global fit result for comparison. Two solid horizontal black lines represent the upper and lower limit for parameter tan β from theoretical constraints. the benchmark point m Φ = 800 GeV, √ λv 2 = 300 GeV. We note that the upper bound on tan β and the lower (negative) bound on cos(β − α) coming from κ g is mainly due to the large contribution from b-quark loop with a enhanced κ b . The overall range is slightly smaller than that obtained from the tree-level only result, shown by region enclosed by the dashed lines. The distorted shape of the global fit results, comparing to the tree-level only results is due to the interplay between both the tree-level contribution and loop corrections. Note that while κ Z can be measured with less than 0.2% precision, it is less constraining comparing to other couplings given the 1/ tan β (tan β) enhanced sensitivities for κ t,c (κ b,τ ) at small (large) tan β region.
To illustrate the dependence on m Φ and λv 2 , which enter the loop corrections, in Fig. 5, we show the 95% C.L. allowed region in the cos(β − α)-tan β plane given CEPC Higgs precision, for m Φ = 800 GeV, 3 TeV, the one-loop level effects almost decouple and the final allowed region is close to the tree-level results. Comparing with the constraints on the cos(β − α)-tan β plane via LHC searches with A → hZ channel as shown in Fig. 3, and the current and HL-LHC Higgs coupling precision measurements [28], the future Higgs factory can constrain the 2HDM parameter space at least an order of magnitude better in the allowed cos(β − α) range.
High precision on the Higgs coupling measurements can also be used to constrain the mass of the heavy Higgs bosons running in the loop. In Fig. 6, we show the 95% C.L. allowed region in the m Φ -tan β plane for  We also show the allowed regions in the m Φ -tan β plane under theoretical considerations in Fig. 6 with the different colors for different choices of cos(β −α). While all ranges of m Φ and tan β are allowed in the alignment limit of cos(β − α) = 0, once cos(β − α) deviates away from 0, large m Φ as well as small and large tan β regions are ruled out by theoretical considerations. Combining both the theoretical constraints and precision Higgs measurements, a constrained region in m Φ -tan β can be obtained for the non-alignment cases.
For √ λv 2 = 300 GeV, larger loop corrections further modify the allowed region in m Φ and tan β. The tt threshold region m Φ ≈ 350 GeV is inaccessible and the range of tan β is shrunk to 0.3 − 1.5 when cos(β − α) varies from 0 to 0.005. For the negative cos(β − α) = −0.005, the allowed region divides to two parts. The part with m Φ ≤ 1000 GeV has a wide range for parameter tan β, while for m Φ > 1000 GeV, 0.4 < tan β < 1.6. Theoretical considerations further limit the range of tan β to be between 0.35 and 3, as shown by the shaded region. For cos(β − α) = ±0.005, m Φ has an upper limit of about 2750 GeV from theoretical considerations. While λv 2 ≡ m 2 Φ − m 2 12 /s β c β is a good parameter to use since it is directly linked to the triple Higgs self-couplings, sometime it is convenient to fix the soft Z 2 breaking parameter m 2 12 instead. The resulting 95% C.L. allowed region in the m Φ -tan β plane is shown in Fig. 7 for m 12 = 0 (left panel) and 300 GeV (right panel). The theoretical constraints as discussed in the previous section are also indicated with the shaded gray regions. They have little dependence on the cos(β − α) value when m 2 12 is kept fixed. For m 12 = 0, m Φ = √ λv 2 is constrained to be less than around 250 GeV. For larger values of m 12 , the rather narrow region in the plane as seen in the right panel indicates a strong correlation between m Φ and tan β for large tan β, approximately scaled as tan β ∼ (m Φ /m 12 ) 2 , which minimizes the corresponding λv 2 value and thus its loop effects. The indirect probe in m Φ via Higgs precision measurements complements the direct search limits at the LHC, especially in the intermediate tan β wedge region where the direct search limits are the most relaxed.

Case with non-degenerate heavy Higgs boson masses
Going beyond the degenerate case, both the Higgs and Z-pole precision observables are sensitive to the mass splittings between the non-SM heavy Higgs bosons. In Fig. 8,   For the case of m H = m H ± (lower panels), the allowed range of ∆m Φ is larger, up to about 400 GeV for √ λv 2 = 0, and up to about 500 GeV for √ λv 2 = 300 GeV. Note that the region for ∆m Φ = 0 corresponds to the situation of cos(β − α) = 0 in Fig. 6, which is much less restrictive than the non-degenerate case ∆m Φ = 0.
In Fig. 9 Figure 9. Constraints on the ∆m A -∆m C plane from individual Higgs coupling measurement (color curves), and the 95% C.L. global fit results (red shaded region), for tan β = 0.2(left), 1 (middle), tan β = 7 (right) under alignment limit, with m H = 800 GeV, √ λv 2 = 300 GeV. For individual coupling constraint, the dashed line represents negative limit, while solid line represents the positive limit. Regions between the solid and dashed curves are the allowed region. For κ γ , region above the line is allowed.
plane from individual Higgs coupling measurements in color curves, and the 95% C.L. global fit results in the red shaded region, for tan β = 0.2 (left panel), 1 (middle panel) and 7 (right panel) under alignment limit with m H = 800 GeV, √ λv 2 = 300 GeV. For each individual coupling constraint with a "±" error bar, the dashed line is for the negative limit, while the solid line is for the positive limit. The range between the two lines is the survival region. Under the alignment limit, κ Z is independent of tan β as apparent in the figure. For Type-II 2HDM, generally speaking, κ b,τ are tan β-enhanced, while κ c is cot β-enhanced. Thus for small tan β, the main constraint on the mass splitting comes from κ c and leads to a small overlapping red region with κ Z as the global fit result of ∆m A ∼ −40 GeV to 0 GeV (left panel). For large tan β, it is due to κ b,τ , resulting in ∆m A ∼ −50 GeV to −250 GeV (right panel). For tan β ∼ 1, constraints from both κ b,τ and κ c are relatively relaxed, leading to a larger allowed region in the mass splittings ∆m A ∼ −250 GeV to 400 GeV (middle panel) mostly due to κ Z . The range of ∆m C is typically between −200 GeV to 100 GeV constrained from κ Z . κ γ mainly involves the charged Higgs loops and only constrains weakly. Note that κ g does not constrain the mass splittings significantly and therefore is not shown in the plots.
In Fig. 10, we present the 95% C.L. allowed region in the ∆m A -∆m C plane, for m H = 800 GeV (left panels) and 2000 GeV (right panels), again under the alignment limit. The upper panels are for 2000 GeV, the allowed ranges of the mass difference are much more relaxed and are almost independent of tan β. For √ λv 2 = 300 GeV, however, the largest ranges for ∆m C,A could be achieved for tan β ∼ 2, for both benchmark choices of m Φ , due to the constraints from individual couplings, as illustrated in Fig. 9. For m H = 2000 GeV, the allowed ranges of the mass difference varies little with 0.5 < tan β < 2, but shrink quickly for larger tan β.  In Fig. 11, we show the 95% C.L. contours in the ∆m A -∆m C plane, focusing on the cos(β − α) dependence given by different color codes, for Higgs (solid curves) and Z-pole precision (dashed curves) constraints individually (left panels), and combined (right panels), with upper rows for m H = 800 GeV, √ λv 2 = 0, middle rows for m H = 800 GeV, √ λv 2 = 300 GeV, and bottom rows for m H = 2000 GeV, √ λv 2 = 0. tan β = 1 is assumed for the plots. For the Higgs precision fit, the alignment limit cos(β − α) = 0 (blue curve) typically gives the largest allowed ranges. Even for small deviation away from the alignment limit, cos(β − α) = ±0.007, ∆m A is constrained to be positive for cos(β − α) = 0.007, and it splits into two branches for cos(β − α) = −0.007. The Z-pole precision measurements force the mass splittings to either ∆m C ∼ 0 or ∆m C ∼ ∆m A , equivalent to m H ± ∼ m H,A . The dependence on cos(β − α) for Z-pole constraints is almost non-noticeable given the small range of cos(β − α) allowed under the current LHC Higgs precision measurements.
Combining both the Higgs and Z-pole precisions (right panels), the range of ∆m C,A are further constrained to be less than about 200 GeV in the alignment limit for m H = 800 GeV, √ λv 2 = 0, with positive (negative) values for the mass splittings preferred for positive (negative) cos(β − α). For √ λv 2 = 300 GeV, loop corrections play a more important role. For cos(β − α) = 0.007, only thin strip of ∆m C ∼ 0 and 0 ∆m A 500 GeV is allowed. For cos(β − α) = −0.007, −250 GeV ∆m C ∼ ∆m A −100 GeV as well as thin slice of ∆m C ∼ 0 for negative ∆m A could be accommodated. For larger m H = 2000 GeV, while the ranges for mass splittings are typically larger under the alignment limit, deviation from the alignment limit leads to tighter constraints due to the suppressed loop contributions.
The Higgs and Z-pole precision measurements at future lepton colliders provide complementary information. While the Z-pole precision is more sensitive to the mass splittings between the charged Higgs boson and the neutral ones (either m H or m A ), the Higgs precision measurements in addition could impose an upper bound on the mass splitting between the neutral ones. Furthermore, the Higgs precision measurements are more sensitive to the parameters cos(β − α), tan β, √ λv 2 and the masses of heavy Higgs bosons.

Comparison between different lepton colliders
In this section, we present a brief comparison for the potential reach of different machines, including CEPC, FCC-ee, and ILC precision shown in Tab. 2 for Z-pole precision and Tab. 3 for Higgs precision. In Fig. 12 for the Z-pole precision, FCC-ee has the best performance because of the higher proposed luminosity at Z-pole. For the combined fit, FCC-ee shows the best constraint, dominanted by the Z-pole effects.

Summary and Conclusions
In this paper, we examined the impacts of the precision measurements of the SM parameters at the proposed Z-factories and Higgs factories on the extended Higgs sector. We first summarized the anticipated accuracies on determining the EW observables at the Z-pole and the Higgs factories in Section 2. Those expectations serve as the general guidances and inputs for the following studies for BSM Higgs sector. We illustrated this by studying in great detail the well-motivated theory, the Type-II 2HDM. Previous works focused on either just the tree-level deviations, or loop corrections under the alignment limit, and with the assumption of degenerate masses of the heavy Higgs bosons. In our analyses, we extended the existing results by including the tree-level and one-loop level effects of non-degenerate Higgs masses. The general formulation, theoretical considerations and the existing constraints to the model parameters were presented in Section 3, see Fig. 1−Fig. 3. The main results of the paper were presented in Section 4, where we performed a global fit to the expected precision measurements in the full model-parameter space. We first set up the global χ 2 -fitting framework. We then illustrated the simple case with degenerate heavy Higgs masses as in Fig. 4 with the expected CEPC precision. We found that in the parameter space of cos(β − α) and tan β, the largest 95% C.L. range of | cos(β − α)| 0.008 could be achieved for tan β around 1, with smaller and larger values of tan β tightly constrained by κ g,c and κ b,τ , respectively. Comparing to the tree-level only results [28], cos(β − α) shifts to negative values for tan β > 1. Smaller heavy Higgs masses and larger λv 2 lead to larger loop corrections, as shown in Fig. 5.
The limits on the heavy Higgs masses also depend on tan β, λv 2 and cos(β − α), as shown in Fig. 6 and alternatively in Fig. 7 varying m 2 12 . While the most relaxed limits can be obtained under the alignment limit with small λv 2 , deviation away from the alignment limit leads to much tighter constraints, especially for allowed range of tan β. The reach seen in the m Φ -tan β plane is complementary to direct non-SM Higgs search limits at the LHC and future pp colliders, especially in the intermediate tan β region when the direct search limits are relaxed.
It is important to explore the extent to which the parametric deviations from the degenerate mass case can be probed by the precision measurements. Fig. 8 showed the allowed deviation for ∆m Φ with the expected CEPC precision and Fig. 9 demonstrated the constraints from the individual decay channels of the SM Higgs boson. As shown in Fig. 10, the Higgs precision measurements alone constrain ∆m A,C to be less than about a few hundred GeV, with tighter constraints achieved for small m H , large λv 2 and small/large values of tan β. Z-pole measurements, on the other hand, constrain the deviation from m H ± ∼ m A,H . We found that the expected accuracies at the Z-pole and at a Higgs factory are quite complementary in constraining mass splittings. While Z-pole precision is more sensitive to the mass splittings between the charged Higgs and the neutral ones (either m H or m A ), Higgs precision measurements in addition could impose an upper bound on the mass splitting between the neutral ones. Combining both Higgs and Z-pole precision measurements, the mass splittings are constrained even further, as shown in Fig. 11, especially when deviating from the alignment limit. Furthermore, Higgs precision measurements are more sensitive to parameters like cos(β − α), tan β, √ λv 2 and the masses of heavy Higgs bosons. We found that except for cancelations in some correlated parameter regions, the allowed ranges are typically tan β ∼ 0.2 − 5, | cos(β − α)| < 0.008, |∆m Φ | < 200 GeV . (5.1) For the sake of illustration, we mostly presented our results using the CEPC precision on Higgs and Z-pole measurements. The comparison among different proposed Higgs factories of CEPC, FCC-ee and ILC are shown in Fig. 12 and Fig. 13. While ILC with different centerof-mass energies has slightly better reach in Higgs precision fit, FCC-ee has slightly better reach in Z-pole precisions.
The precision measurements of the SM parameters at the proposed Z and Higgs factories would significantly advance our understanding of the electroweak physics and shed lights on possible new physics beyond the SM, and could be complementary to the direct searches at the LHC and future hadron colliders.