\section{Introduction}
\label{sec:introduction}
\input{introduction.tex}

\section{Related Work}
\label{sec:relatedwork}
\input{relatedwork.tex}

\section{Methodology}
\label{sec:methodology}
\input{methodology.tex}

\section{Experiments}
\label{sec:experiments}
\input{experiment.tex}

\section{Discussion}\label{conclusion}

We have proposed a simple yet effective attention distillation approach, \ie, SAD, to improve the representation learning of CNN-based lane detection models. SAD~is validated in various models (\ie, ENet, ResNet-18, ResNet-34, and ResNet-101) and achieves consistent performance gains in three popular benchmarks (\ie, TuSimple, CULane and BDD100K), demonstrating the effectiveness of SAD. The results show that SAD~can generally improve the visual attention of different layers in various networks.
It would be interesting to extend this idea to other tasks that demands fine-grained attention to details, such as image saliency detection and image matting.

\noindent
\textbf{Acknowledgement:} This work is supported by SenseTime Group Limited, the General Research Fund sponsored by the Research Grants Council of the Hong Kong SAR (CUHK 14241716), Singapore MOE AcRF Tier 1 (M4012082.020), NTU SUG, and NTU NAP.

{\small
\bibliographystyle{ieee}

\subsection{Dataset}

\begin{table*}[!t]
\caption{Basic information of three lane detection datasets.}
\vskip 0.1cm
\label{dataset_table}
\centering
\small{
\begin{tabular}{c|c|c|c|c|c|c|c}
\hline
Name & \# Frame & Train & Validation & Test & Resolution & Road Type & \# Lane $>$ 5 ? \\
\hline
\hline
TuSimple~\cite{tusimple} & 6, 408 & 3, 268 & 358 & 2, 782 & 1280 $\times$ 720 & highway & $\times$ \\
CULane~\cite{pan2017spatial} & 133, 235 & 88, 880 & 9, 675 & 34, 680 & 1640 $\times$ 590 & urban, rural and highway & $\times$ \\
BDD100K~\cite{yu2018bdd100k} & 80, 000 & 60, 000 & 10, 000 & 10, 000 & 1280 $\times$ 720 & urban, rural and highway &  $\surd$ \\
\hline
\end{tabular}
}
\vspace{-3ex}
\end{table*}

\begin{figure}[t]
  \centering
  \includegraphics[width=1.0\linewidth]{./dataset.pdf}
  \vskip -0.2cm
  \caption{Typical video frames of TuSimple, CULane and BDD100K datasets. Ground-truth lanes are marked in green color.}
  \centering
  \vskip -0.2cm
  \label{fig:dataset}
\end{figure}

\vspace{0.1cm}
\noindent
\textbf{Datasets}.
Figure~\ref{fig:dataset} shows several video frames of three datasets that we use in our experiments. They are TuSimple~\cite{tusimple}, CULane~\cite{pan2017spatial} and BDD100K~\cite{yu2018bdd100k}. TuSimple and CULane are widely used in the literature. Many algorithms~\cite{pan2017spatial, neven2018towards, ghafoorian2018gan} have been tested in TuSimple since it was the largest lane detection dataset before 2018. As to CULane, it contains many challenging driving scenarios like crowded road conditions or roads under poor lighting (see Fig.~\ref{fig:dataset}). BDD100K is originally designed for lane instance classification. However, since there are typically multiple lanes in an image and these lanes are usually very close to each other, using instance segmentation algorithms will yield inferior performance. Therefore, we choose to only detect lanes without differentiating lane instances for BDD100K.
We discuss the details of transforming the original ground truths for our task in the following section on implementation details.
Table \ref{dataset_table} summarizes their details. Note that the last column of Table~\ref{dataset_table} shows that TuSimple and CULane have no more than 5 lanes in a video frame while BDD100K typically contains more than 8 lanes in a video frame. Besides, TuSimple is relatively easy while CULane and BDD100K are more challenging considering the total number of video frames and road types. Note that the original BDD100K dataset provides 100K video frames, in which 70K are used for training, 10K for validation and 20K for testing. However, since the ground-truth labels of the testing partition are not publicly available, we keep the training set unchanged but use the original validation set for testing.
A new validation set is allocated separately from the training set, as shown in Table~\ref{dataset_table}.

\noindent
\textbf{Evaluation metrics}. To facilitate comparisons against previous studies, we follow the literature and use the corresponding evaluation metrics for each particular dataset.

\noindent \textit{1) TuSimple.} We use the official metric (accuracy) as the evaluation criterion. Besides, false positive ($FP$) and false negative ($FN$) are also reported. Accuracy is computed as~\cite{tusimple}: $Accuracy = \frac{N_{pred}}{N_{gt}}$,
where $N_{pred}$ is the number of correctly predicted lane points and $N_{gt}$ is the number of ground-truth lane points.

\noindent \textit{2) CULane.} Following~\cite{pan2017spatial}, to judge whether a lane is correctly detected, we treat each lane as a line with 30 pixel width and compute the intersection-over-union (IoU) between labels and predictions. Predictions whose IoUs are larger than 0.5 are considered as true positives (TP). Then, we use $F_{1}$ measure as the evaluation metric, which is defined as: $F_{1} = \frac{2 \times Precision \times Recall}{Precision + Recall}$,
where $Precision = \frac{TP}{TP + FP}$ and $Recall = \frac{TP}{TP + FN}$.

\noindent \textit{3) BDD100K.} Since there are typically more than 8 lanes in an image, we decide to use pixel accuracy and IoU of lanes to evaluate the performance of different models.

\noindent
\textbf{Implementation details}.
Following~\cite{pan2017spatial}, we resize the images of TuSimple and CULane to 368$\times$640 and 288$\times$800, respectively. As to BDD100K, we resize the image to 360$\times$640 to save memory usage. The lanes of BDD100K are labelled by two lines. Training the networks using the provided labels is tricky. Therefore, based on these two lines, we calculate the center lines as new targets. We dilate ground-truth lanes of the training set of BDD100K as 8 pixels to provide denser targets while keeping these of testing set unchanged (2 pixels). We use SGD~\cite{bottou2010large} to train our models and the learning rate is set to 0.01. Batch size is set as 12 and the total number of training episodes is set as 1800 for TuSimple and 60K for CULane and BDD100K. The cross entropy loss of background pixels is multiplied by 0.4. Loss coefficients $\alpha$, $\beta$, and $\gamma$ are set as 0.1. Since we select lane pixel accuracy and IoU as the evaluation criterion for BDD100K dataset, we alter the original segmentation branch to output binary segmentation maps to facilitate the evaluation on BDD100K. The lane existence prediction branch is also removed for the BDD100K evaluation.

We empirically found that several practical techniques, \ie, data augmentation and IoU loss, can considerably enhance the performance of CNN-based lane detection models. As to data augmentation, we use random rotation, random cropping and horizontal flipping to process the input images.
In our experiments, we apply the same segmentation losses and augmentation strategy to our method, SCNN, ResNet baselines, and deep supervision methods, to ensure a fair comparison. Since the source codes of LaneNet~\cite{neven2018towards} and EL-GAN~\cite{ghafoorian2018gan} are not available, we use their results reported in their papers.

\subsection{Results}

\begin{table}[!t]
\caption{Performance of different algorithms on TuSimple testing set. Here "R-18-SAD~" denotes ResNet-18 + SAD~and we use the same abbreviation in the following sections.}
\label{tusimple_table}
\centering
\small{
\begin{tabular}{c|c|c|c}
\hline
Algorithm & Accuracy & FP & FN \\
\hline \hline
ResNet-18~\cite{he2016deep} & 92.69\% & 0.0948 & 0.0822 \\
ResNet-34~\cite{he2016deep} & 92.84\% & 0.0918 & 0.0796 \\
ENet~\cite{paszke2016enet} & 93.02\% & 0.0886 & 0.0734 \\
LaneNet~\cite{neven2018towards} & 96.38\% & 0.0780 & 0.0244 \\
EL-GAN~\cite{ghafoorian2018gan} & 96.39\% & \textbf{0.0412} & 0.0336 \\
SCNN~\cite{pan2017spatial} & 96.53\% & 0.0617 & \textbf{0.0180} \\
\hline \hline
\textbf{R-18-SAD~(ours)} & 96.02\% & 0.0786 & 0.0451 \\
\textbf{R-34-SAD~(ours)} & 96.24\% & 0.0712 & 0.0344 \\
\textbf{ENet-SAD~(ours)} & \textbf{96.64\%} & 0.0602 & 0.0205 \\
\hline
\end{tabular}
}
\vspace{-2ex}
\end{table}

