Strong first order electroweak phase transition in 2HDM confronting future Z&Higgs factories

The electroweak phase transition can be made first order by extending the Standard Model (SM) Higgs sector with extra scalars. The same new physics can explain the matter-antimatter asymmetry of the universe by supplying an extra source of CP violation and sphaleron processes. In this paper we study the existence of a strong first order electroweak phase transition (SFOEWPT) in the type-I and type-II two Higgs doublet models (2HDM). We focus on how the SFOEWPT requirement constraints the spectrum of non-SM Higgs. Through a parameter space scan, we find that SFOEWPT suggests upper limits on the masses of heavy Higgs $m_{A/H/H^\pm}$, which is less than 1 TeV. High temperature expansion and Higgs vacuum uplifting is used for an analytical understanding of our results. We also study the probe ability on SFOEWPT from Higgs and $Z$-pole precision measurements at the one-loop level at future Higgs&$Z$ factories. And together with theoretical constraints, sizeable loop corrections require $m_A \approx m_{H^\pm}>m_H$ to meet SFOEWPT condition in Type-II 2HDM.


Introduction
The discovery of the Higgs boson in 2012 completes the Standard Model (SM) [1,2], yet there remain observations that cannot be explained by it. One of the most famous puzzles is the baryon asymmetry of the universe (BAU), which sees the visible matter in our universe being dominated by baryons whilst the amount of anti-baryons is negligible. Particle physics models that can successfully explain the BAU need to satisfy the three Sakharov conditions [3]. The SM was once considered as a candidate model [4][5][6], since baryon number conservation can be broken by an electroweak sphaleron process [4,7,8], the CKM matrix provides CP violation, and the electroweak phase transition can induce a departure from equilibrium if the Higgs boson is light enough. However, such an electroweak baryogenesis (EWBG) mechanism in the SM framework turns out to fail, since the CP violation present in the CKM matrix is too small [9] and the measured Higgs mass is too heavy to trigger a strong first order electroweak phase transition (SFOEWPT) [10,11]. Thus for a successful baryogenesis, new physics beyond the SM (BSM) is required to supply a new source of CP violation and a strong out-of-equilibrium process [12,13]. In this work we focus on the latter issue.
It is well known that, compared with other baryogenesis mechanisms, e.g. leptogenesis [105] or the Affleck-Dine mechanism [106], EWBG can be detected at the electroweak scale and part of the parameter space can be covered by current or expected collider experiments. In 2HDMs, in addition to the SM-like Higgs boson h, there are three non-SM Higgs bosons, H/A/H ± . H/A/H ± couple to h and help to build an energy barrier between the symmetric phase and the SU (2) × U (1) broken phase when the temperature of the universe is around the electroweak scale. Then the phase transition, which is tunneling through the energy barrier, can be first order and strong enough. In our study, we will show that, in order for this strong first order phase transition to occur, the masses of H, A and H ± bosons should all be smaller than about 1 TeV, and generally there needs to be a relatively large mass splitting between the heavy Higgs bosons H, A and H ± .
H, A, and H ± bosons with a mass lighter than 1 TeV can be directly produced at the current LHC or future hadron colliders like the HE-LHC [107] or the SPPC [108,109]. Channels like A/H → tt/bb/ττ , H ± → tb, or A → HZ [110][111][112][113] can be used for detection or exclusion. Besides, through mixing and loop effects, the non-SM Higgs bosons also change the predicted value of the oblique parameters S, T and U , and reduce the Higgs couplings κ i = g h ii 2HDM /g h ii SM relative to the SM expectation. Future e + e − colliders like the ILC [114], FCC-ee [115,116] and CEPC [108,109] will copiously produce Z and Higgs bosons, and thus those observables (especially the hZZ coupling) can be measured with unprecedented precision. In this work, we perform a global fit to obtain the parameter space of 2HDMs that simultaneously satisfies a SFOEWPT and the expected measurement precision at future Z and Higgs factories.
The structure of this paper is as follows. In Section 2 we briefly introduce our 2HDM models and calculation methods. In Section 3 we list all relevant measurements that can be used to constrain the parameter space of the type-I and type-II 2HDMs. Section 4 starts with an analytic analysis which helps readers to understand the features of the electroweak phase transition in 2HDMs. Then we study three simplified typical cases, and present the most general scan result. We conclude this work in Section 5.
2 The electroweak phase transition in 2HDMs 2.1 Two Higgs Doublet Models 2HDMs without a Z 2 symmetry generally induce dangerous flavour-violating couplings at tree level. In this work we therefore consider 2HDMs with a soft Z 2 symmetry breaking. The tree-level scalar potential for a 2HDM can be written as: We consider a CP-conserving case, in which all mass parameters m 2 ij and quartic couplings λ i are real. After electroweak symmetry breaking (EWSB), the two SU (2) L Higgs doublets Φ i obtain VEVs v i , and they can be expanded in the component real scalar fields: . (2.2) We further define the ratio of VEVs as tan β ≡ v 2 /v 1 . Two of the three m 2 ij can be replaced by other parameters by imposing conditions that result from minimising the Higgs potential Thus the squared mass matrices of the CP-even, CP-odd, and charged Higgs are: Here λ 345 ≡ λ 3 + λ 4 + λ 5 . After diagonalization, the mass eigenstates are related to the original fields by the rotation matrices: We choose our input parameters to be: The mass of the SM-like Higgs boson m h is fixed to the current central measured value 125.09 GeV [117]. Then the λ i can be re-expressed in terms of these input parameters. Considering the theoretical constraints, including vacuum stability, perturbativity, and unitarity, we introduce sin β cos β , (2.12) following the notation in [118]. Under the assumption of degenerate heavy Higgs masses m H = m A = m H ± , there is no theoretical restriction on the tan β range when √ λv 2 = 0. Type-I and Type-II 2HDMs have different Z 2 parity assignments, and thus the couplings between scalar and other particles have a different dependence on tan β and the mixing angle α. The main difference between the Type I and Type II models is the dependence of the couplings Aff and Hff on the value of tan β. Couplings between A/H and down-type fermions are suppressed by 1 tan β in the Type I model, but are enhanced by tan β in the Type II model. Thus the Type II model is generally more constrained by experiments than the Type I model when tan β is large.
Here we need to emphasize that in the 2HDM we can set the mass of the non-SM Higgs bosons A/H/H ± to an arbitrarily high scale. This is because of the presence of m 2 11,12,22 , with m 2 12 breaking Z 2 symmetry in Eq. 2.1. As can be seen from Eq. 2.5 to Eq. 2.7 , the squared masses of A/H/H ± arise from two types of contribution. One of them involves terms of the form λ i v j v k , which are bounded by perturbative unitarity and thus cannot be too large. Upper-limits on these terms are roughly given by 4πv 2 ≈ (870 GeV) 2 . Another part of Higgs mass squares come from m 2 12 (m 2 11/22 are transformed through Eq. (2.3) and Eq. (2.4)), and these terms can in principle be set to any value without violating theoretical requirements. This makes the search for evidence of 2HDMs an endless game: you can never completely falsify a New Physics model containing hypothetical particles which have no upper limits on their mass.
However, in the following part of this work we will show that the requirement of a SFOEWPT imposes upper limits on the masses of the A/H/H ± bosons, making it possible to fully verify or falsify the idea of EWBG in 2HDMs in the near future.

