Strong first order EWPT & strong gravitational waves in Z3-symmetric singlet scalar extension

The nature of electroweak (EW) phase transition (PT) is of great importance. It may give a clue to the origin of baryon asymmetry if EWPT is strong first order. Although it is a cross over within the standard model (SM), a great many extensions of the SM are capable of altering the nature. Thus, gravitational wave (GW), which is supposed to be relics of strong first order PT, is a good complementary probe to new physics beyond SM (BSM). We in this paper elaborate the patterns of strong first order EWPT in the next to simplest extension to the SM Higgs sector, by introducing a Z3-symmetric singlet scalar. We find that, in the Z3-symmetric limit, the tree level barrier could lead to strong first order EWPT either via three or two-step PT. Moreover, they could produce two sources of GW, despite of the undetectability from the first-step strong first order PT for the near future GW experiments. But the other source with significant supercooling which then gives rise to α∼O0.1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \alpha \sim \mathcal{O}(0.1) $$\end{document} almost can be wholly covered by future space-based GW interferometers such as eLISA, DECIGO and BBO.


I. INTRODUCTION
The early universe must have experienced electroweak phase transition (EWPT), although we do not know its nature, namely whether it is first or second order phase transition.Within the standard model (SM) of particle physics, it must proceed continuously owing to the heaviness of the SM(-like) Higgs boson.However, the SM is commonly believed to be just a low energy effective theory of more fundamental theory, and in the more complete theory EWPT may be first order.Actually, we do have a strong motivation that this may be the case in the context of EW baryogenesis (EWBG) [1,2], which relies on the departure of thermal equilibrium furnished by the first order EWPT.We should further require EWPT to be strong enough so that the baryon asymmetry would not be washed out.This would amount to imposing the following condition: (1.1) with T * being the temperature of EWPT and h * the vacuum expected value (VEV) of the SM Higgs field h at T * .Therefore, SM extensions that could realize strong first order EWPT (SFOEWPT) is of great interest.These extensions could help to build a barrier between the EW vacuum and a metastable vacuum at tree or loop level [2,3].In order to generate a thermal cubic term for h, the latter mechanism usually needs many bosonic degrees of freedom that couple to the Higgs doublet with strength O (1), which thus tends to violate perturbativity near the weak scale.In addition, it suffers from gauge dependence issue [4].On the other hand, the former mechanism utilizing a tree level barrier can avoid these drawbacks.This mechanism is most easily implemented in the extended Higgs sectors by a (supersymmetric) singlet S, containing effective tree-level cubic terms ∼ S 3 + S|H| 2 with H the SM Higgs doublet [5][6][7][8][9][10][11][12][13][14][16][17][18][19][20][21][22][23][24].
If the extended Higgs sector respects some symmetry such as Z 2 , under which S → −S and H → H, an alternative way to the desired tree level barrier is available in the symmetric limit where S does not acquire VEV at the present universe [12,[25][26][27][28][29][30].Such a scenario is associated with multi-step PT's [14,15,31,32,34] which may happen if there is (was) metastable vacua (denoted as Ω meta which breaks some symmetry like Z 2 in this example) except for the desired one (Ω EW ) showing EWSB: The universe may have been once in the intermediate phase Ω meta and then tunneled through a tree level barrier to the phase Ω EW , recovering the Z 2 symmetry. 1Since the barrier is present at tree level, SFOEWPT could be realized simply by reducing the vacuum energy gap between Ω meta and Ω EW , which would lower the critical temperature [14,26].However, very recently Ref. [33] raised the question on easy SFOEWPT in that model.Previous studies did not consider if there is a bounce solution giving S 3 (T ) which satisfies S 3 (T * )/T * ∼ O(140), the condition for bubble nucleation.As noticed in Ref. [33], this condition could rule out a large portion of the parameter space, in particular the most attractive one with weak couplings that leads to the so called nightmare scenario at colliders.As a matter of fact, SFOEWPT can be probed via the Higgs self-coupling measurement [10,22,23,[35][36][37][38][39][40][41][42][43][44][45], however, a sufficiently good precision is unlikely until the next generation of colliders such as the International Linear Collider (ILC) [46] and the Compact Linear Collider (CLIC) [47].The plan of ILC shows that the Higgs self-coupling can be determined with 10% accuracy by upgrading the center-of-mass energy to √ s = 1 TeV [48][49][50]. 1 It is tempting to regard S as the dark matter (DM) candidate.However, it fails, at least being the dominant DM component because the relic density is suppressed owing to the large singlet-Higgs coupling which is required by SFOEWPT [25,26].
Maybe the only available test to the nightmare scenario is the gravitational wave (GW) signal.It is well known that GWs were emitted after the bubble collision during first order PT (FOPT) [51][52][53].The resulting GW spectrum shows a characteristic peak which is related to the PT temperature T * , or the bubble nucleation temperature.Therefore, EWPT, which would have taken place at T * ∼ 100 GeV with the corresponding peak frequency ∼ 10 −5 Hz, may leave imprints at the GW observation experiments [54,55].On the other hand, encouraged by the discovery [56] of GWs by LIGO [57], the next generation detectors: eLISA [58], DECIGO [59] and BBO [60], designed to be sensitive to GW density Ω GW h 2 10 −16 − 10 −10 (depending on frequency ≃ 10 −3 − 10 −1 Hz), will be launched in the near future [54].Thus, SFOEWPT can be examined by the future GW interferometers.Once GW signals with that frequency are observed, it is about to yield deep implications to new physics.
In this paper we study an obvious alternative to Z 2 , the Z 3 -symmetric singlet scalar extension.Theoretically, Z 3 is not inferior to Z 2 at all, both of them frequently used in new physics model building [61][62][63][64][65][66][67][68][69].Although the next simplest extension to the SM Higgs sector, it is strange that such a simple model has not been studied yet in the context of EWPT and GW.Besides, as one of our original but abortive idea, the Z 3 extension is supposed to have one more GW source since the first-step PT is supposed to be first order due to the appearance of cubic term S 3 , and therefore it is a good case in point to demonstrate novel twin-peak GW.Unfortunately, this source is beyond the sensitivity of coming experiments.Still, after employing analytical and numerical methods, we find that the model shows remarkable new features such as the three-step PT, realizing SFOEWPT with strong GW signals which can be well received by eLISA and DECIGO.
The paper is organized as follows: In Section II we introduce the model; In Section III we analyze the vacuum structure at zero temperature as well as its evolution.In Section IV we study the gravitation waves from multi-step PT's in our Z 3 extension of the SM.Then conclusions and discussions are put in the final section.

II. Z 3 SYMMETRIC HIGGS SECTOR WITH A SINGLET SCALAR
Despite of receiving much more attention, a Z 2 discrete symmetry does not take obvious advantage over Z 3 from theoretical viewpoints.On the contrary, we will see that the Higgs sector extended by a Z 3 -symmetric scalar has more interesting features.Besides, this model is one of the very few model which can be studied analytically, at least in the leading order.
On top of the Higgs doublet, the Higgs sector contains an isospin complex singlet scalar S transforming as S → e i2w S with w = π/3 under Z 3 , while the SM fields including the SM Higgs doublet H are neutral under Z 3 .Then the most general renormalizable and Z 3 -symmetric scalar potential V (H, S) is given by Compared to the Z 2 -symmetric model, there is just one more parameter describing the cubic term A s S 3 , and it will give rise to distinguishable difference from the Z 2 -symmetric model.In this paper we do not consider the possibility that S makes the DM candidate [63], because we failed in finding viable parameter space with λ sh ∼ O(0.01) that is necessary to accommodate correct DM phenomenology.But we would like to point out that this Higgs sector could be a part of more complete model where DM is provided by another ingredients of the complete model.For example, one can consider two loop radiative neutrino mass models with DM running in the loop [61,65,70], where DM candidate is well furnished, for instance, a Dirac fermion or even another singlet scalar lighter than S. Since this DM phenomenology involves other parameters and can be irrelevant to PT, we will not enter its details in this paper, keeping in mind that the Z 3 -symmetric extension could be a good simplified model for many BSM for DM and neutrino physics.
Expanding the scalar fields around their classical backgrounds, e.g. S = (s cl +h s +ia s )/ √ 2 and so on, one gets the tree level Higgs potential after dropping the subscripts "cl": 2) The vacuum stability condition reads as At zero temperature T = 0, the SM Higgs parameters are fixed to be λ h ≈ m 2 h /(2v 2 ) = 0.13 up to radiative corrections, with v = 246 GeV and µ 2 h ≈ (88.4 GeV) 2 for m h = 125 GeV.For later use, we give the explicit forms for the field dependent mass squared matrix for the CP-even/odd and three Goldstone bosons: (2.4) where for completeness we also incorporate the leading order thermal masses.In addition, the field dependent masses for the weak gauge bosons and top quarks are given by, for example, Ref. [43].These equations shall furnish the starting point of the loop corrections employed in the program CosmoTransitions [71], which we adopt to study the phase transitions in the model described in this section.

