aainstitutetext:  Physik Department T31, James-Franck-Straße 1, Technische
Universität München, D-85748 Garching, Germanybbinstitutetext:  Theoretical
Physics Department, CERN, 1211 Geneva 23, Switzerlandccinstitutetext:  School
of Physics, Korea Institute for Advanced Study, Seoul 02455,
Koreaddinstitutetext:  College of Physics, Sichuan University, Chengdu 610065,
China

# A multi-component SIMP model with $U(1)_{X}\rightarrow Z_{2}\times Z_{3}$

Soo-Min Choi b Jinsu Kim c Pyungwon Ko d and Jinmian Li<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Multi-component dark matter scenarios are studied in the model with $U(1)_{X}$
dark gauge symmetry that is broken into its product subgroup $Z_{2}\times
Z_{3}$ á la Krauss-Wilczek mechanism. In this setup, there exist two types of
dark matter fields, $X$ and $Y$, distinguished by different $Z_{2}\times
Z_{3}$ charges. The real and imaginary parts of the $Z_{2}$-charged field,
$X_{R}$ and $X_{I}$, get different masses from the $U(1)_{X}$ symmetry
breaking. The field $Y$, which is another dark matter candidate due to the
unbroken $Z_{3}$ symmetry, belongs to the Strongly Interacting Massive
Particle (SIMP)-type dark matter. Both $X_{I}$ and $X_{R}$ may contribute to
$Y$’s $3\rightarrow 2$ annihilation processes, opening a new class of SIMP
models with a local dark gauge symmetry. Depending on the mass difference
between $X_{I}$ and $X_{R}$, we have either two-component or three-component
dark matter scenarios. In particular two- or three-component SIMP scenarios
can be realised not only for small mass difference between $X$ and $Y$, but
also for large mass hierarchy between them, which is a new and unique feature
of the present model. We consider both theoretical and experimental
constraints, and present four case studies of the multi-component dark matter
scenarios.

††preprint:  TUM-HEP-1319/21
CERN-TH-2021-030
KIAS-P21008

##  1 Introduction

For a few decades or so, the Weakly Interacting Massive Particle (WIMP)
paradigm has been one of the mainstreams of dark matter (DM) physics (see, for
example, Refs. Arcadi:2017kky ; Roszkowski:2017nbc  for a recent review).
Within the WIMP scenarios, DM mass could be of the electroweak scale and its
interactions with the Standard Model (SM) particles have the weak-interaction
strength, and thermal relic density could be realised by freeze-out mechanism.
One would have anticipated to observe DM particle by (in)direct detections and
at colliders. However the null results from the LHC searches for new particles
and direct detections of DM with mass in the electroweak scale lead us to
consider more seriously new paradigms beyond the WIMP paradigm.

The Strongly Interacting Massive Particle (SIMP) paradigm Hochberg:2014dra ;
Hochberg:2014kqa , where the DM relic density is dominantly determined by
$3\rightarrow 2$ processes in the dark sector, is one of such new paradigms
and gained a great attention during the last several years; see Refs.
Bernal:2015bla ; Lee:2015gsa ; Choi:2015bya ; Hochberg:2015vrg ; Choi:2016hid
; Choi:2016tkj ; Dey:2016qgf ; Choi:2017mkk ; Choi:2017zww ; Choi:2018iit ;
Hochberg:2018vdo ; Dey:2018yjt ; Choi:2019zeb ; Maity:2019vbo ; Choi:2020ara ;
Katz:2020ywn  for studies on various aspects of SIMP scenarios. DM in the SIMP
scenarios typically has a strong coupling and small mass scale (typically $10$
MeV–$1$ GeV) Hochberg:2014dra . Thus the SIMP DM scenario is well-suited for
explaining the small-scale problems of collision-less vanilla cold-DM paradigm
such as the core-cusp problem and too-big-to-fail problem; see e.g. Ref.
Tulin:2017ara .

As of now, parameter spaces for most of the single-component DM models are
being tightly constrained by various experiments. Single-component DM models
are simple for theoretical ideas on the model building, phenomenological
analyses, and comparison with various experimental data. However, there is no
compelling reason that the DM density of the current Universe is composed of
only one component. For example the visible sector of the Universe contributes
only $\sim 5\%$ to the energy density of the current Universe, and it has a
number of stable or long-lived particles; electron, proton, photon, and three
species of active light neutrinos. Presently the energy density of the dark
Universe in the form of matter is about five times more than the visible
Universe. Therefore it would be an oversimplification to assume that the dark
sector of the Universe consists of a single-component DM. Even if the dark
sector of the current Universe is dominated by a single-component DM, there
could be more DM species in the earlier Universe that would modify the
evolution of the early Universe. They may leave some footprints that can be
observed in the current Universe, e.g. on Cosmic Microwave Background or
matter power spectra.

One important feature that has to be taken into account seriously in DM
physics is that DM particle should be absolutely stable or its lifetime should
be much longer than the age of the Universe. Usually the DM stability is
supposed to be guaranteed by some dark symmetries such as $Z_{2}$. If this
dark symmetry is a global symmetry, it will be violated by gravity effects.
For example, let us consider a real scalar DM model with dark $Z_{2}$ symmetry
under which the real scalar DM field $S$ is $Z_{2}$-odd, namely transforming
like $S\rightarrow-S$, whereas all the SM fields are $Z_{2}$-even. Then the DM
$S$ will be stable at renormalizable level. However there will be
$Z_{2}$-breaking dimension-5 operators such as

$\displaystyle\frac{1}{\Lambda}S(H^{\dagger}H)^{2}\,,\quad\frac{1}{\Lambda}SF_{\mu\nu}F^{\mu\nu}\,,\quad{\it
etc.}$

Then $S$ will decay too fast to be a good cold-DM candidate even for
$\Lambda\sim M_{\rm Pl}$ where $M_{\rm Pl}$ is the Planck mass, unless its
mass is very light $m_{S}\lesssim\mathcal{O}(1)$ eV. The $Z_{2}$ breaking from
dimension-6 operators will not be harmful at all, since one can make DM
lifetime long enough for $\Lambda\ll M_{\rm Pl}$, like the proton lifetime in
the SM with $B$\- and $L$-violating dimension-6 operators. One can evade this
problem of dimension-5 operators that violate dark global symmetry by
promoting the dark global symmetry to a local dark gauge symmetry. In quantum
field theory, these could be achieved by assuming that DM carries a dark
charge that is either exactly conserved or approximately conserved. Finally
the symmetry associated with the exactly conserved charge is implemented into
a local gauge symmetry Baek:2013qwa  (see also Sec. IV B and C of Ref.
Ko:2014nha  for more discussion on this issue). Alternatively, the
approximately conserved dark charge may be an accidental symmetry of an
underlying dark gauge theory, like the $U(1)$ baryon number in the SM being an
accidental symmetry of the SM and being broken only by dimension-6 operators.
Thus, for either absolutely stable or long-lived DM, we end up with a local
gauge theory 111 One loophole in this argument is that the DM particle can be
long-lived enough if it has a very small mass so that its decay channels are
kinematically limited. Axion or light sterile neutrinos are good examples for
such a case. . This strategy was the cornerstone of the SM of particle
physics, and it could be extended to DM model building. In DM models with a
dark gauge symmetry, there are dark gauge bosons and dark Higgs bosons in
addition to the DM particles with various dark charges. Furthermore, there
could be more than one generation of DM particles, as in the SM. There could
be more than one stable or long-lived particles in the dark sector in general.

Therefore, from theoretical point of view, it may be more natural to have a
multi-component DM scenario once the symmetry that is responsible for the DM
stability gets enlarged, e.g. $Z_{N}$ with $N\geq 4$, a product of two
symmetry groups such as $Z_{2}\times Z_{2}^{\prime}$, topologically stable DM,
or dark gauge group with larger than rank-1, etc. Recent studies on the multi-
component DM scenarios within WIMP scenarios include Refs. Ko:2010at ;
Drozd:2011aa ; Aoki:2012ub ; Baek:2013dwa ; Ko:2014bka ; Bian:2014cja ;
Karam:2015jta ; Karam:2016rsz ; Aoki:2016glu ; Ko:2016fcd ;
Bhattacharya:2016ysw ; Ahmed:2017dbb ; Chakraborti:2018lso ; Poulin:2018kap ;
YaserAyazi:2018lrv ; Chakraborti:2018aae ; Aoki:2018gjf ; Bhattacharya:2019fgs
; Chen:2019pnt ; Yaguna:2019cvp .

One of the most interesting multi-component scenarios proposed so far would be
the dark QCD (DQCD) with a light Nambu-Goldstone boson DM (dark pion DM). One
can consider both WIMP Hur:2007uz ; Ko:2008ug ; Ko:2009zz ; Bai:2010qg ;
Ko:2010rj ; Hur:2011sv ; Bai:2013xga ; Hatanaka:2016rek  and SIMP scenarios
Hochberg:2014kqa ; Hochberg:2015vrg ; Choi:2018iit ; Katz:2020ywn  in DQCD
models . In the multi-component SIMP models proposed so far in the literature,
various DM components are related with each other. For example, dark pion DM
with different dark flavours in DQCD models are related with each other by a
dark flavour symmetry. They are not independent with each other when we
consider their origin, and their masses cannot be arbitrarily different. In
the case of the dark pion DM in DQCD, for example, one assumes that the dark
pions are pseudo-Nambu-Goldstone bosons. This picture requires that
$m_{q_{i}}+m_{q_{j}}\ll\Lambda_{\rm dc}$, where $\Lambda_{\rm dc}$ is the dark
confinement scale, which is nothing but the good old partially-conserved-
axialvector-current (PCAC) condition, in order that dark pion DMs are good
Nambu-Goldstone bosons. On the other hand, dark pion DM mass is given by
$m_{\pi}^{2}\sim m_{q}\Lambda_{\rm dc}$, which is bounded from above by
$\Lambda^{2}_{\rm dc}$. Therefore, one cannot imagine a wide separation in the
dark pion DM masses with different flavours. They are bounded by $\Delta
m_{\pi}^{2}\sim\Lambda^{2}_{\rm dc}$. In other words, we can say that they are
not genuine multi-component SIMP models where arbitrary masses are allowed,
including the large mass hierarchy. For genuine multi-component DM models,
different DM components should have different masses (and different lifetimes
in the case of a decaying DM), different (dark) charges, and carry even
different spins.

In this paper we introduce two independent fields $X$ and $Y$ with different
dark charges (and even different spins) from the beginning, and show that two-
or three-component SIMP scenarios can be realised. Therefore the model we
present in the following could be taken as a real multi-component SIMP
scenario. Our study is not based on the effective field theory approach. The
model constructed in this work is renormalizable and gauge invariant under the
SM and the dark gauge symmetry groups, and thus a UV-complete one. Our
conclusion can thus be taken as reliable in the entire energy range
considered.

At this point, it is worthwhile to remind ourselves that WIMP and SIMP are
about mechanisms to reproduce observed DM relic density in a given DM model,
and not about the underlying particle physics models for DM. Within the same
DM model, different mechanisms can be realised depending on the parameter
values and mass scales. Likewise, either single-component or multi-component
WIMP or SIMP scenarios can be realised within the same DM model, depending on
the parameter values and mass scales of various particles in the model. The
model discussed in the following shows how rich the dark sector can be,
depending on the parameter spaces (couplings and mass scales), and how DM
particles are stabilised ($U(1)$ dark charge assignments in our case). The
multi-component SIMP scenarios are realised only in a particular corner of the
parameter space. In a generic parameter space, we expect that single-component
SIMP scenarios or WIMP scenarios would be realised, but these two options
shall not be discussed in detail, since they are not qualitatively different
from the literature of the SIMP/WIMP scenarios, and thus are not the main
interest of this work.

In this work we consider a multi-component DM model with $U(1)_{X}$ dark gauge
symmetry broken into its product subgroup $Z_{2}\times Z_{3}$ á la Krauss-
Wilczek mechanism Krauss:1988zc . There exist two types of DM fields, $X$ and
$Y$, which are distinguished by different $U(1)_{X}$ charges with the
following symmetry breaking pattern: $U(1)_{X}\rightarrow Z_{2}\times Z_{3}$.
The stability of $X$ is guaranteed by the unbroken $Z_{2}$ symmetry, while the
unbroken $Z_{3}$ symmetry ensures the stability of $Y$. Then we can realise a
multi-component SIMP scenario through various DM number-changing scattering
processes such as $YYY\rightarrow XX$. Note that $X$ can be either scalar or
fermion, whereas $Y$ should be scalar. The relic densities of different DM
fields are determined through interactions between the dark fields. For the
$Z_{2}$-charged field $X$, number-changing $2\rightarrow 2$ processes are
responsible for $X$’s relic density. On the other hand, the relic density of
the $Z_{3}$-charged field $Y$ is given by number-changing $3\rightarrow 2$
processes. As the processes responsible for the relic densities of the DM
fields occur within the dark sector and both $X$ and $Y$ participate in the
number-changing processes, we dub our model a multi-component SIMP DM
scenario. The field $X$ receives different masses for its real and imaginary
parts from the vacuum expectation value (VEV) of the dark Higgs, breaking the
$U(1)_{X}$ symmetry. Due to the interactions between $X$ and $Y$, both $X_{I}$
and $X_{R}$ may contribute to $Y$’s $3\rightarrow 2$ processes. It opens a new
class of SIMP models. Depending on the mass gap between $X_{I}$ and $X_{R}$,
two-component or three-component DM scenarios can be realised. One of the most
distinctive features of this model is that the different DM fields $X$ and $Y$
may have vastly different masses, while significantly contributing to the
total DM relic density today. This is in sharp contrast to most of the studies
on multi-component SIMP where each component could coexist when the mass
difference is small enough Choi:2016hid ; Hochberg:2018vdo ; Katz:2020ywn .
For completeness we consider both cases with large mass hierarchy and small
mass difference.

It should be emphasised that the present work, which considers multi-component
SIMP scenarios with a large mass hierarchy, is mostly motivated by theoretical
reasons, as a kind of existence proof of such a case within a theoretically
well-defined DM model. Still there is an interesting phenomenological
consequence. The coupling used to maintain the kinetic equilibrium between the
DM sector and the SM sector may also be used to detect the boosted $Y$ (from
$X$ annihilation) in DM and neutrino detectors. The boosted DM signal is one
of the most important predictions of the multi-component DM scenarios with
large mass hierarchy.

The paper is organised as follows. We first describe our model setup in
Section 2, where we set the notations, explain the $U(1)_{X}$ symmetry
breaking, and identify DM candidates. In Section 3, we study phenomenology of
the multi-component DM scenarios. The evolution of the DM system is described
by four coupled Boltzmann equations, each for $Y$, $X_{I}$, $X_{R}$, and
$Z^{\prime}$. We numerically solve the Boltzmann equations and present two
benchmark cases for the two-component scenario and two benchmark cases for the
three-component scenario. In Section 3.2, we consider theoretical constraints,
such as perturbativity and unitarity, as well as experimental constraints,
including LHC constraints for the invisible Higgs decay, the DM direct
detection, and $Z^{\prime}$ searches. One of the crucial constraints for a
SIMP-type DM is the sufficient release of its kinetic energy. We address this
issue in detail in Section 3.2. Taking into account all of the constraints we
present our results in Sections 3.3 and 3.4. Finally we conclude in Section 4.

##  2 Model

We consider a dark $U(1)_{X}$ gauge symmetry with three different types of
scalars in the dark sector: two DM candidates $X$ and $Y$ plus the dark Higgs
field $\phi$. We assign the $U(1)$ charges of these fields as follows 222 We
can also consider another possibility with spin-1/2 dark fermion $\psi$ with
$q_{\psi}=1/2$, instead of dark scalar $X$. We shall not pursue this case in
detail in this paper, since the qualitative features of multi-component SIMP
scenarios would be the same.:

$\displaystyle(q_{X},q_{Y},q_{\phi})=(1/2,1/3,1)\,.$ (1)

