A new insight into the phase transition in the early Universe with two Higgs doublets

We study the electroweak phase transition in the alignment limit of the CP-conserving two-Higgs-doublet model (2HDM) of Type I and Type II. The effective potential is evaluated at one-loop, where the thermal potential includes Daisy corrections and is reliably approximated by means of a sum of Bessel functions. Both 1-stage and 2-stage electroweak phase transitions are shown to be possible, depending on the pattern of the vacuum development as the Universe cools down. For the 1-stage case focused on in this paper, we analyze the properties of phase transition and discover that the field value of the electroweak symmetry breaking vacuum at the critical temperature at which the first order phase transition occurs is largely correlated with the vacuum depth of the 1-loop potential at zero temperature. We demonstrate that a strong first order electroweak phase transition (SFOEWPT) in the 2HDM is achievable and establish benchmark scenarios leading to different testable signatures at colliders. In addition, we verify that an enhanced triple Higgs coupling (including loop corrections) is a typical feature of the SFOPT driven by the additional doublet. As a result, SFOEWPT might be able to be probed at the LHC and future lepton colliders through Higgs pair production.


Introduction
After the discovery of the 125 GeV Higgs boson [1,2] and the accumulation of LHC data, no evidence of the new physics has been observed yet. Therefore, it is time to inquire whether the Standard Model (SM) of particle physics is actually complete to describe the physics at the electroweak scale. In the meantime, the origin of the baryon asymmetry of the Universe (BAU) is still one of the important open puzzles in particle physics and cosmology. To explain the BAU, the three Sakharov conditions [3] must be fulfilled. The electroweak baryogenesis (EWBG) [4] is a possible mean to account for the generation of an asymmetry (imbalance) between baryons and antibaryons produced in the very early Universe. The success of EWBG requires two crucial ingredients: CP violation and strong first order phase transition (SFOPT), neither of which however can be addressed in the SM framework. First, the SM fails to produce a sufficiently large baryon number due to a shortage of CP violation in the Cabibbo-Kobayashi-Maskawa (CKM) matrix. The other shortcoming of the SM is the absence of departure from thermal equilibrium which could have been realized by a SFOPT: for the observed value of the SM-like Higgs mass this is not accomplished. The phase transition in the early Universe from the symmetric phase to the electroweak symmetry breaking (EWSB) phase actually belongs to a smooth crossover type [5]. It has been derived using lattice computation that the phase transition in the SM can only be strong first order when the Higgs mass is around 70-80 GeV [6][7][8][9]. Therefore, a successful EWBG invokes new physics at the electroweak scale [10]. Theories that go beyond the SM typically have an extended Higgs sector, which may contain the ingredients for a SFOPT as well as new CP-violating interactions as needed for EWBG, usually also producing new signatures at colliders.
The two-Higgs-doublet model (2HDM) is the simplest renormalizable framework to realize EWBG. 1 In this model the scalar potential is extended by an additional SU (2) L doublet, where a charged Higgs together with two additional neutral scalars are introduced. Through their portal interactions with the SM-like Higgs, the finite temperature potential can develop a potential barrier during the Universe cooling down, leading to strengthen the phase transition at the critical temperature. On the other hand, the CP violation can exist either explicitly in these portal couplings or spontaneously via a relative phase between the vacuum expectation values (VEVs) of two doublets [12]. Interestingly, the CP violation (beyond the SM) can be detected indirectly at high precision electric dipole moments (EDMs) experiments. With the recent improvement of the EDMs measurements, the CP violation phases that are needed for the baryon asymmetry are severely constrained [13]. In order to evade the bound one may expect a cancellation arising from the different contributions of EDMs predictions, see Ref. [14][15][16][17][18].
The electroweak phase transition (EWPT) in the 2HDM context has been extensively studied for both the CP conserving case [19][20][21] and with the source of CP-violation [22][23][24][25][26][27]. While the CP phase at zero temperature is supposed to play an insignificant in the EWPT process [19,20,23], the CP-violating phase at finite temperature is found to be important in a recent study [27] where the analysis was performed after taking into account the LHC run-2 constraints. In general, none of the scalar states of the 2HDM resembles a SM-Higgs boson that was observed at the LHC. However, such a SM-like Higgs boson h can arise in the alignment limit, a particularly interesting limit of this model where only one Higgs doublet acquires the total electroweak vev, namely the couplings of H to gauge boson pairs vanish while h possesses SM-like couplings [28][29][30]. In terms of the model parameters, this limit corresponds to sin(β − α) (always positive in our convention) to be close to 1. Driven by the LHC Higgs data, in this paper we focus on the alignment limit (here sin(β − α) ≥ 0.99) of the CP-conserving 2HDMs of Type I and Type II models. We consider the lightest CP-even state h to be the 125 GeV SM-like Higgs observed at the LHC. To proceed the numerical analysis we take the points passing all existing experimental bounds (by the time of paper publication) generated from extensive scans in Ref. [28] and additionally employ the 1-loop improved theoretical constraints, the updated measurements coming from flavor physics and the recent LHC run-2 bounds searching for heavy resonances. Our aim is to identity the parameter space of the 2HDM that can lead to a SFOPT and investigate the implications of a SFOPT required by baryogenesis on the LHC Higgs phenomenology.
It is inspiring to note that the cosmological EWPT can leave signatures of gravitational waves (GW) after the nucleation of the true vacuum bubbles, with typical red-shifted spectrum frequency around O(10 −4 − 10 −2 ) Hz [31], which are detectable in the Evolved Laser Interferometer Space Antenna (eLISA) [32], DECi-hertz Interferometer Gravitational wave Observatory (DECIGO), UltimateDECIGO and Big Bang Observer (BBO) [33]. However, these two effects might be quite incompatible due to an opposite preference occurring in the bubble wall velocity. The baryon asymmetry generation process within EWBG demands a relatively low bubble wall velocity in order to have enough time for the chiral asymmetry generation process to take place, this will later be transformed to the baryon asymmetry by the sphelaron process [10]. Of course, when performing the computation of the BAU in the EWBG mechanism, one should keep in mind that in addition to being subject to large theoretical uncertainties, the detailed calculations of the baryon asymmetry rely on the wall velocity of the bubble generated during the EWPT, see Ref. [34][35][36][37][38][39]. On the contrary, a testable GW signal requires a higher strength of the FOPT and a larger wall velocity. Very intriguingly, the recent development [25] shows that it is possible, although difficult in the 2HDM, to simultaneously accomplish the EWBG and produce the detectable GW signals generated during the EWPT especially through acoustic waves [40,41]. This paper is organized as follows. In Sec. 2 we first briefly review the CP-conserving 2HDMs of Type I and Type II and discuss the status in view of the existing experimental bounds. Next, we describe in Sec. 3 the details of the finite temperature potential and provide a fast numerical handle for the thermal potential. In Sec. 4, the one-stage and two-stage PT are demonstrated and classified. Subsequently, we present in Sec. 5 a useful computational scheme used to single out the one-stage PT and, more importantly, to evaluate the critical temperature T c for the one-stage PT. Having studied the theoretical issues of the model and built the computational tools, we then proceed with the numerical analysis and investigate the properties of the phase transition which are presented in Sec. 6. In particular, the relations between T c and extra Higgs masses as well as the influence of the effective potential at zero temperature on the field value of the electroweak symmetry breaking vacuum are analyzed. In Sec. 7 benchmark scenarios leading to the SFOEWPT are established and their implications for future measurements at colliders are also discussed. Finally, Sec. 8 contains our conclusions and outlook for future studies. In Appendix A, explicit formulas for the thermal mass corrections of the SM gauge bosons are given.

