Gravitational wave signatures from an extended inert doublet dark matter model

TL;DR

This paper explores a dark matter model involving an inert doublet and a singlet scalar, analyzing the resulting gravitational wave signatures from early Universe phase transitions and their detectability with future observatories.

Contribution

It introduces an extended inert doublet dark matter model with a singlet scalar and studies its gravitational wave signatures from electroweak phase transitions.

Findings

First-order electroweak phase transition occurs in the model.

Predicted gravitational wave signals are within the sensitivity of future detectors.

The model provides a viable dark matter candidate compatible with current constraints.

Abstract

We consider a particle dark matter model by extending the scalar sector of the Standard Model by an additional SU(2) scalar doublet which is made "inert" (and stable) by imposing a discrete $Z_2$ symmetry under which the additional scalar doublet is odd (and the SM is even) and it does not develop any vacuum expectation value (VEV). The lightest inert particle (LIP) of this inert doublet model (IDM) can be a viable candidate for Dark Matter. The IDM model is further extended by an additional singlet scalar which is also even under $Z_2$ and develop a VEV on spontaneous symmetry breaking. This additional scalar singlet mixes with SM Higgs and on diagonalisation of the mass matrix two CP even scalar eigenstates are obtained one of which is attributed to the physical Higgs (with mass 125.09 GeV). The LIP is the dark matter candidate in the extended model. For such a particle dark matter model we explore the first-order electroweak phase transition and the consequent production of Gravitational Waves (GW) at that epoch of the early Universe and calculate the intensities and frequencies for such waves. We then investigate the detection possibilities of such GWs at the future spaceborne primordial GW detectors such as eLISA, BBO, ALIA, DECIGO, U-DECIGO and ground-based detector aLIGO.

Tables

Table 1. Table 1 : The chosen four benchmarks points (BPs, BP1-4) to explore the GW production from an extended IDM with an additional real singlet scalar. The relic density and scattering cross-section values for each of the BPs are also mentioned in this Table.

BP	$m_{H_0}$	$m_{h_2}$	$v_s$	$\sin \theta$	${\rho }_1$	${\rho }_3$	${\lambda }_L$	${\lambda }_s$	${\lambda }_2$	${\Omega }_{\text{DM}}{h}^2$	${\sigma }_{\mathrm{SI}}$
	in GeV	in GeV	in GeV		in GeV	in GeV					${\mathrm{cm}}^2$
1	30	100	300	0.01	-3	0.01	0.001	0.0012	0.2	0.1220	9.41$\times {10}^{-48}$
2	68	150	400	0.06	-7	0.2	0.002	0.033	0.031	0.1208	3.69$\times {10}^{-48}$
3	62	250	400	0.01	-5	0.5	0.0003	0.0033	0.1	0.1206	1.70$\times {10}^{-49}$
4	76	200	500	0.03	-1	0.1	0.0016	0.0033	0.01	0.1195	3.56$\times {10}^{-48}$

Table 2. Table 2 : The values for the parameters used to calculate the GW intensity for each of the chosen BPs. See text for details.

BP	$v_n$	$T_c$	$T_n$	$\alpha$	$\frac{\beta }{H}$
	in GeV	in GeV	in GeV
1	226.89	135.68	119.86	0.24	317.86
2	191.03	146.89	132.14	0.25	402.89
3	14.22	130.78	116.99	0.23	692.01
4	209.95	170.92	158.24	0.19	783.65

Equations

$$V=m_1^2{\Phi }_H^{†}{\Phi }_H+m_2^2{\Phi }_I^{†}{\Phi }_I+\frac{1}{2}{m}_s^2{S}^2+{\lambda }_1{\left({\Phi }_H^{†}{\Phi }_H\right)}^2+{\lambda }_2{\left({\Phi }_I^{†}{\Phi }_I\right)}^2+{\lambda }_3\left({\Phi }_H^{†}{\Phi }_H\right)\left({\Phi }_I^{†}{\Phi }_I\right)$$

$$+{\lambda }_4\left({\Phi }_I^{†}{\Phi }_H\right)\left({\Phi }_H^{†}{\Phi }_I\right)+\frac{{\lambda }_5}{2}\left[{\left({\Phi }_I^{†}{\Phi }_H\right)}^2+{\left({\Phi }_H^{†}{\Phi }_I\right)}^2\right]+{\rho }_1\left({\Phi }_H^{†}{\Phi }_H\right)S$$

$$+{\rho }_1^{\prime }\left({\Phi }_I^{†}{\Phi }_I\right)S+{\rho }_2\left({\Phi }_H^{†}{\Phi }_H\right)S^2+{\rho }_2^{\prime }\left({\Phi }_I^{†}{\Phi }_I\right)S^2+\frac{{\rho }_3}{3}{S}^3+\frac{{\rho }_4}{4}{S}^4.$$

$$\phi_{H}=\left(\begin{array}[]{c}G^{+}\\ \dfrac{1}{\sqrt{2}}\left(v+h+iG_{0}\right)\\ \end{array}\right),\hskip 2.84526pt\phi_{I}=\left(\begin{array}[]{c}H^{+}\\ \dfrac{1}{\sqrt{2}}\left(H_{0}+iA_{0}\right)\\ \end{array}\right),\hskip 2.84526ptS=v_{s}+s,$$

$$\dfrac{\partial V}{\partial h}\Big{|}_{h=s=H^{+}=H_{0}=A_{0}=0}=\dfrac{\partial V}{\partial s}\Big{|}_{h=s=H^{+}=H_{0}=A_{0}=0}=0,$$

$$m_1^2+{\lambda }_1{v}^2+{\rho }_1{v}_s+{\rho }_2{v}_s^2=0,$$

$$m_s^2+{\rho }_3{v}_s+{\rho }_4{v}_s^2+\frac{{\rho }_1{v}^2}{2v_s}+{\rho }_2{v}^2=0.$$