Then, the gauge-invariant renormalizable Lagrangian of this model is given by

$\displaystyle\mathcal{L}$ $\displaystyle=\mathcal{L}_{\overline{{\rm
SM}}}-\frac{1}{4}\hat{X}_{\mu\nu}\hat{X}^{\mu\nu}-\frac{\sin\epsilon}{2}\hat{X}_{\mu\nu}\hat{B}^{\mu\nu}+|D_{\mu}H|^{2}+|D_{\mu}X|^{2}+|D_{\mu}Y|^{2}+|D_{\mu}\phi|^{2}-V\,,$
(2)

where $\mathcal{L}_{\overline{{\rm SM}}}$ is the SM Lagrangian without the
Higgs sector, $H$ is the SM Higgs field, and $\epsilon$ is the
$U(1)_{Y}$–$U(1)_{X}$ gauge kinetic mixing parameter. The covariant derivative
of the dark sector is defined as
$D_{\mu}\equiv\partial_{\mu}-iq_{i}g_{X}\hat{X}_{\mu}$, with $g_{X}$ being the
dark gauge coupling and the $U(1)_{X}$-charge assignments of dark fields
$q_{i}$ $(i=X,Y,\phi)$ are defined in Eq. (1).

The scalar potential is given by

$\displaystyle V$
$\displaystyle=m_{X}^{2}|X|^{2}+m_{Y}^{2}|Y|^{2}+\lambda_{\phi}\left(|\phi|^{2}-\frac{v_{\phi}^{2}}{2}\right)^{2}+\lambda\left(|H|^{2}-\frac{v^{2}}{2}\right)^{2}+\lambda_{X}|X|^{4}+\lambda_{Y}|Y|^{4}$
$\displaystyle\quad+\lambda_{XY}|X|^{2}|Y|^{2}+\lambda_{X\phi}|X|^{2}\left(|\phi|^{2}-\frac{v_{\phi}^{2}}{2}\right)+\lambda_{Y\phi}|Y|^{2}\left(|\phi|^{2}-\frac{v_{\phi}^{2}}{2}\right)$
$\displaystyle\quad+\lambda_{XH}|X|^{2}\left(|H|^{2}-\frac{v^{2}}{2}\right)+\lambda_{YH}|Y|^{2}\left(|H|^{2}-\frac{v^{2}}{2}\right)+\lambda_{\phi
H}\left(|\phi|^{2}-\frac{v_{\phi}^{2}}{2}\right)\left(|H|^{2}-\frac{v^{2}}{2}\right)$
$\displaystyle\quad+\left[\mu_{X\phi}X^{\dagger
2}\phi+\lambda^{\prime}_{Y\phi}Y^{\dagger 3}\phi+\text{h.c.}\right]\,,$ (3)

where $v_{\phi}$ and $v$ are VEVs of the dark Higgs and the SM Higgs fields.
After $U(1)_{X}$ symmetry breaking by a nonzero VEV
$\langle\phi\rangle=v_{\phi}/\sqrt{2}$, the last line of Eq. (2) still has
$Z_{2}\times Z_{3}$ discrete gauge symmetries 333 See also Refs. Baek:2014kna
; Baek:2020owl ; Kang:2021oes  and Refs. Ko:2014nha ; Ko:2014loa ; Guo:2015lxa
; Ko:2020gdg  for scalar DM models with $U(1)$ dark gauge symmetry broken into
$Z_{2}$ and $Z_{3}$, respectively.:

$\displaystyle X$ $\displaystyle\rightarrow e^{i\pi}X=-X\,,$ $\displaystyle Y$
$\displaystyle\rightarrow e^{\pm i2\pi/3}Y\,.$

We note that $\mu_{X\phi}$ has the mass dimension one. Expanding the dark
Higgs and the SM Higgs around their VEVs, $v_{\phi}$ and $v$, and writing
$X=(X_{R}+iX_{I})/\sqrt{2}$, the real and the imaginary parts of $X$ get
different masses from the $U(1)_{X}$ symmetry breaking. The mass splitting
between $X_{R}$ and $X_{I}$ is generated by the $\mu_{X\phi}$ term
Baek:2014kna :

$\displaystyle m_{X_{I}}^{2}$
$\displaystyle=m_{X}^{2}-\sqrt{2}\mu_{X\phi}v_{\phi}\,,$ (4) $\displaystyle
m_{X_{R}}^{2}$ $\displaystyle=m_{X}^{2}+\sqrt{2}\mu_{X\phi}v_{\phi}\,.$ (5)

Taking $\mu_{X\phi}>0$ without loss of generality, we see that $X_{I}$ is
always lighter than $X_{R}$.

The mass-squared matrix of the Higgs fields is given by

$\displaystyle
M^{2}=\left(\begin{array}[]{cc}m_{hh}^{2}&m_{hh^{\prime}}^{2}\\\
m_{hh^{\prime}}^{2}&m_{h^{\prime}h^{\prime}}^{2}\end{array}\right)\,,$ (8)

with 444 We denote the physical SM Higgs field and the physical dark Higgs
field by $h$ and $h^{\prime}$, respectively.

$\displaystyle m_{hh}^{2}=2\lambda v^{2}\,,\quad
m_{h^{\prime}h^{\prime}}^{2}=2\lambda_{\phi}v_{\phi}^{2}\,,\quad
m_{hh^{\prime}}^{2}=\lambda_{\phi H}vv_{\phi}\,.$ (9)

The mass eigenvalues are then obtained as follows:

$\displaystyle
m_{1,2}^{2}=\frac{1}{2}\left[(m_{hh}^{2}+m_{h^{\prime}h^{\prime}}^{2})\mp\sqrt{(m_{hh}^{2}-m_{h^{\prime}h^{\prime}}^{2})^{2}+4m_{hh^{\prime}}^{4}}\right]\,.$
(10)

The mass eigenstates, which we denote by $h_{1,2}$, are related to the
interaction eigenstates as

$\displaystyle\left(\begin{array}[]{c}h\\\
h^{\prime}\end{array}\right)=\left(\begin{array}[]{cc}\cos\alpha&-\sin\alpha\\\
\sin\alpha&\cos\alpha\end{array}\right)\left(\begin{array}[]{c}h_{2}\\\
h_{1}\end{array}\right)\,,$ (17)

with the mixing angle

$\displaystyle\alpha=\frac{1}{2}\tan^{-1}\left[\frac{2m_{hh^{\prime}}^{2}}{m_{hh}^{2}-m_{h^{\prime}h^{\prime}}^{2}}\right]\,.$
(18)

With the $U(1)_{Y}$–$U(1)_{X}$ gauge kinetic mixing, the mass matrix of the
gauge boson can be diagonalised as follows. The kinetic mixing can be removed
in terms of $\tilde{B}_{\mu}$ and $\tilde{X}_{\mu}$ defined as Babu:1997st ,

$\displaystyle\left(\begin{array}[]{c}\hat{B}_{\mu}\\\
\hat{X}_{\mu}\end{array}\right)=\left(\begin{array}[]{cc}1&-\tan\epsilon\\\
0&1/\cos\epsilon\end{array}\right)\left(\begin{array}[]{c}\tilde{B}_{\mu}\\\
\tilde{X}_{\mu}\end{array}\right)\,,\qquad\hat{W}_{\mu}=\tilde{W}_{\mu}\,,$
(25)

where we defined $\tilde{W}_{\mu}$ for notational consistency. We then
diagonalise the mass matrix for $\tilde{B}_{\mu}$, $\tilde{X}_{\mu}$, and
$\tilde{W}^{3}_{\mu}$ as

$\displaystyle\left(\begin{array}[]{c}\tilde{B}_{\mu}\\\
\tilde{W}^{3}_{\mu}\\\
\tilde{X}_{\mu}\end{array}\right)=\left(\begin{array}[]{ccc}c_{w}&-s_{w}c_{\zeta}&s_{w}s_{\zeta}\\\
s_{w}&c_{w}c_{\zeta}&-c_{w}s_{\zeta}\\\
0&s_{\zeta}&c_{\zeta}\end{array}\right)\left(\begin{array}[]{c}A_{\mu}\\\
Z_{\mu}\\\ Z^{\prime}_{\mu}\end{array}\right)\,,$ (35)

where $c_{w}\equiv\cos\theta_{w}$, $s_{w}\equiv\sin\theta_{w}$,
$c_{\zeta}\equiv\cos\zeta$, $s_{\zeta}\equiv\sin\zeta$, and $\theta_{w}$ is
the Weinberg angle. Here $\zeta$ is defined as

$\displaystyle\tan 2\zeta\equiv-\frac{m_{\hat{Z}}^{2}s_{w}\sin
2\epsilon}{m_{\hat{X}}^{2}-m_{\hat{Z}}^{2}(c_{\epsilon}^{2}-s_{\epsilon}^{2}s_{w}^{2})}\,,$
(36)

where $m_{\hat{Z}}=\sqrt{(g_{1}^{2}+g_{2}^{2})}v/2$,
$m_{\hat{X}}=g_{X}v_{\phi}$, $c_{\epsilon}\equiv\cos\epsilon$, and
$s_{\epsilon}\equiv\sin\epsilon$. Subsequently, we obtain the following
relations between the original fields and the mass eigenstates:

$\displaystyle\hat{B}_{\mu}$ $\displaystyle=$ $\displaystyle
c_{w}A_{\mu}-(t_{\epsilon}s_{\zeta}+s_{w}c_{\zeta})Z_{\mu}+(s_{w}s_{\zeta}-t_{\epsilon}c_{\zeta})Z_{\mu}^{\prime}\,,$
$\displaystyle\hat{X}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{s_{\zeta}}{c_{\epsilon}}Z_{\mu}+\frac{c_{\zeta}}{c_{\epsilon}}Z_{\mu}^{\prime}\,,$
(37) $\displaystyle\hat{W}^{3}_{\mu}$ $\displaystyle=$ $\displaystyle
s_{w}A_{\mu}+c_{w}c_{\zeta}Z_{\mu}-c_{w}s_{\zeta}Z_{\mu}^{\prime}\,,$

with $t_{\epsilon}\equiv\tan\epsilon$ (similarly for $\zeta$ below). The
masses of the $Z$ and $Z^{\prime}$ gauge bosons are given by

$\displaystyle m_{Z}^{2}=m_{\hat{Z}}^{2}(1+s_{w}t_{\zeta}t_{\epsilon})\,,\quad
m_{Z^{\prime}}^{2}=\frac{m_{\hat{X}}^{2}}{c_{\epsilon}^{2}(1+s_{w}t_{\zeta}t_{\epsilon})}\,.$
(38)

In the following we assume small mixing angles, i.e. $\alpha\ll 1$ and
$\epsilon\ll 1$, hence $h_{2}\approx h$ and $h_{1}\approx h^{\prime}$, and
$m_{Z}^{2}\approx m_{\hat{Z}}^{2}$ and $m_{Z^{\prime}}^{2}\approx
m_{\hat{X}}^{2}$.

In our model, due to the unbroken $Z_{2}\times Z_{3}$ symmetry, fields with
$U(1)_{X}$ charges of 1/2 and 1/3 will be stable, becoming DM candidates. The
$X_{I}$ field with $q_{X}=1/2$ is absolutely stable and makes one of the DM
components. The heavier one $X_{R}$ decays into $X_{I},Z^{\prime(*)}$, and
could be another good DM candidate, if its lifetime is long enough compared
with the age of the Universe. When the mass splitting is so small,
$m_{X_{R}}-m_{X_{I}}\ll{\rm min}\\{2m_{Y},m_{Z^{\prime}}\\}$, that the $X_{R}$
decay channels such as $X_{R}\rightarrow X_{I},Z^{\prime}$ and
$X_{R}\rightarrow X_{I},Y,Y^{*}$ are closed, the $X_{R}$ field again becomes a
good DM candidate (see Appendix A for the three-body decay of $X_{R}$). The
field $Y$ with $q_{Y}=1/3$ is another DM candidate because of the unbroken
$Z_{3}$ symmetry. Then for example

$\displaystyle Y,Y,Y\rightarrow X_{I},X_{I}\,,\quad{\it etc.}$

could be possible through the dark Higgs and/or dark gauge boson exchanges,
which is a new class of SIMP models based on spontaneous breaking of dark
$U(1)_{X}\rightarrow Z_{2}\times Z_{3}$.

We are primarily interested in the setup where the field $Y$ is the lightest
among the DM candidates. To implement the SIMP scenario for $Y$, we assume
that the dark Higgs and the dark gauge boson are heavier than $Y$. Otherwise,
the DM will be pair-annihilated into dark Higgs/dark gauge boson, which is not
of interest in this work. Similarly, if either the mixing between the SM Higgs
boson and dark Higgs or the kinetic mixing between gauge bosons of $U(1)_{Y}$
and $U(1)_{X}$ is sizeable, the DM pair-annihilating into SM particles will be
important once kinematically open. We shall therefore focus on the parameter
spaces that have small mixing between SM Higgs and dark Higgs, as well as
small kinetic mixing in the gauge sector. We note that this setup is also
useful to suppress the DM-nucleon scattering cross section, which is tightly
constrained by the DM direct detection experiments. On the other hand, it is
necessary for $Y$ to lose its kinetic energy through elastic scattering with
SM particles Hochberg:2014dra ; Choi:2015bya . For this, a somewhat light
$Z^{\prime}$ is needed for the SIMP–SM kinetic equilibrium condition as we
discuss in Section 3.2. We focus on the case where $m_{Z^{\prime}}>2m_{Y}$;
thus $Z^{\prime}$ is not stable and will decay. As long as the decay is fast
enough, $Z^{\prime}$ cannot be a DM candidate. For the dark Higgs, we assume
that it is much heavier than $Y$, $X_{I,R}$, and $Z^{\prime}$ such that the
dark Higgs does not play any significant role in the dynamics. The cross
section for the process $X,X\rightarrow Y,Y^{*}$ features a four-point vertex
555 We note that contact interactions between $Z_{2}$-charged field and
$Z_{3}$-charged field are absent at renormalizable level in the fermionic
$Z_{2}$ DM, where the scalar field $X$ is replaced by a spin-1/2 dark fermion
$\psi$.  which, in the large-$m_{h^{\prime}}$ limit, is proportional to
$\lambda_{XY}^{2}$. The relic density of the $X$ field will be diluted through
this process if $\lambda_{XY}$ is large. When the dark Higgs is somewhat
lighter, couplings $\lambda_{X\phi}$ and $\lambda_{Y\phi}$ need to be set
appropriately such that the $X,X\rightarrow Y,Y^{*}$ is suppressed via the
destructive interferences between the contact process and dark Higgs mediated
process. We first assume that $\lambda_{XY}$ is small enough such that this
contact $X$-annihilation process is not efficient in the
large-$m_{h^{\prime}}$ limit. We then comment on the case of destructive
interferences, where $\lambda_{XY}$ may take a larger value, presenting
parameter sets with a lighter dark Higgs and larger dark Higgs-portal
couplings.

In summary the DM candidates are $X_{R}$, $X_{I}$, and $Y$ in our setup. Our
interest is to investigate whether two- or three-component DM scenarios can be
realised for the following four cases: i) $m_{X_{R}}\gg m_{X_{I}}\gg m_{Y}$,
ii) $m_{X_{R}}\gg m_{X_{I}}\sim m_{Y}$, iii) $m_{X_{R}}\sim m_{X_{I}}\gg
m_{Y}$, and iv) $m_{X_{R}}\sim m_{X_{I}}\sim m_{Y}$. The cases i) and ii) will
be the two-component DM scenarios, while iii) and iv) correspond to the three-
component scenarios. We do not consider the inverted mass hierarchy, i.e.
$m_{X_{I}}\lesssim m_{X_{R}}\ll m_{Y}$. On one hand, in order for the
$Z_{3}$-charged field $Y$ to serve as a SIMP DM, the mass should be less than
$\sim$50 MeV, in the large-$m_{h^{\prime}}$ limit, as number-changing
$3\rightarrow 2$ processes are dominant. On the other hand, a DM lighter than
$\sim$5 MeV is heavily constrained by cosmology such as Big Bang
Nucleosynthesis and Cosmic Microwave Background, as explored in Ref.
Sabti:2019mhn . Therefore, we do not consider the inverted mass hierarchy in
this work, and focus on the aforementioned cases i)–iv).