Thermal effective potential
To study the phase transition in the early universe, we need to study the dependence of the free energy density on the order parameter. In our case, the free energy density is the thermal effective potential, and the order parameter is the homogeneous scalar VEV [119]. The thermal effective potential V (φ 1 , φ 2 , T ) at temperature T is composed of four parts: Here V 0 is the tree-level potential of our model, V CW is one-loop Coleman-Weinberg potential, V CT is the counter term, and V T is the thermal correction. The tree-level potential V 0 (φ 1 , φ 2 ) is obtained by replacing the field operators Φ 1 (x) and (2.14) The one-loop Coleman-Weinberg potential V CW (φ 1 , φ 2 ) is given in the MS renormalization scheme by [120]: with the index i running over all massive particles. n i is the degrees of freedom of particle i multiplied by (−1) 2s (s is the spin of particle i), which is -12, -4, 6, 3, 2, 1, 2 and 1 for quarks, leptons, W ± , Z, H ± , G 0 , G ± , and neutral scalars, respectively. c i is 5 6 for gauge bosons, and 3 2 for other particles. m 2 i (φ 1 , φ 2 ) is the mass square of particle i with v 1 and v 2 in its expression being replaced by scalar field value φ 1 and φ 2 . The renormalization scale µ is set to the zero temperature VEV v.
In Eq. 2.11 we choose the scalar masses, mixing angle, and VEV ratio as our input parameters. These parameters are considered as physical parameters. It means that the VEVs are determined by the position of the minimum of the scalar potential, and squared masses are given by the second order partial derivatives of the scalar potential with respect to the scalar fields at the position of the minimum. Adding the Coleman-Weinberg correction will shift both the position of the minimum and the second order partial derivatives of the tree-level potential.
Thus, in order to offset the modification, counter terms V CT (Φ 1 , Φ 2 ) need to be added to the Lagrangian. For a CP-conserving 2HDM, V CT (Φ 1 , Φ 2 ) can be expressed as [52]: Coefficients of counter terms, those δs, need to be fixed by "on-shell" conditions: with ψ i denoting all of the component scalar fields of Φ 1 and Φ 2 . These conditions are evaluated at the minimum of the scalar potential at zero temperature, where After adding these counter terms, our input parameters can be treated as physical parameters which are directly connected to observables.
The thermal correction with ring resummation included is [121,122]: Here, the index i denotes all gauge bosons and scalars, j denotes leptons and quarks, and k denotes scalars and the longitudinal component of gauge bosons. The functions J B,F are two integrals which come from the scalar and fermion thermal corrections respectively: The second line in 2.19 comes from ring resummation, which is used to avoid the infrared divergence that occurs when the scalar mass is much smaller than the temperature.m 2 k (φ 1 , φ 2 , T ) is the thermal Debye mass, an expression for which can be found in the literature [52,122].