The two-Higgs-doublet model
Let us start with a brief review of the tree-level 2HDM at zero temperature. The general 2HDM is obtained by doubling the scalar sector of the SM, two doublets with identical quantum numbers are present. In general, CP violation may be present in the scalar sector and the Yukawa sector contains generic tree-level flavor changing neutral currents (FCNCs) mediated by the neutral scalar states. Here we consider a CP conserving Higgs sector and the absence of tree-level FCNCs. The first condition is obtained by imposing a reality condition on the parameters of the potential, and the second requirement is achieved by imposing a Type I or Type II structure on the Yukawa sector, this is achieved by imposing a softly-broken Z 2 symmetry [42,43].
Denoting by Φ 1 , Φ 2 the two Higgs doublets, the tree-level potential of this model is expressed as, In this basis, the Z 2 symmetry under which Φ 2 → −Φ 2 is manifest in the quartic terms, while it is softly broken by the introduction of the m 2 12 term. In general, m 2 12 and λ 5 are complex. We consider in this work a CP conserving Higgs sector and set all λ i and m 2 12 as real parameters, see [27] for a CP violating study. In this basis, both Higgs doublets have a non-zero vacuum expectation value (vev). We parametrize the degrees of freedom contained in the Higgs doublets as, where v i are the vevs of the two Higgs doublets. At zero temperature one has the relation v 2 (246 GeV) 2 . For convenience, we use the shorthand notation v ≡ v 0T from now on and define v 1 = v cos β and v 2 = v sin β, tan β is therefore the ratio of the two vevs at T = 0.
The mass parameters m 2 11 and m 2 22 in the potential Eq. (2.1) are determined by the potential minimization conditions, here the shorthand notations s β ≡ sin β, c β ≡ cos β and t β ≡ tan β and λ 345 ≡ λ 3 + λ 4 + λ 5 are employed. Though tan β is a physical parameter here, it is still possible to redefine the two doublets and go to a basis where the full vev resides entirely in one of the two doublets: the so-called Higgs basis (H 1 , H 2 ) [44]. This change of basis is generally possible as observed by the invariance of the gauge kinetic terms of the two doublets under a U(2) Higgs flavor transformation. In the Higgs basis H 1 has the full vev v and thus is precisely the SM Higgs doublet. In general however both neutral components mix upon EWSB and a SM-like CPeven mass eigenstate is not automatic. On the contrary, if one of the CP-even eigenstates is parallel to the neutral H 1 direction, this realizes the alignment limit of the 2HDM [45] which the LHC Higgs data appears to favor. In general, in a basis-independent manner, the alignment limit is defined as the presence of a CP-even eigenstate in the vev direction in the scalar field space.
In the electroweak vacuum, the squared mass matrices in the neutral CP-even, CP-odd and charged scalar sectors are respectively given by, (2.6) The CP-even mass eigenstates h and H, with m h ≤ m H , are obtained through the diagonalization of M 2 P , they are expressed in terms of the neutral components of the doublets as, where the mixing angle α is introduced and is expressed in terms of the entries of the mass matrix. Diagonalization of M 2 A leads to a massive CP-odd scalar A and a massless Goldstone boson G 0 , while M 2 ± leads to a charged state H± and a charged massless Goldstone boson G ± . Their tree-level masses read 2 (2.10) 2 We point out that in the review article [46] a factor of 2 is missing in front of λ5 and λ4 + λ5 terms in the formula of m 2 A and m 2 + , respectively. Type I, II Type I Type II Higgs Using Eqs. (2.8)-(2.10) one can inversely solve the potential parameters, λ 1 , . . . , λ 5 in terms of four physical Higgs masses and the CP-even Higgs mixing angle α, supplemented by the Z 2 soft-breaking parameter m 2 12 [45]. This means that the scalar potential can be entirely determined by these seven parameters and therefore allows us to choose them as a set of complete free inputs for the numerical analysis.
As mentioned previously, we imposed a Z 2 symmetry on the potential Eq. (2.1) in order to forbid Higgs-mediated tree-level FCNCs. Out of the four independent realizations of this symmetry in the fermion sector, we study two of them: the so-called Type I model where only Φ 1 couples to fermions and the Type II model where Φ 1 couples to downtype fermions and Φ 2 to up-type fermions, see [47] for details. These particular structures redefine multiplicatively the Higgs couplings to fermions as compared to the SM predictions, we denote as C U,D,V theses multiplicative factors for the up-type fermions, down-type fermions and massive gauge bosons, respectively. The Higgs couplings to massive gauge bosons do not depend on the Z 2 symmetry charges but are directly obtained from gauge symmetry alone. In Table 1 we present these factors for the three physical scalar states of the theory. Important intuition can be gained by re-expressing these factors in terms of (β −α) and β, in particular to understand their behavior in the alignment limit sin(β −α) ≈ 1: (2.14)

Theoretical constraints
For a viable 2HDM scenario, we require here tree-level stability of the potential, which means that Eq. (2.1) has to be bounded from below, requiring In this work we improve the bounds supplemented by the radiative corrected potential, as will be shown in Sec. 3.2. Additional theoretical constraints from S-matrix unitarity and perturbativity are required. Tree-level unitarity 3 imposes bounds on the size of the quartic couplings λ i or various combinations of them [49,50]. Similarly (often less stringent) bounds on λ i may be obtained from perturbativity arguments.

The experimental constraints
Next, we briefly describe the impact of the experimental bounds on the parameter space of the model. First, electroweak precision data (EWPD), essentially the T parameter, constrains the mass difference between m H ± and m A or m H , one of the two neutral states should indeed be approximatively paired with the charged state in order to restore a custodial symmetry of the Higgs sector [51,52]. Second, the recent measurement on BR(B → X s γ) [53] excludes low values of m H ± 580 GeV in the Type II model [54]. As a consequence, the preferred ranges for the scalar masses are pushed above ∼ 400 GeV. Third, LHC measurements of the 125 GeV signal rates put large constraints on the 2HDM parameter space, in particular they tend to favor the alignment limit where the Higgs couplings are similar to the SM ones. To evaluate these constraints, we use Lilith-1.1.3 [55].
Finally, regarding direct searches, we implement the Run-1 and LEP constraints as performed in [28]. A very important search for the Type II model is in the A, H → τ τ channel, either through gluon-fusion or bb-associated production [56,57]. The ATLAS Run-2 constraint is much stronger than the corresponding Run-1 searches, eliminating larger portion of the parameter space at large tan β in particular. For m A < ∼ 350 GeV we only find few scenarios compatible with the experimental constraints in the Type II model 4 . This is both coming from the aforementioned τ τ search, as well as the H → ZA searches for CPodd state down to 60 GeV. This final state has been searched for by the CMS collaboration during Run-1 [58], and leads to severe constrains of the parameter region. The A → Zh channel has been searched for during both LHC Run-1 [58,59] and Run-2 [60] but the resulting constraints have little impact. In Fig. 1 we show the allowed spectra (red pluses) for the two types of models considered here. The points labeled 'no-EWSB' comes from the requirement of proper EWSB at the 1-loop level, which will be extensively discussed in Sec. 3.2. Due to the severe constraints on the mass spectrum of the extra Higgs bosons, these experimental constraints have significant influence on the requirement of a SFOPT as we will see in Sec. 6.
We now move to investigate the possibility of having a first-order phase transition for the surviving sample points. The interesting question is whether the parameter space that LHC Higgs data favors, simultaneously satisfying both theoretical constraints and experimental bounds, can lead to a favorable prediction for a strong first-order phase transition.  no-EWSB Figure 1. The allowed mass spectra (up to 1.2 TeV) of the BSM Higgs bosons confronting with LHC Run-2. In Type II (right), the B-physics constraint on the charged Higgs m H ± ≥ 580 GeV [54] is imposed and nearly excludes the low m A points. Gray points indicate EWSB is not ensured at zero temperature when one-loop effect is included in the Higgs potential.

The effective potential at finite temperature
To study the phase transition we consider the scalar potential of the model at finite temperature. In the standard analysis, the effective potential which is composed of the tree-level potential at zero temperature V 0 (h 1 , h 2 ) derived in Eq. (3.2), the Coleman-Weinberg one-loop effective potential V CW (h 1 , h 2 ) at T = 0 given in Eq. (3.3), the counter-terms V CT given in Eq. (3.11) being chosen to maintain the tree level relations of the parameters in V 0 , and the leading thermal corrections being denoted by V th (h 1 , h 2 , T ). We discuss these terms separately now.