III. MULTI-STEP PHASE TRANSITIONS AND SFOEWPT
Before the start of vacua structure analysis, it is useful to brief the concept of vacuum tunneling at finite temperature, which also plays key roles in gravitational wave radiations.Tunneling from a metastable vacuum to the ground state through a barrier follows the picture of bubble expansion.For a given scalar potential V ( φ, T ) with φ denoting a vector of real scalar fields in the multi dimensional fields space, the (critical) bubble can be found by extremizing the Euclidean action S E (T ), which can be done numerically by the program CosmoTransitions [71].Then, the bubble nucleation rate per unit volume per unit time will be given by Γ(t) = Γ 0 (t) exp[−S E (t)] with the prefactor Γ 0 ∼ T 4 [73].One may also write S E (T ) = S 3 (T )/T where S 3 (T ) is defined as  In order for the nucleated vacuum bubbles to percolate through the whole Universe, the nucleation rate per Hubble volume per Hubble time should reach the unity which determines the FOPT temperature T * .This condition is converted to the very strict relation S 3 (T * )/T * = 4 ln(T * /H * ) ≃ 140 − 150.Such a condition may block some FOPT, especially for the case that the barrier is induced by tree level effects.
A. Vacuum structure: minima in the high temperature expansion

Preliminaries
At T = 0, we are interested in the case where the EWSB but Z 3 -preserving vacuum Ω h ≡ ( h = v, 0) is the ground state, which may be accompanied by a metastable vacuum Ω s ≡ (0, s = 0) or Ω sh ≡ ( h = 0, s = 0).The presence of Ω sh is a new aspect in the Z 3 -symmetric model compared to the Z 2 -symmetric model, and it will make possible three-step PT's in our model.In any case, tree level barrier is indeed furnished, as shown in Fig. 1 To study the vacuum structure and its evolution as T varies, in principle we should find extremes based on the effective potential V eff (s, h; T ) at a given T , which at one loop level contains several terms such as the Coleman-Weinberg (CW) potential and finite temperature corrections [2].But that should render analytical analysis impossible.Therefore, for the purpose of analytic discussion in this Section, we instead work in the high temperature expansion of V eff (s, h; T ), without considering the CW potential 2 .Under this approximation, the tree level effective potential Eq. (2.2) turns out to be the most general form, while the finite temperature corrections are simply encoded in the evolutions of the mass squared parameters s,h the parameters defined at T = 03 and the coefficients c s , c h are calculated to be It is seen that µ 2 s (T ) decreases with T as long as c s = 2(λ sh + 2λ s ) > 0. Note that the linear term A s sT 2 is absent, due to the cancellation between the CP-odd and CP-even contributions, as a result of Z 3 symmetry.