Numerical analysis method
An electroweak phase transition is considered to be strong enough only if the net baryon number generated around the bubble wall is not significantly washed out by the sphaleron process inside the bubble. This condition can be converted to the requirement on the value of "wash out" parameter [123]: Here T c is the critical temperature where a second minimum of V (φ 1 , φ 2 , T ) that breaks SU (2) × U (1) appears, and v c ≡ v 2 1 (T c ) + v 2 2 (T c ) reflects the scale of electroweak symmetry breaking. Here v 1 (T c ) and v 2 (T c ) are the scalar field values which minimize V (φ 1 , φ 2 , T c ). 1 Second order derivatives of V CW suffer from an infrared divergence originating from the massless Goldstone boson when T = 0. This problem can be solved by introducing an IR cut-off mass [50].
The calculation of ξ c suffers from theoretical uncertainties. The first problem is that the ξ c induced by V (φ 1 , φ 2 , T ) is not gauge independent by itself [124][125][126]. Missing higher-order quantum corrections also induce a theoretical uncertainty [127]. For a concrete model, one can use lattice simulations to obtain a reliable value of ξ c [63], but such a non-perturbative calculation is very computationally expensive. Being aware of the theoretical uncertainty in the calculation of ξ c , in this work we relax the criterion of a SFOEWPT to ξ c ≡ vc Tc > 0.9. On the other hand, for a first order phase transition to really happen in the universe, the bubble nucleation rate should be larger than the Hubble expansion rate at the nucleation temperature [33,128]. This requirement can be considered as a further constraint on the 2HDM parameter space. For a conservative estimate, in this work we will not consider a requirement on the bubble nucleation rate.
Analytically, T c and v c can be obtained by solving the following equations: To make (0, 0) and ) also need to be positive definite. However, due to the complicated form of V (φ 1 , φ 2 , T ), solving these equations analytically is quite difficult. Instead, one can search for the critical temperature using a numerical method. There are already public packages which can be used for numerical thermal phase transition analysis, such as CosmoTransitions [129], BSMPT [61], and PhaseTracer [130]. We choose BSMPT for our numerical analysis, since the 2HDM has been implemented in BSMPT as a benchmark model, and BSMPT is written in C++ which helps to save numerical calculation time. In BSMPT, the search for T c is started from a high temperature (the default value is 300 GeV), where the minimum position of V (φ 1 , φ 2 , T ) is (0, 0). Then BSMPT traces the minimum position of V (φ 1 , φ 2 , T ) with decreasing temperature. If BSMPT detects a minimum position jumping (0, 0) ⇒ (v 1 (T ), v 2 (T )) at a certain temperature T , the search stops and the output T is the desired critical temperature T c .
The full thermal phase transition history of the 2HDM could be complicated [59]. Multiple phase transition processes are possible. For baryogenesis, however, only the phase transition that transfers (0, 0) ⇒ (v 1 (T ), v 2 (T )) is relevant. This is because a successful baryogenesis requires the sphaleron rate to be very fast outside the bubble wall, i.e. Γ Sph ∼ (α W T ) 4 . While in the electroweak symmetry breaking phase, the sphaleron rate will be strongly suppressed Thus another phase transition (v 1 (T ), v 2 (T )) ⇒ (w 1 (T ), w 2 (T )) has nothing to do with baryogenesis, because the sphaleron rate outside the bubble will be too low to generate baryon number. We will therefore not take this kind of phase transition into account in this work.

Current and expected bounds
2HDMs are constrained by various theoretical considerations and experimental measurements, such as vacuum stability, perturbativity and unitarity, as well as heavy flavor observations [131], electroweak precision measurements, and LHC Higgs measurements and non-SM Higgs searches [132]. We briefly summarize below the constraints we adopt in the following sections.

Theoretical constraints • Vacuum stability
In order to make the vacuum stable, the scalar potential should be bounded from below [133][134][135][136]: • Perturbativity and unitarity We adopt a general perturbativity condition of |λ i | ≤ 4π, and for the unitarity bound [137][138][139][140][141]: To provide some general insights into the impact of these theoretical constraints, we show in Fig. 1 the allowed regions in the m H − ∆m (left), m H − tan β (middle), and λv 2 − tan β (right) planes, for various fixed values of the other parameters. In the left panel, we take tan β = 3, cos(β − α) = 0, fixing m A = m H ± . Here √ λv 2 = 0, 150, 300, 220, 230 GeV are represented by the red, blue, green, purple, and orange lines, and the region under the lines is allowed by the theoretical constraints. Generally, a larger heavy Higgs mass m H corresponds to a smaller allowed mass splitting ∆m for any specific √ λv 2 . The allowed ∆m also gets smaller when √ λv 2 gets larger, and here there is no region left for √ λv 2 > 232 GeV. In the middle panel with √ λv 2 = 0 GeV, we explore the effect of the parameter cos(β−α). Here, based on the allowed | cos(β −α)| at the current LHC Run-II [142], we take cos(β −α) = ±0.005 (dashed lines), and cos(β −α) = ±0.02 (solid lines) and show the allowed region, which is to the left of the corresponding lines. We fix the mass splitting ∆m = m A/H ± − m H = 200 GeV. Under cos(β − α) = 0, m H < 820 GeV is allowed, independently of tan β. If cos(β − α) = 0, such as the 0.005 region shown by the dashed lines, the allowed regions are Theoretical Constraints: tanβ=3 reduced. As discussed in [118], the allowed regions for opposite-sign cos(β − α) are symmetric around the line tan β = 1.
In the right panel, m H is fixed at 700 GeV, and ∆m = m A/H ± − m H = 0, 100, 200 and 230 are shown. The allowed region is inside of the boundary line. Larger ∆m leads to a smaller allowed λv 2 range, and ∆m > 230 GeV is no longer allowed. For √ λv 2 = 0, there is no restriction on tan β.  164,165], H → hh [166][167][168][169], and A/H → HZ/AZ [170,171]. To investigate the impact on the 2HDM parameter space of the published null results in these searches, we take the cross section times branching fraction limits, σ × BR, from the LHC studies and reinterpret them for our 2HDM model points using the SusHi package [172] to calculate the production cross-section at NNLO level, and the 2HDMC [173] code for Higgs decay branching fractions at tree level.

