# Multicomponent Scalar Dark Matter with an Extended Gauge Sector

Baradhwaj Coleppa<EMAIL_ADDRESS>Indian Institute of Technology
Gandhinagar, Gandhinagar 382 055, India  Kousik Loho<EMAIL_ADDRESS>Indian Institute of Technology Gandhinagar, Gandhinagar 382 055, India  Agnivo
Sarkar<EMAIL_ADDRESS>Regional Centre for Accelerator-based Particle
Physics, Harish-Chandra Research Institute, HBNI, Chhatnag Road, Jhusi,
Prayagraj (Allahabad) 211019, India

###### Abstract

We consider an extension of the Standard Model of particle physics with an
additional $SU(2)$ gauge sector along with an additional scalar bidoublet and
a non-linear sigma field. The neutral components of the bidoublet serve as
dark matter candidates by virtue of the bidoublet being odd under a $Z_{2}$
symmetry. Generic beyond Standard Model constraints like vacuum stability,
invisible decay of higgs, Higgs alignment limit and collider constraints on
heavy gauge bosons restrict the parameter space of this model. In this
multicomponent dark matter scenario, we investigate the interplay between the
annihilation and co-annihilation channels originating from the new gauge
sector as those contribute to the relic abundance. We also inspect the direct
detection constraints on scattering cross-sections of the dark matter
particles with the detector nucleons and present our observations.

HRI-RECAPP-2023-06

##  I Introduction

The Standard Model (SM) of particle physics has been experimentally proven to
be a very successful theory over the last few decades. However, there are a
few limitations of the SM which motivate physicists to look for Beyond the
Standard Model (BSM) scenarios. One such very important drawback with the SM
is that it can not explain the particulate nature of the Dark Matter (DM). DM
consists of more than one fourth of the energy budget of our universe Aghanim
et al. (2020a); Hinshaw et al. (2013). In the quest of solving the DM problem,
one very popular class of models are those that incorporate a Weakly
Interacting Massive Particle (WIMP) in their spectrum Bertone et al. (2005);
Bergstrom (2009); Arcadi et al. (2018). In a typical WIMP model, there is an
extended sector in addition to the SM particles which is weakly coupled to the
SM and protected by a stabilising symmetry. Various scalar extensions
including singlet extensions McDonald (1994); Burgess et al. (2001); Guo and
Wu (2010); Bandyopadhyay et al. (2010); He et al. (2009, 2010); Cline et al.
(2013), two-singlet extensions Abada et al. (2011); Abada and Nasri (2012);
Arhrib and Maniatis (2019); Hamada et al. (2021); Modak et al. (2015);
Maniatis (2021); Bhattacharya et al. (2017a); Maity and Ray (2020); Basak et
al. (2021), inert 2HDM Lopez Honorez et al. (2007); Belyaev et al. (2018),
2HDM with gauge extension Rojas-Abatte et al. (2017) and also various gauge
extensions mostly of different kinds of additional $U(1)$ Bhat and Adhikari
(2020); Nam (2020); Basak and Mondal (2014) with an extended scalar and/or
fermion sector have been studied in this context.

In this work, we extend the gauge sector of SM with an additional $SU(2)$ and
also the scalar sector with an additional scalar bidoublet and a non-linear
sigma field. The SM-like scalar doublet and the sigma field contribute to the
electroweak symmetry breaking whereas the scalar bidoublet provides the DM
candidates protected by an imposed $Z_{2}$ symmetry. An extended gauge sector
introduces co-annihilation channels for the relic calculation as well as
modifies the direct detection constraints - the main goal of the present work
is to study in this specific model context how these phenomenological
considerations affect the parameter space of the model.

The paper is organised as follows: in Sec. II we introduce the model. Then in
Sec. III we discuss the model constraints arising from vacuum stability,
invisible decay of the Higgs particle, Higgs alignment and collider searches
for a heavy gauge boson. In Sec. IV, the DM phenomenology has been explored
while motivating the DM candidates, calculating the relic abundance and the
direct detection constraints before finally presenting our conclusions in Sec.
V. All the relevant couplings are organised in Appendix B after a few comments
on perturbativity in Appendix A. A brief note on the loop level calculations
of Higgs decay to two photons is provided in Appendix C with appropriate loops
relevant to the model under study. The Boltzmann equations and Feynman
diagrams relevant for the relic, are given in Appendix D.

##  II The Model

We consider an $SU(2)$ extension of the gauge sector of the SM. Furthermore,
the scalar sector is extended with an additional scalar bi-doublet odd under a
$Z_{2}$ symmetry and a non-linear sigma field. The gauge sector can now be
written as $SU(2)_{0}\times SU(2)_{1}\times U(1)_{2}$ (where the subscripts
help us to keep track of different gauge charges). We will see that three out
of the four degrees of freedom of the scalar doublet and three degrees of
freedom of the sigma field are eaten up by the six gauge bosons of the gauge
sector, thus giving rise to massive gauge bosons, while the photon remains
massless due to the residual $U(1)_{EM}$ symmetry. The electrically neutral
components of the scalar bidoublet protected by a $Z_{2}$ symmetry furnish the
dark matter candidates of the universe. The scalar sector can be conveniently
expressed in the form

$\Phi_{1}=\begin{pmatrix}G^{+}\\\
\frac{v_{1}+h+iG_{0}}{\sqrt{2}}\end{pmatrix};\qquad\Phi_{2}=\frac{1}{2}\begin{pmatrix}\phi_{0}+i\phi_{3}&i(\phi_{1}-i\phi_{2})\\\
i(\phi_{1}+i\phi_{2})&\phi_{0}-i\phi_{3}\end{pmatrix};\qquad\Sigma=\exp{\frac{i\sigma^{a}\pi^{a}}{F}}.$
(1)

After a few redefinitions, the latter two can be recast as

$\Phi_{2}=\frac{1}{2}\begin{pmatrix}\phi_{0}+i\phi_{3}&i\sqrt{2}\phi^{+}\\\
i\sqrt{2}\phi^{-}&\phi_{0}-i\phi_{3}\end{pmatrix};\qquad\Sigma=\begin{pmatrix}1+\frac{i}{F}\pi^{3}&\frac{i\sqrt{2}}{F}\pi^{+}\\\
\frac{i\sqrt{2}}{F}\pi^{-}&1-\frac{i}{F}\pi^{3}\end{pmatrix},$ (2)

where $\phi^{\pm}=\frac{1}{\sqrt{2}}(\phi_{1}\mp i\phi_{2})$ and similarly for
the $\pi$’s.
The scalars transform under these gauge symmetries in the following manner:

$\Phi_{2}\rightarrow U_{0}^{\dagger}\Phi_{2}U_{1},\qquad\Sigma\rightarrow
U_{1}^{\dagger}\Sigma U_{0},\qquad\Phi_{1}\rightarrow
U_{1}\Phi_{1},U_{2}\Phi_{1}.$ (3)

Thus the $\Phi_{2}$ and $\Sigma$ transform as bidoublets under the
$SU(2)_{0}\times SU(2)_{1}$whereas the $\Phi_{1}$ is a doublet under
$SU(2)_{1}$. The $U_{i}$ signify the unitary transformation matrices
corresponding to the gauge symmetry indexed $i$ and the covariant derivatives
take the form

$\displaystyle D_{\mu}\Phi_{1}$
$\displaystyle=\partial_{\mu}\Phi_{1}-\frac{ig_{1}}{2}W^{a}_{1\mu}\sigma^{a}\Phi_{1}-\frac{ig_{2}}{2}B_{2\mu}\Phi_{1}$
$\displaystyle D_{\mu}\Phi_{2}$
$\displaystyle=\partial_{\mu}\Phi_{2}+\frac{ig_{0}}{2}W^{a}_{0\mu}\sigma^{a}\Phi_{2}-\frac{ig_{1}}{2}\Phi_{2}W^{a}_{1\mu}\sigma^{a}$
$\displaystyle D_{\mu}\Sigma$
$\displaystyle=\partial_{\mu}\Sigma+\frac{ig_{1}}{2}W^{a}_{1\mu}\sigma^{a}\Sigma-\frac{ig_{0}}{2}\Sigma
W^{a}_{0\mu}\sigma^{a}.$

With these definitions in place, the kinetic energy term of the scalar sector
can be written as

$\mathcal{L}\supset(D^{\mu}\Phi_{1})^{\dagger}(D_{\mu}\Phi_{1})+\textrm{Tr}\,[(D^{\mu}\Phi_{2})^{\dagger}(D_{\mu}\Phi_{2})]+\frac{F^{2}}{4}\textrm{Tr}[(D^{\mu}\Sigma)^{\dagger}(D_{\mu}\Sigma)].$
(4)

The $SU(2)_{0}\times SU(2)_{1}\times U(1)_{2}$ symmetry is broken down to
$SU(2)_{L}\times U(1)_{Y}$ as $\Sigma$ develops a vacuum expectation value $F$
and then this symmetry is further broken down to $U(1)_{EM}$ as the $\Phi_{1}$
develops a vacuum expectation value $v$ (= 246 GeV). After symmetry breaking,
collecting the relevant terms, the neutral gauge boson mass matrix can be
written as

$\text{M}^{2}_{N}=\frac{1}{8}\begin{pmatrix}g_{0}^{2}F^{2}&-g_{0}g_{1}F^{2}&0\\\
-g_{0}g_{1}F^{2}&g_{1}^{2}(F^{2}+v^{2})&-g_{1}g_{2}v^{2}\\\
0&-g_{1}g_{2}v^{2}&g_{2}^{2}v^{2}\end{pmatrix}=\frac{g_{0}^{2}F^{2}}{8}\begin{pmatrix}1&-x&0\\\
-x&x^{2}(1+r^{2})&-x^{2}tr^{2}\\\ 0&-x^{2}tr^{2}&x^{2}t^{2}r^{2}\end{pmatrix}$
(5)

where $\frac{g_{1}}{g_{0}}=x\ll 1$, $\frac{g_{2}}{g_{1}}=t$, $\frac{v}{F}=r$.
This matrix can be diagonlized perturbatively in the small parameter $x$. The
eigenvalues corresponding to the eigenstates of the photon, $Z$ and
$Z^{\prime}$ bosons are given by

$m_{\gamma}^{2}=0;\qquad
m_{Z}^{2}=\frac{g_{0}^{2}F^{2}r^{2}x^{2}}{4}(1+t^{2}-x^{2});\qquad
m_{Z^{\prime}}^{2}=\frac{g_{0}^{2}F^{2}}{4}(1+x^{2}),$ (6)

and the mass eigenstates are (writing $t=\frac{s}{c}$)

