Phase Transitions in Twin Higgs Models

We study twin Higgs models at non-zero temperature and discuss cosmological phase transitions as well as their implications on electroweak baryogenesis and gravitational waves. It is shown that the expectation value of the Higgs field at the critical temperature of the electroweak phase transition is much smaller than the critical temperature, which indicates two important facts: (i) the electroweak phase transition cannot be analyzed perturbatively (ii) the electroweak baryogenesis is hardly realized in the typical realizations of twin Higgs models. We also analyze the phase transition associated with the global symmetry breaking, through which the Standard Model Higgs is identified with one of the pseudo-Nambu-Goldstone bosons in terms of its linear realization, with and without supersymmetry. For this phase transition, we show that, only in the supersymmetric case, there are still some parameter spaces, in which the perturbative approach is validated and the phase transition is the first order. We find that the stochastic gravitational wave background is generated through this first order phase transition, but it is impossible to be detected by DECIGO or BBO in the linear realization and the decoupling limit. The detection of stochastic gravitational wave background with the feature of first order phase transition, therefore, will give strong constraints on twin Higgs models.


Introduction
Naturalness of electroweak symmetry breaking (EWSB) has been a guideline for exploring physics beyond the Standard Model (SM). Popular scenarios of physics beyond the SM include supersymmetry (SUSY) and composite Higgs, which are still promising solutions to the (large) hierarchy problem, since they remove the sensitivity of the weak scale to quadratically divergent quantum effects from physics at high energy scales such as the Planck scale and the grand unification scale. However, the discovery of the SM-like Higgs boson and nothing else at the Large Hadron Collider (LHC) poses a problem for naturalness. No new colored particles predicted in these popular scenarios have been observed so far at the LHC, which already leads to fine-tuning in the Higgs potential at sub-percent level. Although we do not know whether nature takes thought for this little hierarchy problem or not, it is interesting to pursue possibilities to ameliorate this fine-tuning and to explore their implications for particle phenomenology and cosmology.
The twin Higgs mechanism [1] is an attractive idea to provide a solution to the little hierarchy problem without introducing new colored states. There are several variations to realize this idea, but every twin Higgs model starts with the assumption that the SM Higgs field can be considered as one of the pseudo-Nambu-Goldstone bosons (pNGBs) arising from spontaneous breaking of a global symmetry G, such as U (4) symmetry, that contains SU (2) A × SU (2) B symmetry in its subgroups, to a smaller group H, such as U (3). Here SU (2) A and the mirror (or twin) SU (2) B are gauged and interchanged under a (approximate) Z 2 symmetry. The SU (2) A gauge symmetry is identified with the SU (2) W symmetry in the SM and spontaneously broken by the vacuum expectation value (VEV) of the Higgs field. By introducing a SU (3) C mirror color symmetry and mirror fermions that are charged under SU (3) C × SU (2) B , quadratic divergence to the Higgs potential coming from the SM colored particles (and SU (2) W gauge bosons) are canceled by the mirror colored particles (and SU (2) B gauge bosons).
The original realization of the twin Higgs idea, which is now called the Mirror twin Higgs model [1], has a mirror copy of all the SM particle content related to the Z 2 symmetry. On the other hand, the Fraternal (minimal) twin Higgs model [2] has a smaller twin particle content, that is, twin W bosons, twin gluons and twin fermions corresponding to the third generation. Other twin particles are not required since the corresponding SM particles give less important contributions to the Higgs potential. In any case, due to the Z 2 symmetry, the quadratic terms of the Higgs potential accidentally preserve the original global symmetry G and the pNGBs associated with G → H breaking are protected from radiative corrections, allowing the natural EWSB. Since every twin partner is not charged under the SM gauge group, this mechanism realizes the so-called neutral naturalness 1 and enables the model to evade stringent LHC bounds. In this mechanism, there still remains the "large hierarchy problem", that is, the hierarchy between the global G → H breaking scale, expected to be up to 5-10 TeV, and more fundamental scales such as the Planck scale or the grand unification scale. It is expected to be addressed by the UV completion such as SUSY [4][5][6][7][8][9][10] or composite Higgs [11][12][13][14][15].
If the twin Higgs mechanism is really realized in nature, it may have significant impacts on cosmology since the models predict new particles in the mirror (or twin) sector and have a rich structure in the Higgs sector. One immediate concern in the Mirror twin Higgs model is the existence of mirror copies of light SM particles. In fact, a twin photon and twin neutrinos give sizable contributions to the radiation energy density, which is strongly disfavored by measurements of Cosmic Microwave Background (CMB) [16] and Big Bang Nucleosynthesis (BBN) [17]. This issue has been recently studied in [18][19][20][21]. The effects of twin baryons on the large scale structure and CMB are also investigated in Ref. [22]. On the contrary, there is a candidate of dark matter in the Mirror twin Higgs model, which has been investigated in Ref. [23]. The Fraternal twin Higgs model does not lead to an extra dark radiation component but still accommodates a dark matter candidate [24][25][26][27].
In addition to the modification in the relatively late-time cosmology described above, we also naturally expect that the thermal history of earlier Universe can differ from the standard one, such as cosmological phase transitions, in the twin Higgs models. For example, in the SM without any extensions, the electroweak phase transition is known to be crossover [28]. Since the Higgs sector is significantly different, it is non-trivial whether the electroweak symmetry is really restored and how the EWSB proceeds even if any. Kilic and Swaminathan addressed the first question and showed that the electroweak symmetry is really restored at high temperature [29]. In this paper, we try to address the second question, that is, how the EWSB proceeds in the twin Higgs models. If the electroweak phase transition is first order, which means that it proceeds through the bubble nucleation, then it is attractive in the cosmological point of view. A first order phase transition generates stochastic gravitational wave (GW) background from bubble collisions [30][31][32][33][34], sound waves [35][36][37][38], and turbulence of the plasma [39][40][41][42][43][44]. The typical peak frequencies of GWs originating from the first order phase transition associated with the EWSB are O(10 −3 ∼ 1) Hz, which are good targets of gravitational wave detectors such as DE-CIGO [45] and BBO [46]. Therefore, if the twin Higgs models generally predict a first order phase transition associated with the EWSB, they can be tested by these observations. Moreover, if the SM Higgs VEV at the critical temperature T C , φ A (T C ), is larger than the critical temperature, φ A (T C )/T C > 1, inside the bubble, the electroweak phase transition is strong first order and the sphaleron decoupling condition is satisfied. (See App. A.1 for the definitions of a first order phase transition and a strong first order phase transition.) A strong first order electroweak phase transition accommodates electroweak baryogenesis [47,48], so that the present baryon asymmetry of the Universe can be explained depending on the model parameters other than the (minimal standard) Higgs sector.
In this paper, we examine how the EWSB and the global symmetry breaking G → H (typically, U (4) → U (3)) proceed in thermal history of the Universe. In particular, we address the question if these phase transitions can be first order in that framework. We find that in the non-supersymmetric case, thermal potential around both the electroweak and global symmetry breaking cannot be analyzed perturbatively, which suggests that both phase transitions are unlikely to be first order and hence we can expect for neither the electroweak baryogenesis nor the generation of stochastic gravitational wave background. Even in the case with supersymmetric UV completion, by limiting ourselves to the linear realization and the decoupling limit where only the mirror stops are added to the nonsupersymmetric model, we find that the EWSB cannot still be analyzed perturbatively and the conclusion is still robust. For the global symmetry breaking, however, we show that, with an appropriate parameter choice, the symmetry breaking can be first order and the associated gravitational wave background is generated, but unfortunately it is too small to be detected by DECIGO or BBO.
We organize this paper as follows. In Sec. 2, we see the electroweak vacuum structure of the non-supersymmetric twin Higgs models as well as supersymmetric ones and discuss the fine-tuning in the model parameters. In Sec. 3, we study twin Higgs models with and without SUSY at non-zero temperature and examine how the EWSB proceeds. In Sec. 4, we examine how the global symmetry breaking proceeds and show that in supersymmetric twin Higgs models the first order phase transition can be realized for appropriate parameter choices but the resultant gravitational wave background is undetectable at planned gravitational wave detectors. Sec. 5 is devoted to our concluding remarks and comments. In App. A, we exhibit the detailed calculations and expression of the thermal potential and the spectrum of the stochastic wave background from a first order phase transition.

