First-order electroweak phase transition in a complex singlet model with $\mathbb{Z}_3$ symmetry

We consider an extension of the Standard Model with a complex singlet scalar, where a global $U(1)$ symmetry is explicitly broken to $\mathbb{Z}_3$ symmetry. We study the two-step electroweak phase transition in the model and find that it can be of first order if the heavy scalar mass falls in the range of $1.2-2$~TeV and the mixing angle $\left | \theta \right |\gtrsim 0.2$ ($11.5^{\circ}$). The Higgs signal strength measurements at the LHC, on the other hand, restrict the mixing angle $\left | \theta \right |\lesssim 0.4$ ($23^{\circ}$). Future colliders including high-luminosity LHC can probe the remaining parameter space of first-order phase transition in this scenario. After the $U(1)$ symmetry breaking, the pseudo-Goldstone boson becomes a dark matter candidate due to a hidden $\mathbb{Z}_2$ symmetry of the model. We find that the pseudo-Goldstone boson can make up a small fraction of the observed dark matter and escape from the constraints of current direct detection. We also show that the stochastic gravitational wave signals from the phase transition are potentially discoverable with future space-based interferometers.


Introduction
The Standard Model (SM) of particle physics has been completed since the discovery of the Higgs boson at the LHC in 2012 [1,2]. However, it is widely believed that new physics is required to explain various phenomena beyond SM, such as the existence of dark matter (DM) and the origin of neutrino mass and mixing. We concern ourselves in this work two prominent particle physics puzzles. One is the origin of baryon number asymmetry in the Universe [3]. In the content of electroweak (EW) baryogenesis [4][5][6][7], the baryon number asymmetry is generated outside the EW broken phase and then captured by the bubbles' expansion in the progress of electroweak phase transition (EWPT). Furthermore, a sufficiently strong first-order EWPT is required to suppress the washout of baryon number asymmetry through sphalerons [7,8]. However, the EWPT in the SM with the observed Higgs boson mass is found to be a crossover [9]. The other is the existence of DM. There is overwhelming evidence from cosmological and astrophysical observations that show about 85% of the matter in the Universe is DM, instead of the normal matter made from SM particles [3,[10][11][12]. The most popular class of DM candidates is that of weakly interacting massive particles (WIMPs), which lies in the mass range of 1 − 1000 GeV [13,14]. These particles decouple from the thermal bath as the early Universe expands and cools and finally reach the appropriate relic density as observed now.
One of the simplest models to trigger a strong first-order EWPT is to introduce a new bosonic degree of freedom to the scalar potential of SM. Various Higgs portal models with (or without) a DM candidate have been widely studied . The extension of SM with a complex singlet scalar for the above-mentioned purpose was proposed and studied in Refs. [15,16], followed by extensive researches on its implications in subsequent works [16][17][18][19]. In the Higgs-portal model with a real singlet scalar as the DM candidate, the scalar is stabilized by a Z 2 symmetry and cannot have a vacuum expectation value (VEV) [44]. DM phenomena related to this model have been studied in Refs. [20][21][22][23][24][25][26][27]. A previous study [45] also shows that such a model cannot trigger a sufficiently strong first-order EWPT except when a large number of singlet scalars (∼ 10) are introduced. In the model with a global U (1) group explicitly broken down to Z 2 symmetry, Ref. [46] shows that all phase transitions leading to the correct vacuum are of the second order. Furthermore, Ref. [47] shows that a large portion of the parameter space for the Z 2 symmetry model will be ruled out by the required bubble nucleation condition, S(T n )/T n ∼ 140.
By now, no obvious evidence for WIMPs has been observed in both DM direct detections and collider searches. In particular, the recent constraints from LUX [48], PandaX-II [49], and XENON1T [50] tell us that the DM interaction cross section with protons and neutrons is extraordinarily tiny, less than about 10 −45 − 10 −46 cm 2 , low enough to rule out most of WIMP DM models. Ref. [51] shows that there is a novel mechanism for suppressing the direct detection cross section if the global U (1) symmetry of Higgs portal model is softly broken to Z 2 symmetry by a mass term. Such a suppression is due to a cancellation in the DM-nucleon scattering amplitude between the SM-like and new Higgs boson contributions at zero momentum transfer. Although the cancellation is spoiled at loop level, the scattering contribution from one loop is trivial and the conclusion remains practically unchanged [52,53].
In this work we are concerned with one of the simplest extensions of SM with a complex singlet scalar. In the model, the U (1) symmetry is softly broken by a cubic term to Z 3 symmetry (see also recent works [42,43]). The real part of the complex scalar acquires a VEV, while the imaginary component (pseudo-Goldstone boson) plays the role of DM candidate due to a residual Z 2 symmetry in the scalar potential. As shown in a previous study [44], such a cubic term could generally induce a potential barrier at tree level and trigger a strong first-order EWPT. We find that in the two-step EWPT scenario and for the heavy scalar mass falling in the range of 1.2 − 2 TeV, the model with a Higgs mixing angle satisfying |θ| 0.2 (11.5 • ) can induce a strong first-order EWPT and the stochastic gravitational wave (GW) signals produced from the phase transition are potentially discoverable by future space-based interferometers, such as LISA [54], DECIGO [55] and BBO [56].
This work is presented as follows. In sections 2 and 3, we describe our model and study detailed properties of the model at tree level. In section 4, we present our search scheme and results for the strength of the EWPT from a comprehensive scan of the parameter space. We take into account the constraints from electroweak precision observables and Higgs searches at the colliders in section 5. The DM phenomenology of the pseudo-Goldstone boson χ is studied in section 6. The discussions of GW production along with the firstorder cosmological phase transition and its detection by future space-based interferometers are given in section 7. We summarize our findings in section 8. Some detailed formulas and parameters are collected in the appendices: appendix A gives the field-dependent mass matrices of the particles at finite temperature; appendix B provides the coefficients of counter-terms in the scalar potential; appendix C gives the partial decay widths of the SM-like Higgs and heavy scalar; and appendix D discusses the sensitivities of space-based interferometers.