##  3 Multi-component Dark Matter Phenomenology

###  3.1 Generalities

Under the assumptions described in the previous section, we have four coupled
Boltzmann equations for describing the evolution of the system consisting of
$\\{Y,X_{I},X_{R},Z^{\prime}\\}$ 666 Since the dark Higgs is assumed to be
sufficiently heavier than the rest of the dark sector fields, it does not play
any significant role in the dynamics. . Schematically the Boltzmann equation
for the species $i$ is given by

$\displaystyle\dot{n}_{i}+3Hn_{i}$
$\displaystyle=-\langle\Gamma\rangle_{i\rightarrow j,k}\left(n_{i}-n_{i}^{\rm
eq}\frac{n_{j}n_{k}}{n_{j}^{\rm eq}n_{k}^{\rm
eq}}\right)-\langle\Gamma\rangle_{i\rightarrow j,k,l}\left(n_{i}-n_{i}^{\rm
eq}\frac{n_{j}n_{k}n_{l}}{n_{j}^{\rm eq}n_{k}^{\rm eq}n_{l}^{\rm eq}}\right)$
$\displaystyle\quad-\langle\sigma v\rangle_{i,j\rightarrow
l,m}\left(n_{i}n_{j}-n_{i}^{\rm eq}n_{j}^{\rm eq}\frac{n_{l}n_{m}}{n_{l}^{\rm
eq}n_{m}^{\rm eq}}\right)$ $\displaystyle\quad-\langle\sigma
v^{2}\rangle_{i,j,k\rightarrow l,m}\left(n_{i}n_{j}n_{k}-n_{i}^{\rm
eq}n_{j}^{\rm eq}n_{k}^{\rm eq}\frac{n_{l}n_{m}}{n_{l}^{\rm eq}n_{m}^{\rm
eq}}\right)\,,$ (39)

where $n_{i}$ ($n_{i}^{\rm eq}$) is the (equilibrium) number density of the
species $i$, $H$ the Hubble parameter, and $\dot{}\equiv d/dt$ with $t$ being
the cosmic time.

Here $\langle\sigma v\rangle$ and $\langle\sigma v^{2}\rangle$ are the
thermally-averaged cross sections, and $\langle\Gamma\rangle$ is the
thermally-averaged decay rate. The thermally-averaged cross section
$\langle\sigma v\rangle$ for the $i,j\rightarrow l,m$ process is given by

$\displaystyle\langle\sigma v\rangle=\frac{1}{S}\frac{1}{F}\sum_{\rm
spins}\frac{1}{64\pi^{2}E_{i}E_{j}(E_{i}+E_{j})}\int
d\Omega|\mathbf{p}_{f}||\mathcal{M}|^{2}\,,$ (40)

where $F$ and $S$ are the spin and symmetry factors, $E_{i,j}$ the energy of
incoming particle, and $\mathbf{p}_{f}$ the 3-momentum of final-state
particle. Similarly, the thermally-averaged cross section $\langle\sigma
v^{2}\rangle$ for the $i,j,k\rightarrow l,m$ process is given by

$\displaystyle\langle\sigma v^{2}\rangle=\frac{1}{S}\frac{1}{F}\sum_{\rm
spins}\frac{1}{128\pi^{2}E_{i}E_{j}E_{k}(E_{i}+E_{j}+E_{k})}\int
d\Omega|\mathbf{p}_{f}||\mathcal{M}|^{2}\,.$ (41)

For the decay rate, see Appendix A. We utilise FeynRules Alloul:2013bka  and
CalcHEP Belyaev:2012qa  to compute the scattering amplitudes. We work in the
nonrelativistic framework and take only $s$-wave contributions.

Instead of $n_{i}$ and $t$, it is convenient to use the yield $Y_{i}$ and the
dimensionless time variable $x$, defined as

$\displaystyle Y_{i}\equiv n_{i}/s(T)\,,\qquad x\equiv m_{i}/T\,,$ (42)

where $s(T)=(2\pi^{2}/45)g_{*s}T^{3}$ is the entropy density with $T$ being
the temperature of the heat bath. Note that, during the radiation-dominated
era,

$\displaystyle dt=\frac{x}{H(m_{i})}dx\,,$ (43)

where $H(m_{i})=\sqrt{(\pi^{2}/90)g_{*}}m_{i}^{2}/M_{\rm P}$. Note also that
$s(T)=s(m_{i})/x^{3}$. Denoting $s(m_{i})\equiv s$ and $H(m_{i})\equiv H$, the
Boltzmann equation can be recast into

$\displaystyle\frac{dY_{i}}{dx}$
$\displaystyle=-\frac{x}{H}\langle\Gamma\rangle_{i\rightarrow
j,k}\left(Y_{i}-Y_{i}^{\rm eq}\frac{Y_{j}Y_{k}}{Y_{j}^{\rm eq}Y_{k}^{\rm
eq}}\right)-\frac{x}{H}\langle\Gamma\rangle_{i\rightarrow
j,k,l}\left(Y_{i}-Y_{i}^{\rm eq}\frac{Y_{j}Y_{k}Y_{l}}{Y_{j}^{\rm
eq}Y_{k}^{\rm eq}Y_{l}^{\rm eq}}\right)$
$\displaystyle\quad-\frac{s}{Hx^{2}}\langle\sigma v\rangle_{i,j\rightarrow
l,m}\left(Y_{i}Y_{j}-Y_{i}^{\rm eq}Y_{j}^{\rm eq}\frac{Y_{l}Y_{m}}{Y_{l}^{\rm
eq}Y_{m}^{\rm eq}}\right)$
$\displaystyle\quad-\frac{s^{2}}{Hx^{5}}\langle\sigma
v^{2}\rangle_{i,j,k\rightarrow l,m}\left(Y_{i}Y_{j}Y_{k}-Y_{i}^{\rm
eq}Y_{j}^{\rm eq}Y_{k}^{\rm eq}\frac{Y_{l}Y_{m}}{Y_{l}^{\rm eq}Y_{m}^{\rm
eq}}\right)\,.$ (44)

The full Boltzmann equations are summarised in Appendix B. We choose
$m_{i}=m_{Y}$.

The DM relic density is then given by

$\displaystyle\Omega_{\rm
DM}h^{2}=\Omega_{y}h^{2}+\Omega_{X_{I}}h^{2}+\Omega_{X_{R}}h^{2}=\frac{s_{0}h^{2}}{\rho_{c}}\left(m_{Y}Y_{y}^{0}+m_{X_{I}}Y_{X_{I}}^{0}+m_{X_{R}}Y_{X_{R}}^{0}\right)\,,$
(45)

where the subscript “0” indicates the present value, $\rho_{c}$ is the
critical energy density, and $h\approx 0.67$ is the scaling factor for the
Hubble parameter. Here $Y_{y}\equiv 2Y_{Y}=2Y_{Y^{*}}$.

The input parameters for the model under consideration are chosen as

$\displaystyle m_{X_{R}}\,,\quad m_{X_{I}}\,,\quad m_{Y}\,,\quad
m_{Z^{\prime}}\,,\quad m_{h^{\prime}}\,,$ $\displaystyle
g_{X}\,,\quad\lambda_{Y}\,,\quad\lambda_{X}\,,\quad\lambda_{XY}\,,\quad\lambda_{Y\phi}\,,\quad\lambda_{X\phi}\,,\quad\lambda^{\prime}_{Y\phi}\,,$
$\displaystyle\epsilon\,,\quad\alpha\,.$

In the following, we fix $m_{Z^{\prime}}=200$ MeV, $m_{h^{\prime}}=30$ GeV,
$\lambda_{X}=0.025$, and $\epsilon=2\times 10^{-4}$. We assume that all the
quartic and portal couplings are non-negative. Furthermore we assume that
$\lambda^{\prime}_{Y\phi}$ is positive. As explained in the previous section,
$\lambda_{XY}$ is required to be small in order to implement multi-component
DM scenarios in the large-$m_{h^{\prime}}$ limit. We first perform the
analysis with $\lambda_{XY}=0$ and give a maximum value of $\lambda_{XY}$ that
do not jeopardise the multi-component scenarios realised by our chosen
benchmark points. We then present another sets of parameter values with a
somewhat smaller $m_{h^{\prime}}$ together with large dark Higgs-portal
couplings, $\lambda_{X\phi}$ and $\lambda_{Y\phi}$. In this case the
destructive interference occurs and $\lambda_{XY}$ may take a larger value.
The mixing angle between the SM Higgs and the dark Higgs is chosen to be
small, $\alpha=10^{-2}$, to evade the invisible Higgs decay constraint.
Consequently, couplings to the SM Higgs field are not important. Before
presenting our results for the multi-component DM scenarios, let us consider
theoretical and experimental constraints.

###  3.2 Constraints

We take into account the following theoretical and experimental constraints.

  * •

Perturbativity
The perturbativity condition imposes the following conditions for the coupling
parameters:

$\displaystyle\lambda_{i}<4\pi\,,\quad\lambda^{\prime}_{Y\phi}<4\pi\,,\quad
g_{X}<4\pi\,,$ (46)

where $i=\\{X,Y,XY,X\phi,Y\phi,XH,YH,\phi H\\}$.

  * •

Unitarity
The unitarity bounds are given by

$\displaystyle|\mathcal{M}|_{i,j\rightarrow l,m}<8\pi\,,$ (47)

for $i,j,l,m=\\{X_{I},X_{R},Y\\}$. The unitarity bound is shown in Fig. 2 in
the ($\Omega_{y}h^{2},g_{X}$) plane.

  * •

Kinetic equilibrium
If the $3\rightarrow 2$ processes are dominant, the kinetic energies of the DM
candidates keep piling up. Thus, the kinetic energies need to be released.
This can be achieved when the DM particles have a way to talk to the SM
thermal bath, so that elastic scattering processes release the kinetic energy.
See Refs. Choi:2019zeb ; Gondolo:2012vh ; Kuflik:2015isi  for details.
Following Ref. Gondolo:2012vh , we compute the momentum relaxation rate
$\gamma$ of a DM particle scattered off the relativistic SM thermal bath
species, namely the electron. We then require the relaxation rate to be larger
than or compatible with the kinetic energy production from the $3\rightarrow
2$ processes, i.e. $\gamma\gtrsim Hx_{\rm KD}^{2}$, until the kinetic
decoupling time $x_{\rm KD}\equiv m_{Y}/T_{\rm KD}$. We set the kinetic
decoupling time to be similar to the freeze-out time, $x_{f}$. In our study,
we use the $Z^{\prime}$-portal interaction to communicate with the
relativistic SM heat bath. Since the relativistic SM degree of freedom is
electron at $T\sim O({\rm MeV})$, we consider the DM–$e$ elastic scattering
process. The most dominant $3\rightarrow 2$ processes are the $Y$-involving
processes, $Y,Y,Y\rightarrow Y,Y^{*}$ and $Y,Y^{*},Y^{*}\rightarrow Y,Y$. Thus
we focus on the kinetic energy release of the field $Y$. From the kinetic
equilibrium condition we see that the $Z^{\prime}$ mass is bounded as follows:

$\displaystyle
m_{Z^{\prime}}\lesssim\left(\frac{31\pi^{3}e^{2}c_{w}^{2}g_{X}^{2}}{1701x_{\rm
KD}^{6}}\right)^{1/4}\left(m_{Y}^{3}M_{\rm P}\right)^{1/4}\sqrt{\epsilon}\,.$
(48)

The kinetic decoupling region is sketched in Fig. 1.

  * •

DM relic density
The current DM relic density is given by Aghanim:2018eyx

$\displaystyle\Omega_{\rm DM}h^{2}=0.1200\pm 0.0012\,.$ (49)

  * •

Invisible Higgs decay
The Higgs invisible decay is bounded Sirunyan:2018owy  as

$\displaystyle{\rm BR}_{h}<0.19\,,$ (50)

at 90% C.L., where ${\rm BR}_{h}$ is branching ratio of the Higgs invisible
decay. We choose a small enough mixing angle $\alpha=10^{-2}$ to satisfy the
bound.

  * •

Direct detection
For the DM direct detection bounds, we consider Xenon10 Essig:2017kqs ,
Xenon1T Aprile:2019jmx ; Aprile:2019xxb , and expected SENSEI-100 bound
Battaglieri:2017aum . For all the benchmark scenarios that we consider (Table
1), the DM-electron cross sections are tiny, and thus safe from the direct
detection constraints. The DM direct detection bounds are depicted in Fig. 2.
The DM-quark scattering bounds are unimportant since the DM masses are less
than 1 GeV in our benchmark scenarios.

  * •

$Z^{\prime}$ searches
We use the BaBar visible/invisible decay of $Z^{\prime}$ constraints
Lees:2014xha ; Lees:2017lec , Belle2 invisible decay of $Z^{\prime}$
constraint Essig:2013vha , and beam dump constraints coming from NA64
Banerjee:2019hmi ; Gninenko:2019qiv , SHiP SHiP:2020noy , E137 Bjorken:1988as
, and Orsay Davier:1989wz ; see also Refs. Andreas:2012mt ; Marsicano:2018krp
. These constraints impose bounds on the kinetic mixing angle $\epsilon$,
depending on the $Z^{\prime}$ mass, and are depicted in Fig. 1.

  * •

DM self-scattering cross section
The SIMP-type DM models naturally have a sizeable self-scattering cross
section that helps to solve the small-scale problems mentioned in Section 1.
We require Spergel:1999mh ; Elbert:2014bma ; Tulin:2017ara ; Chu:2018fzy

$\displaystyle 0.1\;{\rm cm^{2}/g}<\frac{\sigma_{\rm self}}{m_{\rm
DM}}<10\;{\rm cm^{2}/g}\,.$ (51)

The Bullet Cluster Markevitch:2003at ; Clowe:2003tk  imposes a stricter
constraint, $\sigma_{\rm self}/m_{\rm DM}<1\;{\rm cm^{2}/g}$. See also Ref.
Randall:2007ph . A similar bound is obtained from cosmological simulations
with self-interacting DM on the scales of galaxies and galaxy clusters
Rocha:2012jg ; Peter:2012jh .

For the $m_{\rm DM}$, we take an effective DM mass. For the cases i) and iii)
the dominant contribution to $\sigma_{\rm self}$ comes from the $Y$ self-
scattering cross section, as $X_{I}$ and $X_{R}$ are much heavier than $Y$. We
thus choose $m_{\rm DM}=m_{Y}$. For the case of ii) both $X_{I}$ and $Y$ may
contribute; thus we choose $m_{\rm DM}=(m_{Y}+m_{X_{I}})/2$. The case iv) is
more complicated. In this case $X_{R}$, $X_{I}$, and $Y$ all contribute. We
thus consider all the contributions with the effective masses. All the
contributions are weighted by the fractions of the relic density. In summary,