Twin Higgs Models
In the first part of this section, we review the twin Higgs mechanism, which provides a solution to the little hierarchy problem [1,2]. The Higgs mass formulae are also presented. We then discuss the degree of fine-tuning to realize the adequate EWSB in this scenario. In the second part, we describe a supersymmetric realization of the twin Higgs mechanism.

The non-supersymmetric twin Higgs
In the twin Higgs mechanism, the SM Higgs field is identified with pNGBs arising from spontaneous breaking of an approximate U (4) symmetry (explicitly broken by the Yukawa and gauge couplings). 2 Let us consider a linear realization of the mechanism and write a U (4) symmetric potential of a complex scalar field H with the fundamental representation, where λ > 0 is required from the stability of the potential. This potential drives the scalar field H to get a nonzero VEV, f ≡ |H| = m/ √ 2λ. Then, the global U (4) symmetry is broken down to U (3) yielding 7 NGBs. The U (4) symmetry contains the subgroups SU (2) A × SU (2) B and the scalar field can be decomposed as H = (H A , H B ), where H A transforms as a doublet of SU (2) A while H B does as a doublet of SU (2) B . H A is identified with the SM Higgs doublet and the SU (2) A symmetry is regarded as the ordinary SU (2) W gauge symmetry. The SU (2) B symmetry is gauged and becomes the twin SU (2) W . Then, the 6 pNGBs are eaten by the gauge bosons after the symmetry breakings while the remaining one is the observed SM-like Higgs boson h. A physical heavy exotic Higgs h, corresponding to the radial direction, has the mass m h = √ 2λf from Higgs mechanism. As described in the introduction, SU (2) A and the twin SU (2) B are interchanged under a Z 2 symmetry. In the Fraternal twin Higgs model [2], only the Z 2 partners of the third generation of quarks and leptons and the partners of gluons (twin gluons) as well as the twin SU (2) W gauge bosons are introduced. Then, the twin Higgs doublet H B has the following Yukawa coupling similar to the SM top Yukawa coupling, where Q a (a = 1, 2, 3) are twin left-handed top (bottom) quark doublet charged under the twin SU (3) C and t a R are twin right-handed top quarks. y t is the twin top Yukawa coupling whose value is almost the same as the ordinary top Yukawa coupling y t due to the approximate Z 2 symmetry.
The two scalar doublets H A and H B receive quadratically divergent corrections from the top and twin top quarks respectively as well as corrections from the SU (2) A and SU (2) B gauge bosons at one loop. In addition, they receive corrections from the gluons and twin gluons at two-loop level. The quadratically divergent part of their potential is given by where g 2 , g 2 are the SU (2) W and SU (2) W gauge couplings, g 3 , g 3 are the SU (3) C and SU (3) C gauge couplings and Λ is a cutoff scale. The exact Z 2 symmetry leads to y t = y t , g 2 = g 2 , g 3 = g 3 which guarantee that the quadratically divergent part of the potential respects the full U (4) symmetry. Then, the NG nature of the Higgs field is not explicitly broken by the quadratically divergent corrections, addressing the little hierarchy problem. However, the SM Higgs would be exactly massless and inconsistent with our Universe if the U (4) and Z 2 symmetries are exact. Thus we need small breakings of these symmetries. Let us consider the breaking of the U (4) and the Z 2 symmetries to give the appropriate effective Higgs potential. First of all, the gauged SU (2) A × SU (2) B group has already broken the U (4) symmetry explicitly. In addition to the quadratically divergent corrections, this generates logarithmically divergent contributions to the quartic couplings of the form |H A | 4 + |H B | 4 , which do not respect the U (4) symmetry and then contribute to the Higgs boson mass. The explicit Z 2 symmetry breaking is also needed otherwise the hierarchy, v 2 A f 2 , is not fulfilled, which is required to satisfy the constraint from the Higgs coupling measurement. We do not specify a mechanism to generate this breaking in this paper, but just encapsulate the breaking effect in the quadratic and quartic terms of H A . The effective potential of the scalar field H we consider here is then summarized as The first term is the U (4) conserving term coming from the original potential Eq. (2.1) rewritten in terms of H A and H B and the corrections in Eq. (2.3), which determines the U (4) symmetry breaking scale f . The second term that breaks the U (4) symmetry includes the gauge (and top Yukawa) contributions in the Coleman-Weinberg potential. Thus κ 1 will be of order g 4 2 /16π 2 log (Λ/g 2 f ). The third and fourth terms are the Z 2 breaking terms. The third term is induced, e.g., by the quadratic corrections with Z 2 -breaking part in the gauge and matter sector. ρ 1 in the fourth term includes the contribution of the one-loop Coleman-Weinberg potential. In the Fraternal twin Higgs model, the fourth term could arise from the Z 2 breaking effect such as the absence of the U (1) Y gauge symmetry in the twin sector. However, this effect is of order g 4 1 /16π 2 , where g 1 is the U (1) Y gauge coupling constant and tiny. In summary, we take λ, f, κ 1 , σ 1 and ρ 1 to be the model parameters and require σ 1 , κ 1 , ρ 1 < λ so that the second, third and the forth terms in Eq. (2.4) are regarded as perturbations to the first term.
At energies well below the symmetry breaking scale f , we can integrate out the Higgs field H B , which enables us to work with an effective field theory of the SM Higgs field H A . The effective potential of the SM Higgs field can be obtained by setting H B as Using this relation, we find This potential coincides with the SM Higgs potential when the parameters κ 1 , σ 1 are identified with where λ SM ∼ 1 8 is the SM Higgs self-coupling, v A = 246 GeV is the VEV of the Standard Model Higgs field. As denoted above, to satisfy the constraint from the Higgs coupling measurement, the VEV of the Standard Model Higgs field is required to be satisfactorily small compared to the U (4) symmetry breaking scale, that is, v 2 A f 2 . Let us discuss the EWSB conditions (v A 246GeV and m h 125GeV) in the twin Higgs models precisely with the potential (2.4). By expressing the potential (2.4) in terms of the two physical modes φ A and φ B with H A ≡ (0, φ A / √ 2) and H B ≡ (0, φ B / √ 2) and requiring the minimization conditions, ∂V /∂φ A = ∂V /∂φ B = 0, we find the potential minimum given by .
Evaluating the mass matrix ∂V /∂φ i ∂φ j (i, j = A, B) around the potential minimum, we obtain the mass eigenvalues of the system, that is, the SM Higgs boson h and the heavy exotic (global symmetry breaking) Higgs h as [7,8] where the plus sign in front of √ 1 − A corresponds to m 2 h and the negative sign corresponds to m 2 h . With v 2 A /f 2 1, the SM Higgs mass is approximately given by Since we have five parameters f, λ, σ 1 , κ 1 and ρ 1 , after imposing the EWSB conditions v A 246 GeV (2.8) and m h 125GeV (2.9), the system is now described by three parameters. As noted above, we impose the conditions λ > σ 1 , κ 1 , ρ 1 to keep the philosophy of the twin Higgs models. Fig. 1 shows the parameter space that satisfies these conditions for v A /f = 0.223 (f = 1.1 TeV) and 0.123 (f = 2 TeV). We also confirmed that the condition λ > σ 1 is always satisfied. Note that the parameters λ, κ 1 and ρ 1 cannot take arbitrary small values because tiny λ, κ 1 and ρ 1 cannot realize the SM-like Higgs mass. In fact, we can see from Fig. 1 that the smallest values of λ, κ 1 and ρ 1 are roughly given by λ κ 1 ρ 1 0.05. This bound will play important roles when we analyze the dynamics of a phase transition as we will see in Sec. 4. The smallest values of λ, κ 1 , ρ 1 are not sensitive to the breaking scale v A /f and SM-like Higgs mass m h 125GeV. σ 1 > 0 guarantees our assumption of the two-step phase transition as we will see later.
Let us finally examine the fine-tuning in this effective potential. We estimate the degree of tuning by the measure defined in Ref. [49],  Figure 1: The (yellow) region that satisfies λ > ρ 1 , κ 1 and σ 1 > 0 is shown for v A 246 GeV, m h 125GeV, and v A /f = 0.223 (left) and 0.123 (right). The regions above the blue, green and red curves satisfy with the constraints λ > ρ 1 , λ > κ 1 and σ 1 > 0, respectively. effective potential, the set of the observable and parameter that gives the smallest measure is the following one, (2.12) In this calculation, we simply assume the soft breaking scenario, σ 1 ρ 1 , which means that the twin Z 2 symmetry is only broken by the soft term σ 1 f 2 . 3 In order to solve the little hierarchy problem in twin Higgs models, ∆ σ 1 should not take an arbitrary small value. Thus, the symmetry breaking scale f is bounded from above in light of naturalness.

