1 Introduction

The failure of the Standard Model (SM), in describing phenomena like the baryon asymmetry of the universe (BAU) and the dark matter (DM), brings to mind that the SM cannot be considered as a fundamental model. Nevertheless, the discovery of the Higgs boson [1, 2] as the first observed scalar has opened the way to consider the SM as an effective field theory (EFT) and also a window to the Higgs portal. To address the BAU and the DM problems, many theories and models have been proposed beyond the SM such as supersymmetry studies [3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]. Due to the attraction of the Higgs portal, it has been always of interest to investigate the SM extensions which directly challenge the Higgs portal like multi-scalar extensions [18,19,20,21,22,23,24,25,26,27,28,29,30]. The existence of interactions between the Higgs and new scalars makes such models reasonable for explaining the BAU, which needs a strong first-order electroweak phase transition (SFOEWPT), the gravitational waves (GW) produced by an SFOEWPT and the DM. Moreover, such models also have other benefits. First, they are simple and straightforward. Second, they may be renormalized, so no new physics scale is needed. Third, they may be gauge independent, if there exists a barrier in the potential at tree-level [31].

To justify the BAU, there is a need for Baryogenesis to exist [32,33,34,35,36] which itself needs an SFOEWPT, i.e. \(\frac{v_{c}}{T_{c}}\gtrsim 1\) where \(v_{c}\) is the Higgs vacuum expectation value (VeV) at critical temperature \(T_{c}\). This would not happen in the SM, but adding one or more new scalars to the SM potential may lead to an SFOEWPT. With regard to the new potential structure, two different phase transitions (PT) can happen. One of them is one-step PT in which there only exist initial and final phases. Cooling down the universe, it goes through a phase transition and breaks the electroweak symmetry. The other one is two-step (or multi-step) PT in which there also exists an intermediate phase (or more) between initial and final phases [37,38,39,40,41]. The reader is referred to [42,43,44,45,46,47,48,49] for the most recent studies on the EWPT.

During the SFOEWPT, the bubbles with the non-zero VeV nucleates in the plasma. The stochastic GW background arising from the SFOEWPT can be generated by the bubbles collisions and shocks [50,51,52,53,54,55], the sound waves [56,57,58,59], and the Magnetohydrodynamic (MHD) turbulence [60,61,62,63,64] in the plasma. Since the EWPT in the SM is a cross-over instead of strong one, the SM cannot predict the GW produced by the EWPT. So, this is another reason to look for beyond the SM. The recent observations of astrophysical GW [65,66,67,68,69,70,71,72] have brought the hope to detect the GW produced by the EWPT [73,74,75]. The reader is referred to [42, 44,45,46,47,48, 76,77,78,79] for the most recent studies on the GW produced by the EWPT.

As mentioned before, the SM cannot explain DM which existence is well established by cosmological evidence. As the simplest way, this incompetence can be justified by adding one (or more) gauge singlet scalar to the SM. Since the DM should be stable to provide the observed relic density \(\Omega _{c}h^{2}=0.120\pm 0.001\) by \(\mathbf {Planck} \mathbf {2018}\) [80], it is necessary to impose a discrete symmetry on the DM candidate, in present study \(S_{2} \rightarrow -S_{2}\). On the other hand, the global minimum of potential at zero temperature spontaneously breaks this discrete symmetry, so necessarily \(<S_{2}>=0\). The reader is referred to [43,44,45, 47, 49, 81,82,83,84,85,86] for the most recent studies on the DM.

The present work is arranged as follows: in Sect. 2, the most general and renormalizable extension of the SM is presented by adding two scalar sectors \(S_{1}\) and \(S_{2}\) to the usual SM potential.Footnote 1 Assigning a non-zero VeV to \(S_{1}\), the SFOEWPTH can occur in the model. Imposing a \(Z_{2}\) symmetry on \(S_{2}\) makes it a viable candidate for the DM. Also, constraints on the parameter space are discussed. The EWPT, GW and DM are respectively investigated in Sects. 3, 4 and 5. Finally, some conclusions are presented in Sect. 6.

2 The model

The tree-level potential of the model is given by