The tree level potential
Since our model is CP conserving, the classical value of the CP-odd field A is zero and so are the ones for the neutral Goldstone fields. We assume the charged fields do not get VEV during the EWPT process, by taking the classical values for the charged fields to be zero, to strictly respect the U(1) electromagnetic symmetry and therefore ensure the photon massless [46]. 5 The relevant tree level potential V 0 in terms of their classical fields (h 1 , h 2 ) 6 derived from Eq. (2.1) is here we have eliminated m 2 11 and m 2 22 by using the minimization conditions Eq. (2.3).

The Coleman-Weinberg potential at zero temperature
To obtain the radiative corrections of the potential at one-loop level, we use Coleman and Weinberg method [70]. The Coleman-Weinberg (CW) potential in the MS scheme and Landau gauge 7 at 1-loop level has the form: The sum i runs over the contributions from the top fermion, massive W ± , Z bosons, all Higgs bosons and Goldstone bosons 8 ; in the sum s i and n i are the spin and the numbers of degree of freedom for the i-th particle listed in Table 2; Q is a renormalization scale which we fix to Q = v and C i are constants depending on the renormalization scheme. In the MS on-shell scheme employed, C i = 1 2 ( 3 2 ) for the transverse (longitudinal) polarizations of gauge bosons 9 and C i = 3/2 for the particles of other species [72]. Finally, the field- 5 The charge breaking vacuum in multi-Higgs doublet models has been studied in Refs. [61][62][63][64][65][66][67]. Once the U(1) electromagnetic symmetry is broken during the EWPT, the photon acquires mass, which may change the thermal history of the Universe [67]. We leave it to future work. Also, we do not expect the presence of color-breaking vacuum in the process of EWPT since the bosons which actively participate into the evolution of Higgs scalar potential are color neutral. As of our knowledge, color-breaking baryogenesis is achievable in the model with the inclusion of colored bosons (i.e. scalar leptoquarks) [68,69]. 6 To avoid confusion we distinguish the dynamical fields and EW vev in this paper. The classical fields (h1, h2) approach the EW vacuum (v1, v2) at zero temperature. 7 As noted in [71], the VEVs are slightly different in various gauges and the recent study [26] find this effect to be numerically small in the physically interesting regions of parameter space. 8 We ignore the light SM fermions because of the smallness of their masses. In contrast, the inclusion of Goldstone modes is necessary as their masses are non-vanishing for field configurations outside the electroweak vacuum. The photon at zero temperature is strictly massless due to gauge invariance. 9 In most literature Ci = 5/6 is taken for gauge bosons without the distinction between transverse and longitudinal modes. In fact, these two ways of counting are equivalent as the field-dependent mass are identical for both transverse and longitudinal modes at zero temperature. For instance, nZ CZ = 2 × 1/2 + 1 × 3/2 = 3 × 5/6 and nW CW = 2nZ CZ . The mass difference between transverse and longitudinal modes arises from thermal corrections as will see later.
dependent squared massesm 2 i for SM particles include 10 with the corresponding SM Yukawa and gauge couplings being 2m t /v and the ones for scalar bosons are given bŷ where the corresponding matrices M 2 X (X = P, A, ±) are Here the Θ X ij terms listed below are different for X = P, A, ± With V CW being included in the potential, the minimum of the Higgs potential will be slightly shifted, and hence the minimization conditions Eq. (2.3) no longer hold. To maintain these relations, we add the so-called "counter-terms" (CT) [24] 11 , where the relevant coefficients are determined by, 13) which are evaluated at the EW minimum of {h 1 = v 1 , h 2 = v 2 , A = 0} on both sides. As a result, the vevs of h 1 , h 2 and the CP-even mass matrix will not be shifted. One technical difficulty involved at this step arises from the inclusion of the Goldstone bosons in the CW potential. Due to the variation of the scalar field configuration with temperature (which we will see shortly), the Goldstone boson may acquire a non-zero mass at finite temperature, enforcing the inclusion of Goldstone modes in the sum. Nonetheless, in the electroweak vacuum at zero temperature the masses of the Goldstone bosons are vanishing in the Landau gauge, which leads to an infrared (IR) divergence due to the second derivative present in our renormalization conditions Eq. (3.13). This means that renormalizing the Higgs mass at the IR limit is ill-defined [73]. To overcome this divergence, we take a straightforward treatment developed in [24] and impose for Goldstone bosons an IR cut-off at SM Higgs mass, m 2 IR = m 2 h . Although a rigorous prescription used to deal with the Goldstone's IR divergence was developed in [22], Ref. [24] argued that this simple approach can give a good approximation to the exact on-shell renormalization. Practically, in evaluating the derivatives for the CW potential, we remove the Goldstone modes from the sum and add instead the following Goldstone contribution to the right hand of Eq. (3.13) with the replacement for the singular term m 2 G (h 1 , h 2 )| VEV → m 2 IR in the logarithm. Note that the Goldstone bosons have a vanishing contribution to the first derivative evaluated at the vev.
Beyond tree level the true EW vacuum must be preserved when the one-loop corrections are taken into account. This demands that the potential after the inclusion of the CW and counter-terms still form a global minimum at the EW vacuum. As seen in Fig. 2, the CW term (green dotted) often lifts up the potential at the EW vacuum, resulting the local minimum shifting inward or even leading to a false vacuum. On the other hand, the CT effect (blue) drags down the potential at the EW vacuum and thus helps to accomplish a true EW vacuum. As a result, the competition between these two opposite effects determines the existence of a global minimum at the EW vacuum. We present in Fig. 2 two examples where the left one accomplishes a true EW vacuum, while the potential in the right plot has only a local minimum at v. The latter example is phrased 'no-EWSB' in our terminology and such type of points are displayed in Fig. 1. This is an additional important constraint that excludes about 10% (5%) points in the Type I (II) model, in particular for the points with m A ≤ 300 GeV. In the m A < m h /2 regime (termed low-m A scenario), it turns out that EWSB at zero temperature can be achieved as long as at least   one lighter H or H ± is present in the spectrum. For the case where both H and H ± are heavier than ∼ 550 GeV, EWSB would be hardly successful. To understand this, we display in Fig. 3 the one-loop potential at zero temperature (including the CW potential and counter-terms). For simplicity, we assume m H = m H ± , which is typical mass spectra required by the T parameter 12 . The authors of Ref. [74] have shown that low-m A scenario can be phenomenologically alive in the parameter space where the SM-like Higgs h has very small coupling to AA, which leads to tan β 2 or tan β 12 for m H = 600 GeV in the deep alignment limit. The low tan β solution requires a severe tuning in the parameter m 2 12 , and in the allowed range m 2 12 5000 GeV 2 the zero temperature potential (c.f. the 12 The lighter the CP-odd state A is, the stronger the degeneracy between H and H ± should be. red solid line) at EW vacuum v is higher than the one at the origin. Moreover, a proper EW vacuum can be developed as the symmetry soft-breaking parameter m 2 12 increases. This can be achievable for the case of m A ≥ m h /2 where the h → AA decay is kinematically suppressed. On the other hand, the large tan β solution, though possible in Type I model, strongly constrains m 2 12 and tends to lift the potential. Hence, the importance of this class of solution is very marginal and no points were found in our numerical analysis. In addition, tan β 5 in Type II model was already excluded by the CMS bound searching for a light pseudoscalar scalar in the mass range of 20-80 GeV through the bottom-quark associated production and decaying into τ τ final states during Run-1 [75]. As a comparison, we also exhibit the potential at a lower common scale m H = m H ± = 300 GeV. This example is only applicable in Type I model. One can observe that the potential generically reaches a global minimum at the EW vacuum and the depth of this minimum is less sensitive to tan β. This implies that when the new scalars introduced are not heavy, the loop effect is not substantial and thus the potential is largely governed by the tree-level.
We conclude that the requirement of proper EWSB at zero temperature, in synergy with m H ± ≥ 580 GeV required by B-physics measurements [53], entirely exclude the scenario of existing a light pseudoscalar A in Type II model that was delicately studied in Ref. [74]. As will show shortly, these theoretical constraints will play an important role in achieving a strong first-order phase transition.