Supersymmetric twin Higgs models
To address the large hierarchy problem in the twin Higgs scenario, SUSY can provide an attractive solution. Parallel to the case of the ordinary Minimal Supersymmetric Standard Model (MSSM), where Higgs chiral multiplets consist of a pair of doublets, supersymmetric twin Higgs models generally contain four Higgs doublets, (2.13) The chiral multiplets H u , H d are fundamental under the U (4) symmetry and the U (4) multiplets are decomposed into the visible sector fields H A u , H A d and the twin sector fields The superpotential contains an extended version of the ordinary µ-term, Including soft SUSY breaking mass terms, the quadratic part of the U (4) symmetric potential in supersymmetric twin Higgs models is given by (2.14) The quartic part of the U (4) symmetric potential is model dependent. The first term of (2.4) contains the quartic term |H A | 2 |H B | 2 . However, SUSY forbids this type of couplings without further modification of the Higgs sector. There are several proposals to obtain Higgs couplings with twin Higgs fields. Refs. [4][5][6][7] have introduced a massive singlet chiral superfield S with a superpotential SH u H d . The effective theory after integrating out this singlet contains the quartic term |H u H d | 2 . Ref. [8] has considered an additional contribution to the D-term potential from a new U (1) gauge symmetry, under which both the Higgs and the twin Higgs fields are charged. In this paper, we do not go into the details of a specific supersymmetric twin Higgs model, instead, we simply assume the existence of an appropriate U (4) symmetric quartic term and try to extract general features of a supersymmetric twin Higgs scenario. We next consider possible sources of the breaking of the U (4) and the Z 2 symmetries in this scenario. In the non-supersymmetric minimal model, the U (4) symmetry breaking arises only from quantum corrections or from some explicit breaking terms. On the other hand, supersymmetric models have the D-term potential, which breaks the U (4) and the Z 2 symmetries. Here we have assumed the minimal realization of the twin Higgs mechanism, where the twin partner of U (1) Y is not introduced. Unfortunately, the size of the Z 2 breaking in the above D-term potential is insufficient to realize the required hierarchy between the electroweak breaking scale, v A , and the U (4) breaking scale, f . Then, we simply assume the following Z 2 breaking soft mass terms, In order to make the discussion independent of the form of quartic couplings, we take the decoupling limit of the SUSY heavy Higgses and match the theory to the nonsupersymmetric twin Higgs potential. In the decoupling limit, the four Higgs doublets can be written as follows in terms of H A and H B in the non-supersymmetric twin Higgs model, (2.17) Thanks to the approximate Z 2 symmetry, they are almost equal, tan β A tan β B . In the rest of the discussion, we simply assume β A = β B = β and g 2 = g 2 . Note that by taking the decoupling limit (2.17), the structure of the Higgs potential is essentially the same as that of the potential given by (2.4) discussed in the previous subsection. When we require that the supersymmetric twin Higgs potential is matched with Eq. (2.4), we obtain the following relations, cos 2 2β + δκ, where δσ, δκ and δρ represent the radiative corrections. With these expressions, we can evaluate the SM-like Higgs mass from (2.9). Note that it is difficult to realize the SM-like Higgs mass only with the quartic couplings in the D-term potential, κ = g 2 2 8 cos 2 2β and ρ = g 2 1 8 cos 2 2β. We simply assume that there is an additional contribution or a radiative correction to κ and ρ to realize the SM-like Higgs mass. As mentioned in the previous subsection, we impose the EWSB conditions and the conditions λ > σ 1 , κ 1 , ρ 1 to consider the general feature of SUSY twin Higgs models. In Sec. 3 and Sec. 4, we will discuss the order of the phase transitions with imposing these conditions.

