Possibility of multi-step electroweak phase transition in the two Higgs doublet models

We discuss whether a multi-step electroweak phase transition (EWPT) occurs in two Higgs doublet models (2HDMs). The EWPT is related to interesting phenomena such as baryogenesis and a gravitational wave from it. We examine parameter regions in CP-conserving 2HDMs and find certain areas where the multi-step EWPTs occur. The parameter search shows the multi-step EWPT prefers the scalar potential with the approximate $Z_2$ symmetry and a mass hierarchy between the neutral CP-odd and CP-even extra scalar bosons $m_A<m_H$. By contrast, the multi-step EWPT whose first step is strongly first order favors a mass hierarchy $m_A>m_H$. In addition, we compute the Higgs trilinear coupling in the parameter region where the multi-step EWPTs occur, which can be observed at future colliders. We also discuss a multi-peaked gravitational wave from a multi-step EWPT.

In this paper, we study the multi-step EWPTs in the CP-conserving 2HDMs. Because of the absence of the new CP-violating source, the EWBG does not work and we would not discuss it. The study for the CP-violating case remains as future work. The main purpose of this paper is to reveal features of the multi-step PTs. By performing parameter searches, we find certain parameter spaces where the multi-step PT occurs. Furthermore to examine the possibility of verification of the multi-step PT at collider experiments, we compute the deviation of the Higgs trilinear coupling from that in the SM. It is known that the deviation can be large in the 2HDMs [71,72] (see also Ref. [73,74] for recent work), and it would be observed more precisely in future colliders like High-Luminosity Large Hadron Collider (HL-LHC) [75] and International Linear Collider (ILC) [76]. We find that the deviation has a tendency to be large in certain regions when the multi-step PTs occur. In addition, we calculate a two-peaked GW spectrum yielded by a 2-step PT, which can be observed by using LISA, Big Bang Observer (BBO) [77], and Ultimate Deci-Hertz Interferometer Gravitational Wave Observatory (U-DECIGO) [78].
The outline of this paper is as follows: In section II we introduce the generic characteristics of the 2HDMs. Section III is dedicated to give the thermal effective potential. Theoretical constraints considered in our numerical analyses are briefly introduced in Section IV. In section V we show the results of the parameter search for the multi-step PT. Moreover, in section VI, we discuss the predictions for the Higgs trilinear couplings as the collider signatures and for the multi-peaked GW as the cosmological signature for the multi-step PT. Our conclusions are given in section VII.

II. TWO HIGGS DOUBLET MODEL
The tree-level scalar potential of the CP-conserving 2HDMs with a softly broken Z 2 symmetry is written as where Φ i (i = 1, 2) are the SU (2) scalar doublets We here assume that only the neutral CP-even scalar fields have the vacuum expectation values (VEVs) v i , which are real and positive, and satisfy v ≡ v 2 1 + v 2 2 = 246 GeV. The third term with m 2 3 on the right-hand side in Eq. (1) breaks the Z 2 symmetry in the potential softly. The coefficients are taken to be real, although m 2 3 and λ 5 are complex parameters in general. Regarding only the neutral CP-even fields φ i , the Φ i become with A = m 2 3 tan β + λ 1 v 2 cos 2 β, Throughout this paper we take h as the SM-like Higgs boson with m h = 125 GeV. The field-dependent masses of the W boson, the Z boson, and the photon can be written as where g and g are the gauge couplings of SU (2) L and U (1) Y gauge symmetry, respectively. The physical masses of the gauge bosons are derived by taking φ i = v i .
The most general Yukawa term is where Q L and L L are SU (2) L doublets of quarks and leptons, respectively, Y f (f = u, d, l) are the Yukawa matrices of the fermions, and each of Φ f is either Φ 1 or Φ 2 . Since the 2HDMs have the two SU (2) scalar doublets, we assume one of the doublets couples each of the fermions to avoid the tree-level flavor changing neutral current. One of ways to accomplish this is assuming 2HDMs have a Z 2 symmetry. In this case, there are 4 types in the 2HDMs distinguished by the Z 2 charges for each of the fermions as in the Tab. I [79][80][81]. In the Type-I 2HDM, all quarks and charged leptons obtain their masses from the VEV of Φ 2 . In the Type-II 2HDM, the VEV of Φ 2 gives the masses of the up-type quarks, while that of Φ 1 provides those of the down-type quarks and the charged leptons. In the Type-X 2HDM, the charged leptons and quarks obtain their masses from the VEV of Φ 1 and Φ 2 , respectively. In the Type-Y 2HDM, the masses of the down-type quarks are generated by the VEV of Φ 1 , while those of the up-type quarks and the charged leptons are obtained by that of Φ 2 . The field-dependent masses of fermions can be described as where which value assigned to i depends on the types of Yukawa interactions. The physical masses of the fermions are obtained by taking φ i = v i .