The model with Z 3 symmetry
We consider an extension of the SM with just a complex gauge-singlet scalar field S which transforms under a global Z 3 transformation as S → exp(i2π/3)S. Imposing the Z 3 symmmetry associated with S, we can write down the most general renormalizable scalar potential with softly U (1) symmetry breaking: where H denotes the SM Higgs doublet and µ 3 is assumed to be real. The symmetry of global U (1) transformation S → exp(iϑ)S (ϑ is an arbitrary phase) is softly broken by the µ 3 term to the Z 3 symmetry, i.e., the potential remains unchanged only under those transformations with a rotation angle ϑ = 2nπ/3, where n is an integer. Notice that the Hermiticity of the potential implies a symmetry under the transformation S → S * , which turns into a Z 2 symmetry for the imaginary component χ of S (χ → −χ). This ensures the stability of χ and makes it a DM candidate. It is worth mentioning that the usual Z 2 symmetry models do not prohibit a term proportional to |H| 2 S 2 , which could also softly break the U (1) symmetry but at the same time spoil the cancellation mechanism [51]. Such a term is not allowed by the Z 3 symmetry in our model, which shows an advantage over the Z 2 symmetry models. Furthermore, this tree-level potential boasts analytical solutions which may explicitly reveal some of properties of the model.
We note in passing that the breakdown of a discrete symmetry during EWPT in the early universe can potentially lead to the problematic EW-scale cosmic domain walls [57][58][59] whose gravitational effects may result in unacceptable anisotropy in the cosmic microwave background (CMB) radiation [57]. Depending on the stability and evolution of such domain walls, several mechanisms have been proposed to avoid these quandaries [60][61][62]. Details of this issue is beyond the scope of our current study.

Field-dependent masses
The Higgs and singlet scalar can be expanded around their classical backgrounds as where G ± , G 0 , and χ are the Goldstone bosons after spontaneous symmetry breaking. At zero temperature, the VEVs for the two scalars are v ≡ H | T =0 and w ≡ S | T =0 , respectively. We immediately obtain the tree-level potential in terms of the fields h and s Then the field-dependent mass matrix of the scalar bosons is given by The tree-level potential can now be rewritten as where Φ † = (h, s). With an orthogonal rotation from the (h, s) T basis to the physical basis (H, S), one finds the physical masses given by (2.7) The field-dependent mass of pseudo-Goldstone boson χ is Other field-dependent masses of the SM particles are given in appendix A.

Stationary points
We now search for the local minima of the tree-level scalar potential. An interesting scenario is when the potential has two local minima: the EW symmetry broken one located at (v, w) and the EW symmetric one located at (0, w 0 ). The tadpole conditions of the potential are Besides, vacuum stability demands the following bounds on parameters in the potential [63,64] (2.11)

Stationary points along h-axis
The stationary point along the h-axis is given by The condition for a physical vacuum, i.e., µ 2 h /λ h > 0, can be easily satisfied since µ 2 h > 0 holds for most of the parameter space and λ h > 0 is guaranteed by the vacuum stability. However, a zero VEV for the singlet scalar will lead to a vanishing DM mass, which is of no interest to us in this work. To avoid such stationary points from being local minima, we can demand ∂ 2 V /∂s 2 < 0 at these points, giving As we will see below, this condition also ensures a stationary point at (h = 0, s = 0).

Stationary points along s-axis
The stationary points along the s-axis is given by (2.14) The required physical condition is µ 2 3 − 4λ s µ 2 s > 0. We find that for most of the parameter space having sufficiently strong first-order EWPT, the condition µ 2 s < 0 always holds. As we will see below, for the stationary points sitting on the s-axis, the potential located at s + is always lower than the potential located at s − if µ 3 < 0 (i.e., V (0, s + ) < V (0, s − )). On the other hand,

Stationary points at
There are also solutions off the h-axis and s-axis, given by The condition for these solutions to be physical is With vacuum stability condition (2.11), we have λ s − λ 2 m 4λ h > 0. As mentioned above, we can further impose the condition (2.13) on the potential if we demand that the stationary points along the h-axis be unstable. Hence, the condition (2.16) is satisfied for most of the parameter space of interest.
Using Eq. (2.10), we obtain One can show that F < 0 is always true under the assumption of µ 3 < 0. This is in fact equivalent to the condition Therefore, we finally reach the conclusion that with the assumption µ 3 < 0, the tree-level potential located at (h, s + ) is always lower than the one located at (h, s − ). Similarly, We can also prove that F > 0 is always true provided µ 3 > 0. We summarize our conclusion as follows: One can easily verify that these conclusions are also established for the case of local minima along the s-axis. From Eq. (2.8), we see that the pseudo-Goldstone DM mass is To avoid a tachyonic mass for the DM candidate, the signs of µ 3 and w should be opposite. Without loss of generality, the singlet scalar's VEV of is assumed to be positive, and µ 3 should thus have a negative value.

Parameters
Using Eqs. (2.4), (2.9), and (2.10) we have The reason for employing Eqs. (2.9) and (2.10) is to ensure the existence of a local minimum at (v, w) for any choice of parameters. The above parameters can be related to three physical parameters, the masses of two Higgs bosons m H and m S and the mixing angle θ, with the relations where m H = 125 GeV and v = 246 GeV. Thus, we take {w, m S , m χ , θ} as the input parameters of the model.