The thermal effective potential
The finite temperature corrections to the effective potential at one-loop are given by [76] where the functions J B,F are with y ≡ m 2 i (h 1 , h 2 )/T 2 and the upper (lower) sign corresponds to bosonic (fermionic) contributions. The numerical evaluation of this exact integral is very time-consuming (notably for the y < 0 case present for the bosonic degrees of freedom). Thus, computational techniques to reduce the computation time are welcome. A widely used solution is to consider the asymptotic expansions of J B,F . At small y (y 1) 13 , Eq. (3.16) can be approximated by  where a F = π 2 exp( 3 2 − 2γ E ) and a B = 16a F with the Euler constant γ E = 0.5772156649. Whereas at large y, In order to make a quantitive assessment of the approximation precision we plot in Fig. 4 the small/large y approximations as well as the direct numerical evaluation of the integral. (For the evaluation of the latter one we use the NIntegrate function built in Mathematica.) It is clearly seen that the small y approximation (red curve) is valid in the ranges y ∈ (−5, 5) for bosons and y ∈ (0, 5) for fermions, while the large y expansion (blue curve) converges to the exact integral for y > 10 for both functions. A gap is then present between the small and large y approximations in the transition range y ∈ (5, 10). In this situation an interpolation can be introduced to connect smoothly the two approximations. Even though this reduces the deviation of the approximate results from the exact integral to less than 2%, there are still two serious shortcomings. First, this requires a conditional judgement for each state at temperature T to know which approximation should be applied, this largely increases the evaluation time. Second, the above approximations Eqs. (3.17)- (3.19) are only valid for y > 0 as shown in Fig. 4. However, the eigenvalues of the mass matrix of the neutral scalar states can become negative depending on the field configuration. 14 If this happens, Ref. [21] suggests that only the real part of the integral J B should be chosen in the evaluation as the imaginary part is irrelevant in extracting the global minimum. 15 The thermal integrals J B,F given by Eq. (3.16) can be expressed as an infinite sum of modified Bessel functions of the second kind K n (x) with n = 2 [77], with the upper (lower) sign corresponds to bosonic (fermionic) contributions. Our numerical results show that the leading order l = 1 does not provide a good approximation of the full integrals. Instead, inclusion up to l = 5 order in the expansion can match the exact integral very well for both positive and negative y values. Therefore, in this work we take N = 5 in the evaluation of the thermal integrals Eq. (3.20). 16 Fig. 4 also shows that the thermal function is negative for positive y thus dragging the potential down and leading to the formation of two degenerate vacua. As expected, this dragging effect arising from the temperature corrections diminishes as y approaches to the infinity, which corresponds to zero temperature or the decoupling limit. Finally, there is another important part of the thermal corrections to the scalar masses coming from the resummation of ring (or daisy) diagrams [79,80], where M 2 i (h 1 , h 2 , T ) are the thermal Debye masses of the bosons corresponding to the eigenvalues of the full mass matrix which consists of the field dependent mass matrices at T = 0 Eq. (3.9) and the finite temperature correction to the mass function Π X , (X = P, A, ±) given by with the diagonal terms being is the known SM contribution from the SU(2) L and U(1) Y gauge fields and the top quark [79]. It is important to note that the temperature corrections are independent of λ 5 where a possible CP phase can reside. On the other hand, the leading correction to off-diagonal thermal mass is vanishingly small due to Z 2 symmetry imposed in the scalar sector. Moreover, it was argued by [81] that subleading thermal corrections to off-diagonal self-energies are suppressed by additional powers of coupling constants and EW vevs which are usually neglected. Therefore, we shall treat the thermal mass correction Π i as diagonal matrices in the following numerical analysis. The thermal mass corrections of the SM gauge bosons are given in Appendix A.
Historically, there was an alternative algorithm proposed by Parwani in dealing with the thermal corrections [82]. He included the effect of thermal correction from Daisy diagrams by means of substituting ). It is important to note that these two approaches are not physically equivalent and the results produced are quantitively incompatible [21]. They differ in the organization of the perturbative expansion and consistent implementation of higher order terms. The method formulated in Eq. (3.21) restricts the corrections to the thermal masses at one-loop level, whereas Parwani's method inconsistently blends higher-order contributions. Because of this dangerous artifact unrealistically large values of the phase transition strength ξ (defined in Eq. (6.1)) would be obtained. Therefore, we will adopt the former consistent method in the following analysis.

Phase transition: classification
In general, a system may transit from one symmetry phase to another one. Here the electroweak symmetry is broken as the Universe cools down, this is singled as a change in the nature of the global 0-vacuum at high temperature that gets replaced by an electroweak breaking global vacuum at lower temperature. At any given set of parameters, the full effective potential Eq. (3.1) can have several extrema. Our major interest is the global minimum vacuum state, the deepest minimum of the potential. The other extrema can be either saddle points or maxima or local minima of the potential. In studying the thermal phase transition, it is useful to trace the evolution of the extrema as well as calculate the difference in potential depth between the global minimum (called true electroweak (EW) vacuum) and a secondary local minimum.
First, since at very high temperatures electroweak symmetry is not broken, the effective potential has one global minimum, which tends towards the point (h 1 , h 2 ) = (0, 0). We refer to this minimum as the 0-vacuum. As the Universe is cooling down, the parameters that characterize the thermal effects of the model evolve with temperature. This leads to a change of the classical values of h 1 , h 2 fields 17 and thermal phase transition takes place. In general one can classify the thermal phase transition according to the behavior of the vacuum development during the cooling down. For instance, the phase transition may be of first or second order, one-stage or two-stage process. In Fig. 5   origin at T = 163 GeV, and then moves closely along the yellow line until reaching the EW vacuum at zero temperature. This phase transition is called of second order, because the potential minimum shifts continuously while no potential barrier develops during the cooling down. In contrast, the vacuum of the potential displayed in the middle panel is localized in the vicinity of the origin point at high temperature. When the temperature decreases to T 157 GeV an EW vacuum located away from the origin appears, forming two degenerate vacua separated by an energy barrier. This gives rise to the first-order phase transition from the origin to the EW vacuum, which is indicated by the red arrow. In these two examples, the phase transition is termed one-stage. In addition to experiencing only one standard EWSB phase transition, the 2HDM can undergo a two-stage phase transition as the temperature falls as shown in the lower panel. In this mechanism the first stage is a conventional second order PT in which the symmetry is broken, shortly thereafter follows a first order PT. Another remarkable thing is that the ratio of the classical value between the two fields h 2 /h 1 shown in the upper and middle panels has very little dependence on the temperature. However, in general, the value of h 2 (T )/h 1 (T ) is a temperature-dependent parameter and the change in the temperature growth can even be large in magnitude. In particular, the lower panel displays a peculiar behavior of the ratio h 2 (T )/h 1 (T ) as temperature decreases: at first it monotonically increases, resulting in a deviation of the vacuum from the yellow line, then jumps to its zero temperature value (that is tan β) at the transition point and maintains unchanged in the remaining process.