Analysis along the singlet direction
It is illustrative to start the analysis along the singlet direction without taking into account h, where the parameter η ≡ A 2 s /(4λ s µ 2 s ) plays an important role.Under high temperature expansion it evolves as It is straightforward to find out that in different regions of η, the vacuum structures respectively take the forms of: η > 0: There are two minima, separated by the maxima, namely the origin; they are located in the positions For A s > 0, v − is the global minimum, otherwise v + .And we denote the minima as v s .Irrelevant to the sign of A s , the deeper one has negative definite vacuum energy which is a monotonously decreasing function of |η| with the lower limit 8.This case does not admit a very small λ sh since the positive singlet scalar mass requires a large λ sh like the Z 2 case.(Left panel) a smaller λ sh = 0.05 with λ s = 1; (Right panel) a larger λ sh = 0.27 with λ s = 1.In each colored region, we label (using the same color, a convention applied to other figures) all the local minima (with number no lager than two) and their relative orders, says in the green regions there are two minima Ω sh and Ω h with the latter the ground state.
−1 < η < 0: The two minima vanishes and the origin is the only minimum (or only the Higgs doublet acquiring VEV in the two dimensional field space).Note that at very high T where µ 2 s (T ) ∝ T 2 ≫ µ 2 s and thus the potential is always in this region.
η < −1: On top of the origin, there is another minima located at v s,+(−) for A s < (>)0 and the maxima between the two minima is v s,−(+) .Its vacuum energy is This function lies above zero for − 9 8 < η < −1 and crosses zero , the vacuum energy of the symmetric vacuum, at η = − 9 8 and stays below zero for even smaller η.When the two minimums are degenerate, the hight of the barrier separating them is moreover, the distance of two vacua is ∝ |A s |/λ s .Therefore, a larger λ s and a smaller |A s | will help to reduce the barrier and facilitate the bubble nucleation during the transition from the origin to the Z 3 breaking minima.