$${\mu }_h^2=m_1^2+3{\lambda }_1{v}^2+{\rho }_1{v}_s+{\rho }_2{v}_s^2=2{\lambda }_1{v}^2,$$

$${\mu }_s^2=m_s^2+{\rho }_2{v}^2+2{\rho }_3{v}_s+3{\rho }_4{v}_s^2={\rho }_3{v}_s+2{\rho }_4{v}_s^2-\frac{{\rho }_1{v}^2}{2v_s},$$

$$m_{H^{\pm }}^2=m_2^2+\frac{{\lambda }_3{v}^2}{2}+{\rho }_1^{\prime }{v}_s+{\rho }_2^{\prime }{v}_s^2,$$

$$m_{H_0}^2=m_2^2+\left({\lambda }_3+{\lambda }_4+{\lambda }_5\right)\frac{v^2}{2}+{\rho }_1^{\prime }{v}_s+{\rho }_2^{\prime }{v}_s^2,$$

$$m_{A_0}^2=m_2^2+\left({\lambda }_3+{\lambda }_4-{\lambda }_5\right)\frac{v^2}{2}+{\rho }_1^{\prime }{v}_s+{\rho }_2^{\prime }{v}_s^2,$$

$${\mu }_{hs}^2={\rho }_{1}v+2{\rho }_2{v}_{s}v.$$

$$\left(\begin{array}[]{c}h_{1}\\ h_{2}\\ \end{array}\right)=U\left(\begin{array}[]{c}h\\ s\\ \end{array}\right)=\left(\begin{array}[]{cc}\cos\theta&-\sin\theta\\ \sin\theta&\cos\theta\\ \end{array}\right)\left(\begin{array}[]{c}h\\ s\\ \end{array}\right),$$

$$h_1=h\cos \theta -s\sin \theta ,h_2=h\sin \theta +s\cos \theta ,$$

$$\tan \theta =\frac{y}{1+\sqrt{1+y^2}},\text{where}y=\frac{2{\mu }_{hs}^2}{{\mu }_h^2-{\mu }_s^2}.$$

$$m_{h_1,h_2}^2=\frac{\left({\mu }_h^2+{\mu }_s^2\right)}{2}\pm \frac{\left({\mu }_h^2-{\mu }_s^2\right)}{2}\sqrt{1+y^2},$$

$${\lambda }_1=\frac{m_{h_2}^2{\sin }^2\theta +m_{h_1}^2{\cos }^2\theta }{2v^2},$$

$${\rho }_4=\frac{1}{2v_s^2}\left(m_{h_2}^2{\cos }^2\theta +m_{h_1}^2{\sin }^2\theta +\frac{v^2{\rho }_1}{2v_s}-v_s{\rho }_3\right),$$

$${\rho }_2=\frac{m_{h_1}^2-m_{h_2}^2}{4v{v}_s}\sin 2\theta -\frac{{\rho }_1}{2v_s},$$

$$m_1^2=-{\lambda }_1{v}^2-{\rho }_1{v}_s-{\rho }_2{v}_s^2,$$

$$m_2^2=m_{H_0}^2-\frac{v^2}{2}\left({\lambda }_3+{\lambda }_4+{\lambda }_5\right)-{\rho }_1^{\prime }{v}_s-{\rho }_2^{\prime }{v}_s^2,$$

$$m_s^2=-{\rho }_3{v}_s-{\rho }_4{v}_s^2-\frac{{\rho }_1{v}^2}{2v_s}-{\rho }_2{v}^2.$$

$${\lambda }_1,{\lambda }_2,{\rho }_4>0,{\lambda }_3+2\sqrt{{\lambda }_1{\lambda }_2}>0,$$

$${\lambda }_3+{\lambda }_4-|{\lambda }_5|+2\sqrt{{\lambda }_1{\lambda }_2}>0,{\rho }_2+\sqrt{{\lambda }_1{\rho }_4}>0,$$

$${\rho }_2^{\prime }+\sqrt{{\lambda }_2{\rho }_4}>0,$$

$$2{\rho }_2\sqrt{{\lambda }_2}+2{\rho }_2^{\prime }\sqrt{{\lambda }_1}+{\lambda }_3\sqrt{{\rho }_4}+2\left(\sqrt{{\lambda }_1{\lambda }_2{\rho }_4}+\sqrt{\left({\lambda }_3+2\sqrt{{\lambda }_1{\lambda }_2}\right)\left({\rho }_2+\sqrt{{\lambda }_1{\rho }_4}\right)\left({\rho }_2^{\prime }+\sqrt{{\lambda }_2{\rho }_4}\right)}\right)>0,$$

$$\displaystyle 2\rho_{2}\sqrt{\lambda_{2}}+2\rho_{2}^{\prime}\sqrt{\lambda_{1}}+\left(\lambda_{3}+\lambda_{4}-\lambda_{5}\right)\sqrt{\rho_{4}}+2\Big{(}\sqrt{\lambda_{1}\lambda_{2}\rho_{4}}$$

$$\displaystyle+\sqrt{\left(\left(\lambda_{3}+\lambda_{4}-\lambda_{5}\right)+2\sqrt{\lambda_{1}\lambda_{2}}\right)\left(\rho_{2}+\sqrt{\lambda_{1}\rho_{4}}\right)\left(\rho_{2}^{\prime}+\sqrt{\lambda_{2}\rho_{4}}\right)}\Big{)}>0.$$

$$m_{H_0}+m_{A_0}>m_Z,m_{H^{\pm }}>79.3\text{GeV}.$$

$$R_1={\cos }^4\theta \frac{{\Gamma }^{\text{SM}}}{\Gamma },$$

$$\Gamma ={\cos }^2\theta {\Gamma }^{\text{SM}}+{\Gamma }^{\text{inv}},$$

Full text

1 Introduction