$$\begin{aligned} V= & {} -\, m^{2} H^{\dagger }H + \lambda (H^{\dagger }H)^2 + \kappa _{0} S_{1} \nonumber \\&+ 2 (\kappa _{1} S_{1}+\kappa _{2} S_{1}^{2}+\kappa _{3} S_{2}^{2}) H^{\dagger }H \nonumber \\&+ \frac{1}{2} m_{1}^{2} S_{1}^{2} + \frac{\lambda _{1}}{4} S_{1}^{4} + \kappa _{4} S_{1}^{3} \nonumber \\&+ \frac{1}{2} m_{2}^{2} S_{2}^{2} + \frac{\lambda _{2}}{4} S_{2}^{4} + \kappa _{5} S_{1} S_{2}^{2}, \end{aligned}$$
(2.1)

The potential 2.1 is the usual SM potential with two extra gauge singlet scalars and interaction terms which provide Higgs portal between the new scalars and the usual SM particles. H stands for the complex Higgs doublet, \(H=\begin{pmatrix}\chi _{1}+i\chi _{2} \\ \frac{1}{\sqrt{2}}(h+i\chi _{3}) \\ \end{pmatrix}\). \(S_{2}\) stands for the DM imposing \(S_{2}\rightarrow -S_{2}\). Acquiring a non-zero VeV, \(S_{1}\) improves the strength of EWPT. The linear term of \(S_{1}\) can be neglected by a shift in the potential. The \(Z_{2}\) symmetry forbids the existence of linear and cubic terms for \(S_{2}\), so Eq. (2.1) is the most general renormalizable potential which could be made by adding two new scalars. In the unitary gauge at zero temperature, the theoretical fields can be reparameterized in terms of the physical fields,

$$\begin{aligned} H=\begin{pmatrix}0 \\ \frac{1}{\sqrt{2}}(h+v) \\ \end{pmatrix}, S_{1}=s_{1}+\chi ,\quad S_{2}=s_{2}, \end{aligned}$$
(2.2)

where \(v=246.22\,(\mathrm{{GeV}})\) and \(\chi \) are the Higgs and \(S_{1}\) VeV, respectively. Without loss of generality, one can write

$$\begin{aligned} V= & {} -\frac{1}{2} m^{2} h^{2} + \frac{\lambda }{4} h^{4} + (\kappa _{1} s_{1}+\kappa _{2} s_{1}^{2}+\kappa _{3} s_{2}^{2}) h^{2}\nonumber \\&+ \frac{1}{2} m_{1}^{2} s_{1}^{2} + \frac{\lambda _{1}}{4} s_{1}^{4} + \kappa _{4} s_{1}^{3} \nonumber \\&+ \frac{1}{2} m_{2}^{2} s_{2}^{2} + \frac{\lambda _{2}}{4} s_{2}^{4} + \kappa _{5} s_{1} s_{2}^{2}. \end{aligned}$$
(2.3)

In order to have a stable potential, it is required that [22, 87]

$$\begin{aligned} \lambda>0,\quad \lambda _{1}>0,\quad \lambda _{2}>0,\quad \kappa _{2}>-\frac{\sqrt{\lambda \lambda _{1}}}{2},\quad \kappa _{3}>-\frac{\sqrt{\lambda \lambda _{2}}}{2}. \end{aligned}$$
(2.4)

The tadpole equations at \((v,\chi ,0)\) read

$$\begin{aligned} \begin{aligned} m^{2}&=\lambda v^{2} + 2 (\kappa _{1} \chi + \kappa _{2} \chi ^{2}), \\ m_{1}^{2}&=-\lambda _{1} \chi ^{2} - 3 \kappa _{4} \chi - 2 \kappa _{2} v^{2} - \frac{\kappa _{1} v^{2}}{\chi }. \end{aligned} \end{aligned}$$
(2.5)

From the diagonalization of squared-mass matrix and the tadpole equations, one can get

$$\begin{aligned} \begin{aligned} \lambda&= \frac{M_{1}^{2}sin^{2}(\theta )+M_{H}^{2}cos^{2}(\theta )}{2 v^{2}}, \\ \kappa _{2}&= \frac{(M_{H}^{2}-M_{1}^{2})sin(2\theta )}{8v\chi }-\frac{\kappa _{1}}{2\chi }, \\ \lambda _{1}&= \frac{1}{2\chi ^{2}}\left( M_{1}^{2}cos^{2}(\theta )+M_{H}^{2}sin^{2}(\theta )+\frac{\kappa _{1}v^{2}}{\chi }-3 \chi \kappa _{4}\right) , \\ m_{2}^{2}&= M_{2}^{2}-2\kappa _{3}v^{2}-2\kappa _{5}\chi , \end{aligned} \end{aligned}$$
(2.6)