Effective potential
At the one-loop level, the total effective potential is given by where the tree-level potential V (h, s) has been given above, and the other components are discussed below.
At zero temperature, the one-loop corrections to the potential is given using the MS renormalization scheme [65] as where the subscript i = {H, G, S, χ, Z T , Z L , W T , W L , t, b} denote respectively the SM-like Higgs boson H, SM Nambu-Goldstone bosons, heavy scalar S, pseudo-Goldstone DM χ, transverse and longitudinal components of SM gauge bosons, and top and bottom quarks, and N i = {1, 3, 1, 1, 2, 1, 4, 2, −12, −12}. Q is the renormalization scale which is set to be v in this work. The constant C i = 1/2 for gauge boson transverse modes and 3/2 for all the other particles. Due to plasma damping, the validity of the perturbative expansion of the effective potential breaks down at high temperatures. A remedy to this problem is to resum the daisy diagrams to all orders, which results in an additional contribution to the bosonic masses [66]. Thus, we replace the field-dependent bosonic masses at finite temperatures by where the thermal corrections Π i (T ) are given in appendix A. The one-loop potential becomes gauge-dependent when thermal corrections of bosons' masses are introduced [67,68], leading to gauge-dependent critical temperature and GW spectrum produced from phase transition [34,69]. To focus on our topic, in this work we use the Landau gauge with a vanishing gauge-fixing parameter for the potential. At the one-loop level, the finite-temperature contributions to the effective potential are given by [64,70] with the − sign for bosons and + for fermions. To maintain the main properties of the tree-level potential derived above, we add the following counter-terms to the potential at zero temperature The coefficients of the counter-term potential are given in appendix B.

Parameter space for electroweak phase transition
In this section, we scan the parameter space for viable sample points for a sufficiently strong first-order phase transition. We also show the distributions of the model parameters and physical parameters based upon our scan results. T > Tc

Two-step phase transition
A strong first-order EWPT could occur if there is a sufficiently high and wide potential barrier separating the two degenerate vacua of the thermal effective potential at critical temperature. Introducing an extra bosonic degree of freedom could enhance the barrier and thus make the EWPT stronger. In our model, there are two contributions to the barrier in the effective potential [30]: one is the tree-level barrier coming from the cubic term of the tree-level potential; the other one arises from the bosonic thermal corrections to the potential at the one-loop level. To see the latter, one can expand the integrations (3.4) in the high temperature limit, i.e., z ≡ M 2 i (h, s)/T 2 1, and find that there is a term proportional to z 3/2 that leads to terms cubic in both h and s. Most importantly, the barrier term cubic in h is proportional to the Higgs portal coupling λ m , and thus to the mixing angle θ (see Eq. (2.23)) [71]. Therefore a large value of mixing angle could strengthen the EWPT by boosting the potential barrier. The potential barrier is shallow in the SM shallow and its EWPT is confirmed to be a crossover [9].
In this work, we focus on the so-called two-step phase transition as shown in Fig. 1. At very high temperatures, both scalar fields have no VEV. When the temperature drops to a certain value T s , (above the critical temperature T c , at which two degenerate vacua exist concurrently) the scalar potential along the s direction develops a global minimum at (0, w 0 (T s )). We call it the symmetric phase, as shown in the upper panel of Fig. 1. As the temperature further lowers to the critical temperature, another local minima located at (v(T c ), w(T c )), designating the symmetry-broken phase, appears and becomes degenerate with the symmetric phase (0, w 0 (T c )), as shown in the lower left panel of Fig. 1. Meanwhile, a tunneling path, along which is a lowest barrier between the two degenerate minima, opens up. The decrease of potential at the broken phase is much faster than that at the symmetric phase as the temperature approaches zero. As a result, the broken phase moves to (v, w) and becomes a global minimum, as shown by the lower right panel of Fig. 1.

Searching scheme
The strength of the EWPT is measured according to the order parameter is the VEV of SM-like Higgs boson at the critical temperature. For a successful baryogenesis, the first-order EWPT should be strong enough so that the sphaleron process in the broken phase is sufficiently suppressed to avoid baryon asymmetry washout [8]. This gives the conventional criteria for a sufficiently strong EWPT: The precision value of bound should be derived by calculating the sphaleron energy for a specific model. However, the exact values are often very close to 1 and we will neglect such a small difference. Aiming at the search of parameter space giving a sufficiently strong EWPT, we make a random scan of the parameters in the following ranges: In addition to the above enforced restrictions, the parameters are also subject to the constraints discussed in section 2.2, including vacuum stability (2.11) and the conditions to ensure the existence of stationary points at the symmetric phase (0, w 0 ) and the broken phase (v, w). Further constraints come from the requirement of perturbative unitarity [64,72] We also require that V (v, w) < V (0, w 0 ) and V (v, w) < V (0, 0) at zero temperature, so that the broken phase is a global minimum.
To search for the critical temperature where two degenerate vacua coexist, we start from an initial temperature T between a minimum value of temperature T min and a maximum value of temperature T max , we then search between these two temperatures for the local minima of the potential around the positions (0, w 0 ) and (v, w), which can be determined using the analytical formulas (2.14) and (2.15) for given parameters. If the local minimum at the symmetric phase (0, w 0 (T )) is found to be larger (smaller) than the one at the broken phase (v(T ), w(T )), the temperature is increased (decreased) in the next trial. We conclude that no electroweak phase transition for given parameters occurs if the following two cases are met: • Case 1. The local minimum at the symmetric phase (0, w 0 (T max )) is larger than the one at the broken phase (v(T max ), w(T max )).
• Case 2. The local minimum at the symmetric phase (0, w 0 (T min )) is less than the one at the broken phase (v(T min ), w(T min )).
As obvious, the lower and upper temperature limits are critical for our parameters searches. In our preliminary scan of the parameter space, we find that no points with v c /T c 1 can be found at a temperature higher than about 350 GeV. We thus restrict the scan range of temperature to (10 − 350) GeV.

