# Sterile Neutrino Portal Dark Matter from Semi-Production

Ang Liu1 Feng-Lan Shao1<EMAIL_ADDRESS>Zhi-Long Han2
<EMAIL_ADDRESS>Yi Jin2,3 Honglei Li2 1School of Physics and Physical
Engineering, Qufu Normal University, Qufu, Shandong 273165, China
2School of Physics and Technology, University of Jinan, Jinan, Shandong
250022, China
3Guangxi Key Laboratory of Nuclear Physics and Nuclear Technology, Guangxi
Normal University, Guilin, Guangxi 541004, China

###### Abstract

In this paper, we study the feeble sterile neutrino portal dark matter under
the $Z_{3}$ symmetry. The dark sector consists of one fermion singlet $\chi$
and one scalar singlet $\chi$, which transforms as $\chi\to
e^{i2\pi/3}\chi,\phi\to e^{i2\pi/3}\phi$ under the $Z_{3}$ symmetry. Regarding
fermion singlet $\chi$ as the dark matter candidate, the new interaction terms
$y_{\chi}\phi\bar{\chi^{c}}\chi$ and $\mu\phi^{3}/2$ could induce various new
production channels. For instance, when $m_{\phi}>2m_{\chi}$, the pair decay
$\phi\to\chi\chi$ could be the dominant channel, rather than the delayed decay
$\phi\to\chi\nu$. Another appealing scenario is when the dark sector is
initially produced through the scattering process as
$NN\to\chi\chi,NN\to\phi\phi,h\nu\to\chi\phi$, then the semi-production
processes $N\chi\to\phi\phi,N\phi\to\phi\chi,N\chi\to\chi\chi$ could lead to
the exponential growth of dark sector abundances. The phenomenology of sterile
neutrino and the cosmological impact of the dark scalar are also considered in
the $Z_{3}$ symmetric model.

##  I Introduction

The standard model (SM) has made great achievements in particle physics since
its establishment, including but not limited to its outstanding interpretation
of the basic composition of matter and successful prediction of the Higgs
particle ATLAS:2012yve ; CMS:2012qbp . However, there are still some phenomena
that can not be explained by SM, e.g., the origin of tiny neutrino masses and
the nature of dark matter (DM). The former is established by the discovery of
neutrino oscillation Kajita:2016cak ; McDonald:2016ixn , which implies that
neutrino masses are below the eV scale. The latter is indicated by a variety
of evidence, such as the galactic rotation curves, galaxy clusters and large-
scale structure of cosmology Bertone:2016nfn .

A natural idea is seeking a common interpretation of these two problems, which
has been researched extensively Krauss:2002px ; Asaka:2005an ; Ma:2006km ;
Aoki:2008av ; Cai:2017jrq . Traditionally, high scale sterile neutrinos $N$
are introduced to explain the tiny neutrino mass through the type-I seesaw
mechanism Minkowski:1977sc ; Mohapatra:1979ia . If assuming sterile neutrino
has keV-scale mass, it can be regarded as a decaying DM candidate
Dodelson:1993je ; Drewes:2016upu ; Datta:2021elq ; Das:2021nqj . However, the
corresponding parameter space is now tightly constrained by X-ray searches
Ng:2019gch . One pathway to avoid such constraints is imposing additional
symmetry to make the sterile neutrino a stable DM Ma:2006km ; Ma:2007gq . Then
the sterile neutrino becomes the mediator of neutrino mass generation
Hirsch:2013ola .

Despite the requirement of large Yukawa coupling and leptogenesis
Davidson:2002qv  favoring high scale sterile neutrinos, the naturalness
problem suggests that sterile neutrinos should be below $10^{7}$ GeV
Vissani:1997ys . On the other hand, phenomenological studies usually assume
that sterile neutrinos are below the TeV scale in order to be detected at
colliders Han:2006ip ; Abdullahi:2022jlv . In this paper, we also consider
electroweak scale sterile neutrino. Another advantage of the low scale sterile
neutrino is mediating the interaction between the dark matter and SM, which
provides new annihilation or production channels of DM Escudero:2016tzx ;
Escudero:2016ksa ; Campos:2017odj ; Batell:2017rol ; Hall:2019rld ;
Coito:2022kif .

Since particle dark matter was proposed, weakly interacting massive particle
(WIMP) is the most popular candidate Feng:2010gw ; Roszkowski:2017nbc ;
Schumann:2019eaa ; Arcadi:2017kky , which is generated through the freeze-out
mechanism. Many experiments are devoted to searching for it through direct or
indirect ways XENON:2018voc ; XENON:2023sxq ; PandaX-II:2017hlx ;
PandaX-4T:2021bab ; Boveia:2018yeb ; Adam:2021rrw ; AMS:2013fma ; Fermi-
LAT:2015att . Unfortunately, there are no concrete particle DM signals that
have been found so far. An alternative candidate is the feebly interacting
massive particle (FIMP)Hall:2009bx ; Bernal:2017kxu , which is produced via
the freeze-in mechanism. The interaction between FIMP and SM particles is so
weak that it cannot reach the thermal equilibrium state. Consequently, it is
produced non-thermally by the decay or annihilation of some particles in the
early universe.

The feeble sterile neutrino portal DM under the simplest $Z_{2}$ symmetry has
been studied in Refs. Bandyopadhyay:2020qpn ; Cheng:2020gut ;
Falkowski:2017uya ; Liu:2020mxj ; Chang:2021ose . In this work, we attempt to
explore the generation of feeble DM via the sterile neutrino portal with the
$Z_{3}$ symmetry. Within the framework of type-I seesaw, the sterile neutrino
$N$ can provide masses for SM neutrinos via the Yukawa interaction
$y_{\nu}\overline{L}\tilde{H}N$, and couples to the dark sector. The dark
sector contains a fermion singlet $\chi$ and a scalar singlet $\phi$, both of
which transform as $\chi\to e^{i2\pi/3}\chi,\phi\to e^{i2\pi/3}\phi$ under the
exact $Z_{3}$ symmetry. Providing the mass hierarchy of dark particles as
$m_{\chi}\textless m_{\phi}$, then the dark fermion $\chi$ becomes a DM
candidate. The scenario with strong self-interaction dark scalar $\phi$ and DM
produced from the delayed decay $\phi\to\chi\nu$ is studied in Ref.
Ghosh:2023ocl . Different from this previous study, we assume that the dark
scalar $\phi$ is also feeble interacting with SM. Then we perform a
comprehensive investigation of freeze-in production of DM for representative
scenarios. The WIMP scenario of sterile neutrino portal DM has also been
studied in Ref. Bandyopadhyay:2022tsf ; Liu:2023kil .

Compared with the $Z_{2}$ symmetry, the new interactions $\mu\phi^{3}$ and
$y_{\chi}\phi\bar{\chi}^{c}\chi$ in this $Z_{3}$ symmetry will lead to new
viable parameter space for DM. Recently, the semi-production of FIMP DM has
been proposed in Refs. Bringmann:2021tjr ; Hryczuk:2021qtz , which can lead to
the exponential growth of DM abundance. Semi-production of sterile neutrino DM
is then discussed in Ref. Bringmann:2022aim . In this paper, we will show that
the exponential growth of DM via semi-production processes as
$N\chi\to\chi\chi$, $N\chi\to\phi\phi$ and $N\phi\to\phi\chi$ is also possible
in the $Z_{3}$ symmetric model.

The structure of this paper is organized as follows. In Sec. II, we briefly
introduce the sterile neutrino portal DM model with the $Z_{3}$ symmetry. The
evolution of feeble DM relic density for some representative scenarios are
described in Sec. III. Then we analyze the constraints from testable
signatures under certain scenarios in Sec. IV. Finally, discussions and
conclusions are presented in Sec. V.

##  II The Model

| $L$ |  $N$ |  $\chi$ |  $H$ |  $\phi$
---|---|---|---|---|---
$SU(2)_{L}$ | 2 | 1 | 1 | 2 | 1
$U(1)_{Y}$ | $-\frac{1}{2}$ | 0 | 0 | $\frac{1}{2}$ | 0
$Z_{3}$ | 1 | 1 | $\omega$ | 1 | $\omega$
Table 1: Relevant particle contents and the corresponding charge assignments
under the $Z_{3}$ symmetry. Here $\omega\equiv e^{i2\pi/3}$.

The sterile neutrino portal DM further extends the SM, which includes the
sterile neutrino $N$ and a dark sector with a scalar singlet $\phi$ and a
Dirac fermion singlet $\chi$. Among them, $\chi$ is assumed to be the FIMP DM
candidate for illustration. The particle contents and the corresponding charge
assignments are listed in Table 1. The exact $Z_{3}$ symmetry is employed to
ensure the stability of DM $\chi$, under which the dark sector fields $\phi$
and $\chi$ transform non-trivially as $\phi\to e^{i2\pi/3}\phi$ and $\chi\to
e^{i2\pi/3}\phi$ respectively. Yet the sterile neutrino $N$ and SM fields
transform trivially under the $Z_{3}$ symmetry. The scalar potential under the
unbroken $Z_{3}$ symmetry is

$\displaystyle V$ $\displaystyle=$
$\displaystyle-\mu_{H}^{2}H^{\dagger}H+\mu_{\phi}^{2}\phi^{\dagger}\phi+\lambda_{H}(H^{\dagger}H)^{2}+\lambda_{\phi}(\phi^{\dagger}\phi)^{2}+\lambda_{H\phi}(H^{\dagger}H)(\phi^{\dagger}\phi)+\left(\frac{\mu}{2}\phi^{3}+h.c.\right),$
(1)

where $H$ is the standard Higgs doublet. For simplicity, all the parameters
are taken to be real. To guarantee the unbroken $Z_{3}$ symmetry,
$\lambda_{\phi}\textgreater 0$ and $\mu_{\phi}\textgreater 0$ must be
satisfied. After the electroweak symmetry breaking, $h$ and $\phi$ can obtain
physical masses,

$\displaystyle
m_{h}^{2}=-2\mu_{H}^{2},m_{\phi}^{2}=\mu_{\phi}^{2}+\frac{\lambda_{H\phi}v^{2}}{2},$
(2)

where $h$ is identical to the 125 GeV SM Higgs boson and $v=246~{}{\rm GeV}$.
The scalar potential is bounded below with the conditions Belanger:2012zr