$\displaystyle\frac{\sigma_{\rm self}}{m_{\rm DM}}\bigg{|}_{{\it i)},{\it
iii)}}$ $\displaystyle=\;\;\frac{1}{4}\left(\frac{\Omega_{y}}{\Omega_{\rm
DM}}\right)^{2}\frac{1}{m_{Y}}\left(2\sigma_{Y,Y\rightarrow
Y,Y}+\sigma_{Y,Y^{*}\rightarrow Y,Y^{*}}\right)\,,$ (52)
$\displaystyle\frac{\sigma_{\rm self}}{m_{\rm DM}}\bigg{|}_{{\it ii)}}$
$\displaystyle=\;\;\frac{\sigma_{\rm self}}{m_{\rm DM}}\bigg{|}_{{\it i)},{\it
iii)}}+\left(\frac{\Omega_{X_{I}}\Omega_{y}}{\Omega_{\rm
DM}^{2}}\right)\frac{2}{m_{Y}+m_{X_{I}}}\sigma_{X_{I},Y\rightarrow X_{I},Y}$
$\displaystyle\quad+\left(\frac{\Omega_{X_{I}}}{\Omega_{\rm
DM}}\right)^{2}\frac{1}{m_{X_{I}}}\left(\sigma_{X_{I},X_{I}\rightarrow
X_{I},X_{I}}+\sigma_{X_{I},X_{I}\rightarrow Y,Y^{*}}\right)\,,$ (53)
$\displaystyle\frac{\sigma_{\rm self}}{m_{\rm DM}}\bigg{|}_{{\it iv)}}$
$\displaystyle=\;\;\frac{\sigma_{\rm self}}{m_{\rm DM}}\bigg{|}_{{\it
ii)}}+\left(\frac{\Omega_{X_{R}}\Omega_{Y}}{\Omega_{\rm
DM}^{2}}\right)\frac{2}{m_{X_{R}}+m_{Y}}\left(\sigma_{X_{R},Y\rightarrow
X_{R},Y}+\sigma_{X_{R},Y\rightarrow X_{I},Y}\right)$
$\displaystyle\quad+\left(\frac{\Omega_{X_{R}}\Omega_{X_{I}}}{\Omega_{\rm
DM}^{2}}\right)\frac{2}{m_{X_{R}}+m_{X_{I}}}\left(\sigma_{X_{I},X_{R}\rightarrow
X_{I},X_{R}}+\sigma_{X_{I},X_{R}\rightarrow Y,Y^{*}}\right)$
$\displaystyle\quad+\left(\frac{\Omega_{X_{R}}}{\Omega_{\rm
DM}}\right)^{2}\frac{1}{m_{X_{R}}}\left(\sigma_{X_{R},X_{R}\rightarrow
X_{R},X_{R}}+\sigma_{X_{R},X_{R}\rightarrow
X_{I},X_{I}}+\sigma_{X_{R},X_{R}\rightarrow Y,Y^{*}}\right)\,.$ (54)

In Fig. 2, $\sigma_{\rm self}/m_{\rm DM}=0.1$ ${\rm cm^{2}/g}$, $1$ ${\rm
cm^{2}/g}$, and $\sigma_{\rm self}/m_{\rm DM}>10$ ${\rm cm^{2}/g}$ are shown.

Taking all of the listed constraints into account we now present four
benchmark results for the multi-component DM scenarios.

Figure 1: Kinetic decoupling and $Z^{\prime}$-search constraints for the case
i) (left) and cases ii), iii), and iv) (right) in the ($m_{Z^{\prime}}$,
$\epsilon$) plane. See the main text in Section 3.2 for details of the used
constraints. We have the same plot for ii), iii), and iv) since the same
$m_{Y}$, $m_{Z^{\prime}}$, and $\epsilon$ are used. Furthermore, we do not see
a clearly visible difference between the left and right plots since the
difference in $m_{Y}$ is small. The blue points represent our benchmark cases
shown in Table 1.

Figure 2: For each case, the unitarity bound (grey), direct detection bound
(orange), and large self-scattering cross section bound $\sigma_{\rm
self}/m_{\rm DM}>10$ ${\rm cm^{2}/g}$ (red) are shown for different values of
the dark gauge coupling $g_{X}$ and the fraction of the $Y$ relic density
$\Omega_{y}h^{2}/0.12$. The black lines represent DM–electron scattering cross
section values, $10^{-40}$ ${\rm cm^{2}}$ (dashed), $10^{-41}$ ${\rm cm^{2}}$
(dot-dashed), and $10^{-42}$ ${\rm cm^{2}}$ (dotted) from top to bottom,
respectively. The red lines correspond to two different values of the self-
scattering cross section, $0.1$ ${\rm cm^{2}/g}$ (dotted) and $1$ ${\rm
cm^{2}/g}$ (dashed). For details of the used constraints, see the main text in
Section 3.2. The blue points represent our benchmark cases shown in Table 1.
Here we vary only $g_{X}$ and the fraction of the $Y$ relic density, with the
condition $(\Omega_{y}+\Omega_{X_{I}}+\Omega_{X_{R}})h^{2}=0.12$. For the two-
component DM cases i) and ii), $\Omega_{X_{R}}=0$ is chosen, while for the
three-component DM cases iii) and iv), $\Omega_{X_{R}}=\Omega_{X_{I}}$ is
assumed, except for the self-scattering cross section for which we follow the
method outlined in Section 3.2. The other input parameters are fixed (see
Table 1). We note that the correct present-day DM relic density constraint is
not strictly imposed, except at the benchmark points. Thus, it is important to
note that not all the white region is allowed. On the other hand, this
approach allows one to comprehensively understand the various constraints and
where our benchmark cases lie. The full picture requires an intensive
parameter scan by numerically solving the coupled Boltzmann equations which is
beyond the scope of the current work.

###  3.3 Two-component Scenario

When the mass of $X_{R}$ is sufficiently larger than
$m_{X_{I}}+m_{Z^{\prime}}$, $X_{R}$ cannot be a DM candidate as it decays into
$X_{I}$ and $Z^{\prime}$. Thus we have the two-component scenario. We consider
two cases: i) $m_{X_{R}}\gg m_{X_{I}}\gg m_{Y}$ and ii) $m_{X_{R}}\gg
m_{X_{I}}\approx m_{Y}$. The input parameter values are shown in Table 1.

Case |  $m_{X_{R}}$ [MeV] |  $m_{X_{I}}$ [MeV] |  $m_{Y}$ [MeV] | $g_{X}$ | $\lambda_{Y}$ | $\lambda_{Y\phi}$ | $\lambda_{X\phi}$ | $\lambda^{\prime}_{Y\phi}$
---|---|---|---|---|---|---|---|---
i) | 800 | 200 | 50 | 0.85 | 6.41 | 0.9 | 0.9 | 0.255
ii) | 400 | 40 | 37.5 | 1.3 | 6.27 | 1.3 | 0.5 | 0.295
iii) | 150.01 | 150 | 37.5 | 1.3 | 5.54 | 0.05 | 0.05 | 0.271
iv) | 40.001 | 40 | 37.5 | 1.3 | 6.27 | 2.65 | 2.2 | 0.295
Table 1: Input parameter values for the four benchmark cases. Cases i) and ii)
correspond to the two-component scenarios and cases iii) and iv) correspond to
the three-component scenarios. We chose $m_{Z^{\prime}}=200$ MeV,
$m_{h^{\prime}}=30$ GeV, $\lambda_{X}=0.025$, $\alpha=10^{-2}$, and
$\epsilon=2\times 10^{-4}$. Our four benchmark cases are marked as blue points
in Figs. 1 and 2.

We solve the Boltzmann equations given in Appendix B for the given parameters,
and the solutions are shown in Fig. 3. We see that, in the case of
$m_{X_{I}}\gg m_{Y}$, $X_{I}$ freezes out first at $x\approx 5$, followed by
the $Y$ freeze-out at $x\approx 20$. When $m_{X_{I}}\approx m_{Y}$, both
$X_{I}$ and $Y$ freeze out at almost the same time $x\approx 20$. In the case
i), we find that $X_{I}$ ($Y$) constitutes about 23.7% (76.3%) of the total
relic density, while in the case ii), the fractions of $X_{I}$ and $Y$ relics
are respectively 19.7% and 80.3%. Thus we identify these scenarios as two-
component DM scenarios. The self-scattering cross sections are given by $1.49$
${\rm cm^{2}/g}$ and $3.83$ ${\rm cm^{2}/g}$ for cases ${\it i)}$ and ${\it
ii)}$, respectively. We summarise the results in Table 2.

###  3.4 Three-component Scenario

When the mass gap between $X_{R}$ and $X_{I}$ is small enough, $X_{R}$ may
become a good DM candidate. We take the mass splitting to be small enough,
$m_{X_{R}}-m_{X_{I}}\ll{\rm min}\\{2m_{Y},m_{Z^{\prime}}\\}$, such that the
$X_{R}$ decay channels such as $X_{R}\rightarrow X_{I},Z^{\prime}$ and
$X_{R}\rightarrow X_{I},Y,Y^{*}$ are kinematically closed 777 Note that the
decay channel $X_{R}\rightarrow X_{I},\ell,\bar{\ell}$, where $\ell$ is the SM
fermion, can still remain open through a virtual $Z/Z^{\prime}$ exchange. For
the chosen mass gap of $\lesssim 10$ keV, $\ell=\nu$ is possible. However, for
our benchmark cases with such a small gap, the lifetime of $X_{R}$,
$\displaystyle\Gamma_{X_{R}\rightarrow X_{I},\nu,\bar{\nu}}\approx 6.06\times
10^{-46}\;{\rm GeV}\;\left(\frac{\epsilon}{2\times
10^{-4}}\right)^{2}g_{X}^{2}\left(\frac{m_{X_{R}}-m_{X_{I}}}{10\;{\rm
keV}}\right)^{5}\,,$ is longer than the age of the Universe as pointed out in
Ref. Baek:2020owl . We shall thus ignore this decay channel in the following.
Note also that the partial width of the $X_{R}$ decaying into $X_{I}$ and
three photons, as shown in Ref. Harigaya:2020ckz , is around twelve orders of
magnitude smaller than $\Gamma_{X_{R}\rightarrow X_{I},\nu,\bar{\nu}}$. .
Thus, together with $Y$ and $X_{I}$, we have a three-component scenario. We
consider two cases: iii) $m_{X_{R}}\approx m_{X_{I}}\gg m_{Y}$ and iv)
$m_{X_{R}}\approx m_{X_{I}}\approx m_{Y}$. The input parameter values are
shown in Table 1. We chose a large dark gauge coupling,
$g_{X}\sim\mathcal{O}(1)$.

While the decay channels of the field $X_{R}$ are kinematically forbidden, it
is still possible that the field $X_{R}$ may be converted into the field
$X_{I}$ via the conversion process $X_{R},Y\rightarrow X_{I},Y$. We
numerically solve the full Boltzmann equation and found that the mass gap
should be less than $\sim$10 keV, i.e. $m_{X_{R}}-m_{X_{I}}\lesssim 10$ keV to
suppress such a process, thereby achieving three-component scenarios.

The solutions of the Boltzmann equations are shown in Fig. 3. Since
$m_{X_{R}}\approx m_{X_{I}}$ we do not see any clear deviations between the
$X_{R}$ evolution and $X_{I}$ evolution. Similar to the cases i) and ii), we
see that, in the case of $m_{X_{I}}\gg m_{Y}$, $X_{I}$ freezes out first at
$x\approx 5$, followed by the $Y$ freeze-out at $x\approx 20$. When
$m_{X_{I}}\approx m_{Y}$, both $X_{I}$ and $Y$ freeze out at almost the same
time $x\approx 20$. In the case iii), we find that $X_{I}$ and $X_{R}$ ($Y$)
constitute about 20% each (60%) of the total relic density, while in the case
iv), the fractions of $X_{I}$ and $X_{R}$ relics are about 9% each, with 82%
of $Y$ relic. Thus we identify these scenarios as three-component DM
scenarios. The self-scattering cross sections are given by $1.66$ ${\rm
cm^{2}/g}$ and $4.04$ ${\rm cm^{2}/g}$ for cases ${\it iii)}$ and ${\it iv)}$,
respectively. We summarise the results in Table 2.

Figure 3: The solutions of the Boltzmann equations summarised in Appendix B are shown for cases i)–iv). The red, blue, black, and brown solid lines are the yields of $Y$, $X_{I}$, $X_{R}$, and $Z^{\prime}$, respectively; here $Y_{y}\equiv 2Y_{Y}=2Y_{Y^{*}}$. The dashed lines correspond to the equilibrium states. In the cases i) and ii) $Y$ and $X_{I}$ relics become frozen out, indicating two-component DM scenarios. On the other hand, in the cases iii) and iv) all the DM candidates freeze out, and thus we have three-component scenarios. We note that there is no visible difference between the $X_{I}$ relic and $X_{R}$ relic in the cases iii) and iv) due to the small mass gap. A large gauge coupling $g_{X}\sim\mathcal{O}(1)$ is chosen for our benchmark points. In this case, once the mass gap becomes larger than $\sim$10 keV, a significant amount of the relic of $X_{R}$ is converted into the $X_{I}$ relic, mainly through $X_{R},Y\rightarrow X_{I},Y$, becoming a two-component scenario. In all cases we see that $Z^{\prime}$ follows its equilibrium state. Case | $\Omega_{X_{R}}/\Omega_{\rm DM}$ | $\Omega_{X_{I}}/\Omega_{\rm DM}$ | $\Omega_{y}/\Omega_{\rm DM}$ |  $\sigma_{\rm self}/m_{\rm DM}$ (${\rm cm^{2}/g}$) | DM scenario
---|---|---|---|---|---
i) | 0 | 0.237 | 0.763 | 1.49 | Two-component
ii) | 0 | 0.197 | 0.803 | 3.83 | Two-component
iii) | 0.186 | 0.207 | 0.607 | 1.66 | Three-component
iv) | 0.087 | 0.088 | 0.825 | 4.04 | Three-component
Table 2: The fractions of relic density for each DM candidate field and the
self-scattering cross sections are shown for the four benchmark cases. The
input parameters are summarised in Table 1. The total DM relic density is
$\Omega_{\rm DM}h^{2}=0.12$. As expected from the mass gap between $X_{I}$ and
$X_{R}$, cases i) and ii) give rise to two-component DM scenarios, while
three-component DM scenarios are realised in the cases iii) and iv). The self-
scattering cross sections are somewhat larger than 1 ${\rm cm^{2}/g}$ imposed
by the Bullet Cluster constraint Markevitch:2003at ; Clowe:2003tk ;
Randall:2007ph  (see also Refs. Rocha:2012jg ; Peter:2012jh  where the similar
bound is obtained from cosmological simulations with self-interacting DM), but
well within the bound, 10 ${\rm cm^{2}/g}$.

Let us comment on the parameter $\lambda_{XY}$ before we conclude. In our
benchmark cases i)–iv) we chose $\lambda_{XY}$ to be zero so that the contact
interaction between $X$ and $Y$ gets suppressed, helping realizations of the
multi-component DM scenarios in the large-$m_{h^{\prime}}$ limit. For the
chosen parameter values of the benchmark cases, we present the upper bound on
$\lambda_{XY}$ to achieve such multi-component scenarios. We first require
that $\Omega_{Y}/\Omega_{\rm DM}\lesssim 90\%$ be the condition for a multi-
component DM scenario. For the case i), we obtained the maximum value of
$\lambda_{XY}=1.5\times 10^{-4}$ with $\lambda_{Y}=5.45$, in which case the
relic of the $Y$ ($X_{I}$) field consists 90% (10%) of the total relic
density. For the case ii), the maximum value is given by
$\lambda_{XY}=4.7\times 10^{-5}$ with $\lambda_{Y}=4.62$. The $Y$ ($X_{I}$)
relic consists 90% (10%) of the total relic density. For the case iii) (case
iv)), $\lambda_{XY}=4.0\times 10^{-4}$ ($3.8\times 10^{-4}$) is the maximum
value, with $\lambda_{Y}=3.71$ (4.95) that gives $\Omega_{Y}/\Omega_{\rm
DM}=90\%$ and $\Omega_{X_{I,R}}/\Omega_{\rm DM}=5\%$. Note that the value of
$\lambda_{Y}$ is modified accordingly in order to satisfy the correct relic
density of $\Omega_{\rm DM}h^{2}=0.12$. The rest of the input parameters are
kept to be the same as in Table 1. As such, the value of $\lambda_{XY}$ needs
to be small to realise multi-component DM scenarios. This is not a general
statement, however, since the dark Higgs field is assumed to be massive
compared to the other dark fields. If the dark Higgs is chosen to be somewhat
lighter and/or the couplings $\lambda_{X\phi}$ and $\lambda_{Y\phi}$ are
somewhat larger, the value of $\lambda_{XY}$ may become larger, due to
different signs in terms between the contact interaction and the dark Higgs
mediated processes, causing cancellation, i.e. destructive interferences. In
Table 3 we present such cases; we choose $m_{Y}=37.5$ MeV,
$m_{Z^{\prime}}=200$ MeV, $m_{h^{\prime}}=10$ GeV, $g_{X}=1$, $\lambda_{X}=1$,
$\lambda_{X\phi}=\lambda_{Y\phi}=10$, and $\lambda_{Y}=5$ for all the four
cases.