$A_{\mu}=sxW_{0\mu}^{3}+sW^{3}_{1\mu}+cB_{2\mu}$ (7)
$Z_{\mu}=cxW_{0\mu}^{3}+cW_{1\mu}^{3}-sB_{2\mu}$ (8)
$Z_{\mu}^{\prime}=W_{0\mu}^{3}-xW_{1\mu}^{3}.$ (9)

The charged gauge boson mass matrix can similarly be computed to be

$M_{C}^{2}=\frac{g_{0}^{2}F^{2}}{4}.\begin{pmatrix}1&-x\\\
-x&x^{2}(1+r^{2})\end{pmatrix},$ (10)

and the mass eigen values are

$m_{W}^{2}=\frac{g_{0}^{2}F^{2}r^{2}x^{2}}{4}(1-x^{2});\qquad
m_{W^{\prime}}^{2}=\frac{g_{0}^{2}F^{2}}{4}(1+x^{2}).$ (11)

The mass eigenstates of the charged gauge bosons are

$W_{\mu}=xW_{0\mu}+W_{1\mu}$ (12) $W_{\mu}^{\prime}=W_{0\mu}-xW_{1\mu}$ (13)

where the one with the smaller mass is identified as the $W$ boson and the
heavier one as the $W^{\prime}$ boson. It is to be noticed that the massive SM
gauge bosons get their contribution mostly from site $SU(2)_{1}$. However, for
the BSM gauge bosons, $SU(2)_{0}$ is the major contributor. In the limit
$g_{0}\rightarrow\infty$ (or equivalently $x\rightarrow 0$), the
$m_{W^{\prime},Z^{\prime}}$ decouple and the gauge spectrum at low energies
thus reduces to that of the SM as expected. For all practical purposes $x$ can
be parametrized as the ratio of the light and heavy charged gauge boson masses
(to leading order) as

$x^{2}\approx\frac{m_{W}^{2}}{r^{2}m^{2}_{W^{\prime}}}$ (14)

which can be readily seen from Eqns. 11 and 6.

After taking into account the gauge charges of the scalars and the imposed
$Z_{2}$ symmetry, the most general scalar potential can be written as

$\displaystyle\begin{split}V(\Phi_{1},\Phi_{2},\Sigma)&=m_{2}^{2}\,\textrm{Tr}(\Phi_{2}^{\dagger}\Phi_{2})+\lambda_{1}\bigg{[}\Phi_{1}^{\dagger}\Phi_{1}-\frac{v^{2}}{2}\bigg{]}^{2}+\lambda_{2}\bigg{[}\textrm{Tr}(\Phi_{2}^{\dagger}\Phi_{2})\bigg{]}^{2}\\\
&+\lambda_{12}\bigg{[}\Phi_{1}^{\dagger}\Phi_{1}-\frac{v^{2}}{2}\bigg{]}\bigg{[}\textrm{Tr}(\Phi_{2}^{\dagger}\Phi_{2})\bigg{]}+\lambda_{23}F^{2}\bigg{[}\textrm{Tr}(\Sigma\Phi_{2}\Sigma\Phi_{2})\bigg{]}+\tilde{\lambda}_{23}F^{2}\bigg{[}\textrm{Tr}(\Sigma\Phi_{2})\bigg{]}^{2}.\end{split}$
(15)

After spontaneous symmetry breaking the higgs mass ($m_{h}^{2}$) is given by
the expression $2\lambda_{1}v^{2}$ and the charged scalar mass
($m_{\phi^{\pm}}^{2}$) is given by $m_{2}^{2}-\lambda_{23}F^{2}$. The two
neutral components of $\Phi_{2}$ have masses given by the following
expressions:

$m_{\phi_{0}}^{2}=m_{2}^{2}+\lambda_{23}F^{2}+2\tilde{\lambda}_{23}F^{2},~{}~{}~{}~{}~{}m_{\phi_{3}}^{2}=m_{2}^{2}-\lambda_{23}F^{2}.$
(16)

It is imperative to note here that the $Z_{2}$ symmetry of the $\Phi_{2}$ bi-
doublet works as the stabilising symmetry required for the DM under the WIMP
scenario. However, this discrete symmetry can not forbid the number changing
processes among the components of this bidoublet. Thus for a positive
$(\lambda_{23}+\tilde{\lambda}_{23})$, the $\phi_{3}$ component will the
natural DM candidate and for negative values of the coupling parameters the
$\phi_{0}$ component serves the same purpose from the stability argument.
However, if the couplings are such that these two are degenerate, then we have
an interesting possibility of a multicomponent DM scenario where both neutral
components are cosmologically stable DM candidates - we describe how these
will play out in more detail in Sec. IV.

The Yukawa sector is assumed to mimic that of the SM as $\Phi_{1}$ plays a
similar role of the SM Higgs doublet. However, the fermions can still couple
to the heavy gauge bosons due to gauge mixing. All the relevant couplings can
be straightforwardly computed (again in a perturbative expansion in $x$) and
are listed in Appendix B. Specifically, the purely scalar sector couplings are
given in Tables 2 and 3. The scalar couplings with the gauge bosons are listed
in Tables 4,5,6,7,8,9 and 10 and the relevant couplings of the gauge bosons
with the fermions can be found in Tables 11 and 12 and that of the gauge
bosons among themselves in Table 13.

##  III Model Constraints

In this section, we impose all the theoretical and experimental constraints on
the model. We will not only identify all the relevant parameters for carrying
out further analysis, but also assign values to some and constrain the rest
with various experimental observations. Among the many model parameters, the
following are the most relevant for our analysis: $v$, $m_{h}$,
$m_{\phi_{3}}$, $\lambda_{12}$, $\lambda_{23}$, $\tilde{\lambda}_{23}$,
$m_{W^{\prime}}$ (or $m_{Z^{\prime}}$) and $F$ (or equivalently $r$). The
other important parameters can be represented in terms of these above
mentioned ones. For example, $x$ can be represented in terms of
$m_{Z^{\prime}}$ and $F$ following 6 or more conveniently through mass ratio
as in 11 (for a given $v$, $F$ and $r$ can be used interchangeably). We also
set $m_{h}=125$ GeV and $v=246$ GeV in this paper. We now look at four key
constraints on the model parameters coming from vacuum stability, the
branching ratio of invisible decay of the Higgs, the Higgs alignment limit and
collider constraints on the heavy gauge bosons.

###  III.1 Vacuum Stability

It is the quartic part of the scalar potential that contributes to the vacuum
stability conditions Bhattacharyya and Das (2016). After spontaneous symmetry
breaking when $\Sigma$ develops a vacuum expectation value, the quartic part
of the scalar potential in our model can be written in the unitary gauge
$\langle\Sigma\rangle=1$ as,

$V_{4}=\lambda_{1}(\Phi_{1}^{\dagger}\Phi_{1})^{2}+\lambda_{2}[Tr(\Phi_{2}^{\dagger}\Phi_{2})]^{2}+\lambda_{12}(\Phi_{1}^{\dagger}\Phi_{1})[Tr(\Phi_{2}^{\dagger}\Phi_{2})]$
(17)

The vacuum stability conditions can now be found through the standard
procedure described in Bhattacharyya and Das (2016). Let’s assume that
$\Phi_{1}^{\dagger}\Phi_{1}=a$ & $Tr(\Phi_{2}^{\dagger}\Phi_{2})=b$. Then,

$V_{4}=\lambda_{1}a^{2}+\lambda_{2}b^{2}+\lambda_{12}ab=(\sqrt{\lambda_{1}}a-\sqrt{\lambda_{2}}b)^{2}+(2\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{12})ab.$
(18)

  * •

Let us consider first the $b=0$ direction ($a\rightarrow\infty$). We demand
$V_{4}\geq 0$, which implies $\lambda_{1}a^{2}\geq 0$. The condition is then
$\lambda_{1}\geq 0$.

  * •

Next we consider the $a=0$ direction ($b\rightarrow\infty$). Similarly by
demanding $V_{4}$ to be positive, one can derive the condition
$\lambda_{2}\geq 0$.

  * •

Lastly, the direction $a=\sqrt{\frac{\lambda_{2}}{\lambda_{1}}}b$ (where
$a,b\rightarrow\infty$) can be considered. Again demanding $V_{4}\geq 0$ one
can find the condition $\lambda_{12}\geq-2\sqrt{\lambda_{1}\lambda_{2}}$.

Now collecting all the equations together the vacuum stability conditions for
this model are given by

$\lambda_{1}\geq 0,\qquad\lambda_{2}\geq
0\qquad,\textrm{and}\qquad\lambda_{12}\geq-2\sqrt{\lambda_{1}\lambda_{2}}.$
(19)

Another relevant theoretical constraint, the perturbativity bound, is
discussed in Appendix A.

###  III.2 Higgs Invisible Decay

The Higgs invisible decay branching ratio is given by

$BR_{h\rightarrow inv}=\frac{\Gamma_{h\rightarrow inv}}{\Gamma_{h\rightarrow
inv}+\Gamma_{h\rightarrow SM}}.$ (20)

Among ATLAS ATL (2020) and CMS Sirunyan et al. (2019) measurements of this
quantity, the more stringent bound comes from the ATLAS experiment and is 11%.
In Fig. 1, assuming $\phi_{3}$ to be the dark matter candidate, we present
this bound in the parameter space $m_{\phi_{3}}-\lambda_{12}$ for dark matter
masses smaller than half of $m_{h}$ (in order, of course, that $h$ can decay
into a pair of them). There is an upper bound on $\lambda_{12}$ as a function
of $m_{\phi_{3}}$ as indicated by the shaded region in the figure. For higher
DM masses the restriction will be drastically lenient due to off-shell
suppression of the decay width.

Figure 1: The shaded region in the $\lambda_{12}-m_{\phi_{3}}$ space is
disallowed by constraints from Higgs invisible branching ratio.

###  III.3 Higgs Alignment Limit

Given the extended gauge and scalar sector in our model, the standard Higgs
boson couplings to various constituents of the SM are modified due to gauge
mixing. For example, the coupling between Higgs and the pair of $W$ boson
(Table 4 ) will take the form

$\lambda_{hW^{+}_{\mu}W^{-}_{\mu}}=\frac{e^{2}v}{2\sin^{2}\theta_{w}}\left(1+x^{2}\sin^{2}\theta_{w}\right)+h.c.$
(21)