$\lambda_{H}>0,\quad\lambda_{\phi}>0,\quad\lambda_{H\phi}+2\sqrt{\lambda_{H}\lambda_{\phi}}>0.$
(3)

Meanwhile, the estimation of the lifetime of the desired stable vacuum derives
an upper bound on the trilinear coupling, namely $\mu/m_{\phi}\textless
2\sqrt{\lambda_{\phi}}$ Hektor:2019ote . In the following calculation, we take
$\mu=m_{\phi}$ and $\lambda_{\phi}=1$ to meet the above inequality.

The singlet sterile neutrino $N$ not only provides mass for SM neutrinos
through the type-I seesaw mechanism, but also mediates the interaction between
SM and the DM. The new Yukawa interactions and mass terms can be written as

$-\mathcal{L}_{Y}\supset\left(y_{\nu}\overline{L}\widetilde{H}N+y_{N}\phi\bar{\chi}N+\frac{1}{2}m_{N}\overline{N^{c}}N+{\rm
h.c.}\right)+y_{\chi}\phi\bar{\chi^{c}}\chi+m_{\chi}\bar{\chi}\chi,$ (4)

where $\widetilde{H}=i\sigma_{2}H^{\ast}$. The tiny neutrino mass is generated
via the first item, and can be expressed as

$m_{\nu}=-\frac{v^{2}}{2}y_{\nu}~{}m_{N}^{-1}y_{\nu}^{T}.$ (5)

To obtain sub-eV scale light neutrino mass, the Yukawa coupling
$y_{\nu}\lesssim\mathcal{O}(10^{-6})$ is required with electroweak scale
sterile neutrino $N$. In the following studies, we fix $y_{\nu}=10^{-6}$ for
the benchmark points. The seesaw induced mixing angle between the active and
sterile neutrino is then
$\theta=y_{\nu}v/\sqrt{2}m_{N}\lesssim\mathcal{O}(10^{-6})$.

##  III Relic Density

Figure 1: Typical Feynman diagrams for dark sector generation, which also
appear in the $Z_{2}$ symmetric model.

We consider the fermion singlet $\chi$ as the FIMP DM candidate in this paper.
The dark scalar singlet $\phi$ is also assumed feeble interacting with SM, and
is lighter than the sterile neutrino. Meanwhile, the electroweak scale sterile
neutrino $N$ is always in thermal equilibrium via neutrino oscillation
Li:2022bpp  or additional interactions Bandyopadhyay:2022xlp . The generation
of dark scalar $\phi$ is relatively simple, including the Higgs portal
annihilation ${\rm SM}\to\phi\phi$, the sterile neutrino portal direct decay
$N\to\phi\chi$, scattering process $h\nu\to\chi\phi$, $hN\to\chi\phi$, pair
annihilation $NN\to\phi\phi$ and semi-production $N\chi\to\phi\phi$. As for
fermion DM $\chi$, it can be produced through plenty of processes, such as
direct decay $N\to\phi\chi$, delayed decay $\phi\to\chi\nu$, pair decay
$\phi\to\chi\chi$, pair production $NN\to\chi\chi$, semi-production
$N\chi\to\chi\chi$, $N\phi\to\phi\chi$, conversion processes
$\phi\phi\to\chi\chi$ and so on. In addition to the pair decay
$\phi\to\chi\chi$, the semi-production processes $N\chi\to\chi\chi$,
$N\chi\to\phi\phi$ and $N\phi\to\phi\chi$ are new in this $Z_{3}$ symmetric
model. Typical Feynman diagrams for dark sector generation and conversion are
shown in Figures 1 and 2. For simplicity, we neglect those channels with petty
influences of the relic density of the dark sector, e.g. $h\phi\to\phi\phi$,
$h\phi\to\chi\chi$. The relevant Boltzmann equations describing the evolution
of dark sector abundances are given by:

$\displaystyle\frac{dY_{\phi}}{dz}$ $\displaystyle=$
$\displaystyle\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{h\nu\to\chi\phi}\left(Y_{h}^{{\rm eq}}Y_{\nu}^{{\rm
eq}}-\frac{Y_{h}^{{\rm eq}}Y_{\nu}^{{\rm eq}}}{Y_{\chi}^{{\rm
eq}}Y_{\phi}^{{\rm
eq}}}Y_{\chi}Y_{\phi}\right)+\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{Nh\to\chi\phi}\left(Y_{N}^{{\rm eq}}Y_{h}^{{\rm
eq}}-\frac{Y_{N}^{{\rm eq}}Y_{h}^{{\rm eq}}}{Y_{\chi}^{{\rm eq}}Y_{\phi}^{{\rm
eq}}}Y_{\chi}Y_{\phi}\right)$ (6) $\displaystyle+$
$\displaystyle\frac{k}{z^{2}}\left\langle\sigma v\right\rangle_{{\rm
SM}\to\phi\phi}\left(\left(Y_{{\rm SM}}^{{\rm
eq}}\right)^{2}-\left(\frac{Y_{{\rm SM}}^{{\rm eq}}}{Y_{\phi}^{\rm
eq}}\right)^{2}Y_{\phi}^{2}\right)+\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{NN\to\phi\phi}\left(\left(Y_{N}^{{\rm
eq}}\right)^{2}-\left(\frac{Y_{N}^{{\rm eq}}}{Y_{\phi}^{\rm
eq}}\right)^{2}\\!Y_{\phi}^{2}\right)$ $\displaystyle+$ $\displaystyle
k^{\star}z\tilde{\Gamma}_{N\to\phi\chi}\left(Y_{N}^{\rm eq}-\frac{Y_{N}^{\rm
eq}}{Y_{\phi}^{\rm eq}Y_{\chi}^{\rm
eq}}Y_{\phi}Y_{\chi}\right)+\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{N\chi\to\phi\phi}\left(Y_{N}^{{\rm
eq}}Y_{\chi}-\frac{Y_{N}^{\rm eq}Y_{\chi}^{\rm eq}}{(Y_{\phi}^{\rm
eq})^{2}}Y_{\phi}^{2}\right)$ $\displaystyle-$
$\displaystyle\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{\phi\phi\to\chi\chi}\left(Y_{\phi}^{2}-\left(\frac{Y_{\phi}^{\rm
eq}}{Y_{\chi}^{\rm
eq}}\right)^{2}Y_{\chi}^{2}\right)-k^{\star}z\tilde{\Gamma}_{\phi\to\chi\nu}\left(Y_{\phi}-\frac{Y_{\phi}^{\rm
eq}}{Y_{\chi}^{\rm eq}}Y_{\chi}\right)$ $\displaystyle-$ $\displaystyle
k^{\star}z\tilde{\Gamma}_{\phi\to\chi\chi}\left(Y_{\phi}-\frac{Y_{\phi}^{\rm
eq}}{(Y_{\chi}^{\rm eq})^{2}}Y_{\chi}^{2}\right)$
$\displaystyle\frac{dY_{\chi}}{dz}$ $\displaystyle=$
$\displaystyle\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{h\nu\to\chi\phi}\left(Y_{h}^{{\rm eq}}Y_{\nu}^{{\rm
eq}}-\frac{Y_{h}^{{\rm eq}}Y_{\nu}^{{\rm eq}}}{Y_{\chi}^{{\rm
eq}}Y_{\phi}^{{\rm
eq}}}Y_{\chi}Y_{\phi}\right)+\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{Nh\to\chi\phi}\left(Y_{N}^{{\rm eq}}Y_{h}^{{\rm
eq}}-\frac{Y_{N}^{{\rm eq}}Y_{h}^{{\rm eq}}}{Y_{\chi}^{{\rm eq}}Y_{\phi}^{{\rm
eq}}}Y_{\chi}Y_{\phi}\right)$ (7) $\displaystyle+$
$\displaystyle\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{NN\to\chi\chi}\left((Y_{N}^{{\rm eq}})^{2}-\frac{(Y_{N}^{\rm
eq})^{2}}{(Y_{\chi}^{\rm
eq})^{2}}Y_{\chi}^{2}\right)+\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{N\phi\to\phi\chi}\left(Y_{N}^{{\rm
eq}}Y_{\phi}-\frac{Y_{N}^{\rm eq}}{Y_{\chi}^{\rm eq}}Y_{\phi}Y_{\chi}\right)$
$\displaystyle+$ $\displaystyle\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{N\chi\to\chi\chi}\left(Y_{N}^{{\rm
eq}}Y_{\chi}-\frac{Y_{N}^{\rm eq}}{Y_{\chi}^{\rm
eq}}Y_{\chi}^{2}\right)-\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{N\chi\to\phi\phi}\left(Y_{N}^{{\rm
eq}}Y_{\chi}-\frac{Y_{N}^{\rm eq}Y_{\chi}^{\rm eq}}{(Y_{\phi}^{\rm
eq})^{2}}Y_{\phi}^{2}\right)$ $\displaystyle+$
$\displaystyle\frac{k}{z^{2}}\left\langle\sigma
v\right\rangle_{\phi\phi\to\chi\chi}\left(Y_{\phi}^{2}-\left(\frac{Y_{\phi}^{\rm
eq}}{Y_{\chi}^{\rm
eq}}\right)^{2}Y_{\chi}^{2}\right)+k^{\star}z\tilde{\Gamma}_{N\to\phi\chi}\left(Y_{N}^{\rm
eq}-\frac{Y_{N}^{\rm eq}}{Y_{\phi}^{\rm eq}Y_{\chi}^{\rm
eq}}Y_{\phi}Y_{\chi}\right)$ $\displaystyle+$ $\displaystyle
k^{\star}z\tilde{\Gamma}_{\phi\to\chi\nu}\left(Y_{\phi}-\frac{Y_{\phi}^{\rm
eq}}{Y_{\chi}^{\rm
eq}}Y_{\chi}\right)+2k^{\star}z\tilde{\Gamma}_{\phi\to\chi\chi}\left(Y_{\phi}-\frac{Y_{\phi}^{\rm
eq}}{(Y_{\chi}^{\rm eq})^{2}}Y_{\chi}^{2}\right),$