Case | $\frac{m_{X_{R}}}{m_{X_{I}}}$ | $\frac{m_{X_{I}}}{m_{Y}}$ | $\lambda^{\prime}_{Y\phi}$ | $\lambda_{XY}$ | $\frac{\Omega_{X_{R}}}{\Omega_{\rm DM}}$ | $\frac{\Omega_{X_{I}}}{\Omega_{\rm DM}}$ | $\frac{\Omega_{y}}{\Omega_{\rm DM}}$ |  $\frac{\sigma_{\rm self}}{m_{\rm DM}}$ (${\rm cm^{2}/g}$)
---|---|---|---|---|---|---|---|---
i′) | 4 | 4 | 0.2 | 0.02 | 0 | 0.30 | 0.70 | 1.79
ii′) | 10 | 1.07 | 0.29 | 0.03 | 0 | 0.21 | 0.79 | 5.46
iii′) | 1.0001 | 4 | 0.29 | 0.04 | 0.24 | 0.28 | 0.48 | 2.07
iv′) | 1.0001 | 1.07 | 0.29 | 0.04 | 0.09 | 0.10 | 0.81 | 5.74
Table 3:  Input parameters with non-zero $\lambda_{XY}$ and outcomes of DM
relics and the self-scattering cross sections. The rest of the input
parameters are chosen as follows: $m_{Y}=37.5$ MeV, $m_{Z^{\prime}}=200$ MeV,
$m_{h^{\prime}}=10$ GeV, $g_{X}=1$, $\lambda_{X}=0.1$,
$\lambda_{X\phi}=\lambda_{Y\phi}=10$, and $\lambda_{Y}=5$. All the cases are
free from the constraints listed in Section 3.2. We see that the multi-
component DM scenarios are still realised with sizeable $\lambda_{XY}$ values.
This is due to the destructive interference between the $X$–$Y$ contact
interaction and the dark Higgs mediated processes.

##  4 Conclusions

In this paper, we constructed multi-component dark matter scenarios with
$U(1)_{X}$ dark gauge symmetry broken into $Z_{2}\times Z_{3}$. Two
independent types of dark matter fields, namely $X$ and $Y$, distinguished by
their different charges under $Z_{2}\times Z_{3}$, emerge as separately stable
dark matter candidates. The relic densities of different dark matter fields
are governed by interactions between the dark fields within the dark sector
through number-changing $2\rightarrow 2$ and $3\rightarrow 2$ processes.
Unlike many other SIMP models in the literature, $X$ and $Y$ are completely
independent of each other, with different masses and spins in principle. Due
to this reason we call our model a genuine multi-component SIMP dark matter
scenario.

A nonzero vacuum expectation value of the dark Higgs field gives mass
splitting to the real and imaginary parts of the $X$ field. When the mass gap
between the real part $X_{R}$ and imaginary part $X_{I}$ is large, only the
lighter field $X_{I}$ becomes a stable dark matter candidate which, together
with $Y$, gives rise to a two-component dark matter scenario. On the other
hand, once the mass gap is small enough, both $X_{R}$ and $X_{I}$ may
contribute to the dark matter relic density; we thus have a three-component
scenario. In our analysis, the dark gauge coupling $g_{X}$ is chosen to be
$\mathcal{O}(1)$, and we found that the mass gap should be less than $\sim$10
keV, i.e. $m_{X_{R}}-m_{X_{I}}\lesssim 10$ keV, in order to achieve three-
component scenarios.

The unbroken $Z_{3}$ symmetry ensures the stability of the $Y$ field, which
becomes another dark matter candidate of the SIMP-type. Due to the
interactions between $X_{I,R}$ and $Y$, both $X_{I}$ and $X_{R}$ may
contribute to $Y$’s $3\rightarrow 2$ annihilation processes through the
$Z^{\prime}$ and the dark Higgs. These kinds of processes open a new class of
SIMP models.

We presented four choices of the parameter set for the multi-component dark
matter scenarios, two for a two-component dark matter scenario, i)
$m_{X_{R}}\gg m_{X_{I}}\gg m_{Y}$ and ii) $m_{X_{R}}\gg m_{X_{I}}\approx
m_{Y}$, and two for a three-component dark matter scenario, iii)
$m_{X_{R}}\approx m_{X_{I}}\gg m_{Y}$ and iv) $m_{X_{R}}\approx
m_{X_{I}}\approx m_{Y}$, by solving the full coupled Boltzmann equations and
considering both theoretical and experimental constraints. Figures 1 and 2
show a set of constraints. The parameter choices are given in Table 1 and the
solutions of the Boltzmann equations are shown in Fig. 3. For each case we
computed fractions of the relic density. The results, together with values of
the self-scattering cross section, are summarised in Table 2. We found that,
in the large-$m_{h^{\prime}}$ limit, the $\lambda_{XY}$ coupling, which
dilutes the relic density of the $X$ field via the $X$–$Y$ contact
interaction, needs to take a small value to realise a multi-component dark
matter scenario. On the other hand, the $\lambda_{XY}$ coupling may take a
larger value once the dark Higgs field becomes lighter and the dark Higgs
portal couplings, $\lambda_{X\phi}$ and $\lambda_{Y\phi}$, become larger. In
this case destructive interferences occur and the dilution of the $X$ relic is
suppressed even with a larger value of $\lambda_{XY}$. In Table 3 we present
four parameter sets of this case.

Our findings of the multi-component SIMP-type dark matter scenarios are new
and interesting. One of the most distinctive features of our model is that the
different dark matter fields $X$ and $Y$ may have very different mass scales,
while significantly contributing to the total dark matter relic density today,
as demonstrated in the cases i) and iii). In the current work a scalar
$Z_{2}$-charged field is considered. It is certainly possible to have a
fermionic $Z_{2}$-charged field 888 Note that the $Z_{3}$-charged dark matter
cannot be fermionic.  as studied e.g. in Refs. Baek:2014poa ; Ko:2019wxq ;
Baek:2020owl . Furthermore, we investigated four characteristic show-cases for
the multi-component dark matter scenario. In order to fully understand the
model and its phenomenology, an intensive parameter scan is better to be
performed. We shall tackle them in future works.

##  Appendix A Three-body Decay of $X_{R}$

Consider a decay of particle 0 to particles 1, 2, and 3:

$\displaystyle 0\rightarrow 1,2,3\,.$

The evolution of the number density of particle 0 is then governed by the
following Boltzmann equation:

$\displaystyle\frac{dn_{0}}{dt}+3Hn_{0}$ $\displaystyle=\int
d\Pi_{0}d\Pi_{1}d\Pi_{2}d\Pi_{3}(2\pi)^{4}\delta^{(4)}(p_{0}-p_{1}-p_{2}-p_{3})$
$\displaystyle\quad\times\left[f_{1}f_{2}f_{3}(1\pm
f_{0})|\mathcal{M}|^{2}_{1,2,3\rightarrow 0}-f_{0}(1\pm f_{1})(1\pm
f_{2})(1\pm f_{3})|\mathcal{M}|^{2}_{0\rightarrow 1,2,3}\right]$
$\displaystyle\approx\int
d\Pi_{0}d\Pi_{1}d\Pi_{2}d\Pi_{3}(2\pi)^{4}\delta^{(4)}(p_{0}-p_{1}-p_{2}-p_{3})|\mathcal{M}|^{2}_{0\rightarrow
1,2,3}$ $\displaystyle\quad\times\left(f_{1}f_{2}f_{3}-f_{0}\right)\,,$ (55)

where

$\displaystyle d\Pi_{i}\equiv\frac{d^{3}p_{i}}{(2\pi)^{3}2E_{i}}\,,$ (56)

and we assumed $1\pm f_{i}\approx 1$ and CP conservation, i.e.
$|\mathcal{M}|^{2}_{0\rightarrow 1,2,3}=|\mathcal{M}|^{2}_{1,2,3\rightarrow
0}$. Now, using $f_{i}=(n_{i}/n_{i}^{\rm eq})f_{i}^{\rm eq}$ and $f_{0}^{\rm
eq}=f_{1}^{\rm eq}f_{2}^{\rm eq}f_{3}^{\rm eq}$, we obtain

$\displaystyle\frac{dn_{0}}{dt}+3Hn_{0}$ $\displaystyle\approx-\int
d\Pi_{0}d\Pi_{1}d\Pi_{2}d\Pi_{3}(2\pi)^{4}\delta^{(4)}(p_{0}-p_{1}-p_{2}-p_{3})|\mathcal{M}|^{2}_{0\rightarrow
1,2,3}$ $\displaystyle\quad\times\frac{f_{0}^{\rm eq}}{n_{0}^{\rm
eq}}\left(n_{0}-n_{0}^{\rm eq}\frac{n_{1}n_{2}n_{3}}{n_{1}^{\rm eq}n_{2}^{\rm
eq}n_{3}^{\rm eq}}\right)\,.$ (57)

Defining

$\displaystyle\Gamma_{0\rightarrow 1,2,3}\equiv\frac{1}{2m_{0}}\int
d\Pi_{1}d\Pi_{2}d\Pi_{3}(2\pi)^{4}\delta^{(4)}(p_{0}-p_{1}-p_{2}-p_{3})|\mathcal{M}|^{2}_{0\rightarrow
1,2,3}\,,$ (58)

and the thermally-averaged decay rate

$\displaystyle\langle\Gamma\rangle_{0\rightarrow 1,2,3}\equiv\int
d\Pi_{0}2m_{0}\Gamma_{0\rightarrow 1,2,3}\frac{f_{0}^{\rm eq}}{n_{0}^{\rm
eq}}\,,$ (59)

the Boltzmann equation can be written as

$\displaystyle\frac{dn_{0}}{dt}+3Hn_{0}=-\langle\Gamma\rangle_{0\rightarrow
1,2,3}\left(n_{0}-n_{0}^{\rm eq}\frac{n_{1}n_{2}n_{3}}{n_{1}^{\rm
eq}n_{2}^{\rm eq}n_{3}^{\rm eq}}\right)\,.$ (60)

In terms of the yields, the Boltzmann equation can be rewritten as follows:

$\displaystyle\frac{dY_{0}}{dx}=-\frac{x}{H(m)}\langle\Gamma\rangle_{0\rightarrow
1,2,3}\left(Y_{0}-Y_{0}^{\rm eq}\frac{Y_{1}Y_{2}Y_{3}}{Y_{1}^{\rm
eq}Y_{2}^{\rm eq}Y_{3}^{\rm eq}}\right)\,,$ (61)

where $x\equiv m/T$ with $m$ being a mass variable.

For the decay process $X_{R}\rightarrow X_{I}+Y+Y^{*}$, we have

$\displaystyle\frac{dY_{X_{R}}}{dx}=-\frac{x}{H(m)}\langle\Gamma\rangle_{X_{R}\rightarrow
X_{I},Y,Y^{*}}\left(Y_{X_{R}}-Y_{X_{R}}^{\rm
eq}\frac{Y_{X_{I}}Y_{y}^{2}}{Y_{X_{I}}^{\rm eq}(Y_{y}^{\rm
eq})^{2}}\right)\,,$ (62)

where $Y_{y}=Y_{Y}+Y_{Y^{*}}=2Y_{Y}$. The decay rate is given by

$\displaystyle\langle\Gamma\rangle_{X_{R}\rightarrow X_{I},Y,Y^{*}}$
$\displaystyle=\frac{1}{64\pi^{3}m_{X_{R}}}\int_{E_{1}^{\rm min}}^{E_{1}^{\rm
max}}\int_{E_{2}^{\rm min}}^{E_{2}^{\rm
max}}|\mathcal{M}|^{2}_{X_{R}\rightarrow X_{I},Y,Y^{*}}dE_{2}dE_{1}\,,$ (63)
$\displaystyle|\mathcal{M}|^{2}_{X_{R}\rightarrow X_{I},Y,Y^{*}}$
$\displaystyle=\frac{g_{X}^{4}m_{X_{R}}^{2}(E_{1}+2E_{2}-m_{X_{R}})^{2}}{9(m_{X_{I}}^{2}-2m_{X_{R}}E_{1}+m_{X_{R}}^{2}-m_{Z^{\prime}}^{2})^{2}}\,,$
(64)

with

$\displaystyle E_{2}^{\rm min}$
$\displaystyle=\frac{1}{2m_{23}^{2}}\left[(m_{X_{R}}-E_{1})m_{23}^{2}-\sqrt{(E_{1}^{2}-m_{X_{I}}^{2})\lambda(m_{23}^{2},m_{Y}^{2},m_{Y}^{2})}\right]\,,$
(65) $\displaystyle E_{2}^{\rm max}$
$\displaystyle=\frac{1}{2m_{23}^{2}}\left[(m_{X_{R}}-E_{1})m_{23}^{2}+\sqrt{(E_{1}^{2}-m_{X_{I}}^{2})\lambda(m_{23}^{2},m_{Y}^{2},m_{Y}^{2})}\right]\,,$
(66)

and

$\displaystyle E_{1}^{\rm min}=m_{X_{I}}\,,\qquad E_{1}^{\rm
max}=\frac{m_{X_{R}}^{2}+m_{X_{I}}^{2}-4m_{Y}^{2}}{2m_{X_{R}}}\,,$ (67)

where

$\displaystyle m_{23}^{2}$
$\displaystyle=m_{X_{R}}^{2}-2m_{X_{R}}E_{1}+m_{X_{I}}^{2}\,,$ (68)
$\displaystyle\lambda(x,y,z)$ $\displaystyle\equiv
x^{2}+y^{2}+z^{2}-2xy-2xz-2yz\,.$ (69)

##  Appendix B The Full Boltzmann Equations

In this appendix, we present the full Boltzmann equations for $Y$, $X_{R}$,
$X_{I}$, and $Z^{\prime}$. The Boltzmann equation for $Y$ is, using
$Y_{Y}=Y_{Y^{*}}=Y_{y}/2$, given by