\begin{table*}[!t]
\caption{Performance ($F_{1}$-measure) of different algorithms on CULane testing set. For crossroad, only FP is shown. The second column denotes the proportion of each scenario in the testing set.}
\label{culane_table}
\centering
\small{
\begin{tabular}{c|c|c|c|c|c||c|c}
\hline
Category & Proportion & \textbf{ENet-SAD~} & \textbf{R-18-SAD~} & \textbf{R-34-SAD~} & \textbf{R-101-SAD~} & ResNet-101~\cite{he2016deep} & SCNN~\cite{pan2017spatial} \\
\hline \hline
Normal & 27.7\% & 90.1 & 89.8 & 89.9 & \textbf{90.7} & 90.2 & 90.6 \\
Crowded & 23.4\% & 68.8 & 68.1 & 68.5 & \textbf{70.0} & 68.2 & 69.7 \\
Night & 20.3\% & 66.0 & 64.2 & 64.6 & \textbf{66.3} & 65.9 & 66.1 \\
No line & 11.7\% & 41.6 & 42.5 & 42.2 & \textbf{43.5} & 41.7 & 43.4 \\
Shadow & 2.7\% & 65.9 & 67.5 & \textbf{67.7} & 67.0 & 64.6 & 66.9 \\
Arrow & 2.6\% & 84.0 & 83.9 & 83.8 & \textbf{84.4} & 84.0 & 84.1 \\
Dazzle light & 1.4\% & \textbf{60.2} & 59.8 & 59.9 & 59.9 & 59.8 & 58.5 \\
Curve & 1.2\% & 65.7 & 65.5 & \textbf{66.0} & 65.7 & 65.5 & 64.4 \\
Crossroad & 9.0\% & 1998 & 1995 & \textbf{1960} & 2052 & 2183 & 1990 \\
Total & -- & 70.8 & 70.5 & 70.7 & \textbf{71.8} & 70.8 & 71.6 \\
\hline \hline
Runtime (ms) & -- & \textbf{13.4} & 25.3 & 50.5 & 171.2 & 171.2 & 133.5  \\
Parameter (M) & -- & \textbf{0.98} & 12.41 & 22.72 & 52.53 & 52.53 & 20.72  \\
\hline
\end{tabular}
}
\vspace{-3ex}
\end{table*}

\begin{table}[!t]
\caption{Comparative results on BDD100K test set.}
\label{bdd100k_table}
\centering
\small{
\begin{tabular}{c|c|c}
\hline
Algorithm & Accuracy & IoU \\
\hline \hline
ResNet-18~\cite{he2016deep} & 30.66\% & 11.07 \\
ResNet-34~\cite{he2016deep} & 30.92\% & 12.24 \\
ResNet-101~\cite{he2016deep} & 34.45\% & 15.02 \\
ENet~\cite{paszke2016enet} & 34.12\% & 14.64 \\
SCNN~\cite{pan2017spatial} & 35.79\% & 15.84 \\
\hline \hline
\textbf{R-18-SAD~(ours)} & 31.10\% & 13.29 \\
\textbf{R-34-SAD~(ours)} & 32.68\% & 14.56 \\
\textbf{R-101-SAD~(ours)} & 35.56\% & 15.96 \\
\textbf{ENet-SAD~(ours)} & \textbf{36.56\%} & \textbf{16.02} \\
\hline
\end{tabular}
}
\vspace{-2ex}
\end{table}

Tables \ref{tusimple_table}-\ref{bdd100k_table} summarize the performance of our methods, \ie, ResNet-18-SAD, ResNet-34-SAD, and ENet-SAD~against state-of-the-art algorithms on the testing set of TuSimple, CULane and BDD100K datasets.
We also report the runtime and parameter count of different algorithm in Table~\ref{culane_table} so that we can compare the performance with the complexity of the model taken into account. The runtime is recorded using a single GPU (GeForce GTX TITAN X) and the final value of runtime is obtained after averaging the runtime of 100 samples.

It is observed that ENet-SAD~outperforms all baselines in BDD100K while achieving compelling performance in TuSimple and CULane. Considering that ENet-SAD~has 20 $\times$ fewer parameters and runs 10 $\times$ faster compared with SCNN on CULane testing set, the performance strongly suggests the effectiveness of SAD.
It is observed that ResNet-18-SAD and ResNet-34-SAD achieve slightly inferior performance to ENet-SAD~despite their larger model capacity. The is because ResNet-18 and ResNet-34 only use spatial upsampling as the decoder while ENet has a specially designed decoder for the task.
It is noteworthy that SAD also helps given a deeper model. Specifically, we apply SAD~to ResNet-101, and find that it increases the F$_{1}$-measure from 70.8 to 71.8 in CULane and the accuracy increases from 34.45\% to 35.56\% in BDD100K.

We show some qualitative results of our algorithm and some baselines in these three benchmarks. As can be seen in Fig.~\ref{fig:tusimple_result}, ENet-SAD~can detect lanes more precisely than ENet~\cite{paszke2016enet} in TuSimple and CUlane.
As can be seen in Fig. ~\ref{fig:bdd100k_result}, the output probability maps of ENet-SAD~are more compact and contain less noise compared with those of vanilla ENet and SCNN in poor lighting conditions. However, since many images in BDD100K contain more than 8 lanes and are collected in challenging scenarios like severe occlusion and poor lighting conditions, the performance of all algorithms is unsatisfactory and needs further improvement. In general, SAD~can improve the visual attention as well as the detection performance in challenging conditions like crowded roads and poor light conditions.

We also perform experiments that apply SAD and remove the effect of the P1 branch by blocking the gradient of the P1 branch from the main branch. Results show that ENet-SAD (without supervision from P1 branch) can still achieve 96.61$\%$ on TuSimple, 70.8 on CULane and 36.54$\%$ on BDD100K, which means the performance gains come mainly from SAD~itself.

\begin{figure}[t]
  \centering
  \includegraphics[width=\linewidth]{./tusimple_culane_result.pdf}
  \vskip -0.1cm
  \caption{Performance of different algorithms on (a) TuSimple and (b) CULane testing sets. The number below each image denotes the accuracy. Ground-truth lanes are drawn on the input image.}
  \centering
  \vskip -0.25cm
  \label{fig:tusimple_result}
\end{figure}

\begin{figure}[t]
  \centering
  \includegraphics[width=\linewidth]{./bdd100k_result.pdf}
  \vskip -0.2cm
  \caption{Performance of different algorithms on BDD100K testing set. We visualize the probability maps to better showcase the effect of adding self attention distillation. The brightness of the pixel indicates the probability of this pixel belonging to lanes. The number below each image denotes the pixel accuracy of lanes. Ground-truth lanes are drawn on the input image.}
  \centering
  \vskip -0.5cm
  \label{fig:bdd100k_result}
\end{figure}

\subsection{Ablation Study}

We investigate the effects of different factors, e.g., the mimicking path, on the final performance. Besides, we also perform extensive experiments to investigate the timepoint to introduce SAD~in the training process.

\noindent \textbf{Distillation paths of SAD.} We summarize the performance of performing SAD~between different blocks of ENet in Table \ref{attention_ablation_table}. We have a few observations. (1) SAD~works well in the middle and high-level layers. (2) Adding SAD~in low level layers will degrade the performance.
The reason why SAD~does not work in low-level layers is that these layers are originally designated to detect low-level details of the scene. Making them to mimic the attention maps of later layers will inevitably harm their ability of detecting local features since later layers encode more global information.
Besides, we also find that mimicking the attention maps of the neighbouring layer successively brings more performance gains compared with mimicking those of non-adjacent layers ($P_{23}$ + $P_{34}$ outperforms $P_{24}$ + $P_{34}$).
We conjecture that attention maps of neighbouring layers are closer from the semantic perspective compared with those of non-neighbouring layers (see Fig.~\ref{fig:attention_map}).

\noindent \textbf{Backward distillation.} We also tested another distillation scheme that makes higher layers to mimic lower layers. It decreases the performance of ENet from 93.02\% to 91.26\% in TuSimple dataset. This is not surprising as low-level attention maps contain more details and are more noisy. Having higher-level layers to mimic lower layers will inevitably interfere the global information captured in higher layers, hampering the crucial clues for the lane detection task.

\begin{table}[!t]
\caption{Performance of different variants of ENet-SAD~on TuSimple testing set. $P_{ij}$ denotes that the output of the $i$-$th$ block of ENet mimics the output of the $j$-th block.}
\label{attention_ablation_table}
\centering
\small{
\begin{tabular}{c|c||c|c||c|c}
\hline
Path & Accuracy & Path & Accuracy & Path & Accuracy \\
\hline \hline
$P_{12}$ & 91.22\% & $P_{23}$ & 94.72\% & $P_{23}$, $P_{24}$ & 95.38\% \\
$P_{13}$ & 91.36\% & $P_{24}$ & 94.63\% & \textbf{$P_{23}$, $P_{34}$} & \textbf{96.64\%} \\
$P_{14}$ & 91.47\% & $P_{34}$ & 95.29\% & $P_{24}$, $P_{34}$ & 96.52\% \\
\hline
\end{tabular}
}
\vspace{-2ex}
\end{table}

\noindent \textbf{SAD~ v.s. Deep Supervision.} We also compare SAD~with deep supervision~\cite{Xie_2015_ICCV}. Here, deep supervision denotes the algorithm that uses the labels directly as supervision for each layer in the network. More specifically, we use 1x1 convolution and bilinear upsampling to obtain the prediction of intermediate layers and use the cross entropy loss to train the intermediate outputs of the model. We empirically find that adding deep supervision in blocks 2 to 4 obtains the most significant performance gains.
As can be seen in Table~\ref{deep_ablation_table}, SAD~brings more performance gains than deep supervision in all three benchmarks.
We attribute this to the following reasons. Firstly, compared with labels that are considered sparse and rigid, SAD~provides softer attention targets that capture more contextual information that indicate the
scene structure. Distilling information from attention maps of later layers helps previous layers to grasp the contextual signals. Secondly, a SAD path offers a feedback connection from deeper layers to shallower layers. The connection helps facilitate reciprocal learning between successive layers through attention distillation.