where we use the definition $z\equiv m_{\chi}/T$, and $T$ is the temperature.
The parameters $k$ and $k^{\star}$ are defined as $k=\sqrt{\pi
g_{\star}/45}m_{\chi}M_{Pl}$ and
$k^{*}=\sqrt{45/4\pi^{3}g_{\star}}M_{Pl}/m_{\chi}^{2}$ respectively, where
$g_{\star}$ is the effective number of degrees of freedom of the relativistic
species and $M_{Pl}=1.2\times 10^{19}$ GeV is the Planck mass. The thermal
decay width $\tilde{\Gamma}_{i}$ is calculated as
$\Gamma_{i}\mathcal{K}_{1}/\mathcal{K}_{2}$ with $\mathcal{K}_{1,2}$ being the
first and second modified Bessel Function of the second kind. The
corresponding decay widths are given by

$\displaystyle\Gamma_{{N\rightarrow}{\chi\phi}}$ $\displaystyle=$
$\displaystyle\frac{y_{N}^{2}}{16\pi{m_{N}}}{\left(\frac{(m_{N}+m_{\chi})^{2}-m_{\phi}^{2}}{m_{N}^{2}}\right)}{\lambda^{1/2}({m_{N}}^{2},{m_{\phi}}^{2},{m_{\chi}}^{2})},$
(8) $\displaystyle\Gamma_{\phi\to\chi\nu}$ $\displaystyle=$
$\displaystyle\frac{y_{N}^{2}y^{2}_{\nu}v^{2}\,m_{\phi}}{16\pi
m^{2}_{N}}{\left(\frac{m_{\phi}^{2}-m_{\chi}^{2}}{m_{\phi}^{2}}\right)^{2}},$
(9) $\displaystyle\Gamma_{\phi\to\chi\chi}$ $\displaystyle=$
$\displaystyle\frac{y_{\chi}^{2}}{4\pi{m_{\phi}^{2}}}\left(m_{\phi}^{2}-4m_{\chi}^{2}\right)^{3/2},$
(10)

where the kinematic function $\lambda(a,b,c)$ is defined as

$\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc.$ (11)

Moreover, the thermal average cross sections $\langle\sigma v\rangle$ are
calculated numerically by micrOMEGAs Belanger:2013oya . For the feeble dark
sector, the above Boltzmann equations are solved with the initial condition
$Y_{\chi}=Y_{\phi}=0$. To avoid possible double counting of generated on-shell
particles in the $s$-channel, we also apply the real intermediate states
subtraction Kolb:1979qa . In the above Boltzmann equations, the dark sector
distribution functions following the equilibrium behavior are assumed. More
precise calculations involving semi-production processes can be found in
Ref.Bhatia:2023yux .

The various production channels for DM $\chi$ in this $Z_{3}$ symmetric model
heavily depend on the masses of the dark sector and sterile neutrino.
Depending on whether the decays $N\to\phi\chi$ and $\phi\to\chi\chi$ are
kinematically allowed, we classify the mass spectrum into four scenarios,
namely (1): $m_{N}>m_{\phi}+m_{\chi}$ with $m_{\phi}<2m_{\chi}$, (2):
$m_{N}>m_{\phi}+m_{\chi}$ with $m_{\phi}>2m_{\chi}$, (3):
$m_{N}<m_{\phi}+m_{\chi}$ with $m_{\phi}<2m_{\chi}$, (4):
$m_{N}<m_{\phi}+m_{\chi}$ with $m_{\phi}>2m_{\chi}$, where for the latter two
scenarios $m_{\phi}<m_{N}$ is also satisfied. Theoretically, there are also
four scenarios when $m_{\phi}>m_{N}$. By replacing the contribution of
$N\to\phi\chi$ with $\phi\to N\chi$, we find that the results for
$m_{\phi}>m_{N}$ scenarios are quite similar to the $m_{\phi}<m_{N}$
scenarios, so we will not repeat the $m_{\phi}>m_{N}$ scenarios in this paper.

In the following study, we additionally calculate the results under the
$Z_{2}$ symmetry for comparison. Specifically, we give priority to considering
benchmark points under the $Z_{3}$ symmetry to meet the Planck observed relic
density $\Omega_{\text{DM}}h^{2}=0.12$ Planck:2018vyg , whereupon use the
parameters occurring under the $Z_{2}$ symmetry at the same time, i.e.
$\\{m_{\chi},m_{\phi},m_{N},y_{N},y_{\nu},\lambda_{H\phi}\\}$, to calculate
the abundances of dark particles. In addition, the mass of DM is fixed as 100
GeV for illustration.

Figure 2: Feynman diagrams for the conversion processes $\phi\phi\to\chi\chi$,
and various semi-production processes $N\chi\to\chi\chi$, $N\chi\to\phi\phi$,
$N\phi\to\phi\chi$.

###  III.1 Scenario 1

Scenario 1 | $m_{\chi}$ |  $m_{\phi}$ |  $m_{N}$ |  $y_{N}$ |  $y_{\chi}$ |  $y_{\nu}$ |  $\lambda_{H\phi}$ |  $\mu$
---|---|---|---|---|---|---|---|---
$a$ | 100 | 150 | 300 | $10^{-12}$ | $10^{-12}$ | $10^{-6}$ | $2.1\times 10^{-11}$ | 150
$b$ | 100 | 150 | 300 | $10^{-12}$ | $5\times 10^{-4}$ | $10^{-6}$ | $2.1\times 10^{-11}$ | 150
$c$ | 100 | 150 | 300 | $2.8\times 10^{-12}$ | $10^{-12}$ | $10^{-6}$ | $10^{-14}$ | 150
$d$ | 100 | 150 | 300 | $2.8\times 10^{-12}$ | $5\times 10^{-4}$ | $10^{-6}$ | $10^{-14}$ | 150
Table 2: The parameter choices for scenario 1, the units of masses involved
are GeV.

In scenario 1, we consider that the direct decay $N\to\phi\chi$ is opened,
while the pair decay $\phi\to\chi\chi$ is prohibited. The production of dark
scalar can be classified into two kinds of process. One is the SM Higgs portal
through the coupling $\lambda_{H\phi}$, and the other one is the sterile
neutrino portal via the coupling $y_{N}$. Meanwhile, the new Yukawa coupling
$y_{N}$ contributes to the conversion processes as shown in Figure 2. To
illustrate the impact of these conditions, we select four sets of parameters
in Table 2. The corresponding evolution of $Y_{\phi}$ and $Y_{\chi}$ are shown
in Figure 3.

In scenario 1 (a), we choose the Higgs portal coupling $\lambda_{H\phi}$ being
much larger than the sterile neutrino portal coupling $y_{N}$. In this way,
the dark scalar $\phi$ is dominantly generated through the process ${\rm
SM}\to\phi\phi$, and the decay channel $N\to\phi\chi$ is subdominant. Due to
relatively tiny $y_{N}$ and $y_{\chi}$, the DM abundance $Y_{\chi}$ from
direct decay $N\to\phi\chi$ is miserly, meanwhile contributions from other
$2\to 2$ scattering processes are also negligible. With the cross section
$\langle\sigma v\rangle_{{\rm SM}\to\phi\phi}\simeq 1.9\times
10^{-44}~{}\rm{cm^{3}/s}$, the Planck observed DM abundance is generated via
${\rm SM}\to\phi\phi$ followed by the delayed decay $\phi\to\chi\nu$. In
Figure 3 (a), we can see that the evolution of $Y_{\chi}$ and $Y_{\phi}$ are
consistent in the $Z_{2}$ and $Z_{3}$ symmetry all the time, thus $R_{\chi}$
equals one invariably. This is because of the same generation pattern for the
dark sector with only $\phi\to\chi\nu$ allowed in this scenario.

In scenario 1 (b), the value of $y_{\chi}$ is increased to $5\times 10^{-4}$
compared with scenario 1 (a), meanwhile, the other parameters are kept the
same. As shown in Figure 2, there are new $s$-channel and $t$-channel
contributions to the conversion process $\phi\phi\to\chi\chi$ under the
$Z_{3}$ symmetry which do not involve the coupling $y_{N}$. Different from the
$\phi\phi\to\chi\chi$ process, the other conversion processes are suppressed
by the smallness of $y_{N}$. The corresponding cross section $\langle\sigma
v\rangle_{\phi\phi\to\chi\chi}=6.8\times 10^{-29}~{}\rm{cm^{3}/s}$ has been
greatly enhanced for this scenario, which causes the transition of dark scalar
$\phi$ into DM $\chi$. The results are shown in the panel (b) of Figure. 3,
where $Y_{\chi}$ is increased by a factor of 2.5 before $\phi$ decays compared
with the $Z_{2}$ case. According to our calculation, the conversion becomes
significant when $y_{\chi}\gtrsim 10^{-4}$, i.e., the cross section
$\langle\sigma v\rangle_{\phi\phi\to\chi\chi}\gtrsim 2.7\times
10^{-30}~{}\rm{cm^{3}/s}$. The conversion effect leads to the production of DM
$\chi$ earlier than the $Z_{2}$ case. The ratio $R_{\chi}$ remains on a
downward trend until it becomes a constant after $\phi$ totally freeze-in. The
value of this constant is proportional to the conversion rate $\langle\sigma
v\rangle_{\phi\phi\to\chi\chi}$. In this scenario, the dark scalar $\phi$ is
mainly produced via the process ${\rm SM}\to\phi\phi$ as in scenario 1 (a), so
the same amount of abundance $Y_{\phi}$ is expected provided the absence of
conversion $\phi\phi\to\chi\chi$, which leads to a final reduction of
$R_{\chi}$ to one after the scalar decays via $\phi\to\chi\nu$.

Figure 3: The evolution of dark sector abundances $Y_{\chi}$ (green),
$Y_{\phi}$ (blue) and the ratio $R_{\chi}$ (red) in scenario 1. The solid
lines represent the evolution of the dark sector under $Z_{3}$ symmetry, while
the dashed lines are for the $Z_{2}$ symmetry. The ratio $R_{\chi}$ equals
$Y_{\chi}^{(Z_{3})}/Y_{\chi}^{(Z_{2})}$, where $Y_{\chi}^{(Z_{3})}$ and
$Y_{\chi}^{(Z_{2})}$ are the abundance of DM $\chi$ under the $Z_{3}$ and the
$Z_{2}$ symmetry respectively. The orange dot-dashed lines are the Planck
observed relic density for $m_{\chi}=100$ GeV.