Analysis in the (s, h) space
Now let us move to the two-dimensional field space, where new minima, such as the trivial one along the doublet direction Ω h , arise.We start from the vanishing tadpole equations: ).In the limit λ sh → 0, they are nothing but the rescaling of the corresponding parameters in the singlet sector by a factor 4λ h .Then it is straightforward to see that Eq. (3.9) admits solutions having nonvanishing s and h if 1/η S > −1 and moreover λ sh s 2 /2 < µ 2 h .These solutions are denoted as (u s,± , v h ) with the singlet VEV similar to Eq. (3.5), with the sign of v h irrelevant and assigned a positive sign hereafter.The extreme having larger magnitude of singlet VEV has the potential to be the minima Ω sh while the other one is a saddle point.To realize the potential, one further condition should be fulfilled, i.e., the Higgs doublet becomes tachyonic in Ω s , otherwise Ω s is the minima rather than Ω sh .We will add more details about this point elsewhere.More explicitly, the conditions for accommodating Ω sh are summarized as The first inequality sign takes " < " and " > " for a negative and positive µ 2 s , respectively.Obviously, a smaller λ sh can readily satisfy all these inequalities and thus Ω sh is well expected.On the other hand, if λ sh is not very small (says a few 0.1) and moreover λ s is relatively large (∼ 1), the appearance of Ω sh is rare.One can clearly confirm these on the η − m s plane, from the gray and green shaded regions in Fig. 2 and Fig. 3.
The potential energy of Ω sh can be written as the sum of two parts, with one the doublet contribution, while the potential energy of the saddle point has a similar expression with f sh,+ (η S ) → f sh,− (η S ) and Like f − (η) introduced before, f sh,+ (η S ) is a monotonically decreasing function of η S and crosses zero at η S = −9/8.The height of the barrier is ∆E sh = A 4 S /(24λ 3 S λ h )(1 + 1/η S ) 3/2 .Note that the second term in Eq. (3.12), namely E h = −µ 4 h /4λ h , is exactly the vacuum energy of the minima along the doublet direct, Ω h .Therefore, Ω sh stays above Ω h if The first can be readily satisfied for a relatively larger λ sh but smaller λ s , and it does not impose much severer condition other than Eq.(3.11).However, for a smaller λ sh and thus λ S > 0, may be just a narrow strip left; see the green shaded regions in the left panels of Fig. 2 and Fig. 3.

More on Ω s and Ω sh
It is of importance to address the connection between Ω sh and Ω s .Although not proved here 4 , they do not coexist with each other.Actually, for most parameter profiles, the first and the second conditions in Eq. (3.11) suffice the presence of Ω sh .However, there indeed exists parameter space where those two conditions are satisfied while λ sh v 2 s /2 > µ 2 h is also met.Then, in this case Ω s has the priority and Ω sh is not present.In the sense of evolution under temperature, P(T ) ≡ λ sh v 2 s (T )/2 − µ 2 h (T ) somehow can be regarded as an order parameter between the two phases Ω s and Ω sh : When P(T ) crosses zero (from above 0), Ω s transits to Ω sh provided that the other two conditions in Eq. (3.11) are satisfied.In other words, the PT Ω s → Ω sh is second order because it is related to continuous variation of P(T ).However, there is another situation in the region with λ sh v 2 s /2 − µ 2 h < 0. That is, maybe those two conditions are violated and then Ω sh is still not present, only with Ω h left.Its implication is that in this case Ω s may transit to Ω h via second order PT instead of first order.Despite of irrelevant to our study, we did observe such kind of PT in the numerical study.
At T = 0 the existence of Ω h and a metastable Ω s /Ω sh does not guarantee SFOEWPT by that tree level barrier.We have to trace the phase evolution up to high temperature.Anyway, these general analysis still guides us to explore the patterns of PT.
But we keep in mind that higher order terms may change them, in particular the behaviors at the higher temperature.
B. Vacuum structure: evolution/phase transition5 Our universe today is assumed to be in the Ω h phase, but it may have experienced other phases at the earlier stage.Due to the multi minima structure in this model, the patterns of PT are rich and complicated.For our purpose, we wish that there was a PT pattern such as two-step PT Ω 0 → Ω s (Ω sh ) → Ω h or even three-step Ω 0 → Ω s → Ω sh → Ω h .Among them, Ω 0 to Ω sh being the first-step cannot happen in the parameter space of our interest, and the reason will be clear in the footnote 6.We will find that other cases are feasible in the weak coupling region.But the way to achieve them is quite non-trivial and we may get a clue from tracking the evolutions of η(T ) and µ 2 h (T ).At very high T ∞ , the universe is in the symmetric phase Ω 0 because η(T ∞ ) → 0 − and µ 2 h (T ∞ ) < 0. When the universe cooled down, local minima Ω h or/and Ω s appeared when µ 2 h (T h ) and η(T s ) respectively approached 0 and −1.Then their birth temperature are estimated to be respectively.The above estimation on T h is mildy lower than the one calculated in the code, which gives T h ≈ 160 GeV for a weak coupling λ sh < 1.If T h > T s , in the weak coupling region of λ sh 1, Ω 0 would immediately roll down to Ω h at T h .Such a case should be avoided and therefore we should at least impose the condition T s > T h .In practice, this condition should be strengthened to allow the commencement of PT Ω 0 → Ω s , that is T h < T * s with T * s the critical temperature where Ω s became degenerate with Ω s , which can be determined by η(T * s ) = −9/8; see the discussion below Eq. (3.7).However, even given a small enough η, whether Ω 0 → Ω s was completed at the PT temperature T * s ∈ (T h , T * s ) depends on if there was a sufficiently large bubble nucleation rate such that S 3 (T * s )/T * s 140.If not, the universe would be confined in the symmetric phase until the transition to Ω h .Unfortunately, usually T * s is significantly below T * s because Ω 0 → Ω s is hampered by a tree level barrier.Then the strengthened condition T * s > T h is transformed into a bound on the initial η (assumed to be negative), where η * ≡ η(T * s ) < −9/8 but the concrete value cannot be calculated in this paper. 6his raises another possibility which might happen in the η > 0 case.It involves another key temperature T 0 s (again assumed to lie above T h ) at which η(T 0 s ) ≃ 0 and the origin turned into a maxima; below T 0 s the model switched to the η > 0 branch.If Ω 0 was hold above T 0 s , then it would transit to Ω s via second order PT instead of first order like before.We found that some of the two-step PT samples Ω 0 → Ω s → Ω h belongs to this pattern; see the example points labelled as triangle in Fig. 8.
Even Ω 0 → Ω s succeeded, it is likely that subsequently Ω s second order transited to Ω sh rather than to Ω h directly as explained before.For instance, the vacuum energy of Ω h , which is E h = −µ 4 h (T )/4λ h , stayed above that of Ω s , which is E s (T ), or again S 3 (T )/T was too large during the window [T sh , T * s ] with T sh the temperature of appearance of Ω sh ; it is also the cross over temperature of transition Ω s → Ω sh .Of importance, to avoid the former, one confronts with an upper bound on η(T ) derived from f − (η(T )) > −24λ 3 s µ 4 h (T )/λ h A 4 s ≡ −ǫ T .Despite the lack of an explicit expression, practically, on account of the suppression µ h (T ) 4