Recently, with the detection of a gravitational wave (GW) event [1, 2] from a black hole binary merger with the Laser Interferometer Gravitational-Wave Observatory (LIGO) the era of GW-astronomy has begun. The production of GW can also be associated with several other different mechanisms such as topological defects of the domain walls and cosmic strings [3], inflationary quantum fluctuations [4], preheating [5], first-order phase transitions in the early Universe [6, 7], etc. In this work, we explore production of GW from the first-order phase transitions in the early Universe for a proposed dark matter particle physics model that is constructed by simple extension of Standard Model (SM) by an inert doublet and a scalar singlet. Also a strong first-order electroweak phase transition helps to explain electroweak baryogenesis [8]-[11]. It is to be noted that the first-order cosmological phase transition originates from the bubble nucleation of a true vacuum state at a temperature known as the nucleation temperature at which the probability for a single bubble to nucleate within the horizon volume is of the order one. Initially, the bubbles are considered to have all possible shapes with different surface tension and internal pressure. The bubbles with the size just large enough for avoiding collapse are considered as a critical bubble. The bubbles which are smaller than the critical bubble tend to collapse whereas the larger bubbles tend to expand due to the pressure difference between the false and true vacua. During the collisions, the bubbles cannot retain their spherical symmetry which initiates phase transitions and as a result, gravitational waves are produced. The details of the GW production mechanisms are discussed in [12]-[26]. The GWs are produced from the strong first-order electroweak phase transition mainly by the following three mechanisms: bubble collisions [12]-[17], sound waves in the hot plasma [18]-[21] and magnetohydrodynamic turbulence of bubbles [22]-[26] in the early Universe.

The electroweak phase transition as explained by SM of particle physics initiated by spontaneous symmetry breaking through Higgs mechanism when the Higgs acquires a vacuum expectation value (VEV) is a smooth crossover and not a first-order phase transition [27]-[29]. However extension of SM by adding scalar singlet field or a Higgs like doublet can induce strong first-order phase transitions. The phase transitions can be of two types: a) a one-step process, which involves only initial and final phases and b) a two or multi-step processes, which involves one or more intermediate metastable phases along with the initial and final phases [30]-[38]. In order to explore the gravitational wave production from the first-order phase transition, many authors have considered different types of particle dark matter (DM) models [39]-[52]. In this work, we extend the SM by adding an extra Higgs doublet and a real singlet scalar. The added doublet is an inert doublet [53]-[63] in the sense that it does not have any direct coupling with fermion. A $Z_{2}$ symmetry is imposed to make it stable. The lightest stable inert particle is attributed to a viable particle candidate of dark matter. The extra singlet scalar mixes with the SM Higgs. This model has already been discussed in previous works [64, 65]. We demonstrate that this model, in addition to provide a viable particle candidate, induces strong first-order electroweak phase transition. Moreover, introducing a new scalar particle increases the degrees of freedom in the thermal plasma and improves the strength of the electroweak phase transition. In this work, the mass of the second physical scalar $h_{2}$, appearing due to the interaction between singlet scalar and the SM Higgs has been considered heavier as well as lighter than the Higgs mass. We constrain the model parameters by using vacuum stability [66], perturbativity, Large Electron-Positron Collider (LEP)[67] and Large Hadron Collider (LHC) bounds, PLANCK bound on the DM relic density [68] and the limits given by spin independent direct detection experiments like XENON-1T [69], PandaX-II [70], LUX [71] and DarkSide-50 [72]. We also use limits from future generation colliders such as High Luminosity LHC (HL-LHC), International Linear Collider (ILC) and China Electron Positron Collider (CEPC) or TLEP [40]. We choose some benchmark points (BPs) from the allowed parameter space to calculate the GW production due to the first-order phase transition induced by the present model. We also discuss the detectability of such GW by the future space interferometers such as Big Bang Observer (BBO) [73, 74], Evolved Laser Interferometer Space Antenna (eLISA) [73], Advanced Laser Interferometer Antenna (ALIA) [75], DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) [73, 76], Ultimate-DECIGO (U-DECIGO) [77] and ground-based detector Advanced LIGO (aLIGO) [73, 78].

The paper is organised as follows. In Sect. 2, we present the extension of an inert doublet model by introducing a singlet scalar and derive the relations between model parameters. In Sect. 3, we discuss both the theoretical and experimental bounds that we have used to constrain the model parameter space. The calculations of relic density and direct detection cross-section of this extended inert doublet model (IDM) are given in Sect. 4. In this section, we also discuss the viable model parameter space from all the constraints mentioned in Sect. 3. In Sect. 5, we present the finite temperature effective potential to study the electroweak phase transitions in our model. The production mechanisms of GWs from the first-order phase transitions are also furnished in this section. Finally, we summarize and conclude our work in Sect. 6.

2 The Model

In this work we extend the SM of particle physics by an extra Higgs doublet $\Phi_{I}$ and a real singlet scalar $S$. While $\Phi_{I}$ is $Z_{2}$ odd, the SM and the other added scalar singlet is $Z_{2}$ even. The extra doublet does not acquire any VEV, while the added singlet acquires a VEV on spontaneous symmetry breaking and mixes with SM Higgs. The dark matter candidate is the lightest of the two neutral scalars of the inert doublet. The potential of the scalar sector of the model can be expressed as

[TABLE]

As mentioned, after spontaneous symmetry breaking SM Higgs field $\Phi_{H}$ acquires a non zero VEV $v=$246.22 GeV and also the scalar $S$ acquires a VEV $v_{s}$. The SM Higgs, additional Higgs doublet and the scalar particle is represented as

[TABLE]

where $G^{+}$ and $G_{0}$ are the charged and neutral Goldstone bosons respectively. After SSB the Goldstone bosons are absorbed in the gauge bosons ($W^{\pm},Z$) [65]. After minimising the scalar potential represented in Eq. (1) using the conditions

[TABLE]

we obtain

[TABLE]