In scenario 1 (c), we consider the opposite case with $\lambda_{H\phi}\ll
y_{N}$. For $\lambda_{H\phi}=10^{-14}$, the Higgs portal process ${\rm
SM}\to\phi\phi$ is heavily suppressed, so as the other $2\to 2$ scattering
processes with $y_{N}\sim y_{\chi}\sim 10^{-12}$. The direct decay
$N\to\phi\chi$ becomes the dominant contribution of $Y_{\phi}$ and $Y_{\chi}$,
which leads to $Y_{\phi}=Y_{\chi}$ at the beginning. The final abundance of
dark scalar is then converted into DM via the delayed decay $\phi\to\chi\nu$.
In this scenario, the ratio $R_{\chi}$ equals to one all the time as shown in
Figure 3 (c).

In scenario 1 (d), the conversion process $\phi\phi\to\chi\chi$ is also
enhanced with relatively large $y_{\chi}$. Since the abundance $Y_{\phi}$
already equals $Y_{\chi}$ from $N\to\phi\chi$ as in scenario 1 (c), the strong
conversion process does not affect the evolution of the dark sector.

Based on the above results, we can conclude that when the direct decay
$N\to\phi\chi$ is allowed and the delayed decay $\phi\to\chi\nu$ is the only
decay mode of dark scalar, the final DM abundance in the $Z_{3}$ symmetric
model is the same as in the $Z_{2}$ symmetric model, although the conversion
process $\phi\phi\to\chi\chi$ could impact the evolution of DM. So in scenario
1, we can not distinguish the $Z_{3}$ symmetry from the $Z_{2}$ symmetry.

###  III.2 Scenario 2

Scenario 2 | $m_{\chi}$ |  $m_{\phi}$ |  $m_{N}$ |  $y_{N}$ |  $y_{\chi}$ |  $y_{\nu}$ |  $\lambda_{H\phi}$ |  $\mu$
---|---|---|---|---|---|---|---|---
$a$ | 100 | 250 | 400 | $10^{-12}$ | $10^{-12}$ | $10^{-6}$ | $2.0\times 10^{-11}$ | 250
$b$ | 100 | 250 | 400 | $10^{-12}$ | $5\times 10^{-4}$ | $10^{-6}$ | $1.9\times 10^{-11}$ | 250
$c$ | 100 | 250 | 400 | $3.7\times 10^{-12}$ | $10^{-12}$ | $10^{-6}$ | $10^{-14}$ | 250
$d$ | 100 | 250 | 400 | $3.7\times 10^{-12}$ | $5\times 10^{-4}$ | $10^{-6}$ | $10^{-14}$ | 250
Table 3: The parameter choices for the four cases in scenario 2, the units of
masses involved are GeV.

For scenario 2, we increase the mass of the dark scalar to open the pair decay
$\phi\to\chi\chi$, while keeping the decay of $N\to\phi\chi$ allowed. Because
the delayed decay $\phi\to\chi\nu$ is further suppressed by the small mixing
angle $\theta$, the pair decay $\phi\to\chi\chi$ is the dominant mode even
with $y_{N}\simeq y_{\chi}$. Four sets of parameters are chosen in Table 3.
Although the generation mode of dark scalar $\phi$ in scenario 2 is consistent
with the corresponding cases in scenario 1, the final conversion of
$\phi\to\chi$ is significantly different. Figure. 4 shows the corresponding
evolution of dark particles.

Figure 4: Same as Figure. 3, but for scenario 2.

In scenario 2 (a), the contributions from direct decay $N\to\phi\chi$ to the
dark sector abundances are tiny. The dark scalar $\phi$ is dominantly produced
from ${\rm SM}\to\phi\phi$. Correct abundance $Y_{\chi}$ is obtained with
$\langle\sigma v\rangle_{{\rm SM}\to\phi\phi}\simeq 6.3\times
10^{-45}~{}\rm{cm^{3}/s}$ followed by the pair decay $\phi\to\chi\chi$. The
conversion of $\phi\to\chi$ happens much earlier than the $Z_{2}$ symmetric
model due to $\Gamma_{\phi\to\chi\chi}\gg\Gamma_{\phi\to\chi\nu}$. The ratio
$R_{\chi}$ equals one before $\phi$ decays, and quickly increases to 41.6
after $\phi$ decays. Since this pair decay converts one $\phi$ into two
$\chi$, the observed DM abundance $Y_{\chi}^{\text{obs}}$ is realized with
$Y_{\phi}(z=10)=Y_{\chi}^{\text{obs}}/2$ in the $Z_{3}$ symmetric model. In
the $Z_{2}$ symmetric model, the conversion is via the delayed decay
$\phi\to\chi\nu$, which leads to
$Y_{\chi}(z=\infty)=Y_{\phi}(z=10)=Y_{\chi}^{\text{obs}}/2$. So the final
ratio $R_{\chi}$ is two in scenario 2 (a).

In scenario 2 (b), the relatively large $y_{\chi}$ not only enhances the
conversion rate of $\phi\phi\to\chi\chi$, but also increases the decay width
$\Gamma_{\phi\to\chi\chi}$. Our numerical calculation finds that compared with
scenario 2 (a), a slightly smaller $\lambda_{H\phi}$ with $\langle\sigma
v\rangle_{{\rm SM}\to\phi\phi}\simeq 5.7\times 10^{-45}~{}\rm{cm^{3}/s}$ could
satisfy the Planck constraint. Once produced, the dark scalar decays quite
quickly into a DM pair, which results in $Y_{\phi}\ll Y_{\chi}$. The inverse
conversion process and the fast pair decay transform a small part of the dark
sector as $2\chi\to 2\phi\stackrel{{\scriptstyle\rm
decay}}{{\longrightarrow}}4\chi$, which makes the generation of DM more
efficient in this scenario. The ratio $R_{\chi}$ decreases during the
evolution, and finally $R_{\chi}$ reaches about 2.1 in scenario 2 (b).

In scenario 2 (c), the dark sector abundances $Y_{\phi}$ and $Y_{\chi}$ are
initially produced via the direct decay $N\to\phi\chi$. Then the dark scalar
$\phi$ is converted to DM $\chi$ by the pair decay $\phi\to\chi\chi$. The
cascade decay chain is $N\to\phi\chi\to\chi\chi\chi$ in the $Z_{3}$ symmetric
model. Under the $Z_{2}$ symmetry, the decay chain is
$N\to\phi\chi\to\chi\nu\chi$. So as shown in Figure 4 (c), the ratio
$R_{\chi}$ increases to 3 after $\phi$ decays in the $Z_{3}$ symmetric model,
and then decreases to 3/2 after $\phi$ decays in the $Z_{2}$ symmetric model.

In scenario 2 (d), the initial dark sector abundances from $N\to\phi\chi$
decay are much smaller than in scenario 2 (b), so the contribution from the
conversion process $\phi\phi\to\chi\chi$ is too small to make $Y_{\chi}$
exceed obviously even with the same $y_{N}$. Therefore, the increase of
$R_{\chi}$ in the early stage is mainly determined by $\phi\to\chi\chi$. The
final ratio $R_{\chi}$ is also 3/2 in scenario 2 (d).

The new pair decay $\phi\to\chi\chi$ makes the $Z_{3}$ symmetric model
different from the $Z_{2}$ symmetric model. With the same couplings in the
$Z_{3}$ symmetric model, the generated DM abundance in the $Z_{2}$ symmetric
model is always smaller than the observed value. Depending on the dominant
generation process of dark scalar, the ratio $R_{\chi}$ is also different.
When the dark scalar is dominantly produced via the Higgs portal ${\rm
SM}\to\phi\phi$, the final ratio is $R_{\chi}\gtrsim 2$. Meanwhile, if the
dark scalar is generated from direct decay $N\to\phi\chi$, the predicted final
ratio is $R_{\chi}=3/2$. The dark scalar is short-lived in the $Z_{3}$
symmetric model due to the relatively large partial decay width
$\Gamma_{\phi\to\chi\chi}$. Then the tight constraints from cosmology can be
easily satisfied in scenario 2.

###  III.3 Scenario 3

Scenario 3 | $m_{\chi}$ |  $m_{\phi}$ |  $m_{N}$ |  $y_{N}$ |  $y_{\chi}$ |  $y_{\nu}$ |  $\lambda_{H\phi}$ |  $\mu$
---|---|---|---|---|---|---|---|---
$a$ | 100 | 140 | 200 | $10^{-12}$ | $10^{-12}$ | $10^{-6}$ | $2.2\times 10^{-11}$ | 140
$b$ | 100 | 140 | 200 | $10^{-12}$ | $5\times 10^{-4}$ | $10^{-6}$ | $2.2\times 10^{-11}$ | 140
$c$ | 100 | 140 | 200 | $4.5\times 10^{-7}$ | $10^{-12}$ | $10^{-6}$ | $10^{-14}$ | 140
$d$ | 100 | 140 | 200 | $1.5\times 10^{-7}$ | $6.2\times 10^{-1}$ | $10^{-6}$ | $10^{-14}$ | 140
Table 4: The parameter choices for the four cases in scenario 3, the units of
masses involved are GeV.

The sterile neutrino portal coupling $y_{N}$ is at the order of
$\mathcal{O}(10^{-12})$ aiming not to exceed the observed DM relic abundance
from direct decay $N\to\phi\chi$ in the previous two scenarios. In scenario 3,
we consider that both $N\to\phi\chi$ and $\phi\to\chi\chi$ are prohibited
kinematically. Compared to the previous two scenarios, the $2\to 2$ scattering
channels as $NN\to\chi\chi$ and $h\nu\to\chi\phi$ will dominate the production
of $\chi$ at the very beginning in this scenario. Besides the Higgs portal
${\rm SM}\to\phi\phi$ channels, the other scattering processes can also make
considerable contributions to the production of $\phi$. We take four sets of
parameters in Table 4 to illustrate this scenario. In addition, the evolution
of the abundance of dark particles is shown in Figure 5.

Figure 5: Same as Figure. 3, but for scenario 3.

