aainstitutetext: School Of Physical Sciences, Indian Association for the
Cultivation of Science, 2A and 2B, Raja S.C. Mullick Road, Kolkata 700032,
Indiabbinstitutetext: Harish-Chandra Research Institute, A CI of Homi Bhabha
National Institute, Chhatnag Road, Jhunsi, Prayagraj (Allahabad) 211019,
Indiaccinstitutetext: Regional Centre for Accelerator-based Particle Physics,
Harish-Chandra Research Institute,
Prayagraj (Allahabad) 211019, India

# Interplay among gravitational waves, dark matter and collider signals in the
singlet scalar extended type-II seesaw model

Purusottam Ghosh b,c Tathagata Ghosh b,c and Subhojit Roy
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study the prospect of simultaneous explanation of tiny neutrino masses,
dark matter (DM), and the observed baryon asymmetry of the Universe in a
$Z_{3}$-symmetric complex singlet scalar extended type-II seesaw model. The
complex singlet scalar plays the role of DM. Analyzing the thermal history of
the model, we identify the region of the parameter space that can generate a
first-order electroweak phase transition (FOEWPT) in the early Universe, and
the resulting stochastic gravitational waves (GW) can be detected at future
space/ground-based GW experiments. First, we find that light triplet scalars
do favor an FOEWPT. In our study, we choose the type-II seesaw part of the
parameter space in such a way that light triplet scalars, especially the
doubly charged ones, evade the strong bounds from their canonical searches at
the Large Hadron Collider (LHC). However, the relevant part of the parameter
space, where FOEWPT can happen only due to strong SM doublet-triplet
interactions, is in tension with the SM-like Higgs decay to a pair of photons,
which has already excluded the bulk of this parameter space. On the other
hand, the latest spin-independent DM direct detection constraints from
XENON-1T and PANDA-4T eliminate a significant amount of parameter space
relevant for the dark sector assisted FOEWPT scenarios, and it is only
possible when the complex scalar DM is significantly underabundant. In short,
we conclude from our analysis that the absence of new physics at the HL-LHC
and/or various DM experiments in the near future will severely limit the
prospects of detecting a stochastic GW at future GW experiments and will
exclude the possibility of electroweak baryogenesis within this model.

###### Keywords:

Beyond Standard Model, Cosmology of Theories beyond the Standard Model,
Electroweak Phase transition, Gravitational wave, Dark Matter

††preprint: HRI-RECAPP-2022–013

##  1 Introduction

The discovery of the Higgs Boson at the Large Hadron Collider (LHC) around 125
GeV ATLAS:2012yve ; CMS:2012qbp  was a breakthrough moment in our
understanding of the laws of nature, completing the particle content of the
Standard Model (SM). Despite the unparalleled successes of the SM, it cannot
be a complete theory of nature as it fails to explain several phenomena in
nature. Among many other issues, the SM cannot explain either the observed
baryon asymmetry of the Universe (BAU) or tiny neutrino masses, in addition to
not providing a good dark matter (DM) candidate.

In order to generate baryon asymmetry in any model, one needs to satisfy the
three well-known Sakharov conditions Sakharov:1967dj . One of the conditions,
i.e., the out-of-equilibrium condition can be achieved if the electroweak
phase transition (EWPT) in the SM had been a strong first-order one. However,
after the discovery of the 125 GeV Higgs, one can conclude that it is only a
smooth crossover. Therefore, the SM cannot satisfy all the Sakharov’s
conditions. Going beyond the SM (BSM), with the addition of new scalars, one
can modify the scalar potential so that a first-order electroweak phase
transition (FOEWPT) can occur in the early Universe. Then one can generate the
observed BAU via the electroweak baryogenesis (EWBG) mechanism Trodden:1998ym
; Anderson:1991zb ; Huet:1995sh ; Morrissey:2012db .

On the other hand, several astrophysical and cosmological data, including
rotation curves of galaxies, the bullet cluster, gravitational lensing, and
the anisotropy of the cosmic microwave background Hu:2001bc , provide strong
motivation for the presence of DM in the Universe. By analyzing the
anisotropies in the CMB data Hu:2001bc , WMAP WMAP:2006bqn  and PLANCK
Planck:2018vyg  experiments one can find that around a quarter of the energy
density of the Universe is composed of non-baryonic, non-luminous matter with
only known interaction through gravity. Apart from its abundance, which is
precisely measured by the PLANCK experiment to be $\Omega_{\rm
DM}h^{2}=0.120\pm 0.001$ Planck:2018vyg , the nature of DM and its non-
gravitational interactions are still unknown to us.

At the same time, the evidence of neutrino oscillations Kajita:2016cak ;
McDonald:2016ixn  has opened up a new frontier in particle physics, requiring
small neutrino masses, and consequently, offering another evidence for
phenomena beyond the SM. Various seesaw mechanisms, such as type-I
Minkowski:1977sc ; Yanagida:1979as ; Mohapatra:1979ia , type-II
Konetschny:1977bn ; Magg:1980ut ; Lazarides:1980nt ; Schechter:1980gr ;
Cheng:1980qt ; Bilenky:1980cx  etc., are the most popular way to generate tiny
neutrino masses. In this article, we focus on type-II seesaw only, where the
SM scalar sector is extended by a $SU(2)$ triplet ($\Delta$) with hypercharge
$Y=2$. Neutrinos obtain masses after the electroweak symmetry breaking (EWSB)
when the neutral component of the triplet acquires an induced vacuum
expectation value (${\it vev}$).

We extend the type-II seesaw model by an extra complex singlet ($S$) and
further impose a discrete $Z_{3}$ symmetry on the model, due to which $S$
becomes a stable DM candidate. The presence of the triplet and the complex
singlet in the model and their interaction with the SM Higgs doublet can
modify the Higgs potential in such a way that a strong first-order electroweak
phase transition (SFOEWPT) is feasible in the early Universe. Hence, by
proposing this economical extension of the SM, we attempt to alleviate three
major shortcomings of it.

The non-observation of DM in its direct search experiments has severely
constrained the simplest weakly interacting massive particle (WIMP)-like DM
scenarios. For instance, DM direct searches impose stringent limits on the
parameter space of the simplest SM Higgs portal scalar DM Casas:2017jjg ;
Bhattacharya:2017fid . However, the situation can drastically change in the
presence of new DM annihilation channels, and one can evade those constraints.
In this model, DM annihilation due to the presence of portal interaction of DM
with the triplet sector and the $Z_{3}-$symmetry induced DM semi-annihilation
channels allow a significant amount of permissible parameter space for DM that
satisfies the observed DM relic density constraint from the PLANCK experiment
Planck:2018vyg , the DM direct detection (DMDD) limits from XENON-1T
Aprile:2018dbl  and PANDAX-4T PandaX-4T:2021bab , and the DM indirect search
bounds from FERMI-LAT Fermi-LAT:2017bpc  and MAGIC MAGIC:2016xys .

A natural consequence of a first-order phase transition (FOPT) in the early
Universe is the generation of relic gravitational waves (GW) via the
nucleation of bubbles of the broken phase Apreda:2001us ; Grojean:2004xa ;
Weir:2017wfa ; Alves:2018oct ; Alves:2018jsw ; Alves:2019igs ; Alves:2020bpi ;
Chatterjee:2022pxf ; Caprini:2019egz ; Witten:1984rs ; Hogan:1986qda ;
Ellis:2018mja ; Alanne:2019bsm . Such a GW signal can be detected in LISA
LISA:2017pwj , or other proposed ground-based and space-borne experiments,
viz., ALIA Gong:2014mca , Taiji Hu:2017mde , TianQin TianQin:2015yph , aLigo+
Harry:2010zz , Big Bang Observer (BBO) Corbin:2005ny  and Ultimate(U)-DECIGO
Kudoh:2005as , if the amplitude is high enough. These studies receive more
attention after the observation of GW from LIGO and VIRGO collaborations
LIGOScientific:2016aoc ; LIGOScientific:2017vwq ; LIGOScientific:2018mvr ;
LIGOScientific:2020ibl . Note that, the first evidence of stochastic GWs has
been announced very recently by NanoGrav NANOGrav:2023gor  and EPTA
EPTA:2023fyk  collaborations. This discovery provides a new opportunity to
explore physics scenarios beyond the SM. Certain regions of the parameter
space of the present model can be probed via various future GW experiments due
to the production of the GW as a result of an FOEWPT in the scalar sector.

In this paper, we study the interplay among the three sectors mentioned above,
i.e. the production of GW resulting from FOEWPT, the DM prospect, and the LHC
probes of this model. We perform a dedicated scan of the motivated region of
the parameter space. Thereafter, we impose various theoretical constraints,
such as boundedness of the scalar potential from below, electroweak vacuum
stability at zero temperature, perturbativity, unitarity, as well as the
relevant experimental constraints coming from the neutrino sector, flavor
sector, electroweak precision study, dark sector and the heavy Higgs searches
at the LHC. We allow minimal mixing between SM-like Higgs and the triplet-like
neutral scalar to comply with the Higgs signal strength bounds. We further
study the dependencies of the decay width of the SM Higgs into a photon pair
on the model parameters. We then discuss the overall phenomenology of a
complex scalar DM in the presence of the scalar triplet and the SM doublet and
the novel effect of $Z_{3}-$symmetry induced semi-annihilation processes. We
split the study of the FOPT of this model into two parts. First, we examine
the parameter space, where an FOEWPT is facilitated by the triplet in
conjunction with the SM doublet. Next, we study the feasibility of generating
an FOEWPT with the help of the dark sector. Thereafter, we establish a
correlation between the precise measurement of the observed Higgs boson
properties with the detection prospect of the produced GW when the doublet-
triplet interactions are strong. Then we analyze the connection between a
large DM direct detection spin-independent (DMDDSI) cross-section and an
FOEWPT due to significant interaction between DM and the SM Higgs doublet. To
illustrate such correlations, we present a few benchmark points in this work
and discuss the patterns of phase transitions in each case followed by
thorough discussions on the detection prospect at various GW detectors, the
signal-to-noise ratio (SNR) for LISA and possible complementing searches at
the HL-LHC and future DM direct-detection experiments.

The paper is organized as follows: In section 2, we briefly discuss the
theoretical framework of the model, including the interaction pieces of the
Lagrangian, the WIMP DM production mechanism and its possible detection
prospects, the study of EWPT and the production of GW from it. The relevant
theoretical and experimental constraints on the model parameter space are
outlined in section 3. The motivated region of parameter space for this study
has been discussed in section 4. In section 5, we present the results of this
work. Finally, we summarize the outcome of our discussion in section 6.
Various tadpole relations, field-dependent masses of scalar and fermionic
degrees of freedom with thermal correlations (daisy coefficients), various
Feynman diagrams for DM annihilation and various constraints are presented in
the appendix.

##  2 The theoretical framework

###  2.1 The Model

As discussed above, despite many successes of the SM, it cannot explain the
smallness of neutrino masses, the observed matter-antimatter asymmetry and
accommodate a particle DM candidate among other things. To address these three
issues, we consider a model by introducing an additional discrete $Z_{3}$
symmetry and enlarging the particle content of the model by adding a
$SU(2)_{L}$ scalar triplet, $\Delta$, with hypercharge, $Y=2$, and a complex
scalar singlet, $S$, which transforms non-trivially under $Z_{3}$ Yang:2021buu
. The quantum numbers of the BSM scalars in our model along with the SM Higgs
doublet under the extended gauge group $SU(3)_{C}\times SU(2)_{L}\times
U(1)_{Y}\times Z_{3}$ are tabulated in Table 1.

Fields | $\underbrace{SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}}$ $\otimes Z_{3}$
---|---
Complex Scalar DM | $S=\frac{1}{\sqrt{2}}(h_{s}+ia_{s})$ |  1 1 0 $e^{i\frac{2\pi}{3}}$
Scalar Triplet | $\Delta=\left(\begin{matrix}\frac{\Delta^{+}}{\sqrt{2}}&\Delta^{++}\\\ \frac{1}{\sqrt{2}}\big{(}h_{t}+ia_{t}\big{)}&-\frac{\Delta^{+}}{\sqrt{2}}\end{matrix}\right)$ |  1 3 2 1
Higgs doublet | $H=\left(\begin{matrix}G^{+}\\\ \frac{1}{\sqrt{2}}\big{(}h_{d}+ia_{d}\big{)}\end{matrix}\right)$ |  1 2 1 1

Table 1:  Charge assignments of the scalar fields in the model under the gauge
group $\mathcal{G}\equiv\mathcal{G}_{\rm SM}\otimes Z_{3}$ where
$\mathcal{G}_{\rm SM}\equiv SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}$.
Hypercharge ($Y$) of the field is obtained by using the relation:
$Q=I_{3}+\frac{Y}{2}$, where $I_{3}$ is the third component of isospin and $Q$
is the electromagnetic charge.

The part of the Lagrangian relevant for our study in this paper is given
below:

$\displaystyle\mathcal{L}\supset\left(D^{\mu}H\right)^{\dagger}\left(D_{\mu}H\right)+{\rm
Tr}\left[\left(D^{\mu}\Delta\right)^{\dagger}\left(D_{\mu}\Delta\right)\right]-V(\Delta,H)+\mathcal{L}^{\rm
Yukawa}+\mathcal{L}^{\rm DM},$ (1)

where

$\displaystyle D_{\mu}H$ $\displaystyle=$
$\displaystyle\Big{(}\partial_{\mu}-ig_{2}\frac{\sigma^{a}}{2}W_{\mu}^{a}-ig_{1}\frac{Y_{H}}{2}B_{\mu}\Big{)}H,$
$\displaystyle D_{\mu}\Delta$ $\displaystyle=$
$\displaystyle\partial_{\mu}\Delta-
ig_{2}\left[\frac{\sigma^{a}}{2}W_{\mu}^{a},\Delta\right]-{ig_{1}}\frac{Y_{\Delta}}{2}B_{\mu}\Delta~{}.$
(2)

The most general scalar potential involving the SM Higgs doublet and the
scalar triplet can be written as Arhrib:2011uy ; Yang:2021buu :

$\displaystyle V(H,\Delta)$ $\displaystyle=$
$\displaystyle-\mu_{H}^{2}(H^{\dagger}H)+\lambda_{H}(H^{\dagger}H)^{2}~{}$ (3)
$\displaystyle+\mu_{\Delta}^{2}{\rm
Tr}\left[\Delta^{\dagger}\Delta\right]+\lambda_{1}\left(H^{\dagger}H\right){\rm
Tr}\left[\Delta^{\dagger}\Delta\right]+\lambda_{2}\left({\rm
Tr}[\Delta^{\dagger}\Delta]\right)^{2}$ $\displaystyle+\lambda_{3}~{}{\rm
Tr}[\left(\Delta^{\dagger}\Delta\right)^{2}]+\lambda_{4}~{}\left(H^{\dagger}\Delta\Delta^{\dagger}H\right)+\left[\mu\left(H^{T}i\sigma^{2}\Delta^{\dagger}H\right)+h.c.\right]~{}~{}.$

Please note that in the above potential $\mu_{\Delta}^{2}>0$ and the Higgs
triplet cannot trigger any spontaneous symmetry breaking (SSB). So, similar to
the SM, the EWSB happens when the doublet, $H$, develops a vacuum expectation
value (VEV), $\langle{H}\rangle={v_{d}}/{\sqrt{2}}$. However, the presence of
the cubic term $H^{T}i\sigma^{2}\Delta^{\dagger}H$ induces a non-vanishing VEV
for the triplet, $\langle{\Delta}\rangle=v_{t}/\sqrt{2}$, after the EWSB. Then
the scalar fields, $H$ and $\Delta$ can be represented as:

$\displaystyle H=\left(\begin{matrix}G^{+}\\\
\frac{h_{d}+v_{d}+ia_{d}}{\sqrt{2}}\end{matrix}\right)~{},~{}~{}\Delta=\left(\begin{matrix}\frac{\Delta^{+}}{\sqrt{2}}&\Delta^{++}\\\
\frac{h_{t}+v_{t}+ia_{t}}{\sqrt{2}}&-\frac{\Delta^{+}}{\sqrt{2}}\end{matrix}\right).$
(4)

Note that $\sqrt{v_{d}^{2}+2v_{t}^{2}}=v=246$ GeV and in the alignment limit
$v_{t}<<v_{d}$. Minimizing the scalar potential at the vacuums ($v_{d}$ and
$v_{t}$), the bare masses of the scalar potential can be expressed in terms of
other free parameters and can obtain the following conditions:

$\displaystyle\mu_{H}^{2}$ $\displaystyle=$
$\displaystyle-\frac{2M_{\Delta}^{2}v_{t}^{2}}{v_{d}^{2}}+\frac{1}{2}v_{t}^{2}(\lambda_{1}+\lambda_{4})+\lambda_{H}v_{d}^{2}$
$\displaystyle\mu_{\Delta}^{2}$ $\displaystyle=$ $\displaystyle
M_{\Delta}^{2}-v_{t}^{2}(\lambda_{2}+\lambda_{3})-\frac{1}{2}v_{d}^{2}(\lambda_{1}+\lambda_{4})~{}~{}~{}~{}{\rm
with}~{}~{}~{}M_{\Delta}^{2}=\frac{\mu v_{d}^{2}}{\sqrt{2}v_{t}}.$ (5)

The field composition of $V(H,\Delta):$ The two CP even states ($h_{d}$ and
$h_{t}$) are mixed up after the EWSB and give rise to two physical eigenstates
($h^{0}$, $H^{0}$) under the orthogonal transformation. The physical and
interaction states are related as,

$\displaystyle
h^{0}=h_{d}~{}\cos\theta_{t}+h_{t}~{}\sin\theta_{t}~{},~{}~{}~{}~{}H^{0}=-h_{d}~{}\sin\theta_{t}+h_{t}~{}\cos\theta_{t}$
(6)

where $\theta_{t}$ is the mixing angle defined by

$\displaystyle\tan 2\theta_{t}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}\mu
v_{d}-(\lambda_{1}+\lambda_{4})v_{d}~{}v_{t}}{M_{\Delta}^{2}-\frac{1}{4}\lambda_{H}v_{d}^{2}+(\lambda_{2}+\lambda_{3})v_{t}^{2}}~{}.$
(7)