The mass matrix of the scalar sector is obtained by calculating the second-order derivatives $\left(\dfrac{\partial^{2}V}{\partial{h}^{2}},\hskip 2.84526pt\dfrac{\partial^{2}V}{\partial{s}^{2}},\hskip 2.84526pt\dfrac{\partial^{2}V}{\partial{H^{\pm}}^{2}},\hskip 2.84526pt\dfrac{\partial^{2}V}{\partial{H_{0}}^{2}},\hskip 2.84526pt\dfrac{\partial^{2}V}{\partial{A_{0}}^{2}},\hskip 2.84526pt\dfrac{\partial^{2}V}{\partial{h}\partial{s}}\right)$ of the scalar potential (Eq. (1)) and the elements are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

As $h$ and $s$ mix, we diagonalise the mass matrix in $h$, $s$ basis by a unitary matrix $U$

[TABLE]

to obtain two physical mass eigenstates $h_{1}$ and $h_{2}$ as

[TABLE]

where $\theta$ is the mixing angle that can be computed from

[TABLE]

The expressions for the mass eigenstate of the two physical scalars $h_{1}$ and $h_{2}$ are given as

[TABLE]

where the ‘+’ sign is for $h_{1}$ and ‘$-$’ sign is for $h_{2}$. In the present work, $h_{1}$ is attributed to the SM like Higgs boson with mass $m_{h_{1}}=125.09$ GeV [79] and $h_{2}$ is the other scalar with mass $m_{h_{2}}$. Here we consider both the cases when $m_{h_{2}}>m_{h_{1}}$ and when $m_{h_{2}}<m_{h_{1}}$. Considering the coupling $\lambda_{5}$ (in Eq. (1)) to be less than zero, we get $H_{0}$ to be the lightest stable particle and the dark matter candidate in our present work. Using Eqs. (4)-(14) we obtain the following relations

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

3 Constraints

In this section various theoretical and experimental bounds are given. These are used to constrain the model parameter space.

3.1 Theoretical Constraints

Vacuum Stability**

From the vacuum stability conditions the bounds on the couplings are given as [66]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Perturbativity**

All the quartic couplings in the tree-level potential (Eq. (1)) must be less than $4\pi$ to be consistent with the perturbative conditions.

3.2 Experimental Constraints

Collider Constraints**

From the LEP experiment, the bound on the model parameter space is given as **[67]**

[TABLE]

The bounds are also obtained from the LHC experimental results. The signal strength of the SM like Higgs $h_{1}$ in the present model can be expressed as

[TABLE]

where $\Gamma^{\text{SM}}$ and $\Gamma$ are the total SM Higgs decay width and total decay width of SM like Higgs boson of mass 125.09 GeV. The expression of $\Gamma$ can be written as

[TABLE]

where $\Gamma^{\text{inv}}$ is the invisible Higgs decay width. In our case, there are two possible invisible decay channels of $h_{1}$, one of them is $\Gamma^{inv}\left(h_{1}\rightarrow H_{0}H_{0}\right)$ (for $m_{H_{0}}\leq m_{h_{1}}/2$, $m_{H_{0}}$ being the mass of the dark matter particle $H_{0}$) and the other one is $\Gamma^{inv}\left(h_{1}\rightarrow h_{2}h_{2}\right)$ (for $m_{h_{2}}\leq m_{h_{1}}/2$) and they are expressed as

[TABLE]

and

[TABLE]

where the expressions for the couplings $g_{h_{1}H_{0}H_{0}}$ and $g_{h_{1}h_{2}h_{2}}$ are given as

[TABLE]

with

[TABLE]

and

[TABLE]

The invisible decay branching fraction for SM like scalar can be expressed as

[TABLE]

We adopt the bounds on the invisible decay branching fraction for SM scalar to be $\text{Br}_{\text{inv}}\leq 24\%$ **[80]** (for $m_{h_{1}}\geq m_{H_{0}}/2$). The limit on the scalar mixing $\sin\theta\leq 0.4$ **[81]-[83]** and the signal strength of SM Higgs $R_{1}\geq 0.84$ **[84, 85]** is obtained from LHC experimental results. However, study of future collider experiments can further constrain the scalar mixing angle **[40]**. From the study of Ref. **[40]** it is found that High Luminosity LHC (HL-LHC) restricts the mixing angle to be $\cos\theta\geq 0.94$. Moreover bounds from International Linear Collider (ILC) is found to be more stringent. Projected limit on mixing angle is found to be $\cos\theta\simeq 0.98$ for ILC-1 operating at $\sqrt{s}=250$ GeV and $\cos\theta\simeq 0.99$ for ILC-3 at $\sqrt{s}=1$ TeV. The mixing angle is further constrained to $\cos\theta\sim 0.997$ when China Electron Positron Collider (CEPC) or TLEP sensitivity is taken into account. Therefore, we consider $\sin\theta\leq 0.063$ ($\cos\theta\geq 0.998$), consistent with the limits from future collider experiments.

PLANCK constraint on relic density**

The relic density of dark matter candidate $H_{0}$ must satisfy the PLANCK observational limit for this to be a viable candidate of dark matter. The relic density limit given by PLANCK observation is

[TABLE]

where $\Omega_{\text{DM}}$ is the DM relic density normalised by the critical density of the Universe and ${\rm{h}}$ is the Hubble parameter normalised by a value of 100 Km s*-1** Mpc**-1** [68].*

Direct Searches of Dark matter**

Direct detection experiments of dark matter put an upper bound on dark matter nucleon elastic scattering cross-sections for different dark matter masses. In the present work, we consider the results of the following direct detection experiments of dark matter to constrain the model parameter space XENON-1T **[69]**, PandaX-II **[70]**, LUX **[71]** and DarkSide-50 **[72]**.

4 Dark matter phenomenology