where \(M_{H}=126 \,(\mathrm{{GeV}})\), \(M_{1}\), \(M_{2}\) and \(\theta \) are the Higgs mass,Footnote 2 the physical mass of \(S_{1}\), the physical mass of \(S_{2}\) (the DM mass) and the mixing angle, respectively. In Ref. [89], by performing a global fit to the Higgs data from both \(\mathbf {ATLAS}\) and \(\mathbf {CMS}\), the constraint on the mixing angle was given \(|\theta |\le 32.86^{\circ }\) at \(95\%\) confidence level (CL). In Ref. [90], by performing a universal Higgs fit, the upper limit on the mixing angle was given \(|\theta |\le 30.14^{\circ }\) at \(95\%\) CL. In the present work, a Monte Carlo scan is performed over the parameter space with

$$\begin{aligned} \begin{aligned}&5\;\mathrm{{GeV}} \le M_{1} \le 750 \;\mathrm{{GeV}},\quad 5\;\mathrm{{GeV}}\le M_{2} \le 750 \;\mathrm{{GeV}},\\&\quad -23^{\circ }\le \theta \le 23^{\circ }, \\&\quad -80\;\mathrm{{GeV}}\le \kappa _{1}\le 80\;\mathrm{{GeV}},\quad 0.0001\le \kappa _{3}\le 0.1,\\&\quad -80\;\mathrm{{GeV}}\le \kappa _{4}\le 80\;\mathrm{{GeV}}, \\&\quad -80\;\mathrm{{GeV}}\le \kappa _{5}\le 80\;\mathrm{{GeV}},\quad 30 \;\mathrm{{GeV}}\le \chi \le 120 \;\mathrm{{GeV}},\\&\quad 0\le \lambda _{2}\le 4. \end{aligned} \end{aligned}$$
(2.7)

3 Electroweak phase transition

To investigate the EWPT in a model, one needs the finite temperature effective potential given by

$$\begin{aligned} V_{eff}=V_{tree{\text {-}}level}+V_{1-loop}^{T=0}+V_{1-loop}^{T\ne 0}, \end{aligned}$$
(3.1)

where \(V_{tree{\text {-}}level}\), \(V_{1-loop}^{T=0}\) and \(V_{1-loop}^{T\ne 0}\) are the tree-level potential (2.3), the one-loop corrected potential at zero temperature (the so-called Coleman-Weinberg potential) and the one-loop finite temperature corrections, respectively. The last two read

$$\begin{aligned} \begin{aligned}&V_{1-loop}^{T=0}=\pm \frac{1}{64\pi ^2}\sum _{i=h,s_{1},s_{2},W,Z,t} n_i m_i^4\left[ \log \frac{m_i^2}{Q^2}-C_i\right] ,\\&V_{1-loop}^{T\ne 0}=\frac{T^4}{2\pi ^2} \sum _{i=h,s_{1},s_{2},W,Z,t} n_i J_\pm \left[ \frac{m_{i}^{2}}{T^{2}}\right] ,\\&J_\pm \left( \frac{m_{i}^{2}}{T^{2}}\right) = \pm \int _0^\infty dy\, y^2 \log \left( 1\mp e^{-\sqrt{y^2+\frac{m_{i}^{2}}{T^{2}}}}\right) , \end{aligned} \end{aligned}$$
(3.2)

where \(n_{i}\), \(m_{i}\), Q and \(C_{i}\) denote the degrees of freedom, the field-dependent masses, the renormalization scale and the numerical constants, respectively. The degrees of freedom and the numerical constants are respectively given by \((n_{h,s_{1},s_{2}},n_{W},n_{Z},n_{t})=(1,6,3,12)\) and \((C_{W,Z},C_{h,s_{1},s_{2},t})=(5/6,3/2)\). The upper (lower) sign is for bosons (fermions). Assuming the longitudinal gauge bosons polarizations are screened by plasma, thermal masses just contribute to the scalars, so Daisy corrections become small and can be neglected. There are three possibilities to deal with the renormalization scale Q. First one is to add some counter terms to the effective potential (3.1) to make it independent of Q without shifting VeV at zero temperature [91, 92]. Second one is to set Q at a proper scale, like \(Q=160\,(\mathrm{{GeV}})\) the running value of the top mass, \(Q=246.22\,(\mathrm{{GeV}})\) EW scale and \(Q=1(\mathrm{{TeV}})\) for supersymmetry purposes. Third one is to take Q as a free parameter to avoid shifting VeV at zero temperature. Here, the last one is considered.