The corresponding mass eigenvalues of these physical eigenstates are given by

$\displaystyle m_{h^{0}}^{2}$ $\displaystyle=$
$\displaystyle\left(M_{\Delta}^{2}+2v_{t}^{2}(\lambda_{2}+\lambda_{3})\right)\sin^{2}\theta_{t}+2\lambda_{H}v_{d}^{2}\cos^{2}\theta_{t}-\frac{v_{t}\sin
2\theta_{t}\left(2M_{\Delta}^{2}-v_{d}^{2}(\lambda_{1}+\lambda_{4})\right)}{v_{d}}$
$\displaystyle m_{H^{0}}^{2}$ $\displaystyle=$
$\displaystyle\left(M_{\Delta}^{2}+2v_{t}^{2}(\lambda_{2}+\lambda_{3})\right)\cos^{2}\theta_{t}+2\lambda_{H}v_{d}^{2}\sin^{2}\theta_{t}+\frac{v_{t}\sin
2\theta_{t}\left(2M_{\Delta}^{2}-v_{d}^{2}(\lambda_{1}+\lambda_{4})\right)}{v_{d}}.$

With the mass hierarchy $m_{h^{0}}<m_{H^{0}}$, the lighter state, $h^{0}$ acts
like the SM higgs with mass, $m_{h^{0}}=125$ GeV.

Similarly, the mixing between two CP odd states ($a_{d}$, $a_{t}$) leads to
one massless Goldstone mode which is associated with the gauge boson $Z$ in
the SM and one massive pseudo scalar, $A^{0}$ and its mass ($m_{A^{0}}$) is
given by,

$\displaystyle m_{A^{0}}^{2}$ $\displaystyle=$
$\displaystyle\frac{M_{\Delta}^{2}\left(4v_{t}^{2}+v_{d}^{2}\right)}{v_{d}^{2}}~{}~{}.$
(9)

In the singly charged scalar sector, the orthogonal rotation of $G^{\pm}$ and
$\Delta^{\pm}$ fields lead to one massless Goldstone state absorbed by the
longitudinal components of the SM gauge boson $W^{\pm}$ and one massive
charged scalar eigenstate, $H^{\pm}$. The mass of $H^{\pm}$ is given by,

$\displaystyle m_{H^{\pm}}^{2}$ $\displaystyle=$
$\displaystyle\frac{\left(2v_{t}^{2}+v_{d}^{2}\right)\left(4M_{\Delta}^{2}-\lambda_{4}v_{d}^{2}\right)}{4v_{d}^{2}}.$
(10)

The doubly charged scalar ($H^{\pm\pm}\equiv\Delta^{\pm\pm}$) mass is given
by,

$\displaystyle m_{H^{\pm\pm}}^{2}$ $\displaystyle=$ $\displaystyle
M_{\Delta}^{2}-\lambda_{3}v_{t}^{2}-\frac{\lambda_{4}v_{d}^{2}}{2}~{}~{}.$
(11)

Therefore the scalar potential, $V(H,\Delta)$ contains seven massive physical
states: two CP even states ($h^{0},~{}H^{0}$), one CP odd state $(A^{0})$, two
singly charged scalar ($H^{\pm}$) and two doubly charged scalar
($H^{\pm\pm}$). The scalar potential, $V(H,\Delta)$ has seven independent free
parameters as

$\displaystyle\\{M_{\Delta},~{}v_{t},~{}\lambda_{H}(\simeq
0.129),~{}\lambda_{1},~{}\lambda_{2},~{}\lambda_{3},~{}\lambda_{4}\\}$ (12)

whereas the physical masses and mixing angle can be represented in terms of
free parameters using above equations.

In the limiting case, $v_{t}^{2}/v_{d}^{2}<<1$, some interesting mass
relations are followed by:

$\displaystyle\big{(}m_{H^{\pm\pm}}^{2}-m_{H^{\pm}}^{2}\big{)}\approx-\frac{\lambda_{4}v_{d}^{2}}{4}~{};~{}~{}\big{(}m_{H^{\pm}}^{2}-m_{A^{0}}^{2}\big{)}\approx-\frac{\lambda_{4}v_{d}^{2}}{4}~{};~{}~{}m_{H^{0}}\approx
m_{A^{0}}\approx M_{\Delta}.$ (13)

Depending on the sign of the quartic coupling, $\lambda_{4}$, there are two
types of mass hierarchy among different components of the triplet scalar:

  * •

when $\lambda_{4}$ is negative: $m_{H^{\pm\pm}}>m_{H^{\pm}}>m_{H^{0},A^{0}}$,

  * •

when $\lambda_{4}$ is positive: $m_{H^{\pm\pm}}<m_{H^{\pm}}<m_{H^{0},A^{0}}$.

In this work, we particularly focus on the first scenario where $\lambda_{4}$
is negative. The mass difference between $H^{\pm\pm}$ and $H^{\pm}$ defined
here by $\Delta m$ as: $\Delta m=m_{H^{\pm\pm}}-m_{H^{\pm}}$.

Neutrino masses: The Yukawa Lagrangian involving the SM lepton doublets ($L$)
and scalar triplet ($\Delta$) to generate neutrino masses can be written as
FileviezPerez:2008jbu ,

$-\mathcal{L}^{\rm
Yukawa}=(y_{L})_{\alpha\beta}~{}\overline{L_{\alpha}^{c}}~{}i\sigma^{2}\Delta~{}L_{\beta}+h.c.$
(14)

where $y_{L}$ is a $3\times 3$ symmetric complex matrix. The heavy triplet
scalar helps to generate neutrino masses via the familiar type-II seesaw
mechanism. Neutrino masses are generated once the neutral triplet scalar gets
induced non-zero vev $v_{t}$. The light neutrino mass in the flavor basis is
given by FileviezPerez:2008jbu

$(m_{\nu})_{\alpha\beta}=\sqrt{2}(y_{L})_{\alpha\beta}{v_{t}},$ (15)

where $\\{\alpha,\beta\\}$ be the flavour indices and $m_{\nu}$ be the
$3\times 3$ light neutrino mass matrix. In order to obtain the physical
neutrino masses, the mass matrix, $m_{\nu}$ can be diagonalised using the
Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, $U_{\text{PMNS}}$, as

$U_{\text{PMNS}}^{T}~{}m_{\nu}~{}U_{\text{PMNS}}=\text{diag}(m_{1},m_{2},m_{3}),$
(16)

where

$U_{\text{PMNS}}=\begin{pmatrix}c_{12}c_{13}&s_{12}c_{13}&s_{13}e^{-i\delta}\\\
-c_{23}s_{12}-s_{23}s_{13}c_{12}e^{i\delta}&c_{23}c_{12}-s_{23}s_{13}s_{12}e^{i\delta}&s_{23}c_{13}\\\
s_{23}s_{12}-c_{23}s_{13}c_{12}e^{i\delta}&-s_{23}c_{12}-c_{23}s_{13}s_{12}e^{i\delta}&c_{23}c_{13}\end{pmatrix}\begin{pmatrix}1&0&0\\\
0&e^{i\alpha_{2}}&0\\\ 0&0&e^{i\alpha_{3}}\end{pmatrix}.$ (17)

$m_{i}~{}(i=1,2,3)$ is the eigenvalue of the light neutrino mass eigenstate
$\nu_{i}$. Here $c_{ij}=\cos\theta_{ij}$ and $s_{ij}=\sin\theta_{ij}$ (with
$i,j=1,2,3$ and $i\neq j$) are three mixing angles (with
$0\leq\theta_{ij}\leq\pi/2$) in neutrino sector and $\delta$ denotes the Dirac
CP-phase (with $0\leq\delta\leq 2\pi$) responsible for CP violation in
neutrino oscillations. Here $\alpha_{2,3}$ be the Majorana phases confined to
the range $[0,\pi]$. The parameters of neutrino oscillations are related to
the mass-squared differences defined as $\Delta m_{sol}^{2}\equiv\Delta
m_{21}^{2}=m_{2}^{2}-m_{1}^{2}$, $\Delta m_{atm}^{2}\equiv\Delta
m_{31}^{2}=m_{3}^{2}-m_{1}^{2}$ for $m_{3}>m_{2}>m_{1}$ (Normal Hierarchy
(NH)) and $\Delta m_{atm}^{2}\equiv\Delta m_{32}^{2}=m_{3}^{2}-m_{2}^{2}$ for
$m_{2}>m_{1}>m_{3}$ (Inverted Hierarchy (IH)) as well as with the three mixing
angles ($\theta_{ij}$).

The DM Lagrangian: The $Z_{3}$ symmetry invariant Lagrangian for the scalar
singlet, $S$ which act as a stable DM candidate, can be written as
Yang:2021buu

$\displaystyle\mathcal{L}^{\rm DM}=\frac{1}{2}|\partial_{\mu}S|^{2}-V^{\rm
DM},$

where $V^{\rm DM}$ is the scalar potential involving the DM, $S$. The scalar
potential reads

$\displaystyle V^{\rm DM}$ $\displaystyle=$
$\displaystyle\mu_{S}^{2}~{}|S|^{2}+\lambda_{S}|S|^{4}+\frac{\mu_{3}}{3!}\Big{(}S^{3}+{S^{*}}^{3}\Big{)}$
(18) $\displaystyle+\lambda_{SH}H^{\dagger}H(S^{*}S)+\lambda_{S\Delta}~{}{\rm
Tr}\left[\Delta^{\dagger}\Delta\right]\big{(}S^{*}S\big{)}~{}~{}~{}~{}.$

The complex scalar field, $S$ can be a stable DM particle under the so-called
freeze-out mechanism. After EWSB the mass of the DM turns out to be:

$\displaystyle
m_{S}^{2}=\mu_{S}^{2}+\frac{1}{2}\lambda_{SH}v_{d}^{2}+\frac{1}{2}\lambda_{S\Delta}v_{t}^{2}\equiv
m_{\rm DM}^{2}~{}.$ (19)

The free parameters involved in the dark sector are following

$\displaystyle\\{~{}m_{S}(\equiv m_{\rm
DM}),~{}\mu_{3},~{}\lambda_{SH},~{}\lambda_{S\Delta},~{}\lambda_{S}\\}.$ (20)

In the following subsection, we examine in detail how the aforementioned free
parameters play a role in the phenomenology of DM.

###  2.2 The Dark Matter

In this subsection, we discuss the phenomenology of the scalar DM of this
model in the presence of the scalar triplet and the SM scalar doublet. The
singlet-like DM interacts with the doublet ($H$) and the triplet ($\Delta$)
sectors of the model through the scalar portal interactions
$\lambda_{SH}(S^{*}S)(H^{\dagger}H)$ and $\lambda_{S\Delta}(S^{*}S){\rm
Tr}[\Delta^{\dagger}\Delta]$, respectively. These interactions are essential
for the DM annihilation in the early Universe. At the same time, the cubic
term ( $\mu_{3}(S^{3}+{S^{*}}^{3})$) in the dark scalar potential has a novel
feature of DM semi-annihilation. At the early time of the Universe, DM is in
thermal equilibrium with the bath particles via the portal interactions as
well as the cubic interactions in the dark sector. During that period of time,
the interaction rate between DM and thermal bath particles (TBPs)
($\Gamma_{\rm DM-TBP}=n_{S}\langle{\sigma v}_{SS\to{\rm TBPs}}\rangle$)
dominates over the Hubble expansion rate ($\mathcal{H}$). With the expansion
of the Universe, DM freezes out from the thermal bath when $\Gamma_{\rm DM-
TBP}<\mathcal{H}$ yielding the DM density of the current Universe. This
phenomenon is often referred to as the WIMP miracle Kolb:1990vq . In this
setup, the number density of DM is mainly governed by several number changing
processes which are classified as:

$\displaystyle{\rm DM~{}annihilation}:\hskip 8.5359pt$
$\displaystyle(i)~{}S~{}S^{*}~{}\to{\rm SM~{}~{}SM}~{}~{}$
$\displaystyle(ii)~{}S~{}S^{*}~{}\to{\rm X~{}~{}Y};~{}~{}~{}\\{\rm
X,Y\\}=\\{H^{0},A^{0},H^{\pm},H^{\pm\pm}\\}$ $\displaystyle{\rm DM~{}semi-
annihilation}:\hskip 8.5359pt$ $\displaystyle~{}~{}~{}~{}~{}S~{}S~{}\to
S~{}\xi;~{}~{}~{}~{}~{}~{}~{}\\{\xi\\}=\\{h^{0},H^{0}\\}~{},$

where the Feynman diagrams of these number changing processes of DM are shown
in figures 16, 17 and 18, respectively, in appendix D. The evolution of DM
number density ($n_{S}$) in terms of more convenient variables, such as the
co-moving number density ($Y_{S}=\frac{{n_{S}}}{s}$) and dimensionless
parameter $x=\frac{m_{S}}{T}$, is described by the Boltzmann equation(BEQ) in
the following form Hektor:2019ote ; Yang:2021buu :

$\displaystyle\frac{dY_{S}}{dx}$ $\displaystyle=$
$\displaystyle-\frac{s~{}\langle\sigma v\rangle_{S~{}S\to
SM~{}SM}}{\mathcal{H}~{}x}(Y_{S}^{2}-{Y_{S}^{eq}}^{2})~{}\Theta\Big{(}m_{S}-m_{SM}\Big{)}$
(21) $\displaystyle-\frac{s}{\mathcal{H}~{}x}\sum_{X,Y}\langle\sigma
v\rangle_{S~{}S\to
X~{}Y}(Y_{S}^{2}-{Y_{X}^{eq}}{Y_{Y}^{eq}})~{}\Theta\Big{(}2m_{S}-m_{X}-m_{Y}\Big{)}$
$\displaystyle-\frac{1}{2}\frac{s~{}\langle\sigma v\rangle_{S~{}S\to
SSM}}{\mathcal{H}~{}x}(Y_{S}^{2}-Y_{S}Y_{S}^{eq})~{}\Theta\Big{(}m_{S}-m_{SM}\Big{)}~{}.$

Depending on the mass hierarchies among the DM and various scalars of this
model, different number changing processes start to contribute to DM number
density. This fact has been illustrated in the BEQ using the theta function
($\Theta$). The co-moving equilibrium density of $i$th particle is defined as:
$Y_{i}^{eq}=\frac{45}{4\pi^{4}}\frac{g_{i}}{g_{*s}}\Big{(}\frac{m_{i}}{T}\Big{)}^{2}K_{2}[\frac{m_{i}}{T}]$
where $g_{i}$ is the degree of freedom of the $i-$th species and $g_{*s}$ is
the internal relativistic degree of freedom associated with entropy. $K_{2}$
is the modified Bessel function of the second kind. The entropy density and
Hubble expansion rates are defined as: $s=g_{*s}\frac{2\pi^{2}}{45}T^{3}$ and
$\mathcal{H}=\sqrt{\frac{\pi^{2}g_{*}}{90}}\frac{1}{M_{\rm Pl}}T^{2}$,
respectively with $g_{*}=106.7$ and $M_{\rm Pl}=2.4\times 10^{18}$ GeV.
$\langle\sigma v\rangle_{SS\to ij}$ is the thermal average cross-section of
the corresponding number changing process: $SS\to ij$ Kolb:1990vq . Solving
the above BEQ, today’s DM relic density is given by Kolb:1990vq ,

$\displaystyle\Omega_{S}h^{2}$ $\displaystyle=$ $\displaystyle 2.752\times
10^{8}~{}\Big{(}\frac{m_{S}}{\rm GeV}\Big{)}~{}Y_{S}(x\to\infty).$ (22)

The DM abundance ($\Omega_{S}h^{2}$) followed from equations 21 and 22 depends
on the following independent free parameters which are entered in thermal
average cross sections:

$\displaystyle\\{m_{S},~{}\mu_{3},~{}\lambda_{SH},~{}\lambda_{S\Delta}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm
from~{}dark~{}sector},$
$\displaystyle\\{~{}\lambda_{1},~{}\lambda_{2},~{}\lambda_{3},~{}\lambda_{4}~{}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm
from~{}scalar~{}triplet~{}sector.}$ (23)

The dark sector self coupling, $\lambda_{S}$ does not affect the computation
of the relic density. We will discuss the role of the above parameters on DM
density in section 5.1. Note that, the present day DM relic density is
constrained from the Planck observation Planck:2018vyg  which will be
discussed in the subsection 3.1.

$S$$n$$h^{0},H^{0}$$S$$n$ Figure 1: Feynman diagram for spin-independent DM-
nucleon scattering process for the scalar DM.

There are several attempts have been made to detect DM in laboratory
experiments. DMDD is one of them where incoming DM flux scatters with the
nuclei in the target crystals and the recoil of the target nuclei can be
searched for a DM signal. The effective spin-independent DM-Nucleon scattering
cross-section ($\sigma_{\rm DD}^{\rm SI}$) mediated via $t-$channel $CP-$even
scalars ($h^{0}$ and $H^{0}$) exchange diagrams, as shown in figure 1, read as
($q\to 0$):

$\sigma_{\rm DD}^{\rm SI}=\Big{(}\frac{\Omega_{S}h^{2}}{\Omega_{\rm
DM}h^{2}}\Big{)}~{}\Big{(}\frac{1}{4\pi}\bigg{(}\frac{f_{n}\mu_{n}}{m_{S}}\bigg{)}^{2}\bigg{(}\frac{m_{n}}{v_{d}}\bigg{)}^{2}~{}\bigg{[}\frac{\lambda_{D1}}{m_{h^{0}}^{2}}+\frac{\lambda_{D2}}{m_{H^{0}}^{2}}\bigg{]}^{2}\Big{)},$
(24)

where $\Big{(}{\Omega_{S}h^{2}}/{\Omega_{\rm DM}h^{2}}\Big{)}\equiv\xi_{\rm
DM}$ is the fractional DM density. The coupling $\lambda_{D1}$ and
$\lambda_{D2}$ are defined as