Direct searches at LHC Run-II
As a first example, taking the benchmark point cos(β − α) = 0, m 2 12 = m 2 H cos β sin β and m A = m H + = m H + 200 GeV, we show the current collider limits in the m H − tan β plane in Fig. 2, for both the Type I and Type II models. The various channels include H/A → bb H/A → tt (dot-dashed magenta) and 4t production (dashed purple), as well as the exotic decay channel A → HZ (blue). Other decays such as A → Zh and H → hh will only contribute if cos(β − α) deviates from zero at tree level [174]. For the Type-I model (left panel of Fig. 2), the exotic decay A → HZ channel covers m H < 2m t , tan β < 5 totally, and can reach to tan β = 10. Top quarks searches, 4t + A/H → tt, cover m H < 800 GeV for tan β < 0.3, and m H < 650 GeV for tan β < 1.1 A/H → τ τ, γγ then exclude the region m < 350 GeV, tan β < 1 Generally because of the cot β-enhanced Yukawa coupling in Type-I model, only the small tan β region can be explored [132]. In the Type-II 2HDM (right panel), the top quark and H/A → γγ constraints are similar to those for the Type-I model, while the fermionic decays A/H → bb, τ τ, µµ could exclude m H to 800 GeV when tan β > 10 generally. Since the Hbb, and Hτ τ couplings are tan β-enhanced, the A → HZ decay contributes a lot at medium and large tan β regions, tan β > 0.5, m H < 2m t and tan β > 15, m H < 600 GeV Thus m H < 2m t is totally excluded by all channels together in Type-II model, and only 1.5 < tan β < 10 is allowed for m H < 650 GeV, which is important for our later study of the electroweak phase transition.

Higgs and Z pole precision measurements
The SM has been tested with high precision via observables measured at the Z-pole from LEP-I [175] and the LHC [176]. Future lepton colliders will further improve the precision of measurements in the Higgs sector, and we therefore include hypothetical future lepton collider results in our study. In Ref. [177], it was shown that the precision reached by several future e + e − machines, including the CEPC program with an integrated luminosity of 5.6 ab −1 [108,109], the FCC-ee program with 5 ab −1 of integrated luminosity [115,116], and the ILC with various center-of-mass energies [114], is similar. Thus, following the approach adopted in Ref. [177,178], we will explore the CEPC proposals in detail.
In our analyses, we take the S, T, U data at 95% Confidence Level (C.L.) from Table 2 of Ref. [177], and the Higgs precision measurements from Table 3 in the same reference. We use a χ 2 profile-likelihood fit, being the current best-fit central value for current measurements, and 0 for future measurements at the first term for Z sector. The σ ij are the error matrix, σ 2 ij ≡ σ i ρ ij σ j with σ i and correlation matrix ρ ij given in [177]. For the second term about Higgs sector, Higgs precision measurements are used to perform global fit with µ 2HDM i = (σ i ×Br i ) 2HDM /(σ i ×Br i ) SM is the signal strength for various Higgs search channels, σ µ i is the estimated error for each process. The studies [177,178] show that one-loop level electroweak corrections to SM Higgs couplings have probe ability to heavy Higgs with Higgs precision measurements, and thus our study of For future colliders, the various µ obs i are set to be unity in the current analyses, assuming no deviations from the SM observables.
In the following analyses, the overall χ 2 is calculated, and use to determine the allowed parameter region at the 95% C.L. For the one-, and two-parameter fits, the corresponding ∆χ 2 = χ 2 − χ 2 min values at the 95% C.L. are 3.84, and 5.99 respectively.

Flavour constraints
The charged Higgs H ± boson couples to up and down type fermions, and thus observations from flavor physics put strong bounds on its mass and couplings [179]. Among various flavor observations, measurements related to B physics provide the most stringent limits on tan β and m H ± . For example, m H ± < 580 GeV in the Type-II 2HDM has been excluded by the measurement of BR(B → X s γ) [180]. ∆M Bs and BR(B s → µ + µ − ) further exclude m H ± < 1 TeV in the Type-II 2HDM when tan β < 0.7. The region with tan β < 1 and m H ± < 1 TeV in the Type-I 2HDM has been excluded by B → X s γ [180]. In our study, we take these constraints into account.

The Phase Transition of 2HDM
To get a better understanding of the electroweak phase transition in 2HDMs, we will first discuss it in the context of some approximate or limiting cases, focusing on the relationship between the phase transition and the Higgs vacuum uplifting. Then we will consider several benchmark cases, varying one or two parameters to dig into the effects of constraints as well as features of the Higgs potential. Our general results will follow these specific cases.

High Temperature Expansion
Due to the complicated form of the thermal effective potential Eq. (2.13) and its intricate thermal evolution history, it is difficult to tell whether a specific point can successfully trigger a SFOEWPT in the early universe through a simple formula or argument. To simplify the analysis of the phase transition, people generally use a high temperature expansion, limited to the leading terms of the thermal correction functions J B and J F . Then the thermal effective potential can be simplified to a polynomial function of the Higgs field value: Here φ h ≡ cos βφ 1 + sin βφ 2 is the scalar field that breaks the SU (2) × U (1) symmetry at zero temperature. Due to the simple form of Eq. 4.1, we can use the minimization condition and V (0, T c ) = V (v c , T c ) to directly calculate the wash-out parameter: At tree-level, the coefficients µ 2 andλ in Eq. (4.1) are: The coefficients D and E are induced from the leading thermal corrections: T (m 2 (φ h )) 3/2 + · · · (4.4) Here m 2 (φ h ) is the mass square of a massive particle with v 2 in it being replaced by φ 2 Considering the most massive particles in the 2HDM, D and E can be expressed as: In the expression for E, the term E (H/A/H ± ) denotes the contributions from the non-SM Higgs bosons. We cannot explicitly write out the expression for E (H/A/H ± ) because, as we said in Section 2.1, the mass of the H/A/H ± bosons come from two sources. Schematically, the φ h -dependent non-SM Higgs squared masses can be expressed as: Here M 2 = m 2 12 sin β cos β is the scale at which the Z 2 symmetry is broken. α can be A, H, or H ± , and λ α is a linear combination of the λ i (i = 1, 2, 3, 4, 5) parameters. In the alignment limit cos(β − α) = 0, the expressions for λ α are: ) So the non-SM Higgs bosons provide a term in V (φ h , T ) which is not exactly proportional to φ 3 h : We can simplify the above expression in two limiting cases: And so in these two limiting cases: The above expression needs to be multiplied by 2 if α is H ± . Although expression 4.13 is obtained in a limiting case, it helps us to understand which of the input parameters are particularly relevant for a SPOEWPT. When the non-SM Higgs masses are dominated by M 2 , the spectrum tends to be degenerate, and the phase transition strength tends to be reduced as the non-SM Higgs boson masses increase. When the non-SM Higgs masses are dominated by λ α v 2 , the spectrum tends to be split, and the phase transition strength tends to be increased as the non-SM Higgs boson masses increase.