The main idea of the EWPT is that the early universe, which from particle physics point of view may be described by potential (3.1), is in a high phaseFootnote 3 with \(VeV=(<h>,<s1>,<s2>)^{high}\) at high temperatures. Cooling down the universe, a new phase appears with \(VeV=(<h>,<s1>,<s2>)^{low}\). As the universe cools down, the two phases become degenerate at the critical temperature \(T_{c}\). Since the strength of the EWPT is governed by \(\xi =\frac{v_{c}}{T_{c}}\), all that needs to be done is to calculate \(v_{c}\) and \(T_{c}\) from the following conditions:

$$\begin{aligned} \begin{aligned}&\left. \frac{\partial V_{eff}}{\partial h}\right| _{(<h>,<s1>,<s2>)^{high},T=T_{c}} = 0, \\&\left. \frac{\partial V_{eff}}{\partial h}\right| _{(<h>,<s1>,<s2>)^{low},T=T_{c}} = 0,\\&\left. \frac{\partial V_{eff}}{\partial s_{1}}\right| _{(<h>,<s1>,<s2>)^{high},T=T_{c}} = 0, \\&\left. \frac{\partial V_{eff}}{\partial s_{1}}\right| _{(<h>,<s1>,<s2>)^{low},T=T_{c}} = 0,\\&\left. \frac{\partial V_{eff}}{\partial s_{2}}\right| _{(<h>,<s1>,<s2>)^{high},T=T_{c}} = 0, \\&\left. \frac{\partial V_{eff}}{\partial s_{2}}\right| _{(<h>,<s1>,<s2>)^{low},T=T_{c}} = 0,\\&V_{eff}\Big |_{(<h>,<s1>,<s2>)^{high},T=T_{c}} \\&\quad = V_{eff}\Big |_{(<h>,<s1>,<s2>)^{low},T=T_{c}}. \end{aligned} \end{aligned}$$
(3.3)

The last condition guarantees degeneracy and the others guarantee existence of the high and low vacua. There is no analytical solution for the problem, so the calculations are implemented with the CosmoTransitions package [93]. The benchmark points and the corresponding results are presented in Tables 1 and 2, respectively. Here, the exact calculations are performed by CosmoTransitions to get the effective potential, compared to Ref. [87] which used the high temperature expansion. Though, the results of Ref. [87] should be improved for the high temperature expansion case. An extension of the SM with two new scalars was recently studied in Ref. [94], but there are some differences between it and the present work. First, the high temperature expansion was used in [94]. Second, the cubic term \(S_{1}^{3}\), which plays a crucial role in the EWPT as a barrier at tree-level, was not considered in  [94].

Table 1 Benchmark points which provide the SFOEWPT
Full size table

4 Gravitational waves

The SFOEWPT may justify not only the BAU but also the GW signal produced by the EWPT. Actually, the EWPT occurs at a temperature lower than \(T_{c}\), in which the first broken phase bubbles nucleate in the symmetric phase plasma of the early universe. The transition probability is given by \(\Gamma (T)=\Gamma _{0}(T)e^{-S(T)}\) where \(\Gamma _{0}(T)\) is of order \(\mathcal {O}(T^{4})\) and S is the 4-dimensional action of the critical bubbles. For temperatures sufficiently greater than zero, it can be assumed \(S=\frac{S_{3}}{T}\) where the 3-dimensional action is given by

$$\begin{aligned} S_{3} = 4\pi \int {dr r^{2} \left[ \frac{1}{2}\left( \partial _{r} \vec {\phi }\right) ^2 + V_{eff}\right] }. \end{aligned}$$
(4.1)

Here, \(\vec {\phi }=(h,s1,s2)\). The critical bubble profiles, which minimize the action (4.1), can be calculated from the equation of motions. The temperature for a particular configuration, which gives the nucleation probability of order \(\mathcal {O}(1)\), is the nucleation temperature \(T_{n}\).