In this section we furnish the dark matter relic density and direct detection scattering cross-section formulas of the present extended inert doublet model with an additional real singlet scalar. These will be used to compute and constrain the model parameter space.

4.1 Relic Density

In order to calculate the dark matter relic density one needs to solve the Boltzmann equation which can be expressed as **[86]**

[TABLE]

where $n_{H_{0}}$ is the dark matter number density and $n_{H_{0}}^{\text{eq}}$ is the dark matter number density at thermal equilibrium. In Eq. (36), $\langle\sigma v\rangle$ is the thermal average total annihilation cross-section $\sigma$ times the relative velocity $v$ of the dark matter candidate ($H_{0}$) and $H$ is the Hubble parameter. The expressions for $\langle\sigma v\rangle$ at a temperature $T$ can be written as

[TABLE]

where $\sigma(s)$ is the total annihilation cross-section of the dark matter particle $H_{0}$, $\sqrt{s}$ is the centre of mass energy, $K_{1}$ and $K_{2}$ are the first and second-order modified Bessel functions respectively. The relic density of the dark matter candidate $H_{0}$ can be computed using the expressions which are given as

[TABLE]

with

[TABLE]

where $Y_{F}$ is the value of $Y$ at the freeze-out temperature $T_{F}$, $G$ the universal gravitational constant, $g_{*}$ is the degrees of freedom and $T_{0}$ is the temperature at the present epoch.

For the purpose of computations we used the publicly available packages namely FeynRules **[87]** where the extended IDM is implemented and finally micrOMEGAs **[88]** (which includes the package CalcHEP within its framework) to calculate the relic density. In order to calculate the total annihilation cross-sections of the dark matter candidate, the annihilation channels $H_{0}H_{0}\rightarrow W^{+}W^{-},\hskip 2.84526ptZZ,\hskip 2.84526ptf\bar{f}$ (where $W$, $Z$ denote the gauge bosons and $f$ denotes the SM fermions) as well as $H_{0}H_{0}\rightarrow h_{1}h_{1},\hskip 2.84526pth_{2}h_{2},\hskip 2.84526pth_{1}h_{2}$ via s, t and u channels as also the relevant four-point interactions have been considered. The micrOMEGAs **[88]** has been used to estimate the relic-density for dark matter candidate $H_{0}$ considering all the annihilation channels mentioned above.

4.2 Direct detection

The dark matter particle $H_{0}$ interacts with the SM particles via Higgs exchange. The expression for the spin independent elastic scattering cross-section is given as

[TABLE]

where $m_{N}$ is the nucleon mass and $f$ is the nucleon form factor which has been approximated as 0.3 **[89]**. In Eq. (40) the couplings $\lambda_{h_{1}H_{0}H_{0}}$ and $\lambda_{h_{2}H_{0}H_{0}}$ can be written as

[TABLE]

[TABLE]

We calculate the dark matter scattering cross-sections using Eqs. (40)-(42) and then check whether they satisfy the direct detection cross-section bounds obtained from DM direct detection experiments such as XENON-1T **[69]**, PandaX-II **[70]**, LUX **[71]** and DarkSide-50 **[72]**.

4.3 Constraining the parameter space

The parameter space is first constrained with the theoretical bounds given in Sect. 3.1. It is then further constrained by the collider bounds given in Sect. 3.2. The constrained parameter space thus obtained is then used to calculate the relic densities of the dark matter (Sect. 4.1) in this model. The parameter space is finally constrained by comparing these calculated relic densities with PLANCK results (Eq. (35)). We have checked that the scattering cross-sections and their variations with the dark matter masses calculated using this finally constrained parameter space lie below the upper bounds of such variations given by different dark matter direct detection experiments such as XENON-1T **[69]**, PandaX-II **[70]**, LUX **[71]** and DarkSide-50 **[72]**.

In the present work we choose four benchmark points (BPs) from the constrained parameter space and use them to calculate the intensities and frequencies of the gravitational wave that would be produced by this first-order phase transition initiated by the present dark matter model. These are given in Table 1. For the purpose of demonstration, we calculate the relic densities for each of the BPs and show their variations with the mass of the DM candidate $m_{H_{0}}$ in this model. They are given in left panels of Figures 1-4. As mentioned, the relic densities are computed using micrOMEGAs code **[88]**. We also show in the same Figures (right panels) the variations of dark matter nucleon scattering cross-sections as a function of dark matter masses with the same BPs.

The dips in the relic density plots (left panels of Figures 1-4) attribute to the fact that at these dark matter masses the annihilation cross-section suffers a sharpe increase resulting in reduction in relic densities. From the left panels of Figures 1-4, it is evident that when the DM mass ($m_{H_{0}}$) is close to ($m_{h_{1}}/2$) and $m_{h_{2}}/2$, sudden reductions in the relic density occur. This is because at these DM masses annihilation cross-sections increase significantly. For example, for $m_{h_{2}}=100$ GeV the dip in the relic density for the dark matter mass around 62.5 GeV and 50 GeV are due to the resonant increase of annihilation cross-section. In addition, dips in the relic density plot for DM masses close to $m_{W}$ (mass of $W$ boson), $m_{Z}$ (mass of $Z$ boson), $m_{h_{1}}$ (mass of Higgs boson) and $m_{t}$ (mass of top quark) can also be noticed. At these DM masses, new annihilation channels, $H_{0}H_{0}\rightarrow W^{+}W^{-}$, $H_{0}H_{0}\rightarrow ZZ$, $H_{0}H_{0}\rightarrow h_{1}h_{1}$ and $H_{0}H_{0}\rightarrow t\bar{t}$ open up. Besides, when $m_{H_{0}}$ is higher than $m_{W}$, relic density decreases and becomes underabundant due to the increase in the annihilation cross-section. The DM relic abundance condition is again satisfied at the higher values of DM masses ($m_{H_{0}}$). For example, in case of BP1 DM particle reaches required relic abundance at $m_{H_{0}}$=750 GeV.