$\displaystyle\lambda_{D1}$ $\displaystyle=$
$\displaystyle\Big{(}~{}~{}\lambda_{SH}~{}v_{d}\cos\theta_{t}+\lambda_{S\Delta}v_{t}\sin\theta_{t}\Big{)}\cos\theta_{t}~{}\xrightarrow{\theta_{t}\to
0}~{}\simeq\lambda_{SH}~{}v_{d}$ $\displaystyle\lambda_{D2}$ $\displaystyle=$
$\displaystyle\Big{(}-\lambda_{SH}~{}v_{d}\sin\theta_{t}+\lambda_{S\Delta}v_{t}\cos\theta_{t}\Big{)}\sin\theta_{t}\xrightarrow{\theta_{t}\to
0}~{}\simeq 0$ (25)

The reduced mass of DM-nucleon system is defined by $\mu_{n}$
$\big{(}=\frac{m_{S}m_{n}}{m_{S}+m_{n}}\big{)}$ with $m_{n}=0.946$ GeV and the
nucleon form factor $f_{n}=0.28$ Alarcon:2012nr . For the limiting case:
$\theta_{t}\to 0$, $\sigma_{\rm DD}^{\rm SI}$ solely depends on the the Higgs
portal DM coupling, $\lambda_{SH}$ and the DM mass, $m_{S}$. The null results
from the various DM direct search experiments such as XENON-1T Aprile:2018dbl
and PANDA-4T PandaX-4T:2021bab  put a strong constraint on the DMDD cross-
sections which can be further expressed in terms of the model parameters.

At the same time, DM can also be detected at various indirect search
experiments such as the space-based telescopes, Fermi Large Area Telescope
(Fermi-LAT) Fermi-LAT:2017bpc  and the ground-based telescopes, Major
Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) MAGIC:2016xys . These
experiments are looking for the gamma-ray flux which is produced via the
production of the SM particles either through DM annihilation or via decay in
the local Universe. The photons which are emitted from WIMP-like DM lie in the
gamma-ray regime and behave as ideal messengers of DM indirect detection
(DMID). The total gamma-ray flux due to DM annihilation into the SM charge
pairs($\psi\overline{\psi}$) in a specific energy range is given by
MAGIC:2016xys

$\displaystyle\Phi_{\psi\overline{\psi}}$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}\frac{\langle\sigma
v\rangle_{SS\to\psi\overline{\psi}}}{2m_{S}^{2}}\int_{E_{\rm min}}^{E_{\rm
max}}\frac{dN_{\gamma}}{dE_{\gamma}}dE_{\gamma}\underbrace{\int
dx~{}\rho_{S}^{2}\big{(}r(b,l,x)\big{)}}_{J}~{},$ (26)

where $\psi:=\\{\mu^{-},\tau^{-},b,W^{-}\\}$. The $J$ factor contains the
astrophysical information about the DM distribution in the galaxy and
$\frac{dN_{\gamma}}{dE_{\gamma}}$ is the energy spectrum of incoming photons
from DM annihilation. Measuring the gamma-ray flux and using the standard
astrophysical inputs, the indirect search experiments put limits on the
thermal average cross-section of DM annihilation into different SM charge
pairs like
$\mu^{-}\mu^{+},~{}\tau^{-}\tau^{+},~{}b\overline{b},~{}W^{-}W^{+}$. In order
to compare the experimental bound with the theoretically estimated thermal
average cross-section of the corresponding DM annihilation channel, one can
express the above equation 26 as:

$\displaystyle\langle\sigma v\rangle_{SS\to\psi\overline{\psi}}^{\rm eff}$
$\displaystyle=$ $\displaystyle\big{(}\Omega_{S}/\Omega_{\rm
DM}\big{)}^{2}\langle\sigma v\rangle_{SS\to\psi\overline{\psi}}=\xi_{\rm
DM}^{2}\langle\sigma v\rangle_{SS\to\psi\overline{\psi}}~{}~{}.$ (27)

Here the thermal average cross-sections for DMID are scaled with the
fractional DM density as $\xi_{\rm DM}^{2}$. In this setup the thermal average
annihilation cross-sections,$\langle\sigma
v\rangle_{SS\to\psi\overline{\psi}}$ are mostly dominated by $s$-channel SM
Higgs ($h^{0}$) exchange diagram and varies as
$\propto\lambda_{SH}^{2}/m_{S}^{2}$ (for $\sin\theta_{t}\to 0$). Non-
observation of any DM signal in such indirect detection experiments puts a
constraint on the thermal average cross-section of DM annihilation into
different SM charge pairs. Those constraints are further expressed in terms of
the model parameters.

So far, we have talked about the analytical form of DM abundance, DMDDSI
cross-section and DMID cross-section which will help to understand the
numerical results. For the numerical analysis, we have used publicly available
numerical packages. We first implemented the model in ${\tt FeynRules}$
Alloul:2013bka  and then the outputs are fed into ${\tt micrOMEGAs}$
Belanger:2006is  to obtain DM abundance, $\Omega_{S}h^{2}$ and DMDDSI and DMID
cross-sections.

###  2.3 Study of EWPT and the production of GW

In this section, we discuss the theoretical formulation of EWPT in our model
and examine the possibility of the EWPT being a strong first-order one, which
in turn may facilitate EWBG. We also outline the production of stochastic GW
from FOPT.

Note that in the present study of EWPT, we assume that no spontaneous or
explicit ~~$CP$~~ occurs in the Higgs sector. However, one can incorporate a
$Z_{3}$ symmetric $CP$-violating dimension-7 operator,
$y_{t}\,\overline{Q}\tilde{H}\,\big{(}1+ic\frac{S^{3}}{\Lambda^{3}}\big{)}\,t_{R}\,+\,h.c.$,
in this model to generate $CP$\- violation which is essential for EWBG to take
place. We do not perform a detailed calculation of EWBG in this paper since
the focus of the work is to study the correlation/interplay of the prospects
of detecting stochastic GW from an SFOEWPT in an almost unconstrained
parameter space at various future proposed GW experiments, the collider
signals of this parameter space, and the prospects of detecting DM at various
DM direct detection (DMDD) experiments. However, we should mention that the
detailed baryogenesis calculation for an equivalent $Z_{2}$-symmetric
$CP$-violating dimension-6 operator has been performed in the literature
Vaskonen:2016yiu ; Ellis:2022lft . One can follow the same prescription for
our $Z_{3}$ symmetric $CP$-violating dimension-7 operator given above. We
leave such an EWBG calculation for a future study.

####  2.3.1 Effective Higgs potential at finite-temperature

The tree-level scalar potential of the model (see equation 3) relevant for the
study of EWPT, in terms of the $CP$-even Higgs fields {$h_{d}$, $h_{t}$ and
$h_{s}$}, can be written as

$\displaystyle V_{0}(h_{d},h_{t},h_{s})=$
$\displaystyle\,\frac{1}{4}(\lambda_{H}h_{d}^{4}+(\lambda_{2}+\lambda_{3})h_{t}^{4}+\lambda_{S}h_{s}^{4}+(\lambda_{1}+\lambda_{4})h_{d}^{2}h_{t}^{2}+\lambda_{SH}h_{d}^{2}h_{s}^{2}+\lambda_{S\Delta}h_{t}^{2}h_{s}^{2})$
$\displaystyle+\frac{\mu_{3}}{6\sqrt{2}}h_{s}^{3}-\frac{\mu}{\sqrt{2}}h_{d}^{2}h_{t}-\frac{\mu_{H}^{2}}{2}h_{d}^{2}+\frac{\mu_{\Delta}^{2}}{2}h_{t}^{2}+\frac{\mu_{S}^{2}}{2}h_{s}^{2}.$
(28)

The vacuum expectation values (${\it vevs}$) of the scalar fields are assumed
to be real at all temperatures. As we have discussed earlier that at $T=0$,
the minimum of the scalar potential is $\langle h_{d}\rangle\equiv v_{d}$,
$\langle h_{t}\rangle\equiv v_{t}$ and $\langle h_{s}\rangle\equiv 0$ where
$v=\sqrt{v_{d}^{2}+v_{t}^{2}}=246$ GeV. The zero temperature tree-level
potential (as described in equation 2.3.1) gets quantum corrections from all
fields that couple to $h_{d}$, $h_{t}$ and $h_{s}$. These corrections can be
written in terms of the well-known Coleman-Weinberg (CW) Coleman:1973jx  one-
loop potential in the following form:

$\displaystyle V_{\rm
CW}(h_{d},h_{t},h_{s},T=0)=\sum_{i}(-1)^{2s_{i}}n_{i}\frac{m_{i}^{4}(h_{d},h_{t},h_{s})}{64\pi^{2}}\left[\ln\frac{m_{i}^{2}(h_{d},h_{t},h_{s},)}{Q^{2}}-C_{i}\right]\;\
,$ (29)

where ‘$Q$’ is the renormalization scale which we set at $Q=246$ GeV . The sum
‘$i$’ runs over all particles in the model. The constants $C_{i}$ depend on
the choice of the renormalization scheme. In this work, we consider
$\overline{\text{MS}}$ on-shell scheme, whereby for the transverse
polarizations of the gauge bosons $C_{i}=\frac{1}{2}$ and for their
longitudinal polarizations and for all other particles (scalars and fermions)
$C_{i}=\frac{3}{2}$, i.e., $C_{W^{\pm},Z}=5/6$ and for other species
$C_{i}=3/2$. Here, $m_{i}$ is the field-dependent mass of the $i$-th species
(see Appendix C for the details), and $n_{i}$ and $s_{i}$ are their associated
number of degrees of freedom and spin, respectively. The number of degrees of
freedom of all the fields in the model relevant for equation 29 is given
below:

$n_{h_{d},h_{t},h_{s},a_{d},a_{t},a_{s},G^{0}}=1,\,\,n_{H_{\Delta}^{\pm},H_{\Delta}^{\pm\pm},G^{\pm}}=2,\,\,n_{Z}=3,\,\,n_{W^{\pm}}=6,\,\,n_{t}=12\,\,.$
(30)

In this work, we work in the Landau gauge as the ghost bosons decouple in this
gauge and we do not need to consider them in our calculation.

Note that the location of the electroweak minimum in the field space for the
tree-level potential is shifted by the one-loop CW potential. Thus, suitable
counter-terms need to add to the effective potential to ensure that the
minimum of the effective potential coincides with that of the tree-level
potential in the $\overline{\text{MS}}$ renormalization scheme. The added
counter-terms to the potential are parameterized as,

$\displaystyle V_{\rm
CT}=\delta_{\mu_{H}^{2}}h_{d}^{2}+\delta_{\mu_{\Delta}^{2}}h_{t}^{2}+\delta_{\lambda_{H}}h_{d}^{4}+\delta_{\lambda_{2}}h_{t}^{4}+\delta_{\lambda_{S\Delta}}h_{d}^{2}h_{t}^{2}+\delta_{\lambda_{SH}}h_{d}^{2}h_{s}^{2}\;.$
(31)

Various coefficients in $V_{\rm CT}$ (of equation 31) can be found from the
on-shell renormalization conditions at zero temperature:

$\frac{\partial(V_{\text{CW}}+V_{\text{CT}})}{\partial\phi_{{}_{i}}}=0\,\,\,,\,\,\,\\\
\frac{\partial^{2}(V_{\text{CW}}+V_{\text{CT}})}{\partial\phi_{{}_{i}}\partial\phi_{{}_{i}}}=0\,\,\,,$
(32)

where $\phi_{{}_{i}},\phi_{{}_{j}}=\\{h_{d},h_{t},h_{s}\\}$. All the
derivatives are taken at the true EW minima, i.e., $h_{d}=v_{d}$,
$h_{t}=v_{t}$ and $h_{s}=0.$ The various coefficients of the counter-term
potentials are presented in the appendix B. Note that, in the Landau gauge
that we opt for this work, the masses of the Goldstone modes vanish at the
physical minimum. As a result, when physical masses and coupling coefficients
are estimated from derivatives of the loop-corrected effective potential, it
leads to divergences Martin:2014bca ; Elias-Miro:2014pca  and it can be
handled by using an infrared regulator by changing the masses of the Goldstone
modes, $m_{G}^{2}\to m_{G}^{2}+\mu_{\rm IR}^{2}$. In this work, we adopt the
approach which is considered in reference Baum:2020vfl ; Chatterjee:2022pxf
i.e., for numerical calculation it is sufficient to set $\mu_{\rm
IR}^{2}=1\,{\rm GeV}^{2}$.

At finite-temperatures, the effective potential receives additional
corrections and that at the one-loop level are given by Dolan:1973qd ;
Weinberg:1974hy ; Kirzhnits:1976ts

$\displaystyle V_{\rm
th}(h_{d},h_{t},h_{s},T)=\frac{T^{4}}{2\pi^{2}}\,\sum_{i}n_{i}J_{B,F}\left(\frac{m_{i}^{2}(h_{d},h_{t},h_{s})}{T^{2}}\right)\;,$
(33)

where $n_{i}$ are the numbers of degrees of freedom for particles as discussed
earlier. The thermal functions $J_{B(F)}$ are for bosonic (fermionic)
particles defined as

$J_{B/F}(y^{2})=\pm{\rm
Re}\int_{0}^{\infty}x^{2}\ln\left(1\mp\exp^{-\sqrt{x^{2}+y^{2}}}\right){\rm
d}{x}\,,$ (34)

where the lower (upper) sign corresponds to fermions (bosons), $\beta\equiv
1/\text{T}$ and the sum includes all the particles as described in equation
29.

Note that, at high temperature the treatment of the perturbative expansion of
the effective potential no longer remain valid. The quadratically divergent
contributions from the so-called non-zero Matsubara modes need to be re-summed
through the insertion of thermal masses in the one-loop propagator. The
corrections are only to the Matsubara zero-modes, i.e., to the longitudinal
polarization states of the vector bosons and to all the scalars. Since the
thermal contributions to the transverse modes are suppressed due to gauge
symmetry Espinosa:1992kf . To make the expansion reliable, we adopt the well-
known Parwani method Parwani:1991gq  and we denote the thermally improved
field-dependent masses as $M_{i}^{2}$, where
$M_{i}^{2}=\text{eigenvalues}[m_{i}^{2}+\Pi(T^{2})]$ where
$\Pi(T^{2})=c_{ij}T^{2}$ and $c_{ij}$’s are the so-called daisy coefficients.
Note that the gauge symmetries plus the discrete $Z_{3}$-symmetry of the
present model set the off-diagonal terms of the $\Pi(T^{2})$ matrix to zero.
Thus, the Daisy coefficients are listed in equations 68.
Including the Coleman-Weinberg and the thermal corrections, the temperature-
dependent effective potential at one-loop order is given by

$V_{\text{eff}}(T)=V_{0}+V_{\rm CW}(M_{i}^{2})+V_{\rm CT}+V_{\rm
th}(M_{i}^{2},T)\,.$ (35)

The EWPT may now be investigated using this potential, with the minimum of the
potential being tracked as a function of temperature. Note that the total one-
loop effective potential (equation 35) carries explicit gauge dependence.
Since, the important quantities of the study of phase transitions such as the
locations of the extrema of $V_{\text{eff}}(T)$, as well as the ratio
$\phi_{c}(T_{c})/T_{c}$, both are gauge-dependent Dolan:1973qd ;
Nielsen:1975fs ; Fukuda:1975di ; Laine:1994zq ; Baacke:1993aj ; Baacke:1994ix
; Garny:2012cg ; Espinosa:2016nld ; Patel:2011th ; Arunasalam:2021zrs ;
Lofgren:2021ogg . 111 Using the Nielsen identities Nielsen:1975fs ;
Fukuda:1975di , we can get gauge-independent variables of the effective
potential. In addition, the one-loop effective potential of equation 35
explicitly depends on the choice of the renormalization scale ($Q$), which
might have a more significant impact than the gauge uncertainty Laine:2017hdk
.

The production of the stochastic GW from FOPTs will be discussed in the next
subsection. Later, we shall use this to identify a suitable region of
parameter space of the model that can be investigated using various future GW
experiments.

####  2.3.2 Production of stochastic GW from FOPT

In this $Z_{3}$-symmetric singlet scalar extended type-II seesaw model, multi-
step FOPTs may occur in the early Universe, potentially generating a
stochastic background of GW. The amplitude of this stochastic background, in
contrast to GW from a binary system, is a random quantity that is unpolarized,
isotropic and has a Gaussian distribution Allen:1997ad . As a result, the two-
point correlation function can characterize this, and it is proportional to
the power spectral density $\Omega_{\text{GW}}{\rm h}^{2}$. The “cross-
correlation” method between two or more detectors can be used to detect this
kind of stochastic GW Caprini:2015zlo ; Cai:2017cbj ; Caprini:2018mtu ;
Romano:2016dpx ; Christensen:2018iqi .

To determine the GW signals arising from FOPTs and evaluate their evolution
with temperature we require the knowledge of the following set of portal
parameters: $T_{n},\alpha,\beta/H_{n},v_{w}$. Here, the bubble nucleation
temperature, $T_{n}$, is the approximate time at which one bubble per Hubble
volume forms due to a phase transition. A dimensionless quantity, ‘$\alpha$’,
is defined as the ratio of the energy released from the phase transition to
the total radiation energy density at the time when the phase transition is
complete. Thus, it relates to the energy budget of the FOPT and is given by
Espinosa:2010hh

$\alpha=\frac{\rho_{\text{vac}}}{\rho^{*}_{\text{rad}}}=\frac{1}{\rho^{*}_{\text{rad}}}\left[T\frac{\partial\Delta
V(T)}{\partial T}-\Delta V(T)\right]\Bigg{|}_{T_{*}},$ (36)

where $T_{*}=T|_{t_{*}}$. $t_{*}$ is the time when the FOPT completes. In the
absence of substantial reheating effects, $T_{*}\simeq T_{n}$. The total
radiation energy density of the plasma background is
$\rho^{*}_{\text{rad}}=g_{*}\pi^{2}T_{*}^{4}/30$, where $g_{*}$ is the number
of relativistic degrees of freedom at $T=T_{*}$. In this work, we consider
$g_{*}\sim 100$. The difference between the potential energies at the false
minimum and the true minimum is denoted by $\Delta
V(T)=V_{\text{false}}(T)-V_{\text{true}}(T)$. The other two parameters are the
parameter $\beta$, which roughly indicates the inverse time duration of the
phase transition, and the parameter $v_{w}$, which is the bubble wall
velocity.