The electroweak phase transition
As discussed in Sec. 2, twin Higgs models generally accommodate breakings of the two symmetries. One of them is the standard EWSB and another is the breaking of the U (4) symmetry to the U (3) one, through which the SM Higgs field is identified with one of the pNGBs. Let us call the phase transition corresponding to the latter breaking the U (4)-breaking phase transition. In this paper, we analyze not only the electroweak phase transition but also this U (4)-breaking phase transition in cosmology. In this section, we study the order of the electroweak phase transition in twin Higgs models, with and without SUSY, especially.
Before going to the detailed calculation, we would like to discuss the thermal history in the early Universe. At high-temperature, both of φ A and φ B fields are trapped at the origin of the potential due to the thermal mass terms. When the temperature cools down, another minimum different from the origin appears. Below the critical temperature, φ A and φ B fields eventually roll down or tunnel to the true vacuum, and the U (4) symmetry and its subgroup, the SM electroweak symmetry, finally break down. However, we do not know how these two phase transitions proceed. Let us denote the temperatures when φ A and φ B fields acquire their VEVs by T A and T B , respectively. In general, there are three possible trajectories of these two phase transitions, which are schematically described in Fig. 2. The red line (1) shows the trajectory of a two-step phase transition with T B T A , in which φ B field acquires its VEV first and φ A field does later. The blue solid line (2) shows the trajectory of a one-step phase transition with T A ∼ T B , in which the Higgs field rolls (or tunnels) to the true vacuum directly. The green dotted line (3) shows the trajectory of another two-step phase transition with T A T B , in which φ A field acquires its VEV first and φ B field does later. In this paper, we consider the case with T B T A and we call the phase transition at which φ B field acquires its VEV, the U (4)-breaking phase transition. Let us consider the condition under which this case happens. The thermal mass terms for φ A and φ B fields are given by where ζ A and ζ B represent the numerical coefficients depending on the coupling constants. The critical temperatures T A and T B are evaluated by the condition m A (T A ) = m B (T B ) = 0, which yields Taking into account the twin The region with σ 1 > 0 is also shown in Fig. 1. We shall study the strength of the electroweak phase transition. The thermal resummed effective potential for both of H A and H B Higgs fields are calculated at the one-loop order in the same way as Ref. [29]. We take account of the top, twin top quarks, SU (2) W × U (1) Y gauge bosons, and SU (2) W gauge bosons, respectively, because they give dominant contributions to the effective potential. The general expression of the thermal effective potential is summarized in appendix A.1.
Let us calculate the thermal one-loop resummed effective potential Eq. (A.1) for H A and H B in the non-supersymmetric case starting from the effective potential (2.4). Since we take account of not only the H A field but also the H B field, we consider the following background fields, The number of degrees of freedom (d.o.f) and the field dependent masses of SU (2) W ×U (1) Y gauge bosons, SU (2) W gauge bosons, top quark and twin top quark are given by respectively Note that here we considered the Fraternal model where the mirror U (1) gauge fields are absent, but that we expect that the basic results are unchanged even if we include them since the U (1) gauge coupling is tiny. With the supersymmetric completions visible and mirror stops might also contribute, but we do not take account of them by assuming they are sufficiently heavy through the Higgs φ B 's VEV. With this assumption, our conclusion is applicable also to the case with the supersymmetric UV completions. The one-loop effective potential V eff is then given by See the App. A.1 for the details. Here, κ, ρ and σ are the tree-level couplings and do not include the contribution of the one-loop Coleman-Weinberg potential. In addition, we consider the ring diagram contributions denoted by V ring discussed in appendix A.1 to improve the perturbativity. Since the masses of the SU (2) W gauge bosons originating from the VEV of the H B field are much larger than thermal corrections to the masses around the critical temperature, the SU (2) W ring diagram contributions can be neglected. On the other hand, the ring diagram contributions coming from SU (2) W × U (1) Y gauge bosons are not negligible and we need to take them into account. V ring was computed in Ref. [50] and is given by with Here, n i represents the number of d.o.f. for each longitudinal mode. We do not take account of transverse modes of the SU (2) W × U (1) Y gauge bosons because the magnetic masses they receive from the environment are suppressed by g 4 1 or g 4 2 and hence give minor contributions.
In our case, the electroweak phase transition occurs after the U (4)-breaking phase transition. Therefore, during the electroweak phase transition, φ B already gets a non-zero VEV, φ B (T ) = 0. Then, in the same way as Eq. (2.5), we integrate out the φ B (T ) field by setting Here, φ A(B) (T ) represent the temperature dependent VEVs, respectively. It should be noticed that, when we take the T = 0 limit, Eq. (3.19) is reduced to Eq. (2.5). The one-loop resummed effective potential Eq. (A.1) for φ A can be written as where V 0 , V CW , V thermal , and V ring are given by Eq. (3.11), (3.12), (3.13), and (3.14), respectively. For the zero temperature part V 0 +V CW , we set the renormalization conditions given by Neglecting O(φ 6 A ) terms, we obtain the following expression, where the suffix i represents only the SM contribution. Now the system is parameterized only by the U (4)-breaking scale f since the condition (3.19) and renormalization condition Eqs. (3.21) and (3.22) completely fix the other model parameters, κ 1 , σ 1 , ρ 1 , and λ. In Ref. [29], it was shown that the one-loop effective potential obtained by use of the relation Eq. (3.19) exhibits the restoration of the electroweak symmetry at high temperature, which guarantees the presence of the electroweak phase transition. In this paper, we try to clarify the order of the electroweak phase transition. For this purpose, we will first check the validity of perturbative expansion near the critical temperature. As is seen in App. A, the perturbative expansion is valid only when the following condition is satisfied, where T C is the critical temperature of the EWSB. Here the critical temperature T C is defined so that the electroweak symmetry preserving and breaking vacua are degenerate. φ A (T C ) represents the expectation value of φ A for a breaking phase at T C . In Fig. 3, the ratio φ A (T C )/T C is plotted for each U (4) symmetry breaking scale f . We have evaluated the one-loop resummed effective potential given in Eq. (3.20) without resort to the high temperature expansions. It is easily seen that the larger a breaking scale f is, the smaller φ A (T C )/T C is. This fact can be easily understood as follows. The thermal contributions from the twin particles could strengthen the first order nature of the EWSB. However, the twin partners acquire masses proportional to φ B (T ) through the Higgs mechanism.
Thus, a larger breaking scale f leads to larger masses of the twin particles, which easily induces thermal decoupling of twin particles during the electroweak phase transition. This decoupling makes φ A (T C )/T C in our case approach the value in the standard model case. Hence a larger v A /f indeed increases φ A (T C )/T C . However, the largest value of φ A (T C )/T C for f > 2v A required by the constraint of the Higgs coupling measurement is at most 0.2, which is not large enough to satisfy the criteria (3.24). Therefore, we conclude that the higher order effects cannot be neglected and the perturbative expansion is not valid near T C . For the correct analysis, lattice simulations are required. This result has an important implication for the electroweak baryogenesis because it requires the sphaleron decoupling condition φ A (T C )/T C > 1 around the critical temperature. Our results strongly suggest that this condition is hardly satisfied in the Fraternal twin Higgs model as long as the condition (3.19) is valid, and we cannot expect for the implementation of the electroweak baryogenesis. This conclusion remains unchanged even if we go beyond the Fraternal model, as long as the condition (3.19) and the assumption of the trajectory of two-step phase transition are adopted. We do not exclude the possibility to have the strong first order electroweak phase transition once we relax one of these assumptions, which is beyond the scope of the present study.
Finally let us comment on some issues on UV completions. We here do not assume concrete UV physics (SUSY and composite Higgs) in our analysis and analyze the electroweak phase transition by use of effective field theory for the φ A field. As long as its usage is valid, our result is still robust in supersymmetric and composite twin Higgs models. However, it was shown in Refs. [51,52] that the electroweak phase transition can be the strong first order in the composite Higgs scenario. In the setup adopted in Refs. [51,52], the electroweak phase transition and the confinement phase transition, which corresponds to the U (4)-breaking phase transition in twin Higgs models, occurred simultaneously. In addition, the SM-like Higgs field couples with an additional scalar field. Thus this approach does not apply to our consideration.