Parameter distributions
We generate one million random floats uniformly for each of the input parameters, among which about 1.8% are found to be able to trigger a sufficiently strong phase transition while fulfilling the other basic criteria mentioned above. We show the distributions of various physical parameters in Fig. 2, and summarize our observations from the figures as follows: 1. There is a lower bound in the distribution of m S , i.e., m S 500 GeV. This is because in order to trigger a phase transition, the local minimum at the broken phase (v(T ), w(T )) should be larger than the local minimum at the symmetric phase (0, w 0 (T )) when the temperature is higher than critical temperature. To see this, we show in the upper left panel of Fig. 3 , T ) at temperature T = 300 GeV (the distribution of critical temperature can be found in Fig. 4) as a function of m S , with θ = 0.4 and m χ = 300 GeV. The blue, green, and red points in Fig. 3 represent the results with w = 300 GeV, 400 GeV, and 500 GeV, respectively. As shown in the drawing, the potential difference slowly increases with m S and remains negative for m S 500 GeV. At larger values of m S , the potential difference sharply increases to a positive value. Therefore, the lower bound m S 500 GeV is to guarantee a successful phase transition.
2. There is a lower bound in the distribution of |θ|, i.e., |θ| 0.2. Such a relatively large lower bound on the absolute value of the mixing angle is straightforwardly derived from the requirement of a sufficiently strong first-order EWPT. As shown in Eq. (2.23), the mixing angle |θ| directly controls the Higgs-portal interacting strength |λ m | (whose distribution can be found in the left panel of the second row of Fig. 4). the tunneling path to be along a linear combination of the two fields (see Fig. 1) and finally lead to a strong EWPT (see also Ref. [73] for similar conclusions).
3. The lower bound on w and the decrease in the distribution of |θ| at larger values are from the bound of unitarity (last condition of Eq. (4.3)). We plot this bound condition in various planes in Fig. 3, the colored regions are excluded by the constraints. In the upper right panel of Fig. 3, we show the unitarity bounds in the θ-m S plane, fixing m χ = 300 GeV. The values of w in this panel are 100 GeV (green), 200 GeV (yellow), 300 GeV (blue), 400 GeV (light blue), 600 GeV (pink), and 1000 GeV (sandy-brown), respectively. As shown in this panel, when w = 100 GeV, the heavy scalar mass m S is restricted to lie below about 500 GeV. As one increases w, more parameter space opens up. Sufficient sample points for strong EWPT are available for w 200 GeV. These indicate that the unitarity bound on large values of m S can be avoided by increasing w. In the lower left panel of Fig. 3, we plot the unitarity bounds in the wm S plane, with θ = 0.2 (blue), 0.4 (pink), and 0.6 (green), respectively. As shown in this panel, the available parameter space is largely reduced with the increase mixing angle, which explains the decreasing behavior in the |θ| distribution at large values. When m S 500 GeV, nearly all of the parameter space is excluded by the unitarity bound for w 200 GeV. The constraints become independent of m S for w 500 GeV and θ 0.4.

4.
The distribution of pseudo-Goldstone DM mass m χ has a peak around (200−400) GeV and then decreases at larger values. The decreasing behavior of m χ distribution involves a few contributions. One is the unitarity bound; it is easy to check that the unitarity bound becomes tighter for a larger DM mass m χ . Another reason is the condition (2.13), which is to ensure that the stationary point located along the h axis cannot be a local minimum. We show this constraints in the lower right panel of Fig. 3, the colored regions are excluded by the constraints. In this plot, the value of w is fixed at 1000 GeV, the green, yellow, pink, and blue regions represent the bounds with m S = 500 GeV, 700 GeV, 1000 GeV, and 1500 GeV, respectively. We see that for θ 0.2, the parameter space with m χ 1000 GeV is excluded. However, we should note that this bound is imposed only when there is a stationary point along h axis, which requires µ 2 h > 0. From the distribution of µ 2 h shown in the upper left panel of Fig. 4, we find that only about half of the total sample points have a stationary point along h axis. We also note that such a requirement is somehow too stringent because of the loop corrections on the potential. Further constraint may arise from the requirement of d∆V bs (T )/dT 2 > 0 near the critical temperature [63], we find this bound is relatively weak and do not give a detailed account here. 5. We see that within the parameter ranges considered here, the distributions of w and m S do not decrease with increasing values. Such distributions are consistent with a recent global fit performed in Ref. [32].
In Fig. 4 we show the distributions of the various parameters in the scalar potential, which are determined from the four physical parameters above using Eqs.

Experimental constraints
In this section, we consider various constraints of collider experiments on our model, including the electroweak precision observable measurements and Higgs signal strength measurements.