FIG. 6. (Left panel):
A demonstration on the viable three-step EWPT strips on the (η, m s ) plane with λ sh = 0.24, varying λ s = 0.7, 0.8, 0.9, 1.0; for easy identification we plot the region with points labelled by numbers 7, 8..., respectively.The overall picture is consistent with the analytical analysis in the high-temperature expansion.(Right panel): Distributions of two FOPT temperatures, T * s and T * h ; the red dashed line denotes T h , the typical second order PT temperature for Ω 0 → Ω h .thus a small ǫ T at the higher temperature, it gives Combined with Eq. (3.16) we immediately see that the allowed region for η is fairly narrow.In particularly, in the heavy singlet region with A s , µ s ≫ µ h , basically the two-step PT scenario cannot be accommodated in the η < −1 region.And we did not find successful two-step pattern Ω s → Ω h in this region, at least for weak couplings.
If Ω s → Ω h fails, we are left with the three-step PT Ω 0 → Ω s → Ω sh → Ω h .But from the numerical searches it is found that the final-step, in particular in the η < 0 case, tends to suffer a serious bubble nucleation problem and consequently the universe would be confined in the metastable vacuum Ω sh .

C. Numerical samples
We use the code cosmoTransitions [71] for numerical studies on PT in the Z 3 symmetric scalar Higgs sector.Even though only four parameters (λ s , λ sh , A s , µ 2 s ) are introduced, it is still a very time-consuming job to employ a full parameter space scanning.Instead, we shall here focus on the possible PT scenarios indicated by the analytical analysis and demonstrate the typical behaviors of the viable parameter space.
First we consider the η < 0 case.Based on the analysis on vacuum structure, we know that it may allow three-step PT for a relatively large λ sh .But η is supposed to lie within a narrow region to keep Ω sh be the metalstable local minimum instead of global minimum.Whereas the condition to admit the three-step PT is even more strict.For instance, in Fig. 6, we display the feasible regions on the (η, m s ) plane for a given λ sh = 0.24,7 choosing  several values of λ s ; each one corresponds to a slim strip, with a shape similar to the green band shown in the left panel of Fig. 2. It is clearly seen that increasing λ s pushes |η| towards the smaller and narrower region; at the same time, m s is in the heavier region.Eventually, for λ s ≃ 1.1, the allowed region of η is closed.Although we cannot figure out the fundamental reason for that kind of behavior, which must be a complex interplay of several effects, the direct reason is clear from the right panel of Fig. 6.Increasing λ s lowers T * s and it will eventually go below T h , thus shutting down the three-step PT.On the other hand, when λ s becomes fairly small (thus for a much larger v s ), the regions for Ω sh in the (η, m s ) plane will be occupied by that for Ω s .The reason can be traced back to the failure of the third condition in Eq. (3.11).In addition to that, the barrier height between Ω sh and Ω h accordingly increases and thus the tunneling becomes more and more difficult.This compels T * h to approach a fairly low temperature, an obvious trend in the left panel of Fig. 2.These two effects together yield the lower bound on λ s , merely a little bit smaller than 0.7 in this numerical example.
Unlike the above case where only three-step PT is found, in the η > 0 region we find that the two-step PT Ω 0 → Ω s → Ω h can happen, with the first-step either second or first order, depending on the relevant parameters.We describe them one by one in the following: Second order -first oder: This is not surprising since our model basically reduces to the Z 2 -symmetric model in the η → 0 + limit. 8First we study this limit, where the first-step is second order except for a very large λ sh .We obtain similar conclusions in Refs.[26,28,30,33]: The Higgs portal coupling λ sh cannot be made quite small and λ s should also be relatively large; a large λ sh is able to lower T * h thus increasing the SOPT strength, and increasing λ s could further help, shown in Fig. 7.We then investigate the role of A s = 09 , which tends to make Ω s be the global minima; see Fig. 3.For the λ s = 1 example, A s is restricted to be smaller than tens of GeV and thus the resulting deviations as expected are not significant.But it can still increase or decrease T * h with appreciate amount, see the green and blue points in the left panel of Fig. 7.
First order -first oder: For a large λ s = 3, the metastable Ω s can be accommodated for much larger A s ∼ O(100) GeV.That large A s , by contrast, is able to change the nature of transition Ω 0 → Ω s , into the first order type; furthermore, the strength of the second-step can be significantly enhanced and then reopens the smaller λ sh region with λ sh ∼ (a few)×0.1;see the right panel of Fig. 7. From this figure one can find that again the requirement T * s T h yields an upper bound on |A s | 300 GeV in this example.Note that the figures indicate that for a given A s , the allowed region for λ sh is restricted and within this region increasing λ sh could again lead to a lower T * h .Three-step: Three-step PT may be also possible from the analytical analysis and the numerical results confirm this.In fact, we find that it is much easier than the realization in the previous scenario with η < 0, where we have shown fine-tuning the parameters seems unavoidable.Whereas here it is realized in a wide parameter space, including the quite small λ sh region.Since the main features are similar to the previous one, we will not add specific figures for this case.But one may gain some impression from the summarizing Fig. 8.
To have a comprehensive impression on the patterns of PT in the Z 3 -symmetric model, we summarize the parameter region of multi-step PT in Fig. 8, where three-step PT and two-step PT are plotted.The three-step PT case for η < 0 is shown in Fig. 8 (left) with (λ s , m s [GeV]) = (1.0,150) which is corresponding to λ s = 0.9 in Fig. 6 with λ sh = 0.24.On the other hand, we can see the region of the three-step PT case for η > 0 in Fig. 8 (right).Most of the region of SFOEWPT is realized by the two-step PT.For lager A s , the firststep PT significantly becomes the FOPT.Notice that the one-step EWPT becomes the second-order if we assume small λ sh .SFOEWPT with the one-step EWPT along the Higgs doublet direction does not appear for the parameter range in Fig. 8, because it is realized for m s 400 GeV with large λ sh by the non-decoupling thermal loop effects even for A s = 0 GeV as discussed in Refs.[26,30,33,42,43].