\begin{table}[!t]
\caption{Performance of SAD~and deep supervision applied to different base models on TuSimple, CULane and BDD100K testing sets.}
\label{deep_ablation_table}
\centering
\small{
\begin{tabular}{c|c|c|c|c}
\hline
\multirow{2}*{Algorithm} & TuSimple & CULane & \multicolumn{2}{|c}{BDD100K} \\
\cline{2-5}
~ & Accuracy & Total & Accuracy & IoU \\
\hline \hline
ENet & 93.02\% & 68.4 & 34.12\% & 14.64 \\
ENet-Deep & 94.69\% & 69.6 & 35.61\% & 15.38 \\
ENet-SAD & \textbf{96.64\%} & \textbf{70.8} & \textbf{36.56\%} & \textbf{16.02} \\
\hline
R-18 & 92.69\% & 67.9 & 30.66\% & 11.07 \\
R-18-Deep & 94.14\% & 68.8 & 30.95\% & 12.23 \\
R-18-SAD & \textbf{96.02\%} & \textbf{70.5} & \textbf{31.10\%} & \textbf{13.29} \\
\hline
R-34 & 92.84\% & 68.6 & 30.92\% & 12.24 \\
R-34-Deep & 94.52\% & 69.2 & 31.72\% & 13.59 \\
R-34-SAD & \textbf{96.24\%} & \textbf{70.7} & \textbf{32.68\%} & \textbf{14.56} \\
\hline
\end{tabular}
}
\vspace{-3ex}
\end{table}

\noindent \textbf{When to add SAD.} Recall that we assume a half-trained model before we add SAD into the training. Here, we investigate the different timepoints to add SAD. As can be seen in Fig.~\ref{fig:when_add_SAD}, although different timepoints of introducing SAD  achieve almost the same performance in the end, the time to add SAD~ has an effect on the convergence speed of the networks. We attribute the phenomenon to the quality of the target attention maps produced by later layers. In earlier training stage, deeper layers have not been trained well and therefore the distillation targets produced by these layers are of low quality. Introducing SAD at these earlier stages is not as fruitful. Conversely, adding SAD in later training stage would benefit the representation learning of the previous layers.

\begin{figure}[t]
  \centering

  \includegraphics[width=0.8\linewidth]{./when_add_SAD_culane_new.pdf}
  \vskip -0.2cm
  \caption{Performance of adding self attention distillation on the ENet model at different training episodes on the CULane validation set. The number in the legend denotes the episode when self attention distillation is added. "Baseline" denotes the ENet model without self attention distillation.}
  \vskip -0.1cm
  \centering
  \label{fig:when_add_SAD}
\end{figure}

\subsection{Dataset}

\begin{table*}[!t]
\caption{A brief description about three lane detection datasets.}
\label{dataset_table}
\centering
\begin{tabular}{c|c|c|c|c|c|c|c}
\hline
Name & \# Frame & Train & Validation & Test & Resolution & Road Type & \# Lane $>$ 5 ? \\
\hline
TuSimple & 6, 408 & 3, 626 & -- & 2, 782 & 1280 $\times$ 720 & highway & $\times$ \\
\hline
CULane & 133, 235 & 88, 880 & 9, 675 & 34, 680 & 1640 $\times$ 590 & urban, rural and highway & $\times$ \\
\hline
BDD100K & 80, 000 & 70, 000 & -- & 10, 000 & 1280 $\times$ 720 & urban, rural and highway &  $\surd$ \\
\hline
\end{tabular}
\end{table*}

\begin{figure}[t]
  \centering
  \includegraphics[width=1.0\linewidth]{./figure/dataset.pdf}
  \vskip -0.2cm
  \caption{Typical video frames of TuSimple, CULane and BDD100K datasets.}
  \centering
  \label{fig:dataset}
\end{figure}

Table \ref{dataset_table} records the basic information of three lane detection datasets. Note that the last column of Table \ref{dataset_table} shows that TuSimple and CULane have no more than 5 lanes in a video frame while BDD100K typically has more than 8 lanes in a video frame. Besides, TuSimple is relatively easy while CULane and BDD100K are more challenging considering the total number of video frames and road types. Note that original BDD100K dataset provides 100 k video frames, 70 k for training, 10 k for validation and 20 k for testing. However, since they do not provide ground-truth labels for testing set, we keep the training set unchanged but use the original validation set for testing.

\subsection{Evaluation Criterion}

\noindent \textbf{TuSimple:} in TuSimple, we use the official metric (accuracy) as the evaluation criterion. Besides, false positive ($FP$) and false negative ($FN$) are also listed. The following is the equation to compute accuracy~\cite{tusimple}:

\begin{equation}
Accuracy = \frac{N_{pred}}{N_{gt}} ,
\end{equation}

where $N_{pred}$ is the number of correctly predicted lane points and $N_{gt}$ is the number of ground-truth lane points.

\noindent \textbf{CULane and BDD100K:} to judge whether a lane is correctly detected, we treat each lane as a line with fixed pixel width (30 for CULane and 8 for BDD100K) and compute the intersection-over-union (IoU) between labels and predictions. Predictions whose IoUs are larger than 0.5 are considered as true positives (TP). Then, we use $F_{1}-measure$ as the evaluation metric formulated as follows:

\begin{equation}
F_{1}-measure = \frac{2 \times Precision \times Recall}{Precision + Recall} ,
\end{equation}

where $Precision = \frac{TP}{TP + FP}$ and $Recall = \frac{TP}{TP + FN}$.

\subsection{Methods}

We give a brief introduction to several CNN-based lane detection models. Note that the following approaches all adopt the encoder and decoder structure as the main network.

\noindent \textbf{SCNN:} SCNN uses VGG-16 as the encoder, bilinear upsampling as the decoder and performs message passing to better exploit the long, thin prior of lanes. Besides, the lane existence prediction branch is added to provide extra performance gains as well as facilitate the evaluation process.

\noindent \textbf{ResNet-101:} similar to SCNN except that the encoder is ResNet-101 and no message passing is used, it achieves compelling performance in CULane.

\noindent \textbf{LaneNet:} it is mainly comprised of two sub-networks, i.e., LaneNet and H-Net. The former is used to generate the segmentation maps differentiating lane instances and background. The later is utilized to obtain the lane points from original images and segmentation maps to facilitate evaluation.

\noindent \textbf{EL-GAN:} it adopts U-Net~\cite{ronneberger2015u} as the generator (i.e., the encoder and decoder) and DenseNet~\cite{huang2017densely} as the discriminator. The introduction of GAN aims at generating more realistic segmentation maps as well as easing the complex post-processing.

\noindent \textbf{GC-LD:} it uses VGG-16 as the backbone and performs lane and road detection simultaneously. Both geometric prior and inner relationship between two detection tasks are exploited to boost the overall performance of their models.

\subsection{Preliminary}

Lane detection is commonly formulated as a semantic segmentation task. More specifically, given an input image \textbf{X}, the objective is to assign a label $l_{ij}$ ($l_{ij}$ = 1, ..., $N_{c}$) to each pixel $(i, j)$ of \textbf{X}, comprising the segmentation map $s$. Here, $N_{c}$ is the number of classes. The objective is to learn a mapping $\mathcal{F}$: \textbf{X} $\mapsto$ $s$. Recent studies use CNN as $\mathcal{F}$ for end-to-end prediction. The task of lane existence prediction is also introduced to facilitate the evaluation process. We use $b$ to represent the binary labels that indicate the existence of lanes. Then, the function becomes $\mathcal{F}$: \textbf{X} $\mapsto$ ($s$, $b$).

\subsection{Self Attention Distillation}

Apart from training our lane detection network with the aforementioned semantic segmentation and lane existence prediction losses, we aim to perform layer-wise and top-down attention distillation to enhance the representation learning process. The proposed SAD does not require any external supervision or additional labels since attention maps are derived from the network itself.

In general, attention maps can be divided into two categories, \ie, activation-based attention maps~\cite{zagoruyko2016paying} and gradient-based attention maps~\cite{zagoruyko2016paying}. The activation-based attention maps are obtained via processing the activation output of a specific layer while the gradient-based ones are obtained via using the layer's gradient output. In the experiment, we empirically find that activation-based attention distillation yields considerable performance gains while gradient-based attention distillation barely works. Hence, in the following sections we only discuss the activation-based attention distillation.

\noindent \textbf{Activation-based attention distillation.} We use $A_{m} \in R^{C_{m} \times H_{m} \times W_{m}}$ to denote the activation output of the $m$-th layer of the network, where $C_{m}$, $H_{m}$ and $W_{m}$ denote the channel, height and width, respectively. Let $M$ denote the number of layers in the network. The generation of the attention map is equivalent to finding a mapping function $\mathcal{G}$: $R^{C_{m} \times H_{m} \times W_{m}} \to R^{H_{m} \times W_{m}}$. The absolute value of each element in this map represents the importance of this element on the final output. Therefore, this mapping function can be constructed via computing statistics of these values across the channel dimension. More specifically, the following three operations~\cite{zagoruyko2016paying} can serve as the mapping function: $\mathcal{G}_{sum}(A_{m}) = \sum_{i=1}^{C_{m}}| A_{mi}|$, $\mathcal{G}_{sum}^{p}(A_{m}) = \sum_{i=1}^{C_{m}}| A_{mi}|^{p}$ and $\mathcal{G}_{\max}^{p}(A_{m}) = \max_{i=1, C_{m}}| A_{mi}|^{p}$. Here, $p > 1$ and $A_{mi}$ denotes the $i$-th slice of $A_{m}$ in the channel dimension.

\begin{figure}[b]
  \centering
  \includegraphics[width=1.0\linewidth]{./compare_func.pdf}

  \caption{Attention maps of the block 4 of the ENet model using different mapping functions.}
  \centering
  \vskip -0.2cm
  \label{fig:func}
\end{figure}

\begin{figure*}[t]
  \centering
  \includegraphics[width=0.78\linewidth]{./framework.pdf}
  \vskip -0.2cm
  \caption{An instantiation of using SAD. $E_{1}$ $\sim$ $E_{4}$ comprise the encoder of ENet~\cite{paszke2016enet}, $D_{1}$ and $D_{2}$ comprise the decoder of ENet. Following \cite{pan2017spatial}, we add a small network to predict the existence of lanes, denoted as $P_{1}$. AT-GEN is the attention generator.}
  \centering
  \vskip -0.2cm
  \label{fig:pipeline}