In scenario 3 (a), the dark scalar $\phi$ is dominantly produced via ${\rm
SM}\to\phi\phi$. With $\langle\sigma v\rangle_{{\rm SM}\to\phi\phi}\simeq
2.4\times 10^{-44}~{}\rm{cm^{3}/s}$, correct DM relic abundance $Y_{\chi}$ is
obtained by delayed decay $\phi\to\chi\nu$. It is obvious in Figure 5 (a) that
the contribution from scattering to the generation of DM $\chi$ is much lower
than that from $N\to\phi\chi$ decay. Without the contribution from direct
decay $N\to\phi\chi$ to $Y_{\phi}$, a slightly larger $\lambda_{H\phi}$ is
required compared with scenario 1 (a). The ratio $R_{\chi}$ is invariant to
one due to the same transformation process under the two symmetries.

In scenario 3 (b), the conversion process $\phi\phi\to\chi\chi$ is enhanced,
which becomes the dominant production mode of $\chi$. The large conversion
rate leads $R_{\chi}$ to rise to an enormous value $\sim\mathcal{O}(10^{13})$
in the initial time, and then decrease to one with the completion of
$\phi\to\chi\nu$.

In scenario 3 (c), the contribution of ${\rm SM}\to\phi\phi$ can be ignored
due to tiny $\lambda_{H\phi}$. The dark sector is primarily generated by
scattering processes as $NN\to\chi\chi,NN\to\phi\phi$, $h\nu\to\chi\phi$ at
the very beginning. The typical scattering cross sections are $\langle\sigma
v\rangle_{NN\to\chi\chi}\simeq 1.2\times 10^{-48}~{}\rm{cm^{3}/s}$,
$\langle\sigma v\rangle_{NN\to\phi\phi}\simeq 1.8\times
10^{-48}~{}\rm{cm^{3}/s}$ and $\langle\sigma v\rangle_{\chi\phi\to h\nu}\simeq
2.6\times 10^{-47}~{}\rm{cm^{3}/s}$ for the benchmark point. It can be seen
from Figure. 5 (c) that the generated dark abundances from scattering are two
orders of magnitudes lower than the observed value under the $Z_{2}$ symmetry.
Nevertheless, the new semi-production processes $N\chi\to\phi\phi$ and
$N\phi\to\phi\chi$ are enhanced with $\mu=m_{\phi}$ and $y_{N}=4.5\times
10^{-7}$ under the $Z_{3}$ symmetry, which results in the exponential growth
of dark sector abundances. It is worth mentioning that the assumption of
thermal equilibrium of sterile neutrino is important to realize such
exponential growth Bringmann:2021tjr . For the benchmark point, the DM
abundance $Y_{\chi}$ is much larger than the dark scalar abundance $Y_{\phi}$,
so the contribution from delayed decay $\phi\to\chi\nu$ to the total
$Y_{\chi}$ is not obvious. Naturally, the ratio $R_{\chi}$ exponentially
increases to $R_{\chi}^{\text{max}}\simeq 4.6\times 10^{2}$ until the end of
the semi-production processes. Afterwards $R_{\chi}$ is affected by
$\phi\to\chi\nu$, and finally decreases to $2.3\times 10^{2}$ in scenario 3
(c).

In scenario 3 (d), we reduce the value of $y_{N}$, so $Y_{\chi}$ will
eventually fail to satisfy the observed relic density even with the
enhancement by the semi-production processes $N\chi\to\phi\phi$ and
$N\phi\to\phi\chi$ as in scenario 3 (c). On the other hand, $y_{\chi}$ is
taken as a large value $6.2\times 10^{-1}$, which then increases the third
semi-production processes $N\chi\to\chi\chi$ with $\langle\sigma
v\rangle_{N\chi\to\chi\chi}\simeq 2.6\times 10^{-36}~{}\rm{cm^{3}/s}$. The new
semi-production process $N\chi\to\chi\chi$ will cause additional contribution
to the exponential growth of $Y_{\chi}$ to satisfy the Planck constraint.
Meanwhile, the cross section of conversion process $\phi\phi\to\chi\chi$ is
greatly enhanced to about $1.0\times 10^{-22}~{}\rm{cm^{3}/s}$, which makes an
equal amount of $Y_{\phi}$ and $Y_{\chi}$ when $z\lesssim 1$. Afterward, the
conversion process quickly converts the dark scalar into DM. The ratio
$R_{\chi}$ exponentially increases to $R_{\chi}^{\text{max}}\simeq 4.8\times
10^{3}$, and $R_{\chi}$ finally decreases to $2.4\times 10^{3}$.

The former two cases in scenario 3 indicate that when the DM abundance is
dominant by the delayed decay $\phi\to\chi\nu$, the predicted final DM
abundances of $Z_{2}$ and $Z_{3}$ are the same. However, when DM is primarily
generated through the neutrino portal scattering process $NN\to\chi\chi$ and
$h\nu\to\phi\chi$, the semi-production processes $N\chi\to\phi\phi$,
$N\phi\to\phi\chi$ and $N\chi\to\chi\chi$ could lead to the exponential growth
of the dark sector abundances. The latter two cases in scenario 3 have quite
different predictions between the $Z_{2}$ and $Z_{3}$ symmetric models, thus
are useful to distinguish these two models.

###  III.4 Scenario 4

Scenario 4 | $m_{\chi}$ |  $m_{\phi}$ |  $m_{N}$ |  $y_{N}$ |  $y_{\chi}$ |  $y_{\nu}$ |  $\lambda_{H\phi}$ |  $\mu$
---|---|---|---|---|---|---|---|---
$a$ | 100 | 250 | 300 | $10^{-12}$ | $10^{-12}$ | $10^{-6}$ | $2.1\times 10^{-11}$ | 250
$b$ | 100 | 250 | 300 | $10^{-12}$ | $5\times 10^{-4}$ | $10^{-6}$ | $2.0\times 10^{-11}$ | 250
$c$ | 100 | 250 | 300 | $3.2\times 10^{-7}$ | $10^{-12}$ | $10^{-6}$ | $10^{-14}$ | 250
$d$ | 100 | 250 | 300 | $4.6\times 10^{-7}$ | $5\times 10^{-4}$ | $10^{-6}$ | $10^{-14}$ | 250
Table 5: The parameter choices for the four cases in scenario 4, the units of
masses involved are GeV.

Scenario 4 has also opened the pair decay $\phi\to\chi\chi$ in contrast with
scenario 3. Besides the final decay mode of dark scalar $\phi$, the initial
generation channels of the dark sector in scenario 4 are consistent with that
in scenario 3. Table 5 and Figure 6 correspond to the selection of parameters
and the evolution of dark abundances, respectively.

Figure 6: Same as Figure. 3, but for scenario 4.

In scenario 4 (a), the dark scalar $\phi$ is produced via the Higgs portal
${\rm SM}\to\phi\phi$ process. Productions from $2\to 2$ scattering processes
are quite inefficient, and the DM $\chi$ is generated by the fast pair decay
$\phi\to\chi\chi$ under the $Z_{3}$ symmetry. Compared with scenario 3 (a), a
slightly smaller $\lambda_{H\phi}$ is enough to realize the correct DM relic
abundance, which is also due to the pair decay. This decay can lead to the
ratio $R_{\chi}$ increasing to $\mathcal{O}(10^{12})$, and then decreasing to
two finally.

In scenario 4 (b), both the conversion process $\phi\phi\to\chi\chi$ and decay
$\phi\to\chi\chi$ are greatly enhanced. Same as in scenario 2 (b), these two
processes lead to more efficient production of DM than scenario 4 (a), so a
smaller $\lambda_{H\phi}$ in this scenario is enough to produce correct DM
abundance. The ratio $R_{\chi}$ quickly reaches the maximum value of $\sim
10^{14}$, then gradually decreases to 2.2.

In scenario 4 (c), the dark sector abundances are firstly generated by the
$2\to 2$ scattering processes with typical cross section $\langle\sigma
v\rangle_{NN\to\phi\phi}\simeq 2.2\times 10^{-49}~{}\rm{cm^{3}/s}$,
$\langle\sigma v\rangle_{NN\to\chi\chi}\simeq 7.4\times
10^{-50}~{}\rm{cm^{3}/s}$ and $\langle\sigma v\rangle_{\chi\phi\to h\nu}\simeq
1.0\times 10^{-47}~{}\rm{cm^{3}/s}$ for the benchmark point. Then the
relatively large semi-production processes $N\chi\to\phi\phi$ and
$N\phi\to\phi\chi$ exponentially enhance the dark sector abundances. The ratio
$R_{\chi}$ exponentially increases to $7.3\times 10^{2}$, and is further
enlarged by the pair decay $\phi\to\chi\chi$. Finally $R_{\chi}$ decreases to
$5.8\times 10^{2}$ due to the delayed contribution of $\phi\to\chi\nu$ under
the $Z_{2}$ symmetry.

In scenario 4 (d), the large pair decay width $\Gamma_{\phi\to\chi\chi}$ makes
the dark scalar $\phi$ quite short-lived. The produced dark scalar rapidly
decays into the DM pair, rather than taking part in the semi-production
processes $N\chi\to\phi\phi$ and $N\phi\to\phi\chi$, which clearly weakens the
exponential enhancement effect. Therefore, a larger $y_{N}$ is required to
produce the observed DM abundance compared with scenario 4 (c). The ratio
$R_{\chi}$ exponentially increases to $5.8\times 10^{2}$, then decreases to
$2.9\times 10^{2}$ finally.

Similar to scenario 2, the pair decay $\phi\to\chi\chi$ is more efficient in
producing DM abundance in the $Z_{3}$ symmetric model even when the dark
scalar is generated through the Higgs portal ${\rm SM}\to\phi\phi$.
Exponential enhancement by the semi-production processes $N\chi\to\phi\phi$
and $N\phi\to\phi\chi$ are also possible in this scenario. However, the rapid
pair decay $\phi\to\chi\chi$ may weaken the enhancement effect.

##  IV Phenomenology

The sterile neutrino portal FIMP DM model has rich phenomenology Liu:2022cct .
Despite the DM $\chi$ being hard to detect, both the sterile neutrino $N$ and
dark scalar $\phi$ lead to observable signatures. The sterile neutrino $N$ can
be directly produced at colliders Abdullahi:2022jlv . Meanwhile, the neutrino
from delayed decay $\phi\to\chi\nu$ affects the Cosmic Microwave Background
(CMB), the energetic neutrino spectrum and the effective number of
relativistic neutrino species Liu:2022cct .