III. THE EFFECTIVE POTENTIAL AT FINITE TEMPERATURE
A. The one-loop corrected effective potential An EWPT is caused by the temperature change of the effective scalar potential. To study the PT, we consider the thermal effective potential. The one-loop corrected effective potential at the finite temperature V β is where V 0 , V CW , V CT , and V β 1 are the tree-level potential (4), the one-loop level potential at zero temperature (the Coleman-Weinberg potential), the counterterm potential, and the one-loop level potential at the finite temperature, respectively.
The Coleman-Weinberg potential in the MS scheme is written by [82] where k indicates scalar and gauge bosons and fermions, and n k , m k , and µ are the degrees of freedom of each fields, the field-dependent masses of each fields, and the renormalization scale which we set µ = 246 GeV, respectively. The upper (lower) sign corresponds to the bosonic (fermionic) contribution. The corresponding degrees of freedom are n k = 2, 1, 1, 1, 6, 3, 2, 12, 12, and 4 for k = H ± , H, h, A, W, Z, γ, t, b, and τ , respectively. We only consider the fermions which have the non-negligible contributions. The constant c k are equal to 1/2 for transverse gauge bosons and 3/2 for the other particles in the MS scheme. The V CW changes the coordinate of the global minimum of the potential from that of V 0 . We introduce the counterterm potential V CT for fixing the coordinate, the masses and the mixing angles of the scalar fields to be equal to the tree-level ones. Thus, we impose the following five conditions to determine V CT : Following Ref. [16], we set V CT with five parameters δm 2 1 , δm 2 2 , δλ 1 , δλ 2 , and δλ 345 , Hence, the conditions (20) give where . We calculate δm 2 1 , δm 2 2 , δλ 1 , δλ 2 , and δλ 345 numerically and substitute them for V CT . However, there are infrared divergences in the second derivatives of V CW , which are proportional to log m 2 NG where m NG indicate the masses of Nambu-Goldstone (NG) bosons. To avoid these divergences, we use the approximation which is shown in Ref. [83]. In this approximation, m NG are approximated as the mass of the SM-like Higgs boson, i.e. m NG → m h . This approximation is justified because the divergences are only logarithmic, hence the changes of the masses of the NG bosons do not make large differences.
The one-loop thermal contributions to the potential can be written as [84] where T represents the temperature and the upper (lower) sign indicates the bosonic (fermionic) contribution. We calculate the integral in Eq. (23) numerically. The squared-masses of the scalar bosons in Eq. (23) can become negative for certain sets of the coordinate (φ 1 , φ 2 ) and T 4 . In that case, we adopt a method that is discarding the imaginary part of the thermal potential, which is related to the instability of the field configuration [85], and taking only the real part (e.g. Ref. [15]).

B. Resummation
Although V β contains the corrections to the one-loop level, the contributions of higher loop diagrams get larger as the temperature rises. The dominant diagrams at the high temperature are called daisy diagrams [84]. We perform resummation which is the method for taking into account the corrections from the diagrams [86,87]. Although there are two methods for the resummation, we apply the Parwani method [86] 5 . The resummation is achieved by appending the corrections from the scalar and gauge boson polarization tensors in the infrared limit Π B (T ) to the masses of bosons m 2 and inserting these corrected masses to V β 1 in Eq. (23) [89]. The index B represents the species of bosons.
In the 2HDMs, we carry out the resummation concretely as the following. The resummation for scalar fields are performed by adding the contributions of the two-point functions to the mass parameters m 1 and m 2 in the mass matrices Eqs. (7) and (9) [90] where c i are the coefficients of correction terms and determined by Π B (T ) which depends on the types of Yukawa interactions. In the Type-I 2HDM, they can be written by [16] For the other Yukawa types, one can obtain the coefficients by apportioning the Yukawa coupling terms in Eq. (26) to c 1 and c 2 according to Tab. I. And then, we append the correction terms to the non-diagonalized scalar matrices M 2 w ± , M 2 z , and M 2 h in Eqs. (7) and (9) as and obtain the corrected scalar masses by diagonalizing them. For the gauge fields, one can carry out the resummation by appending the contributions to only the longitudinal component of the mass matrices. Following Ref. [90], the corrected masses of the longitudinal W boson can be written by The corrected mass matrix of the longitudinally polarized Z boson and photon in the gauge basis is By diagonalizing it, the corrected masses of the Z boson and the photon are obtained as with

IV. THEORETICAL CONSTRAINTS AND EW-VACUUM STABILITY
For the theoretical constraints on the model, we consider constraints from the boundedness from below (BFB) of V 0 , which is described as [91][92][93][94] the perturbativity, |λ n | < 4π (n = 1, 2, · · · 5), (33) and the tree-level unitarity [95,96]. Furthermore the absolute tree-level stability of the EW vacuum is required [97,98], where the negative m 2 3 is disfavored. In following analyses, we confirm numerically that the EW vacuum is the global minimum in the region for |φ i | ≤ 10 TeV, and remove the cases where φ 2 1 + φ 2 2 = 246 GeV is not satisfied at the global minimum.