The U (4)-breaking phase transition
In this section, we explore the U (4)-breaking phase transition in twin Higgs models with and without supersymmetric completion. For the concrete calculation, we adopt the Fraternal model, but general features would apply to other models. 4 As discussed in Sec. 3, we assume that the U (4)-breaking phase transition occurs first and the electroweak phase transition does next in the following discussion.

The case of the twin Higgs model without UV completion
Let us first consider the twin Higgs model without any UV completions, in the sense that no new particles other than the mirror particles to the SM are involved. The U (4)-breaking phase transition generally depends on UV physics such as SUSY and composite Higgs. However, if new particles in the UV completion are sufficiently heavy during the phase transition, we can safely neglect the effect of these particles. We shall study the strength of the U (4)-breaking phase transition by using the potential (2.4) with this assumption.
In our set up, the Higgs field H A is trapped at the origin of the potential H A = 0 due to the thermal mass term during the U (4)-breaking phase transition. Thus, we take the background fields as and calculate the resummed one-loop potential given by Eq. (A.1) for the field H B . We take account of the twin top and SU (2) W gauge bosons which give dominant contributions to the effective potential. On the other hand, a larger quartic coupling λ > g 2 2 makes the U (4)-breaking phase transition weaker φ B (T C )/T C < g 2 (see Eq. (A.18)) hence we consider a small quartic coupling λ < g 2 2 . Since the quartic coupling λ < g 2 2 is smaller than g 2 and y t , we neglect the H A and H B loop contributions to the effective potential in the following discussion.
With the field dependent masses of the twin top and SU (2) W gauge bosons given by Eqs. (3.7) and (3.9), The one-loop effective potential V eff is expressed as In order to find the ring diagram contribution V ring , we need to evaluate the thermal masses of the SU (2) W gauge bosons denoted by (A.10). We here take account of the one-loop selfenergy of longitudinal modes [54] in which only one generation (third generation) is included for the Fraternal model. The transverse modes receive the magnetic masses, but they are suppressed by the factor of g 4 2 and hence we omit them as is the case of the electroweak phase transition. We then obtain the ring diagram contribution V ring given by When we use the high-temperature expansions given by Eq. (A.6) and Eq. (A.7), the resummed one-loop effective potential takes the following form: where T 2 , (4.10) Thanks to the twin Z 2 symmetry, y t y t and g 2 g 2 , the U (4)-breaking phase transition described by the potential (4.9) is similar to the electroweak phase transition in the SM which has been analyzed by perturbative [55] and non-perturbative [28,56] methods. The most reliable approach to clarify the order of the phase transition is the lattice simulation [28,56]. It was shown that the electroweak phase transition in the SM is the first order when m H 70 − 80 GeV (or λ SM 0.04) is satisfied. One might wonder if the difference between U (4)-breaking sector and the SM sector prevents us from adopting the results of the SM to the U (4)-breaking case. But these differences are negligible for our purpose in the following reasons. First of all, the breaking scale of the U (4)-breaking phase transition, f , is different from that of the electroweak phase transition, v A . However, the order of the electroweak phase transition in the SM depends on the parameter λ SM /g 2 2 [56], but not v A . Thus, we have only to identify λ + κ 1 in our model with λ SM in the SM electroweak phase transition. Second, there is no U (1) Y gauge boson in the Fraternal twin Higgs model. Since the U (1) Y gauge coupling g 1 is tiny compared to the SU (2) W gauge coupling g 2 , we can neglect this effect safely. Indeed, the original paper [56] also does not include U (1) Y and they concluded that the error due to this assumption is small enough. Finally, the coefficient of the thermal mass (4.6) for the Fraternal model 5 differs from that of the SU (2) W gauge bosons in the SM. In the lattice simulation [56], they use the three-dimensional effective Lagrangian obtained by integrating out all fermions and the longitudinal modes of SU (2) W gauge bosons [57][58][59], which affects the values of the parameters λ SM and g 2 . We confirmed that this difference gives only 10 percent changes in the parameters of three-dimensional effective Lagrangian, and hence we can safely neglect it. Therefore, the order of the U (4)-breaking phase transition can be analyzed by use of the result of the electroweak phase transition in the SM. We conclude that the U (4)-breaking phase transition is the first order when λ + κ 1 0.04 is satisfied, thanks to y t y t and g 2 g 2 .
As discussed in Sec. 2.1, the parameters κ 1 and λ are bounded below, λ + κ 1 0.1, due to the EWSB conditions and the conditions λ > σ 1 , κ 1 , ρ 1 as we can see from Fig. 1. Therefore, it cannot satisfy the condition for the first order U (4)-breaking phase transition, λ + κ 1 0.04. We also expect no gravitational wave production because of the absence of a first order phase transition in the case of twin Higgs models without any UV completions. The differences in the Fraternal and Mirror models give minor effects and are within the uncertainties in our estimate. More generally, our conclusion is robust in any models as long as the tree-level potential is given by Eq. (2.4) and there are no additional light degrees of freedom so that the thermal masses for twin SU (2) W gauge bosons do not differ so much.

The case of supersymmetric twin Higgs models
In the previous subsection, we do not consider effects of UV physics such as composite Higgs and SUSY on the U (4)-breaking phase transition. If other fields strongly couple to the Higgs field H B , we cannot apply the argument in the previous subsection. We here consider supersymmetric twin Higgs models and explore the order of the U (4)-breaking phase transition. Especially, since any such models contain twin stops, which are strongly coupled to the Higgs field H B and possibly light at the restored phase in the absence of the Higgs VEV, we focus on the effect of light twin stops. Hereafter, we take the decoupling limit, simply assuming that every supersymmetric partner except for twin stops acquires a large soft mass and decouples with thermal plasma during the U (4)-breaking phase transition. We will show that there is some parameter space where the U (4)-breaking phase transition is the first order and estimate the gravitational wave amplitude generated through this phase transition.
Let us calculate the one-loop resummed effective potential (A.1). Since we take the decoupling limit as explained in section 2.2, the background fields are given by We take account of the left and right-handed twin stops, the twin top quarks and the SU (2) W gauge bosons, which give dominant contributions to the effective potential. We neglect the Higgs loop correction as is the non-supersymmetric case. The tree level potential V 0 in Eq. (A.1) is concretely written as where the second term includes the D-term contribution, κ ⊃ g 2 2 cos 2 2β/8. The thermal one-loop corrections V CW and V thermal are evaluated as follows. The field dependent masses of twin top quarks and SU (2) W gauge bosons are given by Eqs. (3.9) and (3.7). Those of the left and right-handed twin stops can be written as In order to calculate the ring diagram contribution V ring , we need to evaluate thermal masses of the longitudinal mode of the SU (2) W gauge bosons. As a result of the twin Z 2 symmetry, thermal masses of the SU (2) W gauge bosons can be calculated in the same way as the case of the MSSM [54]. The thermal masses of the longitudinal mode of the SU (2) W gauge bosons, the left and right-handed twin stops are given by Then the temperature dependent mass matrix is given by Moreover, in our set up, the twin QCD two-loop contribution is non-negligible compared to the resummed one-loop effective potential because the strong coupling g 3 and the top Yukawa coupling y t are large compared to the other matter couplings. In the MSSM, the sunset diagram, which gives the dominant contribution, is evaluated in Ref. [60]. We adopt it and calculate the two-loop twin QCD contribution as 3T .