From the right panels of Figures 1-4, one observes DM direct detection scattering cross-section ($\sigma_{\text{SI}}$) decreases as the DM mass increases. The nature of these plots can be easily understood from Eq. (40). The right panels of Figures 1-4 show that the calculated $\sigma_{\text{SI}}$ for four BPs are consistent with the bounds given by dark matter direct detection experiments like XENON-1T, LUX, PandaX-II and DarkSide-50. Thus $H_{0}$ can be considered as a viable particle candidate of dark matter with the mass of the order of GeV - TeV.

5 Electroweak Phase transition and Gravitational Waves Production in Extended Inert Doublet Dark Matter Model

In this section we explore the electroweak phase transition and production of GWs from the considered dark matter model.

5.1 Effective Potential

To study the electroweak phase transition (EWPT) in the present model we add the finite temperature correction with the tree-level potential (Eq. (1)). The finite temperature effective potential can be written as **[90]**

[TABLE]

where $V_{1-loop}^{T=0}$ and $V_{1-loop}^{T\neq 0}$ are the one-loop corrected potential at zero temperature and at finite temperature respectively. The one-loop effective potential at zero temperature is given by **[90]**

[TABLE]

where the ‘+’ sign is for bosons and ‘-’ sign is for fermions. Here, $n_{i}$ is the number of degrees of freedom, $m_{i}$ is the field-dependent masses of the particles and $i\equiv(h,H_{0},A_{0},H^{\pm},s,W,Z,t$). The degrees of freedom of these particle species are $n_{W^{\pm}}=6$, $n_{Z}=3$, $n_{t}=12$ and $n_{h,H_{0},A_{0},H^{+},H^{-},s}=1$. The field-dependent masses of the gauge bosons and scalars in terms of $v$ (VEV of $\phi_{H}$), $v_{\phi}$ (VEV of $\phi_{I}$) and $v_{s}$ (VEV of $S$) at $T=0$ are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $g$ and $g^{\prime}$ denote the $SU(2)_{L}$ and $U(1)_{Y}$ gauge couplings of the SM respectively. Note that, for $v_{\phi}=0$, the Eqs. (47) - (51) reduce to the expressions given in Eqs. (5) - (9). We work in the Landau gauge ($\xi=0$), where the Goldstone masses become zero at $T=0$ and also there are no ghost contributions **[91]**. However, the gauge invariant methods ($\xi\neq 0$) **[92]** can also be applied. The quantity $Q$ denotes the renormalisable scale which we take $Q=246.22$ GeV in our calculations. However, the selection of $Q$ is not unique. One can vary $Q$ and check the uncertainty due to that in different thermal parameters **[93, 94]**. We varied $Q$ for a benchmark point and presented the sensitivity for different parameters in Sect. 5.3. In Eq. (44), $C_{i}$ represents a numerical constant. For $W,\hskip 2.84526ptZ$ boson the $C_{W,Z}=5/6$ and for the other particles $C_{h,H_{0},A_{0},H^{\pm},s,t}=3/2$. The one-loop effective potential at finite temperature has the form **[90]**

[TABLE]

where

[TABLE]

We include the daisy resummation procedure in $V_{1-loop}^{T\neq 0}$ **[95]** by adding the temperature dependent terms to the boson masses. Now, the corrected thermal masses **[52]**, **[96]-[98]** are $\mu_{1}^{2}(T)=m_{1}^{2}+c_{1}T^{2}$, $\mu_{2}^{2}(T)=m_{2}^{2}+c_{2}T^{2}$ and $\mu_{3}^{2}(T)=m_{s}^{2}+c_{3}T^{2}$, where

[TABLE]

[TABLE]

[TABLE]

here $y_{t}$ denotes the top Yukawa coupling. In this work we use the CosmoTransitions package **[90]** to compute the finite temperature corrections to the tree-level potential.

5.2 Gravitational Waves Production from Dark Matter

The bubble nucleation of a true vacuum state (from several metastable states) with surface tension and internal pressure at a temperature known as the nucleation temperature is central to the first-order cosmological phase transition. They can produce different sizes small and large the former types of which tend to collapse whereas the latter types tend to expand after attaining the criticality. These are the latter types, the critical bubbles collide with each other and upon losing their spherical symmetries in the process drive the phase transition and eventual emission of gravitational waves.

The bubble nucleation rate per unit volume at a particular temperature can be written as **[99]**

[TABLE]

where $\Gamma_{0}\left(T\right)\propto T^{4}$, $S_{3}\left(T\right)$ is the Euclidean action of the critical bubble. The Euclidean action $S_{3}\left(T\right)$ can be expressed as **[99]**

[TABLE]

where $V_{eff}$ is the effective finite temperature potential represented in Eq. (43). Nucleation of the bubble occurs at the nucleation temperature $T_{n}$ if it satisfies the condition $S_{3}\left(T_{n}\right)/T_{n}\approx 140$ **[90]**.

The gravitational waves are produced from the first-order phase transition mainly by the three mechanisms namely bubble collisions **[12]-[17]**, sound wave **[18]-[21]** and turbulence in the plasma **[22]-[26]**. The total GW intensity $\Omega_{\text{GW}}{\rm h}^{2}$ as a function of frequency can be expressed as the sum of the contributions from the three components **[12]-[26]**

[TABLE]

The component from the bubbles collision $\Omega_{\text{col}}{\rm h}^{2}$ is given by

[TABLE]

with the parameter $\beta$

[TABLE]

where $T_{n}$ is the nucleation temperature and $H_{n}$ is the Hubble parameter at $T_{n}$. A conservative estimate of the bubble wall velocity $v_{w}$ can be expressed as **[16, 45, 48, 100]**

[TABLE]