\end{figure*}

The differences between these mapping functions are depicted in Fig.~\ref{fig:func}. Compared with $\mathcal{G}_{sum}(A_{m})$, $\mathcal{G}_{sum}^{p}(A_{m})$ puts more weights to areas with higher activations. The larger the $p$ is, the more focus is placed on these highly activated areas. Compared with $\mathcal{G}_{\max}^{p}(A_{m})$, $\mathcal{G}_{sum}^{p}(A_{m})$ is less biased since it calculates weights across multiple neurons instead of selecting the maximum value of these neuron activations as the weight. In the experiment, we empirically find that using $\mathcal{G}_{sum}^{2}(.)$ as the mapping function yields the most performance gains.

\noindent \textbf{Adding SAD to training.} The intuition behind SAD~is that the attention maps of previous layers can distil useful contextual information from those of successive layers. Following ~\cite{zagoruyko2016paying}, we also perform spatial softmax operation $\Phi(.)$ on $\mathcal{G}_{sum}^{2}(A_{m})$. Bilinear upsampling $\mathcal{B}(.)$ is added before the softmax operation if the size of original attention maps is different from that of targets. However, different from Sergey \etal ~\cite{zagoruyko2016paying} who perform attention distillation within two networks, the proposed self attention distillation is performed within the network itself.

Adding SAD to an existing network is straight-forward. It is possible to introduce SAD at different timepoint of the training, which could affect the convergence time. We will show an evaluation in the experiment section. Here we assume an ENet half-trained to 40K episodes.
As shown in Fig. \ref{fig:pipeline}, we add an attention generator, abbreviated as AT-GEN, after each $E_2$, $E_3$, and $E_4$ encoder block of ENet. Formally, AT-GEN is represented by a function $\Psi(.) = \Phi(\mathcal{B}(\mathcal{G}_{sum}^{2}(.)))$. A successive layer-wise distillation loss is formulated as follows:
\begin{equation}
\label{eqn:sad_loss}
\begin{split}
\mathcal{L}_\mathrm{distill}(A_{m}, A_{m+1}) =
\sum_{m=1}^{M-1} \mathcal{L}_\mathrm{d} (\Psi(A_{m}), \Psi(A_{m+1})),
\end{split}
\end{equation}
where $\mathcal{L}_\mathrm{d}$ is typically defined as a $L_{2}$ loss and $\Psi(A_{m+1})$ is the target of the distillation loss. In the example shown in Fig. \ref{fig:pipeline}, we have the number of layers $M = 4$. Note that we do not assign different weights to different SAD paths, although this is possible. We found that this uniform scheme works well in our experiments.

The total loss is comprised of four terms:
\begin{equation}
\label{eqn:total_loss}
\begin{split}
\mathcal{L} = & \underbrace{\mathcal{L}_\mathrm{seg}(s, \hat{s}) + \alpha \mathcal{L}_\mathrm{IoU}(s, \hat{s})}_{\text{segmentation loss}} \\
& + \underbrace{\beta \mathcal{L}_\mathrm{exist}(b, \hat{b})}_{\text{existence loss}} + \underbrace{\gamma \mathcal{L}_\mathrm{distill}(A_{m}, A_{m+1})}_{\text{distillation loss}}.
\end{split}
\end{equation}
Here, the first two terms are segmentation losses that comprise of the standard cross entropy loss $\mathcal{L}_\mathrm{seg}$ and the IoU loss $\mathcal{L}_\mathrm{IoU}$. The IoU loss aims at increasing the intersection-over-union between the predicted lane pixels and ground-truth lane pixels. It is formulated as $\mathcal{L}_{IoU} = 1 - \frac{N_{p}}{N_{p} + N_{g} - N_{o}}$, where $N_{p}$ is the number of predicted lane pixels, $N_{g}$ is the number of ground-truth lane pixels and $N_{o}$ is the number of lane pixels in the overlapped areas between predicted lane areas and ground-truth lane areas.
$\mathcal{L}_\mathrm{exist}$ is the binary cross entropy loss. $\hat{s}$ is the segmentation map produced by the network and $\hat{b}$ is the prediction of the existence of lanes. The parameters $\alpha$, $\beta$, and $\gamma$ balance the influence of segmentation losses, existence loss, and distillation loss on the final task.

It is noteworthy that the SAD paths can be generalized to dense connections beyond the example shown here. For instance, we can add block~1 $ \xrightarrow[]{\text{mimic}}$ block~3, block~1 $ \xrightarrow[]{\text{mimic}}$ block~4, and block~ 2 $ \xrightarrow[]{\text{mimic}}$ block~4 in addition to the current paths. In general, the number of possible SAD paths for a network with a depth of $M$ layers is $\frac{M(M-1)}{2}$. We will evaluate this possibility in our experiments.

\noindent \textbf{Visualization of attention maps with and without SAD.}
We investigate the influence of SAD by studying the attention maps of different blocks in ENet with and without SAD. More results will be reported in Section~\ref{sec:experiments}. Both networks with and without SAD are trained up to 60K episodes.
We visualize the attention maps of four existing blocks in ENet. As can be observed in Fig.~\ref{fig:learning}, after adding SAD, the attention maps of ENet become more concentrated on task-relevant objects, \eg, lanes, vehicles and road curbs. This would in turn improve the lane detection accuracy, as we will show in the experiments.

\begin{figure*
  \centering
  \includegraphics[width=0.75\linewidth]{./vis_sad_learning.pdf}

  \caption{Attention maps of ENet with and without self attention distillation. Both networks with and without SAD are trained up to 60K episodes. SAD is applied to ENet at 40K training episodes.}
  \centering
  \vskip -0.2cm
  \label{fig:learning}
\end{figure*}

\subsection{Lane Prediction}
\label{subsec:inference}

The output of the model is not post-processed for TuSimple and BDD100K except CULane. For CULane, in the inference stage, we feed the image into the ENet model. Then the multi-channel probability maps and the lane existence vector are obtained. Following~\cite{pan2017spatial}, the final output is obtained as follows: First, we use a 9 $\times$ 9 kernel to smooth the probability maps. Then, for each lane whose existence probability is larger than 0.5, we search the corresponding probability map every 20 rows for the position with the highest probability value. In the end, we use cubic splines to connect these positions to get the final output.

\subsection{Architecture Design}
\label{subsec:architecture}

The original ENet model is an encoder-decoder structure comprised of $E_{1} \sim E_{4}$, $D_{1}$ and $D_{2}$. Following \cite{pan2017spatial}, we add a small network $P_{1}$ to predict the existence of lanes. The encoder module is shared to save memory space.
Apart from this modification, we also observed some useful techniques to modify ENet for achieving better performance in the lane detection task.
Dilated convolution \cite{yu2015multi} is added to replace the original convolution layers in the lane existence prediction branch to increase the receptive field of the network without increasing the number of parameters. In the original design, the resolution of feature maps of $E_{4}$ is only 36 $\times$ 100 for CULane. This leads to severe loss of information. Hence, we use feature concatenation to fuse the output of $E_{4}$ with that of $E_{3}$ so that the output of the encoder can benefit from information encoded in previous layers.


\section{Introduction}
Recently, with the detection of a gravitational wave (GW) event \cite{LIGOScientific:2016,{LIGOScientific:2017}} 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 \cite{Vilenkin:2000}, inflationary quantum fluctuations \cite{Starobinsky:1979ty}, preheating \cite{Khlebnikov:1997di}, first-order phase transitions in the early Universe \cite{Witten:1984rs, {Hogan:1986qda}}, 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 \cite{Kuzmin:1985}-\cite{Morrissey:2012}. 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 \cite{Kosowsky:1992a}-\cite{Caprini:2009}. The GWs are produced from the strong first-order electroweak phase transition mainly by the following three mechanisms: bubble collisions \cite{Kosowsky:1992a}-\cite{Caprini:2008}, sound waves in the hot plasma \cite{Hindmarsh:2014}-\cite{Hindmarsh:2015} and magnetohydrodynamic turbulence of bubbles \cite{Caprini:2006}-\cite{Caprini:2009} 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 \cite{Kajantie:1996}-\cite{Csikor:1999}. 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 \cite{Land:1992}-\cite{Blinov:2015}. 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 \cite{Kozaczuk:2014kva}-\cite{Huang:2017rzf}.
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 \cite{Ma:2006}-\cite{Lopez:2011} 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 \cite{Banik:2014cfa,{Banik:2014eda}}.
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 \cite{Kannike:2012pe}, perturbativity,  Large Electron-Positron Collider (LEP)\cite{Beringer:2012} and Large Hadron Collider (LHC) bounds, PLANCK bound on the DM relic density \cite{Ade:2013zuv} and the limits given by spin independent direct detection experiments like XENON-1T \cite{Aprile:2015uzo}, PandaX-II  \cite{Tan:2016zwf}, LUX \cite{Akerib:2016vxi} and DarkSide-50 \cite{Marini:2016haq}. 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 \cite{Profumo:2014opa}. 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) \cite{Moore:2015,{Corbin:2006}}, Evolved Laser Interferometer Space Antenna (eLISA) \cite{Moore:2015}, Advanced Laser Interferometer Antenna (ALIA) \cite{Gong:2015}, DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) \cite{Moore:2015,{Musha:2017}}, Ultimate-DECIGO (U-DECIGO) \cite{Kudoh:2006} and ground-based detector Advanced LIGO (aLIGO) \cite{Moore:2015,{aLIGO:2019}}.

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.