V. NUMERICAL RESULTS
In this section, we discuss the parameter space where the multi-step PT occurs. To study the PT, CosmoTransitions [99] is used in the analyses. We also study the region where the strongly first order PT occurs in the multi-step PT. The strength of the PT ξ is defined by where T c is the critical temperature, at which minima degenerate between two phases, and v c is the critical value of φ 2 1 + φ 2 2 at T c . As the criterion for a strong PT, we consider ξ ≥ 1, where the sphaleron processes are suppressed enough in the SU (2) broken phase. We especially focus on the cases that the first step PTs of the 2-step PTs are strongly first order 6 , which we name "the strong 2-step PTs" 7 . In this case, the sphaleron rate is suppressed in the broken phase, and it is expected that v c and T c at the subsequent step PTs are respectively larger and smaller than the previous PT. Therefore, the inequality ξ ≥ 1 is kept and the sphaleron rate is also suppressed at the subsequent step PTs 8 .
In the following analysis, instead of the eight parameters (m 2 1 , m 2 2 , m 2 3 , λ 1−5 ) in V 0 , we take the following set as the input parameters: Here m h = 125 GeV and v = 246 GeV. In order to put restrictions on the range of the other input parameters, we consider the experimental constraints from the EW precision data, the B → X s γ decays, the H ± → τ ν decays, and the coupling measurements of the Higgs boson. The EW precision data can be satisfied by assuming the mass degeneracy between the charged scalar boson and at least one of the extra neutral scalar bosons, m H ± m A or m H , which makes the custodial symmetry recovered and hence the ρ parameter ρ 1 like in the SM [100]. For m H ± , the range m H ± < 590 GeV is excluded from B → X s γ decays in the Type-II and -Y 2HDMs [101], while m H ± ≤ 170 GeV is excluded from H ± → τ ν decays in the Type-X 2HDM [102]. The constraints from the coupling measurements of the Higgs boson [103] show that e.g., is excluded at tan β = 2 (10) in the Type-I 2HDM, and | cos(β − α)| > 0.15 (0.05) is excluded at tan β = 2 (10) in the Type-X 2HDM. In the Type-II and -Y 2HDMs, the constraints are stricter than those in the Type-I and -X 2HDMs. In our analyses, based on the above constraints, we take the range for the input parameters in the Type-I and -X 2HDMs as shown in Tab TeV for the other extra neutral scalar boson. For the mixing angles, we take tan β = 2 − 10, and | cos(β − α)| ≤ 0.25 in the Type-I 2HDM, while the alignment limit, cos(β − α) = 0, in the Type-X 2HDM. In the analyses we focus on the relatively small value for m 3 as 0 ≤ m 3 ≤ 100 GeV, since the strong 2-step PTs which we are interested in prefer to occur for the smaller m 3 and do not occur for m 3 100 GeV (cf. Fig. 6). As depicted in Tab. II, we take every 10 GeV in m A and m H , 0.5 in tan β, 0.05 in cos(β − α) (though we set cos(β − α) = 0 in the Type-X 2HDMs), and 5 GeV in m 3 .
In the Type-II and -Y 2HDMs, we take the same ranges for the input parameters with those in the Type-X 2HDM, but m Φ = m H ± ≥ 590 GeV by the constraint from B → X s γ. In these cases, we have found that the stability of the EW vacuum is not realized because the contributions of the heavy extra scalar fields lift up the potential significantly at the EW vacuum (see Eq. (19)) and the origin (φ 1 , φ 2 ) = (0, 0) becomes the global minimum. Hence we discuss only Type-I and -X 2HDMs hereafter.
In this subsection, we show the results in the Type-I 2HDM with m A = m H ± . Fig. 1 represents the allowed parameter region by the theoretical constraints (the BFB, the perturbativity, and the tree-level unitarity) in the m A vs. m H plane at tan β = 2 (left) and 7 (right). It shows in the case of tan β = 7 the upper limit on m H is lower than that in the case of tan β = 2, as m H 290 (440) GeV for tan β = 7 (2).
The left panels of Fig. 2 exhibit the parameter points where the 1-step and multi-step PTs (left) occur in the m A vs. m H (top), m A vs. tan β (middle), and m A vs. cos(β − α) (bottom) planes. The yellow, blue, and purple points show the results for the 1-step, 2-step, and 3 or more step PTs, respectively 9 . Here in addition to the theoretical constraints considered in Fig. 1, the constraint for the stability of the EW vacuum is further imposed. Compared the top left panel in Fig. 2  Left: Parameter points where the 1-step and multi-step PTs occur in the m A vs. m H (top), m A vs. tan β (middle), and m A vs. cos(β − α) (bottom) planes in the Type-I 2HDM with m A = m H ± . The yellow, blue, and purple points show the results for the 1-step, 2-step, and 3 or more step PTs, respectively. Right: Number of points where the 1-step and multi-step PTs occur as a function of m A − m H (top), tan β (middle) and cos(β − α) (bottom). The 1-step, 2-step, and 3 or more step PTs are colored by yellow, blue, and purple, respectively. The grey dashed lines in the panels represent R multi , which are the ratios of the number of points for the multi-step PTs to that for all PTs. the constraint. In the left panels of Fig. 2, the range of m A where the multi-step PTs occur gets larger, as m H increases, and tan β and cos(β − α) (except for tan β 10 and cos(β − α) 0.25) respectively decreases. In the top left panel of Fig. 2, the parameter region of the multi-step PTs overlaps with that of the 1-step PTs, but for the region where m H 420 GeV with m A 550 GeV or 200-300 GeV, the 3 or more step PTs occur mostly. The right panels in Fig. 2 represent the number of points for the 1-step (yellow), 2-step (blue), 3 or more step PTs (purple), respectively, as a function of m A − m H (top), tan β (middle) and cos(β − α) (bottom). The ratios R multi in the panels, plotted as grey dashed lines, are the ratios of the number of points for the multi-step PTs to that for all PTs, R multi = # of points for the multi-step PTs # of points for all PTs .
We see that, in the top right panel of Fig. 2, the multi-step PTs favor m A − m H < 0, and the ratio R multi 1 is obtained at m A − m H −210 GeV. Hence, when R multi becomes around 1, m H 390 GeV and m A is around 200 GeV. Moreover, the middle and bottom right panels of Fig. 2 show that R multi becomes larger for the smaller tan β and cos(β − α), and reaches about 10% at tan β = 2 and cos(β − α) = −0.25, respectively. For tan β 10 and cos(β − α) 0.25, we can see that R multi are only a few %, respectively, although the allowed ranges of m A for the multi-step PTs are wide in the middle and bottom left panels in Fig. 2.
As expected from Fig. 1, the regions with m H 420 GeV are realized for tan β 2 in this analysis. In such a low tan β case, the B → µ + µ − process gives the constraint on m H ± as e.g., m H ± > 340 (125) GeV at tan β 2 (3) at 95%CL in the Type-I and -X 2HDMs [101]. Therefore the region where the 3 or more step PTs occur mostly with m H 420 GeV and m A (= m H ± ) 200-300 GeV is excluded by the constraint from B → µ + µ − . We have found that even when the constraint is taken into account in this analysis, the multi-step PTs favor the mass hierarchy m A < m H and e.g.,