The collider signatures of sterile neutrino $N$ will be analyzed briefly. The
electroweak scale $N$ can be produced at LHC via the process $pp\to
W\to\ell^{\pm}N$. The cross section of this process is determined by the
mixing angle $\theta$. Lepton number violation signature arises from the decay
$N\to\ell^{\pm}W^{\mp}\to\ell^{\pm}q_{1}\bar{q}_{2}$ Cai:2017mow . When
$m_{N}<m_{W}$, the three-body decay via off-shell $W/Z$ is the dominant
channel, which leads to the displaced vertex signature Drewes:2019fou . In
Figure 7 (a), we summarize the status and future prospect of $N$. By searching
for the displaced vertex signature, a quite large part of the parameter space
with $m_{N}<m_{W}$ can be covered in the future. For our benchmark scenarios,
$y_{\nu}=10^{-6}$ corresponds to $\theta^{2}\sim 10^{-12}$ with a natural
seesaw relation, which is clearly beyond the scope of future sensitivity. Of
course, a larger mixing angle is possible by tuning the parametrization
parameters Barman:2022scg .

Figure 7: Status and future prospect of sterile neutrino $N$ (left). The gray
areas have been excluded by current experiments Blondel:2022qqo . The purple,
red, blue and black dashed lines are the future limits from SHiP SHiP:2018xqw
; Gorbunov:2020rjx , CEPC CEPCStudyGroup:2018ghi , LHC Pascoli:2018heg ;
Izaguirre:2015pga  and FCC-hh Antusch:2016ejd , respectively. The pink line
indicates the seesaw predicted limit. Cosmological constraints of dark scalar
$\phi$ (right). In the right panel, the black dotted line represents the
cosmological constraint discussed in Hambye:2021moy  with $m_{\phi}=100~{}{\rm
GeV}$, the red and purple dotted lines represent the two epochs of CMB and
present, respectively. The circle, triangle, star, and diamond represent
scenarios 1 to 4. Meanwhile, the orange, green, blue and gray samples
represent the four cases (a) to (d) for each scenario.

Then we will focus on the cosmological constraints on $\phi\to\chi\nu$ in
different scenarios under the $Z_{3}$ symmetry. The secondary particles
emitted by the neutrino from delayed decay $\phi\to\chi\nu$ have a great
impact on the CMB anisotropies and spectral distortions. In Figure 7 (b), we
show the corresponding cosmological constraints, where the fractional
abundance $f_{\phi}=\Omega_{\phi}/\Omega_{\rm DM}$,
$\varepsilon=(m^{2}_{\phi}-m^{2}_{\chi})/2m^{2}_{\phi}$ denotes the fraction
of the energy of $\phi$ that has been transferred to neutrinos
Blackadder:2014wpa . According to Table 2 and Figure 3, scenario 1 (a) and 1
(b) predict the same results, thus are overlapped in Figure 7 (b). The same is
true for scenario 1 (c), 1 (d) and scenario 3 (a), 3 (b).

In scenario 1, the typical lifetime of dark scalar $\tau_{\phi}$ is about
$10^{11}\sim 10^{12}$ s with the tiny coupling $y_{N}\sim 10^{-12}$. While
scenarios 1 (a) and 1 (b) are excluded by CMB constraint, scenarios 1 (c) and
1 (d) are still marginally allowed. The DM relic density from $N\to\phi\chi$
has limited the coupling $y_{N}\lesssim 10^{-12}$, so the simplest way to
improve scenario 1 is increasing $y_{\nu}$ to about $\mathcal{O}(10^{-5})$. In
scenario 3, we have $\tau_{\phi}\sim 10^{11}$ s for case (a) and (b),
meanwhile $\tau_{\phi}\sim 10^{1}$ s for case (c) and (d), respectively. The
former two cases are also excluded. Different from scenario 1, increasing
$y_{N}$ for scenarios 3 (a) and (b) is also viable, because the DM relic
density from scattering processes only requires $y_{N}\lesssim 10^{-7}$. As
for scenarios 2 and 4, the branching ratio of $\phi\to\chi\nu$ is much smaller
than that of $\phi\to\chi\chi$, which results in only a tiny part of $\phi$
decaying into neutrinos. Therefore, scenarios 2 and 4 can easily satisfy the
cosmological constraints. In the following discussion, scenario 1 and scenario
3 are considered preferentially as $\phi\to\chi\nu$ is the only decay channel.

Figure 8: The predicted neutrino fluxes at present for scenarios 1 and 3. The
yellow and gray dotted lines are the thermal and nuclear solar neutrino flux
Vitagliano:2019yzm , The black squares and purple triangles represent the
diffuse supernova neutrino background (DSNB) flux measured at the KamLAND
KamLAND:2011bnd  and SK Super-Kamiokande:2013ufi , respectively. The red
points are the atmospheric neutrino data from SK Super-Kamiokande:2015qek .
The orange, green, blue, and gray solid lines correspond to case (a), (b), (c)
and (d) for each scenario, while solid and dot-dashed lines correspond to
$y_{\nu}=10^{-6}$ and $y_{\nu}=10^{-5}$ respectively.

The energetic neutrinos generated by the delayed decay of $\phi$ will be
captured by current neutrino experiments. The neutrino flux at present is
calculated as Bandyopadhyay:2020qpn

$\Phi_{\rm cos}\equiv
E_{\nu}\frac{d\varphi}{dE_{\nu}}=\left(\frac{n_{\phi}}{\tau_{\phi}}\right)\left(\frac{e^{-t(x)/\tau_{\phi}}}{H(x)}\right)\theta^{{}^{\prime}}(x),$
(12)

where $E_{\nu}$ is the observed neutrino energy, $d\varphi/dE_{\nu}$ is the
predicted neutrino flux, $n_{\phi}$ is the number density of $\phi$ if it is
stable, $\theta^{{}^{\prime}}(x)$ is the Heaviside theta function. The cosmic
time $t(x)$ at red-shift $1+x$ and the Hubble parameter $H(x)$ in the standard
cosmology are given by

$\displaystyle t(x)$ $\displaystyle\approx$
$\displaystyle\frac{4}{3H_{0}}\left(\frac{\Omega_{\rm r}^{3/2}}{\Omega_{\rm
m}^{2}}\right)\left(1-\left(1-\frac{\Omega_{\rm m}}{2(1+x)\Omega_{\rm
r}}\right)\sqrt{1+\frac{\Omega_{\rm m}}{(1+x)\Omega_{\rm r}}}\right),$ (13)
$\displaystyle H(x)$ $\displaystyle=$ $\displaystyle
H_{0}\sqrt{\Omega_{\Lambda}+(1+x)^{3}\Omega_{\rm m}+(1+x)^{4}\Omega_{\rm r}},$
(14)

where $x=E_{0}/E_{\nu}-1$ with initial energy
$E_{0}=(m_{\phi}^{2}-m_{\chi}^{2})/2m_{\phi}$, the Hubble constant
$H_{0}=100h~{}\rm{km/s/Mpc}$ with $h=0.6727$ Planck:2018vyg . The dark energy,
matter and radiation fractions are $\Omega_{\Lambda}=0.6846,\Omega_{\rm
m}=0.315$ and $\Omega_{\rm r}=9.265\times 10^{-5}$, respectively.

The neutrino fluxes generated in scenarios 1 and 3 are shown in Figure 8. The
predicted neutrino fluxes with $y_{\nu}=10^{-6}$ for both scenarios are
allowed by current observation. However, such small $y_{\nu}$ may not be
favored by CMB constraint, so we also show the results of $y_{\nu}=10^{-5}$. A
larger $y_{\nu}$ leads to the dark scalar decaying earlier, resulting in less
energetic neutrino flux at present.

The neutrinos generated from $\phi\to\chi\nu$ also increase the effective
number of relativistic neutrino species $N_{\rm eff}$, which can be written as

$N_{\rm
eff}=\frac{7}{8}\left(\frac{11}{4}\right)^{4/3}\left(\frac{\rho_{\nu}}{\rho_{\gamma}}\right)=3\left(\frac{11}{4}\right)^{4/3}\left(\frac{T_{\nu}}{T_{\gamma}}\right)^{4},$
(15)

where $\rho_{\nu}$ and $\rho_{\gamma}$ represent the energy densities of light
neutrinos and photons respectively, $T_{\nu}$ and $T_{\gamma}$ are their
corresponding temperatures. By modifying the evolution equations of $T_{\nu}$
and $T_{\gamma}$ in SM EscuderoAbenza:2020cmq ; Escudero:2018mvt , the
corresponding equations that conform to our model are

$\displaystyle\frac{dT_{\gamma}}{dt}$ $\displaystyle=$
$\displaystyle-\frac{4H\rho_{\gamma}+3H(\rho_{e}+p_{e})+\frac{\delta\rho_{\nu_{e}}}{\delta
t}+2\frac{\delta\rho_{\nu_{\mu}}}{\delta t}-\varepsilon\xi_{\rm
EM}\frac{\rho_{\phi}}{\tau_{\phi}}}{\frac{\partial\rho_{\gamma}}{\partial
T_{\gamma}}+\frac{\partial\rho_{e}}{\partial T_{\gamma}}},$ (16)
$\displaystyle\frac{dT_{\nu}}{dt}$ $\displaystyle=$ $\displaystyle-
HT_{\nu}+\frac{\frac{\delta\rho_{\nu_{e}}}{\delta
t}+2\frac{\delta\rho_{\nu_{\mu}}}{\delta t}+\varepsilon(1-\xi_{\rm
EM})\frac{\rho_{\phi}}{\tau_{\phi}}}{3\frac{\partial\rho_{\nu}}{\partial
T_{\nu}}}.$ (17)

where $\rho_{\gamma,e,\nu}$ denote the energy densities of $\gamma$, $e$ and
$\nu$. $\rho_{\phi}$ expresses the energy density of $\phi$ provided it is
stable. $p_{e}$ is the pressure density of $e$. $\xi_{\rm EM}$ represents the
energy fraction that the neutrinos inject into electromagnetic plasma, which
is assumed to be zero for the selection of $m_{\phi}$ in this work
Hambye:2021moy . The neutrino-electron energy density transfer rate
$\delta\rho_{\nu}/\delta t$ is taken from Refs. Escudero:2018mvt ;
EscuderoAbenza:2020cmq . In addition, we don’t distinguish the flavor of
neutrinos here.