One can notice that the deviation of this coupling from the SM counterpart is
proportional to $x^{2}$. A similar deviation can be found in other SM-like
gauge couplings in Tables 4 and 5. On the other hand the couplings between the
SM Higgs and different fermion pairs in this model maintain the same structure
as SM due to the Yukawa sector of both these models being the same. To analyse
how new physics can modify Higgs properties we look into the branching ratio
of $h\to\gamma\gamma$. Apart from the usual SM fermion and $W$ boson loops,
the heavy $W^{{}^{\prime}}$ boson as well as the charged scalar $\phi^{\pm}$
will contribute to this loop mediated process. After incorporating these
additional contributions we have calculated the branching ratio
$\text{BR}\left(h\to\gamma\gamma\right)$ (for a detailed calculation, see
Appendix C) in this model for a range of charged Higgs and $W^{\prime}$
masses. In Fig. 2, we present the allowed region (shaded region is disallowed)
of the parameter space after recasting ATLAS and CMS measurements - we have
set the ratio between two _vev’s_ $r$ at 0.25 and presented our results for
two different $\lambda_{12}$ benchmark values given by 0.1 and 0.01. It is
found that for both the scenarios, the region $m_{W^{\prime}}>$ 1 TeV and
$m_{\phi^{\pm}}>$ 200 GeV is allowed.

Figure 2: The dark cyan shaded region in the $m_{\phi^{\pm}}$ vs.
$m_{W^{{}^{\prime}}}$ plane is disallowed from the branching ratio
measurements of $h\to\gamma\gamma$ channel.

###  III.4 Collider Constraints on Heavy Gauge Bosons

Due to its extended gauge sector, the model admits heavy charged and neutral
gauge bosons. The $Z^{\prime}$ couples to different SM fields primarily
through gauge mixing. As a result, this $Z^{\prime}$, if produced in
colliders, can be discovered from its decays to different SM final states.
Specifically, the $Z^{\prime}$ can decay into $\ell_{L}\bar{\ell}_{L}$ ($\ell$
represents all possible charged leptons), $\nu_{L}\bar{\nu}_{L}$ and
$q_{L}\bar{q}_{L}$ ($q$ represents all possible SM coloured particles) pairs
where the the corresponding coupling is suppressed by the parameter $x$
_w.r.t_ the SM counterpart. One can notice in Tables 11 and 12 that the
$Z^{\prime}$ specifically couples to the left-handed chiral fermions since it
is induced purely from a gauge mixing in the neutral gauge sector and does not
receive any contribution from the $U(1)_{2}$ (see Eqn. 9). In addition to the
fermion pairs, the $Z^{\prime}$ can decay into two other SM final states $Zh$
and $WW$. Apart from these SM final states the $Z^{\prime}$ boson does couple
to the dark sector particles $\phi^{+}\phi^{-}$ and $\phi_{0}\phi_{3}$. The
three point vertices corresponding to these two decay modes are proportional
to $\frac{1}{x}$ at leading order (see Table 10). As a consequence the decay
width of these two modes are significantly larger than SM final states. In
Table 1 we present the analytic formulas for these decay widths.

Decay Modes | Decay rate expressions
---|---
$\Gamma\left(Z^{{}^{\prime}}\to f\bar{f}\right)$ | $\frac{N_{c}\lambda^{2}_{Z^{{}^{\prime}}f\overline{f}}}{48\pi}m_{Z^{{}^{\prime}}}$
$\Gamma\left(Z^{{}^{\prime}}\to W^{+}W^{-}\right)$ | $\frac{\lambda^{2}_{WWZ^{{}^{\prime}}}m_{Z^{{}^{\prime}}}}{192\pi}\sqrt{1-\frac{m^{2}_{W}}{m^{2}_{Z^{{}^{\prime}}}}}\left(1+16\frac{m^{2}_{W}}{m^{2}_{Z^{{}^{\prime}}}}-68\frac{m^{4}_{W}}{m^{4}_{Z^{{}^{\prime}}}}-48\frac{m^{6}_{W}}{m^{6}_{Z^{{}^{\prime}}}}\right)$
$\Gamma\left(Z^{{}^{\prime}}\to Zh\right)$ | $\frac{\lambda^{2}_{hZZ^{{}^{\prime}}}}{48\pi m_{Z^{{}^{\prime}}}}\left(2+\frac{\left(m^{2}_{Z^{{}^{\prime}}}+m^{2}_{Z}-m^{2}_{h}\right)^{2}}{4m^{2}_{Z}m^{2}_{Z^{{}^{\prime}}}}\right)\sqrt{1+\frac{\left(m^{2}_{h}-m^{2}_{Z}\right)^{2}}{m^{4}_{Z^{{}^{\prime}}}}-2\frac{\left(m^{2}_{h}+m^{2}_{Z}\right)^{2}}{m^{4}_{Z^{{}^{\prime}}}}}$
$\Gamma\left(Z^{{}^{\prime}}\to\phi_{0}\phi_{3}\right)$ | $\frac{\lambda^{2}_{\phi_{0}\phi_{3}Z^{{}^{\prime}}}}{48\pi m_{Z^{{}^{\prime}}}}\left(m^{2}_{Z^{{}^{\prime}}}-2\left(m^{2}_{\phi_{0}}+m^{2}_{\phi_{3}}\right)\right)\sqrt{1+\frac{\left(m^{2}_{\phi_{0}}-m^{2}_{\phi_{3}}\right)}{m^{4}_{Z^{{}^{\prime}}}}-2\frac{\left(m^{2}_{\phi_{0}}+m^{2}_{\phi_{3}}\right)}{m^{2}_{Z^{{}^{\prime}}}}}$
$\Gamma\left(Z^{{}^{\prime}}\to\phi^{+}\phi^{-}\right)$ | $\frac{\lambda^{2}_{\phi^{+}\phi^{-}Z^{{}^{\prime}}}}{192\pi m_{Z^{{}^{\prime}}}}\left(m^{2}_{Z^{{}^{\prime}}}-4m^{2}_{\phi^{\pm}}\right)\sqrt{1-4\frac{m^{2}_{\phi^{\pm}}}{m^{2}_{Z^{{}^{\prime}}}}}$
Table 1: The relevant $Z^{\prime}$ decay rates for the purposes of this
section - the expressions for the couplings given herein can be found in
Appendix B.

Here $N_{c}$ stands for the colour factor (3 and 1 for quarks and leptons
respectively) and the $\lambda$’s denote the relevant cubic vertices. We have
neglected the mass of the fermions in comparison with the $Z^{\prime}$ mass
for the fermionic decay modes. In principle, there exists a cubic vertex of
the form $\lambda_{Z^{{}^{\prime}}WW^{{}^{\prime}}}$. However, in the present
model $m_{Z^{{}^{\prime}}}$ and $m_{W^{{}^{\prime}}}$ are degenerate and the
$Z^{{}^{\prime}}$ decay via $WW^{{}^{\prime}}$ channel is kinematically
forbidden. In Fig. 3 (left) we present the corresponding BRs of these channels
for $m_{Z^{{}^{\prime}}}$ranging from 1 TeV to 5 TeV. As expected, the BR of
the $\phi^{+}\phi^{-}$ and $\phi_{0}\phi_{3}$ modes dominate for the entire
$m_{Z^{{}^{\prime}}}$ range. Among the SM final states the $WW$ has the
maximum and $Zh$ the minimum branching ratio. For the present calculation we
have set $r$ = 0.25 and $m_{\phi_{0}}=m_{\phi_{3}}=m_{\phi^{\pm}}$ = 400 GeV.
In Fig. 3 (right) we present the
$\frac{\Gamma_{Z^{{}^{\prime}}}}{m_{Z^{{}^{\prime}}}}$ in percentage. There is
a visible growth in this ratio at the mass range $m_{Z^{{}^{\prime}}}$ ¿ 2.5
TeV. However, the overall value remains below $\lesssim$ 30$\%$ for the entire
range. This is to emphasize that for the current benchmark choice the
perturbative treatment is applicable. For different parameter choices, one
should explicitly calculate the
$\frac{\Gamma_{Z^{{}^{\prime}}}}{m_{Z^{{}^{\prime}}}}$ and ensure that the
corresponding value is less than unity.

Figure 3: Left: The branching ratio of $Z^{{}^{\prime}}$ for all possible two
body final states. Here we have set $r$ = 0.25 and
$m_{\phi_{0}}=m_{\phi_{3}}=m_{\phi^{\pm}}$ = 400 GeV. The branching ratio for
$Z^{{}^{\prime}}\to qq~{}$ mode (solid orange line) is summed over all six SM
quark final states. Similarly, the $Z^{{}^{\prime}}\to\ell^{+}\ell^{-}$ (blue
solid) and $Z^{{}^{\prime}}\to\nu\nu~{}$ (green dashed) modes are calculated
while considering all three lepton flavours. Right: The ratio
$\frac{\Gamma_{Z^{{}^{\prime}}}}{m_{Z^{{}^{\prime}}}}$ for
$m_{Z^{{}^{\prime}}}$ mass ranging from 1 TeV to 5 TeV. Here
$\Gamma_{Z^{{}^{\prime}}}$ is the total decay width which is computed
considering all possible two body decay modes.

From the above discussion, it is evident that the $Z^{{}^{\prime}}$ in this
model primarily decays via exotic scalars coming from the scalar bidoublet and
that the branching ratios in the SM final states are significantly low.
Nevertheless, it is useful to check the direct detection limits to understand
the allowed parameter space. Before the LHC era, the LEP and Tevatron
colliders have extensively searched for the heavy neutral gauge bosons. The
LEP searches Electroweak (2003) can be divided into two categories: i)
precision measurement around the $Z$ pole which can effectively constraint the
mixing angle between the $Z-Z^{{}^{\prime}}$ and ii) direct searches for the
heavy neutral boson in the $e^{+}e^{-}\to f\bar{f}$ (here $f$ is all possible
SM fermions) at the mass range above $Z$-pole. On other hand, at Tevatron the
$Z^{{}^{\prime}}$ boson was primarily looked for in the $p\bar{p}\to
Z^{{}^{\prime}}X\to\ell^{+}\ell^{-}X$ channel and using the data, both the
CDFAbe et al. (1995, 1997) and D0 Abachi et al. (1996); Abazov et al. (2001)
collaborations have put bounds on $m_{Z^{{}^{\prime}}}$ for different $U(1)$
extensions of the SM. The commonality between the LEP and Tevatron is the
corresponding bounds are strictly dependent on the $Z^{{}^{\prime}}$ couplings
to SM fermions. For a comprehensive analysis on how these results can be
translated in a specific model, the interested reader can consult Ref. Carena
et al. (2004); Bandyopadhyay et al. (2018). At LHC, both the Drell-Yan (DY)
and Vector Boson Fusion (VBF) production mechanisms are considered for the
$Z^{{}^{\prime}}$ searches. After production, the $Z^{{}^{\prime}}$ can decay
to dileptons, $Zh$, and $WW$ final states. For the present study, we will
consider these three final states and will recast the corresponding ATLAS and
CMS bounds in our model. To calculate the signal cross section in our model we
have used the prescription illustrated in Ref. Pappadopulo et al. (2014);
Flórez et al. (2017). We begin with the process $pp\to Z^{{}^{\prime}}Zh$,
where the $Z^{{}^{\prime}}$ is produced via DY and decays into $Zh$. The $Z$
further decays leptonically and the Higgs boson decays to $\bar{b}b$. Both the
ATLAS Aad et al. (2023) and CMS Sirunyan et al. (2021) have searched the
$Z^{{}^{\prime}}$ boson in this mode. The model dependent upper bounds on the
$\sigma\times\text{BR}$ is displayed in Fig. 4 (upper panel left) by red
dashed (ATLAS) and blue dotted curve (CMS) curves. The black solid line
represents the corresponding number for the present model.