\section{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

\begin{equation}\label{eq:1}
\begin{aligned}
V=m_1^2 \Phi_H^\dagger\Phi_H + m_2^2 \Phi_I^\dagger\Phi_I+\dfrac{1}{2}m_s^2 S^2 + \lambda_1 \left(\Phi_H^\dagger\Phi_H \right)^2 + \lambda_2 \left(\Phi_I^\dagger\Phi_I \right)^2 + \lambda_3 \left(\Phi_H^\dagger\Phi_H \right) \left(\Phi_I^\dagger\Phi_I \right) \\+ \lambda_4 \left(\Phi_I^\dagger\Phi_H \right) \left(\Phi_H^\dagger\Phi_I \right) + \dfrac{\lambda_5}{2} \left[\left(\Phi_I^\dagger\Phi_H \right)^2+\left(\Phi_H^\dagger\Phi_I \right)^2 \right] + \rho_1 \left(\Phi_H^\dagger\Phi_H \right) S \hspace{2.7cm} \\+\rho_1^{\prime} \left(\Phi_I^\dagger\Phi_I \right) S + \rho_2 \left(\Phi_H^\dagger\Phi_H \right) S^2+ \rho_2^{\prime} \left(\Phi_I^\dagger\Phi_I \right) S^2 + \dfrac{\rho_3}{3} S^3+ \dfrac{\rho_4}{4} S^4. \hspace{3.5cm}
\end{aligned}
\end{equation}
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
\begin{equation}\label{eq:2}
\phi_H=\left(
\begin{array}{c}
 G^+\\
 \dfrac{1}{\sqrt{2}} \left(v+h+iG_0\right)\\
\end{array}
\right), \hspace{1mm}\phi_I=\left(
\begin{array}{c}
 H^+ \\
 \dfrac{1}{\sqrt{2}} \left(H_0+iA_0 \right)\\
\end{array}
\right), \hspace{1mm}S=v_s+s,
\end{equation}
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$) \cite{Banik:2014eda}.
After minimising the scalar potential represented in Eq. (\ref{eq:1}) using the conditions
\begin{equation}\label{eq:3}
\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,
\end{equation}
we obtain
\begin{equation}\label{eq:4}
\begin{aligned}
m_1^2+\lambda_1v^2+\rho_1v_s+\rho_2v_s^2=0,\\
m_s^2+\rho_3v_s+\rho_4v_s^2+\dfrac{\rho_1v^2}{2v_s}+\rho_2v^2=0.
\end{aligned}
\end{equation}
The mass matrix of the scalar sector is obtained by calculating the second-order derivatives $\left(\dfrac{\partial ^2V}{\partial {h}^2},\hspace{1mm} \dfrac{\partial ^2V}{\partial {s}^2}, \hspace{1mm}\dfrac{\partial ^2V}{\partial {H^{\pm}}^2},\hspace{1mm} \dfrac{\partial ^2V}{\partial {H_0}^2},\hspace{1mm} \dfrac{\partial ^2V}{\partial {A_0}^2},\hspace{1mm} \dfrac{\partial ^2V}{\partial {h}\partial {s}}\right)$ of the scalar potential (Eq. (\ref{eq:1})) and the elements are

\begin{equation}\label{eq:5}
\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,
\end{equation}

\begin{equation}\label{eq:6}
\mu_s^2=m_s^2+\rho_2 v^2+2\rho_3 v_s+3\rho_4 v_s^2=\rho_3v_s+2\rho_4v_s^2-\dfrac{\rho_1v^2}{2v_s},
\end{equation}

\begin{equation}\label{eq:8}
m_{H^\pm}^2=m_2^2+\dfrac{\lambda_3v^2}{2}+\rho_1^\prime v_s+\rho_2^\prime v_s^2,
\end{equation}

\begin{equation}\label{eq:9}
m_{H_0}^2=m_2^2+\left(\lambda_3+\lambda_4+\lambda_5\right)\dfrac{v^2}{2}+\rho_1^\prime v_s+\rho_2^\prime v_s^2,
\end{equation}

\begin{equation}\label{eq:10}
m_{A_0}^2=m_2^2+\left(\lambda_3+\lambda_4-\lambda_5\right)\dfrac{v^2}{2}+\rho_1^\prime v_s+\rho_2^\prime v_s^2,
\end{equation}