$\displaystyle\frac{dY_{y}}{dx}$
$\displaystyle=\frac{2x}{H}\left[\langle\Gamma\rangle_{Z^{\prime}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}-\frac{Y_{Z^{\prime}}^{\rm eq}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)+\langle\Gamma\rangle_{X_{R}\rightarrow
X_{I},Y,Y^{*}}\left(Y_{X_{R}}-\frac{Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}Y_{X_{I}}Y_{y}^{2}\right)\right]$
$\displaystyle\quad+\frac{2s}{Hx^{2}}\Bigg{[}\langle\sigma
v\rangle_{Z^{\prime},Z^{\prime}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}^{2}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)+\langle\sigma
v\rangle_{X_{R},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{X_{R}}^{2}-\frac{(Y_{X_{R}}^{\rm eq})^{2}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)$ $\displaystyle\quad+\langle\sigma
v\rangle_{X_{I},X_{I}\rightarrow
Y,Y^{*}}\left(Y_{X_{I}}^{2}-\frac{(Y_{X_{I}}^{\rm eq})^{2}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)+\langle\sigma v\rangle_{X_{I},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{X_{I}}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)$
$\displaystyle\quad+\frac{1}{2}\langle\sigma v\rangle_{Z^{\prime},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{Z^{\prime}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)\Bigg{]}$
$\displaystyle\quad+\frac{2s^{2}}{Hx^{5}}\Bigg{[}-\frac{1}{8}\langle\sigma
v^{2}\rangle_{Y,Y,Y^{*}\rightarrow Y^{*},Y^{*}}\left(Y_{y}^{3}-Y_{y}^{\rm
eq}Y_{y}^{2}\right)-\frac{1}{8}\langle\sigma v^{2}\rangle_{Y,Y,Y\rightarrow
Y,Y^{*}}\left(Y_{y}^{3}-Y_{y}^{\rm eq}Y_{y}^{2}\right)$
$\displaystyle\quad-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{I},Y,Y\rightarrow
X_{I},Y^{*}}\left(Y_{X_{I}}Y_{y}^{2}-Y_{y}^{\rm eq}Y_{X_{I}}Y_{y}\right)$
$\displaystyle\quad-\frac{3}{8}\langle\sigma v^{2}\rangle_{Y,Y,Y\rightarrow
X_{I},X_{R}}\left(Y_{y}^{3}-\frac{(Y_{y}^{\rm eq})^{3}}{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{I},Y,Y\rightarrow
X_{R},Y^{*}}\left(Y_{X_{I}}Y_{y}^{2}-\frac{Y_{X_{I}}^{\rm eq}Y_{y}^{\rm
eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}Y_{y}\right)$
$\displaystyle\quad+\frac{1}{2}\langle\sigma
v^{2}\rangle_{X_{I},X_{I},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{X_{I}}^{2}Y_{y}-\frac{(Y_{X_{I}}^{\rm
eq})^{2}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y^{*}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}Y_{y}^{2}-\frac{Y_{X_{R}}^{\rm eq}(Y_{y}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$
$\displaystyle\quad-\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y\rightarrow
X_{I},Y^{*}}\left(Y_{X_{R}}Y_{y}^{2}-\frac{Y_{X_{R}}^{\rm eq}Y_{y}^{\rm
eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{y}\right)-\frac{3}{8}\langle\sigma
v^{2}\rangle_{Y,Y,Y\rightarrow X_{R},X_{R}}\left(Y_{y}^{3}-\frac{(Y_{y}^{\rm
eq})^{3}}{(Y_{X_{R}}^{\rm eq})^{2}}Y_{X_{R}}^{2}\right)$
$\displaystyle\quad-\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y\rightarrow
X_{R},Y^{*}}\left(Y_{X_{R}}Y_{y}^{2}-Y_{y}^{\rm
eq}Y_{X_{R}}Y_{y}\right)+\frac{1}{2}\langle\sigma
v^{2}\rangle_{X_{I},X_{R},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{X_{I}}Y_{X_{R}}Y_{y}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)$
$\displaystyle\quad-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{(Y_{X_{I}}^{\rm
eq})^{2}}Y_{X_{I}}^{2}\right)-\frac{3}{8}\langle\sigma
v^{2}\rangle_{Y,Y,Y\rightarrow X_{I},X_{I}}\left(Y_{y}^{3}-\frac{(Y_{y}^{\rm
eq})^{3}}{(Y_{X_{I}}^{\rm eq})^{2}}Y_{X_{I}}^{2}\right)$
$\displaystyle\quad+\frac{1}{2}\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Y^{*}\rightarrow
Y,Y}\left(Y_{Z^{\prime}}^{2}Y_{y}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)+\frac{1}{2}\langle\sigma
v^{2}\rangle_{X_{R},X_{R},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{X_{R}}^{2}Y_{y}-\frac{(Y_{X_{R}}^{\rm
eq})^{2}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)$
$\displaystyle\quad-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y^{*},Y^{*}\rightarrow
Z^{\prime},Y}\left(Y_{Z^{\prime}}Y_{y}^{2}-Y_{y}^{\rm
eq}Y_{Z^{\prime}}Y_{y}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{(Y_{X_{R}}^{\rm eq})^{2}}Y_{X_{R}}^{2}\right)$
$\displaystyle\quad+\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Z^{\prime}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}^{3}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{3}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)+\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm eq}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)$ $\displaystyle\quad+\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{R},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{R}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{X_{R}}^{\rm eq})^{2}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)+\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{I}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{I}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{X_{I}}^{\rm eq})^{2}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)\Bigg{]}\,.$
(70)

The Boltzmann equations for $X_{I}$ and $X_{R}$ are

$\displaystyle\frac{dY_{X_{I}}}{dx}$
$\displaystyle=\frac{x}{H}\Bigg{[}\langle\Gamma\rangle_{X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}-\frac{Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)+\langle\Gamma\rangle_{Z^{\prime}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}-\frac{Y_{Z^{\prime}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad+\langle\Gamma\rangle_{X_{R}\rightarrow
X_{I},Y,Y^{*}}\left(Y_{X_{R}}-\frac{Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}Y_{X_{I}}Y_{y}^{2}\right)\Bigg{]}$
$\displaystyle\quad+\frac{s}{Hx^{2}}\Bigg{[}2\langle\sigma
v\rangle_{Z^{\prime},Z^{\prime}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}^{2}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{(Y_{X_{I}}^{\rm eq})^{2}}Y_{X_{I}}^{2}\right)+2\langle\sigma
v\rangle_{X_{R},X_{R}\rightarrow
X_{I},X_{I}}\left(Y_{X_{R}}^{2}-\frac{(Y_{X_{R}}^{\rm
eq})^{2}}{(Y_{X_{I}}^{\rm eq})^{2}}Y_{X_{I}}^{2}\right)$
$\displaystyle\quad+\langle\sigma v\rangle_{X_{R},Y\rightarrow
X_{I},Y}\left(Y_{X_{R}}Y_{y}-\frac{Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{X_{I}}Y_{y}\right)-2\langle\sigma v\rangle_{X_{I},X_{I}\rightarrow
Z^{\prime},Z^{\prime}}\left(Y_{X_{I}}^{2}-\frac{(Y_{X_{I}}^{\rm
eq})^{2}}{(Y_{Z^{\prime}}^{\rm eq})^{2}}Y_{Z^{\prime}}^{2}\right)$
$\displaystyle\quad-2\langle\sigma v\rangle_{X_{I},X_{I}\rightarrow
Y,Y^{*}}\left(Y_{X_{I}}^{2}-\frac{(Y_{X_{I}}^{\rm eq})^{2}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)-\langle\sigma v\rangle_{X_{I},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{X_{I}}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)\Bigg{]}$
$\displaystyle\quad+\frac{s^{2}}{Hx^{5}}\Bigg{[}-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
Z^{\prime},Z^{\prime}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm
eq}}Y_{Z^{\prime}}^{2}\right)-\langle\sigma
v^{2}\rangle_{X_{I},X_{R},Y\rightarrow
Z^{\prime},Y}\left(Y_{X_{I}}Y_{X_{R}}Y_{y}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{y}\right)$
$\displaystyle\quad+\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)+\frac{1}{4}\langle\sigma
v^{2}\rangle_{Y,Y,Y\rightarrow X_{I},X_{R}}\left(Y_{y}^{3}-\frac{(Y_{y}^{\rm
eq})^{3}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad-\frac{1}{2}\langle\sigma
v^{2}\rangle_{X_{I},Y,Y\rightarrow
X_{R},Y^{*}}\left(Y_{X_{I}}Y_{y}^{2}-\frac{Y_{X_{I}}^{\rm eq}Y_{y}^{\rm
eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}Y_{y}\right)-2\langle\sigma
v^{2}\rangle_{X_{I},X_{I},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{X_{I}}^{2}Y_{y}-\frac{(Y_{X_{I}}^{\rm
eq})^{2}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)$ $\displaystyle\quad+\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{Z^{\prime}}^{2}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)-\langle\sigma
v^{2}\rangle_{X_{I},X_{I},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{I}}^{2}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$
$\displaystyle\quad+\langle\sigma v^{2}\rangle_{X_{R},X_{R},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}^{3}-\frac{(Y_{X_{R}}^{\rm
eq})^{3}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)+\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y^{*}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}Y_{y}^{2}-\frac{Y_{X_{R}}^{\rm eq}(Y_{y}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$ $\displaystyle\quad+\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{X_{I}}^{2}\right)+\frac{1}{2}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y\rightarrow
X_{I},Y^{*}}\left(Y_{X_{R}}Y_{y}^{2}-\frac{Y_{X_{R}}^{\rm eq}Y_{y}^{\rm
eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{y}\right)$
$\displaystyle\quad-\langle\sigma v^{2}\rangle_{X_{I},X_{R},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{X_{I}}Y_{X_{R}}Y_{y}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},X_{I}\rightarrow
Z^{\prime},X_{R}}\left(Y_{Z^{\prime}}^{2}Y_{X_{I}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)$
$\displaystyle\quad-3\langle\sigma v^{2}\rangle_{X_{I},X_{I},X_{I}\rightarrow
Z^{\prime},X_{R}}\left(Y_{X_{I}}^{3}-\frac{(Y_{X_{I}}^{\rm
eq})^{3}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)-\langle\sigma
v^{2}\rangle_{X_{I},X_{R},X_{R}\rightarrow
Z^{\prime},X_{R}}\left(Y_{X_{I}}Y_{X_{R}}^{2}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)$
$\displaystyle\quad+\frac{1}{2}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{(Y_{X_{I}}^{\rm
eq})^{2}}Y_{X_{I}}^{2}\right)+\frac{1}{2}\langle\sigma
v^{2}\rangle_{Y,Y,Y\rightarrow X_{I},X_{I}}\left(Y_{y}^{3}-\frac{(Y_{y}^{\rm
eq})^{3}}{(Y_{X_{I}}^{\rm eq})^{2}}Y_{X_{I}}^{2}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}^{2}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},Y\rightarrow
X_{R},Y}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}Y_{y}\right)$
$\displaystyle\quad-2\langle\sigma v^{2}\rangle_{X_{I},X_{I},Y\rightarrow
Z^{\prime},Y}\left(Y_{X_{I}}^{2}Y_{y}-\frac{(Y_{X_{I}}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{y}\right)+\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Z^{\prime}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}^{3}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{3}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{I}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{I}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)+\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{R},X_{R}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{R}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad+\langle\sigma v^{2}\rangle_{Z^{\prime},X_{R},Y\rightarrow
X_{I},Y}\left(Y_{Z^{\prime}}Y_{X_{R}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{y}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm eq}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)$ $\displaystyle\quad-2\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{I}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}(Y_{X_{I}})^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{X_{I}}^{\rm eq})^{2}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)\Bigg{]}\,,$
(71)

and

$\displaystyle\frac{dY_{X_{R}}}{dx}$
$\displaystyle=\frac{x}{H}\Bigg{[}\langle\Gamma\rangle_{Z^{\prime}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}-\frac{Y_{Z^{\prime}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)-\langle\Gamma\rangle_{X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}-\frac{Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$
$\displaystyle\quad-\langle\Gamma\rangle_{X_{R}\rightarrow
X_{I},Y,Y^{*}}\left(Y_{X_{R}}-\frac{Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}Y_{X_{I}}Y_{y}^{2}\right)\Bigg{]}$
$\displaystyle\quad+\frac{s}{Hx^{2}}\Bigg{[}2\langle\sigma
v\rangle_{Z^{\prime},Z^{\prime}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}^{2}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{(Y_{X_{R}}^{\rm eq})^{2}}Y_{X_{R}}^{2}\right)-2\langle\sigma
v\rangle_{X_{R},X_{R}\rightarrow
Z^{\prime},Z^{\prime}}\left(Y_{X_{R}}^{2}-\frac{(Y_{X_{R}}^{\rm
eq})^{2}}{(Y_{Z^{\prime}}^{\rm eq})^{2}}Y_{Z^{\prime}}^{2}\right)$
$\displaystyle\quad-2\langle\sigma v\rangle_{X_{R},X_{R}\rightarrow
X_{I},X_{I}}\left(Y_{X_{R}}^{2}-\frac{(Y_{X_{R}}^{\rm
eq})^{2}}{(Y_{X_{I}}^{\rm eq})^{2}}Y_{X_{I}}^{2}\right)-2\langle\sigma
v\rangle_{X_{R},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{X_{R}}^{2}-\frac{(Y_{X_{R}}^{\rm eq})^{2}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)$ $\displaystyle\quad-\langle\sigma
v\rangle_{X_{R},Y\rightarrow X_{I},Y}\left(Y_{X_{R}}Y_{y}-\frac{Y_{X_{R}}^{\rm
eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{y}\right)-\langle\sigma
v\rangle_{X_{I},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{X_{I}}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)\Bigg{]}$
$\displaystyle\quad+\frac{s^{2}}{Hx^{5}}\Bigg{[}-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
Z^{\prime},Z^{\prime}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm
eq}}Y_{Z^{\prime}}^{2}\right)-\langle\sigma
v^{2}\rangle_{X_{I},X_{R},Y\rightarrow
Z^{\prime},Y}\left(Y_{X_{I}}Y_{X_{R}}Y_{y}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{y}\right)$
$\displaystyle\quad-2\langle\sigma v^{2}\rangle_{X_{R},X_{R},Y\rightarrow
Z^{\prime},Y}\left(Y_{X_{R}}^{2}Y_{y}-\frac{(Y_{X_{R}}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm
eq}}Y_{Z^{\prime}}Y_{y}\right)+\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)$ $\displaystyle\quad+\frac{1}{4}\langle\sigma
v^{2}\rangle_{Y,Y,Y\rightarrow X_{I},X_{R}}\left(Y_{y}^{3}-\frac{(Y_{y}^{\rm
eq})^{3}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)+\frac{1}{2}\langle\sigma
v^{2}\rangle_{X_{I},Y^{*},Y^{*}\rightarrow
X_{R},Y}\left(Y_{X_{I}}Y_{y}^{2}-\frac{Y_{X_{I}}^{\rm eq}Y_{y}^{\rm
eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}Y_{y}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{Z^{\prime}}^{2}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)-\langle\sigma
v^{2}\rangle_{X_{I},X_{I},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{I}}^{2}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$
$\displaystyle\quad-3\langle\sigma v^{2}\rangle_{X_{R},X_{R},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}^{3}-\frac{(Y_{X_{R}}^{\rm
eq})^{3}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y^{*}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}Y_{y}^{2}-\frac{Y_{X_{R}}^{\rm eq}(Y_{y}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$ $\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{X_{I}}^{2}\right)-\frac{1}{2}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y\rightarrow
X_{I},Y^{*}}\left(Y_{X_{R}}Y_{y}^{2}-\frac{Y_{X_{R}}^{\rm eq}Y_{y}^{\rm
eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{y}\right)$
$\displaystyle\quad+\frac{1}{2}\langle\sigma v^{2}\rangle_{Y,Y,Y\rightarrow
X_{R},X_{R}}\left(Y_{y}^{3}-\frac{(Y_{y}^{\rm eq})^{3}}{(Y_{X_{R}}^{\rm
eq})^{2}}Y_{X_{R}}^{2}\right)-\langle\sigma
v^{2}\rangle_{X_{I},X_{R},Y^{*}\rightarrow
Y,Y}\left(Y_{X_{I}}Y_{X_{R}}Y_{y}-\frac{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)$ $\displaystyle\quad+\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},X_{I}\rightarrow
Z^{\prime},X_{R}}\left(Y_{Z^{\prime}}^{2}Y_{X_{I}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)+\langle\sigma
v^{2}\rangle_{X_{I},X_{I},X_{I}\rightarrow
Z^{\prime},X_{R}}\left(Y_{X_{I}}^{3}-\frac{(Y_{X_{I}}^{\rm
eq})^{3}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)$ $\displaystyle\quad-\langle\sigma
v^{2}\rangle_{X_{I},X_{R},X_{R}\rightarrow
Z^{\prime},X_{R}}\left(Y_{X_{I}}Y_{X_{R}}^{2}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)+\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}^{2}\right)$
$\displaystyle\quad+\langle\sigma v^{2}\rangle_{Z^{\prime},X_{I},Y\rightarrow
X_{R},Y}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}Y_{y}\right)-2\langle\sigma
v^{2}\rangle_{X_{R},X_{R},Y^{*}\rightarrow
Y,Y}\left(Y_{X_{R}}^{2}Y_{y}-\frac{(Y_{X_{R}}^{\rm eq})^{2}}{Y_{y}^{\rm
eq}}Y_{y}^{2}\right)$ $\displaystyle\quad+\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Z^{\prime}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}^{3}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{3}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)+\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{I}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{I}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{R},X_{R}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{R}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{R},Y\rightarrow
X_{I},Y}\left(Y_{Z^{\prime}}Y_{X_{R}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{y}\right)$
$\displaystyle\quad+\frac{1}{2}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{(Y_{X_{R}}^{\rm
eq})^{2}}Y_{X_{R}}^{2}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm eq}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)$ $\displaystyle\quad-2\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{R},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{R}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{X_{R}}^{\rm eq})^{2}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)\Bigg{]}\,.$
(72)