The tunnelling probability from the false vacuum to the true vacuum per unit
time per unit volume is given by Turner:1992tz

$\Gamma(T)\simeq T^{4}\left(\frac{S_{3}\left(T\right)}{2\pi
T}\right)^{3/2}e^{-S_{3}\left(T\right)/T},$ (37)

where the Euclidean action $S_{3}\left(T\right)$ of the background field
($\phi$), in the spherical polar coordinate, is given by

$S_{3}\left(T\right)=4\pi\int dr\hskip
2.84526ptr^{2}\left[\dfrac{1}{2}\left(\partial_{r}\vec{\phi}\right)^{2}+V_{eff}\right]\,\,,$
(38)

with $V_{eff}$ being the finite temperature total effective potential defined
in equation 35, and $\vec{\phi}$ denoting the three components ${\it vev}$ of
the scalar fields. The condition for the formation of a bubble of critical
size can be found by extremizing this Euclidean action. For this study, we use
the publicly available toolbox $\mathbf{CosmoTransitions}$ Wainwright:2011kj
to find bounce solutions of the above Euclidean action. It uses a path
deformation method and this computation is the most technically challenging
part. The bubble nucleation rate ($\Gamma$) pet unit volume at temperature T
is given by $\frac{\Gamma(T)}{V}\propto T^{4}e^{-S_{3}/T}$. The nucleation
temperature is usually determined by solving the following equation of the
nucleation rate,

$\int_{T_{n}}^{\infty}\frac{dT}{T}\frac{\Gamma(T)}{H(T)^{4}}\simeq 1.$ (39)

This equation states that the nucleation probability of a single bubble within
one horizon volume is approximately 1, and it reduces to the condition
$S_{3}(T)/T\approx 140$, solving which one can obtain the nucleation
temperature $T_{n}$ Apreda:2001us . This is the highest temperature for which
$S_{3}/T\lesssim 140$. For all $T>0$, if $S_{3}/T>140$, the corresponding
transition to the true vacuum or nucleation of bubbles of critical size of the
broken phase does not happen. This occurs because either the barrier between
the two local minima is too high or the distance between them is too large in
the field space, resulting in a low tunnelling probability. For $T<T_{c}$, if
this condition is not fulfilled, the system would remain trapped at the
metastable or false minimum. Thus, in spite of the presence of a deeper
minimum of the potential at the zero temperature the Universe remains at the
false minimum. This phenomenon suggests that, rather than just establishing
the existence of a critical temperature for an FOEWPT, it is more important to
investigate the successful nucleation of a bubble of the broken electroweak
phase.

The inverse time duration of the phase transition, $\beta$, is obtained from
the relation

$\beta=-\frac{dS_{3}}{dt}\Bigr{|}_{t_{*}}\simeq\dfrac{\dot{\Gamma}}{\Gamma}=H_{*}T_{*}\frac{d(S_{3}/T)}{dT}\Bigr{|}_{T_{*}},$
(40)

where $H_{*}$ is the Hubble rate at $T_{*}$. For a strong GW signal, the ratio
needs to be low which implies a relatively slow phase transition.

Having defined the portal we need to know two more parameters to calculate the
GW power spectrum coming from an FOPT in the early Universe. When a phase
transition occurs in a thermal plasma, the released energy is shared between
the plasma’s kinetic energy, which causes a bulk motion of the fluid in the
plasma resulting in GW, and heating the plasma. The quantity $\kappa_{v}$ is
the fraction of latent heat energy converted into the bulk motion of the
fluid, which takes the form Espinosa:2010hh ; Chiang:2019oms

$\kappa_{v}\simeq\left[\dfrac{\alpha}{0.73+0.083\sqrt{\alpha}+\alpha}\right]\,\,.$
(41)

Finally, we need to know $\kappa_{\text{turb}}$, which is a fraction of
$\kappa_{v}$ that is used to generate Magneto-Hydrodynamic-Turbulence (MHD) in
the plasma. The value of $\kappa_{\text{turb}}$ is unknown but it is expected
that $\kappa_{\text{turb}}\approx(5\sim 10)\,\kappa_{v}$ Hindmarsh:2015qta ,
and we consider this fractional value to be 0.1 in all our GW calculations.
Now we have all the relevant parameters at our disposal, using which we can
calculate the GW energy density spectrum.

It is generally accepted that during a cosmological FOPT, the GW can be
produced by three distinct processes: bubble wall collisions, long-standing
sound waves in the plasma, and MHD turbulence. For bubble collisions, the GW
is generated by the stress-energy tensor located at the wall, and the
mechanism is referred to as the “envelope approximation”222In fact, the
“Envelope approximation” is actually two approximations. In the first
approximation, the expanding bubble’s stress-energy tensor is believed to be
non-zero only in an infinitesimally thin shell on the bubble surfaces. In the
second approximation, when two bubbles interact, it is assumed that this
stress-energy tensor vanishes instantly, leaving just the ‘envelopes’ of the
bubbles to interact Kosowsky:1992rz ; Kosowsky:1991ua ; Kosowsky:1992vn ;
Caprini:2019egz .. However, for a phase transition proceeding in a thermal
plasma, bubble collisions’ contribution to the total GW energy density is
believed to be negligible Bodeker:2017cim . On the other hand, the bulk motion
of the plasma gives rise to velocity perturbations in it, resulting in the
generation of sound waves in a medium made up of relativistic particles. The
sound waves receive the majority of the energy released during the phase
transition Hindmarsh:2013xza ; Giblin:2013kea ; Giblin:2014qia ;
Hindmarsh:2015qta . This relatively long-living acoustic production of GW is
often regarded as the most dominant one 333Here, we are presuming that a
bubble expanding in plasma can achieve a relativistic terminal velocity. The
energy in the scalar field is negligible in this situation and the most
significant contributions to the signal are predicted to come from the fluid’s
bulk motion in the form of sound waves and/or MHD. This is the non-runaway
bubble in a plasma (NP) scenario Caprini:2015zlo ; Schmitz:2020syl ..
Percolation may also cause turbulence in the plasma, particularly MHD
turbulence since the plasma is completely ionized. This corresponds to the
third source Caprini:2006jb ; Kahniashvili:2008pf ; Kahniashvili:2008pe ;
Kahniashvili:2009mf ; Caprini:2009yp ; Kisslinger:2015hua  of GW production.
Thus, summing the contributions of the above two processes (i.e. only the
contributions from the sound wave and the MHD turbulence) one obtains the
overall GW intensity $\Omega_{\text{GW}}{\rm h}^{2}$ from the SFOPT for a
given frequency $f$,

$\Omega_{\text{GW}}h^{2}\simeq\Omega_{\text{sw}}h^{2}+\Omega_{\text{turb}}h^{2},$
(42)

where $h\approx 0.673$ DES:2017txv .

The contribution to the total GW power spectrum from the sound waves,
$\Omega_{\text{sw}}h^{2}$, in equation 42 can be modelled by the following fit
formula Hindmarsh:2019phv

$\Omega_{\text{SW}}{\rm h}^{2}=2.65\times
10^{-6}\Upsilon(\tau_{SW})\left(\dfrac{\beta}{H_{\star}}\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}}\,\,,$
(43)

where $T_{\star}$ is the temperature just after the end of GW production and
$H_{\star}$ is the Hubble rate at that time. For this analysis, we consider
the approximation that $T_{\star}\approx T_{n}$. The present day peak
frequency $f_{\text{SW}}$ for the sound wave contribution is

$f_{\text{SW}}=1.9\times 10^{-5}\hskip
2.84526pt\text{Hz}\left(\dfrac{1}{v_{w}}\right)\left(\dfrac{\beta}{H_{\star}}\right)\left(\dfrac{T_{n}}{100\hskip
2.84526pt\text{GeV}}\right)\left(\dfrac{g_{*}}{100}\right)^{\frac{1}{6}}\,\,.$
(44)

It has been shown in recent studies Guo:2020grp ; Hindmarsh:2020hop  that
sound waves have a finite lifetime, which causes $\Omega_{\text{sw}}h^{2}$ to
be suppressed. The multiplication factor $\Upsilon(\tau_{SW})$ is considered
in equation 43 to take this suppression into account.

$\Upsilon(\tau_{SW})=1-\frac{1}{\sqrt{1+2\tau_{\text{sw}}H_{\ast}}}.$ (45)

Here, the lifetime $\tau_{\text{sw}}$ is considered as the time scale when the
turbulence develops, approximately given by Pen:2015qta ; Hindmarsh:2017gnf :

$\displaystyle\tau_{\text{sw}}\sim\frac{R_{\ast}}{\bar{U}_{f}},$ (46)

where $R_{\ast}$ is the mean bubble separation, and $\bar{U}_{f}$ is the root-
mean-squared (RMS) fluid velocity. $R_{\ast}$ is related to the duration of
the phase transition parameter, $\beta$, through the relation
$R_{\ast}=(8\pi)^{1/3}v_{w}/\beta$ Hindmarsh:2019phv ; Guo:2020grp . On the
other hand, from hydrodynamic analyses, references Hindmarsh:2019phv ;
Weir:2017wfa  has shown that the RMS fluid velocity is given by,
$\bar{U}_{f}=\sqrt{3\kappa_{v}\alpha/4}$. At
$\tau_{\text{sw}}\rightarrow\infty$, $\Upsilon(\tau_{SW})$ approaches the
asymptotic value $1$.

In contrast, the contribution $\Omega_{\text{turb}}h^{2}$ from the MHD
turbulence in equation 42 can be modeled by Caprini:2015zlo

$\Omega_{\text{turb}}h^{2}=3.35\times
10^{-4}\left(\frac{H_{*}}{\beta}\right)\left(\frac{\kappa_{\text{turb}}\alpha}{1+\alpha}\right)^{\frac{3}{2}}\left(\frac{100}{g_{s}}\right)^{\frac{1}{3}}v_{w}\frac{(f/f_{\text{turb}})^{3}}{[1+(f/f_{\text{turb}})]^{\frac{11}{3}}(1+8\pi
f/h_{\star})},$ (47)

where turbulence produced GW spectrum present day peak frequency
$f_{\text{turb}}$ is given by,

$f_{\text{turb}}=2.7\times
10^{-2}~{}\text{mHz}\frac{1}{v_{w}}\left(\frac{\beta}{H_{*}}\right)\left(\frac{T_{*}}{100~{}\text{GeV}}\right)\left(\frac{g_{s}}{100}\right)^{\frac{1}{6}},$
(48)

with the parameter

$h_{*}=16.5\times 10^{-6}\hskip
2.84526pt\text{Hz}\left(\dfrac{T_{n}}{100\hskip
2.84526pt\text{GeV}}\right)\left(\dfrac{g_{*}}{100}\right)^{\frac{1}{6}}\,\,.$
(49)

Before ending our discussion on the framework of GW production from FOPT in
the early Universe, a few words are necessary on EWBG. For the early Universe
to undergo a successful EWBG via an SFOEWPT, one of the parameters discussed
in this subsection, the bubble wall velocity ($v_{w}$), needs to play an
important role. From the preceding GW calculation, it is clear that a larger
$v_{w}$ is required for a stronger GW production, but to achieve the observed
matter anti-matter asymmetry a small subsonic $v_{w}$ is needed. Consequently,
a significantly large $v_{w}$, which can generate detectable GW signals, is
detrimental to the process of producing the observed baryon asymmetry.
However, recent studies recognise that $v_{w}$ may not be the velocity that
enters EWBG calculation. This is due to the non-trivial plasma velocity
profile surrounding the bubble wall No:2011fi . This will require an in-depth
analysis of the effect of particle transport near the wall and we leave it for
a future study. For the current work, we assume that expanding bubbles reach a
relativistic terminal velocity in the plasma, i.e., $v_{w}\simeq 1$.

##  3 Theoretical and experimental constraints

Having set up the framework for our study, we discuss various theoretical and
phenomenological constraints on the model parameters used in our analysis. The
theoretical constraints coming from the EW vacuum stability at zero
temperature, perturbativity, partial wave unitarity and the phenomenological
constraints from the flavor sector, the neutrino oscillation data and the
electroweak precision observables including the ‘$\rho$’ parameter are well
discussed in the literature Yang:2021buu . For the sake of completeness of
this work, we discuss those in the appendix E and F. Other phenomenological
constraints from the DM experiments and those coming from the neutral and
charged Higgs boson searches at the LHC are discussed below. Such a discussion
is useful to determine the targeted region of parameter space for our work.

###  3.1 DM constraints

Besides neutrino masses, another objective of the model used in this paper is
to explain the observed DM relic density of the Universe. In Section 2.2 we
provided the analytical formulae necessary to calculate the DM abundance,
$\Omega_{S}h^{2}$, and th$\sigma_{\rm DD}^{\rm SI}$) within this model. Here,
we discuss how we study the DM phenomenology of the model numerically. First,
we generated the model files using FeynRules Alloul:2013bka  and then feed the
necessary files into micrOMEGAs(v4.3) Belanger:2006is  and perform a scan over
the model parameters relevant for the dark sector, as listed in equation 23.
The package micrOMEGAs evaluate $\Omega_{S}h^{2}$ and $\sigma_{\rm DD}^{\rm
SI}$ for each benchmark point.

It is well known that the Planck collaboration measured the DM relic abundance
to be $0.120\pm 0.001$, which we take into cognizance in our DM analysis. On
the other hand, the most stringent bounds on $\sigma_{\rm DD}^{\rm SI}$ are
provided by the latest XENON-1T Aprile:2018dbl  and PANDA-4T PandaX-4T:2021bab
results. While imposing these upper limits on $\sigma_{\rm DD}^{\rm SI}$ we
take into account a scale factor ($\xi_{{}_{\text{D}M}}$), which is the ratio
of the computed relic density for the DM particle of the model and the
observed DM relic density of the Universe, i.e.,
$\xi_{{}_{\text{D}M}}=\frac{\Omega_{S}h^{2}}{0.120}$. $\xi_{{}_{\text{D}M}}$
alleviates the strong bounds on $\sigma_{\rm DD}^{\rm SI}$ for a significant
fraction of our scanned points. More details about how the above DM relic
density and DMDD constraints are incorporated in our numerical analysis are
discussed in Section 5.1.

The indirect search of DM followed from the experiments such as Fermi-LAT and
MAGIC also put a constraint on the individual thermal average annihilation
cross-section of DM into SM charge pairs as $\xi_{\rm DM}^{2}\langle\sigma
v\rangle_{SS\to\psi\overline{\psi}}$ where
$\psi:=\\{\mu^{-},\tau^{-},b,W^{-}\\}$ MAGIC:2016xys . When the density of DM
is under-abundant, the scale factor $\xi_{\rm DM}$ further reduces the
effective indirect search cross-section and alleviates bounds on the
corresponding annihilation process. Like the direct search, same portal
coupling, $\lambda_{SH}$ is also involved in the indirect search cross-
section. The direct search of DM puts a strong constraint on the portal
coupling, $\lambda_{SH}$, which leads to a small indirect search cross-section
well below the existing upper bound of the corresponding process. Therefore in
our discussion of available DM parameter space in Section 5.1, we ignore the
indirect search constraints, which are automatically satisfied once direct
search constraints take into account.

###  3.2 LHC constraints

Finally, we discuss the phenomenology of the Higgs sector at the LHC in this
subsection. The production of the triplet scalars at the LHC via $s$-channel
exchange of $Z/\gamma^{*}$ and $W^{\pm}$ lead to various final states of the
SM particles from their prompt decays. The wide array of phenomenological
consequences of type-II seesaw at the LHC has been thoroughly studied in the
literature Huitu:1996su ; Chakrabarti:1998qy ; Chun:2003ej ; Akeroyd:2005gt ;
Garayoa:2007fw ; Kadastik:2007yd ; Akeroyd:2007zv ; FileviezPerez:2008jbu ;
delAguila:2008cj ; Akeroyd:2009hb ; Melfo:2011nx ; Aoki:2011pz ;
Akeroyd:2011zza ; Chiang:2012dk ; Chun:2012jw ; Akeroyd:2012nd ; Chun:2012zu ;
Dev:2013ff ; Banerjee:2013hxa ; delAguila:2013mia ; Chun:2013vma ;
Kanemura:2013vxa ; Kanemura:2014goa ; Kanemura:2014ipa ; kang:2014jia ;
Han:2015hba ; Han:2015sca ; Das:2016bir ; Babu:2016rcr ; Mitra:2016wpr ;
Cai:2017mow ; Ghosh:2017pxl ; Crivellin:2018ahj ; Du:2018eaw ; Dev:2018kpa ;
Antusch:2018svb ; Aboubrahim:2018tpf ; deMelo:2019asm ; Primulando:2019evb ;
Padhan:2019jlc ; Chun:2019hce ; Ashanujjaman:2021txz ; Mandal:2022zmy ;
Dutta:2014dba . Numerous collider searches have been carried out by the CMS
and ATLAS collaborations to hunt for the same at the LHC ATLAS:2012hi ;
Chatrchyan:2012ya ; ATLAS:2014kca ; Khachatryan:2014sta ; CMS:2016cpz ;
CMS:2017pet ; Aaboud:2017qph ; CMS:2017fhs ; Aaboud:2018qcu ; Aad:2021lzu .
However, there has not been any evidence of an excess over the SM background
expectations so far. Thus, these searches strongly constrain the parameter
space relevant to the type-II seesaw model.