Higgs Vacuum Uplifting
Another method that can help us to understand which parameters are important for SFOEWPT, is to calculate the depth of the zero temperature Higgs potential [55]. For a shallow Higgs potential, it is easier to develop an energy barrier between the symmetric phase and the broken phase than for a deep Higgs potential, when the temperature is high. Thus generally speaking, there is an inverse relation between the phase transition strength and the depth of the vacuum energy. We follow the notation at Ref. [55] and define the SM vacuum energy density as F SM 0 . The value of F SM 0 is about −1.25 × 10 8 GeV 4 . The vacuum energy density of the 2HDM is denoted by F 0 . We can further define a dimensionless parameter: ∆F 0 /|F SM 0 | > 0 means that the 2HDM vacuum energy is uplifted from the SM value, whilst ∆F 0 /|F SM 0 | cannot exceed 1, otherwise the zero temperature vacuum will be unstable. The numerical results in [55] show a positive correlation between ξ c and the parameter ∆F 0 /|F SM 0 |. However, we find that the relationship is only valid for m H ≤ 500 GeV, the range Ref [55] explored, and the parameters may become negatively-correlated for large m H . To illustrate this, here we refine their analysis by considering a benchmark case: sin β cos β = 0 (to meet the theoretical constraints, as in the right panel of Fig. 1).  The one-loop level Higgs vacuum uplifting in the alignment limit cos(β − α) = 0 has been given in Ref [55]: To illustrate the idea underlying Higgs vacuum uplifting, here we display the whole shape of the zero temperature Higgs potential. In the left panel of Fig. 3, we present the zero temperature Higgs potential along the φ h ≡ cos βφ 1 + sin βφ 2 direction, with m H = 200, 400, 600, 800 GeV represented by red, orange, green and blue lines respectively . The SM Higgs potential is also shown with black dashed line for comparison. It is clear that as m H increases, the height of the minimum point of the Higgs potential continues to rise, and the shape of the Higgs potential becomes shallower. For large m H , F 0 > 0, generating an unstable vacuum. Thus for a stable vacuum m H cannot be too large.
To find the relationship between ξ c and ∆F 0 /|F SM 0 |, in the right panel of Fig. 3 we present both ∆F 0 /|F SM 0 | and ξ c as functions of m H . In the plot, the left y axis is for ∆F 0 /|F SM 0 | with the red dashed line representing the relationship with m H . While ξ c , the right y axis, is shown by the solid green line. Here ξ c is calculated numerically from the package BSMPT.
As the red dashed line, it is clear that there is a linear relationship between ∆F 0 /|F SM 0 | and m H , similar as the left panel. But as the green line shows, ξ c is not monotonically dependent on m H and gets the maximum value around m H = 500 GeV. This result can be understood by our high temperature expansion analysis. Generally as m H , equal to M 2 in our scenario to meet theoretical constraints, becomes too large, the non-SM Higgs mass is dominated by M 2 and E (H/A/H ± ) get smaller as Eq. (4.13). Thus the phase transition strength becomes weaker as m H increases from Eq. (4.7) and Eq. (4.1) .
Since ∆F 0 /|F SM 0 | always gets larger when m H grows, while ξ c gets larger at first (m H < 500 GeV here), and then gets smaller, we can conclude ∆F 0 /|F SM 0 | is not monotonically correlated with ξ c . This conclusion is different to the previous study [46].
In order to get a more robust relationship between ∆F 0 /|F SM 0 | and ξ c , as well as exploring the mass splitting effects ∆m = m A/H ± − m H , we extend the benchmark case by including different mass splittings between the non-SM Higgs bosons: The reason for us to consider different mass splittings is that the mass splitting between different non-SM Higgs bosons is roughly proportional to the size of the couplings λ i . Generally speaking, the greater the couplings λ i , the easier it is for the non-SM Higgs bosons to change the shape of the Higgs thermal potential from Eq. To understand our scan results, we need to invoke the analysis we performed in the last subsection. In our scenario, we have the following relationships between different parameters: Thus, following the discussion we presented in the last subsection, if the value of ∆m is fixed and m H is not too large, the phase transition strength will increase as m H and ∆F 0 /|F SM 0 | increase. But if m H becomes too large and dominates m A/H ± , the phase transition strength will decrease as m H increases, until the vacuum becomes unstable, i.e. ∆F 0 = |F SM 0 |. In the right panel of Fig. 4, we therefore observe that ξ c first rises as ∆F 0 = |F SM 0 | increases (equivalent to m H increasing), and then ξ c decreases as ∆F 0 = |F SM 0 | (and m H ) continues to increase. For the right panel of Fig. 4, depending on the mass splitting and the phase transition features, we can divide the parameter space into three regions: 1. The small mass splitting region, with mass splitting ∆m < 160GeV. In this case, the Higgs vacuum energy cannot be uplifted too high, which means that these points are safe from vacuum stability bounds, and m H can vary from 200GeV to 1TeV. Due to the small value of ∆m, however, λ A/H ± v 2 is too small to satisfy the requirement of a SFOEWPT.
2. The medium mass splitting region, with mass splitting ∆m ∈ [160GeV, 230GeV]. In this case, most of the parameter space is still safe from the vacuum stability constraint.
When m H is not too large, it helps to enhance the phase transition strength. As m H grows to dominate the mass expression of m A/H ± , ξ c rapidly decreases. We also observe that ξ c rises firstly and then falls as ∆F 0 /|F SM 0 | increases. The middle region of ∆F 0 /|F SM 0 | is favored by the existence of a SFOEWPT.
3. The large mass splitting region, with mass splitting ∆m > 230GeV. In this region ∆F 0 /|F SM 0 | starts from a value that is larger than 0.4, and quickly touches the vacuum stability bound ∆F 0 /|F SM 0 | = 1 as m H increases. This means that, before m H has increased to be able to dominate m A/H ± , the vacuum is already unstable. We therefore observe that ξ c increases nearly monotonically as ∆F 0 /|F SM 0 | increasing.
Through the above discussion, it is clear that the upper limits on the non-SM Higgs boson masses in the 2HDM come from vacuum stability (when the mass splitting is large), or the SFOEWPT requirement (when the mass splitting is medium-large). Without the SFOEWPT requirement, the non-SM Higgs bosons can be arbitrarily heavy without violating vacuum stability, provided the mass splitting between them is small enough.
On the other hand, the black dotted lines in the right panel of Fig. 4 clearly show the relationship between ξ c and m H . We found that ξ c is a monotonically increasing function of ∆F 0 /|F SM 0 | when m H < 500 GeV , such as m H = 250, 450 GeV in the right panel of Fig. 4. But with larger mass, the phase transition strength ξ c gets smaller, and when m H > 850 GeV, ξ c can no longer reach 0.9. To avoid the unstable vacuum, larger m H needs smaller ∆m as in Fig. 3. Therefore, a too large m H will result in a too small λ A/H ± v 2 , which could not generate a SFOEWPT. In Table 1