Table 2 The values of the VeV of the high and the low phases, \(T_{c}\) and the strength of the SFOEWPT
Full size table

The GW may be produced by the collision of the bubbles at some temperature \(T_{*}\), it is usually assumed \(T_{*}=T_{n}\). Supposing that the friction force is not enough to prevent the bubbles from running away, the GW signal is given by

$$\begin{aligned} \Omega _{GW} h^{2}\simeq \Omega _{col} h^{2}+\Omega _{sw} h^{2}+\Omega _{turb} h^{2}. \end{aligned}$$
(4.2)

As seen, the GW signal is given by the sum of bubbles collision, sound wave and turbulence in the plasma which respectively read [53, 55, 56, 59, 64, 73, 95, 96]

$$\begin{aligned} \Omega _{col} h^{2}= & {} 1.67 \times 10^{-5} \left( \frac{\beta }{H}\right) ^{-2} \frac{0.11\, v_{b}^{3}}{0.42+v_{b}^{2}} \nonumber \\&\quad \times \left( \frac{\kappa \, \alpha }{1+\alpha }\right) ^{2} \left( \frac{g_{*}}{100}\right) ^{-\frac{1}{3}}\frac{3.8\, \left( \frac{f}{f_{col}}\right) ^{2.8}}{1+2.8\, \left( \frac{f}{f_{col}}\right) ^{3.8}},\nonumber \\ \Omega _{sw} h^{2}= & {} 2.65 \times 10^{-6} \left( \frac{\beta }{H}\right) ^{-1} v_{b} \left( \frac{\kappa _{v}\, \alpha }{1+\alpha }\right) ^{2} \nonumber \\&\quad \times \left( \frac{g_{*}}{100}\right) ^{-\frac{1}{3}} \left( \frac{f}{f_{sw}}\right) ^{3} \left( \frac{7}{4+3\left( \frac{f}{f_{sw}}\right) ^{2}}\right) ^{\frac{7}{2}},\nonumber \\ \Omega _{turb} h^{2}= & {} 3.35 \times 10^{-4} \left( \frac{\beta }{H}\right) ^{-1} v_{b} \left( \frac{\epsilon \, \kappa _{v}\, \alpha }{1+\alpha }\right) ^{\frac{3}{2}} \nonumber \\&\qquad \times \left( \frac{g_{*}}{100}\right) ^{-\frac{1}{3}} \frac{\left( \frac{f}{f_{turb}}\right) ^{3} \left( 1+\frac{f}{f_{turb}}\right) ^{-\frac{11}{3}} }{1+\frac{8\pi f}{h_{*}}}, \end{aligned}$$
(4.3)

with

$$\begin{aligned} \begin{aligned}&f_{col} = 16.5\times 10^{-6}\, \mathrm {Hz}\, \left( \frac{0.62}{v_{b}^{2}-0.1\,v_{b}+1.8}\right) \\&\quad \qquad \times \left( \frac{\beta }{H}\right) \left( \frac{T_{n}}{100\, \mathrm {GeV}}\right) \left( \frac{g_{*}}{100}\right) ^{\frac{1}{6}},\\&f_{sw} = 1.9\times 10^{-5}\, \mathrm {Hz}\, \left( \frac{1}{v_{b}}\right) \left( \frac{\beta }{H}\right) \left( \frac{T_{n}}{100\, \mathrm {GeV}}\right) \left( \frac{g_{*}}{100}\right) ^{\frac{1}{6}},\\&f_{turb} = 2.7\times 10^{-5} \mathrm {Hz} \left( \frac{1}{v_{b}}\right) \left( \frac{\beta }{H}\right) \left( \frac{T_{n}}{100\, \mathrm {GeV}}\right) \left( \frac{g_{*}}{100}\right) ^{\frac{1}{6}},\\&h_{*} = 16.5\times 10^{-6}\, \mathrm {Hz}\, \left( \frac{T_{n}}{100\, \mathrm {GeV}}\right) \left( \frac{g_{*}}{100}\right) ^{\frac{1}{6}},\\&\kappa =1-\frac{\alpha _{\infty }}{\alpha },\\&\kappa _{v}=\frac{\alpha _{\infty }}{\alpha }\left( \frac{\alpha _{\infty }}{0.73+0.083\sqrt{\alpha _{\infty }}+\alpha _{\infty }}\right) ,\\&\alpha _{\infty }=\frac{30}{24\, \pi ^2 g_{*}} \left( \frac{v_{n}}{T_{n}}\right) ^{2} \left( 6\left( \frac{M_{W}}{v}\right) ^{2} \right. \\&\quad \qquad \left. + 3\left( \frac{M_{Z}}{v}\right) ^{2} + 6\left( \frac{M_{top}}{v}\right) ^{2} \right) . \end{aligned} \end{aligned}$$
(4.4)
Table 3 The values of the VeV of the high and the low phases, \(T_{n}\), \(\alpha \) and \(\beta /H\)
Full size table