IV. STRONG GRAVITATIONAL WAVE FROM PT WITH SUPERCOOLING
FOPT due to a tree level barrier typically gives rise to significant supercooling, which results in sizable free energy release during PT [54].Hence, strong GW is well expected in our Z 3 symmetric model.As a characteristic signal in our model giving three-step PT with two FOPT, two sources of GW are furnished.However, it turns out that only the second FOPT, namely the SOEWPT is detectable in the near future.Actually, this PT pattern is totally within the sensitivities of upcoming GW experiments such as eLIAS, DECIGO and BBO.
A. α & β description on GW from SOPT The GW spectrum from SOPT is very complicated, depending on quite a few details of the bubble dynamics, which is another very complicated job depending on the details of SOPT.But the complications can be parameterized by several parameters, with the most crucial two, α and β, which capture the main features of SOPT dynamics and largely determine the features of GW spectrum.We will follow the conventions in Ref. [54].
The parameter α is the total energy budget of SOPT normalized by the radiative energy with g * (= 108.75) being the relativistic degrees of freedom in the plasma at the PT temperature T * .The liberated latent heat ǫ = −(∆V + T ∂V /∂T )| T * , with ∆V the vacuum energy gap between two vacua.For SOPT which typically has a significant supercooling, the latent heat actually is the vacuum energy.For most SFOEWPT parameter region, we have α ≪ 1 owing to the smallness of ∆V .Another parameter β is defined to be the variation of action with respect to time at T * : with the Hubble constant H * ≡ 1.66 √ g * T 2 * /m pl .So, its inverse characterizes the duration of PT (τ ∼ 1/β) thus the GW peak frequency; usually it is much shorter than the Hubble time scale 1/H * .It is convenient to introduce the dimensionless parameter β ≡ β/H * ≫ 1.
The stochastic GW background in the linear approximation receive three contributions (In principle, in the case with two SOPT there are two GW sources thus six contributions.): where three terms stand for relics originating from bubble collision, sound waves and magnetohydrodynamics (MHD) turbulence in the plasma, respectively.Their structures can be factorized as with a denoting the subscripts in Eq. (4.3).Concretely, • Bubble collision: In the envelope approximation, c col = 1.67 [75] with v b being the velocity of the bubble wall.κ col is the fraction of latent heat deposited in a thin shell close to the PT front.The shape factor is defined as where f col is the red-shifted peak frequency, with h * = 1.65 × 10 −2 mHz T * 100 GeV g * 100 the value of the inverse Hubble time at T * redshifted to today.An analytical treatment to this source was made in Ref. [72].
• Sound waves & MHD turbulence: Both are due to bulk motion, giving c = 2.65 × 10 , respectively [76,77].The enhancement H * /β is traced back to the longer lasting time in producing GWs than the previous case.κ sw ≈ α/0.73 (for v b ≃ 1) is the fraction of latent heat transferred to the bulk motion of the fluid, while κ turb ≈ ǫ turb κ sw with ǫ turb ∼ 5 − 10% the fraction of bulk motion that is turbulent [77].The shape factors for these two cases are with the corresponding redshifted peak frequencies given by f sw[turb] = 1.14 [1.64] These values obtained from simulations suffer from uncertainties and can only be fully trusted in certain regions of (α, v b ) [54]; for instance, the expression Ω sw h 2 is safely reliable only for α 0.1 [77], but we may still use it for the larger α.