Case1: alignment limit with fixed mass splitting
Following the previous approximate analysis of the electroweak phase transition, we now investigate a series of benchmark cases, starting with, Here we take After these theoretical and experimental constraints, the allowed parameter region is approximately located around m H ∈ (380, 830) GeV, and tan β ∈ (1, 10). The colored region m H ∈ (380, 700) GeV shows the parameter space which can generate a SFOEWPT, with dashed lines indicating the phase transition strength ξ c . We can see that, generally, the strength ξ c gets its maximal value around m H = 500 GeV, which is discussed in the right panel of Fig. 3. The green dash-dotted lines show ∆F 0 /|F SM 0 | = 0.42, 0.53, 0.63, which grows with larger m H and is independent of tan β. We can therefore again conclude that the SFOEWPT strength is not monotonically dependent on m H or ∆F 0 /|F SM 0 |. Finally there is a black band region round m H = 700 GeV, which means that the phase transition strength ξ c < 0.9. Beside the black band region, there is a grey region which is allowed by various constraints, but ξ c in this region has no value. This is because the phase transition in this region is not first order, and thus we can not find the critical temperature and calculate ξ c . We have also checked that Higgs and Z-pole precision measurements give no constraints in this case since cos(β − α) = 0 and ∆m = 200 GeV. The allowed regions are divided into three parts, the colorful region with ξ c > 0.9, the light grey region (mostly above the colorful region) with ξ c < 0.9, and the dark grey region in which a SFOEWPT cannot occur.
Again here √ λv 2 = 0 is set to avoid the constraints on the parameter tan β. In Eq. (4.2), we show the constraints arising from the requirement of a SFOEWPT and other observables in the m H − m A/H ± plane of the Type-II 2HDM. For the various heavy Higgs search channels, only A → HZ gives a visible constraint (shown by the blue region), which can exclude the region with m H < 350 GeV, m A/H ± < 800 GeV. B-physics constraints, shown by the hatched cyan dashed line, exclude m H ± < 580 GeV. Since here we have m A = m H ± and cos(β − α) = 0, the Higgs and Z-pole precision constraints are satisfied automatically. On the other hand, the theoretical constraints, indicated by hatched black lines, give a strong limit on the mass splitting range, roughly ∆m = m A/H ± − m H ∈ (−50, 200) GeV.
The allowed regions are divided into three parts, the colorful region with ξ c > 0.9, the light grey region which is mostly above the colorful region with ∆m = m A/H ± − m H ≈ 200 GeV.) with ξ c < 0.9, and the dark grey region in which a phase transition cannot occur. From the colored region, we find that, both a too large or too small ∆m will not allow for a SFOEWPT. As discussed in Fig. 4, for a too small ∆m, the Higgs vacuum energy cannot be uplifted high enough to generate a phase transition, while too large a value of ∆m will result in an unstable potential F 0 = |F SM 0 |, where the the potential at second EW minimal is higher than the it at the origin. This is also responsible for the upper limit on m H , as the analysis around Eq. (4.18) shows, since too small a value of m H ∆m cannot generate a proper barrier for a SFOEWPT.