Finally, the Boltzmann equation for $Z^{\prime}$ is given by

$\displaystyle\frac{dY_{Z^{\prime}}}{dx}$
$\displaystyle=\frac{x}{H}\Bigg{[}-\langle\Gamma\rangle_{Z^{\prime}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}-\frac{Y_{Z^{\prime}}^{\rm eq}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)-\langle\Gamma\rangle_{Z^{\prime}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}-\frac{Y_{Z^{\prime}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad+\langle\Gamma\rangle_{X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}-\frac{Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)\Bigg{]}$
$\displaystyle\quad+\frac{s}{Hx^{2}}\Bigg{[}-2\langle\sigma
v\rangle_{Z^{\prime},Z^{\prime}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}^{2}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{(Y_{X_{I}}^{\rm eq})^{2}}Y_{X_{I}}^{2}\right)-2\langle\sigma
v\rangle_{Z^{\prime},Z^{\prime}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}^{2}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{(Y_{X_{R}}^{\rm eq})^{2}}Y_{X_{R}}^{2}\right)$
$\displaystyle\quad-2\langle\sigma v\rangle_{Z^{\prime},Z^{\prime}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}^{2}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)+2\langle\sigma
v\rangle_{X_{R},X_{R}\rightarrow
Z^{\prime},Z^{\prime}}\left(Y_{X_{R}}^{2}-\frac{(Y_{X_{R}}^{\rm
eq})^{2}}{(Y_{Z^{\prime}}^{\rm eq})^{2}}Y_{Z^{\prime}}^{2}\right)$
$\displaystyle\quad-\langle\sigma v\rangle_{Z^{\prime},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{Z^{\prime}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)+2\langle\sigma
v\rangle_{X_{I},X_{I}\rightarrow
Z^{\prime},Z^{\prime}}\left(Y_{X_{I}}^{2}-\frac{(Y_{X_{I}}^{\rm
eq})^{2}}{(Y_{Z^{\prime}}^{\rm eq})^{2}}Y_{Z^{\prime}}^{2}\right)\Bigg{]}$
$\displaystyle\quad+\frac{s^{2}}{Hx^{5}}\Bigg{[}\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
Z^{\prime},Z^{\prime}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm
eq}}Y_{Z^{\prime}}^{2}\right)+\langle\sigma
v^{2}\rangle_{X_{I},X_{R},Y\rightarrow
Z^{\prime},Y}\left(Y_{X_{I}}Y_{X_{R}}Y_{y}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{y}\right)$
$\displaystyle\quad+\langle\sigma v^{2}\rangle_{X_{R},X_{R},Y\rightarrow
Z^{\prime},Y}\left(Y_{X_{R}}^{2}Y_{y}-\frac{(Y_{X_{R}}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm
eq}}Y_{Z^{\prime}}Y_{y}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)$ $\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{Z^{\prime}}^{2}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)+\langle\sigma
v^{2}\rangle_{X_{I},X_{I},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{I}}^{2}Y_{X_{R}}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$
$\displaystyle\quad+\langle\sigma v^{2}\rangle_{X_{R},X_{R},X_{R}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}^{3}-\frac{(Y_{X_{R}}^{\rm
eq})^{3}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)+\frac{1}{4}\langle\sigma
v^{2}\rangle_{X_{R},Y,Y^{*}\rightarrow
Z^{\prime},X_{I}}\left(Y_{X_{R}}Y_{y}^{2}-\frac{Y_{X_{R}}^{\rm eq}(Y_{y}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{I}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{I}}\right)$ $\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}^{2}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},X_{I}\rightarrow
Z^{\prime},X_{R}}\left(Y_{Z^{\prime}}^{2}Y_{X_{I}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)$
$\displaystyle\quad+\langle\sigma v^{2}\rangle_{X_{I},X_{I},X_{I}\rightarrow
Z^{\prime},X_{R}}\left(Y_{X_{I}}^{3}-\frac{(Y_{X_{I}}^{\rm
eq})^{3}}{Y_{Z^{\prime}}^{\rm eq}Y_{X_{R}}^{\rm
eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)+\langle\sigma
v^{2}\rangle_{X_{I},X_{R},X_{R}\rightarrow
Z^{\prime},X_{R}}\left(Y_{X_{I}}Y_{X_{R}}^{2}-\frac{Y_{X_{I}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{X_{R}}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Y\rightarrow
Z^{\prime},Y}\left(Y_{Z^{\prime}}^{2}Y_{y}-Y_{Z^{\prime}}^{\rm
eq}Y_{Z^{\prime}}Y_{y}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{I},X_{I}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{(Y_{X_{I}}^{\rm eq})^{2}}Y_{X_{I}}^{2}\right)$
$\displaystyle\quad-\langle\sigma v^{2}\rangle_{Z^{\prime},X_{I},Y\rightarrow
X_{I},Y}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{y}-Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}Y_{y}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}^{2}\right)$
$\displaystyle\quad-\langle\sigma v^{2}\rangle_{Z^{\prime},X_{I},Y\rightarrow
X_{R},Y}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm eq}}Y_{X_{R}}Y_{y}\right)-2\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Y\rightarrow
Y^{*},Y^{*}}\left(Y_{Z^{\prime}}^{2}Y_{y}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{2}}{Y_{y}^{\rm eq}}Y_{y}^{2}\right)$ $\displaystyle\quad+\langle\sigma
v^{2}\rangle_{X_{I},X_{I},Y\rightarrow
Z^{\prime},Y}\left(Y_{X_{I}}^{2}Y_{y}-\frac{(Y_{X_{I}}^{\rm
eq})^{2}}{Y_{Z^{\prime}}^{\rm eq}}Y_{Z^{\prime}}Y_{y}\right)-3\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Z^{\prime}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}^{3}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{3}}{Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{I}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{I}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}}{Y_{X_{R}}^{\rm
eq}}Y_{X_{I}}Y_{X_{R}}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{R},X_{R}\rightarrow
X_{I},X_{R}}\left(Y_{Z^{\prime}}Y_{X_{R}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm eq}}Y_{X_{I}}Y_{X_{R}}\right)$
$\displaystyle\quad-\langle\sigma v^{2}\rangle_{Z^{\prime},X_{R},Y\rightarrow
X_{I},Y}\left(Y_{Z^{\prime}}Y_{X_{R}}Y_{y}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}^{\rm eq}}{Y_{X_{I}}^{\rm
eq}}Y_{X_{I}}Y_{y}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
X_{R},X_{R}}\left(Y_{Z^{\prime}}Y_{y}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{y}^{\rm eq})^{2}}{(Y_{X_{R}}^{\rm eq})^{2}}Y_{X_{R}}^{2}\right)$
$\displaystyle\quad-\langle\sigma v^{2}\rangle_{Z^{\prime},X_{R},Y\rightarrow
X_{R},Y}\left(Y_{Z^{\prime}}Y_{X_{R}}Y_{y}-Y_{Z^{\prime}}^{\rm
eq}Y_{X_{R}}Y_{y}\right)-\frac{1}{2}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y\rightarrow
Y,Y}\left(Y_{Z^{\prime}}Y_{y}^{2}-Y_{Z^{\prime}}^{\rm eq}Y_{y}^{2}\right)$
$\displaystyle\quad-3\langle\sigma
v^{2}\rangle_{Z^{\prime},Z^{\prime},Z^{\prime}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}^{3}-\frac{(Y_{Z^{\prime}}^{\rm
eq})^{3}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)-\frac{1}{4}\langle\sigma
v^{2}\rangle_{Z^{\prime},Y,Y^{*}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{y}^{2}-Y_{Z^{\prime}}^{\rm eq}Y_{y}^{2}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{I}}Y_{X_{R}}-\frac{Y_{Z^{\prime}}^{\rm
eq}Y_{X_{I}}^{\rm eq}Y_{X_{R}}^{\rm eq}}{(Y_{y}^{\rm
eq})^{2}}Y_{y}^{2}\right)-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{R},X_{R}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{R}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{X_{R}}^{\rm eq})^{2}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)$
$\displaystyle\quad-\langle\sigma
v^{2}\rangle_{Z^{\prime},X_{I},X_{I}\rightarrow
Y,Y^{*}}\left(Y_{Z^{\prime}}Y_{X_{I}}^{2}-\frac{Y_{Z^{\prime}}^{\rm
eq}(Y_{X_{I}}^{\rm eq})^{2}}{(Y_{y}^{\rm eq})^{2}}Y_{y}^{2}\right)\Bigg{]}\,.$
(73)

###### Acknowledgements.

The work is supported in part by the DFG Collaborative Research Centre
“Neutrinos and Dark Matter in Astro- and Particle Physics” under Grant No. SFB
1258 (SMC), by KIAS Individual Grants, Grant No. PG021403 (PK), and by
National Research Foundation of Korea (NRF) Grant No. NRF-2019R1A2C3005009
(PK), funded by the Korea government (MSIT), and the Fundamental Research
Funds for the Central Universities, by the National Science Foundation of
China under Grant No. 11905149 (JL). The authors thank Chih-Ting Lu and Shu-Yu
Ho for mentioning on the $\lambda_{XY}$ term.