Numerical procedures: T c evaluation scheme
The dynamics of the EWPT is governed by the effective potential at finite temperature Eq. (3.1) in our model. For purposes of analyzing the temperature evolution of the potential involving both h 1 and h 2 , it is convenient to work with a polar coordinate representation of the classical fields h 1 (T ) and h 2 (T ). To that end, we define h(T ) and θ(T ) via The tree-level potential Eq. (3.2) in the (h, θ) plane becomes 4) and the remaining parts of the effective potential are much more involved and hence not shown here. When a first order PT takes place, a local minimum with h = 0 develops and becomes degenerate with the symmetric minimum h = 0 as the temperature decreases, this defines the critical temperature T c , and the two minima are separated by a potential barrier. Therefore, the evaluation of T c is of great importance in studying the EWPT and its cosmological consequences. A straightforward approach is to decrease the temperature by small steps and make a potential plot (like Fig. 5) at each step. Then the global minimum of the potential (starting from the EW vacuum at zero temperature) can be followed stepby-step and the critical point is found once the minimum displays a jump rather than a smooth transition. Obviously, this graphic method is feasible only for benchmark points but is barely applicable for extensive scan due to its non-numerical nature. To date several numerical methods have been developed. In Ref. [19], going from zero temperature to higher temperature, the critical point is taken to be the last one for which the minimum lied below the origin. This approach is no able to resolve the 2-step phase transitions where the potential experiences a second order PT prior to the first order PT, giving rise to a vacuum shift from the origin at higher temperature. To overcome this problem, the authors of [21] used advanced numerical algorithms to search for the global minimum of the effective potential in the (h, θ) plane for each temperature. We employ a method consisting of the following procedures: First, we deal with points for which the ratio h 2 /h 1 is (approximately) temperatureindependent, that is θ(T ) = β. The effective potential Eq. (3.1) reduces to a function of two parameter -temperature T and T -dependent field norm h(T ). In this case, one can easily determine the critical temperature T c and the field norm v c at which the potential reaches a minimum by solving the equations In searching for the solution of the above equations, we require a difference between the potential at the minimum and its value at the origin smaller than 10 −10 GeV 4 . As a consequence, the solution for v c would be a value close to zero if there were no degenerate minima of the effective potential present in the process of temperature drop. This means that below a certain small value of v c we do not expect a decent probability of achieving a first-order phase transition. Instead, very likely such points lead to a second-order phase transition. For this reason we employ a technical cut v c > 1 GeV in order to remove these points. Next, we are going to deal with the points that exhibit an explicit temperature dependence for the ratio of two fields. This type of points often lead to a 2-stage phase transition [83,84], as illustrated in the last row of Fig. 5. There must exist a global minimum for which tan θ(T ) = tan β is not obeyed at a certain temperature or within a small temperature interval. In this situation, (v c , T c ) obtained as a solution of Eq. (5.5) is not the critical vev and temperature where the phase transition occurs because the true vacuum is no longer located at the origin. Searching for the global minimum should be performed not along the tan β line but on a two-dimensional (h, θ) space. We employ an algorithm which uses the steepest descent method to find the global minimum of the effective potential. At T c the searched minimum is then compared with the value of the effective potential evaluated at v c and the one with the lower value is chosen as the candidate for the global minimum. For the general 2-stage PT, tracing the (temperature) evolution of the global minimum on a 2D plane is inevitable. Here, we discard points leading to a 2-stage PT and focus on the scenarios featuring a 1-stage PT. In the following analysis, we only retain parameter points with T c ≤ 300 GeV. 18