the region would be constrained by the extra Higgs boson search
H → AZ at the LHC [104][105][106]. It is generally more severe for the low tan β and the cos(β − α) closer to zero [107,108]. We leave the detailed analyses including the constraints from such extra Higgs boson searches for future work.
The left panels of which are shown by green dashed lines. It is notable that in the top right panel of Fig. 3 the strong 2-step PTs favor m A − m H > 0 and R st2 becomes R st2 = 100% at m A − m H 210 GeV. Moreover, the strong 2-step PTs are likely to occur at the small tan β and cos(β − α), respectively, as shown in the middle and bottom right panels of Fig. 3. The parameter regions of the strong 2-step PTs do not receive the constraint from B → µ + µ − , while some of them with m A − m H > m Z could be constrained from the A → HZ decay at the LHC [107,108].
The multi-step PTs tend to occur for the smaller tan β and cos(β − α). Next, we investigate the correlations between these parameters and m 2 2 that is an important parameter to determine the path of PT. As we will see later in Fig. 5, the VEVs after the first step in the multi-step PTs have a tendency to be located mainly along the φ 2 axis. Therefore we expect that m 2 2 should be negative and have large enough magnitude to make the multi-step PTs occur because such m 2 makes the potential decrease in the direction of the φ 2 axis. From Eq. (6), m 2 2 can be written as In parameter region. Hence, as m H | sin α| increases, the negative m 2 2 with the large magnitude can be obtained. As shown in Fig. 1, the larger m H is allowed for the smaller tan β. On the other hand, | sin α| increases as cos(β − α) gets smaller in our parameter space. Thus, the minimum value of m 2 2 decreases as tan β and cos(β − α) get smaller, respectively. Meanwhile, at cos(β − α) 0.2, sin α is possible to be zero. In this case, the second and last terms in Eq. (38) vanish and the value of m H does not affect m 2 2 . However, in the region where cos(β − α) 0.25 and tan β 10, | sin α| can be large to some extent and the contribution of the last term in Eq. (38) recovers, which leads to the negative m 2 2 with the slightly large magnitude. This case is presented in Fig. 4. As described above, when m 2 2 is negative with the large magnitude, the first step in the multistep PT tends to occur along the φ 2 axis. In order to see more clearly, we show in Fig. 5 the VEVs after each step of the 2-step (left) and 3-step (right) PTs in the φ 2 vs. φ 1 plane. The dark-green points indicate the VEVs after the first step PTs and the green points show the ones after the second step PTs. In the right panel of Fig. 5, the light-green points indicate the VEVs after the third step PTs. The EW vacuum, where the relation φ 2 1 + φ 2 2 = 246 GeV is satisfied, is located on the pink line. One can see that the dark-green points have a tendency to be along the axes, and the most of VEVs after the last step PTs are located in the direction of the EW vacuum. In the left panel of Fig. 5, the VEVs after the first step PTs extend to in the directions of the φ 1 or φ 2 axis, even though many of them extend to the φ 2 axis. We find that those cases have the negative  Additionally, it is found that the VEVs after the first step PTs are located on or near the φ 1 axis when tan β 2, cos(β − α) is large and m 3 is near zero. On the other hand, for the 3-step PTs in the right panel of Fig. 5, the VEVs after the first step PTs favor going along the φ 2 axis rather than φ 1 . The magnitudes of those VEVs after the first step PTs are not so large, therefore it would be difficult for the first step PT of the 3-step PTs to become the strongly first order. Similar results to above are also seen in the Type-I 2HDM with m H = m H ± and the Type-X 2HDM with m A or m H = m H ± .
Finally, the left panels of Fig. 6 show the strength of PT ξ of the first step PT as a function of the m 3 10 . Considering the tree-level potential V 0 in Eq. (4), the large m 2 3 makes the potential decrease in the region far from the axes and the magnitude of the negative m 2 2 small (cf. Eq. (38)). These make the directions of the first step PT toward the region far from axes, so that it is difficult for the multi-step PTs whose first step PTs occur along the axes to happen. From the analysis, we have found that the maximum magnitude of the VEVs after the first step of the multi-step PTs are gradually larger as m 3 gets smaller. Therefore, ξ of the first step PT has a tendency to be large for the smaller m 3 . We have also found that the larger R multi and R st2 are obtained for the smaller m 3 as in the right panels of Fig. 6. Note that no parameter points where the strong 2-step PTs occur were found in m 3 > ∼ 90 GeV. We summarize in Tabs. III and IV in Appendix. B, the values or ranges of input parameters where the ratios R multi and R st2 have the maximum values, respectively, for Type-I and -X 2HDMs.   Fig. 8. We can also see that they happen only in tan β 4 and cos(β − α) 0 from the right panels of the second and third lines. In addition, the small m 3 is favored when the strong 2-step PTs happen in the bottom right panel. Some parameter points in the above region are excluded by the constraint from B → µ + µ − , e.g. most of the points for the multi-step PTs with m H 330-340 GeV. Besides, we have clarified that the ratio R st2 is the largest in m A − m H > 0, the small tan β, cos(β − α), and m 3 , respectively, as shown in Tab. IV of Appendix B. These tendencies are not changed even if we consider the constraint from B → µ + µ − . In the cases of the Type-X 2HDMs, we take the alignment limit cos(β − α) = 0. PTs increases as m 3 decreases and reaches around 1.2. From our analyses, R multi in the Type-X with m A = m H ± has the largest value as 21% at m A − m H −130 GeV. It also gets the maximum value at tan β 2 and m 3 0, respectively. We have found that the region for the multi-step PTs with m A − m H −130 GeV is realized for tan β 2 in this case, so that such a region is excluded by the constraints from both the B → µ + µ − and H → AZ [108]. In the allowed region by B → µ + µ − constraint, we have also found that the R multi has the larger value for the m A − m H < 0, e.g. R multi = 9% for m A − m H −90 GeV, and the smaller tan β, respectively.
On the other hand, the strong 2-step PTs only occur in the narrow region with m A > m H as shown in the top right panel of Fig. 9. Moreover, the bottom right panel of Fig. 9 shows that the strong 2-step PTs happen only in m 3 15 GeV and they predict ξ 1-1.3. Although the points for the strong 2-step PTs do not receive the constraint from B → µ + µ − , those with m H 150 GeV would be excluded by the A → HZ search [108], while those with m H 330 GeV remain.
As above, the parameter region where the multi-step PTs occur in this case gets narrow com- In this case we have found that there are no parameter points where the more than 3-step PTs occur, therefore the purple points show the points for only the 3-step PTs. Although the largest value of m A for the multi-step PTs is larger than that in the Type-X with m A = m H ± , the other tendencies in Fig. 10 are similar. The maximum value of ξ for the multi-step PTs reaches near 2 in the bottom left panel. We have found that R multi is the largest at m A −m H −210 GeV, tan β 2, and m 3 0, respectively (cf. Tab. III of Appendix B). However, the region with the large magnitude of the negative m A − m H , which is found at tan β 2, is excluded by the constraint from H → AZ decay [108]. On the other hand, in the region where the H → AZ channel does not open, we find that R multi is obtained for Meanwhile, from the top right panel in Fig. 10, same as the other cases, the strong 2-step PTs only occur when the mass hierarchy m A > m H exists with tan β 2. Additionally, they happen when m 3 is small as m 3 30 GeV in the bottom right panel of Fig. 10. The constraint from B → µ + µ − excludes the part of the region where the multi-step PTs happen, e.g. 290 GeV m H 340 GeV, hence the region for the strong 2-step PTs is excluded.
We have confirmed that the result for cos(β − α) = 0 in the Type-I 2HDM with m H = m H ± have same tendency with the ones in the Type-X 2HDM with m H = m H ± .
To summarize briefly, the region where the multi-step PTs is likely to occur is where m A − m H is negative with large magnitude, tan β is small, cos(β − α) is negative and small in the Type-I 2HDMs (it is fixed at zero in the Type-X 2HDMs), and m 3 is small, respectively. Different from the feature of the multi-step PTs, the strong multi-step PTs occur only when the mass hierarchy m A > m H exists, while the tendencies for the other parameters are similar. Finally, we show the results of two specific cases. Fig. 11 shows the parameter points where the 1-step, 2-step, strong 2-step, and 3 or more step PTs occur in the m A vs. m H plane in the Type-I 2HDM with m A = m H ± . The other input parameters are fixed as tan β = 2, cos(β − α) = −0.2 (left panel) or 0 (right panel), and m 3 = 0. We can see the regions for the 1-step and multi-step PTs are almost divided. In addition, R multi in the left panel of Fig. 11 is larger than that in the right panel, which implies the multi-step PTs favor the negative values of cos(β − α). We can also find that the strong 2-step PTs occur only with the mass hierarchy m A > m H . The above features are also seen in the Type-I 2HDM with m H = m H ± . Note that the right panel of Fig. 11 has similar tendencies with the result in the Type-X 2HDM with m A = m H ± as described before. Taking into account the constraint from B → µµ decays, the region of m A (= m H ± ) 340 GeV is excluded. In the survival parameter space in the left panel of Fig. 11, the multi-step PTs occur mostly for m H > ∼ 300 GeV.