Here in this work, we use the above mentioned expression (Eq. (64)) for calculating the bubble wall velocity. However, in some literatures, $v_{w}$ is taken to be 1 for simplicity (**[49, 50, 101]**). A more detailed discussion on the choice of $v_{w}$ can be found in **[102]** which results in $v_{w}$ $\gtrsim$ 0.2. The parameter $\kappa$ in Eq. (62) represents the fraction of latent heat deposited in a thin shell which is given by

[TABLE]

with **[49, 103]**

[TABLE]

In the above, $v_{n}$ is the vacuum expectation value of Higgs at $T_{n}$ and $m_{W}$, $m_{Z}$ and $m_{t}$ are the masses of W, Z and top quarks respectively. The parameter $\alpha$ is defined as the ratio of vacuum energy density $\rho_{\text{vac}}$ released by the electroweak phase transition to the background energy density of the plasma $\rho_{*}^{\text{rad}}$ at $T_{n}$. The expression of $\alpha$ has the form

[TABLE]

with

[TABLE]

and

[TABLE]

The parameter $\rho_{vac}$ can also be calculated using the trace anomaly **[21]**, where the same requires an additional factor of $1/4$ in front of the term $T\dfrac{dV}{dT}$ (Eq. (68)). The quantity $f_{\text{col}}$ in Eq. (62) is the peak frequency produced by the bubble collisions which takes the form

[TABLE]

The sound wave (SW) component of the gravitational wave (Eq. (61)) is given by

[TABLE]

where $\kappa_{v}$ is the faction of latent heat transformed into the bulk motion of the fluid which has the following form

[TABLE]

In Eq. (71) $f_{\text{SW}}$ denotes the peak frequency produced by the sound wave mechanisms which takes the form

[TABLE]

The component from the turbulence in the plasma $\Omega_{\text{turb}}{\rm h}^{2}$ is given by

[TABLE]

where $\epsilon=0.1$ and $f_{\text{turb}}$ denotes the peak frequency produced by the turbulence mechanism which can be written as

[TABLE]

In Eq. (74) the parameter $h_{*}$ has the following form

[TABLE]

We repeat our analysis by considering an updated expression of $\kappa_{v}$ **[104, 105]**, where the prefactor $\dfrac{\alpha_{\infty}}{\alpha}$ is taken to be 1 and furnish our results in Sect. 5.3. The updated expression of $\kappa_{v}$ is given by **[104, 105]**

[TABLE]

In this work, Eqs. (61)-(77) are used for the computation of gravitational wave intensity.

5.3 Calculations and Results

As mentioned earlier the computation of GW intensity from the first-order phase transition in the present particle DM model (inert doublet and a scalar singlet extended SM) has been performed for four chosen benchmark points for the model parameters given in Table 1.

We calculate gravitational wave intensity from the model and compare it with the sensitivity curves of different GW detectors such as BBO, eLISA, ALIA, DECIGO, U-DECIGO and aLIGO. The GW intensity depends mainly on factors like strength of the first-order phase transition (the parameter $\alpha$), the time-scale of the phase transition (the parameter $1/\beta$), bubble wall velocity $v_{w}$, nucleation temperature $T_{n}$ and Higgs VEV $v_{n}$ at the nucleation temperature $T_{n}$. In order to calculate the GW intensity, we first calculate the transition temperature of the first-order phase transition. In calculating this, the finite temperature effective potential (Eq. (43)) is first computed. For these calculations we use a publicly available package namely the Cosmotransition package **[90]**. The tree-level potential (Eq. (1)) serves as an input to this package and provides the parameters related to the phase transition. The GW intensity is estimated by using Eqs. (61)-(77). In this work, we have selected four BPs (Table 1) such that they satisfy all constraints mentioned in Sect. 3. The relic abundance $\Omega_{\text{DM}}\text{h}^{2}$ and the direct detection scattering cross-section $\sigma_{\text{SI}}$ obtained from each of the four BPs are also given in Table 1.

In Figure 5 we show the phase transition properties for BP1. The left and right panels of Figure 5 show the phase structure of the model and the tunneling profile as a function of the bubble radius respectively. As can be seen in the left panel of Figure 5, there exist two transition temperatures at $T_{n}=119.86$ GeV and $T_{n}=226.55$ GeV. We only consider the lower of the two as the phase transition temperature because the low nucleation temperature is more sensitive to probe the GW signal. An electroweak phase transition occurs when the temperature of the Universe drops which results in a separation of potential between a high phase and a low phase by a potential barrier. In this case, a first-order phase transition occurs at the nucleation temperature $T_{n}=119.86$ GeV. The phase transition properties for BP2, BP3 and BP4 are also studied but no significant differences from BP1 are noted. The computed values of the parameters ($v_{n},\hskip 2.84526ptT_{c},\hskip 2.84526ptT_{n},\hskip 2.84526pt\alpha,\hskip 2.84526pt\beta/H$) corresponding to each of the four BPs to be uesd for the calculations of GW intensity are furnished in Table 2.

From Table 2, one can see that the nucleation temperature $T_{n}$ is smaller than the critical temperature $T_{c}$ (the temperature at which there exist two degenerate minima separated by a potential barrier) for each of the BPs. Since, the selection of $Q$ is not unique, we have varied $Q$ form 68 GeV to 270 GeV and found that the parameters have the following uncertainty: 17.6 %, 14.8 % and 0.9 % for $T_{c}$, $T_{n}$ and $v_{n}$ respectively. This indicates that the thermal parameters can vary significantly with $Q$.