B. Excellent prospects of Z 3 -symmetric model at GW detectors
As one of the main original motivation to study the Z 3 -symmetric singlet scalar extension, we expected that there would be two sources of GW and therefore a distinguishable twin-peak GW could show up.Qualitatively it is true and we indeed find that two SOPTs can come from both three and two-step PTs, thus contributing to two GW sources.Unfortunately, quantitatively it is disappointing.
We display the results on the (α, β) plane in the Fig. 9 and Fig. 10, with the experimental sensitivities of eLISA [54,78] and DECIGO [59] labelled by the shaded regions.The sensitivity regions of four eLISA detector configurations described in Table I in Ref. [54] are denoted by "C1", "C2", "C3" and "C4".The expected sensitivities for the future DE-CIGO stages are labeled by "Correlation", "1 cluster" and "Pre" following Ref.[59].The transition temperature T * depends on the model parameters (see, Fig. 6 and Fig. 7) and the velocity of the bubble wall v b is uncertain.Although the experimental sensitivities on the (α, β) depend on T * and v b , we take T * = 50 GeV and v b = 0.95 as a reference for the purpose of illustration.It is seen that typically one needs α O(0.01) for the near future detection.However, the first source from Ω 0 → Ω s turns out to be undetectable since it always gives α 0.01.On the other hand, the other source, in particular in the three-step PT case, is very promising and almost the whole parameter space can be covered.One of the main reasons causing this difference is that the first-step happened at a relatively high temperature T * s 160 GeV, which typically is rather higher than the EWPT temperature T * h 100 GeV; recalling that α ∝ 1/T 4 , thus the first source is suppressed.A lower T * h also Regardless of the wall velocity issue, there is another serious problem which probably renders the three-step PT scenario irrelevant to EWBG.During EWPT, the phase outside the bubble, Ω sh , already breaks the EW symmetry and therefore the baryogensis process, which should be effective outside the bubble, actually was ineffective due to the suppression on sphaleron rate by e −E sph (T * h )/T * h ∼ e − 4π g 2 h * /T * h ≪ 1.However, we cannot completely excludes this case, since the idea of locally recovering EW symmetry used in the supersonic EWBG might work here, and it deserves a specific study elsewhere.