A. Higgs trilinear couplings
In this section, to research the possibility of the verification of the multi-step PT by collider experiments, we discuss the Higgs trilinear coupling λ hhh . The coupling λ hhh is derived by calculating the third derivative of the effective potential with respect to the SM-like Higgs fields at the EW vacuum as The trilinear coupling corrected by the leading 1-loop contribution of the top quarks in the SM is written by where N c is the color number of the top quarks. We determine the deviation of the Higgs trilinear coupling from that in the SM as When δλ hhh is equal to zero, the coupling has the same value as in the SM. In the following, we analyze δλ hhh by the Type-I and -X 2HDMs. The current limits on the Higgs trilinear coupling from Higgs pair production are −4.2 < δλ hhh < 10.9 (at 95% CL) from ATLAS [109]. At the future measurement, like the HL-LHC, the limit could reach an accuracy of about 50-60% with 3 ab −1 data [75], while the ILC operating at 500 GeV has the possibility to measure δλ hhh with 27% of precision [76]. 300 GeV. In other words, compared with the results of the 1-step PTs, the values of δλ hhh for the multi-step PTs have a tendency to be large at the same value of m A . The upper right panel shows the smaller tan β is, the larger the maximum value of δλ hhh is for tan β 8. In the lower left panel, when the multi-step PTs occur, the maximum value of δλ hhh becomes larger as cos(β −α) gets smaller except for cos(β −α) 0.25, especially in the negative values of cos(β − α). Note that there are parameter points that have the negative deviations. They are in | cos(β − α)| 0.1 and m 3 50 GeV found out by our analysis. The region where δλ hhh 2.5 for the multi-step PTs is 490 GeV m A 560 GeV, tan β 2, and cos(β − α) −0.2 in our parameter space. In the region, we find that m H is in 320 GeV m H 410 GeV.
The lower right panel of Fig. 12 represents δλ hhh for the 2-step and strong 2-step PTs as a function of m A . It indicates δλ hhh have the possibility of being large as 0.5-2.5 when the strong 2-step PTs happen. When δλ hhh is 2.5 with the strong 2-step PTs, we have found m A 510 GeV, m H 320 GeV, tan β 2, and cos(β − α) −0.25. The deviations δλ hhh 0.5-2.5 would be tested at future colliders such as the HL-LHC and the ILC. The deviations δλ hhh for the Type-I 2HDM with m H = m H ± are shown in Fig. 13. In the upper left panel, the parameter points where the multi-step PTs occur are located on the upper side in m A 400 GeV. The behavior of the predictions for the multi-step PTs in Fig. 13 are similar to the ones in the Type-I 2HDM with m A = m H ± . From panels except for the lower right panel in Fig. 13, we find that the region where the multi-step PTs occur with δλ hhh 2.5 is 590 GeV m A 640 GeV, tan β 2, and cos(β − α) −0.2. In the region, the range of m H is 370 GeV m H 420 GeV. Additionally, in the lower right panel of Fig. 13, the range of δλ hhh for the strong 2-step PTs occur is about 0.5-2.0, while the largest value of such δλ hhh is slightly smaller than that in the Type-I 2HDM with m A = m H ± . When δλ hhh is 2 with the strong 2-step PTs, we have found 600 GeV m A 610 GeV, 300 GeV m H 360 GeV, tan β 2, and cos(β − α) −0.15. Although the constraint from B → µ + µ − excludes the part of the above region (especially the region with m H 330-340 GeV as described in Section V A 2), it still remains and such δλ hhh would be tested at the future collider experiments. at m A 400 GeV, tan β 2, and m 3 20 GeV, and 310 GeV m H 340 GeV. Although the largest value is smaller than that in the Type-I 2HDMs, it can be accessed at the future collider experiments. The lower left panel shows that the values of δλ hhh for the multi-step PTs converge to around 0.7 as m 3 gets larger. Such dependence of δλ hhh on m 3 is not seen in the Type-I 2HDMs. Moreover, the predicted values of δλ hhh stay positive in all regions, although the negative δλ hhh are also predicted in the Type-I 2HDMs. These differences between the Type-I and Type-X are mainly due to the range of cos(β − α). The lower right panel in Fig. 14 shows that the range where the strong 2-step PTs occur is 0.5 δλ hhh 1.2. When δλ hhh is 1.2 with the strong 2-step PTs, we have found m A 400 GeV, m H 340 GeV, and tan β 2. The deviations δλ hhh 0.5-1.2 would be explored at the future colliders like the HL-LHC and the ILC.
The deviations δλ hhh in the Type-X 2HDM with m H = m H ± are shown in Fig. 15. We obtain similar features to the ones in the Type-X 2HDM with m A = m H ± , except that the region where the multi-step PTs occur gets broad upward. From Fig. 15 except for the lower right panel (also from the top left panel in Fig. 10), we see that δλ hhh 1.2 for the multi-step PTs is predicted for 480 GeV m A 510 GeV, tan β 2, m 3  be tested at the future collider experiments. In the lower right panel in Fig. 15, we see that the strong 2-step PTs give δλ hhh 1.2 around m A 490 GeV. However, all regions where the strong 2-step PTs occur is excluded by the constraint from B → µ + µ − as mentioned in Section V B 2. Fig. 16 shows the predictions for δλ hhh in the two cases in Fig. 11. In the left (right) panel of Fig. 16, we take tan β = 2, cos(β − α) = −0.2 (0), and m 3 = 0 in the Type-I 2HDM with m A = m H ± . Compared with the same value of m A , δλ hhh for the multi-step PTs have a tendency to be larger than that for the 1-step PTs. We can also see the largest value of δλ hhh where the multi-step PTs occur at cos(β − α) = −0.2 (left) is greater than that at cos(β − α) = 0 (right).
Meanwhile, δλ hhh for the strong 2-step PTs are relatively large as about δλ hhh 1-2 (left) and 1 (right), respectively. The regions of the strong 2-step PTs do not receive the constraint from B → µµ decays since m A (= m H ± ) > 340 GeV.