Case3: alignment limit with m H =700 GeV
In our previous case studies, we always had the simple assumption of m A = m H ± to satisfy the oblique constraints from Z-pole measurements, and also to simplify the parameter space. Here to study the general mass splitting region, we take another benchmark case, m A ∈ (500, 1200) GeV, m H ± ∈ (500, 1200) GeV, m H = 700 GeV, cos(β − α) = 0, √ λv 2 = 0, tan β = 3.
√ λv 2 = 0 is once more set to avoid the constraints on the parameter tan β, and we take m H = 700 GeV as an example. In Eq. (4.3), we show the electroweak phase transition and other constraints in the plane of m A −m H ± in the Type-II 2HDM. The theoretical constraints are now particularly important, as the region with hatched black lines acts as a boundary on the allowed parameter space. The lower limits on both m A and m H ± are approximately 670 GeV, while the upper limits are 970 and 930 GeV respectively. This is because, once there is a large mass splitting, λ 1−5 will be enlarged [178]. For the various new physics search channels, only the oblique constraints make an effect here. As the hatched blue dashed lines show, the allowed regions are around either m A = m H ± or m H ± = m H = 700 GeV.
The allowed regions are divided into three parts, the colorful region with ξ c > 0.9 allowing a SFOEWPT, the light grey region (mostly above the colorful region) with ξ c < 0.9, and the dark grey region without a first order phase transition. To understand the features here, we also have green dash-dotted lines for ∆F 0 /|F SM 0 |, which gets large when m A , m H ± increases. We also note that, to get a proper vacuum energy uplifting, at least one of m A or m H ± should be large. For instance, F 0 /|F SM 0 | = 0.4 requires m H ± = 900 GeV when m A = m H , or m A = 950 GeV when m A = m H ± , or m A ≈ m H ± ≈ 870 GeV. The region with ξ c > 0.9 is located at F 0 /|F SM 0 | ∈ (0.37, 0.63). The large mass limit comes from F 0 /|F SM 0 | → 1, where the vacuum is not stable, while the small mass limit comes from Eq. (4.18), where there is only limited vacuum uplifting and a barrier to generating a SFOEWPT.