(4.26)
In this expression, the high-temperature expansion [61] and mass-averaging approximation [62] are used. According to the discussions in Refs. [63,64], their usage is justified for our purpose. It should be noticed that this negative logarithmic dependence of φ B in Eq. (4.26) gives an additional contribution to the potential barrier between the origin and another minimum. Without taking this contribution into account, we would underestimate φ B (T C )/T C .  φ B (T n )/T n α β/H(T n ) 682 1 7 × 10 −3 7 × 10 4 Table 1: Parameters T n , φ B (T n )/T n , α and β/H(T n ) for the evaluation of the spectrum of gravitational wave background with the benchmark point λ = 0.05, g 3 = 1, κ 1 = 0.05, M stop = 0, X t = 0, v A /f = 0.123 and tan β = 10.
As discussed in App. A.1, in order to have the first order phase transition and gravitational wave production, φ B (T C )/T C g 2 is required. This ratio gets larger for a smaller λ + κ 1 (see App. A.1.). As discussed in Sec. 2 (see Fig. 1), we have the conditions λ > 0.05 and κ 1 > 0.05, from the requirements λ > ρ 1 , κ 1 and m h 125 GeV. Thus hereafter we take λ 0.05 and κ 1 0.05 as the benchmark point. For simplicity we require the quartic coupling κ 1 is dominated by the D-term, κ ( g 2 2 /8) cos 2 2β, so that tan β 10. The value of the twin QCD coupling constant g 3 can be somewhat different from the value of the visible QCD coupling constant g 3 because the exact Z 2 symmetry is not necessary from the view point of naturalness [2,18]. Here we simply set the twin QCD coupling to be g 3 = 1. The change of the value of g 3 allowed by naturalness leads to a 10% effect for φ B (T C )/T C . In addition, we take X t = 0 in our evaluation for the following reason. A non-zero X t tends to induce unwanted color-breaking vacua. In order to avoid the appearance of such vacua, larger soft masses are required, which reduces the ratio between the effective mass and the cubic term. Thus, with a non-zero X t , φ B (T C )/T C will get smaller compared to the case with a vanishing X t . Now the ratio φ B (T C )/T C is determined by the twin stop soft parameters and the U (4)-breaking scale f . Figure 4 shows the ratio φ B (T C )/T C as the function of the left and right-handed twin stop (common) soft masses | m 2 The renormalization scale is set to be µ = T . We can see that the ratio φ B (T C )/T C takes the maximal value for the massless limit of the light twin stop M stop 0, which is roughly 0.9. For other choices of the ratio v A /f 0.1, required from the point of view of naturalness, we confirmed that φ B (T C )/T C 0.9 > g 2 for M stop 0. Thus, for this parameter choice, the phase transition is the first order, which leads to the generation of the gravitational wave. Note that here we admit the strong violation of the Z 2 symmetry in the soft stop mass, but we assume that the Z 2 symmetry is hold for tan β otherwise we cannot have the Mexican-hat type U (4)-breaking potential. Now let us evaluate the spectrum of gravitational wave background generated in this model. For this purpose, we need to estimate the nucleation temperature T n , the latent heat density α and the duration of the phase transition β (see App. A.2 for the detailed definition). They can be obtained by solving the bounce equations for the thermal resummed effective potential V (φ B , T ) = V 0 + V CW + V thermal + V ring + V (2) thermal . Table 1 shows the values of these parameters for our benchmark points, λ = 0.05, κ 1 = 0.05, g 3 = 1, tan β = 10, X t = 0, and v A /f = 0.123. Figure 5 shows the spectrum of the gravitational wave background for our benchmark point (see App. A.2 for the formalism to calculate it). The most dominant source of the gravitational wave for our benchmark point is found to be the sound wave of the plasma bulk motion after the bubble collision, Ω gw h 2 Ω sw h 2 given by (A.33). The peak frequency is around O(10)Hz and the peak amplitude of gravitational wave is around O(10 −19 ) due to the large β/H(T n ) 7 × 10 4 and small α 7 × 10 −3 . We can easily see that it is well below the sensitivities of DECIGO and BBO. It was pointed out in Ref. [65] that the formula (A.33) overestimates the gravitational wave amplitude in the region of large β/H (typically, β/H 10 2 for α > 5 × 10 −3 ). Thus it should be emphasized that the gravitational wave amplitude in Fig. 5 is an upper bound and we expect that it will be much weaker in reality.
It is nontrivial whether our benchmark point, which gives the maximal ratio φ B (T C )/T C , gives the maximal amplitude of the gravitational wave background. We numerically confirmed that it is approximately maximal for our benchmark point. Concretely, • For λ, κ 1 and M stop , we confirmed that smaller λ + κ 1 and M stop give larger gravitational wave amplitude. Since we restrict them as λ > 0.05, κ 1 > 0.05, and M stop > 0, 6 our benchmark point gives the maximal amplitude.
• The peak amplitude of gravitational wave, Ω peak gw h 2 , does not depend on the breaking scale f . We can write the effective potential as where φ, f and M stop are parameters normalized by the temperature, φ ≡ φ/T, f ≡  6 Note that when we allow the negative twin stop soft masses, the gravitational wave amplitude would be larger. In this case, however, we have to take account of the SU (3) C breaking minimum hence we do not consider such a scenario in this paper. tional wave, Ω peak gw h 2 , only depends on the α and β/H parameters at T = T n hence we get Ω peak gw h 2 = Ω peak gw h 2 ( φ, f , M stop )| T =Tn . The nucleation temperature is roughly given by T n T B T C , where T B and T C are given in Sec. 3 and App. A.1, respectively. From the expression (3.2), we can easily find f /T n f /T B = const. In addition, from the expression (A.18), the fraction φ B (T n )/T n φ B (T C )/T C does not depend on the breaking scale f (the quartic coupling ξ is less sensitive to the change of T n ). Thus, when we vary the breaking scale f with M stop = 0, the peak amplitude of gravitational wave does not change. On the other hand, the peak frequency ν peak is proportional to the nucleation temperature, T n , hence a smaller f leads to a lower peak frequency due to the lower nucleation temperature. We numerically confirm this behavior.
• We numerically confirm that a smaller tan β makes the gravitational wave amplitude larger. However, a smaller tan β leads to a larger up-type Higgs-top Yukawa coupling, Y t = y t / sin β. Here we impose the perturbative condition of the Yukawa coupling Y 2 t /(4π) 1 at the electroweak scale. This condition gives tan β 0.28. For tan β = 0.28, the peak amplitude of gravitational wave is larger than that of tan β = 10 by merely around factor 10.
Note also that the change of the value of the twin QCD coupling constant allowed by naturalness affects the amplitude of gravitational wave by around factor 10 at most, and hence this effect does not change our result significantly. Thus, we conclude that, even if we take the effect of a light twin stop into account, it is almost impossible to generate gravitational wave background detectable by DECIGO or BBO.
Finally, we would like to give some comments. We have assumed that φ B acquires the VEV first and φ A does later. In order to verify this assumption, we have calculated the thermal resummed effective potential V (φ A , φ B , T n ) for both of the Higgs fields φ A and φ B when φ B gets the VEV at T = T n . We numerically confirmed that the potential minimum appears only in the φ B direction at T n given in TABLE. 1. Therefore, the assumption of two-step phase transition is validated.
The resummed effective potential at finite temperature depends on a gauge-fixing parameter. In our calculation, we adopted the Landau gauge. The effect of gauge dependence is discussed in, e.g., Ref. [66,67]. According to Ref. [67], the uncertainty due to gauge choice is roughly one or two order magnitude for Ω gw h 2 . Even when we take this uncertainty into account, the gravitational wave amplitude shown in Fig. 5 still does not reach the detectable regions by DECIGO and BBO. Therefore, our conclusion is still robust.