Figure 4: The direct search limits in the channels $\text{DY}\to
Z^{{}^{\prime}}\to Zh$ (upper panel left), $\text{VBF}\to Z^{{}^{\prime}}jj\to
WWjj$ (upper panel right) and $\text{DY}\to\ell^{+}\ell^{-}$ (lower panel).
The red-dashed, blue-dotted and black solid line denote the ATLAS, CMS and the
model specific curves respectively. The arXiv IDs which are written beside the
ATLAS and CMS highlight the corresponding articles Sirunyan et al. (2021); Aad
et al. (2023); Aaboud et al. (2018); Tumasyan et al. (2022); Aaboud et al.
(2017); Sirunyan et al. (2018) which are used in this calculation.

In the case of VBF mechanism the $Z^{{}^{\prime}}$ boson can decay via $WW$
mode. After that, one of the $W$ bosons decays hadronically and the other one
decays leptonically. In Fig. 4 (upper panel right), we display the
corresponding ATLAS, CMS and model specific $\sigma\times\text{BR}$ values via
red dashed, blue dashed and black solid line respectively. We have also
translated the LHC limits in dilepton final state where the $Z^{{}^{\prime}}$
is produced via DY mechanism. The corresponding limits Aaboud et al. (2017);
Sirunyan et al. (2018) are displayed in Fig. 4 (lower panel). Considering
these search channels one can see that the most stringent bound can be
obtained from the $Zh$ mode and to satisfy current collider limits one have to
set $m_{Z^{{}^{\prime}}}\gtrsim$ 2.5 TeV for the specific benchmark points
that we have chosen. Apart from the heavy neutral gauge bosons the LHC also
looked for the additional charged gauge bosons. However the number of searches
is comparatively lower as well as less stringent in the latter case.

##  IV Dark Matter Phenomenology

In this section we consider our model described in Sec. II as a plausible
solution to the dark matter problem. In what follows, we discuss the choice of
the DM candidate, calculate the relic abundance for a multicomponent scenario,
and address the direct detection constraints.

###  IV.1 Choice of Dark Matter Candidate

For a generic parameter choices in this model, the couplings mentioned in
Table 10 allow the decays
$\phi_{0}\rightarrow\phi^{\pm}W^{\prime\mp},\phi_{3}Z^{\prime}$ for a positive
$\lambda_{23}+\tilde{\lambda}_{23}$, and
$\phi_{3}\rightarrow\phi_{0}Z^{\prime}$ for negative values of
$\lambda_{23}+\tilde{\lambda}_{23}$. Such decays, despite being suppressed by
the heavy gauge boson mass, would still make one of the scalars unstable. In
other words, constraining such decays with the age of universe would push the
heavy gauge boson masses to very high values essentially decoupling them from
SM gauge sector for all practical purposes. While this is certainly a
reasonable path to take, in this work our aim to to treat the multi-component
DM scenario while maintaining a meaningful way of discussing the gauge sector
in the context of present experiments and constraints. Thus, we work under the
mass degeneracy limit of the $Z_{2}$ odd sector that makes the neutral
components cosmologically stable DM candidates. Referring to Eqn. 16, we see
that this corresponds to the particular parameter choice
$\tilde{\lambda}_{23}=-\lambda_{23}$. The charged components will completely
annihilate via photon couplings. For the mass degenerate scenario, all
components of the $Z_{2}$ odd sector are stable at tree level in the context
of collider production as well. However, in the absence of phenomenologically
interesting higher dimensional terms in our scalar Lagrangian, loop level
effects can give rise to a sub-GeV mass non-degeneracy between the charged and
neutral components Cirelli et al. (2006) leading to charged scalar decays via
MeV scale charged pions as in reference Belyaev et al. (2018), thus avoiding
any constraints pertaining to the production of an otherwise stable charged
particle in colliders which could potentially disallow charged scalar masses
(and thus the DM masses) upto a few TeV.

###  IV.2 Relic Abundance

In the early universe all the particles were in equilibrium in the high
temperature cosmic soup. Over time as the universe expanded and the
temperature of the universe started to drop, the interaction rate of WIMPs
with the other particles of the cosmic plasma became comparable to the Hubble
parameter. It became difficult for the WIMPs to find another WIMP to
annihilate with and thus they were frozen out remaining in the universe as a
relic which one can compare with the observed DM relic abundance. The value of
relic abundance as measured by Planck satellite is given by $\Omega
h^{2}=0.120\pm 0.001$ Aghanim et al. (2020b). There are two DM candidates in
this model as argued in Sec. IV.1. The dark matter candidates can annihilate
mainly via higgs mediated channels to SM states, for example $b\bar{b}$,
$\tau\bar{\tau}$, $W^{+}W^{-}$, $ZZ$ and $hh$, as allowed by the kinematics in
various DM mass range. The Feynman diagrams are given in Appendix D. Among the
individual contributions of the DM candidates to the relic, the contribution
from $\phi_{3}$ is much less than that of $\phi_{0}$ for the DM mass range
starting from $m_{W}$ upto around 300 GeV. This is due to the fact that in
this model there exists a $\phi_{3}\phi_{3}W^{+}W^{-}$ coupling (see Table 6)
which is absent for $\phi_{0}$, and thus $\phi_{3}$ annihilates more through
that channel leaving less relic in the universe. This coupling comes into
effect around a DM mass of $m_{W}$ and finally goes away as the co-
annihilation channels start to dominate at a DM mass of around 300 GeV.111The
exact DM mass value can be different than 300 GeV depending upon relevant
couplings and the heavy gauge boson mass. As we move into the higher DM mass
range, co-annihilation among $\phi_{0}$ and $\phi_{3}$ via $Z^{\prime}$
mediated diagrams into $WW$, $Zh$ and quark and lepton pairs start to dictate
the relic abundance. All these properties are evident in Fig. 5 where the
relic abundance contribution of each dark matter candidate as well as the
total contribution have been displayed as a function of dark matter mass for a
mass degenerate DM scenario. The $Z^{\prime}$ mass is fixed at 5 TeV as a
benchmark choice and the relic varies with scalar sector parameter
$\lambda_{12}$ in the annihilation dominated low DM mass region whereas at
high DM mass region the $\lambda_{12}$ dependence is absent since the
$Z^{\prime}$ mediated co-annihilation starts dictating the relic. There are
two peaks in the velocity averaged cross section due to the Breit-Wigner
resonance which are related to the mass of the two mediators $h$ and
$Z^{\prime}$. One can also observe the sudden decline in the $\phi_{0}$ relic
curve around $m_{W}$ as the Higgs mediated $WW$ channel becomes kinematically
accessible and similarly for the Higgs mediated $ZZ$ and $hh$ channels. Both
$\phi_{0}$ and $\phi_{3}$ annihilation cross sections are dominated by
$\bar{b}b$ channel until the $WW$ becomes kinematically accessible. However,
for $\phi_{0}$ annihilation a combination of $WW$ and $hh$ channels dominate
at higher DM masses in contrast to $\phi_{3}$ annihilation where the $WW$
channel completely takes over. The co-annihilation is also dominated by the
$WW$ channel once it becomes kinematically available. The region above the
black line is disallowed by overabundance as per Planck data. The relevant
Boltzmann equations are given in Appendix D.

Figure 5: Relic abundance as a function of dark matter mass for various
couplings (lagrangian parameters) in the multi-component dark matter scenario
with annihilation via $h$ portal and co-annihilation via $Z^{\prime}$ portal
for $m_{Z^{\prime}}=5$ TeV.

Next, we turn to the total relic abundance as a function of DM mass in Fig. 6
but this time for various values of $Z^{\prime}$ mass and two benchmark values
of $\lambda_{12}$ to complement our results in Fig. 5. For the higher DM mass
range not much difference can be noticed in this logarithmic scale plot for
various $m_{Z^{\prime}}$ values apart from the shift in the Breit-Wigner
resonance. However, in the lower DM mass range there is strong competition
from the Higgs mediated channels. A heavier $Z^{\prime}$ mediated process
generates less cross section compared to the Higgs mediated processes in low
DM mass region. However, a lighter $Z^{\prime}$ mediated process can produce
comparable cross section to the higgs mediated ones and thus contribute to the
relic significantly even in the low DM mass range as evident in Fig. 6
(right).

Figure 6: Relic abundance as a function of dark matter mass for various
$Z^{\prime}$ masses in the multi-component dark matter scenario.

###  IV.3 Direct Detection Constraints

Direct detection experiments like Xenon1T Aprile et al. (2018, 2022),
PandaX-4T Meng et al. (2021) and LUX-ZEPLIN (LZ) Aalbers et al. (2022)
constrain the parameter space by non-observation as they provide upper limits
on the scattering cross-section of the DM candidates with the detector nuclei
in the non-relativistic limit. The expression222We have assumed that the quark
contribution to each of the nucleons are approximately equal. for spin
independent scattering cross section of $\phi_{0}$ with the nucleons of the
detector material is given by Basak et al. (2021)

$\sigma_{N,0}^{SI}=\frac{\mu_{0}^{2}}{4\pi
m_{\phi_{0}}^{2}}\bigg{[}\frac{\lambda_{h\phi_{0}\phi_{0}}}{m_{h}^{2}}\bigg{(}\sum_{q=u,d,s}\lambda_{hqq}\frac{m_{N}}{m_{q}}f_{Tq}^{N}+\frac{2}{27}f_{TG}^{N}\sum_{q=c,b,t}\lambda_{hqq}\frac{m_{N}}{m_{q}}\bigg{)}\bigg{]}^{2},$
(22)

where, $\mu_{0}=\frac{m_{\phi_{0}}m_{N}}{m_{\phi_{0}}+m_{N}}$. The matrix
elements ($f^{N}_{Ti}$) are given in reference Ellis et al. (2000).

In addition to the Higgs mediated $\phi_{0}N\rightarrow\phi_{0}N$ scattering
process mentioned above, this model also exhibits a $Z^{\prime}$ mediated
$\phi_{0}N\rightarrow\phi_{3}N$ process. In the mass degeneracy limit of the
DM candidates, both the contributions are to be added together to arrive at
the combined cross-section. The relevant couplings
$\lambda_{\phi_{0}\phi_{3}Z^{\prime}}=-\lambda_{\phi_{3}\phi_{0}Z^{\prime}}$
and
$\lambda_{Z^{\prime}\bar{u}_{L}u_{L}}=-\lambda_{Z^{\prime}\bar{d}_{L}d_{L}}$
are given in Tables 10 and 11. In the center of mass frame the cross section
is given by