Electroweak precision observables
We first consider the constraints arising from electroweak precision observables (EWPOs) [74,75]. We introduced two scalars to the SM, one of them is a real singlet scalar which can mix with the SM Higgs boson. Both the SM Higgs boson and the real singlet scalar will contribute to the SM gauge boson self-energies at loop-level and finally induce corrections to the oblique parameters S, T , and U . The other one is a pseudo-Goldstone boson χ which becomes a DM candidate due to an accidental Z 2 symmetry hidden in our Z 3 symmetry model. The pseudo-Goldstone boson do not mix with SM particles and thus do not affect the oblique parameters. Constraints from the EWPO's can alter the distribution of m χ via the parameters m S and θ. The contributions of the heavy scalar to the oblique parameters ∆O i ≡ O i − O SM i (here i denotes T , S, or U ) can be quantified as follows [28,44]: where m Z and m W are the SM weak gauge boson masses, the cosine of the weak mixing angle c 2 The loop functions f T (x) and f S (x) are given by [32,76] In order to analyze the impacts of the EWPO's on the strong EWPT parameter distributions, we follow the procedure given in Refs. [44,73,77,78] by defining a ∆χ 2 as where ∆O exp i denotes the experimental measurements of the deviations of oblique parameters from its SM reference values. We take the most recent analysis from the Gfitter Group [79] ∆S = 0.04 ± 0.11, ∆T = 0.09 ± 0.14, ∆U = −0.02 ± 0.11. (5.4) The covariance matrix Σ 2 ij ≡ σ i ρ ij σ j involves σ i as the various errors given in Eqs. is the correlation matrix of the experiment. The electroweak observables are governed by only two parameters: the heavy scalar mass m S and the mixing angle θ. We consider the singlet scalar extended models to be consistent with the EWPO's if the oblique parameters lie within the 95% confidence level (CL) ellipsoid, which corresponds to taking ∆χ 2 ewpo ≤ 5.99 for the models with given m S and θ.
We plot the contours of ∆χ 2 ewpo in Fig. 5 and find that in the m S -θ space, the region with m S 1500 GeV and θ 1.0 (to the right of the pink curve) is excluded at the 95% CL. The constraint is seen to be not very strong.

Higgs signal strengths
In our extension of the SM with a complex singlet scalar, the SM-like Higgs boson coupling strengths to other SM particles are modified by a common factor of cos θ. This leads to the productions of XX (where X can be a SM gauge boson or fermion) via the SM-like Higgs boson being suppressed by a factor cos 2 θ, provided that no additional decay channel is permitted. Note that the same factor of cos 2 θ is involved in all major production channels (gluon fusion, vector boson fusion and tth). To take into account the published exclusion bounds from Higgs searches at the LEP, Tevatron and LHC experiments, we make use of the HiggsBounds_v4.3.1 [80] package, which calculates the predicted signal rates for the search channels considered in the experimental data. More information on search channels and experimental data used in HiggsBounds can be found in Ref. [80]. The HiggsBounds package determines whether a point in the model parameter space is excluded at 95% CL by comparing the predicted signal rates against the expected and observed cross section limits from the direct Higgs searches [80].
Distinct signal strengths, defined as the production rate times the decay branching fraction relative to the SM expectation, i.e., µ i ≡ (σ × BR) i /(σ × BR) SM i , in various decay channels including γγ, W W * , ZZ * , bb, and τ + τ − have already been measured with high precision at the LHC. From these signal strengths, one can obtain information on the couplings of the Higgs boson to SM particles and derive constraints on the extension models. For the model considered in this work and in the narrow width approximation, the signal strength of H is [28] where Γ SM H denotes the total decay widths of the SM-like Higgs boson with mass being set to m H , and Γ H→XX is the width of H decaying to a pair of X (= χ, H), which can be found in appendix C. We see that the signal strengths of H is suppressed by two factors: cos 2 θ and the presence of new decay channels. Even when the new decay channels are kinematically forbidden, the signal strength is still reduced by cos 2 θ. This means that a generic signature of the mixing of the SM Higgs boson with an extra singlet scalar boson can be derived from a reduced signal of the Higgs bosons at the LHC [28].
We use the HiggsSignals_v1.4.0 package [81] to estimate the χ 2 of a given model and assess which sample points are allowed by the signal strength measurements. The package assumes a Gaussian probability distribution and uses the peak-centered method for the calculation of χ 2 = χ 2 µ + χ 2 m , where χ 2 µ is evaluated by comparing the signal strength measurements for the peak to the model-predicted signal strengths and χ 2 m is evaluated by comparing the model-predicted Higgs boson mass and the observed one if a mass measurement is also available [81]. The signal strength measurements used in the HiggsSignals analysis are summarized in table 1.
We find that χ 2 min , the minimum of χ 2 , is obtained in the model with the mixing angle θ = 0, which means no mixing for the SM Higgs boson. Totally 89 observables are used in each fit, and χ 2 min = 102.02. We conclude that a sample point in the parameter space is not excluded by the experimental results at 95% CL if ∆χ 2 hs = χ 2 − χ 2 min < 9.49 (for four free parameters).
We plot ∆χ 2 hs as a function of the mixing angle in Fig. 6, with the pink, blue, and green curves denoting the results for (m S , m χ ) = (800, 300) GeV, (800, 50) GeV, and (50, 50) GeV, respectively. As discussed above, when m S , m χ > m H /2, the SM Higgs signal strength is reduced by a factor of cos 2 θ, and thus ∆χ 2 hs increases with |θ|. When m χ < m H /2, the SM Higgs invisible decay H → χχ is kinematically allowed, leading to an additional suppression in the SM Higgs signal strength, as is indicated by the blue curve. One can also see that the opening of new decay channels will play a dominant role in suppressing the SM Higgs signal strength. Furthermore, if m S < m H /2 the decay channel H → SS opens up and dominates ∆χ 2 hs for |θ| 0.3, as shown by the green curve. In summary, the constraints from Higgs signal strength measurements can be divided into two cases:   the Higgs signal strength scales as cos 2 θ. The mixing angle of the SM Higgs boson with an extra singlet scalar should be 0.4 (23 • ) at 95% CL.
• Case 2. If there is at least one extra particle is lighter than half of the SM Higgs boson mass, the Higgs signal strength will receive an additional suppression, the constraint on the mixing angle becomes more rigorous than Case 1.  Figure 7. Distributions of the input parameters. The blue histogram represents those samples that can trigger a strong EWPT, the green histogram represents those samples that further survive the experimental constraints.