Discussion and conclusion
We have investigated the dynamics of the electroweak phase transition and the phase transition associated with global U (4) breaking in twin Higgs models with and without supersymmetric completion. In Sec. 3, we found that the electroweak phase transition in twin Higgs models cannot be analyzed perturbatively as long as the effective potential is given by (2.4) and (3.19). It does not satisfy the condition of a strong first order phase transition, and hence we cannot expect for the realization of the electroweak baryogenesis as well as the generation of gravitational wave background. In Sec. 4.1, we considered the U (4)-breaking phase transition in twin Higgs models without any UV completions such as composite Higgs and SUSY. We confirmed that the U (4)-breaking phase transition is the first order only when λ + κ 1 0.04 is satisfied. However, as discussed in Sec. 2.2, we obtained the relation λ + κ 1 > 0.1 in order to realize the adequate EWSB and the conditions λ > σ 1 , κ 1 , ρ 1 . Thus, the U (4)-breaking phase transition cannot be the first order, and we expect that there is no gravitational wave production. In Sec. 4.2, we considered the U (4)-breaking phase transition with supersymmetric UV completions in the decoupling limit where only the effect of light twin stops is taken into account. We calculated the resummed effective potential including the dominant two-loop twin QCD contribution. Then, we confirmed that the U (4)-breaking phase transition can be analyzed perturbatively only when the light twin stop masses with M soft 0 are realized. We calculated the largest possible gravitational wave amplitude within the parameters for which the EWSB conditions and the conditions λ > σ 1 , κ 1 , ρ 1 are satisfied. However, we found that the gravitational wave amplitude cannot reach the detectable regions by DECIGO and BBO.
We conclude that it is impossible to produce large enough amplitude of gravitational wave to be detected by DECIGO or BBO in twin Higgs models, under our assumptions such as taking the decoupling limit, the perturbative conditions λ > σ 1 , κ 1 , ρ 1 and the trajectory of two-step phase transition. We need to give a comment. If there is an additional field strongly coupled to the Higgs fields H A and H B , the dynamics of the electroweak phase transition and the U (4)-breaking phase transition will be changed due to the additional contribution to the effective potential. For example, as mentioned in Sec. 2.2, there is a singlet scalar field coupled to the Higgs field H A and H B in F-term twin Higgs models. If such a singlet scalar field is sufficiently light during the U (4)-breaking phase transition, the situation might be dramatically changed. In this paper, we do not consider such specific cases because we are mostly interested in giving model independent predictions.

A Finite temperature effective potential and phase transition
In this appendix, we give the details of the calculations used for evaluating thermal potential and stochastic gravitational wave background from the first order phase transition in Secs. 3 and 4. We also give the criteria to judge the validity of the perturbative calculation for thermal potential.
A.1 Thermal effective potential and validity of perturbation theory Here we give the way to calculate the thermal potential used in Secs. 3 and 4. We mainly follow the discussions in Ref. [68].
The thermal effective potential is divided into three parts and can be schematically written as Here V 0 , V CW and ∆V th represent the tree level potential, the one-loop Coleman-Weinberg potential and the thermal contributions respectively. Since we shall see that the perturbative calculation is not necessarily justified at high-temperature, we will later take into account higher order effects partially to improve the perturbativity. The Coleman-Weinberg potential V CW is given by where n i is the number of degrees of freedom of a particle i, m i (φ) represents the field dependent mass of the particle i, and µ is a renormalization scale, and (−) F i gives 1 for bosons and −1 for fermions. Here we adopt the DR renormalization scheme. Thermal contributions to the effective potential include the one-loop effective potential given by where i runs the particle species and the suffixes B and F represent Boson and Fermion contributions, respectively. We here adopt the imaginary time formalism. For later use, we note that at high temperature m 2 (φ)/T 2 1, they are approximated as where γ E is the Euler's constant. However, it will be immediately seen that this perturbative expansion breaks down at high temperature. The quadratic divergent contributions to the self-energy from the n-loop diagram, often called as the ring diagram or daisy diagram [62], behave as [69] where a is a constant determined by the coupling constants, which are the expanding parameters in the zero-temperature perturbative calculations. Thus, for aT 2 /m 2 (φ) 1, the perturbative calculation is not valid especially in calculating the critical temperature at the phase transition.
By taking a closer look at the structure of the divergences, we can see that they come from the Matsubara zero mode of bosonic particles. Thus this problem is relaxed by "resumming" the ring diagrams of bosonic particles where we replace the mass of the bosonic particle m Bi (φ) in the one-loop thermal potential by the dressed one, where Π i (T ) is the one-loop self energy of the bosonic particle corresponding to the ring diagrams. Here c denotes the contribution of gauge and Yukawa couplings. This is equivalent to adding to the thermal potential so that After the resummation, the n-loop quadratically divergent diagram behaves as Thus for a < 1, the condition for the perturbative expansions to be validated is improved as aT m(φ) 1. (A.14) When non-Abelian gauge fields are involved, we need to take into account another subtle issue. Although the transverse modes of the gauge fields are massless at the one-loop perturbative calculation, it is known that through the non-perturbative process it receives the so-called magnetic mass, ∼ g 4 T 2 , with g being the gauge coupling. Then, with a similar discussion given above, the higher loop of non-Abelian gauge bosons will give the contributions with the powers of g 2 T /m(φ) [70,71] and the perturbation breaks down at high temperature [72], In this case, even the resummed effective potential (A.1) is not reliable and the dynamics of phase transition should be analyzed by lattice simulations. Since we expect the parameter a given above is at most unity, we conclude that the resummed effective potential is valid for γ 2 ≡ g 2 T /m(φ) gT /φ < 1 when m(φ) gφ.
Let us now give our criteria for a first order phase transition and a "strong" first order phase transition to occur. Starting from the one-loop effective potential (A.1), we can approximate it as Here M 2 (T ), E and ξ(T ) represent a temperature dependent mass, a numerical coefficient depending on coupling constants, and a temperature dependent self-coupling, respectively. Note that the coefficient E comes from the loops from bosonic particles (see Eqs. (A.6) and (A.7)). The thermal potential (A.16) can have two minima. One is at the origin, while the other is not, depending on the temperature and other model parameters. As the temperature decreases, the two minima can get degenerated. We define the temperature at which the two minima degenerate as the critical temperature, T C . At that temperature, the effective potential (A.16) is written as where φ(T C ) = 0 is the other minimum at T C . Or we can write Below the critical temperature, the minimum other than the origin is energetically favored and hence tunneling from the origin (symmetric phase) to the other minimum (broken phase) can occur and the symmetry breaks down. We have seen that at the small field values, including the origin, the resummed effective potential (A.1) or its approximated one (A.16) is not reliable. However, we here give the criteria that the perturbative calculation is allowed to use and the symmetry breaking is first order if the potential minimum in the broken phase satisfies the condition that the perturbative calculation is valid, gT C /φ(T C ) < 1, since the tunneling rate is determined mainly by the information of the potential around the broken phase but not the symmetric phase. Moreover, we define the phase transition as strong first order if is satisfied. The difference between first order phase transition (gT C /φ(T C ) < 1) and strong first order phase transition (T C /φ(T C ) < 1) is important when we consider electroweak baryogenesis because the sphaleron decoupling condition is given by Eq. (A.19).
On the other hand, this difference is not important when we discuss the gravitational wave background generated by first order phase transition. Since a first order phase transition proceeds through bubble nucleation, the production of gravitational wave background requires a first order phase transition, not a strong first order phase transition as we will see later. From (A.18), we can see that a strong first order phase transition takes place if the self-coupling ξ(T C ) is small enough and the cubic prefactor E is large enough. This is because the parameter E determines the height of the barrier between the origin and the other minimum. Since the cubic term comes from the bosonic loop contribution, bosons strongly coupled to φ are needed for a strong first order phase transition. Before closing this subsection, let us comment on the effect of the ring diagram on the strength of the phase transition. In the expression of the ring diagram contribution (A.11), the thermal field dependent mass m(φ, T ) is roughly given by Eq. (A.10). After the resummation, if the thermal mass is much larger than the zero-temperature part, that is, m 2 B (φ) Π at T = T C , (m 2 ) 3 2 behaves like a constant term Π 3 2 (T C ) which does not give the potential barrier. This effect makes φ(T C ) small hence the resummation generally makes the phase transition weaker.