V. CONCLUSIONS AND DISCUSSIONS
Gravitational waves, which are supposed to be relics of strong first order PT, is a good probe to new physics beyond SM complementary to collider searches.Within SM, EWPT is second order, while EWBG requires SFOEWPT.Such a situation inspires a great many extensions to the SM Higgs sector where a potential barrier is created during EWPT by, e.g., non-thermal tree level effects.One simple implement of this idea is in the mixing Higgs-singlet models which includes a sizable cubic term like S|H| 2 [9][10][11][12][16][17][18][19][20][21][22][23], and then a doublet-singlet mixing could contribute to enhance the strength of the FOPT.As a result, such kinds of models can be tested by the synergy between the measurements of various Higgs boson couplings at future collider experiments and the observation of GWs at future spacebased interferometers as discussed in Refs.[22,23].In another implementation imposing unbroken discrete symmetry like Z 2 [12,25,26,[28][29][30]33], multi-step PT could utilize a tree level barrier.But generically the absence of mixing renders the tests at colliders very difficult without taking enough large λ sh coupling as discussed in Refs.[26,30,33,42,43].In this paper, we have focused on such the nightmare scenario.
Following the second line, in this paper we studied in great detail the patterns of SFOEWPT in the next to simplest extension to the SM Higgs sector, the Z 3 -symmetric extension to the Higgs sector.It involves one more term, the cubic term for the singlet (as a comparison, the Z 2 model can be a limit of turning off the corresponding coupling A s ), that gives rise to remarkable difference in PT than the Z 2 -symmetric model.It can not only significantly enhance the potential barrier thus the PT strength, but also lead to three-step PT through the intermediate phase Ω sh in contrast with the Z 2 -symmetric model.SFOEWPT is expected to be fairly strong due to the significant supercooling.Especially, the three-step PT produces two sources of GW.Despite of the undetectability from the first-step in the near future, the other source with significant supercooling and thus leads to α ∼ O(0.1), basically can be completely covered by eLISA and DECIGO.Additionally, there were also studies on FOPT thus GWs associated with the singlet scalar, without taking into account EWPT [83,84].
We end up the paper with some open questions.First, we encounter some new phenomenologies/problems in multi-step SFOEWPT, which needs further understanding.Second, if the three-step PT in our model can give sizable baryon number asymmetry is unclear and should be studied elsewhere, although it would depend on the specific models for baryogenesis involving a new source of CP violation.In the last, the SOPT phenomenologies are quite rich in this model, and if we give up the Z 3 -symmetric ground state, there are even more intriguing scenarios.For instance, the pattern Ω 0 → Ω h → Ω sh is in the bulk parameter and has important consequence to dark matter.We leave this part to future publication.
FIG.1.Examples for two vacua separated by a tree level barrier: Ω h with a metastable Ω sh (left panel) or metastable Ω s (right panel).The contours are equipotential curves, rescaled by large numbers which are of no importance here.The red stars and hearts stand for minima and saddle points, respectively.The dashed lines schematically show the tunneling paths from Ω sh/s → Ω h .

FIG. 7 .
FIG.7.First order phase transition temperature as the function of λ sh for the two-step PT (Ω 0 → Ω s → Ω h ) in the η > 0 region.For the second-first order PT pattern (left panel), we show the Z 2 -like case with A s = 0 GeV (black lines) and the deviations from non-zero A s by fixing λ sh (green and blue lines), for four cases (λ s , m s [GeV]) = (1, 100), (1, 150), (3, 100),(3, 150).The first-first order PT pattern (right panel) arises in the larger λ s region which allows a larger A s ; for (λ s , m s [GeV]) = (3, 150), we show A s [GeV] = 100, 200, 300 (blue lines).For each dashed line, the upper and the lower ends denote T * s and T * h , respectively.In these plots we just keep the points which give SFOEWPT.

FIG. 8 .
FIG.8.Global picture of multi-step PT in the (A s , λ sh ) plane for (λ s , m s [GeV]) = (0.9, 150) (left panel) and (1.0, 100) (right panel).PT of three-step (red, circle), two-step (green, triangle for the second-first order PT or star for the first-first order PT) and one-step (blue, square) are plotted.Filled plots satisfy the condition of SFOEWPT in Eq. (1.1).In η < 0 region, the three-step PT can happen only in a very narrow space, consistent with Fig.6.
1.Examples for two vacua separated by a tree level barrier: Ω h with a metastable Ω sh (left panel) or metastable Ω s (right panel).The contours are equipotential curves, rescaled by large numbers which are of no importance here.The red stars and hearts stand for minima and saddle points, respectively.The dashed lines schematically show the tunneling paths from Ω sh/s → Ω h .