\(v_{n}\) and \(g_{*}\) are the Higgs VeV and the number of the relativistic degrees of freedom at \(T_{n}\), respectively. Here, \(\epsilon =0.1\), and \(g_{*}\) is read from the MicrOMEGAs package [97, 98]. Still, there are three important parameters which should be defined. One of them is the bubble wall velocity, since assumed that the bubbles run away, given by \(v_{b} \simeq 1\). The two others, \(\alpha \) and \(\beta \), are given as follows

$$\begin{aligned} \begin{aligned}&\alpha = \left. \frac{\rho _{vac}}{\rho ^{*}}\right| _{T_{n}},\\&\beta = \left. \left[ H\, T\, \frac{d}{dT}\left( \frac{S_{3}}{T}\right) \right] \right| _{T_{n}}, \end{aligned} \end{aligned}$$
(4.5)

where \(\rho _{vac}=\left( V_{eff}^{high}-T dV_{eff}^{high}/dT \right) -\left( V_{eff}^{low}-T dV_{eff}^{low}/dT \right) \), \(\rho ^{*}=g_{*}\pi ^{2}T_{n}^{4}/30\) and \(H_{n}\) are the latent heat (vacuum energy density) released by the EWPT, the background energy density of the plasma and Hubble parameter at \(T_{n}\), respectively. Using the CosmoTransitions package [93], the parameters \(\alpha \), \(\beta /H\), \(v_{n}\) and \(T_{n}\) are calculated and presented in Table 3. In Fig. 1, the GW signals are plotted versus frequency for the benchmark points of Table 1. To check if the GW signals for the benchmark points 1 fall within the sensitivity of GW detectors, the sensitivity curves of \(\mathbf {eLISA}\), \(\mathbf {ALIA}\), \(\mathbf {DECIGO}\) and \(\mathbf {BBO}\) detectorsFootnote 4 are also plotted in the Fig. 1. As seen from the Fig. 1, the dashed blue line corresponding to the GW signal for the BM7 may be detected by \(\mathbf {N2A1M5L6}\) and \(\mathbf {N2A5M5L6}\) configurations of \(\mathbf {eLISA}\) and \(\mathbf {BBO}\) detectors. The dashed red and yellow lines corresponding, respectively, to the GW signal for the BM4 and the BM6 may be detected by \(\mathbf {N2A5M5L6}\) configuration of \(\mathbf {eLISA}\) and \(\mathbf {BBO}\) detectors. The dashed orange line corresponding to the BM2 may be detected by \(\mathbf {DECIGO}\) and \(\mathbf {BBO}\) detectors. The dashed green, cyan and purple lines corresponding, respectively, to the GW signal for BM1, BM5 and BM8 cannot be detected by the mentioned detectors. The GW signal for the BM3 is not big enough to be shown at the scale of the Fig. 1.

Fig. 1
figure 1

The dashed blue, red, yellow, orange, green, cyan and purple lines represent the GW signal for BM7, BM4, BM6, BM2, BM1, BM5 and BM8, respectively. The solid black lines represent the sensitivity curves and are labeled by the name of the detectors. For eLISA, the sensitivity curves are labeled by the name of the configuration

Full size image
Fig. 2
figure 2

Phase transition properties of BM2: the subfigure a presents \(S_{3}/T\) versus temperature, the dashed horizontal red line shows \(S_{3}/T=140\) where nucleation occurs. The subfigure b presents the tunneling profile as a function of radius. The subfigure c presents the norms of high (green line) and low (blue line) phases as functions of temperature, the dashed vertical red line shows the nucleation temperature. The subfigure d presents the contour levels of the potential at the nucleation temperature \(T_{n}=41.61 \,(\mathrm{{GeV}})\), the dashed black line shows the tunneling path

Full size image