Properties of the first order EWPT
The strength of the phase transition is quantified as the ratio of the norm of the neutral fields to the temperature at the critical point, Here v c = h 1 2 + h 2 2 + A 2 , 19 in general, represents the value of the norm of all scalar fields involved at the broken vacuum at critical temperature T c . Note that when interpreting the ratio as the strength of the electroweak phase transition, one should be aware of its gauge dependence [76,[85][86][87]. In order to ensure that a baryon number generated during the phase transition is not washed out, a strong first-order phase transition is demanded and occurs if ξ ≥ 1 [88]. 20 Before presenting the main results, we discuss the specific features of the parameter space compatible with the theoretical and experimental constraints and at the same time leads to first order and second order phase transition. We will show results for both Type I and Type II models.

first order vs. second order phase transition
It has been shown in Fig. 5 that both first order and second order phase transition can take place in the 2HDM. Whether first order or second order PT is developed depends on the mass spectrum among the three extra Higgs bosons, which is directly related to the five quartic couplings λ i and the soft symmetry breaking parameter m 2 12 through Eqs. (2.8)-(2.10). Thus, it would be very interesting and useful if one can divide the entire model parameter space into different sectors where distinct dynamics of vacuum evolution leading to first order and second order PT take place. An initial attempt along this direction was made in [90] in accordance with the general geometric analysis of [91]. In [90] the authors introduced several discriminators in terms of certain combinations of λ i and succeeded in dividing into four sectors which do not overlap. However, the analysis conducted in [90] is oversimplified, only the effect of thermal mass corrections was included. When considering the full effective potential Eq. (3.1), such division may be highly difficult or even impossible, which is reflected in Fig. 6 where we map our first order (red) and second order (blue) PT points in the 2D space of model parameters and none of the parameters exclusively distinguish the two types of PT points. 21 As expected, λ 1 and λ 2 have marginal influence since both of them only enter into the masses of two neutral CP-even scalars. On the contrary, λ 3,4,5 can be potentially used as discriminators as they affect the masses of 19 Since we restrict ourselves to a CP-conserving model, the global minimum has A = 0. In general there may exist local minima that are CP-violating, while a recent study [21] found that it is always vanishes up to numerical fluctuations at both T = 0 and T = Tc. 20 The choice of the washout factor is subject to additional uncertainties. It was argued that the EW sphaleron is not affected much if extra degrees of freedom are SM-gauge singlets [89] but the situation in the presence of an additional doublet is unclear yet. As a more conservative choice, other criterion such as ξ ≥ 0.7 was also taken in other works. 21 We examined that none of the discriminators defined in [90] can effectively isolate the first order (red) and second order (blue) PT points when considering the full effective potential. three extra Higgs bosons simultaneously. For instance, λ 5 is bounded from -6 to 6 (2) in Type I (II) model and a small value of λ 5 tends to induce a second order PT unless the sum λ 345 is negatively large. Another important observation is that for a given value of tan β, larger m 12 , allowing for larger m H , favours a second order PT. This points to the fact that the phase transition in the theory degrades to the SM case when the new scalars reside in a decoupled sector, as expected intuitively. In reverse, it has the implication that first order PT is more probable for a small or modest value of m 12 when m H is fixed. To illustrate this, we evaluate the phase transition properties in the process of slowly varying m 12 , assuming the alignment limit and a common mass scale among the three BSM states M = 600 GeV for simplicity. The situation is shown in Fig. 7, where red and green curves represent tan β = 1 and 1.5, respectively. This plot can be used to track the evolution of the critical vev v c : it starts from zero (in the second order PT stage) at large m 12 to a non-zero value (in the first order PT stage). The jump from the second order PT to the first order PT is indicated by a dashed line with an arrow. Notably, a severe fine-tuning on m 12 is required for a successful first order PT and the v c value approaches the EW vacuum at smaller m 12 . This interesting behavior is explicitly illustrated in Fig 8 which gives, for tan β = 1, the 1-loop potential curve at zero temperature (left) and the finite temperature effective potential evaluated at the critical temperature (right) for various values of m 12 . As m 12 decreases, thermal effects generate a higher potential barrier and simultaneously push the degenerate  vacuum towards the EW vev v, giving rise to a growth in ξ (owing to the small fluctuation on T c in the stage of the first order PT). On the other hand, a smaller contribution from the m 2 12 term to the tree-level and 1-loop potential at zero temperature will remove the potential barrier. For example, the SM potential V ∼ λh 4 when the mass term µ 2 → 0. Consequently, the desired vacuum disappears, resulting in a terminal value of v c near v, as we will also see in Fig. 9. Furthermore, the effect of increasing tan β on the phase properties is also visible by comparing the red and green curves. For a larger tan β and the same mass spectrum, the first order PT is realised at a lower value of m 12 and in the meanwhile the 'no-EWSB' situation takes place at a smaller value of v c .
As also seen in Fig. 6, most of our points have tan β close to one, which agree well with the findings of previous studies [19]. Yet we would like to clarify that such preference is absolutely not the consequence of requiring a (strong) first order PT. The underlying reason is that in the vicinity of tan β 1, a large range of m 2 12 satisfying the theoretical constraints outlined in Sec. 2.1 is allowed. 22 Oppositely, m 2 12 is strongly constrained in the high tan β region and a fine-tuning is required, which will greatly increase the difficulty of accumulating the points by means of random scan. Numerically, very limited range of tan β is allowed for large m 2 12 .

Properties of the first order EWPT
We now turn to discuss the general properties of the first-order PT accomplished in the 2HDM. The crucial parameters of the phase transition include the critical temperature T c , the field value v c at T c and their ratio ξ = v c /T c which is used as a measure of the strength of the EWPT. In Fig. 9 we display in the (T c , ξ) plane the points consistent with all theoretical constraints on the potential and up-to-date LHC limits at Run-2. Three black contours from top to down correspond to v c = 250, 135, 50 GeV. We first discuss the 22 The correlation between tan β and m 2 12 were discussed in details in Ref. [92].
impact of extra scalars in the spectrum. Suppose all extra scalars are heavy (i.e, above 800 GeV) and thus their masses are highly degenerate required by the EWPD (see Fig. 1), then the sector consisting of the new scalars decouple from the SM Higgs and the dynamics of phase transition behaves like the SM. Of course, the strength of EWPT is not closely related to the masses of any of additional Higgs bosons but more directly linked to the mass splittings among them, which can be explicitly visualized in Fig. 14 presented later. A general tendency observed is that v c is more constrained as T c decreases. In the extreme case of T c 100 GeV, the thermal effect, while still playing the role of lifting the effective potential and forming two degenerate minima, is too weak to compete with the zero-temperature loop corrections to the potential. As a result, the critical classical field value is mostly localized around v, which makes it slowly vary with respect to the temperature change. Nonetheless, v c shown in Fig. 9 does not exceed the zero temperature EW vacuum value v owing to the EW vacuum run away ('no-EWSB' bound) as sketched in Fig. 7, implying that the PT strength ξ necessarily improves at low T c . More quantitively, this leads to a maximum PT strength ξ 5 at T c = 50 GeV, and, on the other hand, implies an upper bound on T c at 250(350) GeV for ξ ≥ 1(0.7) 23 . In addition, we observe that a lower bound on T c for each value of v c . For the value of the critical classical field v c being slightly away from the EW vacuum v, the lower bounds on the critical temperature would be around T c 100 GeV in Type I and the lower bound on T c in Type II model is slightly raised due to the lack of m A ≤ 350 GeV points. We stress that this is an useful finding that one can utilize to greatly optimize the algorithm for the evaluation of T c . Last, we point out that the extremum, if coexisting in the vicinity of v c 135 GeV, often develops to a local maximum (corresponding to a barrier) rather than a local minimum of the potential, which causes a narrow gap dividing the displayed points into two parts.
An explicit dependence of the critical temperature T c on the mass spectrum of the three extra Higgs bosons can be visualized in Fig. 10, where we display all points that pass the applied constraints as in Fig. 9 and additionally fulfill a strong first order EW phase transition (i.e., ξ ≥ 1).
Having explored these SFOPT behaviors, we shall investigate the relation between critical classical field values, critical temperatures and different contributions to the effective potential in the model. While the thermal contribution is crucial in controlling the process of vacuum tunneling, lots of attempts have been made to describe the properties of the phase (i.e. v c and T c ) from the effective potential at zero temperature. A recent progress was reported in Ref. [93] (within the framework of the CP-conserving 2HDM) that the strength of the phase transition is dominantly controlled by the value of F 0 , the depth of the 1-loop potential at zero temperature between the symmetry unbroken vacuum h = 0 and the symmetry broken vacuum h = v which corresponds to, in our notation, is the full 1-loop potential at zero temperature. Using the normalized depth ∆F 0 defined in [93] one can derive an upper bound that definitely guarantees the PT to be strong, for example, ∆F 0 /F SM 0 −0.34 necessarily leads to ξ ≥ 1 in the 2HDM. This, of course, can be used as an empirical test to assess the strength of the phase transition. However, a strong first order PT is still possible even though this upper bound is overflowed, in this situation the thermal potential plays a more important role for the thermal evolution of the system. Therefore, while appreciating the advantage of this approach in simplifying the phase transition study, which allows to find regions of the parameter space where a SFOPT could be achieved, we expect a deeper comprehension by investigating not the strength ξ itself, which is not an intrinsic property of the phase transition, but the characteristic quantities derived from the phase dynamics: v c and T c . Interestingly, we find that the magnitude of v c increases towards the zero temperature VEV with the decrease of the vacuum depth ∆V 1-loop 0 (v) independent of the value of T c . This is illustrated in Fig 11 and is one of the nontrivial outcomes of this work. It is naively true that v c v when |∆V 1-loop 0 (v)| 0, which implies that the thermal effects in the presence of extra scalars enhance the value of the effective potential at the SU(2) symmetry broken vacuum and almost do not shift the symmetry broken vacuum at the critical temperature. As expected, as the vacuum depth |∆V 1-loop 0 (v)| increases, v c decreases towards the classical field value of h = 0, which results in a smaller value of ξ for a given T c . In the meanwhile, we emphasize that the precise evaluation of v c (and T c ), of course, requires the inclusion of the temperature-dependent part in the potential. The critical temperature T c is supposed to be more related to the thermal corrections to the effective potential, as demonstrated in the lower panels of Fig. 11.
In addition to the non-thermal loop effect discussed above, the thermal effect in the presence of extra scalars is another promising source driving the SFOPT. In the 2HDM, extra BSM bosonic states are present in the plasma and induce the additional contribution to the thermal mass through the quartic couplings (λ 1,2,3,4 ), see Eq.  Figure 11. Properties of first order PT in the 2HDM. The z-axis in the upper and lower plots represents the vacuum depth of the zero temperature potential |∆V 1-loop 0 (v)| and the thermal potential in the broken vacuum at the critical temperature V T (v c , T c ), respectively. Only T c 300 GeV points are retained. a proper cancellation between their masses and couplings is satisfied, an energy barrier can be generated so that the PT becomes strongly first order [94]. In order to see the importance of the thermal effect, we estimate the thermal masses for three extra Higgs bosons Eq. (3.22) at critical temperature T c and present in Fig. 12 the ratio normalising the zero-temperature masses (the measured masses) as a function of the PT strength ξ. Clearly, the ratio for the three states have a large variation around 1 on both sides, which means their thermal corrections can be either constructive or destructive even for the SFOPT (ξ ≥ 1). In particular, this ratio for the CP-odd A state in Type I model can be up to ∼ 20 owing to the presence of the extremely light A. While the thermal correction tends to suppress the m A and m H ± at T c , the preference over the enhancement on the H (relative to H ± ) is still visible. The importance of the thermal mass maximizes at ξ 1 and becomes less significant as ξ further grows.
Recall that the SU(2) custodial symmetry is not severely broken at zero temperature due to the T parameter in the EWPD which forces small mass difference for |m H ± − m A | or |m H ± − m H | or both. One may be curious whether this symmetry is broken at finite temperature. This is especially interesting when such symmetry plays a crucial role in selecting a particular region of parameter space. In general, the thermal correction to the field dependent masses might results in a shift of the symmetry of the model at finite temperature. The particular case of interest is the Z 2 symmetry cases studied in Refs. [95][96][97] where the Z 2 symmetry is preserved at T = 0 but spontaneously broken at T = 0. To examine if the effect of thermal corrections leads to a shift of the SU(2) custodial symmetry in our model, we estimate the ratio of the thermal mass for two neutral states with respect to that for the charged state at critical temperature T c . The result is illustrated by Fig. 13 where one can observe that the points displayed are well aligned either m 1 with about 10-20% departure, indicating a large violation of the SU(2) custodial symmetry is not possible at finite temperature during the SFOEWPT in the 2HDM. 7 Strong first order EWPT and the implications for future measurements at colliders 7.1 Typical mass spectra and discovery channels at LHC As seen from Fig. 9, a SFOPT is possible in both Type I and Type II models. Then one may wonder what is the LHC Higgs phenomenology associated with a SFOPT. To answer this question, in Fig. 14 we present in the m A versus m H ± (upper) and m H versus m H ± (lower) planes all points that pass the applied constraints as in Fig. 9 and additionally realize a SFOEWPT (i.e., ξ ≥ 0.7). The values of ξ and the mass difference |m A − m H | are indicated in color scale in the upper and lower panels, respectively. We emphasize again that the EWPD, essentially the T parameter, force the mass differences between the charged Higgs boson and at least one of the extra neutral Higgs bosons to be small and strongly favor mass spectra where the masses of all new scalars are close to each other, in the decoupling limit in particular. This severe constraint on the mass spectra for the non-SM Higgs bosons leads to five benchmark scenarios achieving a SFOPT in Type I model.
They are summarized in Table 4 where the characteristic mass spectra and the main decay modes of the heavier neutral Higgs boson (H or A) with an estimate on the branching ratio are given in each scenario.
We start with the analysis in Type I model. First, the most widely studied mass configuration includes a pseudoscalar A with mass within the range of 400 − 600 GeV accompanied with m H ≈ m H ± 200 GeV [19,20]. In this case, m 2 12 must be relatively small since large m 2 12 tends to reduce the strength of the phase transition. This leads to a special relation among the quartic couplings λ 1,2,3 0 and λ 4 −λ 5 5, meaning that the strength of the phase transition is mainly governed by λ 4 and λ 5 , see also [19]. Dictated by symmetry argument, one can image that the mass spectrum consisting of a light CP-odd state and two highly mass degenerate H and H ± can also lead to a SFOPT, which is reflected by the existence of a bulk of red points at the upper left corner (i.e, m H ± 400 − 600 GeV, m A 200 GeV) in Fig. 14. The situation of the model parameters is opposite due to the flip of mass hierarchy among the three BSM Higgs states. To be specific, m 12 is large as a consequence of large m H , and λ 1,2,3 0 and λ 4 −λ 5 −5. Likewise, a SFOPT (ξ ≥ 1) can be also realized provided that m A and m H ± are close to each other, while both having a large gap relative to m H . Strictly speaking, such condition provides two possibilities for the mass spectra: i) m A m H ± 600 GeV and m H 200 GeV and ii) m A m H ± 200 GeV and m H 600 GeV, which correspond to two isolated red-orange points densely distributed along the diagonal line in the lower panel plot. Deduced from Eqs (2.9) and (2.10) the mass degeneracy between A and H ± states in this scenario restrict λ 4 λ 5 , while an additional coupling λ 3 participates into the potential evolution and influences the phase transition. Apart from these four scenarios that are visible in the low panel plot, the upper left plot in Fig. 14 demonstrates an additional possible scenario that is compatible with ξ ≥ 1 where all three non-SM-like Higgs bosons have similar mass scales at 300 − 350 GeV. This scenario was unfortunately ignored [19,20] or paid less attention [21]. 24 It is also worth noting that in this highly degenerate scenario none of λ i couplings can be close to zero if the first order PT takes 24 One might indeed have believed that large mass splitting among the non-SM Higgs bosons are a necessary condition for the requirement of a strong first-order EWPT.