Results
We show our results in Fig. 7 after combining all the collider experimental constraints discussed above. We find that 8104 (in the green histogram) out of 18047 sample points (in the blue histogram) that can trigger a sufficiently strong EWPT survive the experimental bounds. As discussed above, the most stringent constraints come from the Higgs signal strength measurements.
As shown by the blue histogram, all of the m S and most of the m χ sample points that induce a strong EWPT distribute at masses larger than m H /2 = 62.5 GeV. According to the above discussions, we see that there is a universal constraint on the mixing angle, |θ| 0.4, from the Higgs signal strength measurements. Here we have chosen the scan range of m S ≥ 150 GeV, apparently, the bound on the mixing angle can be much stronger if m S and/or m χ is lighter than m H /2. Anyway, the constraint |θ| 0.4 is not dependent on the assumed parameter ranges. Hence, as shown by the lower right panel of Fig. 7, there is a hard cut-off around ±0.4 in the distribution of θ.
The peak around 200 − 400 GeV in the distribution of m χ is now removed, and it tends to a flater distribution in the range 100 − 1500 GeV. While there seems to be no preferred range in the w distribution by the signal strength constraints, the sample points with m S 1.2 TeV are strongly disfavored, as shown in the top two panels. For one thing, the signal strength bounds require a mixing angle |θ| 0.4. Yet samples with large values of mixing angle are mainly associated with a lighter scalar mass m S due to the perturbative unitarity. For another reason, the scalar mass m S should be large enough to induce a strong EWPT when the mixing angle is small.
We summarize our main conclusions obtained in this section below: • The Higgs signal strength measurements give a universal constraint on the mixing angle |θ| 0.4 (23 • ).
• The mass of heavy scalar S should be larger than 1.2 TeV from the combined constraints of Higgs signal strength measurements and perturbative unitarity and the requirement of strong EWPT.
• Our analysis supports |θ| 0.2 (11.5 • ) for the scalar mass m S 2 TeV, which is the requirement of a sufficiently strong EWPT. However, this conclusion depends on the scanning range of scalar mass m S pre-assumed in this study. A strong EWPT for the mixing angle less than 0.2 might be available if the scanned heavy scalar mass extends beyond 2 TeV (we leave this for future studies).

Dark matter phenomenology
In this section, we discuss constraints on the model from the observed DM relic density and null direct search result.

Dark matter relic density
As mentioned above, the pseudo-Goldstone boson χ from the spontaneous symmetry breaking has a Z 2 symmetry, ensuring the stability of the pseudo-Goldstone boson as a DM candidate. In the standard freeze-out scenario, the DM particles are in chemical equilibrium with the other SM particles via annihilation-production reaction in the early Universe. The DM population becomes nonrelativistic and the annihilations take over the thermal productions with the adiabatic expansion of the Universe. At around the temperature that the reaction rate falls below the Universe expansion rate, the DM particles begin to decouple from the thermal bath. The evolution of the DM number density is described by the Boltzmann equation where the abundance Y (T ) denotes the ratio of the DM number density n χ to the entropy density, M pl = 1.22 × 10 19 GeV is the Planck mass, g * is the effective number of relativistic degrees of freedom, Y eq (T ) is the thermal equilibrium abundance, and σv rel is the relativistic thermally averaged annihilation cross section. The resulting DM relic density is given by where Y 0 is the abundance of DM in the present Universe. In our numerical analysis, we make use of the MicrOMEGAs_5.0.4 package [87] to calculate the DM relic density.
In the left panel of Fig. 8, we plot the DM relic density h 2 Ω DM as a function of the DM mass m χ , with the values of w and m S fixed at 1000 GeV and 800 GeV, respectively. In our model, the DM annihilation to the SM particles are mediated by the SM-like Higgs boson H and heavy scalar S. When m χ m H /2, its annihilation process is kinematically suppressed, leading to a large value of DM relic density. The resonant DM annihilation occurs at m χ m H /2 and m S /2, which would result in a sharp decrease of the DM relic density, as shown by the two dips in the curves. In the right panel of Fig. 8, we plot the contours of DM relic density in the m S -w plane. We see that the DM relic density becomes larger as w increases.
In the left panel of Fig. 9, we calculate the DM relic density of the sample points. The blue scatter points represent the samples that have a sufficiently strong EWPT and the green scatter points are the samples further survive all the experimental constraints. The black line denotes the DM relic density given by the Planck satellite's observation of the CMB radiation [88] h 2 Ω obs DM = 0.1188 ± 0.0010.
It is seen that most of the samples have h 2 Ω DM much below the observed result. On one hand, these sample parameters can evade from a over-closed Universe. On the other hand, only a small fraction of DM consists of the pseudo-Goldstone boson from spontaneous symmetry breaking in our model.

Direct detection
For DM direct detection, we use the MicrOMEGAs package to compute the spin-independent DM-nucleon elastic scattering cross section σ SI . As shown above, for the parameter space considered here, only a small fraction of DM is made up of the pseudo-Goldstone bosons. Effective spin-independent DM-nucleon elastic scattering cross section as a function of the DM mass. The black and red curves denote the upper limits on scattering cross section from LUX [48] and XENON1T [50] experiments, respectively. In both plots, the blue scatter points represent the samples that have a sufficiently strong EWPT and the green scatter points are those further surviving the experimental constraints.
To compare with the experimental upper bounds, we have to scale the scattering cross section asσ We simultaneously calculate the DM relic density h 2 Ω DM and scattering cross section σ SI with the help of MicrOMEGAs, then we obtain the effective scattering cross sectionσ SI using Eq. (6.4). We plot our results in the right panel of Fig. 9. Again, the blue scatter points represent the samples that have a sufficiently strong EWPT, and the green scatter points are the samples further surviving all the collider constraints. The black curve denotes the upper limit on DM-nucleon elastic scattering cross section from the LUX [48] experiment, and the red curve from the XENON1T [50] experiment. They are the most stringent constraints on the spin-independent DM-nucleon scattering cross section up to date. The scattering cross section in our model is suppressed by both a small mixing angle θ and a large value of w. It can be seen that our scanned sample points have sustained the most stringent constraints from DM direct detection experiments with most of the scattering cross sections being much below the upper limits from the direct detections.