Figure 2: [Left] Branching ratios of doubly charged scalar, $H^{\pm\pm}$ as a
function of $v_{t}$ for $m_{H^{\pm\pm}}=300$ GeV and $\Delta m=15$ GeV.
[Right] Phase diagram for $H^{\pm\pm}$ decays to different modes is shown in
$\Delta m-v_{t}$ plane by different contour lines: the leptonic decay mode by
blue dotted line, the gauge boson decay mode by orange dotted line and the
cascade decay mode by the blue line. Each contour line corresponds to $99\%$
of the branching ratio for the corresponding decay channel. The mass of
$H^{\pm\pm}$ is kept fixed as $m_{H^{\pm\pm}}=300$ GeV. The cascade dominated
region ($Br(H^{\pm\pm}\to{W^{\pm}}^{*}H^{\pm})>99\%$), which we are interested
in our discussion is shown by the shaded orange region.

The most tantalizing collider signature of a type-II seesaw motivated model is
the production and decay of $H^{\pm\pm}$. It is well known that $H^{\pm\pm}$
can decay to $\ell^{\pm}\ell^{\pm}$ ($W^{\pm}W^{\pm})$ for low (relatively
large) values of $v_{t}$, when the triplet scalars are degenerate in mass. If
one introduces a mass-splitting between the triplet scalars, additionally, the
$H^{\pm\pm}\rightarrow H^{\pm}W^{\pm}$ decay channel opens up for $\Delta m>0$
scenarios. However, as we have seen in Section F.3, the EW precision
observables limit the mass-splitting, $\Delta m$, to be less than 40 GeV.
Hence, this restricts $H^{\pm\pm}$ to decay to either $H^{\pm}\ell^{\pm}\nu$
and $H^{\pm}q\bar{q^{\prime}}$ via an off-shell $W^{\pm*}$. For $\Delta
m\lesssim\mathcal{O}(2)$ GeV, when quarks cease to remain as free particles,
then $H^{\pm\pm}$ decay to $H^{\pm}\pi^{\pm}$ has to be considered also.
Similar behavior can be observed for decays of $H^{\pm}$ as well, again for
$\Delta m>0$, leading to a cascade decay of $H^{\pm\pm}$ in turn. In fact, for
$v_{t}\sim 10^{-3}-10^{-6}$ GeV, the $H^{\pm\pm}(H^{\pm})\rightarrow
H^{\pm}(H^{0}/A^{0})W^{\pm*}$ decay becomes the dominant decay channel for
$\Delta m\gtrsim 10$ GeV. For a detailed discussion on the various decay
channels of $H^{\pm\pm}(H^{\pm})$ we refer interested readers to references
FileviezPerez:2008jbu ; Ghosh:2018drw . To illustrate the above points we show
BFs of $H^{\pm\pm}$ to $\ell^{\pm}\ell^{\pm}$ ($W^{\pm}W^{\pm})$, and the
cascade channel ($H^{\pm}W^{\pm*}$) for $m_{H^{\pm\pm}}=300$ GeV with $\Delta
m=15$ GeV in the left panel of figure 2. In contrast, we display the
Br$(H^{\pm\pm}\rightarrow H^{\pm}W^{\pm*})>99\%$ region by the shaded area in
the right panel of figure 2 as a function of $v_{t}$ as well as $\Delta m$.

In the $\Delta m=0$ case, the ATLAS collaboration places the most stringent
bound of $m_{H^{\pm\pm}}\gtrsim 870$ GeV, assuming
Br($H^{\pm\pm}\rightarrow\mu^{\pm}\mu^{\pm}$) to be 100$\%$ Aaboud:2017qph ,
and using 36.1 fb-1 data at $\sqrt{s}=13$ TeV. On the other hand, ATLAS again
excludes $m_{H^{\pm\pm}}\gtrsim 350$ GeV in the $W^{\pm\pm}$ final state with
139 fb-1 of data Aad:2021lzu . The first bound is applicable in the small
$v_{t}$ regime, while the second is for relatively large $v_{t}$ values. It is
clear from the discussion in the previous paragraph that the above limits are
not applicable to the whole $\Delta m-m_{H^{\pm\pm}}$ plane for the range of
$v_{t}$, where the cascade decay have a significant BF.

Let us first consider $\Delta m<0$ scenarios. For a large mass-splitting
($\mathcal{O}(10)$ GeV) and a moderate $v_{t}$, cascade decays of
$H^{0}/A^{0}\to H^{\pm}W^{\mp*}$ and $H^{\pm}\to H^{\pm\pm}W^{\mp*}$, dominate
over other decays. In this regime, $H^{\pm\pm}$ is the lightest triplet
scalar, and decays to either $\ell^{\pm}\ell^{\pm}$ or $W^{\pm}W^{\pm}$
depending on the choice of $v_{t}$. Hence, the effective pair production
cross-section of $H^{\pm\pm}$ is augmented by cascade decays of $H^{0},A^{0}$
and $H^{\pm}$. So, one can expect the bounds on $m_{H^{\pm\pm}}$ to be strong
for $\Delta m<$ 0 scenarios compared to the $\Delta m=0$ case. In a recent
analysis, reference Ashanujjaman:2021txz  has shown by recasting various ATLAS
and CMS analyses that for $\Delta m=-10$ ($-30$) GeV the present exclusion
limit for $m_{H^{\pm\pm}}$ is 1115 (1076) GeV with $v_{t}\sim 10^{-5}-10^{-6}$
GeV. In the same paper, the authors have improved the bound on
$m_{H^{\pm\pm}}$ for $\Delta m=0$ to 955 GeV (420 GeV) for small (relatively
large) $v_{t}$.

On the other hand, in $\Delta m>0$ scenarios, $H^{\pm\pm}$ is the heaviest
triplet member, and again for a mass-splitting of $\mathcal{O}(10)$ GeV and a
moderate $v_{t}$, $H^{\pm\pm}$ will mostly undergo the cascade decay,
$H^{\pm\pm}\rightarrow W^{\pm*}H^{\pm}\rightarrow
W^{\pm*}W^{\pm*}\,H^{0}/A^{0}$. $H^{0}$ and $A^{0}$ mostly decay to invisible
neutrinos for $v_{t}<$ $10^{-4}$ GeV. $H^{0}$ can also decay to $b\bar{b}$ or
$\tau^{+}\tau^{-}$ if its mixing angle with the SM Higgs is sufficiently
large. The EW precision bound of $\Delta m<40$ GeV ensures that the leptons
and jets arising from $H^{\pm\pm}$ cascade decay will be soft and may not
cross the $p_{T}^{\ell/j}$ thresholds required in canonical
$\ell^{\pm}\ell^{\pm}$ and $W^{\pm}W^{\pm}$ searches for $H^{\pm\pm}$.
Reference Ashanujjaman:2021txz  shows that even with 3 ab-1 of integrated
luminosity data, the LHC will be unable to constrain $m_{H^{\pm\pm}}\gtrsim$
200 GeV for $\Delta m=$ 10 GeV (30 GeV), with $v_{t}$ around $10^{-3}$ GeV to
$10^{-5}$ GeV ($10^{-3}$ GeV to $10^{-6}$ GeV). This scenario is akin to
compressed supersymmetry spectra and a dedicated search strategy is required
to develop to probe this region Baer:2014kya ; Dutta:2012xe ; Dutta:2014jda ;
Ajaib:2015yma ; Dutta:2017nqv .

In addition, the measurement of the properties of the observed 125 GeV Higgs
boson can also restrict some of our model parameters. As the observed Higgs
boson’s couplings to the SM particles are found to be almost SM-like, it will
not allow a significant mixing between the neutral $CP$-even components of the
triplet and the SM doublet. In order to satisfy the Higgs signal strength
bounds, the mixing angle $\sin\theta_{t}$, defined in equation 7, should be
below 0.1 ATLAS:2016neq ; CMS:2018uag . The observed Higgs can decay to a pair
of DM scalars if $m_{S}<m_{h^{0}}/2$. Hence, we impose the limit
Br($h^{0}\rightarrow$ invisible) $<$ 0.15 on our model parameters CMS:2018yfx
.

Furthermore, the decay width of the SM Higgs into a photon pair can be
affected by the presence of $H^{\pm}$ and $H^{\pm\pm}$ in the loops.
$\Gamma(h^{0}\rightarrow\gamma\gamma)$ can both be enhanced or reduced
compared to the SM value depending on the sign of quartic couplings
$\lambda_{1}$ and $\lambda_{4}$. Since the partial decay width of the SM into
photons is tiny, the total decay width of $h^{0}$ barely changes. The signal
strength parameter of $h^{0}\rightarrow\gamma\gamma$ channel is defined as,

$\mu_{\gamma\gamma}=\frac{\Gamma^{\text{N}P}[h^{0}\rightarrow\gamma\gamma]}{\Gamma^{\text{S}M}[h^{0}\rightarrow\gamma\gamma]}$
(50)

where $\Gamma^{\text{SM (NP)}}[h^{0}\rightarrow\gamma\gamma]$ denotes the
decay rate without (with) the inclusion of new physics. We consider that the
production cross-section of $h^{0}$ will remain the same when the effect of NP
is included since $\cos\theta_{t}\rightarrow 1$. For more details on
$\Gamma(h^{0}\rightarrow\gamma\gamma)$ within the type-II seesaw model we
refer the reader to references Arhrib:2011vc ; Gunion:1989we . The current
limits on $\mu_{\gamma\gamma}$ from ATLAS and CMS are $1.04^{+0.10}_{-0.09}$
ATLAS:2022tnm  and $1.12{\pm 0.09}$ CMS:2021kom , respectively.

##  4 Choice of parameter space

###  4.1 Implications from Higgs searches

As discussed in subsection 3.2, LHC direct searches for $H^{\pm\pm}$ will not
be able to limit the triplet-like scalar masses even at 200 GeV for a moderate
$v_{t}$ and relatively large $\Delta m$ ($>10$ GeV) Ashanujjaman:2021txz .
Depending on the value of $v_{t}$ and mixing with $h^{0}$, in this scenario,
$H^{0}/A^{0}$ either decays into $\nu\nu/b\bar{b}/\tau^{+}\tau^{-}$ or
$h^{0}h^{0},ZZ/h^{0}Z$. Henceforth, in this work, we focus on 10 GeV $<\Delta
m<$ 40 GeV and $10^{-3}$ GeV $<v_{t}<$ $10^{-6}$ GeV region of parameter space
of the type-II seesaw part of the model. On the other hand, the rich scalar
sector of this model can generate stochastic GW signals from FOPT in the early
Universe. As the LHC fails to constrain this region from the direct searches
of triplet-like scalars, we try to probe this region of parameter space from
GW experiments.

For $\Delta m>0$, $\lambda_{4}$ needs to be $-$ve (see equation 13). Thus, we
consider $\lambda_{4}<0$ region of parameter space in our analysis. For
successful FOPT to take place in the Higgs sector, the two most important
quartic couplings are $\lambda_{1}$ and $\lambda_{4}$. They are responsible
for mixing the doublet and the triplet fields (see quation 2.3.1). We vary the
absolute value of these couplings from 0 to 4. Other quartic couplings of the
triplet-sector $\lambda_{2}$ and $\lambda_{3}$ do not significantly affect
FOPT. Thus, we set them to fixed values for the rest of our study. Since we
focus on light triplet-like scalars, we require $m_{H^{\pm\pm}}$ within 200
GeV to 600 GeV, and for that purpose, we vary $M_{\Delta}$ from 100 GeV to 800
GeV. In addition, for the above choice of our $v_{t}$ and $m_{H^{\pm\pm}}$
values, our benchmark points pass the flavour constraints discussed in
subsection F.2 as well.

The precision measurements of $h^{0}$ couplings at the LHC restrict the mixing
between the doublet-like and the triplet-like $CP$-even neutral scalars to
small values. Therefore, we limit our scan to tiny $\sin\theta_{t}$ values to
allow minimum mixing and maintain the SM-like nature of $h^{0}$. However,
loop-induced $h^{0}$ decays, such as $h^{0}\rightarrow\gamma\gamma$, can
deviate significantly from the observed limits even at small $\sin\theta_{t}$
values. We discuss the variation of the signal strength, $\mu_{\gamma\gamma}$
of the $h^{0}\rightarrow\gamma\gamma$ channel in the above-discussed region of
parameter space below.

Figure 3: Variation of the signal strength $\mu_{\gamma\gamma}$ in the
$\lambda_{1}-\lambda_{4}$ plane. The other two parameters, $M_{\Delta}$ and
$v_{t}$, of the triplet sector that can alter the signal strength are fixed at
300 GeV and $10^{-5}$ GeV, respectively.

At the one-loop level, the top quark, along with the charged weak gauge bosons
and the charged triplet-like Higgs states contribute to the decay of
$h^{0}\rightarrow\gamma\gamma$. The deviation of the signal strength of this
decay channel $\mu_{\gamma\gamma}$, defined in equation 50, from 1 can become
a signature of new physics beyond the SM. We present that the variation of
$\mu_{\gamma\gamma}$ for $\lambda_{1}>0$, $\lambda_{4}<0$, $M_{\Delta}=300$
GeV and $v_{t}=10^{-4}$ GeV in figure 3. The variation of $\lambda_{1}$ and
$\lambda_{4}$ are presented in the $x$-axis and $y$-axis, respectively,
whereas the variation of $\mu_{\gamma\gamma}$ is indicated by palette colors.
The figure also takes into consideration the theoretical constraints. Due to
constraints from vacuum stability that are discussed in section
E.1,444Particularly from the inequality relations:
$(\lambda_{1}+\lambda_{4})+2\sqrt{\lambda_{H}(\lambda_{2}+\lambda_{3})}\geq 0$
and
$(2\lambda_{1}+\lambda_{4})+4\sqrt{\lambda_{H}(\lambda_{2}+\frac{\lambda_{3}}{2})}\geq
0$ the white empty section of figure 3 is excluded. We note that due to a
destructive interference between the SM and BSM contributions to
$\mu_{\gamma\gamma}$, it is mostly less than 1 for $\lambda_{4}<0$. The
deviation of $\mu_{\gamma\gamma}$ from 1 increases with the increase of the
absolute values of $\lambda_{1}$ and $\lambda_{4}$. Thus the precision study
of $\mu_{\gamma\gamma}$ at the LHC and other proposed collider experiments can
exclude a significant amount of the parameter space of our interest in this
paper. We shall discuss this issue further in section 5.

###  4.2 Implications from DM searches

In contrast to the Higgs sector discussed above, in the dark sector, we vary
$m_{S}$, $\lambda_{SH}$, and $\lambda_{S\Delta}$ within a certain range to
satisfy all DM constraints. $\lambda_{SH}$ coupling has a significant impact
on FOPT as it mixes the singlet $h_{s}$ field with the $h_{d}$-field (see
equation 2.3.1). Thus, the possibility of a strong FOPT is expected to
increase with a large $\lambda_{SH}$. However, as we have already discussed in
section 2.2, increasing the same coupling increases $\sigma_{\rm DD}^{\rm SI}$
as well. Hence, large $\lambda_{SH}$ values are ruled out by DMDD limits. We
shall discuss this issue in more detail in section 5.3. Another dark sector
parameter, $\mu_{3}$, also impacts both DM and the FOPT. This coupling not
only can produce the necessary barrier for FOPT but also contributes to DM
semi-annihilation processes. We vary this parameter in our study satisfying
the constraints $\mu_{3}/m_{S}\lesssim 2\sqrt{\lambda_{S}}$ from the stability
of the global minimum, which is discussed in subsection E.1.

To summarize, the ranges of input parameters of the model we consider for the
scan are listed in table 2. The other parameters are fixed at
$\lambda_{H}=0.129$, $\lambda_{2}=0.1$, $\lambda_{3}=0$, $\lambda_{S}=0.1$.

Varying parameters | $\lambda_{1}$ | $\lambda_{4}$ |  $M_{\Delta}$ (GeV) |  $v_{t}$ (GeV) |  $m_{S}$ (GeV) |  $|\lambda_{SH}|$ |  $|\lambda_{S\Delta}|$ |  $\mu_{3}$ (GeV) |  $|\lambda_{S}|$
---|---|---|---|---|---|---|---|---|---
Ranges | 0 – 4 |  $-4$ – 0 | 100 – 800 |  $10^{-6}$ – $10^{-3}$ | 0 – 1000 | $<4$ | $<6$ | $<1000$ | $<6$
Table 2: Ranges of various model parameters employed for scanning the present
model parameter space.

##  5 Results

In this section, we investigate the DM phenomenology, the thermal history
regarding an FOEWPT, the associated production of stochastic GW and the
interplay among the DM, GW and the collider physics in this model. In the
subsection 5.1, we analyse the dependencies of the DM phenomenology in terms
of both DM relic density and DMDDSI cross-section on the triplet-sector scalar
masses. We study the variation of various model parameters of the dark sector
in the estimation of relic density and DMDDSI cross-section. Implications of
the various model parameters on FOEWPT have been studied in two parts. The
impact of the triplet-sector parameters on the FOEWPT, the production of GW
and the interplay with LHC have been studied in the subsection 5.2. Two
benchmark scenarios are presented to illustrate the effect of the triplet-
sector parameters on the FOEWPT. In the subsection 5.3, we study the
possibilities of FOEWPT by altering the shape of the SM Higgs potential mainly
from the dark sector. The connection between the DM observables and an FOPT
has been discussed. A benchmark scenario that exhibits a two-step phase
transition in the early Universe has been presented to illustrate such a
scenario.

###  5.1 DM Parameter space

In this section, we discuss the DM parameter space in terms of model
parameters which is consistent with DM relic density observed by the PLANCK
experiment and the upper bound on DM-nucleon spin-independent cross-section
provided by XENON-1T and PANDA-4T experiments. Before going to the details of
the parameter space scan, we attempt to comprehend the variation of DM relic
density with DM mass, $m_{S}$ $(\equiv m_{\rm DM})$ and other relevant
parameters of the model such as scalar portal couplings $\lambda_{SH}$,
$\lambda_{S\Delta}$ and trilinear coupling $\mu_{3}$. The abundance of DM also
significantly depends on the masses of triplet scalar, i.e.,
$m_{H^{0}},~{}m_{A^{0}},~{}m_{H^{\pm}}$ and $m_{H^{\pm\pm}}$ which can be
expressed as the independent parameters as mentioned in equation 12. It is
difficult to establish the significance of all the relevant parameters that
vary simultaneously. Therefore, we fix the triplet sector parameters (in
equation 12), which simplifies the setup and allows us to investigate DM
phenomenology. For DM discussion, we consider one benchmark point from the
cascade region of the triplet sector (as discussed in section 4.2) which is
given by:

$\displaystyle\\{M_{\Delta}=367.7~{}{\rm GeV},~{}v_{t}=3.3\times
10^{-4}~{}{\rm
GeV},~{}\lambda_{1}=3.92,~{}\lambda_{2}=0.1,~{}\lambda_{3}=0,~{}\lambda_{4}=-0.989\\}$

which corresponds to the physical masses and mixing angle followed from
equations LABEL:MCPE, 9, 10 and 11 as:

$\displaystyle\\{m_{H^{0}}=367.7,~{}m_{A^{0}}=367.7,~{}m_{H^{\pm}}=387.5,~{}m_{H^{\pm\pm}}=406.4~{}{\rm
GeV},$ $\displaystyle~{}\mu=0.00104,~{}\sin\theta_{t}=-1.03\times
10^{-6}~{}\\}$ (52)

It is worth mentioning here that one can also consider another set of
parameters as mentioned above, for which the underlying physics remains the
same. The only thing that will change is the allowed region of DM mass
depending on the triplet scalar mass. Therefore the phenomenology of DM
depends on the following independent parameters in the dark sector along with
the above-fixed parameters

$\displaystyle\\{m_{S},~{}\mu_{3}/m_{S}=1,~{}\lambda_{SH},~{}\lambda_{S\Delta}\\}.$
(53)

Note that for DM discussion we consider here $\mu_{3}/m_{S}=1$ which is
responsible for the semi-annihilation process, $S~{}S\to S~{}h^{0}$. The
larger values of $\mu_{3}/m_{S}$ lead to a larger semi-annihilation
contribution to DM abundance. However, a stable EW vacuum sets an upper bound
on the trilinear coupling, $\mu_{3}\leq 2m_{S}$ as discussed earlier.

Figure 4: Variation of relic density as a function of DM mass for different
ranges $\lambda_{S\Delta}$ : $0.01\leq\lambda_{S\Delta}\leq 0.04$ (red),
$0.04<\lambda_{S\Delta}\leq 0.1$ (blue), $0.1<\lambda_{S\Delta}\leq 0.3$
(green) and $0.3<\lambda_{S\Delta}\leq 0.5$ (orange). We illustrate two
different choices of $\lambda_{SH}:$ $\\{0.01-0.02\\}$ (left) and
$\\{0.10-0.11\\}$ (right).

In Fig.4, we show the variation of DM abundance ($\Omega_{S}h^{2}$) as a
function of $m_{S}$ for two different choices of $\lambda_{SH}$:
$0.01\leq\lambda_{SH}\leq 0.02$ (left) and $0.10\leq\lambda_{SH}\leq 0.11$
(right). Different colored patches indicate different ranges of DM portal
coupling with triplet, $\lambda_{S\Delta}$ : $0.01\leq\lambda_{S\Delta}\leq
0.04$ (red), $0.04<\lambda_{S\Delta}\leq 0.1$ (blue),
$0.1<\lambda_{S\Delta}\leq 0.3$ (green) and $0.3<\lambda_{S\Delta}\leq 0.5$
(orange) are considered for each plot. The black dotted horizontal lines
indicate the observed DM relic density from the PLANCK data Planck:2018vyg ,
i.e., $\Omega_{\text{DM}}h^{2}=0.120$. With the increase in DM mass, different
number-changing processes open up after crossing some threshold values
resulting in a drop in relic density.

$\bullet$ $m_{S}<m_{h^{0}}$: When DM mass, $m_{S}<m_{h^{0}}$, the DM abundance
is mostly dominated by the $S~{}S\to{\rm SM~{}SM}$ number changing process
which is mediated by both the CP even Higgses, $h^{0}$ and $H^{0}$. And the
density of DM varies as:

$\displaystyle\Omega_{S}h^{2}\propto\frac{1}{\langle\sigma v\rangle_{SS\to{\rm
SM~{}SM}}}.$ (54)

where for the small $\sin\theta_{t}$ limit, $\langle\sigma v\rangle_{SS\to{\rm
SM~{}SM}}\propto{\lambda_{SH}^{2}}/{m_{S}^{2}}$ and almost independent of
$\lambda_{S\Delta}$ which is shown from the above Figs.4. The annihilation
contribution will get suppressed with the increase of DM mass and hence relic
density increases. For a fixed DM mass, with an increase of $\lambda_{\phi
H}$, relic density decreases as it is shown in the right panel figure. A sharp
drop in DM density due to the SM Higgs, ${h^{0}}$ resonance near $m_{S}\sim
m_{h^{0}}/2$. Beyond the Higgs($h^{0}$) pole, $m_{S}>m_{W},m_{Z}$, DM
annihilates into a pair gauge final states resulting in a large enhancement in
$\langle\sigma v\rangle$, which leads to less density.

$\bullet$ $m_{h^{0}}\leq m_{S}<m_{H^{0}}$: A new annihilation channel,
$S~{}S\to S~{}h^{0}$ (semi-annihilation)starts contributing when the DM mass
becomes heavier than Higgs mass ($m_{S}>m_{h^{0}}$) thanks to the trilinear
coupling, $\mu_{3}$ present in this scenario. Therefore the relic density of
DM in this region is followed as:

$\displaystyle\Omega_{S}h^{2}\propto\frac{1}{\langle\sigma v\rangle_{SS\to{\rm
SM~{}SM}}+\langle\sigma v\rangle_{SS\to Sh^{0}}}.$ (55)

The second term in the above equation depends on $\mu_{3}$ and $\lambda_{SH}$
whereas the first term only depends on $\lambda_{SH}$. The semi-annihilation
contributions dominate over the standard annihilation to SM for a larger value
of $\mu_{3}$ which will help to evade the direct search bound on
$\lambda_{SH}$. With the increase of $\lambda_{SH}$, both standard
annihilation and semi-annihilation contributions are significantly enhanced
and hence relic density decreases as it is depicted from the left and right
panel of Fig.4. Note here that the effect of semi-annihilation becomes
pronounced where the propagator suppression is minimal near $m_{S}\gtrsim
m_{h^{0}}$. There is no variation of relic density with $\lambda_{S\Delta}$
because both the contributions are almost insensitive to the change of
$\lambda_{S\Delta}$ for the small $\sin\theta_{t}$ limit. There is no drop in
relic density near $m_{S}\sim m_{H^{0}}/2$ due to the second heavy Higgs,
$H^{0}$. This is because of the heavy Higgs mediated diagrams, $S~{}S\to{\rm
SM~{}SM}$, are strongly suppressed by small $\sin\theta_{t}$ and $v_{t}$.

$\bullet$ $m_{S}\gtrsim m_{H^{0}}$: The new annihilation processes start to
contribute to DM relic when DM mass is larger than the masses of the triplet
states i.e. $2m_{S}>m_{X}+m_{Y}$ and the relic density of DM turns out to be

$\displaystyle\Omega_{S}h^{2}\propto\frac{1}{\langle\sigma v\rangle_{SS\to{\rm
SM~{}SM}}+\langle\sigma v\rangle_{SS\to Sh^{0}}+\langle\sigma v\rangle_{SS\to
X~{}Y}}~{},~{}\\{X,Y\\}:=\\{H^{0},A^{0},H^{\pm},H^{\pm\pm}\\}.$

The Higgs mediated s and t channel diagrams for $SS\to XY$ process are
suppressed by a small mixing angle, $\sin\theta_{t}$ and small $v_{t}$. But
the four-point contact diagrams for the $SS\to XY$ process enhance the cross-
section significantly with the increase of $\lambda_{S\Delta}$. Therefore the
effect of new annihilation to DM relic density becomes more prominent with the
increase of $\lambda_{S\Delta}$. Since this coupling is almost insensitive to
$\sigma_{\rm DD}^{\rm SI}$ for the small mixing angle, one can play with it to
satisfy observed DM density. For the smaller values of $\lambda_{SH}$ and
$\mu_{3}(=m_{S})$, the most dominant contribution to DM relic is coming from
the new annihilation channels, $S~{}S\to X~{}Y$ whereas the other two
contributions, $S~{}S\to{\rm SM,SM}$ and $S~{}S\to S~{}h^{0}$ are sub-
dominate. With the increase of $\lambda_{S\Delta}$, $\langle\sigma
v\rangle_{S~{}S\to X,Y}$ increases and hence relic density decreases which are
shown from the left panel of Fig.4 for fixed choices of
$\lambda_{SH}:\\{0.01-0.02\\}$ and $\mu_{3}/m_{S}=1$. In the right panel of
Fig.4, we considered comparably large values of $\lambda_{SH}:\\{0.10-0.11\\}$
which further increases the $S~{}S\to{\rm SM~{}SM}$ and $S~{}S\to S~{}h^{0}$
contributions and hence total effective annihilation cross-section increases.
As a result, relic density decreases compared with the left panel figure. Note
that a new DM semi-annihilation, $SS\to SH^{0}$ is kinematically allowed here
but it has almost negligible contribution to relic density as the cross-
section is suppressed by the small mixing angle, $\sin\theta_{t}$ and small
$v_{t}$.

Figure 5: [Left]SI spin-independent DM-nucleon scattering cross-section,
($\sigma^{\rm SI}_{\rm DD}$) for relic density allowed(PLANCK) parameter space
as a function of DM mass($m_{S}$). For comparison, the current upper bound on
the SI DM-nucleon cross-section from XENON 1T Aprile:2018dbl  and PANDAX 4T
PandaX-4T:2021bab  data are shown in the same plane. The Neutrino Floor is
represented by a shaded orange region. [Right] Relic(PLANCK)$+$DD(PANDAX-4T)
allowed parameter space is shown in $m_{S}-\lambda_{SH}$ plane for different
ranges of $\lambda_{S\Delta}$ shown by different color patches. Note that the
relic density of DM, $S$ considered here $100\%$ ($\xi_{\rm DM}=1$) of
observed DM density.

We shall now move to the parameter space scan. To get the allowed region of
parameter space, we perform a three-dimensional random scan over the following
parameters keeping other parameters fixed as discussed above.

$\displaystyle\\{m_{S}:~{}\\{30-1000\\}{\rm
GeV},~{}~{}\lambda_{SH}:\\{0.001-0.3\\},~{}~{}\lambda_{S\Delta}:\\{0.01-0.5\\}\\}~{}.$
(57)

In the left panel of Fig.5, we show the relic density allowed parameter space
(with $\xi_{\rm DM}=1$) emerging from the random scan in the plane of
$\sigma_{\rm DD}^{\rm SI}$ versus DM mass, $m_{S}$. As discussed in section
2.2, for a small mixing limit($\sin\theta_{t}\to 0$), the DD cross-section is
proportional to $\lambda_{SH}^{2}$ and almost independent of
$\lambda_{S\Delta}$. Therefore with an increase of $\lambda_{SH}$, the DD
cross-section increases. The latest upper bounds on $\sigma_{\rm DD}^{\rm SI}$
against DM mass from XENON 1T (purple dashed line) and PANDAX 4T (black dashed
line) are shown in the same plane. The parameter space above the dashed line
can be disfavoured from the non-observation of the DM signal of the
corresponding direct search experiments. We find that DM mass below the mass
of the lightest triplet state ($m_{S}<m_{H^{0}}$), the DM-SM Higgs portal
coupling, $\lambda_{SH}$ is required comparably large to satisfy observed
relic density which corresponds to large DD cross-section lies above the
experimental upper limit. As a result, most of the region($m_{S}<m_{H^{0}}$)
is excluded from non-observation of DM signal apart from Higgs poles,
$m_{S}\sim m_{h^{0}}/2$. Once the triplet final states open up, the DM portal
coupling with a triplet, $\lambda_{S\Delta}$ takes up the major role in
controlling DM relic density. With the help of $\lambda_{S\Delta}$ which is
insensitive to DD-cross-section, the pressure on $\lambda_{SH}$ can be reduced
to produce correct relic density which can evade the direct search bound. This
phenomenon can be observed from the right panel of Fig.5. As a result, the
region above the DM mass, $m_{S}>m_{H^{0}}$ is allowed from both relic and
direct search constraints with the help of large $\lambda_{S\Delta}(\gtrsim
0.1)$.

###  5.2 Implications of the triplet sector parameters on SFOEWPT

We have already discussed how the discovery of a 125 GeV Higgs mass in the SM
demonstrates that electroweak breaking occurs through a smooth cross-over
transition. However, it can be an FOPT in the presence of additional scalars,
going Beyond-the SM. A further benefit of FOEWPT is that it can supply one of
the necessary conditions for the explanations of BAU. In this section, we
illustrate the possibility of FOPT in the electroweak sector in our choice of
parameter region of the current model. To explain the observed matter-
antimatter asymmetry via the mechanism of EWBG, the FOPT in the electroweak
sector needs to satisfy another additional condition, i.e.,

$\xi_{n}=\frac{v_{n}}{T_{n}}>1.$ (58)

$v_{n}$ is the ${\it vev}$ in the electroweak minimum at the nucleation
temperature $T_{n}$. The strength of the FOPT is quantified via $\xi_{n}$.
This criterion is required to prevent the wash-out of the generated baryon
asymmetry after the EWPT. We present four plots of the FOPT-allowed scan
results in figure 6 to illustrate the correlations among the model parameters
of the triplet-sector in connection with the FOPT and its strength, $\xi_{n}$.
Here, We fix the various couplings of dark sectors to small values to ensure
that their impact on the FOPT is minimal and we vary $\lambda_{1}$,
$\lambda_{4}$, $v_{t}$ and $M_{\Delta}$ only in the range as described in
table 2. Thus, here we want to study the impact of the triplet-sector
parameters on the FOPT in the $h_{d}$-field direction. In our scan results
$v_{n}\sim v_{h}$, since the triplet ${\it vev}$ always remains very small.

Figure 6: Scattered points that exhibits FOPT along the $h_{d}-$field
direction in the $\lambda_{1}-\lambda_{4}$ plane with $M_{\Delta}$ (strength
of the phase transition, $\xi_{n}$) indicated by the palette in the top left
(bottom right) plot. The variation of the same points in the
$m_{H^{\pm\pm}}-\lambda_{1}$ ($m_{H^{\pm\pm}}-\lambda_{4}$) plane with the
palette indicating the strength of the phase transition $\xi_{n}$ is shown in
the top right (bottom left) plot.

Figure 7: Left: Scatter points in the $(\lambda_{1}+\lambda_{4})-M_{\Delta}$
plane that satisfies FOPT along the $h_{d}-$field direction with the palettes
indicating the strength of the phase transition, $\xi_{n}$. Right: Similar to
the scatter spots on the left plot but in the
$(\lambda_{1}+\lambda_{4})-T_{n}$ plane with showing the variation of
$M_{\Delta}$ via the palette.

The plot on top, left of figure 6 shows the FOPT allowed region of
$\lambda_{1}$ and $\lambda_{4}$ with the $M_{\Delta}$ being indicated by the
palette-color. It shows that FOPT prefers relatively larger $\lambda_{1}$. We
have already discussed that in our choice of parameter region ($\Delta m>0$),
$\lambda_{4}$ has to be $-$ve and from the theoretical constraints
(particularly constraints from the stability) that we have discussed in
section E, the absolute value of $\lambda_{4}$ cannot be too large compared
with $\lambda_{1}$. As a result, the permitted range of $\lambda_{4}$ rises
with larger $\lambda_{1}$. Color variation reveals that $M_{\Delta}$ needs to
be on the smaller side at the relatively smaller values of $\lambda_{1}$. From
equation 11, $m_{H{{}^{\pm}\pm}}\sim M_{\Delta}$ at smaller $\lambda_{4}$.
Thus, FOPT demands relatively light triplet-like scalars at smaller quartic
coupling $\lambda_{1}$. At the relatively larger triplet-like scalar masses
(larger $M_{\Delta}$), FOPT demands larger $\lambda_{1}$. On the other hand,
for a given $\lambda_{1}$, the decrease of the absolute value of $\lambda_{4}$
demands increasing $M_{\Delta}$ for FOPT.

The shape of the SM Higgs potential gets modification from the quartic terms
that are proportional with $\lambda_{1}$ and $\lambda_{4}$ in equation 2.3.1.
Thus, the phase transition pattern in the Higgs sector strongly depends on
these quartic couplings $\lambda_{1}$ and $\lambda_{4}$. In the plot on top,
right and bottom, left of figure 6, we present the variation of
$m_{H^{\pm\pm}}$ with $\lambda_{1}$ and $\lambda_{4}$, respectively, from our
scanned points that satisfy FOPT. The palette-color indicates the associated
strength of the FOPT, $\xi_{n}$, values. For a fixed $m_{H^{\pm\pm}}$, color
variation shows that the strength of the FOPT increases with $\lambda_{1}$.
Thus, larger $\lambda_{1}$ is favoured for an SFOEWPT, which is required for
EWBG. In addition, for a fixed value of the quartic coupling $\lambda_{1}$,
lighter $m_{H^{\pm\pm}}$ is preferred for an SFOEWPT. Thus, larger quartic
coupling $\lambda_{1}$ and lower triplet-like scalars increase the strength of
the FOPT. The other plot indicates that for a given $m_{H^{\pm\pm}}$, the
strength of the FOPT increases with the decrease of the absolute value of
$\lambda_{4}$. The dependence of $\lambda_{1}$ and $\lambda_{4}$ on $\xi_{n}$
can clearly be understood from the bottom and the right plot of figure 6. It
shows that the increase of $\lambda_{1}$ and the decrease of the absolute
value of $\lambda_{4}$ increases the strength of the phase transition.

The effective dependence of these quartic couplings on the phase transition
pattern can be quantified in terms of the sum of these two quartic couplings,
i.e., $\lambda_{1}+\lambda_{4}$. We present the variation of
$\lambda_{1}+\lambda_{4}$ with $M_{\Delta}$ with $\xi_{n}$ being indicated by
the palette-color in the left plot of figure 7. The variation of colors
indicates that the strength of FOPT increases with the increase of the
effective quartic coupling ($\lambda_{1}+\lambda_{4}$) and decrease of
$M_{\Delta}$ (or, the triplet-like scalar masses).