place.
On the other hand, the allowed mass spectrum that is compatible with ξ ≥ 1 in Type II model is quite simple. As explained in Sec. 2.2, the combination effect from Bphysics observables and EWPD pushes m H ± 580 GeV and simultaneously raises the mass scale for at least one of the extra Higgs bosons. Consequently, many scenarios available in Type I model are eliminated, resulting in an allowed mass spectrum that leads to a strong first-order EWPT being quite restrained: m A m H ± ≈ 600 GeV and a large positive mass gap between m H ± and m H : m H ± − m H > ∼ 300 GeV.
Generally speaking, requiring a SFOPT forces down the mass scale for the new scalars and the preferred ranges for all the scalar masses below 600 GeV, which coincidently approaches to the current lower bound on the charged Higgs mass strongly constrained by the latest measurement of B → X s γ. This means that future improvement on B-physics observables may decisively rule out the success of SFOPT in the Type II 2HDM. Of course, Fig. 14 also informs us that weak first order PT (under the criterion of ξ 0.7) would still be possible even if no additional Higgs bosons were discovered below 1 TeV.
Finally, we briefly discuss the prospects of testing the EWPT at the colliders in accordance with the mass spectrum provided above. In the alignment limit sin(β − α) ≈ 1 we consider, the coupling g hAZ is vanishingly small but the coupling g HAZ ∝ sin(β − α) is enhanced. Hence, the branching ratios for A → ZH and H → ZA as long as kinematically allowed can be substantially large depending on the mass spectrum in the model. These results point towards the observation of the A → ZH and/or H → ZA decay channels would be "smoking gun" signatures of 2HDMs with a SFOEWPT. 25 LHC search prospects for the former decay have been analyzed and proposed as a promising EWPT benchmark scenario in [20], while the collider analysis looking at both decays was performed in Ref. [98] but not specifically aiming at the EWPT. In Fig. 15 we show the 13 TeV cross sections at the LHC for these two channels in the gluon-fusion production mode. In all cases, a cross section above the pb level can be achieved for the scenarios realizing a SFOPT. Although these signatures are characteristic ones in most of the 2HDM scenarios discussed above (see Table 4), no strong correlation in these channels is found between ξ and the corresponding cross sections, which means that there is no guarantee to observe these decays in colliders. We leave a detailed collider analysis to future studies.
Searching for a new scalar resonance is performed at the LHC mostly through its decay into SM particles. These decay channels include H → ZZ → 4 , H, A → γγ, τ τ, tt. For the purpose of testing the EWPT, it would be very useful to find channels with strong correlation to the ξ value. The one served as an example here is the gluon-fusion production cross section of A and H in the τ τ decay channel, which is shown in Fig. 16. In general, the gluon-fusion cross section in Type I model is considerably small, so there is very little hope to ever observe A or H in this channel. An exception occurs in the very light CP-odd A region with cross-section as large as the level of 10 − 100 pb [74]. Moreover in this region, m A ≤ 60 GeV, a few points with large ξ values are observed, which could be excluded by 25 Although our results confirm the results in the earlier literature [20], more importantly, we clarify that the decay A → ZH is not a unique "smoking gun" signature of SFOPT in the 2HDM of Type I model. This conclusion is also supported by another recent study [21]. the upcoming experimental searches in that channel. In Type II the situation is different, for a given scalar mass the achievable cross-sections have a lower bound. The large ξ points are located at low m H 350 GeV and have reasonably large cross-sections just below the current experimental upper limit. In short, we estimate that a factor of 4 improvement in the search sensitivity, which is very likely to be reached, would either see an exciting signal or eliminate these points, as a result, the first order PT with strength ξ > 3 can be fully tested at the LHC.