$\sigma^{SI}_{N,0}=\frac{\mu^{2}}{32\pi
sm_{\phi}^{2}}\int_{0}^{\pi}\frac{p_{f}}{p_{i}}|\mathcal{M}|^{2}\sin\theta^{\prime}d\theta^{\prime},$
(23)

where $s$ is the Mandelstam variable, $p_{i}$ and $p_{f}$ are the initial and
final momenta respectively and $\theta^{\prime}$ is the relative angle between
them. We have assumed $m_{\phi_{0}}=m_{\phi_{3}}=m_{\phi}$ and the reduced
mass $\mu$ is given by $\mu=\frac{m_{\phi}m_{N}}{m_{\phi}+m_{N}}$. The
amplitude can be written as

$\displaystyle\begin{split}|\mathcal{M}|^{2}=&\frac{8\lambda_{\phi_{0}\phi_{3}Z^{\prime}}^{2}\lambda_{Z^{\prime}\bar{u}_{L}u_{L}}^{2}f_{N}^{2}}{m_{Z^{\prime}}^{4}}\bigg{[}(E_{3}E_{2}+p_{i}p_{f}\cos\theta^{\prime}+E_{1}E_{2}+p_{i}^{2})(E_{3}E_{4}+p_{f}^{2}+E_{1}E_{4}+p_{i}p_{f}\cos\theta^{\prime})\\\
&-(m_{\phi}^{2}+E_{1}E_{3}-p_{i}p_{f}\cos\theta^{\prime})(E_{2}E_{4}-p_{i}p_{f}\cos\theta^{\prime}-m_{N}^{2})\bigg{]},\end{split}$
(24)

where $f_{N}$ is the sum of the matrix elements and $E_{i}$ are the respective
energies of the individual initial and final states of the scattering process
in consideration.

In the multicomponent dark matter scenario, the quantity of relevance to
realise the constraints from direct detection experiments by non-observation
is given by the following effective cross-section Bhattacharya et al. (2017b);
Herrero-Garcia et al. (2017); Bhattacharya et al. (2020),

$(\sigma_{N,i}^{SI})^{eff}=\frac{\Omega_{i}}{\Omega_{tot}}\sigma_{N,i}^{SI},$
(25)

where the DM-nucleon scattering cross-section is modified with each DM
candidate’s respective contribution to the relic abundance and
$\sigma_{N,i}^{SI}$ is the summation of the contributions from equation 22 and
23. Figure 7 describes the experimental constraints coming from Xenon1T,
PandaX-4T and LZ experiments on the effective cross-section of the DM
candidate $\phi_{0}$ as a function of DM mass where LZ is clearly the most
restricting of them all. A similar set of constraints can be obtained for the
$\phi_{3}$ as well following the same procedure.

Figure 7: Constraints from direct detection experiments on the effective
cross-section of dark matter candidate $\phi_{0}$ for benchmark values
$\lambda_{12}=10^{-2}$(left) and $\lambda_{12}=10^{-3}$(right). The
constraints from Xenon1T(2022), PandaX-4T (2021), LUX-ZEPLIN (2022) and also
the Neutrino Floor have been shown in the figures.

It is evident from Fig. 7 that for a $Z^{\prime}$ mass of 2.5 TeV (motivated
by the lower bound from Fig. 4), much of the DM mass parameter space is
already excluded upto around 2 TeV for both the benchmark values of coupling
parameter $\lambda_{12}=10^{-2}$ and $10^{-3}$. Keeping every other benchmark
value same and modifying only the $Z^{\prime}$ mass to 5 TeV brings the DM
mass exclusion region down to $\sim 100-200$ GeV range and very interestingly
opens up some parameter space in the lower DM mass region of around $10-15$
GeV for $\lambda_{12}=10^{-3}$. A more careful look at the figures makes it
apparent that for $\lambda_{12}=10^{-2}$, the Higgs mediated diagram plays an
important role upto the DM masses around the weak scale and beyond that the
heavy gauge boson mediated diagram takes over. However, the contribution of
the Higgs mediated diagram remains subdominant to the one mediated by the
heavy gauge boson throughout the DM mass range for $\lambda_{12}=10^{-3}$. As
expected, the cross section is more suppressed by a higher value of the
$Z^{\prime}$ mass compared to a lower value. The slightly raised plateau like
shape of the curves in the mid-range of DM masses depicted here are the result
of the fractional contributions of the individual DM components to the relic.
One also has to consider the neutrino floor Billard et al. (2022) and the
neutrino fog O’Hare (2021) in this regard as for cross-sections too small the
DM scattering will be difficult to tell apart due to the presence of neutrino
scattering. Similar results are expected for $\phi_{3}$ upto the relic
fraction as is realized in Fig. 8. The rise (dip) in the curves are indicative
of the substantial (negligible) contribution of $\phi_{3}$ to the relic in
that DM mass region.

Figure 8: Same as figure 7 but for $\phi_{3}$.

There is also in principle a scope for inelastic scattering between the two DM
components by softly breaking their mass degeneracy. However, that would lead
to some problems with the stability of the DM components and this
multicomponent DM scenario will effectively reduce to a single component
scenario as discussed in Sec. IV.1.

##  V Conclusion

In this work we have extended the SM with an SU(2) gauge sector and the scalar
sector with a bidoublet and a non-linear sigma field. The neutral components
of the $Z_{2}$ protected bidoublet are suitable DM candidates. In the DM mass
range starting from a few GeV upto around half of the Higgs mass, the
constraint from invisible decay of higgs becomes important and restricts the
most relevant scalar sector coupling parameter $\lambda_{12}$ with an upper
bound of $\mathcal{O}(10^{-3}-10^{-2})$. However, the sensitivity of direct
detection experiments provide a much more stringent bound for a typical scalar
extension under the WIMP scenario. In many scalar extensions, one has to
resort to the Breit-Wigner resonance to account for the DM relic abundance
while simultaneously satisfying the direct detection constraints. In our
multicomponent DM model the extended gauge sector mitigates those strict
constraints to a great extent. The gauge sector couplings provide additional
channels for the DM annihilation and co-annihilation thus increasing the
velocity averaged cross-section and reducing the relic to avoid overabundance.
In the similar coupling range this model can also be constrained with the
direct detection effective cross-section values from the experimental
observations. Furthermore, we found that a $Z^{\prime}$ mediated scattering
diagram can dominate a scalar mediated diagram in a sizeable region of
parameter space. Thus, models with scalar dark matter within an extended gauge
sector context can thus satisfy all experimental constraints and remain
interesting prospects for future study and detection.

###### Acknowledgements.

BC acknowledges support from the Department of Science and Technology, India,
under Grant CRG/2020/004171. KL thanks Oscar Zapata for helpful comments. AS
thanks Tarak Nath Maity for helpful discussions.

##  Appendix A Perturbativity

The perturbativity of a model typically means that the loop level
contributions towards the scattering amplitude should be lower than the tree
level contribution. It must be ensured that the available couplings in the
model are sufficiently small such that one can employ perturbation theory
reliably to calculate the $s$-matrix. In this model the cubic vertices between
the $Z^{\prime}$ boson and $\phi_{0}\phi_{3}$ or $\phi^{+}\phi^{-}$ scalars
are proportional to $\frac{1}{x}$. As a result, for a very small values of $x$
the decay width $\Gamma(Z^{{}^{\prime}}\to\phi_{0}\phi_{3}/\phi^{+}\phi^{-})$
becomes significantly large. For those values of $x$, the perturbative
treatment is not applicable. From Eq. 14 one can see
$x=\sqrt{\frac{M^{2}_{W}}{r^{2}M^{2}_{Z^{{}^{\prime}}}}}$. This relation
suggest that the parameter $x$ is dependent on two variables $r$ and
$M_{Z^{{}^{\prime}}}$. In fig. 9, we show the variation of the parameter $x$
in $M_{Z^{{}^{\prime}}}$ vs $r$ plane.

Figure 9: The variation of $x$ in $M_{Z^{{}^{\prime}}}$ versus $r$ plane. For
details please see the text.

The black line, blue dashed and red dotted lines signify the contours of $x$
equals to 0.05, 0.1 and 0.5 respectively. From the plot one can notice for the
$r\to 0$ the value of $x$ enhances for the entire range of
$M_{Z^{{}^{\prime}}}$. With these numbers in hand one can evaluate
$\frac{\Gamma_{Z^{{}^{\prime}}}}{M_{Z^{{}^{\prime}}}}$ for the $r$ ranging
from $10^{-1}$ to 1 and $M_{Z^{{}^{\prime}}}$ ranging from 1.25 TeV to 5 TeV.
In fig. 10, the allowed region is given by
$\frac{\Gamma_{Z^{{}^{\prime}}}}{M_{Z^{{}^{\prime}}}}<$ 1, where one can
safely consider perturbativity.

Figure 10: The perturbativity limit is well respected in all regions of
interest. The dark cyan shaded portion represents the disallowed region of
parameter space.

##  Appendix B Couplings