## References

  * (1) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo and F. S. Queiroz, “The waning of the WIMP? A review of models, searches, and constraints,” Eur. Phys. J. C 78, no.3, 203 (2018) [arXiv:1703.07364 [hep-ph]].
  * (2) L. Roszkowski, E. M. Sessolo and S. Trojanowski, “WIMP dark matter candidates and searches—current status and future prospects,” Rept. Prog. Phys. 81, no.6, 066201 (2018) [arXiv:1707.06277 [hep-ph]].
  * (3) Y. Hochberg, E. Kuflik, T. Volansky and J. G. Wacker, “Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 113, 171301 (2014) [arXiv:1402.5143 [hep-ph]].
  * (4) Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky and J. G. Wacker, “Model for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 115, no.2, 021301 (2015) [arXiv:1411.3727 [hep-ph]].
  * (5) N. Bernal, C. Garcia-Cely and R. Rosenfeld, “WIMP and SIMP Dark Matter from the Spontaneous Breaking of a Global Group,” JCAP 04, 012 (2015) [arXiv:1501.01973 [hep-ph]].
  * (6) H. M. Lee and M. S. Seo, “Communication with SIMP dark mesons via Z’ -portal,” Phys. Lett. B 748, 316-322 (2015) [arXiv:1504.00745 [hep-ph]].
  * (7) S. M. Choi and H. M. Lee, “SIMP dark matter with gauged Z3 symmetry,” JHEP 09, 063 (2015) [arXiv:1505.00960 [hep-ph]].
  * (8) Y. Hochberg, E. Kuflik and H. Murayama, “SIMP Spectroscopy,” JHEP 05, 090 (2016) [arXiv:1512.07917 [hep-ph]].
  * (9) S. M. Choi and H. M. Lee, “Resonant SIMP dark matter,” Phys. Lett. B 758, 47-53 (2016) [arXiv:1601.03566 [hep-ph]].
  * (10) S. M. Choi, Y. J. Kang and H. M. Lee, “On thermal production of self-interacting dark matter,” JHEP 12, 099 (2016) [arXiv:1610.04748 [hep-ph]].
  * (11) U. K. Dey, T. N. Maity and T. S. Ray, “Light Dark Matter through Assisted Annihilation,” JCAP 03, 045 (2017) [arXiv:1612.09074 [hep-ph]].
  * (12) S. M. Choi, H. M. Lee and M. S. Seo, “Cosmic abundances of SIMP dark matter,” JHEP 04, 154 (2017) [arXiv:1702.07860 [hep-ph]].
  * (13) S. M. Choi, Y. Hochberg, E. Kuflik, H. M. Lee, Y. Mambrini, H. Murayama and M. Pierre, “Vector SIMP dark matter,” JHEP 10, 162 (2017) [arXiv:1707.01434 [hep-ph]].
  * (14) S. M. Choi, H. M. Lee, P. Ko and A. Natale, “Resolving phenomenological problems with strongly-interacting-massive-particle models with dark vector resonances,” Phys. Rev. D 98, no.1, 015034 (2018) [arXiv:1801.07726 [hep-ph]].
  * (15) Y. Hochberg, E. Kuflik and H. Murayama, “Twin Higgs model with strongly interacting massive particle dark matter,” Phys. Rev. D 99, no.1, 015005 (2019) [arXiv:1805.09345 [hep-ph]].
  * (16) U. K. Dey, T. N. Maity and T. S. Ray, “Boosting Assisted Annihilation for a Cosmologically Safe MeV Scale Dark Matter,” Phys. Rev. D 99, no.9, 095025 (2019) [arXiv:1812.11418 [hep-ph]].
  * (17) S. M. Choi, H. M. Lee, Y. Mambrini and M. Pierre, “Vector SIMP dark matter with approximate custodial symmetry,” JHEP 07, 049 (2019) [arXiv:1904.04109 [hep-ph]].
  * (18) T. N. Maity and T. S. Ray, “Resonant Assisted Annihilation,” JCAP 11, 033 (2019) [arXiv:1907.08262 [hep-ph]].
  * (19) S. M. Choi, J. Kim, H. M. Lee and B. Zhu, “Connecting between inflation and dark matter in models with gauged Z3 symmetry,” JHEP 06, 135 (2020) [arXiv:2003.11823 [hep-ph]].
  * (20) A. Katz, E. Salvioni and B. Shakya, “Split SIMPs with Decays,” JHEP 10, 049 (2020) [arXiv:2006.15148 [hep-ph]].
  * (21) S. Tulin and H. B. Yu, “Dark Matter Self-interactions and Small Scale Structure,” Phys. Rept. 730, 1-57 (2018) [arXiv:1705.02358 [hep-ph]].
  * (22) S. Baek, P. Ko and W. I. Park, “Singlet Portal Extensions of the Standard Seesaw Models to a Dark Sector with Local Dark Symmetry,” JHEP 07, 013 (2013) [arXiv:1303.4280 [hep-ph]].
  * (23) P. Ko and Y. Tang, “Self-interacting scalar dark matter with local $Z_{3}$ symmetry,” JCAP 05 (2014), 047 [arXiv:1402.6449 [hep-ph]].
  * (24) P. Ko and Y. Omura, “Supersymmetric U(1)B x U(1)L model with leptophilic and leptophobic cold dark matters,” Phys. Lett. B 701, 363-366 (2011) [arXiv:1012.4679 [hep-ph]].
  * (25) A. Drozd, B. Grzadkowski and J. Wudka, “Multi-Scalar-Singlet Extension of the Standard Model - the Case for Dark Matter and an Invisible Higgs Boson,” JHEP 04, 006 (2012) [erratum: JHEP 11, 130 (2014)] [arXiv:1112.2582 [hep-ph]].
  * (26) M. Aoki, M. Duerr, J. Kubo and H. Takano, “Multi-Component Dark Matter Systems and Their Observation Prospects,” Phys. Rev. D 86, 076015 (2012) [arXiv:1207.3318 [hep-ph]].
  * (27) S. Baek, P. Ko and W. I. Park, “Hidden sector monopole, vector dark matter and dark radiation with Higgs portal,” JCAP 10, 067 (2014) [arXiv:1311.1035 [hep-ph]].
  * (28) P. Ko and Y. Tang, “$\nu$$\Lambda$MDM: A model for sterile neutrino and dark matter reconciles cosmological and neutrino oscillation data after BICEP2,” Phys. Lett. B 739, 62-67 (2014) [arXiv:1404.0236 [hep-ph]].
  * (29) L. Bian, T. Li, J. Shu and X. C. Wang, “Two component dark matter with multi-Higgs portals,” JHEP 03, 126 (2015) [arXiv:1412.5443 [hep-ph]].
  * (30) A. Karam and K. Tamvakis, “Dark matter and neutrino masses from a scale-invariant multi-Higgs portal,” Phys. Rev. D 92, no.7, 075010 (2015) [arXiv:1508.03031 [hep-ph]].
  * (31) A. Karam and K. Tamvakis, “Dark Matter from a Classically Scale-Invariant $SU(3)_{X}$,” Phys. Rev. D 94, no.5, 055004 (2016) [arXiv:1607.01001 [hep-ph]].
  * (32) S. Bhattacharya, P. Poulose and P. Ghosh, “Multipartite Interacting Scalar Dark Matter in the light of updated LUX data,” JCAP 04, 043 (2017) [arXiv:1607.08461 [hep-ph]].
  * (33) P. Ko and Y. Tang, “Residual Non-Abelian Dark Matter and Dark Radiation,” Phys. Lett. B 768, 12-17 (2017) [arXiv:1609.02307 [hep-ph]].
  * (34) M. Aoki and T. Toma, “Implications of Two-component Dark Matter Induced by Forbidden Channels and Thermal Freeze-out,” JCAP 01, 042 (2017) [arXiv:1611.06746 [hep-ph]].
  * (35) A. Ahmed, M. Duch, B. Grzadkowski and M. Iglicki, “Multi-Component Dark Matter: the vector and fermion case,” Eur. Phys. J. C 78, no.11, 905 (2018) [arXiv:1710.01853 [hep-ph]].
  * (36) M. Aoki and T. Toma, “Boosted Self-interacting Dark Matter in a Multi-component Dark Matter Model,” JCAP 10, 020 (2018) [arXiv:1806.09154 [hep-ph]].
  * (37) S. Chakraborti and P. Poulose, “Interplay of Scalar and Fermionic Components in a Multi-component Dark Matter Scenario,” Eur. Phys. J. C 79, no.5, 420 (2019) [arXiv:1808.01979 [hep-ph]].
  * (38) A. Poulin and S. Godfrey, “Multicomponent dark matter from a hidden gauged SU(3),” Phys. Rev. D 99, no.7, 076008 (2019) [arXiv:1808.04901 [hep-ph]].
  * (39) S. Yaser Ayazi and A. Mohamadnejad, “Scale-Invariant Two Component Dark Matter,” Eur. Phys. J. C 79, no.2, 140 (2019) [arXiv:1808.08706 [hep-ph]].
  * (40) S. Chakraborti, A. Dutta Banik and R. Islam, “Probing Multicomponent Extension of Inert Doublet Model with a Vector Dark Matter,” Eur. Phys. J. C 79, no.8, 662 (2019) [arXiv:1810.05595 [hep-ph]].
  * (41) S. Bhattacharya, P. Ghosh, A. K. Saha and A. Sil, “Two component dark matter with inert Higgs doublet: neutrino mass, high scale validity and collider searches,” JHEP 03, 090 (2020) [arXiv:1905.12583 [hep-ph]].
  * (42) C. R. Chen, Y. X. Lin, C. S. Nugroho, R. Ramos, Y. L. S. Tsai and T. C. Yuan, “Complex scalar dark matter in the gauged two-Higgs-doublet model,” Phys. Rev. D 101, no.3, 035037 (2020) [arXiv:1910.13138 [hep-ph]].
  * (43) C. E. Yaguna and Ó. Zapata, “Multi-component scalar dark matter from a $Z_{N}$ symmetry: a systematic analysis,” JHEP 03, 109 (2020) [arXiv:1911.05515 [hep-ph]].
  * (44) T. Hur, D. W. Jung, P. Ko and J. Y. Lee, “Electroweak symmetry breaking and cold dark matter from strongly interacting hidden sector,” Phys. Lett. B 696, 262-265 (2011) [arXiv:0709.1218 [hep-ph]].
  * (45) P. Ko, “Electroweak symmetry breaking and cold dark matter from hidden sector technicolor,” Int. J. Mod. Phys. A 23, 3348-3351 (2008) [arXiv:0801.4284 [hep-ph]].
  * (46) P. Ko, “Electroweak symmetry breaking and cold dark matter from hidden sector strong interaction,” AIP Conf. Proc. 1178, no.1, 37-43 (2009)
  * (47) Y. Bai and R. J. Hill, “Weakly Interacting Stable Pions,” Phys. Rev. D 82, 111701 (2010) [arXiv:1005.0008 [hep-ph]].
  * (48) P. Ko, “Electroweak symmetry breaking and cold dark matter from strongly interacting hidden sector,” PoS ICHEP2010, 436 (2010) [arXiv:1012.0103 [hep-ph]].
  * (49) T. Hur and P. Ko, “Scale invariant extension of the standard model with strongly interacting hidden sector,” Phys. Rev. Lett. 106, 141802 (2011) [arXiv:1103.2571 [hep-ph]].
  * (50) Y. Bai and P. Schwaller, “Scale of dark QCD,” Phys. Rev. D 89, no.6, 063522 (2014) [arXiv:1306.4676 [hep-ph]].
  * (51) H. Hatanaka, D. W. Jung and P. Ko, “AdS/QCD approach to the scale-invariant extension of the standard model with a strongly interacting hidden sector,” JHEP 08, 094 (2016) [arXiv:1606.02969 [hep-ph]].
  * (52) L. M. Krauss and F. Wilczek, “Discrete Gauge Symmetry in Continuum Theories,” Phys. Rev. Lett. 62, 1221 (1989).
  * (53) S. Baek, P. Ko and W. I. Park, “Local $Z_{2}$ scalar dark matter model confronting galactic ${\mathrm{G}eV}$-scale $\gamma$-ray,” Phys. Lett. B 747, 255-259 (2015) [arXiv:1407.6588 [hep-ph]].
  * (54) S. Baek, J. Kim and P. Ko, “XENON1T excess in local $Z_{2}$ DM models with light dark sector,” Phys. Lett. B 810 (2020), 135848 [arXiv:2006.16876 [hep-ph]].
  * (55) D. W. Kang, P. Ko and C. T. Lu, “Exploring properties of long-lived particles in inelastic dark matter models at Belle II,” [arXiv:2101.02503 [hep-ph]].
  * (56) P. Ko and Y. Tang, “Galactic center $\gamma$-ray excess in hidden sector DM models with dark gauge symmetries: local $Z_{3}$ symmetry as an example,” JCAP 01 (2015), 023 [arXiv:1407.5492 [hep-ph]].
  * (57) J. Guo, Z. Kang, P. Ko and Y. Orikasa, “Accidental dark matter: Case in the scale invariant local B-L model,” Phys. Rev. D 91, no.11, 115017 (2015) [arXiv:1502.00508 [hep-ph]].
  * (58) P. Ko and Y. Tang, “Semi-annihilating $Z_{3}$ Dark Matter for XENON1T Excess,” [arXiv:2006.15822 [hep-ph]].
  * (59) K. S. Babu, C. F. Kolda and J. March-Russell, “Implications of generalized Z - Z-prime mixing,” Phys. Rev. D 57, 6788-6792 (1998) [arXiv:hep-ph/9710441 [hep-ph]].
  * (60) N. Sabti, J. Alvey, M. Escudero, M. Fairbairn and D. Blas, “Refined Bounds on MeV-scale Thermal Dark Sectors from BBN and the CMB,” JCAP 01, 004 (2020) [arXiv:1910.01649 [hep-ph]].
  * (61) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, “FeynRules 2.0 - A complete toolbox for tree-level phenomenology,” Comput. Phys. Commun. 185, 2250-2300 (2014) [arXiv:1310.1921 [hep-ph]].
  * (62) A. Belyaev, N. D. Christensen and A. Pukhov, “CalcHEP 3.4 for collider physics within and beyond the Standard Model,” Comput. Phys. Commun. 184, 1729-1769 (2013) [arXiv:1207.6082 [hep-ph]].
  * (63) P. Gondolo, J. Hisano and K. Kadota, “The Effect of quark interactions on dark matter kinetic decoupling and the mass of the smallest dark halos,” Phys. Rev. D 86, 083523 (2012) [arXiv:1205.1914 [hep-ph]].
  * (64) E. Kuflik, M. Perelstein, N. R. L. Lorier and Y. D. Tsai, “Elastically Decoupling Dark Matter,” Phys. Rev. Lett. 116, no.22, 221302 (2016) [arXiv:1512.04545 [hep-ph]].
  * (65) N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
  * (66) A. M. Sirunyan et al. [CMS], “Search for invisible decays of a Higgs boson produced through vector boson fusion in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Phys. Lett. B 793, 520-551 (2019) [arXiv:1809.05937 [hep-ex]].
  * (67) R. Essig, T. Volansky and T. T. Yu, “New Constraints and Prospects for sub-GeV Dark Matter Scattering off Electrons in Xenon,” Phys. Rev. D 96, no.4, 043017 (2017) [arXiv:1703.00910 [hep-ph]].
  * (68) E. Aprile et al. [XENON], “Search for Light Dark Matter Interactions Enhanced by the Migdal Effect or Bremsstrahlung in XENON1T,” Phys. Rev. Lett. 123, no.24, 241803 (2019) [arXiv:1907.12771 [hep-ex]].
  * (69) E. Aprile et al. [XENON], “Light Dark Matter Search with Ionization Signals in XENON1T,” Phys. Rev. Lett. 123, no.25, 251801 (2019) [arXiv:1907.11485 [hep-ex]].
  * (70) M. Battaglieri, A. Belloni, A. Chou, P. Cushman, B. Echenard, R. Essig, J. Estrada, J. L. Feng, B. Flaugher and P. J. Fox, et al. “US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report,” [arXiv:1707.04591 [hep-ph]].
  * (71) J. P. Lees et al. [BaBar], “Search for a Dark Photon in $e^{+}e^{-}$ Collisions at BaBar,” Phys. Rev. Lett. 113, no.20, 201801 (2014) [arXiv:1406.2980 [hep-ex]].
  * (72) J. P. Lees et al. [BaBar], “Search for Invisible Decays of a Dark Photon Produced in ${e}^{+}{e}^{-}$ Collisions at BaBar,” Phys. Rev. Lett. 119, no.13, 131804 (2017) [arXiv:1702.03327 [hep-ex]].
  * (73) R. Essig, J. Mardon, M. Papucci, T. Volansky and Y. M. Zhong, “Constraining Light Dark Matter with Low-Energy $e^{+}e^{-}$ Colliders,” JHEP 11, 167 (2013) [arXiv:1309.5084 [hep-ph]].
  * (74) D. Banerjee et al. [NA64], “Improved limits on a hypothetical X(16.7) boson and a dark photon decaying into $e^{+}e^{-}$ pairs,” Phys. Rev. D 101, no.7, 071101 (2020) [arXiv:1912.11389 [hep-ex]].
  * (75) S. N. Gninenko, D. V. Kirpichnikov, M. M. Kirsanov and N. V. Krasnikov, “Combined search for light dark matter with electron and muon beams at NA64,” Phys. Lett. B 796, 117-122 (2019) [arXiv:1903.07899 [hep-ph]].
  * (76) C. Ahdida et al. [SHiP], “Sensitivity of the SHiP experiment to light dark matter,” [arXiv:2010.11057 [hep-ex]].
  * (77) J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu, L. W. Mo, T. A. Nunamaker and P. Rassmann, “Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump,” Phys. Rev. D 38, 3375 (1988).
  * (78) M. Davier and H. Nguyen Ngoc, “An Unambiguous Search for a Light Higgs Boson,” Phys. Lett. B 229, 150-155 (1989).
  * (79) S. Andreas, C. Niebuhr and A. Ringwald, “New Limits on Hidden Photons from Past Electron Beam Dumps,” Phys. Rev. D 86, 095019 (2012) [arXiv:1209.6083 [hep-ph]].
  * (80) L. Marsicano, M. Battaglieri, M. Bondi’, C. D. R. Carvajal, A. Celentano, M. De Napoli, R. De Vita, E. Nardi, M. Raggi and P. Valente, “Dark photon production through positron annihilation in beam-dump experiments,” Phys. Rev. D 98, no.1, 015031 (2018) [arXiv:1802.03794 [hep-ex]].
  * (81) D. N. Spergel and P. J. Steinhardt, “Observational evidence for selfinteracting cold dark matter,” Phys. Rev. Lett. 84, 3760-3763 (2000) [arXiv:astro-ph/9909386 [astro-ph]].
  * (82) O. D. Elbert, J. S. Bullock, S. Garrison-Kimmel, M. Rocha, J. Oñorbe and A. H. G. Peter, “Core formation in dwarf haloes with self-interacting dark matter: no fine-tuning necessary,” Mon. Not. Roy. Astron. Soc. 453, no.1, 29-37 (2015) [arXiv:1412.1477 [astro-ph.GA]].
  * (83) X. Chu, C. Garcia-Cely and H. Murayama, “Velocity Dependence from Resonant Self-Interacting Dark Matter,” Phys. Rev. Lett. 122, no.7, 071103 (2019) [arXiv:1810.04709 [hep-ph]].
  * (84) M. Markevitch, A. H. Gonzalez, D. Clowe, A. Vikhlinin, L. David, W. Forman, C. Jones, S. Murray and W. Tucker, “Direct constraints on the dark matter self-interaction cross-section from the merging galaxy cluster 1E0657-56,” Astrophys. J. 606, 819-824 (2004) [arXiv:astro-ph/0309303 [astro-ph]].
  * (85) D. Clowe, A. Gonzalez and M. Markevitch, “Weak lensing mass reconstruction of the interacting cluster 1E0657-558: Direct evidence for the existence of dark matter,” Astrophys. J. 604, 596-603 (2004) [arXiv:astro-ph/0312273 [astro-ph]].
  * (86) S. W. Randall, M. Markevitch, D. Clowe, A. H. Gonzalez and M. Bradac, “Constraints on the Self-Interaction Cross-Section of Dark Matter from Numerical Simulations of the Merging Galaxy Cluster 1E 0657-56,” Astrophys. J. 679, 1173-1180 (2008) [arXiv:0704.0261 [astro-ph]].
  * (87) M. Rocha, A. H. G. Peter, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, J. Onorbe and L. A. Moustakas, “Cosmological Simulations with Self-Interacting Dark Matter I: Constant Density Cores and Substructure,” Mon. Not. Roy. Astron. Soc. 430, 81-104 (2013) [arXiv:1208.3025 [astro-ph.CO]].
  * (88) A. H. G. Peter, M. Rocha, J. S. Bullock and M. Kaplinghat, “Cosmological Simulations with Self-Interacting Dark Matter II: Halo Shapes vs. Observations,” Mon. Not. Roy. Astron. Soc. 430, 105 (2013) [arXiv:1208.3026 [astro-ph.CO]].
  * (89) K. Harigaya, Y. Nakai and M. Suzuki, Phys. Lett. B 809, 135729 (2020) [arXiv:2006.11938 [hep-ph]].
  * (90) S. Baek, P. Ko and W. I. Park, “The 3.5 keV X-ray line signature from annihilating and decaying dark matter in Weinberg model,” [arXiv:1405.3730 [hep-ph]].
  * (91) P. Ko, T. Matsui and Y. L. Tang, “Dark matter bound state formation in fermionic Z2 DM model with light dark photon and dark Higgs boson,” JHEP 10 (2020), 082 [arXiv:1910.04311 [hep-ph]].