According to the Tables 2 and 3, it seems that BM2 is a special point. The value of \(\beta /H\) is large at this point, while, the nucleation temperature is not very close to the critical temperature. At the same time, \(T_{n}\) is low and \(\alpha \) is large.Footnote 5 To clarify the situation of BM2, the phase transition properties of BM2 are shown in Fig. 2. As seen from the Fig. 2a, the slope of \(S_{3}/T\) increases around \(T_{n}\) which indicates the parameter \(\beta /H\) is large. The physics of this situation can be described by the tunneling profile, the norm of phases as a function of temperature, and the contour levels of the potential with the tunneling path. As seen from Fig. 2b, d, the center of bubble is far away from the stable vacuum. Also from Fig. 2c, it is seen that the transition occurs at a temperature where the unstable vacuum is close to disappearance. The values of potential at high and low phases are, respectively, \(V^{high}_{eff}=-91583128.19 (\mathrm{{GeV}}^{4})\) and \(V^{low}_{eff}=-101840540.47 (\mathrm{{GeV}}^{4})\), which give the pressure difference \(\Delta p = 10257412.28 (\mathrm{{GeV}}^{4})\). The barrier location is at \((h, s_{1}, s_{2})=(19.84, 212.19, 0)\) with \(V_{eff}=-91582001.71 (\mathrm{{GeV}}^{4})\) which gives the barrier height \(\Delta V_{barrier \, height}=1126.48 (\mathrm{{GeV}}^{4})\). Clearly, the barrier height is very small, \(\Delta V_{barrier \, height}/\Delta p = 0.0001\). Due to the reasons given above, the bubbles are extremely thick walled. Since the barrier height is very small, the transition duration is very short, accordingly, the parameter \(\beta /H\) is quite large. This extremely thick walled case is similar to the second-order phase transition in which \(\beta /H \rightarrow \infty \) and there is no barrier. Moreover, there is another interesting note for BM2. Due to the cubic term \(s_{1}^{3}\), it is expected that the model has a sizable barrier at tree-level like the supercooled scenario discussed in [78], but this is not the case for BM2. At this point, the model mimics the behavior of supercooled phase transitions with the supercooling parameter \((T_{c}-T_{n})/T_{c}=0.31\), though, the transition is short-lived.

5 Dark matter

As mentioned prior, imposing the \(Z_{2}\) symmetry on \(S_{2}\) makes it a viable candidate for the DM. Considering the freeze-out formalism, the DM relic density abundance can be calculated by solving the Boltzmann equation,

$$\begin{aligned} \frac{dn}{dt}=-3 H n-<\sigma v> (n^{2}-n_{eq}^{2}), \end{aligned}$$
(5.1)

where n, H, \(<\sigma v>\) are the number of the DM particles, the Hubble parameter and the thermally-averaged cross section for the DM annihilation, respectively. It is customary to rewrite the Boltzmann equation in terms of \(Y=n/s\), where s is the total entropy density of the universe, the result is [100]

$$\begin{aligned} \frac{dY}{dx}= & {} -\left( \frac{45\,G\,g_{*}}{\pi }\right) ^{-\frac{1}{2}}\frac{M\,h_{*}}{x^{2}}\left( 1+\frac{1}{3}\frac{T}{h_{*}}\frac{dh_{*}}{dT}\right) \nonumber \\< & {} \sigma v>(Y^{2}-Y_{eq}^{2}), \end{aligned}$$
(5.2)

where \(x=M/T\), M is the DM mass. \(h_{*}\) is the effective degree of freedom for the entropy densities. The DM relic density abundance reads,

$$\begin{aligned} \Omega _{DM} h^{2}\simeq (2.79\pm 0.05) \times 10^{8} \left( \frac{M}{\mathrm{{GeV}}}\right) Y(0). \end{aligned}$$
(5.3)

It is assumed the usual SM particles only interact with Higgs in this model, so the annihilation channels for DM via the Higgs portal s-channel are \(s_{2}s_{2}\rightarrow W^{+}W^{-},ZZ,f\bar{f}\). Also, there exists \(s_{2}s_{2}\rightarrow \phi _{i}\phi _{j}\) (with \(\phi _{i(j)}=h,s_{1}\) and \(i(j)=1,2\)) via s, t and u channels and four-point interactions.