B. Gravitational waves from multi-step PT
The first order PT at the EW scale is the source of GW whose typical spectrum has a peak frequency. Therefore if the first order PT occurs multiple times in a multi-step EWPT, the multipeaked GW can be observed in the space-based interferometers. In this subsection, we study such a possibility in the case of the 2-step PT.
The GW spectrum is characterized by two parameters α GW andβ GW at the nuclear temperature T n . Here T n is determined by the condition that one bubble nucleates per Hubble radius S 3 /T n 140 where S 3 is the O(3) symmetric action. The α GW is given by which is the ratio of the latent heat (T n ) to the radiation density ρ rad (T n ) = g * (π 2 T 4 n )/30 where g * is 110.75 in the 2HDMs. The latent heat in the first order PT is calculated as where ∆V is the difference between the effective potential of two phases before and after the PT.
On the other hand,β GW is defined asβ GW ≡ β GW /H n , where H n is the Hubble parameter at T n and β GW is the inverse time duration of the PT There are three contributions to the GW spectrum at a first order PT: Here h is the dimensionless Hubble parameter, f is the frequency of the GW at present, Ω ϕ is the scalar field contribution from collisions of bubble walls [110][111][112][113][114][115], Ω sw is the contribution from sound waves surrounding the bubble walls [116][117][118][119] and Ω turb is the contribution from magnetohydrodynamic (MHD) turbulence in plasma [120][121][122][123][124][125]. Each contribution is given by α GW andβ GW with the velocity of bubble wall v w and the κ ϕ , κ sw , and κ turb which are the fraction of vacuum energy, respectively, converted into gradient energy of scalar field, bulk motion of the fluid, and MHD turbulence. Numerical simulations and analytic estimates of the individual contributions lead to the following formula: • Scalar field contribution Ω ϕ [115] : where the peak frequency is Hz. • Sound-wave contribution Ω sw [119] : where the peak frequency is Hz.
We assume the bubble wall velocity as v w = 1 for simplicity and set [113,126] and κ turb ≈ 0.1κ sw [119]. We compute the GW spectrums from a strong 2-step PT where both the first and second step PTs are first order. The following parameter set in the Type-I 2HDM is chosen as a benchmark point : In the left and right panels in Fig. 17, the paths of the first and second PT in the strong 2-step PT are respectively shown by the black line in the φ 2 vs. φ 1 plane. The contour plots of the effective potential at T n are also given in Fig. 17. The path of the first step PT runs almost along the φ 2 axis from the origin to (φ 1 , φ 2 ) (3 GeV, 115 GeV) at T n 49.9 GeV. The path of the second step PT goes from (φ 1 , φ 2 ) (5 GeV, 126 GeV) to (φ 1 , φ 2 ) (93 GeV, 219 GeV), which is in the direction of the EW vacuum, at T n = 42.9 GeV. The strengths of the first and second step PT are respectively ξ = 2.1 and 4.2, then both of them satisfy the criterion ξ ≥ 1. The values of (α GW ,β GW ) are (8.1×10 −2 , 8.5×10 3 ) for the first step and (0.16, 1.9×10 2 ) for the second step. The GW spectrums h 2 Ω GW from these PTs are shown in Fig. 18. The observable areas by the future space-based interferometers such as LISA [30][31][32], DECIGO [127,128], BBO [77], U-DECIGO [78], Taiji [129,130], and TianQin [131,132] are also presented. The navy and blue lines represent the GW spectrums from the first and second step PTs which have the peak frequencies around 0.1 Hz and 2 × 10 −3 Hz, respectively. The superposed GW spectrum is shown by the red line. We can see that it has a double peak, which can be observed by BBO or U-DECIGO 11 . Additionally, the deviation of the Higgs trilinear coupling δλ hhh is 2.2 for the benchmark point. Such δλ hhh has the possibility to be measured at the HL-LHC and the ILC. Therefore, the signature of the strong 2-step PT at the benchmark point may be observed in the experiments of both GW and colliders. With a combination of these signatures, it might be possible to identify whether the strong 2-step PT occurred in the early universe.