Gravitational waves
A first-order cosmological phase transition can only occur in the presence of a scalar effective potential barrier separating the symmetry-broken and -unbroken vacua. Although the probability of tunneling from the metastable minimum to the stable one via the instantons is very tiny, the decay of the false vacuum can proceed through thermal fluctuations which help overcome the potential barrier. The tunneling rate per unit volume and time element is approximately given by [45,89] Γ(T ) = A(T )e −S 3 /T , (7.1) where A(T ) [S 3 /(2πT )] 3/2 T 4 and S 3 denotes the three-dimensional on-shell Euclidean action of a instanton. Due to the supercooling effect, the onset of the bubble nucleations can be delayed to a temperature T n smaller than the critical temperature T c . The nucleation temperature T n is defined to be at which the probability of nucleating one bubble per horizon volume is of order one, i.e., p(T ) ∼ 1, where the probability of bubble nucleations per Hubble volume is defined as In a radiation dominated Universe, the Hubble parameter is given by H 2 = 8π 3 g * T 4 /(90M 3 pl ) and g * 110. The condition p(T ) ∼ 1 guarantees the percolation of bubbles in the early Universe and can be translated into the following criterion for determining the nucleation temperature [45] For the phase transition at a characteristic temperature T ∼ O(100 GeV), the condition above is well approximated by S 3 (T n )/T n 140. The successful bubble nucleations at EW scale are guaranteed by Eq. (7.3), which requires a sufficiently large bubble nucleation rate to overcome the expansion rate. On one hand, a sufficiently strong EWPT ensures that the washout of baryon asymmetry through sphalerons is suppressed. On the other hand, a successful bubble nucleation is the requirement of triggering baryogenesis in the EW broken phase. The latter is typically a more stringent requirement on the model.
The Euclidean action for a spherical bubble configuration can be written as By extremizing the Euclidean action, we obtain the following differential equation If Φ(r) represents the profile of a particle in the potential V , then Eq. (7.5) can be treated as the classical equation of motion, which can be solved by the traditional overshooting/undershooting method [89]. In this work, we employ the CosmoTransitions 2.0.2 pack-age [90] to perform the numerical calculations of the bubble profile and Euclidean action. Afterwards, we use Eq. (7.3) to determine the nucleation temperature T n .

Gravitational wave parameters
It has been shown that the stochastic gravitational waves (GW's) produced from a cosmological phase transition can be fully characterized by the knowledge of two primary parameters [91], which are defined as and where T * is the GW generation temperature, ρ rad = π 2 g * T 4 /30 is the radiation energy density in the plasma, and the latent heat associated with the phase transition is given by is the potential difference between the broken phase and the symmetric phase at temperature T . Therefore, the parameter α is related to the maximum available energy budget for gravitational wave emissions. The parameter β represents the rate of time variation of the nucleation rate, whose inverse gives the duration of the bubble nucleation. Consequently, β/H * defines the characteristic frequency of the GW spectrum produced from the phase transition.
In addition to the GW parameters α and β/H * and the nucleation temperature T n , the GW spectrum also depends on the bubble wall velocity v w , which is the expanding speed of the true vacuum. It has been pointed out in Refs. [92,93] that it is the relative velocity between the wall and the plasma in the front (v + ) rather than v w that should be used in the calculations of EW baryogenesis. For a strong EWPT, the condition v + v w could be achieved, making it possible to generate a viable EW baryogenesis and a loud GW signal in the same scenario. In this work, the bubble velocity v w is simply assigned to be around 1 for the calculations of GW spectra.
We show the calculation results of parameters α and β/H * in the left panel of Fig. 10, the colored bar indicates the nucleation temperature T n . The distribution of T n are given in the right panel of Fig. 10. We find that 8564 out of 18047 sample points satisfy the bubble nucleation condition in Eq. (7.3). Note that in the calculations of β/H * with Eq. (7.7), we have assumed that reheating is not significant for typical transitions. In this case, the temperature for GW generation T * is approximately equivalent to the nucleation temperature T n , i.e., T * T n . A strong supercooling could not only enhance the strength of the phase transition, but also change the evolution of the Universe since the expansion of the Universe would be dominated by vacuum energy in the supercooled phase, instead of radiation. In this case, there is a lower bound on the temperature of the transition to ensure the successful bubble percolation and completion of the EWPT [94].

Gravitational wave spectrum
In what follows we review the three processes that are involved in the production of GW's during a first-order phase transition (see Refs. [95,96] and references therein for details): • Collisions of bubble walls and shocks in the plasma. GW's produced from this process depends only on the dynamics of the scalar field. The "envelope approximation" is used in the numerical simulations to estimate the GW spectrum, given by [97]  • Sound waves in the plasma generated subsequently after the bubble collisions. Numerical simulations indicate that the durations of sound waves and turbulence as active sources of GW's are typically much longer than the collisions of the bubble walls. This process contributes a GW spectrum desribed by [98] • Turbulence in the plasma formation after the bubble collisions. Simulations show that only a small fraction ∼ 5 − 10% of the bulk motion from the bubble walls is converted into turbulence. However, GW's produced from this process could play a dominant role at high frequencies, as the GW signals from sound waves decay much faster. The GW spectrum from turbulence can be parameterized as [99] h 2 Ω turb (f ) = 3.35 × 10 −4 H * β κ turb α 1 + α There is a critical phase transition strength α ∞ for the phase transition, which is estimated as where n i is equal to N i for bosons and N i /2 for fermions, with N i already given in section 3.
is the difference of the field-dependent squared masses between the symmetric and broken phases. According to α ∞ , the phase transition can be divided into two cases [95]: • Case 1. Non-runaway bubbles, α < α ∞ . In this case, the bubble walls will reach a terminal velocity and the latent energy transferred into the scalar field is negligible, i.e., κ col 0. The efficiency factor for the sound wave contribution is then given by The efficiency factor for turbulence κ turb is related to κ sw by κ turb = κ sw , where we take = 0.1 in this work.
A recent study [100] suggests that although bubbles may runaway in certain cases, most of their energy is dissipated into the surrounding plasma and very little energy is deposited in the bubble walls. This in turn leads to a negligible contribution to the GW spectrum from bubble collisions. Hence, we only take the GW spectra produced by sound waves and turbulence into account in our calculations.

Space-based interferometers
The frequentist approach is normally used for the experimental investigation of the stochastic GW signals from EWPT, where the detectability of the signals is measured by the corresponding signal-to-noise ratio (SNR) [37,95] where h 2 Ω exp denotes the sensitivity of a GW experiment, N is the number of independent observatories of the experiment, and T obs is the duration of the mission in units of year. The peak frequency of GW spectrum produced from the EWPT is red-shifted to around the milli-Hertz band, which falls right within the range of future space-based GW interferometers. The planned space-based GW experiments considered in this work include the LISA [54] interferometer as well as the proposed successors B-DECIGO [101], DECIGO [55], and BBO [56]. For the auto-correlated experiments LISA and B-DECIGO, N = 1, while N = 2 for the cross-correlated experiments DECIGO and BBO. Following Ref. [37], we assume a mission duration of T obs = 4 years for all of the experiments. The SNR threshold value ρ thr , above which the GW signal is detectable for the experiment, and the other detector parameters are summarized in table 2. The experimental sensitivities are summarized in appendix D.
In Fig. 11, the GW experimental sensitivities are shown in the α-β/H * plane by the colored regions, in which the SNR of a given GW observatory exceeds its detection threshold. We have assumed a nucleation temperature T n = 200 GeV for the determination of these regions. The scatter points are the samples that both generate a strong first-order EWPT and survive the collider searches and DM experiments. The colored bar shows values of the mixing angle θ. The samples with θ > 0 (θ < 0) are plotted in the left (right) panel of Fig. 11). As illustrated by this figure, a considerable portion of sample points with a mixing angle |θ| in the range 0. 25 Table 2. A summary of parameters and assumptions used for the planned space-based interferometers. the LISA experiment. We thus expect that the EWPT scenario depicted in this work will be tested by the future space-based interferometers.