Figure 9: The evolution of $\Delta N_{\rm eff}$ for scenarios 1 and 3. The
calculations are started at $T_{\gamma}=T_{\nu}=10$ MeV with the corresponding
initial time $t_{0}=\frac{1}{2H}|_{T=10~{}\rm MeV}$. The purple and red dashed
lines represent the constraints of $\Delta N_{\rm eff}$ from current Planck
Planck:2018vyg  and future CMB S4 Abazajian:2019eic , respectively. The
orange, green, blue, and gray solid lines correspond to case (a), (b), (c) and
(d) for each scenario, while solid and dot-dashed lines correspond to
$y_{\nu}=10^{-6}$ and $y_{\nu}=10^{-5}$ respectively.

The evolution of $\Delta N_{\rm eff}$ for scenario 1 and 3 are shown in Figure
9, here $\Delta N_{\rm eff}\equiv N_{\rm eff}-N^{\rm SM}_{\rm eff}$ with
$N^{\rm SM}_{\rm eff}=3.045$ Mangano:2005cc ; Grohs:2015tfy ; deSalas:2016ztq
. For scenario 1, results with $y_{\nu}=10^{-6}$ are not favored by current
Planck observation. When $y_{\nu}=10^{-5}$, scenarios 1 (a) and 1 (b) predict
$\Delta N_{\text{eff}}\simeq 0.14$, which is allowed by current Planck limit
and can be further confirmed by future CMB S4. Scenario 1 (c) and 1 (d) with
$y_{\nu}=10^{-5}$ predict $\Delta N_{\text{eff}}\simeq 0.03$, which is beyond
the future sensitivity. For scenario 3, cases (a) and (b) with
$y_{\nu}=10^{-6}$ are excluded by Planck, but are still allowed with
$y_{\nu}=10^{-5}$. Case (c) and (d) predict vanishing small $\Delta
N_{\text{eff}}$ even with $y_{\nu}=10^{-6}$, thus are also hard to detect.

##  V Discussion and Conclusion

The feeble sterile neutrino portal DM with $Z_{3}$ symmetry is studied in this
paper. Besides the sterile neutrino $N$, a dark sector with one fermion
singlet $\chi$ and one scalar singlet $\phi$ is also introduced. The dark
sector $\phi$ and $\chi$ are charged under a $Z_{3}$ symmetry. In addition to
the well-studied sterile neutrino portal Yukawa coupling
$y_{N}\phi\bar{\chi}N$ and Higgs portal coupling
$\lambda_{H\phi}(H^{\dagger}H)(\phi^{\dagger}\phi)$ in the $Z_{2}$ symmetric
model, the $Z_{3}$ symmetry further allows the dark sector Yukawa interaction
$y_{\chi}\phi\bar{\chi^{c}}\chi$ and dark scalar self-interaction
$\mu\phi^{3}/2$. Provided the fermion singlet $\chi$ as the FIMP DM candidate,
the latter two terms could generate new production channels for DM in the
$Z_{3}$ symmetric model.

Because various production channels depend on the mass spectrum, we consider
four specific scenarios to illustrate the evolution of the dark sector. We
find that the dominant production and decay mode of dark scalar $\phi$ has a
great effect on the evolution of DM. When the delayed decay $\phi\to\chi\nu$
is the only decay mode of $\phi$, the dark scalar generated from the Higgs
portal process ${\rm SM}\to\phi\phi$ (as in scenario 1 (a), 1 (b), 3 (a), 3
(b)) or from direct decay $N\to\phi\chi$ (as in scenario 1 (c), 1 (d)) will
lead to the same final DM abundance for both $Z_{2}$ and $Z_{3}$ symmetry,
although the conversion process $\phi\phi\to\chi\chi$ could alert the
evolution of DM in the $Z_{3}$ symmetric model. For natural seesaw required
Yukawa coupling $y_{\nu}=10^{-6}$, these involved scenarios with the only
delayed decay $\phi\to\chi\nu$ mode are already excluded by cosmological
constraints. We show that increasing $y_{\nu}=10^{-5}$ is sufficient to
satisfy all current constraints.

When the pair decay $\phi\to\chi\chi$ is kinematically allowed, it becomes the
dominant decay mode of dark scalar, since the delayed decay $\phi\to\chi\nu$
is heavily suppressed by the tiny mixing angle $\theta\sim 10^{-6}$ in our
analysis. This pair decay $\phi\to\chi\chi$ only appears in the $Z_{3}$
symmetric model, thus definitely leads to a difference between the two kinds
of symmetric models. When the dark scalar is dominantly produced from the
Higgs portal process ${\rm SM}\to\phi\phi$ (as in scenario 2 (a), 2 (b), 4
(a), 4 (b)), the final DM abundance in the $Z_{3}$ symmetry is about twice as
large as it in the $Z_{2}$ symmetry. Meanwhile, if the dark scalar is
generated from the direct decay $N\to\phi\chi$ (as in scenario 2 (c), 2 (d)),
the DM relic abundance ratio of $Z_{3}$ symmetry to $Z_{2}$ symmetry is three
to two. In short, the pair decay is more efficient in producing DM. With a
suppressed branching ratio of $\phi\to\chi\nu$, these scenarios are easily to
avoid the cosmological constraints.

The most interesting scenario is when the dark sector is primarily generated
by the scattering processes as $NN\to\chi\chi,NN\to\phi\phi,h\nu\to\chi\phi$
(as in scenario 3 (c), 3 (d), 4 (c), 4 (d)). Then the semi-production process
$N\chi\to\phi\phi,N\phi\to\phi\chi,N\chi\to\chi\chi$ could lead to the
exponential growth of dark sector abundances in the $Z_{3}$ symmetric model.
Compared with the $Z_{2}$ symmetric model, the final DM abundance of such
scenarios could be enhanced by two to three orders of magnitudes. Our
benchmark points also indicate that the generation of DM $\chi$ is much more
efficient than the dark scalar, which results in a tiny fractional abundance
$f_{\phi}$. Meanwhile, the relatively large Yukawa coupling
$y_{N}\sim\mathcal{O}(10^{-7})$ significantly reduces the lifetime of the dark
scalar $\phi$. These two aspects make such scenarios hard to probe via the
cosmological observables, even when $\phi\to\chi\nu$ is the only decay mode.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China
under Grant No. 11975011, 11805081 and 11635009, Natural Science Foundation of
Shandong Province under Grant No. ZR2019QA021 and ZR2022MA056, the Open
Project of Guangxi Key Laboratory of Nuclear Physics and Nuclear Technology
under Grant No. NLK2021-07.