VII. CONCLUSIONS
In this paper, we have studied the parameter regions where the multi-step and strong 2-step EWPTs occur by scanning the parameter spaces in the CP-conserving Type-I and Type-X 2HDMs with m A or m H = m H ± . In the analyses, we have focused on the small m 3 as 0 ≤ m 3 ≤ 100 GeV. is negative with large magnitude, (ii) tan β is small, (iii) cos(β − α) is negative and small, (iv) m 3 is small. The features (ii), (iii), and (iv) are preferred for the negative m 2 2 with large magnitude, which can yield the minimum point along the φ 2 axis. By contrast, the strong 2-step PTs occur only when the mass hierarchy m A > m H exists in our parameter search, while they have similar features as (ii), (iii), and (iv). On the other hand, the VEVs after the first step of the multi-step PTs have a tendency to be located along the φ 2 (or φ 1 ) axis, and the VEVs after the last step PTs likely to lie in the direction of the EW vacuum.
As the possible physical signatures for the multi-step PTs in the collider experiments, we have investigated the deviation of the Higgs trilinear coupling from that in the SM δλ hhh . The maximum value of δλ hhh increases as tan β and cos(β − α) (which is zero in the Type-X 2HDMs) becomes smaller respectively in the case where the multi-step PTs occur. Compared with the results of the 1-step PTs at the same value of m A , the values of δλ hhh for the multi-step PTs have a tendency to be large. In particular, when the strong 2-step PTs happen, δλ hhh are larger than about 0.5 and the largest value of δλ hhh in the Type-I 2HDMs reach over 2. Such deviations would be measured at future colliders like the HL-LHC and the ILC. As the signatures observed by the space-based interferometers, we have computed the GW spectrums from the strong 2-step PT where the first order PT occurs twice. The superposed GW spectrum has the possibility to have a double peak and be observed by the future observers as BBO and U-DECIGO. The multi-step EWPT might be confirmed by combining the information obtained from the future collider and GW experiments.    can cure the contributions from the negative squared-masses.