The parameter space is constrained by the direct detection DM searches. To do this, one needs to calculate the spin-independent cross section for DM-nucleon scattering,Footnote 6 and compares the result with the \(\mathbf {XENON1T}\) \(\mathbf {2018}\) experiment data [101]. The spin-independent cross section is given by

$$\begin{aligned} \sigma _{SI}=\frac{4 M_{s_{2}}^{2} M_{N}^{2}}{\pi (M_{s{2}}+M_{N})^{2}}\Big | \mathcal {M}_{s_{2}-N} \Big |^{2}, \end{aligned}$$
(5.4)

where \(M_{s_{2}}\), \(M_{N}\) and \(\mathcal {M}_{s_{2}-N}\) are the DM mass, the nucleus mass and the scattering amplitude at low energy limit, respectively. \(\mathcal {M}_{s_{2}-N}\) is related to \(\mathcal {M}_{s_{2}-quark}\), so, calculating effective Lagrangian coefficients and nucleon form factors, \(\mathcal {M}_{s_{2}-N}\) can be obtained from \(\mathcal {M}_{s_{2}-quark}\). Here, the model is implemented in SARAH [102,103,104], the model spectrum is obtained by SPheno [105, 106] and the DM properties are studied by MicrOMEGAs [97, 98]. The results are presented in Table 4.

Table 4 The values of the DM relic abundance and the spin-independent cross sections
Full size table

As seen in Table 4, the relic density of all benchmark points is compatible with \(\mathbf {Planck}\) \(\mathbf {2018}\) data which reports \(\Omega _{c}h^{2}=0.120\pm 0.001\).Footnote 7 Moreover, the results fit with \(\mathbf {XENON1T}\) \(\mathbf {2018}\) experiment which gives an upper limit, less than \(\mathbf {LUX}\) \(\mathbf {2017}\) [107] and \(\mathbf {PandaX}\)-\(\mathbf {II}\) \(\mathbf {2017}\) [108] reports, on the DM-nucleon spin-independent elastic scattering cross section. In the DM study, there are two differences with Ref. [94]. The first is the \(s_{1}^{3}\) interaction which gives a significant contribution to the DM annihilation through \(s_{1}\) s-channel, and consequently to relic density. The second is the spin-independent cross section which was taken to be zero in Ref. [94], but the more realistic case like here is to have a non-zero DM-nucleon cross section, if weakly interacting massive particles (WIMP) constitute the DM. This is the main idea behind \(\mathbf {LUX}\), \(\mathbf {PandaX}\)-\(\mathbf {II}\) and \(\mathbf {XENON1T}\) experiments.

6 Conclusions

The main goal of this work has been to investigate the EWPT, the GW and the DM issues in an extension of the SM by adding two scalar degrees of freedom. To reach the goal, it has been assumed that one of the new scalars has a non-zero VeV to assist the phase transition and the other has no VeV to be a viable DM candidate. It has been seen if one takes the most general renormalizable form of the potential, the model can represent all the signals together. As seen from Tables 2 and 3, the model can have phase transitions from strong (\(\xi \sim 1\)) to very strong (\(\xi \sim 4\)). From Fig. 1, the model presents the GW signals from the frequency range of \(10^{-5}\)(Hz) to 10(Hz) which are detectable by \(\mathbf {eLISA}\), \(\mathbf {BBO}\) and \(\mathbf {DECIGO}\). From Table 4, the model provides the DM signals which are in agreement with the \(\mathbf {Planck}\) \(\mathbf {2018}\) data and the \(\mathbf {XENON1T}\) \(\mathbf {2018}\) experiment. It is seen that the DM candidate may be quite massive with a mass greater than 100(GeV) which belongs to the extremely cold DM; although, since the model has a rich parameter space, the lighter DMs might be found by performing a Monte Carlo simulation via a computer cluster. With all of these, it can be concluded that the SFOEWPT, the GW and the DM signals can successfully be described by the present model as an extension of the SM with two additional real gauge singlet scalars. As a final note, it has been assumed that the GW production from bubble collisions follows from thin-wall and envelope approximations which is usual in the literature. In this assumption, only uncollided parts of the bubbles are taken into account as the GW sources. Recently, it has been shown that the GW production from bubble collisions is analytically solvable [109, 110]. Also, the possibility of using GWs and collider experiments to constrain the EWPT has been discussed in [111]. It is left for future work to study the GW signals of the present model using these recent studies.