$\lambda_{abc}$ | Couplings
---|---
$hhh$ | $\lambda_{1}v$
$h\phi_{0}\phi_{0}$ | $\frac{\lambda_{12}v}{2}$
$h\phi_{3}\phi_{3}$ | $\frac{\lambda_{12}v}{2}$
$h\phi^{+}\phi^{-}$ | $\lambda_{12}v$
Table 2: Cubic scalar couplings. $\lambda_{abcd}$ | Couplings
---|---
$hhhh$ | $\frac{\lambda_{1}}{4}$
$hh\phi_{0}\phi_{0}$ | $\frac{\lambda_{12}}{4}$
$hh\phi_{3}\phi_{3}$ | $\frac{\lambda_{12}}{4}$
$hh\phi^{+}\phi^{-}$ | $\frac{\lambda_{12}}{2}$
$\phi_{0}\phi_{0}\phi_{0}\phi_{0}$ | $\frac{\lambda_{2}}{4}$
$\phi_{3}\phi_{3}\phi_{3}\phi_{3}$ | $\frac{\lambda_{2}}{4}$
$\phi^{+}\phi^{+}\phi^{-}\phi^{-}$ | $\lambda_{2}$
$\phi_{0}\phi_{0}\phi_{3}\phi_{3}$ | $\frac{\lambda_{2}}{2}$
$\phi^{0}\phi^{0}\phi^{+}\phi^{-}$ | $\lambda_{2}$
$\phi^{3}\phi^{3}\phi^{+}\phi^{-}$ | $\lambda_{2}$
Table 3: Quartic scalar couplings. $\lambda_{abc}$ | Couplings | Conjugate
---|---|---
$hZZ$ | $\frac{e^{2}v}{4}\sec^{2}\theta(\operatorname{cosec}^{2}\theta+x^{2})$ |
$hZ^{\prime}Z^{\prime}$ | $\frac{e^{2}v}{4}x^{2}\operatorname{cosec}^{2}\theta$ |
$hZZ^{\prime}$ | $-\frac{e^{2}v}{2}x\operatorname{cosec}^{2}\theta\sec\theta$ |
$hW^{+}W^{-}$ | $\frac{e^{2}v}{2}(\operatorname{cosec}^{2}\theta+x^{2})$ |
$hW^{\prime+}W^{\prime-}$ | $\frac{e^{2}v}{2}x^{2}\operatorname{cosec}^{2}\theta$ |
$hW^{\prime+}W^{-}$ | $-\frac{e^{2}v}{2}x\operatorname{cosec}^{2}\theta$ | $hW^{+}W^{\prime-}$
Table 4: Cubic higgs-gauge couplings. $\lambda_{abcd}$ | Couplings | Conjugate
---|---|---
$hhZZ$ | $\frac{e^{2}}{8}\sec^{2}\theta(\operatorname{cosec}^{2}\theta+x^{2})$ |
$hhZ^{\prime}Z^{\prime}$ | $\frac{e^{2}}{8}x^{2}\operatorname{cosec}^{2}\theta$ |
$hhZZ^{\prime}$ | $-\frac{e^{2}}{4}x\operatorname{cosec}^{2}\theta\sec\theta$ |
$hhW^{+}W^{-}$ | $\frac{e^{2}}{4}(\operatorname{cosec}^{2}\theta+x^{2})$ |
$hhW^{\prime+}W^{\prime-}$ | $\frac{e^{2}}{4}x^{2}\operatorname{cosec}^{2}\theta$ |
$hhW^{\prime+}W^{-}$ | $-\frac{e^{2}}{4}x\operatorname{cosec}^{2}\theta$ | $hhW^{+}W^{\prime-}$
Table 5: Quartic higgs-gauge couplings. $\lambda_{abcd}$ | Couplings | Conjugate
---|---|---
$\phi_{0}\phi_{0}Z^{\prime}Z^{\prime}$ | $\frac{e^{2}\operatorname{cosec}^{2}\theta}{8x^{2}}+\frac{e^{2}}{8}(3+2\cot^{2}\theta)+\frac{e^{2}}{8}(3+\cot^{2}\theta)x^{2}$ |
$\phi_{0}\phi_{0}W^{\prime+}W^{\prime-}$ | $\frac{e^{2}\operatorname{cosec}^{2}\theta}{4x^{2}}+\frac{e^{2}}{4}(3+2\cot^{2}\theta)+\frac{e^{2}}{4}(3+\cot^{2}\theta)x^{2}$ |
$\phi_{3}\phi_{3}W^{+}W^{-}$ | $e^{2}(\operatorname{cosec}^{2}\theta+x^{2})$ |
$\phi_{3}\phi_{3}W^{\prime+}W^{\prime-}$ | $\frac{e^{2}}{4x^{2}}\operatorname{cosec}^{2}\theta-\frac{e^{2}}{4}(\operatorname{cosec}^{2}\theta+\cot^{2}\theta)+\frac{e^{2}}{4}(\cot^{2}\theta-1)x^{2}$ |
$\phi_{3}\phi_{3}W^{\prime+}W^{-}$ | $\frac{e^{2}}{2x}\operatorname{cosec}^{2}\theta-\frac{e^{2}}{2}(\cot^{2}\theta)x$ | $\phi_{3}\phi_{3}W^{+}W^{\prime-}$
$\phi_{0}\phi_{3}W^{\prime+}W^{-}$ | $\frac{ie^{2}}{2x}\operatorname{cosec}^{2}\theta+\frac{ie^{2}}{2}(\operatorname{cosec}^{2}\theta+1)x$ | $-\phi_{0}\phi_{3}W^{+}W^{\prime-}$
Table 6: Quartic gauge couplings for DM annihilation. $\lambda_{abcd}$ | Couplings | Conjugate
---|---|---
$\phi^{+}\phi_{0}W^{-}Z^{\prime}$ | $-\frac{ie^{2}}{2x}\operatorname{cosec}^{2}\theta-\frac{ie^{2}}{2}(\operatorname{cosec}^{2}\theta+1)x$ | $-\phi^{-}\phi_{0}W^{+}Z^{\prime}$
$\phi^{+}\phi_{3}W^{-}A$ | $-e^{2}\operatorname{cosec}\theta-e^{2}(\sin\theta)x^{2}$ | $\phi^{-}\phi_{3}W^{+}A$
$\phi^{+}\phi_{3}W^{-}Z$ | $-e^{2}\cos\theta(\operatorname{cosec}^{2}\theta+x^{2})$ | $\phi^{-}\phi_{3}W^{+}Z$
$\phi^{+}\phi_{3}W^{-}Z^{\prime}$ | $-\frac{e^{2}}{2x}\operatorname{cosec}^{2}\theta+\frac{e^{2}}{2}(\cot^{2}\theta)x$ | $\phi^{-}\phi_{3}W^{+}Z^{\prime}$
$\phi^{+}\phi_{0}W^{\prime-}A$ | $\frac{ie^{2}}{2x}\operatorname{cosec}\theta+\frac{ie^{2}}{2}(\operatorname{cosec}\theta+\sin\theta)$ | $-\phi^{-}\phi_{0}W^{\prime+}A$
$\phi^{+}\phi_{0}W^{\prime-}Z$ | $\frac{ie^{2}}{2x}\cot\theta\operatorname{cosec}\theta+\frac{ie^{2}}{2}\cos\theta(2+\cot^{2}\theta)x$ | $-\phi^{-}\phi_{0}W^{\prime+}Z$
$\phi^{+}\phi_{3}W^{\prime-}A$ | $-\frac{e^{2}}{2x}\operatorname{cosec}\theta+\frac{e^{2}}{2}(\cos\theta\cot\theta)x$ | $\phi^{-}\phi_{3}W^{\prime+}A$
$\phi^{+}\phi_{3}W^{\prime-}Z$ | $-\frac{e^{2}}{2x}\cot\theta\operatorname{cosec}\theta+\frac{e^{2}}{2}(\cos\theta\cot^{2}\theta)x$ | $\phi^{-}\phi_{3}W^{\prime+}Z$
$\phi^{+}\phi_{3}W^{\prime-}Z^{\prime}$ | $e^{2}(\operatorname{cosec}^{2}\theta+x^{2})$ | $\phi^{-}\phi_{3}W^{\prime+}Z^{\prime}$
Table 7: Quartic couplings for charged scalar decay. $\lambda_{abcd}$ | Couplings | Conjugate
---|---|---
$\phi^{+}\phi^{-}AA$ | $e^{2}+e^{2}\sin^{2}\theta x^{2}$ |
$\phi^{+}\phi^{-}ZZ$ | $e^{2}\cot^{2}\theta+e^{2}\cos^{2}\theta x^{2}$ |
$\phi^{+}\phi^{-}Z^{\prime}Z^{\prime}$ | $\frac{e^{2}}{4x^{2}}\operatorname{cosec}^{2}\theta-\frac{e^{2}}{4}(1+2\cot^{2}\theta)-\frac{e^{2}}{4}(1-\cot^{2}\theta)x^{2}$ |
$\phi^{+}\phi^{-}AZ$ | $2e^{2}\cot\theta+2e^{2}(\sin\theta\cos\theta)x^{2}$ |
$\phi^{+}\phi^{-}AZ^{\prime}$ | $\frac{e^{2}}{x}\operatorname{cosec}\theta-e^{2}(\cos\theta\cot\theta)x$ |
$\phi^{+}\phi^{-}ZZ^{\prime}$ | $\frac{e^{2}}{x}\cot\theta\operatorname{cosec}\theta-e^{2}(\cos\theta\cot^{2}\theta)x$ |
$\phi^{+}\phi^{-}W^{+}W^{-}$ | $e^{2}(\operatorname{cosec}^{2}\theta+x^{2})$ |
$\phi^{+}\phi^{-}W^{\prime+}W^{\prime-}$ | $\frac{e^{2}}{2x^{2}}\operatorname{cosec}^{2}\theta+\frac{e^{2}}{2}+\frac{e^{2}}{2}(\operatorname{cosec}^{2}\theta)x^{2}$ |
$\phi^{+}\phi^{-}W^{\prime+}W^{-}$ | $\frac{e^{2}}{2x}\operatorname{cosec}^{2}\theta-\frac{e^{2}}{2}(\cot^{2}\theta)x$ | $\phi^{+}\phi^{-}W^{+}W^{\prime-}$
Table 8: Quartic couplings for charged scalar annihilation. $\lambda_{abcd}$ | Couplings | Conjugate
---|---|---
$\phi^{+}\phi^{+}W^{-}W^{-}$ | $-\frac{e^{2}}{2}(\operatorname{cosec}^{2}\theta+x^{2})$ | $\phi^{-}\phi^{-}W^{+}W^{+}$
$\phi^{+}\phi^{+}W^{\prime-}W^{\prime-}$ | $\frac{e^{2}}{2}(\operatorname{cosec}^{2}\theta+x^{2})$ | $\phi^{-}\phi^{-}W^{\prime+}W^{\prime+}$
$\phi^{+}\phi^{+}W^{\prime-}W^{-}$ | $-\frac{e^{2}}{2x}\operatorname{cosec}^{2}\theta+\frac{e^{2}}{2}(\cot^{2}\theta)x$ | $\phi^{-}\phi^{-}W^{\prime+}W^{+}$
Table 9: A few more couplings between the scalar and gauge sectors. $\lambda_{abcd}$ | Couplings | Conjugate
---|---|---
$(\partial\phi_{0})\phi_{3}Z^{\prime}$ | $-\frac{e}{2x}\operatorname{cosec}\theta-\frac{e}{4\sec\theta}(2\cot\theta+3\tan\theta)x$ |
$(\partial\phi_{0})\phi^{+}W^{\prime-}$ | $-\frac{e}{2x}\operatorname{cosec}\theta-\frac{e}{4\sec\theta}(2\cot\theta+3\tan\theta)x$ | $(\partial\phi_{0})\phi^{-}W^{\prime+}$
$(\partial\phi_{3})\phi_{0}Z^{\prime}$ | $\frac{e}{2x}\operatorname{cosec}\theta+\frac{e}{4\sec\theta}(2\cot\theta+3\tan\theta)x$ |
$(\partial\phi_{3})\phi^{+}W^{-}$ | $-ie\operatorname{cosec}\theta-\frac{ie}{2}(\sin\theta)x^{2}$ | $-(\partial\phi_{3})\phi^{-}W^{+}$
$(\partial\phi_{3})\phi^{+}W^{\prime-}$ | $-\frac{ie}{2x}\operatorname{cosec}\theta+\frac{ie}{4\sec\theta}(2\cot\theta+\tan\theta)x$ | $-(\partial\phi_{3})\phi^{-}W^{\prime+}$
$(\partial\phi^{+})\phi_{0}W^{\prime-}$ | $\frac{e}{2x}\operatorname{cosec}\theta+\frac{e}{4\sec\theta}(2\cot\theta+3\tan\theta)x$ | $(\partial\phi^{-})\phi_{0}W^{\prime+}$
$(\partial\phi^{+})\phi_{3}W^{-}$ | $ie\operatorname{cosec}\theta+\frac{ie}{2}(\sin\theta)x^{2}$ | $-(\partial\phi^{-})\phi_{3}W^{+}$
$(\partial\phi^{+})\phi_{3}W^{\prime-}$ | $\frac{ie}{2x}\operatorname{cosec}\theta-\frac{ie}{4\sec\theta}(2\cot\theta+\tan\theta)x$ | $-(\partial\phi^{-})\phi_{3}W^{\prime+}$
$(\partial\phi^{+})\phi^{-}A$ | $-ie-\frac{ie}{2}(\sin^{2}\theta)x^{2}$ | $-(\partial\phi^{-})\phi^{+}A$
$(\partial\phi^{+})\phi^{-}Z$ | $-ie\cot\theta-\frac{ie}{2}(\sin\theta\cos\theta)x^{2}$ | $-(\partial\phi^{-})\phi^{+}Z$
$(\partial\phi^{+})\phi^{-}Z^{\prime}$ | $-\frac{ie}{2x}\operatorname{cosec}\theta+\frac{ie}{4\sec\theta}(2\cot\theta+\tan\theta)x$ | $-(\partial\phi^{-})\phi^{+}Z^{\prime}$
Table 10: Cubic Couplings of gauge fields and BSM scalars involving derivatives. $\lambda_{abc}$ | Couplings
---|---
$Z\bar{u}_{L}u_{L}$ | $\frac{g_{1}\cos\theta}{2}-g_{2}\tilde{Y}_{u}\sin\theta$
$Z^{\prime}\bar{u}_{L}u_{L}$ | $-\frac{g_{1}x}{2}$
$A\bar{u}_{L}u_{L}$ | $\frac{g_{1}\sin\theta}{2}+g_{2}\cos\theta\tilde{Y}_{u}$
$Z\bar{d}_{L}d_{L}$ | $-\frac{g_{1}\cos\theta}{2}-g_{2}\tilde{Y}_{d}\sin\theta$
$Z^{\prime}\bar{d}_{L}d_{L}$ | $\frac{g_{1}x}{2}$
$A\bar{d}_{L}d_{L}$ | $-\frac{g_{1}\sin\theta}{2}+g_{2}\cos\theta\tilde{Y}_{d}$
$Z\bar{u}_{R}u_{R}$ | $-g_{2}\sin\theta\tilde{Y}_{u}$
$Z^{\prime}\bar{u}_{R}u_{R}$ | $0$
$A\bar{u}_{R}u_{R}$ | $g_{2}\cos\theta\tilde{Y}_{u}$
$Z\bar{d}_{R}d_{R}$ | $-g_{2}\sin\theta\tilde{Y}_{d}$
$Z^{\prime}\bar{d}_{R}d_{R}$ | $0$
$A\bar{d}_{R}d_{R}$ | $g_{2}\cos\theta\tilde{Y}_{d}$
$Wud$ | $\frac{g_{1}}{\sqrt{2}}V_{CKM}$
$W^{\prime}ud$ | $-\frac{g_{1}x}{\sqrt{2}}V_{CKM}$
Table 11: Cubic Couplings of gauge fields with quarks. $\lambda_{abc}$ | Couplings
---|---
$Z\bar{e}_{R}e_{R}$ | $-g_{2}\sin\theta\tilde{Y}_{e}$
$Z^{\prime}\bar{e}_{R}e_{R}$ | $0$
$A\bar{e}_{R}e_{R}$ | $g_{2}\cos\theta\tilde{Y}_{e}$
$Z\bar{\nu}_{L}\nu_{L}$ | $\frac{g_{1}\cos\theta}{2}-g_{2}\sin\theta\tilde{Y}_{l}$
$Z^{\prime}\bar{\nu}_{L}\nu_{L}$ | $-\frac{g_{1}x}{2}$
$A\bar{\nu}_{L}\nu_{L}$ | $\frac{g_{1}\sin\theta}{2}+g_{2}\cos\theta\tilde{Y}_{l}$
$W\bar{e}_{L}\nu_{L}$ | $\frac{g_{1}}{\sqrt{2}}$
$W^{\prime}\bar{e}_{L}\nu_{L}$ | $-\frac{g_{1}x}{\sqrt{2}}$
Table 12: Cubic Couplings of gauge fields with leptons. $\lambda_{abc}$ | Couplings
---|---
$WW\gamma$ | $e\left(1+x^{2}\right)$
$W^{{}^{\prime}}W^{{}^{\prime}}\gamma$ | $e\left(1+x^{2}\right)$
$WW^{{}^{\prime}}\gamma$ | 0
$WWZ$ | $e\cot\theta_{w}\left(1+x^{2}\right)$
$W^{{}^{\prime}}W^{{}^{\prime}}Z^{{}^{\prime}}$ | $\frac{e}{x\sin\theta_{w}}\left(1-x^{4}\right)$
$WWZ^{{}^{\prime}}$ | $-\frac{ex}{\sin\theta_{w}}\left(1+x\right)$
Table 13: Relevant couplings of gauge bosons among themselves.