General results
During the last section, we presented three benchmark cases to discuss the effects of the heavy Higgs masses on the existence of a SFOEWPT in the alignment limit, as well as the influence of a variety of theoretical and current experimental constraints up to the one-loop level.
In this section, we present a more general study of Type-I and Type-II 2HDMs. At the same time, we will explore the impact of future results from Higgs factories, presented in Section 3.3, taking the CEPC precision measurements as an example.
Our parameter scan regions for both Type-I and Type-II are : |α| < π 2 , tan β ∈ (0.2, 50), m A ∈ (10, 1500) GeV , m H ± ∈ (10, 1500) GeV , We perform a random parameter scan in the above parameter region, with the total number of samples exceeding 1 billion, for both Type-I and Type-II models.
In Fig. 8 we show the scan results for the Type-II 2HDM. The grey scatter points are the regions allowed by B physics, theoretical constraints, heavy Higgs direct searches and SM Higgs precision measurements at the current LHC Run-II, and constraints from EW oblique operators. The green points are a subset of the grey ones, which can generate a SFOEWPT, and the red points are further required to meet the constraints from future Higgs precision measurements at CEPC. Compared to Case 1 (Section 4.2), which assumed the alignment limit and set m H ± = m A , here we could divide the whole allowed region into 4 classes, • Class A: Regions with m H < 350 GeV. Here the region has m H ± ≈ m A > m H , and the mass splitting is about (300,500) GeV to meet the constraint m H ± > 580 GeV. Generally √ λv 2 ≈ 0 to allow for such a large mass splitting and tan β is within the region selected by the theoretical constraints shown in Fig. 1. This region can also be divided into two subgroups based on sign(κ b ). When sign(κ b ) = +, m H < 200 GeV, tan β ∈ (5, 10) can escape the constraints from the H → τ τ channel as in the right panel of Fig. 2. At the same time, the large mass splitting m A − m H > 450 GeV weakens the constraint from the A → HZ channel [132]. Another subgroup is sign(κ b ) = −, the so-called wrong-sign Yukawa coupling region with cos(β − α) ≈ 2/ tan β. Here m H can reach 350 GeV, cos(β − α) ∈ (0.2, 0.4), and LHC direct searches require tan β < 10 [142]. Because of the large mass splitting in this region, ∆F 0 /|F SM 0 | is too large to produce a stable vacuum.
• Class B: Regions with 5 < tan β < 12 for m H ≈ 450 GeV. This region is also a wrongsign Yukawa coupling region with sign(κ b ) = −. Generally √ λv 2 ≈ 0 to meet theoretical constraints, and m H ± ≈ m A = m H + 140 GeV with tan β < 12 to meet constraints from the A → HZ and A/H → τ τ channels (see Fig. 2). m A/H ± − m H > 140 GeV to meet B physics constraints, while a larger mass splitting is not allowed by theoretical constraints even though √ λv 2 ≈ 0. As the right panel of Fig. 4 shows, because m H = 450 GeV and ∆m = 140 GeV, ∆F 0 /|F SM 0 | < 0.2 is too small, thus the vacuum uplifting is too small, and there is no SFOEWPT here.
• Class C: Regions with 5 < tan β < 45 for 600 < m H < 700 GeV. Here m H ± = m H with m A < m h = 125 GeV. Again it is a wrong-sign Yukawa coupling region with √ λv 2 ≈ 0. The lower limit of m H comes from B physics and EW oblique constraints, and the upper limit comes from theoretical constraints m H − m A < 650 GeV. In the region, ∆F 0 /|F SM 0 | < 0, thus there is no chance to generate a SFOEWPT.
• Class D: The main allowed region with m H > 350 GeV. The region is similar to the white allowed region in the right panel of Fig. 2. Compared to Case 1 with m H ± = m A = m H = 200 GeV in the alignment limit, here the allowed grey region by current LHC Run-II has no upper limit on m H anymore from theoretical constraints when all parameters are free. When m H < 900 GeV, 1 < tan β < 10 is required to satisfy the constraints from H/A → τ τ , top searches and B physics. When m H > 900 GeV, tan β can take a larger value as the constraining power of the H/A → τ τ channel gets weaker. In this region, there are a number of points with ξ c > 0.9, as shown by green points.
We can see the green parameter space has an upper limit of about 900 GeV. For points that also satisfy CEPC constraints as the red points, the parameter space has an upper limit of about 800 GeV. The right panel of Fig. 8 shows the scan results in the plane of ∆m A = m A − m H and ∆m C = m A − m H ± , allowing us to analyze the Class D parameter space. Here the general structure is ∆m C ≈ ∆m A or ∆m C = 0 because of Z-pole oblique constraints.
For the green points from Class D that satisfy LHC Run-II constraints whilst producing a SFOEWPT, there are mainly three regions. For Class D1, ∆m C ≈ ∆m A ∈ (100, 350) GeV. The region has m H ∈ (350, 600) GeV, sin β cos β ≈ (650 GeV) 2 , tan β ≈ 1. In the left panel of Fig. 9, we show the ∆F 0 |F SM 0 | contours in the plane of ∆m A − m C . We can see that when m H = 900 GeV and M = 650 GeV, ∆m C ≈ ∆m A ∈ (−200, −50) GeV, which results in ∆F 0 |F SM 0 | ∈ (0.37, 0.6). This is one of the essential conditions for a SFOEWPT. In this region, √ λv 2 ∈ (500, 600) GeV, thus theoretical constraints impose tan β ≈ 1 as shown in Fig. 1.  Fig. 9, we show the ∆F 0 |F SM 0 | contours for m H = 700 GeV with dash-dotted lines. Class D2 and Class D3, which are allowed by current LHC indirect Higgs precision measurements and direct heavy Higgs searches, will be excluded by Higgs precision observables at the CEPC. This is because the large √ λv 2 in the two regions will lead to large one-loop level corrections to the SM-like Higgs couplings [118,177,178], and large mass splittings around tan β = 1 are not allowed by precise measurements of the Higgs couplings.
We show our general scan results for the Type-I 2HDM in Fig. 10. The allowed grey, green and red points here cover a larger area than for the Type-II model, which mainly comes from heavy Higgs direct search constraints on the large tan β region. As the benchmark case shown in Fig. 2 shows, there is no constraint on tan β > 2 when m H > 2m t in the Type-I 2HDM because all Hff couplings are reduced as tan β increases.
At the same time, there is also a larger range for cos(β − α) at tan β > 2 compared to the Type-II 2HDM [142]. Thus terms involving cos(β − α) will also become important, and from Ref. [55] we can get, ∆F 0 | general =∆F 0 | cos(β−α)=0 + 1 128π 2 cos(β − α) sin(β − α)(tan β − 1 tan β ) here M 2 = m 2 12 sin β cos β . Because of this additional term, once there is sizable tan β, cos(β − α), the allowed ∆m A , ∆m C to generate the proper ∆F 0 /|F SM 0 | range will be a little different to that in the Type-II case. In other words, the allowed parameter space in the Type-I model is larger than that in the Type-II model.
Generally speaking, compared to the Type-II 2HDM, the upper limit of m H allowed by a SFOEWPT in the Type-I model can still reach to 900 GeV. In the Type-II model, such points have ∆m A,C < 0, and are excluded by Higgs and Z-pole precision measurements. But larger tan β values allow larger mass splittings between the heavy Higgs bosons [177,178], and thus in the Type-I model m H → 900 GeV still satisfy these precision measurements. Similarly the regions with ∆m C ≈ 0, ∆m A < 0 or ∆m C ≈ 0, ∆m A > 0 which are not allowed in the Type-II model can still generate a SFOEWPT in the Type-I model.

Conclusion
In this work, we have revisited the existence of a strong first order electroweak phase transition (SFOEWPT) in the Type-I and Type-II 2HDMs. Using both numerical and analytical analysis methods, we pointed out that ∆F 0 /|F SM 0 | is not monotonically related to ξ c as shown in Fig. 3 and Fig. 4. This conclusion is different to that of a previous study [46].
We also found, SFOEWPT suggests the non-SM Higgs bosons, H/A/H ± , have upper limits on their mass as our benchmark Case 1 Section 4.2 and general scan results Fig. 8 and Fig. 10. This limits comes from the combined requirements of vacuum stability at zero temperature and λ H/A/H ± v 2 corrections term at high temperature.
By combining current bounds from LHC direct and indirect Higgs searches, current electroweak precision measurements, flavour physics, and anticipated precision measurements at the future CEPC Z and Higgs factory, we have shown that the requirement of a SFOEWPT puts strong constraints on the mass spectrum of H/A/H ± : For the type-I 2HDM: sin β cos β . In Type-I 2HDM, because of allowed large tan β region from Higgs precision measurements, Class D2 and D3 are still allowed.
Both Type-I and Type-II requires a sizable mass splitting between different heavy non-SM Higgs. And the suggested upper limits of m A/H/H ± is 900 GeV at current stage, and 800 GeV after including Higgs and Z-pole precisions at CEPC. Such a constrained spectrum points out a clear direction for direct searches at the LHC and future colliders.