\begin{equation}\label{eq:7}
\mu_{hs}^2=\rho_1v+2\rho_2v_sv.
\end{equation}
As $h$ and $s$ mix, we diagonalise the mass matrix in $h$, $s$ basis by a unitary matrix $U$
\begin{equation}\label{eq:11}
\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),
\end{equation}
to obtain two physical mass eigenstates $h_1$ and $h_2$ as
\begin{equation}\label{eq:12}
h_1=h\hspace{1mm}\cos\theta-s\hspace{1mm}\sin\theta, \hspace{1mm}h_2 =h\hspace{1mm}\sin\theta +s\hspace{1mm}\cos\theta,
\end{equation}
where $\theta$ is the mixing angle that can be computed from
\begin{equation}\label{eq:13}
\tan\theta=\dfrac{y}{1+\sqrt{1+y^2}}, \hspace{2mm} \text{where} \hspace{1mm} y=\dfrac{2 \mu_{hs}^2}{\mu_{h}^2-\mu_s^2} \hspace{1mm}.
\end{equation}
The expressions for the mass eigenstate of the two physical scalars $h_1$ and $h_2$ are given as
\begin{equation}\label{eq:14}
m_{h_1, h_2}^2=\dfrac{\left(\mu_{h}^2+\mu_s^2\right)}{2}\pm \dfrac{\left(\mu_{h}^2-\mu_s^2\right)}{2} \sqrt{1+y^2},
\end{equation}
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 \cite{Patrignani:2016xqp} 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. (\ref{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. (\ref{eq:4})-(\ref{eq:14}) we obtain the following relations
\begin{equation}\label{eq:15}
\lambda_1=\dfrac{m_{h_2}^2\sin ^2\theta+m_{h_1}^2 \cos ^2\theta}{2 v^2},
\end{equation}

\begin{equation}\label{eq:16}
\rho_4=\dfrac{1}{2v_{s}^2}\left(m_{h_2}^2\cos ^2\theta+m_{h_1}^2 \sin ^2\theta + \dfrac{v^2\rho_1}{2v_s}-v_s\rho_3\right),
\end{equation}

\begin{equation}\label{eq:17}
\rho_2=\dfrac{m_{h_1}^2-m_{h_2}^2}{4 v\hspace{.1mm} v_{s}}\sin {2\theta}-\dfrac{\rho_1}{2v_s},
\end{equation}

\begin{equation}\label{eq:18}
m_1^2=-\lambda_1v^2-\rho_1v_s-\rho_2v_s^2,
\end{equation}

\begin{equation}\label{eq:19}
m_2^2=m_{H_0}^2-\dfrac{v^2}{2}\left(\lambda_3+\lambda_4+\lambda_5\right)-\rho_1^\prime v_s -\rho_2^\prime v_s^2,
\end{equation}

\begin{equation}\label{eq:20}
m_s^2=-\rho_3v_s-\rho_4v_s^2-\dfrac{\rho_1 v^2}{2v_s}-\rho_2 v^2.
\end{equation}

\section{Constraints}
In this section various theoretical and experimental bounds are given. These are used to constrain the model parameter space.
\subsection{Theoretical Constraints}
\noindent \underline{\it Vacuum Stability}

From the vacuum stability conditions the bounds on the couplings are given as \cite{Kannike:2012pe}
\begin{equation}\label{eq:21}
\lambda_1, \lambda_2, \rho_4>0, \hspace{1mm} \lambda_3+2\sqrt{\lambda_1 \lambda_2}>0,
\end{equation}

\begin{equation}\label{eq:22}
\lambda_3+\lambda_4-|\lambda_5|+2\sqrt{\lambda_1 \lambda_2}>0, \hspace{1mm}\rho_2+\sqrt{\lambda_1 \rho_4}>0,
\end{equation}

\begin{equation}\label{eq:23}
\rho_2^\prime+\sqrt{\lambda_2 \rho_4}>0,
\end{equation}

\begin{equation}\label{eq:24}
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,
\end{equation}

\begin{equation}\label{eq:25}
\begin{aligned}
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} \hspace{3cm}\\+\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. \hspace{1cm}
\end{aligned}
\end{equation}

\noindent \underline{\it Perturbativity}

All the quartic couplings in the tree-level potential (Eq. (\ref{eq:1})) must be less than $4\pi$ to be consistent with the perturbative conditions.
\subsection{Experimental Constraints}
\noindent \underline{\it Collider Constraints}

From the LEP experiment, the bound on the model parameter space is given as  \cite{Beringer:2012}
\begin{equation}\label{eq:26}
m_{H_0}+m_{A_0}>m_Z, \hspace{1mm} m_{H^\pm}>79.3 \hspace{1mm} {\text{GeV}}.
\end{equation}

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
\begin{equation}\label{eq:27}
R_1=\cos ^4\theta\dfrac{\Gamma ^{\text{SM}}}{\Gamma },
\end{equation}
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
\begin{equation}\label{eq:28}
\Gamma =\cos ^2\theta \hspace{1mm} \Gamma ^{\text{SM}}+\Gamma ^{\text{inv}},
\end{equation}
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
 \begin{equation}\label{eq:29}
\Gamma^{inv}\left(h_1\rightarrow H_0 H_0\right)=\dfrac{\left(g_{h_1 H_0 H_0}\right)^2}{16\pi m_{h_1}}\left(1-\dfrac{4m_{H_0}^2}{m_{h_1}^2}\right)^{1/2},
\end{equation}
and
\begin{equation}\label{eq:30}
\Gamma^{inv}\left(h_1\rightarrow h_2 h_2\right)=\dfrac{\left(g_{h_1 h_2 h_2}\right)^2}{16\pi m_{h_1}}\left(1-\dfrac{4m_{h_2}^2}{m_{h_1}^2}\right)^{1/2},
\end{equation}
where the expressions for the couplings $g_{h_1 H_0 H_0}$ and $g_{h_1 h_2 h_2}$ are given as
\begin{equation}\label{eq:69}
g_{h_1 H_0 H_0} = \left( \dfrac{\lambda_L}{2} \cos\theta - \dfrac{\lambda_s}{2} \sin\theta \right)v=\lambda_{h_1 H_0 H_0} v,
\end{equation}
with
\begin{equation}\label{eq:40}
\lambda_L=\lambda_3+\lambda_4+\lambda_5 \hspace{1mm}, \hspace{1mm} \lambda_s=\dfrac{\rho_1^\prime+2\rho_2^\prime v_s}{v},
\end{equation}
and
\begin{equation}\label{eq:71}
\begin{aligned}
g_{h_1 h_2 h_2} = \dfrac{1}{2}(6\lambda_1 v\cos\theta \sin^2\theta + 2\rho_1 \cos^2\theta \sin\theta -\rho_1 \sin^3\theta + 2v\rho_2 \cos^3\theta -4v\rho_2 \cos\theta\sin^2\theta \\+4v_s\rho_2 \cos^2\theta\sin\theta - 2v_s \rho_2 \sin^3\theta - 2\rho_3 \cos^2\theta\sin\theta -6v_s \rho_4 \cos^2\theta\sin\theta). \hspace{1cm}
\end{aligned}
\end{equation}

The invisible decay branching fraction for SM like scalar can be expressed as
\begin{equation}\label{eq:31}
\text{Br}_{\text{inv}}=\dfrac{\Gamma^{\text{inv}}}{\Gamma}.
\end{equation}
We adopt the bounds on the invisible decay branching fraction for SM scalar to be $\text{Br}_{\text{inv}}\leq24\%$ \cite{CERN:2016} (for $m_{h_1}\ge m_{H_0}/2$).
The limit on the  scalar mixing $\sin\theta\le0.4$ \cite{Robens:2015}-\cite{Dupuis:2016} and the signal strength of SM Higgs $R_1\ge0.84$ \cite{Khachatryan:2015, {Aad:2016}}
is obtained from LHC experimental results. However, study of future collider experiments can
further constrain  the  scalar mixing angle \cite{Profumo:2014opa}. From the study of Ref.~\cite{Profumo:2014opa}
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\sim0.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.

\noindent \underline{\it 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
\begin{equation}\label{eq:32}
0.1172\leq \Omega _{\text{DM}}{\rm h}^2\leq 0.1226.
\end{equation}
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}$ \cite{Ade:2013zuv}.

\noindent \underline{\it 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
\cite{Aprile:2015uzo}, PandaX-II  \cite{Tan:2016zwf}, LUX \cite{Akerib:2016vxi} and DarkSide-50 \cite{Marini:2016haq}.

\section{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.
\subsection{Relic Density}
In order to calculate the dark matter relic density one needs to solve the Boltzmann equation which can be expressed as  \cite{Kolb:1990vq}
\begin{equation}\label{eq:33}
\frac{dn_{H_0 }}{dt}+3Hn_{H_0}=-\langle\sigma v\rangle\left[n_{H_0 }^2-\left(n_{H_0 }^{\text{eq}}\right){}^2\right],
\end{equation}
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. (\ref{eq:33}), $\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
\begin{equation}\label{eq:34}
\left<\sigma v\right>=\dfrac{1}{8m_{H_0}^4 TK_2^2\left(\frac{m_{H_0 }}{T}\right)}\int _{4 m_{H_0 }^2}^{\infty }ds\left(s-4 m_{H_0 }^2\right)\sqrt{s}K_1\left(\frac{\sqrt{s}}{T}\right) \sigma (s),
\end{equation}
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
\begin{equation}\label{eq:35}
\Omega_{\rm DM} {\rm h}^2 =2.755\times 10^8
\left(\dfrac{m_{H_0}}{\rm{GeV}}\right) Y_0,
\end{equation}
with
\begin{equation}\label{eq:36}
\dfrac{1}{Y_0}=\dfrac{1}{Y_F}+\left(\dfrac{45 G}{\pi}\right)^{-\frac{1}{2}} \int_{T_0}^{T_F} g_{*}^{1/2} \langle\sigma {\rm v}\rangle dT,
\end{equation}
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 \cite{Alloul:2013bka} where the extended IDM is implemented and finally micrOMEGAs \cite{Belanger:2013oya} (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^-,\hspace{1mm} ZZ,\hspace{1mm} f\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,\hspace{1mm}h_2 h_2,\hspace{1mm} h_1 h_2$ via s, t and u channels as also the relevant four-point interactions have been considered. The micrOMEGAs \cite{Belanger:2013oya} has been used to estimate the relic-density for dark matter candidate $H_0$ considering all the annihilation channels mentioned above.

\subsection{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
\begin{equation}\label{eq:37}
\sigma_{\text{SI}}=\dfrac{m_N^4}{\pi \left( m_{H_0}+m_N\right)^2} f^2 \left(\dfrac{\lambda_{h_1 H_0 H_0}\cos\theta}{m_{h_1}^2}+\dfrac{\lambda_{h_2 H_0 H_0}\sin\theta}{m_{h_2}^2}\right)^2,
\end{equation}
where $m_N$ is the nucleon mass and $f$ is the nucleon form factor which has been approximated as 0.3 \cite{cline:2015}. In Eq. (\ref{eq:37}) the couplings $\lambda_{h_1 H_0 H_0}$ and $\lambda_{h_2 H_0 H_0}$ can be written as

\begin{equation}\label{eq:38}
\lambda_{h_1 H_0 H_0} = \left( \dfrac{\lambda_L}{2} \cos\theta - \dfrac{\lambda_s}{2} \sin\theta \right),
\end{equation}

\begin{equation}\label{eq:39}
\lambda_{h_2 H_0 H_0} = \left( \dfrac{\lambda_L}{2} \sin\theta + \dfrac{\lambda_s}{2} \cos\theta \right).
\end{equation}

We calculate the dark matter scattering cross-sections using Eqs. (\ref{eq:37})-(\ref{eq:39}) and then check whether they satisfy the direct detection cross-section bounds obtained from DM direct detection experiments such as XENON-1T
\cite{Aprile:2015uzo}, PandaX-II  \cite{Tan:2016zwf}, LUX \cite{Akerib:2016vxi} and DarkSide-50 \cite{Marini:2016haq}.
\subsection{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. (\ref{eq:32})).
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
\cite{Aprile:2015uzo}, PandaX-II  \cite{Tan:2016zwf}, LUX \cite{Akerib:2016vxi} and DarkSide-50 \cite{Marini:2016haq}.

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 \ref{t1}. 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 \ref{fig:2}-\ref{fig:4}. As mentioned, the relic densities are computed using micrOMEGAs code \cite{Belanger:2013oya}. 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.
\begin{table}[H]
\centering
\footnotesize
\begin{tabular}{|l|c|c|c|c|c|c|c|c|c|c|c|r|}
\hline
BP&$m_{H_{0}}$&$m_{h_{2}}$&$v_{s}$&$\sin\theta$&$\rho_{1}$&$\rho_{3}$&$\lambda_{L}$&$\lambda_{s}$&$\lambda_{2}$&$\Omega_{\text{DM}}{\rm h}^2$&$\sigma_{{\rm{SI}}}$\\
 &in GeV&in GeV&in GeV&&in GeV&in GeV&&&&&$\rm{cm^2}$\\
 \hline
1&30&100&300&0.01&-3&0.01&0.001&0.0012&0.2&0.1220&9.41$\times 10^{-48}$\\
 \hline
2&68&150&400&0.06&-7&0.2&0.002&0.033&0.031&0.1208&3.69$\times 10^{-48}$\\
 \hline
3&62&250&400&0.01&-5&0.5&0.0003&0.0033&0.1&0.1206&1.70$\times 10^{-49}$\\
\hline
4&76&200&500&0.03&-1&0.1&0.0016&0.0033&0.01&0.1195&3.56$\times 10^{-48}$\\
\hline
\end{tabular}
\caption{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.}\label{t1}
\end{table}
The dips in the relic density plots (left panels of Figures \ref{fig:2}-\ref{fig: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 \ref{fig:2}-\ref{fig: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 \ref{fig:2}-\ref{fig: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. (\ref{eq:37}). The right panels of Figures \ref{fig:2}-\ref{fig: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.
\begin{figure}[H]
\centering
\includegraphics[width=7cm,height=5.3cm]{relicfbp1.png}
\includegraphics[width=7cm,height=5.3cm]{directbp1.png}
 \caption{
Variation of DM relic density (left panel) and variation of DM direct detection scattering cross-section (right panel) as a function of DM mass $m_{H_0}$ for BP1. The shadowed area in the left panel bounded by two horizontal parallel lines shows the PLANCK bounds and the red line represents the calculated relic densities. DM direct detection scattering cross-section (right panel), calculated for different DM masses $m_{H_0}$, is overplotted with the bounds given by dark matter direct detection experiments such as XENON-1T, LUX, PandaX-II and DarkSide-50.}

\label{fig:2}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=8cm,height=7cm]{relicfbp2new.png}
\includegraphics[width=8cm,height=7cm]{directbp2new.png}
 \caption{Same as in Figure \ref{fig:2} but for BP2.}
\label{fig:3}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=8cm,height=7cm]{relicfbp3.png}
\includegraphics[width=8cm,height=7cm]{directbp3.png}
 \caption{Same as in Figure \ref{fig:2} but for BP3.}
\label{fig:3new}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=8cm,height=7cm]{relicfbp4new.png}
\includegraphics[width=8cm,height=7cm]{directbp4new.png}
 \caption{Same as in Figure \ref{fig:2} but for BP4.}
\label{fig:4}
\end{figure}

\section{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.
\subsection{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. (\ref{eq:1})). The finite temperature effective potential can be written as \cite{Wainwright:2012}
\begin{equation}\label{eq:41}
V_{\text{eff}}=V_{\text{}tree-level}+V_{1-loop}^{T=0}+V_{1-loop}^{T\neq0},
\end{equation}
where $V_{1-loop}^{T=0}$ and $V_{1-loop}^{T\neq0}$ 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 \cite{Wainwright:2012}
\begin{equation}\label{eq:42}
V_{1-loop}^{T=0}=\pm \dfrac{1}{64\pi^2} \sum_{i} n_i m_i^4 \left[ \log\dfrac{m_i^2}{Q^2}-C_i\right],
\end{equation}
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
\begin{equation}\label{eq:5b}
m_W^2=\dfrac{1}{4}g^2\left(v^2+v_\phi^2\right)\,\,,
\end{equation}

\begin{equation}\label{eq:5c}
m_Z^2=\dfrac{1}{4}\left(g^2+{g^{\prime}}^2\right)\left(v^2+v_\phi^2\right)\,\,,
\end{equation}

\begin{equation}\label{eq:5a}
\mu_{h}^2=m_1^2+3\lambda_1 v^2+ \dfrac{v_\phi^2}{2}\lambda_L+\rho_1 v_s+\rho_2 v_s^2\,\,,
\end{equation}

\begin{equation}\label{eq:6a}
\mu_s^2=m_s^2+\rho_2 v^2+2\rho_3 v_s+3\rho_4 v_s^2+\rho_2^{\prime}v_\phi^2\,\,,
\end{equation}

\begin{equation}\label{eq:8a}
m_{H^\pm}^2=m_2^2+\lambda_2 v_\phi^2+\dfrac{\lambda_3v^2}{2}+\rho_1^\prime v_s+\rho_2^\prime v_s^2\,\,,
\end{equation}

\begin{equation}\label{eq:9a}
m_{H_0}^2=m_2^2+3\lambda_2 v_\phi^2+\dfrac{v^2}{2}\lambda_L+\rho_1^\prime v_s+\rho_2^\prime v_s^2\,\,,
\end{equation}

\begin{equation}\label{eq:10a}
m_{A_0}^2=m_2^2+\lambda_2 v_\phi^2+\left(\lambda_3+\lambda_4-\lambda_5\right)\dfrac{v^2}{2}+\rho_1^\prime v_s+\rho_2^\prime v_s^2\,\,,
\end{equation}

\begin{equation}\label{eq:8b}
m_{G^\pm}^2=m_1^2+\lambda_1 v^2+\dfrac{\lambda_3 v_\phi^2}{2}+\rho_1 v_s+\rho_2 v_s^2\,\,,
\end{equation}

\begin{equation}\label{eq:5b}
m_{G_0}^2=m_1^2+\lambda_1 v^2+ \left(\lambda_3+\lambda_4-\lambda_5\right)\dfrac{v_\phi^2}{2}+\rho_1 v_s+\rho_2 v_s^2\,\,,
\end{equation}
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. (\ref{eq:5a}) - (\ref{eq:10a}) reduce to the expressions given in Eqs. (\ref{eq:5}) - (\ref{eq:10}).
We work in the Landau gauge ($\xi=0$), where the Goldstone masses become zero at $T=0$ and also there are no ghost contributions \cite{Basler:2016obg}. However, the gauge invariant methods ($\xi\neq0$) \cite{Patel:2011th} 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 \cite{Kainulainen:2019kyp,{Laine:2000rm}}. We varied $Q$ for a benchmark point and presented the sensitivity for different parameters in Sect. 5.3. In Eq. (\ref{eq:42}), $C_i$ represents a numerical constant. For $W, \hspace{1mm} Z$ 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 \cite{Wainwright:2012}
\begin{equation}\label{eq:43}
V_{1-loop}^{T\neq0}=\dfrac{T^2}{2\pi^2} \sum_{i} n_i J_{\pm}\left[\dfrac{m_i^2}{T^2}\right],
\end{equation}
where
\begin{equation}\label{eq:44}
J_{\pm}\left(\dfrac{m_i^2}{T^2}\right)=\pm \int_0^{\infty} dy \hspace{1mm} y^2 \log\left(1\mp e^{-\sqrt{y^2+\dfrac{m_i^2}{T^2}}}\right).
\end{equation}
We include the daisy resummation procedure in $V_{1-loop}^{T\neq0}$ \cite{Arnold} by adding the temperature dependent terms to the boson masses. Now, the corrected thermal masses \cite{Huang:2017rzf}, \cite{Blinov:2015vma}-\cite{Gil:2012ya} 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
\begin{equation}
c_1=\dfrac{6\lambda_1+2\lambda_3+\lambda_4}{12}+\dfrac{3g^2+{g^{\prime}}^2}{16}+\dfrac{1}{6}\rho_2+\dfrac{y_t^2}{4}\,\,,
\end{equation}

\begin{equation}
c_2=\dfrac{6\lambda_2+2\lambda_3+\lambda_4}{12}+\dfrac{3g^2+{g^{\prime}}^2}{16}+\dfrac{1}{6}\rho_2^{\prime}\,\,,
\end{equation}

\begin{equation}
c_3=\dfrac{1}{8}\rho_4+\dfrac{1}{6}\rho_2+\dfrac{1}{6}\rho_2^{\prime}\,\,,
\end{equation}
here $y_t$ denotes the top Yukawa coupling.
In this work we use the CosmoTransitions package \cite{Wainwright:2012} to compute the finite temperature corrections to the tree-level potential.

\subsection{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 \cite{Linde:1983}
\begin{equation}\label{eq:45}
\Gamma=\Gamma_0\left(T \right) e^{-S_3\left(T \right)/T}
\end{equation}
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 \cite{Linde:1983}
\begin{equation}\label{eq:46}
S_{3}=4\pi \int dr \hspace{1mm}r^{2} \left[ \dfrac{1}{2} \left(\partial_{r} \vec{\phi} \right)^2 +V_{eff}\right],
\end{equation}
where
$V_{eff}$ is the effective finite temperature potential represented in Eq. (\ref{eq:41}). 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$ \cite{Wainwright:2012}.

The gravitational waves are produced from the first-order phase transition mainly by the three mechanisms namely bubble collisions \cite{Kosowsky:1992a}-\cite{Caprini:2008}, sound wave \cite{Hindmarsh:2014}-\cite{Hindmarsh:2015} and turbulence in the plasma \cite{Caprini:2006}-\cite{Caprini:2009}. 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 \cite{Kosowsky:1992a}-\cite{Caprini:2009}
\begin{equation}\label{eq:47}
\Omega_{\text{GW}}{\rm h}^2=\Omega_{\text{col}}{\rm h}^2+ \Omega_{\text{SW}}{\rm h}^2+ \Omega_{\text{turb}}{\rm h}^2.
\end{equation}
The component from the bubbles collision $\Omega_{\text{col}}{\rm h}^2$ is given by
\begin{equation}\label{eq:48}
\Omega_{\text{col}}{\rm h}^2=1.67\times 10^{-5} \left(\dfrac{\beta}{H} \right) ^{-2} \dfrac{0.11 v_{w}^3}{0.42+v_{w}^2} \left(\dfrac{\kappa \alpha}{1+\alpha}\right)^2 \left(\dfrac{g_*}{100}\right)^{-\frac{1}{3}}\dfrac{3.8 \left(\dfrac{f}{f_{col}}\right)^{2.8}}{1+2.8 \left(\dfrac{f}{f_{\text{col}}}\right)^{3.8}}\,\,\,,
\end{equation}
with the parameter $\beta$
\begin{equation}\label{eq:49}
\beta=\left[H T \dfrac{d}{dT}\left( \dfrac{S_3}{T}\right) \right]\bigg|_{T_n},
\end{equation}
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 \cite{Kamionkowski:1994, {Chao:2017vrq},{Dev:2019njv},{Steinhardt:1982}}
\begin{equation}\label{eq:50}
v_w=\dfrac{1/\sqrt{3}+\sqrt{\alpha^2+2\alpha/3}}{1+\alpha}.
\end{equation}
Here in this work, we use the above mentioned expression (Eq. (\ref{eq:50})) for calculating the bubble wall velocity. However, in some literatures, $v_w$ is taken to be 1 for simplicity (\cite{Shajiee:2018jdq, {Mohamadnejad:2019vzg}, {Chala:2019rfk}}). A more detailed discussion on the choice of $v_w$ can be found in \cite{Kozaczuk:2015owa} which results in $v_w$ $\gtrsim$ 0.2. The parameter $\kappa$ in Eq. (\ref{eq:48}) represents the fraction of latent heat deposited in a thin shell which is given by
\begin{equation}\label{eq:51}
\kappa=1-\dfrac{\alpha_{\infty}}{\alpha},
\end{equation}
with \cite{Shajiee:2018jdq,{Caprini:2015zlo}}
\begin{equation}\label{eq:52}
\alpha_{\infty}=\dfrac{30}{24\pi^2 g_{*}} \left(\dfrac{v_n}{T_n} \right)^2 \left[6 \left( \dfrac{m_W}{v}\right)^2 +3\left( \dfrac{m_Z}{v}\right)^2 +6\left( \dfrac{m_{t}}{v}\right)^2\right].
\end{equation}
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
\begin{equation}\label{eq:53}
\alpha=\left[\dfrac{\rho_{\text{vac}}}{\rho^*_{\text{rad}}}\right]\bigg|_{T_n}.
\end{equation}
with
\begin{equation}\label{eq:54}
\rho_{\text{vac}}=\left[\left(V_{\text{eff}}^{\text{high}}-T\dfrac{dV_{\text{eff}}^{\text{high}}}{dT} \right)-\left(V_{\text{eff}}^{\text{low}}-T\dfrac{dV_{\text{eff}}^{\text{low}}}{dT} \right)\right],
\end{equation}
and
\begin{equation}\label{eq:55}
\rho^*_{\text{rad}}=\dfrac{g_* \pi^2 T_n^4}{30}.
\end{equation}
The parameter $\rho_{vac}$ can also be calculated using the trace anomaly \cite{Hindmarsh:2015}, where the same requires an additional factor of $1/4$ in front of the term $T\dfrac{dV}{dT}$ (Eq. (\ref{eq:54})).
The quantity $f_\text{col}$ in Eq. (\ref{eq:48}) is the peak frequency produced by the bubble collisions which takes the form
\begin{equation}\label{eq:56}
f_{\text{col}}=16.5\times10^{-6}\hspace{1mm} \text{Hz} \left( \dfrac{0.62}{v_{w}^2-0.1 v_w+1.8}\right)\left(\dfrac{\beta}{H} \right) \left(\dfrac{T_n}{100 \hspace{1mm} \text{GeV}} \right) \left(\dfrac{g_*}{100}\right)^{\frac{1}{6}}.
\end{equation}
The sound wave (SW) component of the gravitational wave (Eq. (\ref{eq:47})) is given by
\begin{equation}\label{eq:57}
\Omega_{\text{SW}}{\rm h}^2=2.65\times 10^{-6} \left(\dfrac{\beta}{H} \right) ^{-1} v_{w} \left(\dfrac{\kappa_{v} \alpha}{1+\alpha}\right)^2 \left(\dfrac{g_*}{100}\right)^{-\frac{1}{3}}\left(\dfrac{f}{f_{\text{SW}}}\right)^{3} \left[\dfrac{7}{4+3 \left(\dfrac{f}{f_{\text{SW}}}\right)^{2}}\right]^{\frac{7}{2}},
\end{equation}
where $\kappa_v$ is the faction of latent heat transformed into the bulk motion of the fluid which has the following form
\begin{equation}\label{eq:58}
\kappa_v=\dfrac{\alpha_{\infty}}{\alpha}\left[ \dfrac{\alpha_{\infty}}{0.73+0.083\sqrt{\alpha_{\infty}}+\alpha_{\infty}}\right].
\end{equation}
In Eq. (\ref{eq:57}) $f_{\text{SW}}$ denotes the peak frequency produced by the sound wave mechanisms which takes the form

\begin{equation}\label{eq:59}
f_{\text{SW}}=1.9\times10^{-5}\hspace{1mm} \text{Hz} \left( \dfrac{1}{v_{w}}\right)\left(\dfrac{\beta}{H} \right) \left(\dfrac{T_n}{100 \hspace{1mm} \text{GeV}} \right) \left(\dfrac{g_*}{100}\right)^{\frac{1}{6}}.
\end{equation}
The component from the turbulence in the plasma $\Omega_{\text{turb}}{\rm h}^2$ is given by
\begin{equation}\label{eq:60}
\Omega_{\text{turb}}{\rm h}^2=3.35\times 10^{-4} \left(\dfrac{\beta}{H} \right) ^{-1} v_{w} \left(\dfrac{\epsilon \kappa_v \alpha}{1+\alpha}\right)^{\frac{3}{2}} \left(\dfrac{g_*}{100}\right)^{-\frac{1}{3}} \dfrac{\left(\dfrac{f}{f_{\text{turb}}}\right)^{3}\left( 1+\dfrac{f}{f_{\text{turb}}}\right)^{-\frac{11}{3}}}{\left(1+\dfrac{8\pi f}{h_{*}}\right)},
\end{equation}
where $\epsilon=0.1$ and $f_{\text{turb}}$ denotes the peak frequency produced by the turbulence mechanism which can be written as
\begin{equation}\label{eq:61}
f_{\text{turb}}=2.7\times10^{-5}\hspace{1mm} \text{Hz} \left( \dfrac{1}{v_{w}}\right)\left(\dfrac{\beta}{H} \right) \left(\dfrac{T_n}{100 \hspace{1mm} \text{GeV}} \right) \left(\dfrac{g_*}{100}\right)^{\frac{1}{6}}.
\end{equation}
In Eq. (\ref{eq:60}) the parameter $h_*$ has the following form
\begin{equation}\label{eq:62}
h_{*}=16.5\times10^{-6}\hspace{1mm} \text{Hz} \left(\dfrac{T_n}{100 \hspace{1mm} \text{GeV}} \right) \left(\dfrac{g_*}{100}\right)^{\frac{1}{6}}.
\end{equation}

We repeat our analysis by considering an updated expression of $\kappa_v$ \cite{Ellis:2018mja,{Espinosa:2010hh}}, 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 \cite{Ellis:2018mja,{Espinosa:2010hh}}

\begin{equation}\label{eq:62a}
\kappa_v\simeq\left[ \dfrac{\alpha_{\infty}}{0.73+0.083\sqrt{\alpha_{\infty}}+\alpha_{\infty}}\right].
\end{equation}
In this work, Eqs. (\ref{eq:47})-(\ref{eq:62a}) are used for the computation of gravitational wave intensity.

\subsection{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. (\ref{eq:41})) is first computed. For these calculations we use a publicly available package namely the Cosmotransition package \cite{Wainwright:2012}. The tree-level potential (Eq. (\ref{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. (\ref{eq:47})-(\ref{eq:62a}). In this work, we have selected four BPs (Table \ref{t1}) 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 \ref{t1}.

In Figure \ref{fig:5} we show the phase transition properties for BP1. The left and right panels of Figure \ref{fig: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 \ref{fig: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, \hspace{1mm}T_c,\hspace{1mm}T_n,\hspace{1mm}\alpha,\hspace{1mm}\beta/H$) corresponding to each of the four BPs to be uesd for the calculations of GW intensity are furnished in Table \ref{t2}.
\begin{table}[H]
\centering
\small
\begin{tabular}{|l|c|c|c|c|c|c|c|c|r|}
\hline
BP&$v_n$&$T_c$&$T_{n}$&$\alpha$&$\dfrac{\beta}{H}$\\
 &in GeV&in GeV&in GeV&&\\
 \hline
1&226.89&135.68&119.86&0.24&317.86\\
\hline
2&191.03&146.89&132.14&0.25&402.89\\
\hline
3&14.22&130.78&116.99&0.23&692.01\\
\hline
4&209.95&170.92&158.24&0.19&783.65\\
\hline
\end{tabular}
\caption{The values for the parameters used to calculate the GW intensity for each of the chosen BPs. See text for details.}\label{t2}
\end{table}
From Table \ref{t2}, 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 \ref{t2} (corresponding to the BPs given in Table \ref{t1}) the GW intensities and frequencies from first-order phase transitions are now calculated from Eqs. (\ref{eq:47})-(\ref{eq:62a}) for each of the sets. In order to test whether the sound wave contribution to the overall GW signal is significant, the suppression factor ($H R_*/\bar{U}_f$, where $R_*$ is the mean bubble separation and $\bar{U}_f$ is the root-mean-
square (RMS) fluid velocity) is calculated \cite{Caprini:2015zlo,{Ellis:2018mja},{Ellis:2019oqb}}. 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 \ref{fig:6}, where the GW intensities for different BPs are calculated using the expressions of $\kappa_v$ following Eq. \ref{eq:58} and \ref{eq:62a} 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 \ref{fig:6} for comparison. One of the ways to graphically represent the sensitivity of a GW detector is called the power-law-integrated sensitivity curve \cite{Thrane:2013oya} 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. \cite{Moore:2015}. 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 \ref{fig: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 \ref{fig: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 \ref{fig:6} and Table \ref{t2} 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. (\ref{eq:62a}) over Eq. (\ref{eq:58}) introduces slight increase in the GW intensity shown in the right panel of Figure \ref{fig: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.
\begin{figure}[H]
\centering
\includegraphics[width=8cm,height=5.5cm]{bp1transition.pdf}
\includegraphics[width=8cm,height=5.5cm]{tunnelingprofilebp1.png}
\caption{Phase transition properties for the BP1. The left panel is for the position of the minima as a function of temperature and the right panel is for the tunneling profile as a function of bubble radius.}
\label{fig:5}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=8.55cm,height=8cm]{graviintenfnew.png}
\includegraphics[width=8.55cm,height=8cm]{graviintenfneweq.png}
 \caption{Comparison of the chosen four BPs with the sensitivity curves of N1A1M2L4, N2A2M5L4, N2A1M5L6 and N2A5M5L6 configurations of eLISA, ALIA, BBO, DECIGO, U-DECIGO, aLIGO and LISA detectors. The left and right panels show the change in the GW intensity vs frequency using the expression for $\kappa_v$ following Eq. (\ref{eq:58}) and Eq. (\ref{eq:62a}) respectively.}
\label{fig:6}
\end{figure}

\section{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. (\ref{eq:49})) 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.

\vspace{5mm}
\noindent {\bf Acknowledgments}

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