Appendix B: Tables for number analyses
We show in Tabs. III, and IV, the values or ranges of input parameters where the ratios R multi (36), and R st2 (37) have the maximum values, respectively, for Type-I and -X 2HDMs. Note that we omit the results for the Type-X 2HDMs in Tab. IV because the number of points for the strong 2-step PTs in these cases are not large enough to consider the dependencies of the ratios on the input parameters. We see that e.g., in Tab. III the multi-step PTs favor m A < m H in all four cases. Similar tendencies in Tabs. III, and IV are seen even if we consider the constraint from B → µ + µ − .
Although we do not discuss the strong 1-step PTs in Section V, we also show the results of the number analyses for them because they are still important in the context of baryogenesis. Fig. 20 gives the number of points for the 1-step (yellow) and strong 1-step (orange) PTs, as a function of m A (top), m H (middle) and m A − m H (bottom) in the Type-I 2HDM with m A = m H ± . The red dashed lines represent R st1 , which are the ratios of the number of points for the strong 1-step PTs to that for the 1-step PTs, We see in the middle left panel that R st1 peaks at m A − m H 220 GeV and reaches close to 1. Hence, the strong 1-step PTs favor the mass hierarchy m A > m H (This is consistent with the results in Ref. [14]). This is the same feature as the one which the strong 2-step PTs have, shown in the top right panel of Fig. 3.
We show in Tab. V, the values or ranges of input parameters where the ratios R st1 (B1) have the maximum values, respectively, for Type-I and -X 2HDMs. Similar tendencies in Tabs. III, IV, and V are seen even if we consider the constraint from B → µ + µ − .