##  Appendix C $h\to\gamma\gamma$ Decay Width

The CP-even Higgs boson can decay into di-photon channel through loop mediated
processes. In case of SM, the $W$ boson loop and different fermion loops
participate in this process. For the present model, in addition to these loop
diagrams the heavy charged gauge boson $W^{{}^{\prime}}$ and the charged
scalar $\phi^{\pm}$ can also engage in the loop induced processes. In Fig. 11
we show the all the Feynman diagrams that can contribute to these process.
Ref. Gunion et al. (2000) have calculated this decay width for MSSM scenario
and adopting their formalism one can reproduce the width for the current
model. The analytic expression for this decay width is given in Eq. 26

$\Gamma\left(h\to\gamma\gamma\right)=\frac{\alpha^{2}g^{2}}{1024\pi^{3}}\frac{m^{3}_{h}}{m^{2}_{W}}\left|\sum_{i}\mathcal{I}^{i}_{h}\right|^{2}$
(26)

Figure 11: All possible loop mediated processes that can contribute to
$h\to\gamma\gamma$ decay mode.

where $\mathcal{I}^{i}_{h}$ represent the individual loop contribution from
different particles that flowing in the loop. In Eq. 27, we present the
explicit expression for $\mathcal{I}^{i}_{h}$.

$\begin{split}\mathcal{I}^{f}_{h}&=N_{c}Q^{2}_{f}R^{h}_{f}F^{h}_{\frac{1}{2}}\left(\tau_{f}\right)\\\
\mathcal{I}^{W/W^{{}^{\prime}}}_{h}&=R^{h}_{W/W^{{}^{\prime}}}F_{1}\left(\tau_{W/W^{{}^{\prime}}}\right)\\\
\mathcal{I}^{\phi^{\pm}}_{h}&=R^{h}_{\phi^{\pm}}\frac{M^{2}_{W}}{M^{2}_{\phi^{\pm}}}F_{0}\left(\tau_{\phi^{\pm}}\right)\end{split}$
(27)

Here $N_{c}$ and $Q_{f}$ represents the color factor and electromagnetic
charges for each SM fermions. The $R^{h}_{i}$s are the relative coupling
strength of the three point vertices $\lambda_{hii}$ (where $i$ can vary based
on the particle flowing in the loop) _w.r.t_ corresponding SM coupling. For
the present model $R^{h}_{f}=1$ as the fermion sector is identical to the SM.
In case of gauge bosons and the charged scalar the relative coupling strength
would be $R^{h}_{W}=\left(1+x^{2}\sin^{2}\theta_{W}\right)$,
$R^{h}_{W^{{}^{\prime}}}=x^{2}$ and
$R^{h}_{\phi^{\pm}}=\frac{2\lambda_{ab}}{g^{2}}$ respectively. To get the
sense of these values one should the check the corresponding vertices that are
present in appendix B. The $F_{j}\left(\tau_{i}\right)$ are the form factors
corresponding to the gauge, fermion and scalar loops which can take following
form.

$\begin{split}F_{0}&=\tau\left[1-\tau f\left(\tau\right)\right]\\\
F_{\frac{1}{2}}&=-2\tau\left[1+\left(1-\tau\right)f\left(\tau\right)\right]\\\
F_{1}&=2+3\tau+3\tau\left(2-\tau\right)f\left(\tau\right)\end{split}$ (28)

where $\tau_{i}=\frac{4m^{2}_{i}}{m^{2}_{h}}$ and depending on the value of
$\tau$ the $f\left(\tau\right)$ can take following form

$\begin{split}f\left(\tau\right)&=\left[\sin^{-1}\left(\sqrt{\frac{1}{\tau}}\right)\right]^{2}~{}~{}~{}\text{if}~{}~{}~{}\tau\geq
1\\\ \text{or}\\\
f\left(\tau\right)&=-\frac{1}{4}\left[\text{ln}\left(\frac{1+\sqrt{1-\tau}}{1+\sqrt{1-\tau}}\right)-i\pi\right]^{2}~{}~{}~{}\text{if}~{}~{}\tau\leq
1\end{split}$ (29)

##  Appendix D Boltzmann Equations

The coupled Boltzmann equations for the multicomponent DM are given by the
following equations :

$\frac{dn_{0}}{dt}+3Hn_{0}=-\langle\sigma
v\rangle_{00}(n_{0}^{2}-n_{0_{eq}}^{2})-\langle\sigma
v\rangle_{03}(n_{0}n_{3}-n_{0_{eq}}n_{3_{eq}})$ (30)
$\frac{dn_{3}}{dt}+3Hn_{3}=-\langle\sigma
v\rangle_{33}(n_{3}^{2}-n_{3_{eq}}^{2})-\langle\sigma
v\rangle_{03}(n_{0}n_{3}-n_{0_{eq}}n_{3_{eq}})$ (31)

where $n_{i}$ is the number density of $\phi_{i}$ and $\langle\sigma
v\rangle_{ij}$ is the velocity averaged cross section of annihilation/co-
annihilation of $\phi_{i}$ and $\phi_{j}$. The expression for the cross
sections can be found in many references including Guo and Wu (2010); Basak
and Mondal (2014); Nam (2020). The corresponding Feynman diagrams for
$\phi_{0}$ annihilation are given in figure 12. The $\phi_{3}$ annihilation
channels are same as the $\phi_{0}$ channels except the additional channel of
figure 13. The Feynman diagrams corresponding to the co-annihilation channels
are given in figure 14.