Summary and conclusion
In this work, we have considered the extension of SM with a complex singlet scalar field S. A global U (1) symmetry associated with S in the scalar potential is softly broken by its cubic terms to a Z 3 symmetry. The real part of the complex scalar field develops a VEV and mixes with the SM Higgs boson after the electroweak symmetry breaking, while the imaginary component, the pseudo-Goldstone boson χ, becomes a DM candidate due to an accidental Z 2 symmetry hidden in the scalar potential. We focus on the two-step phase transition scenario, in which the location of global minimum of the potential moves as (0, 0) → (0, w 0 (T )) → (v(T ), w(T )) with the expansion of the Universe. We then search for the parameter space for sample points where a sufficiently strong EWPT can occur, taking into account the constraints from collider experiments and DM searches. The conclusions we have obtained are summarized as follows: • The requirement of a sufficiently strong EWPT demands |θ| 0.2 (11.5 • ) for m S 2 TeV. The Higgs signal strength measurements, on the other hand, give the constraint |θ| 0.4 (23 • ).
• To trigger a successful phase transition, the real scalar mass m S should be larger than about 500 GeV. The constraints of Higgs signal strength measurements further pushes up the lower bound of the scalar mass: m S 1.2 TeV.
• A small fraction of DM could be made of the pseudo-Goldstone boson χ, while the constraints from current DM direct searches are satisfied.
• GW signals from the first-order phase transition in the model with the mixing angle θ in the range of ∼ 0.25−0.4 can be detectable using future space-based interferometers, such as LISA, DECIGO, and BBO.
We concentrate in this work on the two-step phase transition scenario as described in section 4.1. Other types of multi-step phase transitions can be found in recent works [33,36,71]. The requirement of a large mixing angle to trigger a sufficiently strong EWPT is crucial in the two-step phase transition. The one-step phase transition, (0, 0) → (v(T ), w(T )), can be stronger, and for models with a mixing angle 0.2, the phase transition may be detectable in planned space-based interferometers [36]. We also note that in our analysis, we have restricted the heavy scalar mass m S to be in the range of 150 − 2000 GeV (m S > m H ), the strong EWPT with a lower scalar mass, m S < m H , can be found in Refs. [71,102]. The longitudinal components of W and Z receive thermal corrections from the daisy diagrams, and their masses in terms of W + µ , W − µ , W 3 µ , B 0 µ basis can be written as The field-dependent masses of top and bottom quarks are

C Decay width
Here we give the partial width formulas for some decays of SM-like Higgs H and heavy scalar S:

D Sensitivity of space-based interferometers
Here we summarize the experimental noise power spectral density used in this work. More details can be found in Ref. [37] and references therein. The LISA interferometer is planned to launch in 2034, with the noise power spectral density well approximated by [54,103] S LISA (f ) = 10 3L 2 P OMS (f ) + 2 1 + cos 2 f f * P acc (f ) (2πf ) 4 1 + 6 10 (D.8) The experimental frequency range, duration of the mission and other parameters for the future GW interferometers are summarized in table 2.