Another very important quantity of the FOPT is the temperature at which the PT
starts to occur, i.e., the nucleation temperature ($T_{n}$). We have discussed
in section 2.3.2 that the production of GW peak amplitude and the peak
frequency from an FOPT depend on $T_{n}$. Thus, studying the variations of
$T_{n}$ with the various model parameters is important to identify the region
of parameter space for a strong GW signal from FOPT. Keeping this in mind, We
show the variation of $T_{n}$ with respect to the effective quartic coupling
($\lambda_{1}+\lambda_{4}$) in the right plot of figure 7. The variation of
$M_{\Delta}$ is presented in palette color. It can be seen that at relatively
small effective quartic coupling, $T_{n}$ is relatively large. Also, low
effective quartic coupling demands relatively low $M_{\Delta}$. On the other
hand, at relatively larger effective quartic coupling, $T_{n}$ can go down to
even below 50 GeV at moderate values of $M_{\Delta}$.

Finally, after studying all these plots in figures 6 and 7, we can point out
that stronger FOEWPT demands relatively larger effective quartic couplings in
the potential and relatively low triplet-like scalars in the parameter space.
However, various Higgs searches at the LHC can constrain this parameter space.
In the following subsection, we discuss the connection between the strength of
the FOPT and the production of GW signals and the interplay between the
detection of GW and the LHC searches in this region of parameter space.

Figure 8: Scatter points in the $(\lambda_{1}+\lambda_{4})-M_{\Delta}$ plane
that exhibits FOPT along the $h_{d}-$field direction. The variation of
$\alpha$ ($log(\beta/H_{n})$) is shown in the left (right) plot via the
palette.

####  5.2.1 Production of GW and the interplay with LHC

In section 2.3.2, we have discussed the possibility of producing stochastic GW
background from the cosmological FOPT in the early Universe. This can be
detected in the various future proposed GW space/ground-based detectors. The
important portal parameters that control the GW intensity are $\alpha$,
$\beta/H_{n}$, $T_{n}$ and $v_{w}$. We set $v_{w}=1$ for this work. we have
already discussed the variation of $T_{n}$ in the the bottom, right plot of
figure 6. In this subsection, we discuss the variation of the other two main
portal parameters $\alpha$ and $\beta/H_{n}$, which control the GW spectrum.

In figure 8, we present two plots to show the variation of $\alpha$ (top,
left) and $\beta/H_{n}$ (top, right) at $T=T_{n}$ in colors in the
$(\lambda_{1}+\lambda_{4})-M_{\Delta}$ plane. $\alpha$ ($\beta/H_{n}$)
increases (decreases) with increasing ($\lambda_{1}+\lambda_{4}$) and
decreasing $M_{\Delta}$. These variations have a direct connection with
$\xi_{n}$. Stronger FOPT corresponds to larger $\alpha$ and lower
$\beta/H_{n}$ at $T=T_{n}$. From the discussion in section 2.3.2, it can be
found that the magnitude of the peak of the GW intensity is proportional to
$\alpha$ and inversely proportional to $\beta/H_{n}$ for fixed $v_{w}$ and
$T_{n}$ (see equations 43 and 47). These dependencies can be understood from
their physical definitions. A larger $\alpha$ corresponds to the more energy
transfer from the plasma to the form of GW and a smaller $\beta/H_{n}$ implies
a longer phase transition period. Thus larger $\alpha$ and small $\beta/H_{n}$
enhance the GW intensity. Therefore, one can expect an increase in the GW
intensity with $\xi_{n}$. As a result, the parameter space with larger
effective quartic coupling and smaller triplet-like scalars would produce
stronger GW spectrum intensity. In addition to this, as we discuss in section
2.3.2, it can be found that the value of $\beta/H_{n}$ is crucial to set the
peak frequency range (see equations 44 and 48) which is also important for the
detection purpose at the various future proposed GW experiments.

Figure 9: Left: Parameter points of the scan in the $log(\beta/H_{n})-\alpha$
plane with $T_{n}$ indicating via the palette. Right: Scatter plots of GW peak
amplitude ($\Omega_{\text{GW}}h^{2}$) vs the peak frequency
($f_{\text{peak}}$) in the unit of Hz considering only the contribution coming
from the sound waves mechanism. Variation of the signal strength,
$\mu_{\gamma\gamma}$, is indicated via the palette.

Therefore, in the left plot of figure 9, we present the variation of $\alpha$
and $log(\beta/H_{n})$ at $T=T_{n}$. $T_{n}$ is presented in palette color.
The color variation shows that relatively lower $T_{n}$ corresponds to higher
$\alpha$. On the other hand, $\beta/H_{n}$ measures the inverse duration of
the phase transition. In some situations, the phase transition takes longer to
start, corresponding to relatively lower $T_{n}$ and a larger gap between
$T_{c}$ and $T_{n}$. These scenarios correspond to lower $\beta/H_{n}$. On the
contrary, even at relatively low temperatures, the phase transition can happen
quickly, corresponding to higher $\beta/H_{n}$. In the plot, the purple points
(lower $T_{n}$) is across the wide range of $\beta/H_{n}$.

Variations of these portal parameters $\alpha$, $\beta/H_{n}$, $T_{n}$
correspond to a wide range of GW peak amplitude and peak frequency. We have
already discussed in section 2.3.2 that the sound wave contribution mainly
dominates the GW peak amplitude and the peak frequency (see equations 44 and
43). In the right plot of figure 9, we present the peak amplitude
($\Omega_{\text{GW}}h^{2}$) vs the peak frequency ($f_{\text{peak}}$) in the
unit of Hz considering only the contribution coming from the sound waves
mechanism. The future proposed GW experiments have different sensitive regions
in the intensity amplitude and the frequency range. Thus, studying the
variation of the peak frequency and the peak amplitude of the produced GW is
crucial for detecting the spectrum in future proposed experiments. The
projected sensitivity of those experiments is presented in figures 11 and 15.

We have already discussed that the regions of parameter space with lower
triplet-like scalar masses and correspondingly higher effective quartic
couplings can generate larger GW intensity signals. In section 4.1, we show
the variation of $\mu_{\gamma\gamma}$ with the various parameters of the model
like $\lambda_{1},\lambda_{4}$. It is expected that at relatively larger
effective quartic couplings and light charged triplet-like scalars,
$\mu_{\gamma\gamma}$ deviates significantly from 1. Therefore, we present the
variation of $\mu_{\gamma\gamma}$ in the palette-color in the bottom, right
plot of figure 9. The plot shows that the $\mu_{\gamma\gamma}$ significantly
diverges from 1 for larger GW peak amplitudes. We find that most of the region
of parameter space of FOPT, which can be probed via the GW experiments,
$\mu_{\gamma\gamma}$ is more than $3\sigma$ away from the current ATLAS and
CMS limits, whereas few points lie within the $3\sigma$ limit from the latest
observations. Therefore, most of the region of parameter space of this
scenario that lies within the sensitivity of various GW detectors are already
excluded from the latest precision study of SM-like Higgs boson, particularly
$h^{0}\rightarrow\gamma\gamma$ decays channel. The region that lies within the
$3\sigma-$limit entire region is expected to be tested from the precision
study of SM-like Higgs boson in the HL-LHC and/or future colliders. This
scenario has an interesting interplay between LHC physics and various GW
detectors regarding the possibility of FOPT in the present model. The
possibilities of detecting GW at the future proposed GW detectors due to FOPT
from the current model would be severely limited in the absence of new physics
at the HL-LHC.

In the above discussions, we talk about the phase transition that always
happens in the $h_{d}$-field direction. We have only discussed how different
model parameters from the triplet sector of this model affect the phase
transition quantities and the correlations the production of stochastic GW
signals and the LHC searches. We find that such correlations excluded most of
the regions of the triplet sector that can produce GW as a result of an
FOEWPT. On the other hand, the model parameters from the dark sector,
particularly $\lambda_{SH}$ and $m_{S}$, can play a crucial role in altering
the SM Higgs potential in favour of FOEWPT. Before looking at the influence of
the dark sector parameters on the FOPT in the next section, we first provide
two benchmark points from the above scan results to illustrate the effect of
the triplet-sector parameters on FOEWPT where the dark sector couplings are
small. The implications of the dark sector parameters on the FOPT and the
relationships between the generated GW and other DM experimental,
observational constraints are then discussed in section 5.3.

Input/Observables |  BP1 |  BP2
---|---|---
$\lambda_{1}$ | 3.92 | 3.59
$\lambda_{4}$ | $-0.989$ | $-0.56$
$M_{\Delta}$ | 367.7 | 366.0
$v_{t}$ (GeV) | $3.3\times 10^{-4}$ | $6.1\times 10^{-4}$
$\mu_{3}$ (GeV) | $-540.9$ | 441.1
$\lambda_{SH}$ | 0.0053 | 0.0013
$\lambda_{ST}$ | 0.913 | 0.02
$m_{{}_{\text{D}M}}$ | 432.2 | 420.0
$m_{H^{++}}$ (GeV) | 406.4 | 388.3
$m_{H^{+}}$ (GeV) | 387.5 | 377.3
$m_{H^{0}\sim A^{0}}$ (GeV) | 367.7 | 365.9
$m_{h}$ (GeV) | 125 | 125
$\sin\theta_{t}$ (GeV) | $-1.03\times 10^{-6}$ | $5.69\times 10^{-7}$
$\Omega h^{2}$ | 0.008 | 0.120
$\xi_{\text{DM}}\,\sigma_{\rm DD}^{\rm SI}$ (cm2) | $(9.17)9.31\times 10^{-50}$ | $8.81(8.99)\times 10^{-50}$
SNR (LISA) | 16.7 | $<<1$
Table 3: Set-I benchmark scenarios with relatively larger (smaller) portal
couplings between the triplet (singlet) and the SM Higgs doublet that exhibit
an FOEWPT. Shown are the various model input parameters, relevant masses,
mixing, DM relic density, DMDDSI cross-section and the signal-to-noise ratio
of the produced GW as a result of an FOEWPT in the LISA experiment. The
details of the EWPT of these benchmark scenarios are presented in tables 4 and
5.

####  5.2.2 Benchmark scenarios (Set-I)

In this subsection, we present two benchmark points, BP1 and BP2, in table 3
to illustrate the effect of the triplet-sector parameters on the FOEWPT. In
the triplet-sector, $v_{t}$ is around $10^{-4}$ GeV, $H^{\pm\pm}$ around 400
GeV and for BP1 (BP2) $\Delta m$ is around 20 GeV (10 GeV). In the dark
sectors, $m_{s}>m_{H^{0},A^{0},H^{\pm},H^{\pm\pm}}$ and the DM annihilation
cross-section was large at the early Universe as a result of the opening of
the various annihilation channels of the triple-sector. Since
$\lambda_{S\Delta}$ is smaller in BP2 than in BP1, the DM relic density is
lower in BP1. Lower values of $\lambda_{SH}$ in BP2 result in a substantially
low $\sigma_{\rm DD}^{\rm SI}$ compared to BP1. However, relic density scaling
sets $\sigma_{\rm DD}^{\rm SI}$ value at the same range. For both BP1 and BP2,
$\sigma_{\rm DD}^{\rm SI}$ is below the neutrino floor. Therefore, it is
difficult to probe from the DMDD experiments. However, a significant amount of
parameter space left will be probed (corresponds to relatively larger
$\lambda_{SH}$) from the future DMDD experiments.

BM No |  $T_{i}$ (GeV) | ${\\{h_{d},h_{t},h_{s}\\}}_{\text{false}}$ | $\xrightarrow[\text{type}]{\text{Transition}}$ |  ${\\{h_{d},h_{t},h_{s}\\}}_{\text{true}}$ (GeV) | $\gamma_{{}_{\rm{EW}}}$
---|---|---|---|---|---
BP1 | $T_{c}$ | 76.9 | {0, 0, 0} | FO | {230.1, 0, 0} |
$T_{n}$ | 58.1 | {0, 0, 0} | ,, | {241.3, 0, 0} | 4.15
BP2 | $T_{c}$ | 115.0 | {0, 0, 0} | ,, | {120.9, 0, 0} |
$T_{n}$ | 114.4 | {0, 0, 0} | ,, | {130.4, 0, 0} | 1.14
Table 4: Phase transition characteristics of the benchmark scenarios BP1 and
BP2 presented in table 3. Values of $T_{c}$, $T_{n}$, the corresponding field
values at the false and true phases and the strength of the phase transition,
($\gamma_{{}_{\rm{EW}}}$), along the $SU(2)$-direction are presented for each
benchmark scenario. ‘FO’ corresponds that the phase transition is first-order
type.

Figure 10: Phase flows in benchmark scenarios BP1 and BP2. Individual color
corresponds to a specific minimum of the potential (phase) while each line
indicates the evolution of the phase in the $h_{d}-$field direction as a
function of temperature. The arrows reflect the direction of transition from a
false vacuum to a true vacuum, as calculated at $T_{c}$ and $T_{n}$ along the
$h_{d}$-field. See the text for the details.

The values of $T_{c}$ and $T_{n}$, the strength of the phase transition
$\xi_{n}$ and the field values of the phases at $T_{c}$ and $T_{n}$ are
presented in Table 4. The calculation of $T_{c}$ and $T_{n}$ indicates that
these EWPTs are a one-step transition. In BP1, large supercooling is possible
since $T_{n}$ and $T_{c}$ have a significant difference. This makes the
strength of the transition relatively large. On the other hand, for BP2, the
$T_{n}$ and $T_{c}$ are very close to each other. Also, $T_{n}$ in BP1 is much
smaller than in BP2.

We present the evolution of the symmetric and broken phases with temperature
in figure 10 for these two benchmark points. Each line (red and black color)
shows the field values at a particular minimum as a function of temperature.
The black line represents the symmetric phase, whereas the red line indicates
the broken phase. In addition, the green arrow indicates that an FOPT may
happen in the direction of the arrow since at that temperature and the two
phases connected by the arrow are degenerate. The blue color arrow shows the
transition direction at the nucleation temperature along which phase
transition actually starts. In contrast to BP2, where $T_{c}$ and $T_{n}$ are
pretty near, BP1 exhibits a significant separation between them, resulting in
a stronger phase transition than BP2. The evolution of the broken phase for
both benchmark points indicates that the $h_{d}$ value (red line) finally
evolves to 246 GeV at $T=0$. Thus, the Universe finally evolves to the correct
EW minima, i.e., $v=\sqrt{v_{d}^{2}+v_{t}^{2}}=246$ GeV. Note that the
contribution of the triplet ${\it vev}$ to $v$ is minimal as $v_{t}\sim
10^{-4}$ GeV at $T=0$.

BM No. | $T_{n}$ (GeV) | $\alpha$ | $\beta/H_{n}$
---|---|---|---
BP1 | 58.1 | $0.325$ | $33.58$
BP2 | 114.4 | 0.021 | 26689.5
Table 5: Values of the parameters $T_{n}$, $\alpha$ and $\beta/H_{n}$ (that
control the GW intensity) for the benchmark points BP1 and BP2, presented in
table 3. Figure 11: GW energy density spectrum with respect to frequency for
the benchmark scenarios BP1 and BP2 illustrated against the experimental
sensitivity curves of the future proposed GW detectors such as LISA, SKA,
TianQin, Taiji, aLigo+, BBO, NANOGrav, and U-DECIGO. The solid lines indicate
the overall GW energy density produced for the specific benchmark scenario,
while the broken lines in various colours reflect the individual contributions
from sound waves and turbulence. See the text for details.

Various key parameter values ($\alpha$, $\beta$) pertaining to the GW spectra
generated from the FOPTs for these two benchmark scenarios (BP1 and BP2) are
presented in table 5. The corresponding GW (frequency) spectra are presented
in figure 11 using equations 42 to 49. For each phase transition process, the
contributions from sound waves and turbulence are shown in different colors
with broken lines. The GW peak amplitude and the peak frequency are mostly
dominated by sound contributions. These GW spectra are further compared with
the projected sensitivity of some ground- and space-based GW detectors, viz.,
LISA LISA:2017pwj , ALIA Gong:2014mca , Taiji Hu:2017mde , TianQin
TianQin:2015yph , aLigo+ Harry:2010zz , SKA Carilli:2004nx ; Janssen:2014dka ;
Weltman:2018zrl , NANOGrav McLaughlin:2013ira , Big Bang Observer (BBO)
Corbin:2005ny  and Ultimate(U)-DECIGO Kudoh:2005as  in Fig. 11. Note that, a
significant amount of supercooling happens in BP1 and the corresponding
nucleation temperature is relatively low compared with BP2. This also enhances
the $\alpha$ value for BP1 compared with BP2. Also, the nucleation takes a
much longer time for BP1 compared with BP2. The $\beta$ parameter, which
indicates the inverse of the duration of the phase transition, is smaller in
BP1 than in BP2. Thus, as the FOPT happens much stronger in BP1, the GW
spectrum peak is expected to be much larger in BP1 than in BP2. This behaviour
can be observed in Fig. 11. Note that the GW intensity obtained for BP1 lies
within the sensitivities of LISA, Taiji, ALIA, BBO and UDECIGO, while BP2 lies
only within the sensitivity of UDECIGO and marginally touches the BBO
sensitivity.

The quantity known as the signal-to-noise ratio (SNR) is used to measure the
detection of the GW signal in different experiments. It is defined as
Caprini:2015zlo

$\small{\text{SNR}=\sqrt{\delta\times\mathcal{T}\int_{f_{min}}^{f_{max}}df\bigg{[}\frac{h^{2}\Omega_{\text{GW}}(f)}{h^{2}\Omega_{\text{exp}}(f)}\bigg{]}^{2}},}$
(59)