Triple Higgs couplings and the implications of the future measurements
The scenarios that lead to the first order PT in the model have a mass spectrum below the TeV scale, as can be seen in Fig. 14 and Table 4. The presence of additional scalars that couple to the SM-like Higgs h can modify the triple Higgs coupling hhh at both tree-level and loop-level and thus leads to the deviation with respect to its SM value g SM tree hhh . Moreover, such deviation can be significant near the alignment limit provided being away from the decoupling limit [28,99]. We examine both the tree-level coupling and the one after the inclusion of the one-loop corrections. They are computed by taking the third derivative of the tree level potential V 0 and the one-loop potential V 0 +V CW +V CT with respect to h, respectively and shown in the upper and lower panels of Fig. 17   observe that the triple SM-like Higgs self-coupling g hhh in favor of the highly strong PT (i.e, ξ 3) is close to its SM value g SM tree hhh , while large deviation (mostly suppression) of g hhh from g SM tree hhh is possible for the weakly strong PT (i.e, ξ 1.5). Another transparent observation is that the hhh coupling at tree-level cannot be enhanced in Type I (for m H 600 GeV) and Type II models, see Ref. [28] for analytical understanding of these features. However, we stress that this conclusion will be dramatically changed when the one-loop corrections to the hhh coupling are taken into account. As shown in the lower panel plots, the coupling g hhh at one-loop level are absolutely enhanced in both models and the largest normalized coupling g hhh /g SM tree hhh can be about 2.5, corresponding to ∼150% enhancement. This allow us to conclude that the strong PT (ξ ≥ 1) in the 2HDM typically induces the enhancement on the hhh coupling. Next, we would like to quantitatively explore the relation between the phase transition strength and the content of the derivation the triple Higgs coupling. In general, the loop-level hhh coupling exhibits a larger deviation with increased strength of the phase transition. Whereas, the tree-level hhh coupling shown in the upper panel plots does not display such a proportionality behavior. This dramatic change implies that the loop corrections coming from the CW potential and counter-terms are important in general when the phase transition is of strong first order and can even be dominant over the tree-level contribution in the case of the extremely strong phase transition. It is also apparent in Fig. 17  enhancement on the hhh coupling. In contrast, the hhh coupling normalized to the SM tree value can vary from ∼ 1 to 2.5 for the weakly strong PT of ξ = 1 − 2. This means that large triple Higgs coupling hhh is a necessary but not sufficient condition of realizing the highly strong PT. For instance, if the deviation is smaller than 100%, then possibility of the highly strong PT (ξ 3 (2.5) in Type I (II)) will be eliminated. As a result, the size of the triple Higgs coupling hhh derives an upper bound on the achievable value of ξ. In some sense, this is phenomenologically useful because we have built a connection between the phase transition involving the thermal contribution and a measurable observable at zero temperature. Therefore, the measurement of the triple Higgs coupling could be an indirect approach of probing the phase transition at colliders.
Experimentally, the deviation of the triple Higgs coupling can be detected at both lepton colliders (i.e., ILC [100], CEPC [101] and FCC-ee [102] ) and hadron colliders such as LHC and SppC [103]. At hadron colliders, the resonant Higgs pair production is promising while special attention needs to be paid when the heavier CP-even state H produces a destructive interference with the SM top box diagram process [104]. Upon the sensitivity of 50% supposed to be achieved at the HL-LHC, a large amount of the (nearly entire) parameter space in Type I (II) model leading to strong PT can be probed through the di-Higgs production into bbγγ and bbW + W − channels in the ultimate operation of LHC Run-2 [105][106][107]. In our case, g hhh has the same sign as the SM value and hence results in the destructive interference between the s-channel h-mediator triangle diagram and the top box diagrams of the gg → hh production process. This implies increasing g hhh will decrease the production cross section [108]. Previous studies demonstrated that when g hhh 2.45g SM tree hhh an exact cancellation between these two diagrams is accomplished at the threshold of the di-Higgs invariant mass m hh = 2m t [104,109,110]. Due to the low acceptance at LHC for large g hhh , a cut m hh < 2m t is imposed [104,109]. MVA analysis of Ref. [109] shows that, for the parameter space leading to the SFOPT (presented in Fig. 17), the observation significance in the bbγγ channel with the integrated luminosity of 3 ab −1 at 14 TeV would decrease from 10 to 4 in both Type I and Type II models. In measuring the triple Higgs coupling hhh the lepton machines are typically more powerful, using the Higgs associated process e + e − → Z * → Zh * (hh). The better designed sensitivities at the CEPC [101], FCC-ee [102] and ILC1000 are roughly 20-30%. This indicates that almost the full parameter spaces that are compatible with ξ ≥ 1, particularly for m H 500 GeV, are within the future detection reach.
The other Higgs self-coupling of interest is the Hhh coupling g Hhh , which is also relevant to the Higgs pair production through the s-channel H mediator triangle diagram. The Hhh coupling at one-loop as a function of m H is depicted in Fig. 18. In contrast to the hhh coupling, the one-loop corrections to Hhh coupling are vanishingly small near the alignment limit. We thus do not show the tree-level result. It is important to mention that the Hhh coupling can be significant even in the alignment limit, which can be observed in Fig. 18. For instance, the Hhh coupling is about ±(30 − 50)% of the tree-level SM hhh coupling for the highly strong PT (ξ ≥ 3) and can be even comparable with or larger than g SM tree hhh as the PT is weakly strong (ξ ≈ 1 − 2). Notably, the obtained Hhh coupling g Hhh in the successful SFOEWPT scenarios can have either the same sign as or the opposite sign to the coupling g hhh . The consequence of the sign flip of the g Hhh will affect the s-channel H-resonant triangle diagram contribution to the gg → hh process, whose amplitude is proportional to the product of g Hhh and g Htt = C H U y t , resulting in a change on the m hh lineshape due to the interference between the triangle diagram of the signal and the continuum top box diagram. When the interference is destructive, special attention needs to be paid [104]. The study of Ref. [111] indicates that most of our SFOPT points can be detected at 5σ significance provided that g Hhh × g Htt > 300 GeV at 14 TeV with integrated luminosity of 3 ab −1 using the bbγγ channel.

Conclusions and Outlook
Taking into account theoretical and up-to-date experimental constraints, we studied the electroweak phase transition in the framework of the CP-conserving 2HDM of Type I and Type II models near the alignment limit. The thermal potential was expressed in terms of modified Bessel functions, which allows for a fast numerical evaluation and high precision compared to the simpler high/low temperature approximations. While both 1-stage and 2-stage phase transitions were shown to be realized within the 2HDM, in this paper we focused on scenarios leading to 1-stage phase transitions at electroweak scale, for which first order and the second order phase transitions are distinguished.
We analyzed the properties of the first order phase transition, observing that the field value of the electroweak symmetry breaking vacuum at the critical temperature is strongly related to the vacuum depth of the 1-loop potential at zero temperature, while the critical temperature reflecting the size of the thermal effect is characterized by the temperaturedependent potential. In general, the critical temperature T c tends to be higher as the BSM states becomes heavier, and on the other hand T c can be down to ∼ 100 GeV when at least one light BSM Higgs bosons present in the spectrum. We have also observed that the thermal correction to the mass is important in driving a SFOPT.
The strength of the transition, a key property for the electroweak baryogenesis mechanism, depends largely on the allowed mass spectrum. Requiring a SFOPT with ξ ≥ 1 forces down the mass scale for the new scalars and the preferred ranges for all the scalar masses below 600 GeV. We demonstrate that SFOPT (i.e., ξ ≥ 1) required for baryogenesis is possible in both Type I and Type II models. In Type I model, SFOPT is achievable in the parameter space where a large mass splitting is present between two neutral Higgs bosons such as m H m A and m A m H . In either case, the charged Higgs mass is close to either m H or m A required by the EWPD. The mass spectrum among the extra Higgs bosons in the Type II model is, on the contrary, strongly constrained due to flavor observables, which push the mass of the charged Higgs above ∼ 600 GeV. As a result, scenarios leading to a SFOPT in Type II are m H ± m H m A and m H ± m A m H . In view of large mass splitting between H and A, both pp → H → ZA and pp → A → ZH can be "smoking gun" collider signatures related to a SFOPT in the 2HDMs as the cross sections via gluonfusion production in these two channels predicted for SFOPT points are typically up to ∼ 1 pb. In addition to large mass splitting, SFOPT can also take place in Type I even if all the masses of the three extra Higgs bosons (A, H and H ± ) are degenerate around 350 GeV. Such scenario leads to potentially testable consequences through the A → Zh decay channel at colliders.
Following the analysis of the benchmark scenarios, we investigated the implications of a SFOEWPT on the LHC Higgs phenomenology. Various characteristic collider signatures at the 13 TeV LHC have been identified, among which the gluon-fusion production cross section of A and H in the τ τ decay channel displays a correlation with the PT strength ξ.
It turns out that new physics searches at collider machines can provide an indirect channel to examine the EWPT scenarios. Finally, we verify that an enhancement on the triple Higgs coupling hhh (including loop corrections) is a typical signature of the SFOPT driven by the additional doublet. The PT with larger strength is associated with larger deviation of the loop-level triple Higgs coupling hhh with respect to the SM value, which can help to enhance an energy barrier. Meanwhile, we notice that the other triple Higgs coupling g Hhh can also be comparable with the triple Higgs coupling in the SM for SFOPT so that the search for the heavy neutral Higgs H through the gg → hh process is possible for small tan β since the top Yukawa coupling of the H is proportional to cot β.
We leave for future work the interplay of gravitational waves signals and testable colliders signatures for SFOPT benchmark scenarios presented in this paper. This success would build a link between early Universe cosmology and collider detection, which could provide additional constraints in the allowed parameter space of the 2HDM. We believe that such connection will have a significant physical value and serves as a useful guide for collider search strategies.