Using the computed numerical values for each of the four sets given in Table 2 (corresponding to the BPs given in Table 1) the GW intensities and frequencies from first-order phase transitions are now calculated from Eqs. (61)-(77) for each of the sets. In order to test whether the sound wave contribution to the overall GW signal is significant, the suppression factor ($HR_{*}/\bar{U}_{f}$, where $R_{*}$ is the mean bubble separation and $\bar{U}_{f}$ is the root-mean- square (RMS) fluid velocity) is calculated **[103, 104, 106]**. In our model, the sound wave contribution of the overall GW amplitude turns out to be 1.76 ($>$* 1) for BP1. This implies that the sound wave contribution, in this case, survives more than a Hubble time (long-lasting sound waves) and hence contribute significantly. The applicability of this issue can be tested through a dedicated numerical lattice simulation. The variations of GW intensities with the frequencies for different BPs are shown in left and right panels of Figure 6, where the GW intensities for different BPs are calculated using the expressions of $\kappa_{v}$ following Eq. 72 and 77 respectively. The sensitivity curves of GW detectors (BBO, eLISA, ALIA, DECIGO, U-DECIGO, aLIGO and LISA) are also shown in both of the panels of Figure 6 for comparison. One of the ways to graphically represent the sensitivity of a GW detector is called the power-law-integrated sensitivity curve [107] which calculates power across all frequencies in the sensitivity band by integrating the noise-weighted signal over frequency. However, we have presented our results in such a way that the sensitivity curves for different GW detectors can directly be compared with the stochastic backgrounds and contribution from other types of sources by following Ref. [73]. The future generation space-based GW detectors such as eLISA, BBO, ALIA, and DECIGO are expected to explore the frequency range from millihertz to decihertz. These space missions complement the effort of future generation ground-based GW detector aLIGO by providing the low frequency coverage with much higher sensitivity. This is because of the longer arms of the space-based detectors and absence of the seismic noise. The combined effort of the future space and ground-based telescopes, operational at different energies (eLISA with four different configurations such as N1A1M2L4, N2A2M5L4, N2A1M5L6 and N2A5M5L6 - millihertz band, ALIA, BBO, DECIGO - decihertz band, and aLIGO - decahertz band) will provide unprecedented detection of GWs. In addition, the recently proposed variant of DECIGO based on Fabry-Perot interferometry, is expected to significantly reduce the signal-to-noise ratio of GW backgrounds. As a result of which the best sensitivity in the frequency range 0.1 hertz to 10 hertz will be around $\Omega_{\text{GW}}\text{h}^{2}=10^{-16}-10^{-15}$. It can be seen from Figure 6 that for different BPs, GW intensities attain peaks at different frequencies. The peak of the GW intensity for BP1, BP2, BP3, and BP4 appears at frequencies $1.7\times 10^{-3}$ Hz, $2.3\times 10^{-3}$ Hz, $3.5\times 10^{-3}$ Hz and $5.4\times 10^{-3}$ Hz respectively. We obtained higher GW intensity for BP1 as compared to the other BPs. As seen from Figure 6 the GW intensities for all the BPs (BP1, BP2, BP3 and BP4) lie within the sensitivity curves of the N2A5M5L6 configuration of eLISA, BBO and U-DECIGO. However, BP1 is special because it shows higher intensity than the rest of the cases. It may be inferred from Figure 6 and Table 2 that the GW intensity mainly depends on $\beta$. For BP1, the value of the parameter $\beta$ is the smallest and the corresponding GW intensity is the highest along with the lowest peak frequency at $1.7\times 10^{-3}$. The choice of Eq. (77) over Eq. (72) introduces slight increase in the GW intensity shown in the right panel of Figure 6 as compared to the left panel. Small increase in the falling part of the GW spectrum is observed for three BPs namely BP1, BP2 and BP4 are attributed to the marginal increase of the GW intensity due to the sound wave contribution. However, the primary peak intensities for different BPs remain the same for both of the panels.*

6 Summary and Conclusions

In this work, we have explored the possible production of GWs from the first-order phase transition of the early Universe from a viable particle dark matter model and its detectability with the future GW detectors (BBO, eLISA, ALIA, DECIGO, U-DECIGO and aLIGO). We discussed the first-order phase transition that may be initiated by the present extended SM where the scalar sector of the SM is modified by adding an additional inert scalar doublet and a real scalar singlet. A discrete $Z_{2}$ symmetry has been imposed on the scalar potential which enables the lightest stable inert particle to be a viable particle candidate of dark matter. Due to the imposition of the $Z_{2}$ symmetry, the inert particle does not interact with the SM particles and also it does not acquire any VEV. The added real singlet scalar acquires a VEV on spontaneous symmetry breaking and it mixes with the SM Higgs which resulted in two new physical scalar eigenstates $h_{1}$ and $h_{2}$. Here we considered $h_{1}$ as the SM like Higgs boson and the other scalar $h_{2}$ as the added physical scalar. We constrain the model parameters using the conditions for vacuum stability, perturbativity, LEP and LHC bounds, dark matter relic density as measured by PLANCK and the bounds from the direct detection experiments. The future collider sensitivity is also used to constrain the scalar mixing angle and check for viability of the model with DM phenomenology and Gravitational Wave experiments. From the model parameter space thus constrained, we choose four benchmark points from the constrained parameter space for the computation of GW intensity and explore the GW production from the first-order phase transition of the tree-level potential. We include the finite temperature corrections of the tree-level potential. We calculate the GW intensity for four BPs. It has been found that the GW intensity increases as $\beta$ (Eq. (63)) decreases. In addition, the lower value of $\beta$ also lowers the frequency at which the maximum GW intensity is produced. We then explore whether the GW results for the four BPs lie within the sensitivity range of the GW detectors and found that the GW signals are detectable by the following future generation detectors namely BBO, U-DECIGO and eLISA (configuration - N2A5M5L6) within the detectable range. In this work, we show that the extended inert doublet model under consideration can explain the dark matter as well as the strong first-order electroweak phase transition. Thus we found that our model constrained by limits from future collider experiments, is consistent with DM searches and also in agreement with the sensitivity of GW detectors. The results also imply that the detection of GW signals by future spaceborne detectors could also be useful to ascertain the nature of particle dark matter candidates.

Acknowledgments**

The authors would like to thank A. Dutta Banik and A. Biswas for useful discussions.