where $\delta$ corresponds to the number of independent channels for cross-
correlated detectors (to determine the stochastic origin of GW). It is 2 for
BBO, U-DECIGO and 1 for LISA. The duration of the experimental mission in
years is defined by $\mathcal{T}$. In this work, we consider $\mathcal{T}=5$.
The effective power spectral density of the experiment’s strain noise is
indicated by $\Omega_{\text{exp}}(f)\,h^{2}$. For the detection prospects of
these experiments, the SNR values need to cross the threshold value which
depends on the configuration details of the experiment. Such as, the
recommended threshold number is around 50 for a four-link LISA setup, however,
a six-link design allows for a much lower value of around 10 Caprini:2015zlo .
In this work, we calculate the SNR value only for LISA for these two benchmark
points. For BP1, it is 16.7 and for BP2, the SNR value is way below 1. It is
expected that the SNR value of BP2 will be substantially lower because its GW
signal power spectrum does not fall within the LISA sensitivity curves.
However, other experiments like U-DECIGO can detect this type of benchmark
scenario.

Note that in both the benchmark scenarios, $\mu_{\gamma\gamma}$ is more than
the $3\sigma$ away from the latest LHC limits. Thus, as we discussed in the
description of the right plot of figure 9, FOEWPT preferred region of
parameter space of the triplet sector is mostly excluded from the precision
study of the SM Higgs boson di-photon decay modes. Since the FOEWPT is crucial
for the EWBG, we investigate its viability in the current model by changing
the parameters of the dark sector. In the next section, we will examine such
possibilities and determine the consequences of the dark sector parameters on
the FOEWPT as well as the establishing connections between the produced GW and
other DM experimental constraints.

###  5.3 FOPT from the dark sector

In the previous discussions of this article, we set the dark sector portal
couplings $\lambda_{SH}$ and $\lambda_{S\Delta}$ at smaller values and the
phase transition along the $h_{d}-$field direction is mostly controlled by the
triplet sector parameters of the model, i.e., $\lambda_{1}$, $\lambda_{4}$ and
$M_{\Delta}$. However, larger DM Higgs portal coupling, particularly
$\lambda_{SH}$ as the $h_{d}-$field couples with the $h_{s}-$field through
this portal coupling in the dark sector, can alter the Higgs potential in such
a way that FOEWPT occurs. In addition to this, another intriguing possibility
is the FOPT in the dark sector. In this case, $h_{s}-$field changes from zero
${\it vev}$ at high temperature to non-zero ${\it vev}$ at the electroweak
temperature scale, and then it evolves to zero ${\it vev}$ at $T=0$. Moreover,
a two-step FOPT in the scalar potential is even another intriguing
possibility. In this category, the $h_{s}-$field acquires non-zero ${\it vev}$
at a certain temperature in the early Universe through an FOPT and
subsequently, relatively at a lower temperature another phase transition
occurs where the $h_{d}-$field develops a non-zero ${\it vev}$ whereas the
${\it vev}$ of the $h_{s}-$field goes from a non-zero value to a zero value.
It is crucial to remember that among these three forms of phase transitions,
the one-step phase transition along the $h_{s}-$direction is capable of
generating GW signal but is ineffective for the electroweak baryogenesis.

To study the phase transition along the $h_{s}-$field direction, we first fix
the triplet sector parameters the same as the benchmark scenario BP1 and vary
the dark sector parameters $m_{S}$, $\lambda_{SH}$, $\lambda_{S\Delta}$,
$\lambda_{S}$ and $\mu_{3}$. The DM phenomenology due to the variation of
these DM input parameters for the same triplet-sector input parameters has
already been discussed in section 5.1. In figure 12, we present the variation
of $m_{s}$ and $\lambda_{SH}$ ($\lambda_{S\Delta}$) in the left (right) plot
that exhibits FOPT along the $h_{s}-$field direction. The strength of the
phase transition along the $h_{s}-$field direction, i.e.,
$\xi_{s}=\frac{<h_{s}>}{T_{n}}$ where $<h_{s}>$ is the ${\it vev}$ of the
$h_{s}$ field at the broken phase at $T=T_{n}$, is presented via the palette
color.

In the left plot, the variation of $m_{s}$ and $\lambda_{SH}$ reveals that for
the phase transition along the $h_{s}-$field direction, $\lambda_{SH}$ must be
large for larger $m_{s}$ and lower for smaller $m_{s}$. The region of
parameter space with $m_{s}\gtrsim 425$ GeV is almost disfavoured for FOPT for
$\lambda_{SH}\lesssim 6.0$. On the other hand, the right plot indicates that
the phase transition does not have a preferential bias over the choice of
$\lambda_{S\Delta}$. Equation 19 elucidates such connection between $m_{s}$
and $\lambda_{SH}$. The phase transition along the singlet direction generally
prefers relatively low bare mass squared term ($\mu_{s}^{2}$) in the potential
Ghorbani:2018yfr . Therefore, the parameter space that corresponds to an FOPT
along the $h_{s}$ direction, the maximum possible value of $m_{s}$ is mainly
controlled by the term $\lambda_{SH}v_{d}^{2}$ and it is almost independent of
$\lambda_{S\Delta}$ as $v_{t}<<1$. The color variation in the left plot
reveals that the strength of the phase transition is maximum at the outer
edge, corresponding to the parameter space consisting of very small bare mass
$\mu_{s}$.

The variations of these parameters will also have effects on the observables
of the dark sector, such as the $\Omega_{S}h^{2}$ and $\sigma_{\rm DD}^{\rm
SI}$, and make these scenarios more sensitive to the various DM experiments.
It is therefore anticipated that the production of GW as a result of this type
of FOPT in the dark sector will be connected with DM observables. In the
following subsection, we discuss the relationship between the DM phenomenology
and FOPT along the $h_{s}-$field direction. Note that, in the previous
discussion, we found a similar type of correlation between the FOEWPT-induced
GW intensity and the precision measurement of the signal strength,
$\mu_{\gamma\gamma}$ of the $h^{0}\rightarrow\gamma\gamma$ channel at the LHC.

Figure 12: Scatter plots in the $m_{S}-\lambda_{SH}$
($m_{S}-\lambda_{S\Delta}$) plane in the left (right) plot that exhibits FOPT
along the $h_{s}-$field direction. The variation of the strength of the phase
transition $\xi_{S}$ is shown via the palette.

####  5.3.1 Correlations among the DM observables and an FOPT

In this section, we now turn our attention to the link between the DM
observables and FOPT along the $h_{s}-$field direction. There is an effect on
the various DM observables from the parameter space that exhibits such a phase
transition.

The left plot figure 13 is otherwise similar with the left plot of figure 12
but here the points satisfy the Planck experimental constraints on the relic
density, i.e., $\Omega_{s}h^{2}\leq 0.120$. The relic density of the singlet
scalar field is indicated by the palette color in the plot. Note that, as we
described earlier $m_{s}<h^{0}$, $\lambda_{SH}$ mostly controls the
$\Omega_{s}h^{2}$ whereas for $m_{s}>h^{0}$, the $\mu_{3}$ proportional semi-
annihilation process starts to contribute to DM density. $\lambda_{S\Delta}$
starts to contribute significantly when $m_{s}$ is larger than the triplet-
like scalar masses which is in this case around $\sim$ 400 GeV.

The $CP$-even scalars $h^{0}$ and $H^{0}$ contribute $\sigma_{\rm DD}^{\rm
SI}$. It has been discussed in section 2.2 that for smaller values of $v_{t}$
and $\sin\theta_{t}$, $\sigma_{\rm DD}^{\rm SI}$ mostly depends on the SM-like
Higgs boson medicated process and it increases with $\lambda_{SH}$. Such
enhancement of $\sigma_{\rm DD}^{\rm SI}$ with $\lambda_{SH}$ excludes a
significant amount of parameter space from the latest constraints from
XENON-1T and PANDA-4T. On the other hand, if the DM mass is larger than the
triple-like scalars then larger $\lambda_{S\Delta}$ decreases the relic
density without affecting the $\sigma_{\rm DD}^{\rm SI}$. With this, at a
relatively larger $\lambda_{S\Delta}$ value, the DM could be underabundant.
The small scaling factor due to the low relic density actually reduces the
effective $\sigma_{\rm DD}^{\rm SI}$ value. Such scenarios allow the parameter
space with larger $\lambda_{SH}$ which is otherwise excluded from the latest
constraints.

Figure 13: Left: Scatter points that exhibits FOPT along the $h_{s}-$field
direction and $\Omega_{S}h^{2}<0.120$ in the $m_{S}-\lambda_{SH}$ plane.
Right: Similar to the scatter spots on the left figure, but comply with the
latest DMDDSI experimental limits. The variation of $\Omega_{S}h^{2}$
($\sigma_{\rm DD}^{\rm SI}$) is shown via the palette in the left (right)
plot.

In the right plot of figure 13, we present the points of the left plot that
pass the latest $\sigma_{\rm DD}^{\rm SI}$ constraint from the XENON-1T and
PANDA-4T experiment. The effective $\sigma_{\rm DD}^{\rm SI}$ is presented by
the palette color. Note that due to the requirement of larger $\lambda_{SH}$
for the FOPT a significant amount of the parameter space is excluded from the
$\sigma_{\rm DD}^{\rm SI}$ constraints. Note that for $m_{s}>h^{0}$, the
$\mu_{3}$ dependent semi-annihilation process starts to contribute and because
of that the relic density drops significantly which helps the effective
$\sigma_{\rm DD}^{\rm SI}$ rates to remain within the exclusion limit result
of a significant downward scaling. Because of that, we find some parameter
space for $m_{s}$ after the SM-like Higgs Boson threshold value. Therefore, it
is important to note that the term $\mu_{3}$ has an important role in the dark
sector and the FOPT sector. As it generates the semi-annihilation processes
for the DM annihilation, also it can alter the potential shape and generate a
barrier between the minima which is necessary for an FOPT. In the scenario for
$m_{s}$ larger than the triplet-like scalars, the term proportional to
$\lambda_{S\Delta}$ starts to contribute to the DM annihilation, which can
drastically reduce the relic density. However, in this scenario, as the
triplet-like scalar masses are relatively large $\sim 400$ GeV, this effect
does not start to occur significantly in the FOPT preferred allowed parameter
space as $m_{s}$ is bounded up to $\sim 425$ GeV for $\lambda_{SH}\lesssim
6.0$.

The discussion of this subsection points out the strong correlations between
the FOPT-favoured parameter space and the DM observables. Such correlations
rule out a significant amount of parameter space of having an FOPT along the
$h_{s}-$field direction. In this case, as we consider the triplet-like scalars
a bit heavier, approximately $\sim$400 GeV, the contribution from the
$\lambda_{S\Delta}h_{d}^{2}h_{s}^{2}$ interaction term in the potential on the
relic density estimation significantly start to contribute for comparatively
higher $m_{s}$. In the scenario with relatively smaller triplet-like scalars,
the same interaction becomes significant at relatively small $m_{s}$ that can
open up much more allowed parameter points satisfying the latest DMDDSI cross-
section limits from the XENON-1T and PANDA-4T experiments due to downward
scaling of the underabundant DM. In the next subsection, we choose a benchmark
scenario from that category where the triplet-like scalars are relatively low,
around $\sim$200 GeV, and exhibit FOEWPT in the early Universe and discuss the
such correlation between the FOPT, production of GW and the DM observables is
also discussed.

####  5.3.2 Benchmark scenario (Set-II)

Input/Observables |  BP3
---|---
$\lambda_{1}$ | 0.17
$\lambda_{4}$ | $-0.30$
$M_{\Delta}$ | 193.2
$v_{t}$ (GeV) | $4.5\times 10^{-4}$
$\mu_{3}$ (GeV) | -100.8
$\lambda_{SH}$ | 1.2
$\lambda_{ST}$ | 6.9
$\lambda_{ST}$ | 2.17
$m_{{}_{\text{D}M}}$ | 252.2
$m_{H^{++}}$ (GeV) | 215.3
$m_{H^{+}}$ (GeV) | 204.5
$m_{H^{0}\sim A^{0}}$ (GeV) | 193.2
$m_{h}$ (GeV) | 125
$\sin\theta_{t}$ (GeV) | $8.0\times 10^{-6}$
$\Omega h^{2}$ | $4.3\times 10^{-5}$
$\xi_{\text{{DM}}}\,\sigma_{\rm DD}^{\rm SI}$ (cm2) | $7.1\times 10^{-47}$
SNR (LISA) | $<<1$
Table 6: Set-II benchmark scenario with relatively smaller (larger) portal
couplings between the triplet (singlet) and the SM Higgs doublet that exhibits
a two-step EWPT in the early Universe. Similar to table 3, various model input
parameters, masses, mixing, DM relic density, DMDDSI cross-section and the
signal-to-noise ratio of the produced stochastic GW signal in the LISA
experiment are shown. The details of the phase transition are shown in tables
7 and 8.

In this section, we present a benchmark point (BP3) exhibiting a two-step FOPT
in the early Universe. We consider a relatively light triplet-like scalar of
around 200 GeV and much smaller $\lambda_{1}$ and $|\lambda_{4}|$ values such
that $\mu_{\gamma\gamma}$ is within the 1$\sigma$ (2$\sigma$)-limit of ATLAS
(CMS) latest results. Scenarios like this demand larger quartic couplings
(particularly $\lambda_{SH}$) in the dark sector and smaller DM mass for an
FOPT in the early Universe. In the discussion in section 5.2, we consider very
low quartic couplings of the dark sector, i.e., much smaller $\lambda_{SH}$
and $\lambda_{S\Delta}$. This can be seen in BP1 and BP2 in table 6. Because
of such consideration, we do not find points with smaller quartic coupling
points in the triplet-sector, i.e., parameter points with smaller
$\lambda_{1}$ and $|\lambda_{4}|$, in section 5.2 that can generate an FOEWPT
in the early Universe.

Compared to the scenarios discussed in section 5.2, in the BP3 $\lambda_{SH}$
is relatively large$\sim 1.2$, which facilitates successful FOPT in the scalar
potential. However, $\sigma_{\rm DD}^{\rm SI}$ is expected to be much larger
as a result of a large coupling. To comply with the latest experimental
constraints from XENON-1T and PANDA-4T we consider $m_{s}$ larger than the
triplet-like scalar masses and much larger $\lambda_{S\Delta}$ $\sim 6.9$,
such that the $S~{}S^{*}~{}\to{\rm X~{}~{}Y};~{}~{}~{}\\{\rm
X,Y\\}=\\{H,A^{0},H^{\pm},H^{\pm\pm}\\}$ annihilation processes open up and
contribute significantly to the annihilation cross-section so that DM becomes
underabundant. Such a lower relic density corresponds to a much smaller
scaling fraction $\xi_{\text{DM}}$ which helps the computed $\sigma_{\rm
DD}^{\rm SI}$ rate to comply with the latest experimental stringent limits.

Note that the value of $m_{S}$ is also crucial for the FOPT and cannot be too
large. Due to the comparatively smaller triplet-like scalar masses $\sim$ 200
GeV in this benchmark point, we can keep $m_{S}$ within the range of $\sim
250$ GeV and also take the benefit of the downward relic density scaling
factor on the DMDDSI experimental constraints. This aids in keeping the DMDDSI
cross-section of BP3 below the most recent constraints from the XENON-1T and
PANDA-4T experiments.

BM No |  $T_{i}$ (GeV) | ${\\{h_{d},h_{t},h_{s}\\}}_{\text{false}}$ | $\xrightarrow[\text{type}]{\text{Transition}}$ |  ${\\{h_{d},h_{t},h_{s}\\}}_{\text{true}}$ (GeV) | $\gamma_{{}_{\rm{EW}}}$
---|---|---|---|---|---
BP3 | $T_{c}$ | 163.2 | {0, 0, 0} | FO | {0, 0, 14.8} |
145.7 | {0, 0, 63.6} | ,, | {78.9, 0, 0} |
$T_{n}$ | 163.1 | {0, 0, 0} | ,, | {0, 0, 15.1} |
116.3 | {0, 0, 86.7} | ,, | {197.3, 0, 0} | 1.7
Table 7: Phase transition characteristics of the benchmark point BP3 presented
in table 6. Values of $T_{c}$, $T_{n}$, the corresponding field values at the
false and true phases and the strength of the phase transition,
($\gamma_{{}_{\rm{EW}}}$), along the $SU(2)$-direction are presented. ‘FO’
denotes to FOPT.

Figure 14: Phase flows in benchmark point BP3. Individual color corresponds to
a specific minimum of the potential (phase) while the each line indicates the
evolution of the phase in the $h_{d}-$field direction in the left and
$h_{s}-$field direction in the right as a function of temperature. The arrows
reflect the direction of transition from a false vacuum to a true vacuum, as
calculated at $T_{c}$ and $T_{n}$ along the $h_{d}$-field, whereas an orange
bullet indicates that the corresponding field value does not change much
throughout that transition. See the text for the details.

In contrast with BP1 and BP2, the phase transition pattern is quite different
in BP3. The values of $T_{c}$ and $T_{n}$ for the BP3 benchmark scenario are
presented in table 7 and it suggests that the transition is a two-step FOPT.
The evolution of the field and the various phases with temperature are
presented in figure 14. Note that three types of phase exist. ‘phase-I’
denotes the trivial phase, i.e., all the field values are always zero. On the
other hand, ‘phase-II’ corresponds to the minimum where only the $h_{s}-$
field acquires a non-zero ${\it vev}$ and other field ${\it vevs}$ remain
zero. ‘phase-III’ corresponds to that minimum of the potential where only
$h_{d}-$ field obtains a non-zero ${\it vev}$ and all other field ${\it vevs}$
remain at zero. In the right and left plots of figure 14, the trivial minimum
$\\{0,0,0\\}$ at the high temperature is denoted by the black line (phase-I)
where the broken phases in the direction of $h_{s}$ and $h_{d}$ are presented
in the red color (phase-II) and magenta color (phase-III). These lines
correspond to the field values of a particular minimum as a function of
temperature. The arrows in green color correspond to the temperature and field
values of two phases where they degenerate and a possible FOPT can happen
along the direction of the arrow. The phase at the end of the arrow indicates
that it has a deeper minimum than the other phase for $T<T_{c}$. The
transitions finally start to occur when a bubble of the true minimum finally