A.2 Phase Transition and Gravitational Waves
In this appendix, we review how a first order phase transition proceeds and the stochastic gravitational wave background generated from it is evaluated.
A first order phase transition occurs as a result of true vacuum bubble nucleations. This is understood as quantum or thermal tunneling from a false vacuum to a true vacuum that is separated by a potential barrier. The tunneling rate or the bubble nucleation rate Γ(T ) per unit volume and unit time at finite temperature is evaluated as [73]  Here r is the radial coordinate in the three dimensional polar coordinate system and φ False is the field value of the false vacuum.
The time or the temperature of the phase transition is characterized by the nucleation time t n or the temperature T n , defined as a temperature when the nucleation probability inside one Hubble volume H 3 (T ) gets unity, Since the dominant contribution in the integral (A.24) comes from that around t ∼ t n or T ∼ T n , it can be approximated as Γ(T n ) H 4 (T n ) = 1, (A. 25) which can be used to determine T n . Since the EWSB takes place at T n ∼ O(100) GeV, the bounce action at the time of bubble nucleation is roughly given by Generally speaking, the nucleation temperature is lower than the critical temperature T n < T C . In order to determine the nucleation temperature as well as the bubble profile accurately, we need to solve the equation of motion (A.22) with the boundary condition (A.23) numerically. We here adopt a method dubbed as the under/over-shooting method, developed in Ref. [74]. Now let us give the expressions of the spectrum of the gravitational background from the first order phase transition. Since the broken phase is energetically favored, the nucleated bubbles expand, and collide each other, and finally the whole Universe settles down to the true vacuum. Since the bubble collisions as well as the plasma bulk motion induced by the bubble dynamics are highly inhomogeneous and violent process, gravitational waves are emitted through such processes.
The spectrum of the gravitational wave is determined by the (initial) kinetic energies of the bubbles and the duration of the phase transition. The former is provided by the latent heat density ∆ρ = ρ(φ False , T n ) − ρ(φ True , T n ), where ρ(φ, T ) is the thermodynamic internal energy, but not the potential energy. By identifying the effective thermal potential with the free energy, F(φ, T ) = V (φ, T ), we obtain We parameterize the kinetic energy of bubbles by a dimensionless parameter α representing the ratio between the latent heat density and the radiation energy density, α = ∆ρ ρ rad , (A. 29) where ρ rad = g * π 2 T 4 * /30 denotes the energy density of radiation. Here T * T n is the temperature at which the gravitational waves are emitted. The duration of the phase transition is characterized by the parameter β, defined by Γ(t) Γ 0 e βt , (A. 30) with Γ 0 being a constant. β is expressed in terms of the bounce action as (A. 31) It has been argued that not only the bubble collision or the scalar field dynamics, but also the plasma dynamics caused by the bubble dynamics source the gravitational waves [30][31][32][33][34]. It is indeed found to be the dominant contribution to the gravitational wave background since due to the interaction between the scalar field bubble wall and the plasma, the energy originally carried by bubble walls is quickly taken away to the plasma bulk motion. According to the popular convention, we further classify it into the sound waves in the plasma described in the linear regime, which are generated by the bubble motion and generate gravitational waves around the bubble collision, and the turbulence of plasma bulk motion further developed in the non-linear regime after the bubble collision. Then the total contribution can be schematically written as Ω gw h 2 = Ω bubble h 2 + Ω sw h 2 + Ω tur h 2 Ω sw h 2 + Ω tur h 2 , (A.32) where Ω bubble , Ω sw and Ω tur denote the contributions from the bubble collisions, sound waves and turbulence of the plasma, respectively. For the contributions from sound waves, we adopt the expressions in Ref. [37,41,75,76] as, so that it gives the maximal estimate for the amplitude of the gravitational wave background. In our setup the latent heat density is small, α = O(10 −3∼−2 ). If the bubble wall velocity is smaller and in the deflagration regime, the gravitational wave background is much smaller. It should be also noted that the formula of the gravitational wave coming from the sound waves, (A.33), does not necessarily work and is likely to overestimate the gravitational wave amplitude for large β/H (typically, β/H > 10 2 for α > 5 × 10 −3 ) [65]. Thus, to be precise, our estimate based on the formula (A.33) should be regarded as the upper bound of the gravitational wave amplitude.
Note that by adopting the formula with the envelope approximation in Ref. [78] the contributions from bubble collisions turned out to be subdominant in our setup. 7 Thus we safely omitted the contributions from the bubble collisions in our plots.