Figure 12: Relevant $\phi_{0}$ annihilation channels for relic abundance
calculation. Figure 13: The extra $\phi_{3}$ annihilation channel for relic
abundance calculation compared to $\phi_{0}$ annihilation channels.

Figure 14: Relevant co-annihilation channels for relic abundance calculation.

These equations have been solved to find the relic Kolb and Turner (1990);
Srednicki et al. (1988); Edsjo and Gondolo (1997) by defining an effective
velocity averaged cross section in a similar way as in ref. Ghosh et al.
(2022) considering appropriate channels of production.

## References

  * Aghanim et al. (2020a) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020a), [Erratum: Astron.Astrophys. 652, C4 (2021)], eprint 1807.06209.
  * Hinshaw et al. (2013) G. Hinshaw et al. (WMAP), Astrophys. J. Suppl. 208, 19 (2013), eprint 1212.5226.
  * Bertone et al. (2005) G. Bertone, D. Hooper, and J. Silk, Phys.Rept. 405, 279 (2005), eprint hep-ph/0404175.
  * Bergstrom (2009) L. Bergstrom, New J.Phys. 11, 105006 (2009), eprint 0903.4849.
  * Arcadi et al. (2018) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, Eur. Phys. J. C 78, 203 (2018), eprint 1703.07364.
  * McDonald (1994) J. McDonald, Phys.Rev. D50, 3637 (1994), eprint hep-ph/0702143.
  * Burgess et al. (2001) C. Burgess, M. Pospelov, and T. ter Veldhuis, Nucl.Phys. B619, 709 (2001), eprint hep-ph/0011335.
  * Guo and Wu (2010) W.-L. Guo and Y.-L. Wu, JHEP 1010, 083 (2010), eprint 1006.2518.
  * Bandyopadhyay et al. (2010) A. Bandyopadhyay, S. Chakraborty, A. Ghosal, and D. Majumdar, JHEP 1011, 065 (2010), eprint 1003.0809.
  * He et al. (2009) X.-G. He, T. Li, X.-Q. Li, J. Tandean, and H.-C. Tsai, Phys. Rev. D 79, 023521 (2009), eprint 0811.0658.
  * He et al. (2010) X.-G. He, T. Li, X.-Q. Li, J. Tandean, and H.-C. Tsai, Phys. Lett. B 688, 332 (2010), eprint 0912.4722.
  * Cline et al. (2013) J. M. Cline, K. Kainulainen, P. Scott, and C. Weniger, Phys. Rev. D 88, 055025 (2013), [Erratum: Phys.Rev.D 92, 039906 (2015)], eprint 1306.4710.
  * Abada et al. (2011) A. Abada, D. Ghaffor, and S. Nasri, Phys. Rev. D 83, 095021 (2011), eprint 1101.0365.
  * Abada and Nasri (2012) A. Abada and S. Nasri, Phys. Rev. D 85, 075009 (2012), eprint 1201.1413.
  * Arhrib and Maniatis (2019) A. Arhrib and M. Maniatis, Phys. Lett. B 796, 15 (2019), eprint 1807.03554.
  * Hamada et al. (2021) Y. Hamada, H. Kawai, K.-y. Oda, and K. Yagyu, JHEP 01, 087 (2021), eprint 2008.08700.
  * Modak et al. (2015) K. P. Modak, D. Majumdar, and S. Rakshit, JCAP 03, 011 (2015), eprint 1312.7488.
  * Maniatis (2021) M. Maniatis, Phys. Rev. D 103, 015010 (2021), eprint 2005.13443.
  * Bhattacharya et al. (2017a) S. Bhattacharya, P. Ghosh, T. N. Maity, and T. S. Ray, JHEP 10, 088 (2017a), eprint 1706.04699.
  * Maity and Ray (2020) T. N. Maity and T. S. Ray, Phys. Rev. D 101, 103013 (2020), eprint 1908.10343.
  * Basak et al. (2021) T. Basak, B. Coleppa, and K. Loho, JHEP 06, 104 (2021), eprint 2105.09044.
  * Lopez Honorez et al. (2007) L. Lopez Honorez, E. Nezri, J. F. Oliver, and M. H. G. Tytgat, JCAP 02, 028 (2007), eprint hep-ph/0612275.
  * Belyaev et al. (2018) A. Belyaev, G. Cacciapaglia, I. P. Ivanov, F. Rojas-Abatte, and M. Thomas, Phys. Rev. D 97, 035011 (2018), eprint 1612.00511.
  * Rojas-Abatte et al. (2017) F. Rojas-Abatte, M. L. Mora, J. Urbina, and A. R. Zerwekh, Phys. Rev. D 96, 095025 (2017), eprint 1707.04543.
  * Bhat and Adhikari (2020) I. A. Bhat and R. Adhikari, Phys. Rev. D 101, 075030 (2020), eprint 1906.10185.
  * Nam (2020) C. H. Nam, Eur. Phys. J. C 80, 1114 (2020), eprint 2011.11207.
  * Basak and Mondal (2014) T. Basak and T. Mondal, Phys. Rev. D 89, 063527 (2014), eprint 1308.0023.
  * Bhattacharyya and Das (2016) G. Bhattacharyya and D. Das, Pramana 87, 40 (2016), eprint 1507.06424.
  * ATL (2020) Tech. Rep. ATLAS-CONF-2020-052, CERN, Geneva (2020), URL http://cds.cern.ch/record/2743055.
  * Sirunyan et al. (2019) A. M. Sirunyan et al. (CMS), Phys. Lett. B 793, 520 (2019), eprint 1809.05937.
  * Electroweak (2003) t. S. Electroweak (LEP, ALEPH, DELPHI, L3, OPAL, LEP Electroweak Working Group, SLD Electroweak Group, SLD Heavy Flavor Group) (2003), eprint hep-ex/0312023.
  * Abe et al. (1995) F. Abe et al. (CDF), Phys. Rev. D 51, R949 (1995).
  * Abe et al. (1997) F. Abe et al. (CDF), Phys. Rev. Lett. 79, 2192 (1997).
  * Abachi et al. (1996) S. Abachi et al. (D0), Phys. Lett. B 385, 471 (1996).
  * Abazov et al. (2001) V. M. Abazov et al. (D0), Phys. Rev. Lett. 87, 061802 (2001), eprint hep-ex/0102048.
  * Carena et al. (2004) M. Carena, A. Daleo, B. A. Dobrescu, and T. M. P. Tait, Phys. Rev. D 70, 093009 (2004), eprint hep-ph/0408098.
  * Bandyopadhyay et al. (2018) T. Bandyopadhyay, G. Bhattacharyya, D. Das, and A. Raychaudhuri, Phys. Rev. D 98, 035027 (2018), eprint 1803.07989.
  * Pappadopulo et al. (2014) D. Pappadopulo, A. Thamm, R. Torre, and A. Wulzer, JHEP 09, 060 (2014), eprint 1402.4431.
  * Flórez et al. (2017) A. Flórez, A. Gurrola, W. Johns, Y. D. Oh, P. Sheldon, D. Teague, and T. Weiler, Phys. Lett. B 767, 126 (2017), eprint 1609.09765.
  * Aad et al. (2023) G. Aad et al. (ATLAS), JHEP 06, 016 (2023), eprint 2207.00230.
  * Sirunyan et al. (2021) A. M. Sirunyan et al. (CMS), Eur. Phys. J. C 81, 688 (2021), eprint 2102.08198.
  * Aaboud et al. (2018) M. Aaboud et al. (ATLAS), JHEP 03, 042 (2018), eprint 1710.07235.
  * Tumasyan et al. (2022) A. Tumasyan et al. (CMS), Phys. Rev. D 105, 032008 (2022), eprint 2109.06055.
  * Aaboud et al. (2017) M. Aaboud et al. (ATLAS), JHEP 10, 182 (2017), eprint 1707.02424.
  * Sirunyan et al. (2018) A. M. Sirunyan et al. (CMS), JHEP 06, 120 (2018), eprint 1803.06292.
  * Cirelli et al. (2006) M. Cirelli, N. Fornengo, and A. Strumia, Nucl. Phys. B 753, 178 (2006), eprint hep-ph/0512090.
  * Aghanim et al. (2020b) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020b), eprint 1807.06209.
  * Aprile et al. (2018) E. Aprile et al. (XENON), Phys. Rev. Lett. 121, 111302 (2018), eprint 1805.12562.
  * Aprile et al. (2022) E. Aprile et al. (XENON) (2022), eprint 2210.07591.
  * Meng et al. (2021) Y. Meng et al. (PandaX-4T), Phys. Rev. Lett. 127, 261802 (2021), eprint 2107.13438.
  * Aalbers et al. (2022) J. Aalbers et al. (LZ) (2022), eprint 2207.03764.
  * Ellis et al. (2000) J. R. Ellis, A. Ferstl, and K. A. Olive, Phys.Lett. B481, 304 (2000), eprint hep-ph/0001005.
  * Bhattacharya et al. (2017b) S. Bhattacharya, P. Poulose, and P. Ghosh, JCAP 04, 043 (2017b), eprint 1607.08461.
  * Herrero-Garcia et al. (2017) J. Herrero-Garcia, A. Scaffidi, M. White, and A. G. Williams, JCAP 11, 021 (2017), eprint 1709.01945.
  * Bhattacharya et al. (2020) S. Bhattacharya, N. Chakrabarty, R. Roshan, and A. Sil, JCAP 04, 013 (2020), eprint 1910.00612.
  * Billard et al. (2022) J. Billard et al., Rept. Prog. Phys. 85, 056201 (2022), eprint 2104.07634.
  * O’Hare (2021) C. A. J. O’Hare, Phys. Rev. Lett. 127, 251802 (2021), eprint 2109.03116.
  * Gunion et al. (2000) J. F. Gunion, H. E. Haber, G. L. Kane, and S. Dawson, Front.Phys. 80, 1 (2000).
  * Kolb and Turner (1990) E. W. Kolb and M. S. Turner, _The Early Universe_ , vol. 69 (1990), ISBN 978-0-201-62674-2.
  * Srednicki et al. (1988) M. Srednicki, R. Watkins, and K. A. Olive, Nucl.Phys. B310, 693 (1988).
  * Edsjo and Gondolo (1997) J. Edsjo and P. Gondolo, Phys. Rev. D 56, 1879 (1997), eprint hep-ph/9704361.
  * Ghosh et al. (2022) P. Ghosh, S. Mahapatra, N. Narendra, and N. Sahu, Phys. Rev. D 106, 015001 (2022), eprint 2107.11951.