## References

  * (1) G. Aad et al. [ATLAS], Phys. Lett. B 716, 1-29 (2012) [arXiv:1207.7214 [hep-ex]].
  * (2) S. Chatrchyan et al. [CMS], Phys. Lett. B 716, 30-61 (2012) [arXiv:1207.7235 [hep-ex]].
  * (3) T. Kajita, Rev. Mod. Phys. 88, no.3, 030501 (2016)
  * (4) A. B. McDonald, Rev. Mod. Phys. 88, no.3, 030502 (2016)
  * (5) G. Bertone and D. Hooper, Rev. Mod. Phys. 90, no.4, 045002 (2018) [arXiv:1605.04909 [astro-ph.CO]].
  * (6) L. M. Krauss, S. Nasri and M. Trodden, Phys. Rev. D 67, 085002 (2003) [arXiv:hep-ph/0210389 [hep-ph]].
  * (7) T. Asaka, S. Blanche and M. Shaposhnikov, Phys. Lett. B 631, 151-156 (2005) [arXiv:hep-ph/0503065 [hep-ph]].
  * (8) E. Ma, Phys. Rev. D 73, 077301 (2006) [arXiv:hep-ph/0601225 [hep-ph]].
  * (9) M. Aoki, S. Kanemura and O. Seto, Phys. Rev. Lett. 102, 051805 (2009) [arXiv:0807.0361 [hep-ph]].
  * (10) Y. Cai, J. Herrero-García, M. A. Schmidt, A. Vicente and R. R. Volkas, Front. in Phys. 5, 63 (2017) [arXiv:1706.08524 [hep-ph]].
  * (11) P. Minkowski, Phys. Lett. B 67, 421-428 (1977)
  * (12) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980)
  * (13) S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72, 17-20 (1994) [arXiv:hep-ph/9303287 [hep-ph]].
  * (14) M. Drewes, T. Lasserre, A. Merle, S. Mertens, R. Adhikari, M. Agostini, N. A. Ky, T. Araki, M. Archidiacono and M. Bahr, et al. JCAP 01, 025 (2017) [arXiv:1602.04816 [hep-ph]].
  * (15) A. Datta, R. Roshan and A. Sil, Phys. Rev. Lett. 127, no.23, 231801 (2021) [arXiv:2104.02030 [hep-ph]].
  * (16) A. Das, S. Goswami, V. K. N. and T. K. Poddar, [arXiv:2104.13986 [hep-ph]].
  * (17) K. C. Y. Ng, B. M. Roach, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos and D. R. Wik, Phys. Rev. D 99, 083005 (2019) [arXiv:1901.01262 [astro-ph.HE]].
  * (18) E. Ma, Phys. Lett. B 662, 49-52 (2008) [arXiv:0708.3371 [hep-ph]].
  * (19) M. Hirsch, R. A. Lineros, S. Morisi, J. Palacio, N. Rojas and J. W. F. Valle, JHEP 10, 149 (2013) [arXiv:1307.8134 [hep-ph]].
  * (20) S. Davidson and A. Ibarra, Phys. Lett. B 535, 25-32 (2002) [arXiv:hep-ph/0202239 [hep-ph]].
  * (21) F. Vissani, Phys. Rev. D 57, 7027-7030 (1998) [arXiv:hep-ph/9709409 [hep-ph]].
  * (22) T. Han and B. Zhang, Phys. Rev. Lett. 97, 171804 (2006) [arXiv:hep-ph/0604064 [hep-ph]].
  * (23) A. M. Abdullahi, P. B. Alzas, B. Batell, J. Beacham, A. Boyarsky, S. Carbajal, A. Chatterjee, J. I. Crespo-Anadon, F. F. Deppisch and A. De Roeck, et al. J. Phys. G 50, no.2, 020501 (2023) [arXiv:2203.08039 [hep-ph]].
  * (24) M. Escudero, N. Rius and V. Sanz, JHEP 02, 045 (2017) [arXiv:1606.01258 [hep-ph]].
  * (25) M. Escudero, N. Rius and V. Sanz, Eur. Phys. J. C 77, no.6, 397 (2017) [arXiv:1607.02373 [hep-ph]].
  * (26) M. D. Campos, F. S. Queiroz, C. E. Yaguna and C. Weniger, JCAP 07, 016 (2017) [arXiv:1702.06145 [hep-ph]].
  * (27) B. Batell, T. Han and B. Shams Es Haghi, Phys. Rev. D 97, no.9, 095020 (2018) [arXiv:1704.08708 [hep-ph]].
  * (28) E. Hall, T. Konstandin, R. McGehee and H. Murayama, Phys. Rev. D 107, no.5, 055011 (2023) [arXiv:1911.12342 [hep-ph]].
  * (29) L. Coito, C. Faubel, J. Herrero-García, A. Santamaria and A. Titov, JHEP 08, 085 (2022) [arXiv:2203.01946 [hep-ph]].
  * (30) J. L. Feng, Ann. Rev. Astron. Astrophys. 48, 495-545 (2010) [arXiv:1003.0904 [astro-ph.CO]].
  * (31) L. Roszkowski, E. M. Sessolo and S. Trojanowski, Rept. Prog. Phys. 81, no.6, 066201 (2018) [arXiv:1707.06277 [hep-ph]].
  * (32) M. Schumann, J. Phys. G 46, no.10, 103003 (2019) [arXiv:1903.03026 [astro-ph.CO]].
  * (33) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo and F. S. Queiroz, Eur. Phys. J. C 78, no.3, 203 (2018) [arXiv:1703.07364 [hep-ph]].
  * (34) E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018) [arXiv:1805.12562 [astro-ph.CO]].
  * (35) E. Aprile et al. [XENON], [arXiv:2303.14729 [hep-ex]].
  * (36) X. Cui et al. [PandaX-II], Phys. Rev. Lett. 119, no.18, 181302 (2017) [arXiv:1708.06917 [astro-ph.CO]].
  * (37) Y. Meng et al. [PandaX-4T], Phys. Rev. Lett. 127, no.26, 261802 (2021) [arXiv:2107.13438 [hep-ex]].
  * (38) A. Boveia and C. Doglioni, Ann. Rev. Nucl. Part. Sci. 68, 429-459 (2018) [arXiv:1810.12238 [hep-ex]].
  * (39) W. Adam and I. Vivarelli, Int. J. Mod. Phys. A 37, no.02, 2130022 (2022) [arXiv:2111.10180 [hep-ex]].
  * (40) M. Aguilar et al. [AMS], Phys. Rev. Lett. 110, 141102 (2013)
  * (41) M. Ackermann et al. [Fermi-LAT], Phys. Rev. Lett. 115, no.23, 231301 (2015) [arXiv:1503.02641 [astro-ph.HE]].
  * (42) L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West, JHEP 03, 080 (2010) [arXiv:0911.1120 [hep-ph]].
  * (43) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, Int. J. Mod. Phys. A 32, no.27, 1730023 (2017) [arXiv:1706.07442 [hep-ph]].
  * (44) P. Bandyopadhyay, E. J. Chun and R. Mandal, JCAP 08, 019 (2020) [arXiv:2005.13933 [hep-ph]].
  * (45) Y. Cheng and W. Liao, Phys. Lett. B 815, 136118 (2021) [arXiv:2012.01875 [hep-ph]].
  * (46) A. Falkowski, E. Kuflik, N. Levi and T. Volansky, Phys. Rev. D 99, no.1, 015022 (2019) [arXiv:1712.07652 [hep-ph]].
  * (47) A. Liu, Z. L. Han, Y. Jin and F. X. Yang, Phys. Rev. D 101, no.9, 095005 (2020) [arXiv:2001.04085 [hep-ph]].
  * (48) Z. F. Chang, Z. X. Chen, J. S. Xu and Z. L. Han, JCAP 06, 006 (2021) [arXiv:2104.02364 [hep-ph]].
  * (49) D. K. Ghosh, P. Ghosh and S. Jeesun, [arXiv:2301.13754 [hep-ph]].
  * (50) P. Bandyopadhyay, D. Choudhury and D. Sachdeva, Phys. Rev. D 107, no.1, 015020 (2023) [arXiv:2206.05811 [hep-ph]].
  * (51) A. Liu, Z. L. Han, Y. Jin and H. Li, [arXiv:2306.14091 [hep-ph]].
  * (52) T. Bringmann, P. F. Depta, M. Hufnagel, J. T. Ruderman and K. Schmidt-Hoberg, Phys. Rev. Lett. 127, no.19, 19 (2021) [arXiv:2103.16572 [hep-ph]].
  * (53) A. Hryczuk and M. Laletin, JHEP 06, 026 (2021) [arXiv:2104.05684 [hep-ph]].
  * (54) T. Bringmann, P. F. Depta, M. Hufnagel, J. Kersten, J. T. Ruderman and K. Schmidt-Hoberg, Phys. Rev. D 107, no.7, L071702 (2023) [arXiv:2206.10630 [hep-ph]].
  * (55) G. Belanger, K. Kannike, A. Pukhov and M. Raidal, JCAP 01, 022 (2013) [arXiv:1211.1014 [hep-ph]].
  * (56) A. Hektor, A. Hryczuk and K. Kannike, JHEP 03, 204 (2019) [arXiv:1901.08074 [hep-ph]].
  * (57) S. P. Li and X. J. Xu, JCAP 06, 047 (2023) [arXiv:2212.09109 [hep-ph]].
  * (58) P. Bandyopadhyay, M. Mitra, R. Padhan, A. Roy and M. Spannowsky, JHEP 05, 182 (2022) [arXiv:2201.09203 [hep-ph]].
  * (59) G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, Comput. Phys. Commun. 185, 960-985 (2014) [arXiv:1305.0237 [hep-ph]].
  * (60) E. W. Kolb and S. Wolfram, Nucl. Phys. B 172, 224 (1980) [erratum: Nucl. Phys. B 195, 542 (1982)]
  * (61) D. Bhatia, [arXiv:2308.09801 [hep-ph]].
  * (62) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (63) A. Liu, F. L. Shao, Z. L. Han, Y. Jin and H. Li, [arXiv:2212.10043 [hep-ph]].
  * (64) Y. Cai, T. Han, T. Li and R. Ruiz, Front. in Phys. 6, 40 (2018) [arXiv:1711.02180 [hep-ph]].
  * (65) M. Drewes and J. Hajer, JHEP 02, 070 (2020) [arXiv:1903.06100 [hep-ph]].
  * (66) A. Blondel, C. B. Verhaaren, J. Alimena, M. Bauer, P. Azzi, R. Ruiz, M. Neubert, O. Mikulenko, M. Ovchynnikov and M. Drewes, et al. Front. in Phys. 10, 967881 (2022) [arXiv:2203.05502 [hep-ex]].
  * (67) C. Ahdida et al. [SHiP], JHEP 04, 077 (2019) [arXiv:1811.00930 [hep-ph]].
  * (68) D. Gorbunov, I. Krasnov, Y. Kudenko and S. Suvorov, Phys. Lett. B 810, 135817 (2020) [arXiv:2004.07974 [hep-ph]].
  * (69) J. B. Guimarães da Costa et al. [CEPC Study Group], [arXiv:1811.10545 [hep-ex]].
  * (70) S. Pascoli, R. Ruiz and C. Weiland, JHEP 06, 049 (2019) [arXiv:1812.08750 [hep-ph]].
  * (71) E. Izaguirre and B. Shuve, Phys. Rev. D 91, no.9, 093010 (2015) [arXiv:1504.02470 [hep-ph]].
  * (72) S. Antusch, E. Cazzato and O. Fischer, Int. J. Mod. Phys. A 32, no.14, 1750078 (2017) [arXiv:1612.02728 [hep-ph]].
  * (73) B. Barman, P. S. Bhupal Dev and A. Ghoshal, [arXiv:2210.07739 [hep-ph]].
  * (74) T. Hambye, M. Hufnagel and M. Lucca, JCAP 05, no.05, 033 (2022) [arXiv:2112.09137 [hep-ph]].
  * (75) E. Vitagliano, I. Tamborra and G. Raffelt, Rev. Mod. Phys. 92, 45006 (2020) [arXiv:1910.11878 [astro-ph.HE]].
  * (76) A. Gando et al. [KamLAND], Astrophys. J. 745, 193 (2012) [arXiv:1105.3516 [astro-ph.HE]].
  * (77) H. Zhang et al. [Super-Kamiokande], Astropart. Phys. 60, 41-46 (2015) [arXiv:1311.3738 [hep-ex]].
  * (78) E. Richard et al. [Super-Kamiokande], Phys. Rev. D 94, no.5, 052001 (2016) [arXiv:1510.08127 [hep-ex]].
  * (79) G. Blackadder and S. M. Koushiappas, Phys. Rev. D 90, no.10, 103527 (2014) [arXiv:1410.0683 [astro-ph.CO]].
  * (80) K. Abazajian, G. Addison, P. Adshead, Z. Ahmed, S. W. Allen, D. Alonso, M. Alvarez, A. Anderson, K. S. Arnold and C. Baccigalupi, et al. [arXiv:1907.04473 [astro-ph.IM]].
  * (81) M. Escudero Abenza, JCAP 05, 048 (2020) [arXiv:2001.04466 [hep-ph]].
  * (82) M. Escudero, JCAP 02, 007 (2019) [arXiv:1812.05605 [hep-ph]].
  * (83) G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P. D. Serpico, Nucl. Phys. B 729, 221-234 (2005) [arXiv:hep-ph/0506164 [hep-ph]].
  * (84) E. Grohs, G. M. Fuller, C. T. Kishimoto, M. W. Paris and A. Vlasenko, Phys. Rev. D 93, no.8, 083522 (2016) [arXiv:1512.02205 [astro-ph.CO]].
  * (85) P. F. de Salas and S. Pastor, JCAP 07, 051 (2016) [arXiv:1606.06986 [hep-ph]].