# Common Origin of Dark Matter and Leptogenesis in $U(1)_{B-L}$

Ang Liu1 Feng-Lan Shao1<EMAIL_ADDRESS>Zhi-Long Han2
<EMAIL_ADDRESS>Yi Jin2 Honglei Li2<EMAIL_ADDRESS>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

###### Abstract

In this paper, we investigate the common parameter space of dark matter and
leptogenesis in the $U(1)_{B-L}$ symmetry. This model involves a complex
scalar $\phi$, sterile neutrinos $N$, and Majorana dark matter $\chi$, where
only dark matter $\chi$ is charged under the $Z_{2}$ symmetry. Masses of $N$
and $\chi$ are generated via the Yukawa interactions to $\phi$ after breaking
of the $U(1)_{B-L}$ symmetry. TeV scale sterile neutrinos $N$ are responsible
for the generation of baryon asymmetry through the resonance leptogenesis
mechanism. The new particles in the $U(1)_{B-L}$ have a significant impact on
the dilution of $N$, thus on leptogenesis. Meanwhile, the annihilation
processes of dark matter $\chi$ are almost identical to that of $N$, which
indicates that both leptogenesis and dark matter are closely related to
satisfying the observed results simultaneously. Under various theoretical and
experimental constraints, the viable common parameter space of dark matter and
leptogenesis is obtained for both global and local $U(1)_{B-L}$ symmetry.

##  I Introduction

The standard model (SM) is one of the most successful theories, which provides
excellent interpretations for particle physics. However, some sharp questions
cannot be explained by SM. For instance, observations of neutrino oscillations
Super-Kamiokande:1998kpq ; SNO:2002tuh  indicate that neutrinos have tiny
masses under the constraint from cosmology $\sum m_{\nu}<0.12$ eV
Planck:2018vyg . Another issue is that the amount of CP violation in SM is too
small to account for the observed baryon asymmetry Fong:2012buy . Meanwhile,
various cosmological and astrophysical observations require the existence of
particle dark matter (DM) Bertone:2004pz ; Cirelli:2024ssz . These provide
concrete evidence for new physics beyond the SM. One appealing way is seeking
the common origin of these new physics Asaka:2005pn ; Ma:2006fn ;
Falkowski:2011xh ; Falkowski:2017uya ; Hugle:2018qbw ; Chianese:2019epo ;
Liu:2020mxj .

The famous type-I seesaw mechanism Minkowski:1977sc ; Mohapatra:1979ia  can
naturally explain tiny neutrino masses by introducing sterile neutrinos $N$.
Meanwhile, as described by the leptogenesis mechanism Fukugita:1986hr ;
Davidson:2008bu ; Blanchet:2012bk , the lepton asymmetry could be created by
the out-of-equilibrium decays of $N$ with ${\rm CP}$-asymmetry parameter in
the early universe, afterwards which can be converted to the baryon asymmetry
through the $(B+L)$-violating sphaleron transitions Manton:1983nd ;
Klinkhamer:1984di ; Kuzmin:1985mm ; Khlebnikov:1988sr ; Harvey:1990qw . But
the accompanying question is the mass scale of $N$, namely the seesaw scale,
which can be widely distributed between eV and $10^{14}$ GeV from
Mohapatra:2005wg ; Drewes:2013gca . Typically, the standard thermal
leptogenesis with hierarchical sterile neutrinos requires $m_{N}\gtrsim
10^{9}$ GeV Davidson:2002qv , which is not testable at current or even future
colliders. One viable pathway for TeV scale leptogenesis is the resonant
mechanism with nearly degenerate sterile neutrinos Flanz:1996fb ;
Pilaftsis:1997dr ; Pilaftsis:1997jf ; Pilaftsis:2003gt ; Dev:2017wwc , where
the TeV scale sterile neutrino could lead to a rich phenomenon
Deppisch:2015qwa ; Cai:2017mow ; deGouvea:2013zba ; Alekhin:2015byh .

Another drawback of canonical type-I seesaw is the lack of DM candidates. A
keV scale sterile neutrino could serve as a decaying DM Dodelson:1993je ;
Datta:2021elq , but it is tightly constrained by observations Ng:2019gch . A
stable sterile neutrino DM is possible with discrete symmetry Ma:2006km ,
which usually suffers tight constraints from lepton flavor violation
Kubo:2006yx . Instead of being DM, the sterile neutrino is an ideal messenger
between the dark sector and SM Pospelov:2007mp ; Gonzalez-Macias:2016vxy ;
Escudero:2016ksa ; Batell:2017cmf ; Blennow:2019fhy ; Ballett:2019pyw ;
Biswas:2021kio ; Coy:2021sse ; Borah:2021pet ; Fu:2021uoo ; Das:2022oyx ;
Coito:2022kif ; Biswas:2022vkq ; Liu:2023kil , which can lead to naturally
suppressed DM-nucleon scattering cross section Escudero:2016ksa  and
observable indirect detection signature Batell:2017rol . Currently, the most
popular DM candidate is a weakly interacting massive particle (WIMP) generated
through the freeze-out mechanism. Therefore, many experiments are dedicated to
searching for evidence of their existence XENON:2018voc ; XENON:2023sxq ;
PandaX-II:2017hlx ; PandaX-4T:2021bab ; Boveia:2018yeb ; Adam:2021rrw ;
AMS:2013fma ; Fermi-LAT:2015att .

In seeking for the common origin of neutrino mass, baryon asymmetry and DM, a
meaningful thing to consider is what symmetry follows the new physics. One
natural option is the $B-L$ symmetry Marshak:1979fm ; Mohapatra:1980qe . The
phenomenology of WIMP DM and leptogenesis in the framework of $U(1)_{B-L}$
extended type-I seesaw are separately analyzed in Ref. Okada:2010wd ;
Escudero:2016tzx ; Okada:2018ktp  and Ref. Iso:2010mv ; Heeck:2016oda ;
Dev:2017xry , respectively. Since both DM relic density and baryon asymmetry
via leptogenesis are produced around the TeV scale, they can be closely
correlated by the sterile neutrinos and other new particles in $U(1)_{B-L}$.
In this paper, we further explore the common parameter space for DM and
leptogenesis in the type-I seesaw with $U(1)_{B-L}$ symmetry.

The $U(1)_{B-L}$ extended type-I seesaw introduces sterile neutrinos $N$ to
generate tiny neutrino masses. TeV scale $m_{N}$ is required to produce the
baryon asymmetry through the resonant leptogenesis. A dark fermion $\chi$,
which is odd under the exact $Z_{2}$ symmetry, is regarded as the DM
candidate. In the minimal $B-L$ model, the DM $\chi$ is actually one of the
sterile neutrinos Okada:2010wd . A complex scalar $\phi$ breaks the $B-L$
symmetry, whilst brings masses to $N$ and $\chi$. In this framework, the
annihilation processes of $N$ and $\chi$ mainly affect the final baryon
asymmetry and relic density, which are directly related via the conversion
process $NN\leftrightarrow\chi\chi$. We will demonstrate the common parameter
space for DM and leptogenesis, which can simultaneously satisfy the
corresponding observation values under various constraints.

The structure of this paper is organized as follows. In Sec. II, we briefly
introduce the framework of the model and certain constraints. The calculation
of lepton asymmetry and relic density of DM as well as the analysis of some
phenomenological constraints are revealed in Sec. III for the global
$U(1)_{B-L}$ scenario, and in Sec. IV for the local one. Finally, the
conclusion is presented in Sec. V.

##  II The Framework

| $q_{L}$ | $u_{R}$ | $d_{R}$ | $L_{L}$ | $e_{R}$ |  $N$ |  $\chi$ |  $H$ |  $\phi$
---|---|---|---|---|---|---|---|---|---
$SU(2)_{L}$ | 2 | 1 | 1 | 2 | 1 | 1 | 1 | 2 | 1
$U(1)_{Y}$ | $+\frac{1}{6}$ | $+\frac{2}{3}$ | $-\frac{1}{3}$ | $-\frac{1}{2}$ | -1 | 0 | 0 | $\frac{1}{2}$ | $0$
$U(1)_{B-L}$ | $+\frac{1}{3}$ | $+\frac{1}{3}$ | $+\frac{1}{3}$ | -1 | -1 | -1 | -1 | 0 | $2$
$Z_{2}$ | $+$ | $+$ | $+$ | $+$ | $+$ | $+$ | $-$ | $+$ | $+$
Table 1: Relevant particle contents and corresponding charge assignments.

This model extends complex scalar $\phi$, sterile neutrinos $N$, and DM $\chi$
beyond the SM. The particle contents and the corresponding charge assignments
are listed in Table 1. Regarding the descriptions of Refs. Garcia-Cely:2013wda
; Garcia-Cely:2013nin ; Escudero:2016tzx , the scalar $\phi$ will acquire a
non-zero vacuum expectation value (VEV), which spontaneously breaks the
$U(1)_{B-L}$ symmetry. Ultimately, only $\chi$ is odd under the $Z_{2}$
symmetry, which is naturally identified as the stable DM candidate.

The scalar potential 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),$
(1)

where $H$ is the standard Higgs doublet. $\mu_{H}$ and $\mu_{\phi}$ satisfy
the relationship
$\mu_{H}^{2}=\lambda_{H}v_{H}^{2}+\frac{1}{2}\lambda_{H\phi}v_{\phi}^{2}$ and
$\mu_{\phi}^{2}=\lambda_{\phi}v_{\phi}^{2}+\frac{1}{2}\lambda_{H\phi}v_{H}^{2}$
under the condition of the minimization of the scalar potential
Escudero:2016tzx .

The relevant Yukawa interaction can be written as Escudero:2016tzx

$\displaystyle\mathcal{L}$ $\displaystyle\supset$
$\displaystyle-\frac{\lambda_{N}}{\sqrt{2}}\phi\bar{N}P_{R}N-\frac{\lambda_{\chi}}{\sqrt{2}}\phi\bar{\chi}P_{R}\chi-
y_{\nu}\bar{L}_{L}HP_{R}N+h.c..$ (2)

The tiny neutrino mass is generated via the third term in Eq. (2), and can be
expressed as

$m_{\nu}=-\frac{v_{H}^{2}}{2}y_{\nu}~{}m_{N}^{-1}y_{\nu}^{T}.$ (3)

Here, $v_{H}=246~{}{\rm GeV}$ is the VEV of SM Higgs doublet $H$. Following
the Casas-Ibarra parametrization, the Yukawa matrix $y_{\nu}$ can be expressed
as Casas:2001sr

$y_{\nu}=\frac{\sqrt{2}}{v_{H}}\sqrt{\hat{m}_{N}}R\sqrt{\hat{m}_{\nu}}U_{\rm
PMNS},$ (4)

where $\hat{m}_{N}$ and $\hat{m}_{\nu}$ are the diagonalized mass matrix for
sterile and light neutrinos. $U_{\rm PMNS}$ is the neutrino mixing matrix. It
is possible to achieve relatively large components of $y_{\nu}$ by adjusting
the complex angle of the orthogonal matrix $R$ Barman:2022scg . But our focus
is not this, the parameter related to leptogenesis is the effective neutrino
mass $\tilde{m}\equiv v_{H}^{2}(y_{\nu}^{\dagger}y_{\nu})/m_{N}$, which
determines the decay width of $N\to HL$.

The normal thermal leptogenesis requires $m_{N}\gtrsim 10^{9}$ GeV to generate
a sufficient amount of the baryon asymmetry Davidson:2002qv , which is far
beyond the typical WIMP scale. However, the resonant leptogenesis can solve
this conflict by TeV scale $m_{N}$ Flanz:1996fb ; Pilaftsis:1997dr ;
Pilaftsis:1997jf ; Pilaftsis:2003gt ; Dev:2017wwc . The CP-asymmetry parameter
$\varepsilon_{\rm CP}$ can be enhanced by the self-energy corrections for
nearly degenerate sterile neutrinos Iso:2010mv . Furthermore, the condition
$\Delta m_{N21}^{2}=m_{N2}^{2}-m_{N1}^{2}\simeq m_{N2}\Gamma_{N1}\ll
m_{N2}^{2}$ could result in order unity of $\varepsilon_{\rm CP}$. In the
subsequent calculations, $\varepsilon_{\rm CP}$ will be used as an input
parameter.

###  II.1 The Global $U(1)_{B-L}$ Scenario

In the global $U(1)_{B-L}$ scenario, the complex scalar $\phi$ can be defined
as $\phi=(v_{\phi}+\tilde{\rho}+i\eta)/\sqrt{2}$ with the CP-even real scalar
$\tilde{\rho}$ and CP-odd massless real scalar $\eta$ named Majoron, which is
the Goldstone boson of the spontaneous breaking of the global $U(1)_{B-L}$
symmetry Chikashige:1980ui . After the spontaneous symmetry breaking, there
are two CP-even scalars $h$ and $\rho$, whose masses are

$\displaystyle
m_{h}^{2}=2\lambda_{H}v_{H}^{2}\cos^{2}\theta+2\lambda_{\phi}v_{\phi}^{2}\sin^{2}\theta-\lambda_{H\phi}v_{H}v_{\phi}\sin
2\theta,$ (5) $\displaystyle
m_{\rho}^{2}=2\lambda_{H}v_{H}^{2}\sin^{2}\theta+2\lambda_{\phi}v_{\phi}^{2}\cos^{2}\theta+\lambda_{H\phi}v_{H}v_{\phi}\sin
2\theta.$ (6)

Here $h$ is identical to the 125 GeV SM Higgs boson. $\theta$ is the mixing
angle between two real CP-even scalars $h$ and $\rho$. Then quartic couplings
can be expressed as functions of masses and mixing angle,

$\displaystyle\lambda_{H}$ $\displaystyle=$
$\displaystyle\frac{m_{h}^{2}\cos^{2}\theta+m_{\rho}^{2}\sin^{2}\theta}{2v_{H}^{2}},$
(7) $\displaystyle\lambda_{\phi}$ $\displaystyle=$
$\displaystyle\frac{m_{h}^{2}\sin^{2}\theta+m_{\rho}^{2}\cos^{2}\theta}{2v_{\phi}^{2}},$
(8) $\displaystyle\lambda_{H\phi}$ $\displaystyle=$
$\displaystyle\frac{(m_{\rho}^{2}-m_{h}^{2})\sin 2\theta}{2v_{H}v_{\phi}},$
(9)

where the stability of the scalar potential require
$\lambda_{H\phi}^{2}<4\lambda_{H}\lambda_{\phi}$. This condition can be easily
fulfilled as long as $m_{\rho,h}^{2}>0$. The chiral fermions will obtain
Majorana masses after $U(1)_{B-L}$ symmetry breaking. In the global
$U(1)_{B-L}$ symmetry scenario, the Lagrangian of Eq. (2) becomes

$\displaystyle\mathcal{L}_{\rm global}$ $\displaystyle=$
$\displaystyle-\frac{m_{N}}{2}\bar{N}N-\frac{\lambda_{N}}{2}(-h\sin\theta+\rho\cos\theta)\bar{N}N-\frac{\lambda_{N}}{2}i\eta\bar{N}\gamma^{5}N$
(10)
$\displaystyle-\frac{m_{\chi}}{2}\bar{\chi}\chi-\frac{\lambda_{\chi}}{2}(-h\sin\theta+\rho\cos\theta)\bar{\chi}\chi-\frac{\lambda_{\chi}}{2}i\eta\bar{\chi}\gamma^{5}\chi$
$\displaystyle-(y_{\nu}\bar{L}_{L}HP_{R}N+h.c.),$

where $m_{N}=\lambda_{N}v_{\phi}$ and $m_{\chi}=\lambda_{\chi}v_{\phi}$. The
free parameters are
$\left\\{m_{\rho},m_{\eta},m_{N},m_{\chi},\tilde{m},v_{\phi},\theta,\varepsilon_{\rm
CP}\right\\}$ in the global scenario, then all required couplings
$\lambda_{\phi},\lambda_{H\phi},\lambda_{N},\lambda_{\chi}$ can be expressed
as functions of these parameters. Since different couplings have dissimilar
powers in the renormalization group equation Dev:2015vjd , we use the general
perturbation limits as $\lambda_{N,\chi}<\sqrt{4\pi}$ and
$\lambda_{\phi}<4\pi$ Dev:2017xry .

###  II.2 The Local $U(1)_{B-L}$ Scenario

Different from previous studies Garcia-Cely:2013wda ; Garcia-Cely:2013nin ;
Escudero:2016tzx , we also consider the local scenario. In this case, the
massless Goldstone mode is eaten by the gauge boson $Z^{\prime}$, which
acquires a mass $m_{Z^{\prime}}=2g_{B-L}v_{\phi}$. Here $g_{B-L}$ is the gauge
coupling associated with the $U(1)_{B-L}$ gauge group. The relevant
interactions can be rewritten as

$\displaystyle\mathcal{L}_{\rm local}$ $\displaystyle=$
$\displaystyle-\frac{m_{N}}{2}\bar{N}N-\frac{\lambda_{N}}{2}(-h\sin\theta+\rho\cos\theta)\bar{N}N-\frac{g_{N}^{A}}{2}Z^{\prime}_{\mu}\bar{N}\gamma^{\mu}\gamma^{5}N$
(11) $\displaystyle-$
$\displaystyle\frac{m_{\chi}}{2}\bar{\chi}\chi-\frac{\lambda_{\chi}}{2}(-h\sin\theta+\rho\cos\theta)\bar{\chi}\chi-\frac{g_{\chi}^{A}}{2}Z^{\prime}_{\mu}\bar{\chi}\gamma^{\mu}\gamma^{5}\chi$
$\displaystyle-$ $\displaystyle
g_{f}^{V}Z^{\prime}_{\mu}\bar{f}\gamma^{\mu}f-(y_{\nu}\bar{L}_{L}HP_{R}N+h.c.),$

where $g_{i}^{A,V}=g_{B-L}Q_{i}$ with $Q_{i}$ the the $U(1)_{B-L}$ charge of
particle $i$. $f$ denotes the SM fermions. We need to clarify here that the
vectorial coupling $g_{N,\chi}^{V}$ vanishes for Majorana particles $N$ and
$\chi$, as well as the axial coupling $g_{f}^{A}$ with the disappearing
$U(1)_{B-L}$ charge of SM Higgs doublet $H$ Kahlhoefer:2015bea . In the local
scenario, the free input parameters are
$\left\\{m_{\rho},m_{Z^{\prime}},m_{N},m_{\chi},\tilde{m},v_{\phi},\theta,\varepsilon_{\rm
CP}\right\\}$. Slightly different from the global scenario,
$m_{\rho,N,\chi}\leq\sqrt{4\pi}v_{\phi}$ needs to be satisfied to ensure the
perturbative couplings, i.e., $\lambda_{N,\chi}<\sqrt{4\pi}$ and
$\lambda_{\phi}<2\pi$ Kahlhoefer:2015bea ; Duerr:2016tmh .

##  III The Global $U(1)_{B-L}$ Scenario

###  III.1 Relic Density and Leptogenesis

In the following calculation, we assume that the new scalars $\rho$ and $\eta$
invariably stay in thermal equilibrium by coupling to $H$. Additionally, the
CP-odd scalar $\eta$ with $m_{\eta}=0$ is employed throughout this work for
simplicity. The relevant Feynman diagrams are shown in Fig. 1. Additionally,
almost all contributions of processes $NN(\chi\chi)\to\rho\to{\rm SM}{\rm SM}$
with large $m_{N,\chi}$ come from the boson final states of $V(V=W,Z,h)$.
Therefore, we use $NN(\chi\chi)\to VV$ to represent this process below.
Following Refs. Plumacher:1996kc ; Iso:2010mv ; Dev:2017xry , we can write the
relevant Boltzmann equations as

$\displaystyle\frac{dY_{N}}{dz}$ $\displaystyle=$
$\displaystyle-\frac{z}{sH(m_{N})}\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}-1\right)(\gamma_{N\to HL}+2\gamma_{NL\to qt}+4\gamma_{Nt\to qL})$ (12)
$\displaystyle-$
$\displaystyle\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}\right)^{2}-1\right)(2\gamma_{NN\to\rho\rho,\eta\eta}+2\gamma_{NN\to\rho\eta,h\eta}+2\gamma_{NN\to
VV})$ $\displaystyle+$
$\displaystyle\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{\chi}}{Y_{\chi}^{\rm
eq}}\right)^{2}-\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}\right)^{2}\right)2\gamma_{\chi\chi\to NN},$
$\displaystyle\frac{dY_{\chi}}{dz}$ $\displaystyle=$
$\displaystyle-\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{\chi}}{Y_{\chi}^{\rm
eq}}\right)^{2}-1\right)(2\gamma_{\chi\chi\to\rho\rho,\eta\eta}+2\gamma_{\chi\chi\to\rho\eta,h\eta}+2\gamma_{\chi\chi\to
VV})$ (13) $\displaystyle-$
$\displaystyle\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{\chi}}{Y_{\chi}^{\rm
eq}}\right)^{2}-\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}\right)^{2}\right)2\gamma_{\chi\chi\to NN},$
$\displaystyle\frac{dY_{B-L}}{dz}$ $\displaystyle=$
$\displaystyle\frac{z}{sH(m_{N})}\left(\varepsilon_{{\rm
CP}}\left(\frac{Y_{N}}{Y_{N}^{\rm eq}}-1\right)-\frac{Y_{B-L}}{2Y_{L}^{{\rm
eq}}}\right)\gamma_{N\to HL}$ (14) $\displaystyle-$
$\displaystyle\frac{z}{sH(m_{N})}\frac{Y_{B-L}}{Y_{L}^{{\rm
eq}}}\left(\frac{Y_{N}}{Y_{N}^{{\rm eq}}}\gamma_{NL\to qt}+2\gamma_{Nt\to
qL}\right),$

where $z=m_{N}/T$, the entropy density $s=2\pi^{2}g_{s}T^{3}/(45z^{3})$. The
Hubble expansion rate at temperature $T=m_{N}$ is
$H(m_{N})=\sqrt{4\pi^{3}g_{*}/45}m_{N}^{2}/m_{pl}$ with the Planck mass
$m_{pl}=1.22\times 10^{19}~{}{\rm GeV}$. Here we take both the number of
relativistic degrees of freedom for the entropy density $g_{s}$ and energy
density $g_{\star}$ as 106.75.

Figure 1: The dominant Feynman diagrams for $N$ and $\chi$ annihilation
involving new particles in $U(1)_{B-L}$. Among them, subfigure (a)-(c) are
universal for the global and the local scenarios, while the latter needs to
replace $\eta$ with $Z^{\prime}$. Meanwhile, mixing can lead to processes such
as $NN(\chi\chi)\to\rho,h\to{\rm SMSM}$ in both scenarios, yet the $\rho$
mediated one is dominant for TeV scale $N,\chi$ as shown in panel (d). The
$Z^{\prime}$ mediated channel in panel (d) only appears in the local scenario.

In the above Boltzmann equations, the contributions of the annihilation
processes mediated by $\rho$ to the final SM particles are taken into account
with a small mixing angle $\theta\simeq 0.05$, which is favored by DM direct
detection. Moreover, the lepton asymmetry could be converted to the baryon
asymmetry through the sphaleron processes, i.e., $Y_{B}=\frac{28}{79}Y_{B-L}$.
The DM relic density can be written as
$\Omega_{\chi}h^{2}=m_{\chi}s_{0}Y_{\chi}h^{2}/\rho_{c}$ with the entropy
density of today $s_{0}=2891.2~{}\rm cm^{-3}$ and the critical density of the
universe $\rho_{c}=1.0538\times 10^{-5}h^{2}~{}\rm GeV/cm^{3}$.

The abundance at thermally equilibrium of $N,\chi,L$ can be expressed as

$\displaystyle Y_{N}^{{\rm
eq}}=\frac{45z^{2}}{2\pi^{4}g_{s}}\mathcal{K}_{2}(z),~{}Y_{\chi}^{{\rm
eq}}=\frac{45z^{2}m_{\chi}^{2}}{2\pi^{4}g_{s}m_{N}^{2}}\mathcal{K}_{2}\left(\frac{m_{\chi}}{m_{N}}z\right),~{}Y_{L}^{{\rm
eq}}=\frac{3m_{N}^{3}}{2s\pi^{2}}\zeta(3),$ (15)

where $\mathcal{K}_{2}$ and $\zeta$ are the second modified Bessel function of
the second kind and the Riemann zeta function, respectively. The thermalized
interaction rates for decay and scatter processes are Plumacher:1996kc

$\displaystyle\gamma_{1\to 2+3}$ $\displaystyle=$ $\displaystyle n_{1}^{{\rm
eq}}\frac{\mathcal{K}_{1}}{\mathcal{K}_{2}}\Gamma_{1\to 2+3},$ (16)
$\displaystyle\gamma_{1+2\to 3+4}$ $\displaystyle=$
$\displaystyle\frac{T}{64\pi^{2}}\int_{(m_{1}+m_{2})^{2}}^{\infty}ds^{\prime}\hat{\sigma}(1+2\to
3+4)\sqrt{s^{\prime}}\mathcal{K}_{1}\left(\frac{\sqrt{s^{\prime}}}{T}\right)$
(17)

where $n^{{\rm eq}}$ is the number density in thermal equilibrium, $\Gamma$ is
the decay width. $T$, $s^{\prime}$ and $\hat{\sigma}$ represent the
temperature, the squared center-of-mass energy and the reduced cross section,
respectively. The reduced cross sections of
$NN\to\rho\rho,\eta\eta,\rho\eta,h\eta,VV$ and
$\chi\chi\to\rho\rho,\eta\eta,\rho\eta,h\eta,VV,NN$ can be found in Appendix
A.

The decay width of $N\to HL$ can be expressed as

$\displaystyle\Gamma_{N\to HL}$ $\displaystyle=$
$\displaystyle\frac{m_{N}^{2}\tilde{m}}{8\pi v_{H}^{2}}.$ (18)

Here we fix the effective neutrino mass $\tilde{m}\simeq\sqrt{\Delta m_{\rm
atm}^{2}}\simeq 50~{}\rm meV$ throughout this paper, which is in a natural-
seesaw regime. In addition, the specific representation of $\gamma_{NL\to qt}$
and $\gamma_{Nt\to qL}$ could be found in Ref. Plumacher:1996kc , and some
washout scattering processes with tiny contributions have been overlooked
Borah:2021mri , e.g., $HL\leftrightarrow\bar{H}\bar{L}$.

Figure 2: The evolutions of various abundances $Y_{i}$ with $m_{N}=1000~{}{\rm
GeV}$, $m_{\rho}=500~{}{\rm GeV}$, $\varepsilon_{{\rm CP}}=0.1$, and
$\theta=0.05$ in the global scenario. Subfigures (a)-(d) correspond to four
cases with mass ratios $m_{N}/m_{\chi}=0.5,1,1.5$, and 2 respectively. The
solid black and red lines represent the abundances evolution of $N$ and $\chi$
in thermal equilibrium. The dotted lines in purple and orange are the
observations of baryon asymmetry and DM Planck:2018vyg , namely $Y_{B}^{\rm
obs}\simeq 8.7\times 10^{-11}$ and $\Omega_{{\rm DM}}h^{2}\simeq 0.12$. Solid,
dashed, and dot-dashed lines in green and blue represent $Y_{B-L}$ and
$Y_{\chi}$ with $v_{\phi}=1000~{}{\rm GeV},2000~{}{\rm GeV},3000~{}{\rm GeV}$,
respectively.

The evolutions of various abundances for some benchmarks are shown in Fig. 2.
In these four cases, the evolution of $Y_{N}$ is slightly departure from
thermal equilibrium, so $Y_{N}^{{\rm eq}}$ in all panels represent
approximately the evolution of $Y_{N}$. In panel (a) of Fig. 2, $Y_{B-L}$
needs $v_{\phi}$ to be slightly greater than 1000 GeV to satisfy the
observation of baryon asymmetry. But this result is too small for
observational $Y_{\chi}$, which requires $v_{\phi}$ around 2000 GeV. Both
leptogenesis and DM cannot be satisfied simultaneously by changing $v_{\phi}$
alone, because $Y_{\chi}$ and $Y_{B-L}$ are coherently rising with increasing
$v_{\phi}$, which is directly caused by the decreasing of $\lambda_{\chi}$ and
$\lambda_{N}$. The simplest way to realize common origin is varying the CP-
asymmetry parameter $\varepsilon_{{\rm CP}}$, which means that
$\varepsilon_{{\rm CP}}$ should be less than 0.1 with $v_{\phi}\sim 2000$ GeV
in case (a). In the following three subfigures (b)-(d) of Fig. 2, $Y_{B-L}$
does not show significant changes with the decrease of $m_{\chi}$, which means
that the transition process $NN\to\chi\chi$ is not the dominant contribution
for the dilution of $N$. In contrast, $Y_{\chi}$ continuously increases as
$m_{\chi}$ decreases for fixed $v_{\phi}$. Therefore, an alternative pathway
to satisfy both leptogenesis and DM is modifying $m_{\chi}$, e.g., case (c)
with $m_{N}/m_{\chi}=1.5$ and $v_{\phi}\sim 1000$ GeV. For lighter DM as in
case (d), $v_{\phi}<1000$ GeV is required to satisfy the observed relic
density, meanwhile $\varepsilon_{{\rm CP}}$ should be larger than 0.1 for
leptogenesis.

Figure 3: $Y_{i}$ as a function of $m_{\rho}$ at present in the global
scenario. $m_{N}$ is fixed as 1000 GeV, $\varepsilon_{{\rm CP}}=0.1$, and
$\theta=0.05$. Subfigures (a)-(d) correspond to four cases with mass ratios
$m_{N}/m_{\chi}=0.5,1,1.5$, and 2 respectively. The purple, orange, blue, and
green curves are the same as in Fig. 2. The gray area indicates that
$\lambda_{\phi}$ has exceeded the perturbative limit $\lambda_{\phi}<4\pi$ for
$v_{\phi}=1000$ GeV, while the limits for $v_{\phi}=2000$ GeV and 3000 GeV are
outside the scope of $m_{\rho}$.

We then illustrate the impact of $m_{\rho}$ in Fig. 3. The correct $Y_{B-L}$
and $Y_{\chi}$ occur at the resonance region $m_{\rho}\simeq 2m_{N}$ and
$m_{\rho}\simeq 2m_{\chi}$ separately for $v_{\phi}=3000$ GeV and
$m_{N}/m_{\chi}=0.5$, which is clear that leptogenesis and DM can not be
satisfied simultaneously in this case. When $v_{\phi}=2000$ GeV, there are
three solutions for correct $Y_{\chi}$, with two near the resonance region and
the third one lower than the resonance region. By slightly modifying
$\varepsilon_{{\rm CP}}$, it is possible to obtain the correct $Y_{B-L}$ and
$Y_{\chi}$ at the same time. There will be three values of $\rho$ matching the
observed $Y_{B-L}$ when $v_{\phi}=1000$ GeV. However excessive dilution of
$\chi$ caused by the stronger annihilation processes can only successfully
explain the observed $Y_{\chi}$ when $m_{\rho}\gg m_{N,\chi}$, which is not
allowed by perturbation of $\lambda_{\phi}$. Therefore, we mainly consider
$v_{\phi}>1000$ GeV in the following studies. In panels (b), (c), and (d),
$Y_{B-L}$ can be hardly changed due to the tiny impact of $NN\to\chi\chi$ on
the evolution of $N$. Meanwhile, the DM abundance $Y_{\chi}$ increases as
$m_{\chi}$ decreases, also with the resonance region more concentrated. By
comparing these four cases, we prefer the easiest approach to satisfy both
leptogenesis and DM simultaneously by varying $\varepsilon_{{\rm CP}}$,
especially for light $m_{\rho}$ region. Since both $Y_{B-L}$ and $Y_{\chi}$
favor the resonance region, $m_{N}\sim m_{\chi}$ is required to match the
observed values.

###  III.2 Perturbation and CP-asymmetry

Figure 4: Perturbation constraints and exceeding $\varepsilon_{{\rm CP}}$ on
$m_{\chi}-m_{\rho}$ parameter space in the global scenario. (a)-(d) correspond
to four cases with mass ratios $m_{N}/m_{\chi}=0.5,1,1.5$, and 2 respectively.
The blue stars and black crosses express the regions excluded by perturbation
constraints and exceeding $\varepsilon_{{\rm CP}}$, respectively. The
parameter space composed of lime squares is disallowed by the present direct
detection LZ LZ:2022lsv  limit. The remaining allowed samples are colored by
$\varepsilon_{{\rm CP}}$ (yellow to red).

To obtain the allowed parameter space for the common origin of leptogenesis
and DM, we then explore the following scope

$\displaystyle\begin{aligned} m_{\chi}\in[400,10000]~{}{\rm
GeV},m_{\rho}\in[100,10000]~{}{\rm GeV},v_{\phi}\in[1000,10000]~{}{\rm
GeV},\varepsilon_{{\rm CP}}\in[10^{-6},1].\end{aligned}$ (19)

The mixing angle $\theta$ is important for DM direct detection and collider
signatures. For simplicity, we fix $\theta=0.05$ for illustration. Similar to
the above benchmark scenarios, masses of sterile neutrinos are determined by
the ratio $m_{N}/m_{\chi}=0.5,1,1.5$, and 2 for illustration, because the too
large mass difference between $m_{N}$ and $m_{\chi}$ is not favored for the
common origin. Here we do not consider the lighter DM mass $m_{\chi}<400$ GeV,
which may disrupt the decay process $N\to HL$ and even sphaleron transitions
for the scenarios as $m_{N}/m_{\chi}=0.5$. Meanwhile, $v_{\phi}<1000$ GeV
usually leads to the DM relic density under the experimental value due to too
large coupling. During the scan, the abundances of DM and baryon are required
within the $3\sigma$ range of the Planck observed result Planck:2018vyg ,
i.e., $\Omega_{\chi}h^{2}\in[0.117,0.123],~{}Y_{B}\in[8.6,8.8]\times
10^{-11}$.

The first thing to discuss is the perturbation constraints related to the
couplings, namely, $\lambda_{N,\chi}<\sqrt{4\pi}$ and $\lambda_{\phi}<4\pi$
Dev:2017xry . For certain samples, the dilution effects of $Y_{B-L}$ are so
strong that $\varepsilon_{{\rm CP}}>1$ is required, which can not be realized
even by resonance leptogenesis. The scanning results are shown in Fig. 4. As
shown in panel (a), the parameter space with $m_{\rho}\gtrsim 200$ GeV except
for the DM resonance region $m_{\rho}\simeq 2m_{\chi}$ with $\theta=0.05$, is
excluded by direct detection limit of LZ experiment LZ:2022lsv . Furthermore,
samples in the resonance regions with $m_{\rho}\gtrsim 3000$ GeV will be
excluded by LZ due to smaller $v_{\phi}$. It should be noted that the LZ
exclusion region mainly depends on the value of $\theta$, which will satisfy
the current LZ limit when $\theta\lesssim 0.01$. See more detailed discussion
in Sec. III.3. As we fixed $\theta=0.05$ for all four scenarios in Fig. 4, the
LZ exclusion regions are near the same.

Excessive inhibition of $Y_{B-L}$ results in the need for a larger
$\varepsilon_{{\rm CP}}$ to meet the observed values, sometimes even greater
than 1 in the resonance region of $m_{\rho}\simeq m_{\chi}\simeq 2m_{N}$ for
scenario $m_{N}/m_{\chi}=0.5$. In addition, some parameter space with
$m_{\rho}>2m_{\chi}$ will be disallowed by perturbation constraints, namely
$\lambda_{\phi}\geqslant 4\pi$. Certainly, it is inevitable that the resonance
effect of DM will appear around $m_{\rho}\simeq 2m_{\chi}\simeq 4m_{N}$, where
$\varepsilon_{{\rm CP}}$ will exhibit relatively smaller values
$\varepsilon_{{\rm CP}}\lesssim 10^{-3}$. This is because the correct DM relic
density is satisfied with relatively smaller $y_{\chi}$, thus larger
$v_{\phi}$, which leads to a larger value of $Y_{B-L}$ in the non-resonance
region for leptogenesis. So the baryon asymmetry can be met by simply reducing
$\varepsilon_{{\rm CP}}$. In the special scenario with $m_{N}/m_{\chi}=1$, the
overlapping resonance regions of leptogenesis and DM results in
$\varepsilon_{{\rm CP}}\sim\mathcal{O}(0.1)$ to meet the observations.

Besides the resonance region, DM and leptogenesis can also be simultaneously
satisfied for light scalar as $m_{\rho}\lesssim 200$ GeV. Within this light
scalar region, the required $\varepsilon_{{\rm CP}}$ increases for larger
$m_{N}$. In the scenario with $m_{N}/m_{\chi}=0.5$, $\varepsilon_{{\rm
CP}}\lesssim 1.7\times 10^{-1}$ should be satisfied when $m_{\chi}<10^{4}$
GeV. Increasing the mass ratio $m_{N}/m_{\chi}$ leads to larger required
values of $\varepsilon_{{\rm CP}}$ for certain $m_{\chi}$. As shown in panel
(b) of Fig. 4, $m_{\chi}\gtrsim 7400$ GeV is excluded by $\varepsilon_{{\rm
CP}}>1$ due to too large $m_{N}$ in the light scalar region. The upper bound
on $m_{\chi}$ from $\varepsilon_{{\rm CP}}>1$ decreases as $m_{N}/m_{\chi}$
increases, e.g., $m_{\chi}\gtrsim 3400$ GeV for $m_{N}/m_{\chi}=1.5$ and
$m_{\chi}\gtrsim 2000$ GeV for $m_{N}/m_{\chi}=2$ are excluded. In scenario
$m_{N}/m_{\chi}=2$, the perturbation constraint on $\lambda_{N}$ can also
exclude heavy mass regions, but it is clearly less stringent than the limit
from $\varepsilon_{{\rm CP}}>1$.

###  III.3 Collider Signature and Direct Detection

Figure 5: Predictions for the branching ratio of Higgs invisible decay
BR${}_{\text{inv}}$ (a) and spin-independent cross section $\sigma_{\rm SI}$
(b). The dark gray region in panel (a) is excluded by ATLAS search of Higgs
invisible decay ATLAS:2020kdi , and the black dashed line represents the
predicted result of future CEPC Tan:2020ufz . The lime, blue, and purple
curves stand for $v_{\phi}=1000,2000$, and 3000 GeV, respectively. In
subfigure (b), $v_{\phi}$ and $m_{\rho}$ are fixed as 2000 GeV and 1000 GeV.
The black solid and dashed lines represent the present LZ:2022lsv  and future
LZ:2015kxe  LZ limits. The bottom dot-dashed line is the lower limitation on
the discovery of WIMP in direct detection experiments due to the neutrino
backgrounds Billard:2013qya . The lime, blue, and purple curves express
$\theta=0.1,0.05$, and 0.01, respectively.

In this specific model, the SM Higgs $h$ can decay into $\rho\rho$,
$\eta\eta$, $\chi\chi$, $\nu N$, and $NN$ if kinetically allowed. See Ref.
Escudero:2016tzx  for a detailed discussion. On the other hand, the singlet
scalar $\rho$ is also promising at colliders Dev:2017xry ; Robens:2016xkb ,
which depends on the mixing angle $\theta$. In this paper, we have fixed
$\theta=0.05$, which is allowed by current experimental limits Robens:2022cun
.

In the global scenario, one distinct particle is the Majoron $\eta$, which is
invisible at colliders. With the parameter space in Eq. (19), we focus on the
invisible decay $h\to\eta\eta$, since this channel is always allowed. The
corresponding branching ratio has been constrained by the ATLAS experiment as
ATLAS:2020kdi ,

$\rm BR_{\rm
inv}=\frac{\Gamma_{h\to\eta\eta}}{\Gamma_{h\to\eta\eta}+\Gamma_{\rm
SM}}<0.11,$ (20)

where $\Gamma_{\rm SM}=4$ MeV is the standard Higgs width. The theoretical
Higgs invisible decay width into Majoron is given by Garcia-Cely:2013wda ;
Garcia-Cely:2013nin ; Escudero:2016tzx ,

$\Gamma_{h\to\eta\eta}=\frac{\sin^{2}\theta~{}m_{h}^{3}}{32\pi v_{\phi}^{2}}.$
(21)

The elastic scattering of DM $\chi$ on nuclei mediated by $\rho$ and SM Higgs
$h$ through mixing are tightly constrained by the DM direct detection
experiments LZ:2022lsv ; LZ:2015kxe . The spin-independent scattering cross
section is calculated as Garcia-Cely:2013wda ; Garcia-Cely:2013nin ;
Escudero:2016tzx

$\sigma_{\rm SI}=\frac{C^{2}(\sin 2\theta)^{2}m_{n}^{4}m_{\chi}^{4}}{4\pi
v_{H}^{2}v_{\phi}^{2}(m_{n}+m_{\chi})^{2}}\left(\frac{1}{m_{h}^{2}}-\frac{1}{m_{\rho}^{2}}\right)^{2},$
(22)

where $C\simeq 0.3$ is a nucleon matrix element dependent constant
Belanger:2013oya , and $m_{n}$ represents the nucleon mass.

It can be seen from Eq. (21) and Eq. (22) that the Higgs invisible decay width
and elastic scattering cross section sensitively rely on the mixing angle
$\theta$. Theoretical predictions for the branching ratio of Higgs invisible
decay and scattering cross section for benchmark cases are shown in Fig. 5. In
panel (a), the constraint ability of ATLAS on $\theta$ weakens as $v_{\phi}$
increases. Typically, $\theta\lesssim\mathcal{O}(0.1)$ is required to satisfy
the experimental limit. In the future, the CEPC could probe
BR${}_{\text{inv}}\gtrsim 2.6\times 10^{-3}$ Tan:2020ufz , which corresponds
to $\theta\gtrsim 0.02$ for TeV scale $v_{\phi}$. Meanwhile, the direct
detection limit on $\sigma_{\rm SI}$ has already set more stringent
constraints than Higgs invisible decay. As shown in panel (b) of Fig. 5,
$\theta=0.1$ is excluded by the current LZ limit, although it is allowed by
BR${}_{\text{inv}}$. In the following studies, we will fix $\theta=0.05$ for
illustration. A much smaller $\theta\lesssim 0.01$ is easy to satisfy Higgs
invisible and direct detection limits but is also hard to be detected at
future experiments.

Figure 6: Branching ratio of Higgs invisible decay in the global scenario. We
have fixed $\theta=0.05$ in the calculation. The legends and markers of
samples in the four subfigures are consistent with those in Fig. 4. The solid
and dashed lines correspond to the experimental limits on Higgs invisible
decay of current ATLAS ATLAS:2020kdi  and future CEPC Tan:2020ufz .

Fig. 6 displays the theoretical branching ratio of Higgs invisible decay for
the four scenarios with $m_{N}/m_{\chi}=0.5,1,1.5,2$ within the parameter
space in Eq. (19). By fixing $\theta=0.05$, the predicted branching ratio of
Higgs invisible decay is less than $10^{-2}$, which is one order of magnitude
smaller than current ATLAS limitation. The future CEPC experiment could probe
a large part of the allowed samples. Under the direct detection LZ limit,
allowed samples are also separated into two regions. For light scalar $\rho$
in the non-resonance region (mainly yellow samples as shown in Fig. 4), i.e.,
$m_{\rho}\lesssim 200$ GeV, there is a strong correlation between BRinv and DM
mass $m_{\chi}$. Therefore, we may determine the DM mass by precise
measurement of BRinv once light $\rho$ is discovered. In the resonance region
(mainly red samples as shown in Fig. 4), the predicted BRinv is smaller than
those in the non-resonance region for a certain value of $m_{\chi}$. And only
$m_{\chi}\lesssim 1000$ GeV could lead to observable BRinv at CEPC. Although
the theoretical BRinv is in the range of $[2\times 10^{-4},10^{-2}]$ for all
the four scenarios, allowed samples under current constraints are different,
especially for heavy DM in the non-resonance region. For instance, scenario
$m_{N}/m_{\chi}=2$ predicts BR${}_{\rm inv}\in[3\times 10^{-3},10^{-2}]$ in
the non-resonance region, since $m_{\chi}\gtrsim 2000$ GeV is already excluded
by $\varepsilon_{{\rm CP}}>1$ and non-perturbative.

Figure 7: Spin-independent scattering cross section in the global scenario
with $\theta=0.05$. Legends and markers appearing in these subfigures have the
same meaning as that in Fig. 4. The black solid, dashed and dot-dashed lines
represent the present LZ limit LZ:2022lsv , future LZ limit LZ:2015kxe , and
neutrino background Billard:2013qya .

The scanning results of the scattering cross sections are displayed in Fig. 7,
which are tightly constrained by the LZ limit with $\theta=0.05$. For example,
the present LZ limit at most excludes the lime samples with a cross section
larger than $10^{-46}\rm{cm}^{2}$ when $m_{\chi}\sim 400$ GeV. The future LZ
experiment can further push this limit down by one order of magnitude.
Therefore, most samples with $\theta=0.05$ are within future experimental
reach. With relatively heavy $m_{\rho}$, we find that the samples in the
resonance region (mainly red samples) usually predict detectable direct
detection cross section. On the other hand, the samples in the non-resonance
region (mainly yellow samples) could lead to suppressed scattering cross
section due to cancellation for $m_{\rho}\sim m_{h}$ as in Eq. (22).

As shown in Fig. 5, $\theta$ plays a crucial role in the elastic scattering
cross section. For excessively large $\theta$, such as $\theta=0.1$, we have
checked that almost all samples will be excluded by the present LZ limit.
However, the constraints from the direct detection experiment are extremely
relaxed when $\theta$ takes a smaller value like 0.01. When $m_{N}/m_{\chi}$
increase, as shown in panels (b), (c), and (d), an increasing number of
samples that escape the constraints of direct detection experiments will be
excluded by unitarity constraints and exceeding $\varepsilon_{{\rm CP}}$,
although the exclusion range of LZ has not changed much.

###  III.4 Indirect Detection

The annihilation of DM into SM particles provides an alternative pathway to
probe DM. In the global scenario, DM $\chi$ can directly annihilate into the
final states of $NN$, $\rho\rho$, $\rho\eta$, $\eta\eta$, and $VV$, which is
shown in Fig. 1. At present, the process $\chi\chi\to\rho\rho,\eta\eta$ , and
$VV$ are $p$-wave suppressed by DM relative velocity, thus can not be detected
by indirect detection. So we focus on the $s$-wave processes $\chi\chi\to
NN,\rho\eta$. Once kinetically allowed, the $\chi\chi\to\rho\eta$ channel is
the dominant one Escudero:2016tzx , therefore allowed samples in the non-
resonance region with light $m_{\rho}$ are dominant by this channel. On the
other hand, the $\chi\chi\to NN$ channel becomes dominant only in the
resonance region when $m_{\rho}\simeq 2m_{\chi}$. Meanwhile, the $\chi\chi\to
NN$ channel is kinetically forbidden in scenarios (b), (c), and (d) with
$m_{N}\geq m_{\chi}$. In Fig. 8, we show the theoretical $\left<\sigma
v\right>$ and various experimental limits, where $\left<\sigma v\right>$ is
the thermally averaged cross section for DM annihilation times velocity at
present.

Figure 8: Constraints from the indirect detection experiments in the global
scenario. The solid lines stand for the upper limit on $\left<\sigma v\right>$
for pair annihilation of DM into neutrinos by ANTARES Gozzini:2023nka . The
dot-dashed lines in panel (a) express the bound on $NN$ final state from
Fermi-LAT and H.E.S.S. Campos:2017odj . The dashed lines represent the
projected upper limits of $W^{+}W^{-}$ final states from the Cherenkov
Telescope Array (CTA) CTA:2020qlo ; Mangipudi:2021ivm . The definitions of the
blue squares, black crosses, and lime stars are the same as those in Fig. 4.
The remaining allowed samples are colored by $\varepsilon_{{\rm CP}}$ while
shaped by different final products.

In subfigure (a) of Fig. 8, the traditional WIMP $\left<\sigma v\right>$
namely $2\times 10^{-26}~{}\rm cm^{3}/s$ is mostly composed of channel
$\chi\chi\to\rho\eta$. Assuming $\theta=0.05$, we find that $\rho\to\eta\eta$
is the dominant decay mode of $\rho$. The pseudo scalar $\eta$ could have keV
scale mass due to quantum gravity effects Lattanzi:2014mia , which leads to
the cascade decay $\eta\to\nu\nu$ in the global $U(1)_{B-L}$ scenario. In this
way, the process $\chi\chi\to\rho\eta$ will generate neutrinos via cascade
decays of $\rho$ and $\eta$. Currently, the experimental limit on the
annihilation cross section of DM into neutrinos is at the order of
$10^{-24}~{}{\rm cm}^{3}/s$, which has no limitation on theoretical parameter
space. The sub-dominant decay mode $\rho\to W^{+}W^{-}$ has a branching ratio
of 0.1 for $\theta=0.05$. Therefore, $\left<\sigma v\right>$ of
$\chi\chi\to\rho\eta\to W^{+}W^{-}\eta$ is roughly located near $10^{-27}\rm
cm^{3}/s$, which is one order of magnitude smaller than the projected upper
limit for $W^{+}W^{-}$ final state.

Within the resonance region, both $p$-wave $\chi\chi\to\eta\eta,VV$ and
$s$-wave $\chi\chi\to NN$ channels contribute to the relic density of DM. At
present, the cross sections of $p$-wave $\chi\chi\to\eta\eta,VV$ channels are
velocity suppressed, resulting in the $\chi\chi\to NN$ the dominant channel.
So the cross section of $\chi\chi\to NN$ is typically smaller than the WIMP
value $\left<\sigma v\right>\sim 2\times 10^{-26}~{}\rm cm^{3}/s$. With
$\theta=0.05$, we also find that the samples with $\left<\sigma
v\right>\gtrsim 2\times 10^{-28}\rm cm^{3}/s$ are excluded dominantly by LZ.
Therefore, the $\chi\chi\to NN$ channel is also beyond the reach of indirect
detection.

Different from scenario (a), smaller cross section $\left<\sigma v\right>$ in
the resonance region for scenarios (b), (c) and (d) are dominated by the
$s$-wave process $\chi\chi\to h\eta$, which is suppressed by the mixing angle
$\theta=0.05$. So the the surviving samples with $\left<\sigma
v\right>\lesssim\mathcal{O}(10^{-29})~{}\rm cm^{3}/s$ are difficult to detect
by current or future experiments. We also find that there are few samples are
dominant by the $\chi\chi\to W^{+}W^{-}$ process when $v_{\phi}$ is very
large, namely $\lambda_{H\phi}$ is much greater than $\lambda_{\phi}$. For the
non-resonance region, the dominant annihilation mode is $\chi\chi\to\rho\eta$,
which leads to similar results as scenario (a) for indirect detection. It
should be noted that our indirect detection results are different from the
previous study in Ref. Escudero:2016tzx , mainly because of the non-vanishing
mixing $\theta$. In the simplified case with $\theta=0$, the processes into SM
final states as $\chi\chi\to h\eta,W^{+}W^{-}$ mediated by $\rho$ are absent
as in Ref. Escudero:2016tzx , then the $p$-wave $\chi\chi\to\eta\eta$ becomes
the dominant contribution in the resonance region.

##  IV The Local $U(1)_{B-L}$ Scenario

###  IV.1 Relic Density and Leptogenesis

For the local scenario, the various dominant annihilation processes of $N$ and
$\chi$ correspond to Feynman diagrams in Fig.1 by replacing $\eta$ with
$Z^{\prime}$ in panel (a)-(c). Whereupon the Boltzmann equations Eq.12 can be
rewritten as

$\displaystyle\frac{dY_{N}}{dz}$ $\displaystyle=$
$\displaystyle-\frac{z}{sH(m_{N})}\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}-1\right)(\gamma_{N\to HL}+2\gamma_{NL\to qt}+4\gamma_{Nt\to qL})$ (23)
$\displaystyle-$
$\displaystyle\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}\right)^{2}-1\right)(2\gamma_{NN\to\rho\rho,Z^{\prime}Z^{\prime}}+2\gamma_{NN\to\rho
Z^{\prime},hZ^{\prime}}+2\gamma_{NN\to{\rm SM}{\rm SM}})$ $\displaystyle+$
$\displaystyle\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{\chi}}{Y_{\chi}^{\rm
eq}}\right)^{2}-\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}\right)^{2}\right)2\gamma_{\chi\chi\to NN},$
$\displaystyle\frac{dY_{\chi}}{dz}$ $\displaystyle=$
$\displaystyle-\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{\chi}}{Y_{\chi}^{\rm
eq}}\right)^{2}-1\right)(2\gamma_{\chi\chi\to\rho\rho,Z^{\prime}Z^{\prime}}+2\gamma_{\chi\chi\to\rho
Z^{\prime},hZ^{\prime}}+2\gamma_{\chi\chi\to{\rm SM}{\rm SM}})$ (24)
$\displaystyle-$
$\displaystyle\frac{z}{sH(m_{N})}\left(\left(\frac{Y_{\chi}}{Y_{\chi}^{\rm
eq}}\right)^{2}-\left(\frac{Y_{N}}{Y_{N}^{\rm
eq}}\right)^{2}\right)2\gamma_{\chi\chi\to NN},$
$\displaystyle\frac{dY_{B-L}}{dz}$ $\displaystyle=$
$\displaystyle\frac{z}{sH(m_{N})}\left(\varepsilon_{{\rm
CP}}\left(\frac{Y_{N}}{Y_{N}^{\rm eq}}-1\right)-\frac{Y_{B-L}}{2Y_{L}^{{\rm
eq}}}\right)\gamma_{N\to HL}$ (25) $\displaystyle-$
$\displaystyle\frac{z}{sH(m_{N})}\frac{Y_{B-L}}{Y_{L}^{{\rm
eq}}}\left(\frac{Y_{N}}{Y_{N}^{{\rm eq}}}\gamma_{NL\to qt}+2\gamma_{Nt\to
qL}\right),$

where the notations are the same as those in Eq.12. It should be emphasized
that for TeV scale $N$ and $\chi$ in Fig.1 (d), the SM fermion final states
are dominantly mediated by $Z^{\prime}$, while the SM boson final states are
dominantly mediated by $\rho$. For convenience to distinguish different
contributions, we therefore use $NN(\chi\chi)\to VV(V=W,Z,h)$ and
$NN(\chi\chi)\to f\bar{f}$ to represent the $\rho$ and $Z^{\prime}$ mediated
channels respectively. The new reduced cross sections arising here, i.e.,
$NN\to\rho Z^{\prime},hZ^{\prime},Z^{\prime}Z^{\prime},{\rm SM}{\rm SM}$ and
$\chi\chi\to\rho Z^{\prime},hZ^{\prime},Z^{\prime}Z^{\prime},{\rm SM}{\rm
SM},NN$ have specific expressions shown in Appendix A.

Different from the global scenario, there are various tight constraints on the
new gauge boson $Z^{\prime}$ in the local scenario. Searches of $Z^{\prime}$
in the dilepton final states have excluded $m_{Z^{\prime}}\lesssim 5.15$ TeV
with 140 fb-1 data at LHC ATLAS:2019erb ; CMS:2021ctt . Meanwhile, electroweak
precision test limits from LEP require
$m_{Z^{\prime}}/g_{B-L}=2v_{\phi}\gtrsim 7$ TeV Cacciapaglia:2006pk . In the
following study, we fix $m_{Z^{\prime}}=7$ TeV for illustration.

Figure 9: The evolutions of abundances $Y_{i}$ with $m_{N}=1000~{}{\rm GeV}$,
$m_{\rho}=500~{}{\rm GeV}$, $m_{Z^{\prime}}=7000~{}{\rm GeV}$, and
$\varepsilon_{{\rm CP}}=0.1$ in the local scenario. Subfigures (a)-(d)
correspond to the four cases with mass ratios $m_{N}/m_{\chi}$= 0.5,1,1.5 and
2. The markers of the curves are consistent with that in Fig.2 except for the
blue and green lines, in which solid, dashed and dot-dashed correspond to
$v_{\phi}=4000~{}{\rm GeV}$, $5000~{}{\rm GeV}$, and $6000~{}{\rm GeV}$,
respectively.

We first take some benchmarks to demonstrate the evolution of abundances of
various particles in Fig. 9. In panel (a) for leptogenesis, the dominant
$NN\to\rho\eta$ channel in the global scenario is replaced by the new process
$NN\to\rho Z^{\prime}$. However, the contribution of $NN\to\rho Z^{\prime}$ to
relic density is negligible, since $2m_{N}$ is much smaller than
$m_{\rho}+m_{Z}^{\prime}$ for the benchmark scenarios. We find that the
process $NN\to f\bar{f}$ plays a major contribution, while $NN\to\rho\rho$ is
subdominant. It is clear that decreasing $v_{\phi}$ results in a smaller
$Y_{B-L}$ for a fixed value of $\varepsilon_{{\rm CP}}$. To meet the
observation value of baryon asymmetry, $v_{\phi}$ typically needs to be
smaller than $1000$ GeV for $\varepsilon_{{\rm CP}}=0.1$, which is already
excluded by LEP-II Cacciapaglia:2006pk ; ALEPH:2013dgf  in the local scenario.
So the only viable method is reducing $\varepsilon_{{\rm CP}}$. For DM $\chi$,
the dominant annihilation channel is $\chi\chi\to f\bar{f}$, as
$\chi\chi\to\rho Z^{\prime},Z^{\prime}Z^{\prime}$ are forbidden kinematically
for the benchmark points.

The results for increasing $m_{N}/m_{\chi}$ are shown in panels (b)-(d) of
Fig. 9. Although the reduced cross section of $NN\to\chi\chi$ should become
larger, $Y_{B-L}$ has hardly changed since this channel is not the dominant
contribution. Meanwhile, the decrease of $m_{\chi}$ leads to the increase of
$Y_{\chi}$. Under the LEP-II constraint $v_{\phi}\gtrsim 3500$ GeV, the DM
abundances are much larger than the Planck observation in the non-resonance
region of the benchmark points. Therefore, we can not find the common
parameter space of leptogenesis and DM in the non-resonance region for the
special scenario illustrated. We need a more efficient pathway to produce DM.

Figure 10: $Y_{i}$ as a function of $m_{\rho}$ at present in the local
scenario. $m_{N}$ and $m_{Z^{\prime}}$ are fixed as 1000 and 7000 GeV. The
markers of the curves can be seen in Fig. 3 except for blue and green lines,
in which the solid, dashed, and dot-dashed represent $v_{\phi}=4000~{}{\rm
GeV}$, $5000~{}{\rm GeV}$, and $6000~{}{\rm GeV}$, respectively.

As shown in Fig. 3, the abundances of DM can be greatly reduced in the
resonance region for the global scenario, which portends that leptogenesis and
DM may also coexist in this region for the local scenario. Then, Fig. 10 shows
the influences of $m_{\rho}$ on abundances $Y_{i}$, where the resonance effect
is obvious. In the resonance region at $m_{\rho}\simeq 2m_{N}$, the dominant
channel is $NN\to\eta\eta/VV$ for the global scenario, while it is substituted
for $NN\to VV$ in the local scenario. It should be noted that the $\rho$
portal channels depend on the mixing angle $\theta$, where we use
$\theta=0.05$ in the calculations. In the decoupling limit with $\theta=0$,
the resonance effect will be absent, because $\rho$ can not mediate $NN\to{\rm
SM}{\rm SM}$ process in this special case for the benchmark points with
$m_{N}/m_{\chi}=0.5$. Outside of the resonance region, the dominant channel is
$NN\to\rho\rho$ for light $m_{\rho}$, and then becomes $NN\to f\bar{f}$ when
$m_{\rho}>m_{N}$.

Besides the resonance happens at $m_{\rho}\simeq 2m_{\chi}$, the results of DM
are similar to the sterile neutrino. In panel (a) of Fig. 10 with DM heavier
than sterile neutrino, the $\chi\chi\to NN$ channel is always allowed, which
also contributes to DM annihilation at the $\rho$ resonance. As for panel
(b)-(d), $Y_{B-L}$ shows an obvious resonance effect at $m_{\rho}\simeq
2m_{N}$ under the influence of $NN\to VV$ and $NN\to\chi\chi$. Therefore, the
symmetric $m_{\rho}$ in this region can successfully explain baryon asymmetry
without modifying $\varepsilon_{{\rm CP}}$. In these three cases, DM is not
heavier than sterile neutrino, which indicates that the impact of $\chi\chi\to
NN$ can be ignored even in the resonance region. In addition, the overall
increase of $Y_{\chi}$ can not change the fact that the correct DM abundance
can only be obtained at position $m_{\rho}\simeq 2m_{\chi}$. These results
confirm the condition that to simultaneously satisfy $Y_{B-L}$ and $Y_{\chi}$,
$m_{\rho}$ could appear in the DM resonance region as $m_{\rho}\simeq
2m_{\chi}$, while changing $\varepsilon_{{\rm CP}}$.

As already shown in Fig. 1, DM and sterile neutrino $N$ also annihilate via
the new gauge boson $Z^{\prime}$. Therefore, correct DM abundance can be
obtained at the $Z^{\prime}$ resonance region as $m_{Z^{\prime}}\simeq
2m_{\chi}$ Blanco:2019hah . Varying $m_{Z^{\prime}}$ would lead to quite
similar results as varying $m_{\rho}$ shown in Fig. 10. In this paper, we have
fixed $m_{Z^{\prime}}=7000$ GeV, which means that the $Z^{\prime}$ resonance
appears at DM mass $m_{\chi}\simeq m_{Z^{\prime}}/2=3500$ GeV. This result is
clearly shown in Fig. 11.

###  IV.2 Perturbation and CP-asymmetry

Similar to the global scenario, we also perform a scan in searching for the
common origin of leptogenesis and DM in the parameter ranges as

$\displaystyle\begin{aligned} m_{\chi}\in[400,10000]~{}{\rm
GeV},m_{\rho}\in[100,10000]~{}{\rm GeV},v_{\phi}\in[3500,20000]~{}{\rm
GeV},\varepsilon_{{\rm CP}}\in[10^{-6},1].\end{aligned}$ (26)

During the scan, we have fixed $m_{Z^{\prime}}=7000$ GeV and $\theta=0.05$ for
illustration. The range of $v_{\phi}$ in the local scenario is different from
the global scenario because $v_{\phi}\lesssim 3500$ GeV is already excluded by
LEP-II Cacciapaglia:2006pk ; ALEPH:2013dgf . The gauge coupling is determined
by $g_{B-L}=m_{Z^{\prime}}/(2v_{\phi})$, which is smaller than one in the
above range of $v_{\phi}$.

In the local scenario, the perturbation constraints are
$\lambda_{N,\chi}<\sqrt{4\pi}$ and $\lambda_{\phi}<2\pi$ Kahlhoefer:2015bea ;
Duerr:2016tmh , meanwhile the CP-asymmetry is required to satisfy
$\varepsilon_{\rm CP}<1$. The scanning results for the four specific cases
with mass ratios $m_{N}/m_{\chi}=0.5,1,1.5,2$ are displayed in Fig. 11. Panel
(a) is the case with light sterile neutrino, where all samples satisfy the
perturbation constraints. A few samples with $\varepsilon_{{\rm CP}}>1$ occur
at the $\rho$ or $Z^{\prime}$ resonance due to over dilution. The $\rho$
portal resonance is attributed to $NN\to VV$ with $m_{\rho}\simeq
2m_{N}=m_{\chi}\sim 3500$ GeV, while the $Z^{\prime}$ portal resonance is
dominant by $NN\to f\bar{f}$ with $m_{Z^{\prime}}\simeq 2m_{N}=m_{\chi}\sim
7000$ GeV. For the DM, there are three viable regions for the common origin in
case (a), i.e., the $\rho$ resonance, the $Z^{\prime}$ resonance, and the
heavy DM $m_{\chi}\gtrsim 6000$ GeV. Under the strict constraints of the
present direct detection bound of LZ LZ:2022lsv , $m_{\rho}\lesssim 200$ GeV
is required in the heavy DM region, where the $s$-wave $\chi\chi\to\rho
Z^{\prime}$ is the dominant annihilation channel. Finally, these allowed
samples are accompanied by relatively tiny $\varepsilon_{{\rm CP}}$ in the
resonance region, approximately $\in[10^{-5},10^{-1}]$. But in the non-
resonance region, $\varepsilon_{{\rm CP}}$ has much larger values, namely
$\varepsilon_{{\rm CP}}\in[10^{-1},1]$.

Figure 11: Perturbation and exceeding $\varepsilon_{{\rm CP}}$ constraints in
the local scenario. (a)-(d) correspond to four cases with mass ratios
$m_{N}/m_{\chi}=0.5,1,1.5,2$. The markers of the samples are the same as those
in Fig. 4.

Increasing the mass ratio $m_{N}/m_{\chi}$ leads to different results, which
are shown in panels (b)-(d) of Fig. 11. The $\rho$ and $Z^{\prime}$ resonance
regions for DM are still viable, where allowed samples have a larger
$\varepsilon_{{\rm CP}}$ compared to panel (a). Meanwhile, the area excluded
by perturbation constraints expands with $\lambda_{N}\geqslant\sqrt{4\pi}$ as
$m_{N}$ increases. Due to excessive dilution of $N$, the positions of $\rho$
and $Z^{\prime}$ resonance related to $N$ are disfavored by $\varepsilon_{{\rm
CP}}<1$. However, the $\varepsilon_{{\rm CP}}>1$ excluded samples at resonance
vary for different mass ratio $m_{N}/m_{\chi}$. For case (b) with
$m_{N}=m_{\chi}$, the resonance regions for DM are the same as $N$, so a lot
of samples in these resonance regions are also excluded by $\varepsilon_{{\rm
CP}}>1$. The $Z^{\prime}$ resonance related to $N$ excludes
$m_{\chi}=2m_{N}/3\simeq m_{Z^{\prime}}/3$ for case (c) and
$m_{\chi}=m_{N}/2\simeq m_{Z^{\prime}}/4$ for case (d) where
$m_{Z^{\prime}}=7000$ GeV. The $\rho$ resonance related to $N$ in principle
excludes $m_{\rho}\simeq 2m_{N}=3m_{\chi}\simeq 3m_{Z^{\prime}}/2$ for case
(c) and $m_{\rho}\simeq 2m_{N}=4m_{\chi}\simeq 2m_{Z^{\prime}}$ for case (d),
however the latter case is out the scanned range of $m_{\rho}$. In panels (c)
and (d) of Fig. 11, we also find that the samples in the non-resonance region
are fully excluded by the limitations of perturbation and $\varepsilon_{{\rm
CP}}$. The sterile neutrino heavier than the DM is totally not allowed in the
non-resonance region for the local scenario, which is different from the
global scenario.

###  IV.3 Collider Signal and Direct detection

Figure 12: The current LEP-II Cacciapaglia:2006pk ; ALEPH:2013dgf  and
prospective ILC limits LCCPhysicsWorkingGroup:2019fvj  in the local scenario.
(a)-(d) correspond to four cases with mass ratios
$m_{N}/m_{\chi}=0.5,1,1.5,2$. The black dashed line represents the prospective
ILC limit with $\sqrt{s}=250$ GeV. The markers of the samples are the same as
that in Fig. 4.

The collider signature of gauged $U(1)_{B-L}$ is extensively studied
Basso:2008iv ; Iso:2009nw . In the local scenario, the DM $\chi$ could
contribute to Higgs invisible decay via $h\to\chi\chi$ when
$m_{\chi}<m_{h}/2$, but it is forbidden kinematically with the parameter space
in Eqn. (26). The golden channel for discovering the gauge boson $Z^{\prime}$
is the dilepton channel $pp\to Z^{\prime}\to\ell^{+}\ell^{-}$, which has
excluded $m_{Z^{\prime}}\lesssim 5.15$ TeV at LHC ATLAS:2019erb ; CMS:2021ctt
. The gauge boson $Z^{\prime}$ can also mediate the pair production of sterile
neutrino as $pp\to Z^{\prime}\to NN$ FileviezPerez:2009hdc , which has
excluded the region with $m_{N}\lesssim 1.4$ TeV and $m_{Z^{\prime}}\lesssim
4.4$ TeV CMS:2023ooo . The choice of $m_{Z^{\prime}}=7000$ GeV with TeV scale
$m_{N}$ satisfies these current experimental limits. In the future, the 250
GeV ILC could exclude $v_{\phi}\lesssim 22$ TeV via the process $e^{+}e^{-}\to
f\bar{f}$ LCCPhysicsWorkingGroup:2019fvj ; KA:2023dyz .

The scanned results in the $m_{\chi}-v_{\phi}$ plane are shown in Fig. 12. In
panel (a), it is obvious that the allowed samples at the $\rho$ resonance
satisfy $v_{\phi}\lesssim 11$ TeV. By comparison, $v_{\phi}$ in the
$Z^{\prime}$ resonance can reach a larger value, but it will not exceed $20$
TeV. Meanwhile, the larger the $v_{\phi}$ is, the smaller the CP-asymmetry
$\varepsilon_{\rm CP}$ is required, especially clear in the $Z^{\prime}$
resonance region. For 10 TeV scale heavy DM in the non-resonance region,
survived samples require relatively larger coupling thus smaller $v_{\phi}$ to
satisfy observed relic abundance, which is less than $5500$ GeV here.
Certainly, samples with smaller $v_{\phi}$ tend to have larger scattering
cross section, hence are also more tightly constrained by LZ. Based on the
above results, we can conclude that the future 250 GeV ILC is enough to
completely exclude all the samples within the parameter space in Eqn. (26),
even at the resonance region.

In panel (b) of Fig. 12, we find that the excluded samples at the $Z^{\prime}$
resonance by the constraint of exceeding $\varepsilon_{{\rm CP}}>1$
concentrate in the sharp peak, which leads to the allowed samples satisfying
$v_{\phi}\lesssim 8000$ GeV. Meanwhile, the few survived samples via the
$\rho$ resonance indicate $v_{\phi}\lesssim 6000$ GeV. In panels (c)-(d), as
mentioned in Fig. 11, only some samples at regions of $\rho$ and $Z^{\prime}$
resonance are left under various current limits. But they can not avoid the
constraints of future 250 GeV ILC either. Except for the $Z^{\prime}$
resonance related to $N$, it is clear that a lighter DM $\chi$, i.e., a
lighter sterile neutrino $N$, requires smaller $\varepsilon_{\rm CP}$ for
samples via the $\rho$ resonance.

Figure 13: Spin-independent scattering cross sections obtained through the
Higgs portal in Eq. (22) (blue line), through tree-level $Z^{\prime}$ in Eq.
(27) (purple line) and loop-level $Z^{\prime}$ in Eq. (28) (lime line). The
parameters $\theta$, $v_{\phi}$, and $m_{\rho}$ are fixed as 0.05, 5000 GeV,
and 1000 GeV, respectively. The black solid, dashed and dot-dashed lines
represent the present LZ limit LZ:2022lsv , future LZ:2015kxe  LZ limit, and
the neutrino floor Billard:2013qya .

We then consider the direct detection limits. The elastic scattering cross
sections of $\chi$ can be divided into two kinds according to the mediators.
The Higgs portal $\rho$ and $h$ mediated one are the same as the scenario
discussed in the global scenario, which heavily depends on mixing angle
$\theta$. The expression of the spin-independent scattering cross section can
be seen in Eq. (22). The other one is mediated by $Z^{\prime}$. The spin-
independent scattering cross section at the tree-level can be written
asHooper:2014fda

$\displaystyle\sigma_{\rm
SI}^{\text{tree}}=\frac{2m_{\chi}^{2}m_{n}^{4}(g_{\chi}^{A})^{2}v^{2}}{\pi
m_{Z^{\prime}}^{4}(m_{\chi}+m_{n})^{4}A^{2}}\left[Z(2g_{u}^{V}+g_{d}^{V})+(A-Z)(g_{u}^{V}+2g_{d}^{V})\right]^{2},$
(27)

where $v$ is the non-relativistic velocity of DM at present. $Z$ and $A$ are
the atomic number and atomic mass of the target, respectively. Since $v\sim
10^{-3}$, the scattering cross section through tree-level $Z^{\prime}$ is
heavily suppressed. The velocity unsuppressed scattering via $Z^{\prime}$
appears at the one-loop level. The corresponding spin-independent scattering
cross section is calculated as Kahlhoefer:2015bea

$\displaystyle\sigma_{\rm
SI}^{\text{loop}}=\frac{C^{2}m_{n}^{4}m_{\chi}^{4}(g_{\chi}^{A})^{4}(g_{q}^{V})^{4}}{\pi^{5}(m_{\chi}+m_{n})^{2}m_{\rho}^{4}m_{Z^{\prime}}^{4}},$
(28)

where $q$ in $g_{q}^{V}$ represents quarks.

Figure 14: Constraints from the direct detection experiments in the local
scenario. $\theta$ is also fixed as $0.05$.

The numerical results for a specific scenario are shown in Fig.13. The blue
curve for the Higgs portal scattering cross section is calculated with mixing
$\theta=0.05$, which can be tested by future LZ experiment. For
$m_{\chi}\lesssim 1000$ GeV, the tree-level $Z^{\prime}$ mediated result is
larger than the loop-level one. And the loop induced scattering via
$Z^{\prime}$ becomes larger when DM $\chi$ is heavier than about 1000 GeV. The
scattering cross sections mediated by $Z^{\prime}$ at both tree and loop-level
are smaller than $10^{-54}~{}\rm cm^{2}$, which are over eight orders of
magnitudes smaller than the Higgs portal result. Therefore, we only consider
the spin-independent scattering cross section through the Higgs portal for
simplicity during the scan.

The scanned results of spin-independent scattering cross sections are shown in
Fig.14, where $\theta$ is fixed as 0.05. Except for different excluded samples
under constraints from perturbation and $\varepsilon_{\rm CP}>1$, the
theoretical values of $\sigma_{\rm SI}$ are similar for the four scenarios for
different mass ratios of $m_{N}/m_{\chi}$. The current LZ limit tends to
exclude DM heavier than 3000 GeV, while light DM is safe under the current
limit. This is because $\sigma_{\rm SI}$ is proportional to $m_{\chi}^{4}$
according to Eqn. (22). For certain $v_{\phi}$ determined by DM relic
abundance, the heavier DM thus leads to a larger scattering cross section. For
DM via the $\rho$ portal, allowed samples predict $\sigma_{\rm SI}\gtrsim
3\times 10^{-48}~{}\rm cm^{2}$, which is over the neutrino background. A few
$\rho$ portal samples with $m_{\chi}\lesssim 1500$ GeV would escape the future
LZ limit, but they are within the reach of future DARWIN experiment
DARWIN:2016hyl . However, in panel (b) with $m_{N}/m_{\chi}=1$, future LZ is
likely to exclude all allowed samples at the $\rho$ resonance because of
exceeding $\varepsilon_{{\rm CP}}$. On the other hand, the $Z^{\prime}$
resonance and non-resonance samples could have relatively small $\sigma_{\rm
SI}$, mainly due to the cancellation effect with $m_{\rho}\simeq m_{h}$. Once
again, the above direction detection results heavily depend on the mixing
angle $\theta=0.05$. A smaller mixing as $\theta=0.01$ would result in most of
the samples under the LZ limit.

###  IV.4 Indirect Detection

The theoretical DM annihilation channels in the local scenario are also shown
in Fig. 1. Different from the global scenario, the dominant DM annihilation
channels at present heavily depend on the mass spectrum of the new particles.
Due to the Majorana nature of DM $\chi$, the DM annihilation channels into SM
final states via $\rho$ and $Z^{\prime}$ are $p$-wave suppressed. So most of
the samples in the $\rho$ or $Z^{\prime}$ resonance region are beyond the
reach of indirect detection. For samples in the non-resonance region, the
dominant annihilation channel is the $s$-wave $\chi\chi\to\rho Z^{\prime}$,
which is hopeful to be detected.

Figure 15: Constraints from the indirect detection experiments in the local
scenario. The newly added $\tau^{+}\tau^{-}$ dashed line represents the
projected upper limits of $\tau^{+}\tau^{-}$ final states from the Cherenkov
Telescope Array CTA:2020qlo ; Mangipudi:2021ivm . Other labels are the same as
that in Fig. 8. The allowed samples are marked with heterogeneous shapes to
represent different dominant final states and colored by $\varepsilon_{{\rm
CP}}$.

The results of $\left<\sigma v\right>$ at present with a variety of bounds are
shown in Fig. 15. In panel (a), we first illustrate the case with light
sterile neutrino. At the time of DM freeze out, both $\chi\chi\to\text{SM SM}$
and $\chi\chi\to NN$ via $\rho$ or $Z^{\prime}$ contribute to DM annihilation.
As the former channel is $p$-wave suppressed at present, the dominant
annihilation channel becomes $\chi\chi\to NN$. With a subdominant contribution
to relic density, the $\chi\chi\to NN$ channel has an annihilation cross
section $\left<\sigma v\right>$ less than $10^{-27}~{}\rm{cm^{3}/s}$, which is
far below the constraints of Fermi-LAT and H.E.S.S. on the $NN$ final states
Campos:2017odj . The small annihilation cross section $\left<\sigma v\right>$
of $\chi\chi\to NN$ also makes it beyond the scope of future CTA sensitivity.
Close to the sharp $Z^{\prime}$ resonance, there are a few samples that
dominantly annihilate into $f\bar{f}$ at present. However, their annihilation
cross section is at the order of $\mathcal{O}(10^{-31})~{}\rm{cm^{3}/s}$,
which is also difficult to be detected.

In panels (b), (c), and (d) of Fig. 15, the $s$-wave $\chi\chi\to NN$ channel
is kinematically disallowed at present due to $m_{N}\geq m_{\chi}$. Therefore,
the permitted samples at the $\rho$ resonance dominantly annihilate into
$W^{+}W^{-}$, while samples at the $Z^{\prime}$ resonance annihilate into
$f\bar{f}$. Both of these two channels are $p$-wave, so the corresponding
current annihilation cross sections are less than
$\mathcal{O}(10^{-29})~{}\rm{cm^{3}/s}$, which is far below the experimental
sensitivity.

The promising annihilation channel is the $s$-wave process $\chi\chi\to\rho
Z^{\prime}$ in the non-resonance region, which has a typical annihilation
cross section of $2\times 10^{-26}~{}\rm{cm^{3}/s}$. Although part of such
samples may suffer tight constraints from direct detection, they can always
escape this limit by lowing $\theta$ while keeping $\left<\sigma v\right>$ of
$\chi\chi\to\rho Z^{\prime}$ unchanged. Different from the global scenario,
the additional scalar singlet $\rho$ does not decay into the invisible
Majoron. From Fig. 11, we know that $\rho\lesssim 200$ GeV when $\theta=0.05$.
In this way, $\rho\to W^{+}W^{-}$ becomes the dominant decay mode, since
$\rho\to Z^{\prime}Z^{\prime}$ and $\rho\to NN$ are kinematically disallowed.
Meanwhile, the new gauge boson mainly decays into fermion pair as
$Z^{\prime}\to f\bar{f}$. The detailed calculation of the gamma-ray spectrum
in the full dark matter annihilation channel of $\chi\chi\to\rho Z^{\prime}$
is beyond the scope of this paper. For illustration, we consider the future
CTA limit of the individual $W^{+}W^{-}$ and $\tau^{+}\tau^{-}$ final states.
It is clear that the samples with $m_{\chi}\gtrsim 5800$ GeV are promising in
the future CTA experiment for cases (a) and (b) of Fig. 15. Although the
$\chi\chi\to\rho Z^{\prime}$ channel could also lead to an observable
signature at CTA for case (c) and (d) with sterile neutrino heavier than DM,
such samples are already excluded by perturbation and exceeding
$\varepsilon_{{\rm CP}}$.

Furthermore, close to the $Z^{\prime}$ resonance with $2m_{\chi}\gtrsim
m_{Z^{\prime}}+m_{\rho}$ for light $m_{\rho}$, although some samples are
generated through $\chi\chi\to f\bar{f}$ via $Z^{\prime}$ at freeze out, the
current annihilation of them are dominated by the $s$-wave process
$\chi\chi\to\rho Z^{\prime}$. Due to the subdominant contribution to the DM
abundance, the annihilation cross section of these samples is typically less
than $10^{-26}~{}\rm{cm^{3}/s}$, which is below the CTA $W^{+}W^{-}$
sensitivity. Meanwhile, a few samples are above the $\tau^{+}\tau^{-}$
sensitivity, so these samples might still be promising when considering the
full decays of $\chi\chi\to\rho Z^{\prime}$. On the other hand, when
$2m_{\chi}<m_{Z^{\prime}}+m_{\rho}$ at the $Z^{\prime}$ resonance, the current
dominant channels would be the $s$-wave process $\chi\chi\to hZ^{\prime}$.
However, the corresponding annihilation cross section of certain samples are
suppressed by the mixing $\theta=0.05$, which results in $\left<\sigma
v\right>\lesssim 10^{-28}~{}\rm{cm^{3}/s}$. So such samples are clearly beyond
the experimental reach.

##  V Conclusion

This $U(1)_{B-L}$ model further extends the SM with sterile neutrinos $N$,
Majorana dark matter $\chi$, and complex scalar $\phi$, which are charged
under the $U(1)_{B-L}$ symmetry. An additional $Z_{2}$ symmetry is also
employed to keep the stability of DM $\chi$. After the spontaneous breaking of
$U(1)_{B-L}$ symmetry, both $N$ and $\chi$ obtain masses through the Yukawa
interactions related to $\phi$. Besides generating the tiny neutrino masses
via type-I seesaw, the out-of-equilibrium decays of TeV scale $N$ create the
baryon asymmetry through the resonance leptogenesis. Furthermore, the dilution
of $N$ has a significant impact on leptogenesis with new particles in the
$U(1)_{B-L}$. Since the evolution of $N$ and $\chi$ involve quite similar
processes in the calculation of baryon asymmetry and DM relic density, the
focus of this paper is to explore the common parameter space of these two
aspects under various constraints. Both the global $U(1)_{B-L}$ and local
$U(1)_{B-L}$ scenario are studied in this work.

In the global $U(1)_{B-L}$ scenario, we have a scalar singlet $\rho$ and a
(nearly) massless Majoron $\eta$ after the symmetry breaking. The evolution of
$N$ and $\chi$ are dominantly determined by the pair annihilation processes
$NN,\chi\chi\to\rho\rho,\rho\eta,\eta\eta$. Meanwhile, the conversion process
$\chi\chi\leftrightarrow NN$ mediated by $\rho$ and $\eta$ could also impact.
We have chosen four scenarios, namely $m_{N}/m_{\chi}=0.5,1,1.5,2$, to
illustrate the phenomenology. In addition, the mixing angle $\theta$ between
$\rho$ and Higgs $h$ is important for DM direct detection, which is assigned
to an appropriate value of $0.05$. Under various constraints, correct DM
abundance and baryon asymmetry can be realized at the resonance region of
$m_{\rho}\simeq 2m_{\chi}$ or the non-resonance region of $m_{\rho}\lesssim
200$ GeV. In the resonance region, the required CP-asymmetry
$\varepsilon_{{\rm CP}}$ can be as small as $\mathcal{O}(10^{-5})$, while the
non-resonance samples usually need $\varepsilon_{{\rm CP}}\gtrsim 10^{-2}$.
For the non-resonance region, the perturbation and exceeding
$\varepsilon_{{\rm CP}}$ constraints become tighter when $m_{N}/m_{\chi}$
increases. Although the Higgs invisible decay induced by $h\to\eta\eta$ can
not be detected by LHC, the future CEPC may have a positive signature. With
$\theta=0.05$, most of the samples could lead to detectable DM-nucleon
scattering cross section. For DM annihilation at present, the $s$-wave process
$\chi\chi\to\rho\eta$ mainly decays into neutrinos, which is far below
experimental limits. Due to additional contribution from $p$-wave
$\chi\chi\to\eta\eta$ at freeze out, the $\chi\chi\to NN$ channel typically
has an annihilation cross section smaller than $2\times 10^{-26}~{}\rm
cm^{3}/s$, which leads to most samples beyond the reach of CTA. So the global
$U(1)_{B-L}$ scenario is usually unpromising in indirect detection
experiments.

In the local $U(1)_{B-L}$ scenario, the Majoron $\eta$ is replaced by the
gauge boson $Z^{\prime}$. We have fixed $m_{Z^{\prime}}=7000$ GeV in this
study to respect current experimental limits. The processes as
$NN,\chi\chi\to{\rm SM}{\rm SM}$ via $\rho$ and $Z^{\prime}$ become important.
Four scenarios with $m_{N}/m_{\chi}=0.5,1,1.5,2$ are also chosen for
illustration. Under various experimental limits, the common origin of DM and
leptogenesis is viable at the $\rho$ resonance $m_{\rho}\simeq 2m_{\chi}$ with
$\varepsilon_{{\rm CP}}\gtrsim\mathcal{O}(10^{-5})$ or the $Z^{\prime}$
resonance $m_{Z^{\prime}}\simeq 2m_{\chi}$ with $\varepsilon_{{\rm
CP}}\gtrsim\mathcal{O}(10^{-3})$. For the non-resonance region, we find that
it is only possible for $m_{\chi}\gtrsim 6000$ GeV and $m_{N}/m_{\chi}\leq 1$
with $\varepsilon_{{\rm CP}}\gtrsim\mathcal{O}(10^{-1})$. When the sterile
neutrino $N$ is heavier than the DM, the non-resonance samples are excluded by
exceeding $\varepsilon_{{\rm CP}}$. The existence of the new gauge boson
$Z^{\prime}$ makes this scenario more promising at colliders. In the future,
the 250 GeV ILC could probe all the allowed parameter space by precise
measurement of the $e^{+}e^{-}\to f\bar{f}$ process. The DM-nucleon scattering
is dominant by the Higgs portal processes, so we fix $\theta=0.05$ for
illustration as in the global scenario. In this way, the future direct
detection experiment would test a large portion of the allowed parameter
space. At both resonance regions, the $\chi\chi\to{\rm SM}{\rm SM}$ processes
are $p$-wave suppressed at present due to the Majorana nature of $\chi$, thus
are hard to detect. The promising case for the future CTA experiment is the
non-resonance region with $\chi\chi\to\rho Z^{\prime}$, where $\rho\to
W^{+}W^{-}$ and $Z^{\prime}\to f\bar{f}$ are the dominant decay modes.

To summarize, both the global and local $U(1)_{B-L}$ scenario favor the
resonance region for common origin of DM and leptogenesis. However, the
resonance region is hard to test at indirect detection experiments for both
scenarios. On the other hand, correct relic abundance and baryon asymmetry can
also be simultaneously satisfied in the non-resonance region. For the global
$U(1)_{B-L}$, electroweak scale DM is viable, but is still hard to be detected
at indirect detection. For the local $U(1)_{B-L}$, DM must be over TeV scale,
but is within the future CTA limit. The Majoron $\eta$ could induce testable
Higgs invisible decay in the global scenario, while the gauge boson
$Z^{\prime}$ could be detected at colliders in the local scenario. These
differences can be used to distinguish the global and local scenarios.

## Acknowledgments

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

##  Appendix A The reduced cross sections

The reduced cross sections of annihilation processes involved in the Boltzmann
equations are shown in this appendix. We use CalcHEP Belyaev:2012qa
throughout to calculate them. In the global $U(1)_{B-L}$ scenario, our results
match with the Ref. Dev:2017xry . The reduced cross sections of sterile
neutrino pair annihilation processes $NN\to\rho\rho,\rho\eta,\eta\eta$ are

$\displaystyle\hat{\sigma}(NN\to\rho\rho)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{N}^{4}}{32\pi}\bigg{(}\frac{9\beta_{\rho\rho}(x_{N}-4)r_{\rho
N}^{2}}{(x_{N}-r_{\rho N})^{2}}-\frac{24\beta_{\rho\rho}x_{N}r_{\rho
N}}{x_{N}|x_{N}-r_{\rho N}|}+\frac{12r_{\rho N}(x_{N}+2r_{\rho
N}-8)}{x_{N}|x_{N}-r_{\rho N}|}$ (29) $\displaystyle\times$
$\displaystyle\log\left(\frac{(1-\beta_{\rho\rho})x_{N}-2r_{\rho
N}}{(1+\beta_{\rho\rho})x_{N}-2r_{\rho
N}}\right)-4\beta_{\rho\rho}-\frac{8\beta_{\rho\rho}(r_{\rho
N}-4)^{2}}{(x_{N}-2r_{\rho
N})^{2}-\beta_{\rho\rho}^{2}x_{N}^{2}}-\frac{2}{x_{N}(x_{N}-2r_{\rho N})}$
$\displaystyle\times$ $\displaystyle\Big{(}x_{N}^{2}-2(r_{\rho
N}-4)(2x_{N}-3r_{\rho
N}-4)\Big{)}\log\left(\frac{(1-\beta_{\rho\rho})x_{N}-2r_{\rho
N}}{(1+\beta_{\rho\rho})x_{N}-2r_{\rho N}}\right)\bigg{)},$
$\displaystyle\hat{\sigma}(NN\to\rho\eta)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{N}^{4}}{16\pi
x_{N}^{2}}\left(\beta_{\rho\eta}r_{\rho N}^{2}x_{N}-2(x_{N}-r_{\rho
N})(x_{N}-2r_{\rho N})\log\left(\frac{(1-\beta_{\rho\eta})x_{N}-r_{\rho
N}}{(1+\beta_{\rho\eta})x_{N}-r_{\rho N}}\right)\right),$ (30)
$\displaystyle\hat{\sigma}(NN\to\eta\eta)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{N}^{4}}{32\pi(x_{N}-r_{\rho
N})^{2}}\left(\beta_{\eta\eta}x_{N}(r_{\rho N}^{2}-4x_{N})-2(x_{N}^{2}-r_{\rho
N}^{2})\log\left(\frac{1-\beta_{\eta\eta}}{1+\beta_{\eta\eta}}\right)\right),$
(31)

where $x_{i}=s^{\prime}/m_{i}^{2}$, $r_{ij}=m_{i}^{2}/m_{j}^{2}$. In the above
equations, we have considered a vanishing mixing $\theta$ as Ref. Dev:2017xry
to obtain the simplified expressions. The Yukawa coupling $\lambda_{N}$ is
evaluated as $\lambda_{N}=m_{N}/v_{\phi}$ in the calculation. For $N$ pair
annihilation, the $\beta_{ij}$ function is defined as

$\displaystyle\beta_{ij}=x_{N}^{-1}\sqrt{(1-4x_{N}^{-1})(x_{N}^{2}+r_{iN}^{2}+r_{jN}^{2}-2x_{N}r_{iN}-2x_{N}r_{jN}-2r_{iN}r_{jN})}$
(32)

On the other hand, the annihilation channels into SM final states depend on
the mixing $\theta$. The $NN\to VV$ ($V=W,Z,h$) process is dominant by the
$\rho$ portal for TeV scale $N$, which is calculated as

$\displaystyle\hat{\sigma}(NN\to VV)$ $\displaystyle=$
$\displaystyle\frac{1}{32\pi v_{\phi}^{2}(x_{N}-r_{\rho
N})^{2}}\bigg{(}\beta_{WW}(x_{N}-4)\cos^{2}\theta\Big{(}\cos
2\theta(\lambda_{H\phi}^{2}v_{\phi}^{2}-4\lambda_{H}^{2}v_{H}^{2})$ (33)
$\displaystyle+$
$\displaystyle\frac{g_{e}^{2}\sin^{2}\theta}{\sin^{4}\theta_{w}}(2g_{e}^{2}v_{H}^{2}-m_{W}^{2}\sin^{2}\theta_{w}-2m_{N}^{2}\sin^{2}\theta_{w}x_{N})+\lambda_{H\phi}^{2}v_{\phi}^{2}$
$\displaystyle+$ $\displaystyle
4\lambda_{H}^{2}v_{H}^{2}+4\lambda_{H\phi}\lambda_{H}v_{\phi}v_{H}\sin
2\theta\Big{)}$ $\displaystyle+$
$\displaystyle\frac{\beta_{ZZ}(x_{N}-4)\cos^{2}\theta}{2}\Big{(}\cos
2\theta(\lambda_{H\phi}^{2}v_{\phi}^{2}-4\lambda_{H}^{2}v_{H}^{2})+\frac{g_{e}^{2}\sin^{2}\theta}{\sin^{4}\theta_{w}\cos^{4}\theta_{w}}$
$\displaystyle\times$
$\displaystyle(2g_{e}^{2}v_{H}^{2}-m_{Z}^{2}\sin^{2}\theta_{w}\cos^{2}\theta_{w}-2m_{N}^{2}\sin^{2}\theta_{w}\cos^{2}\theta_{w}x_{N})+\lambda_{H\phi}^{2}v_{\phi}^{2}$
$\displaystyle+$ $\displaystyle
4\lambda_{H}^{2}v_{H}^{2}+4\lambda_{H\phi}\lambda_{H}v_{\phi}v_{H}\sin
2\theta\Big{)}$ $\displaystyle+$
$\displaystyle\beta_{hh}(x_{N}-4)\cos^{2}\theta\Big{(}\lambda_{H\phi}v_{\phi}\cos^{3}\theta+\lambda_{H\phi}v_{H}\sin^{3}\theta+2(3\lambda_{H}-\lambda_{H\phi})v_{H}$
$\displaystyle\times$
$\displaystyle\cos^{2}\theta\sin\theta+2(3\lambda_{H}-\lambda_{H\phi})v_{\phi}\sin^{2}\theta\cos\theta\Big{)}^{2}\bigg{)},$

where $g_{e}$ and $\theta_{w}$ are the electric coupling constant and Weinberg
angle, respectively. At present, the $s$-wave $NN\to h\eta$ channel is also
important, which is calculated as

$\displaystyle\hat{\sigma}(NN\to h\eta)$ $\displaystyle=$
$\displaystyle\frac{\lambda_{N}^{2}}{16\pi
x_{N}^{2}m_{N}^{2}}\bigg{(}4\lambda_{H\phi}m_{N}v_{H}\lambda_{N}\sin\theta\cos\theta(r_{hN}-x_{N})\log\left(\frac{(1-\beta_{h\eta})x_{N}-r_{hN}}{(1+\beta_{h\eta})x_{N}-r_{hN}}\right)$
(34) $\displaystyle+$ $\displaystyle
2m_{N}^{2}\sin^{2}\theta\bigg{(}2\beta_{h\eta}\lambda_{\phi}^{2}x_{N}+\lambda_{N}^{2}(r_{hN}-x_{N})(4\lambda_{\phi}+x_{N}\lambda_{N}^{2})$
$\displaystyle\times$
$\displaystyle\log\left(\frac{(1-\beta_{h\eta})x_{N}-r_{hN}}{(1+\beta_{h\eta})x_{N}-r_{hN}}\right)\bigg{)}+\beta_{h\eta}\lambda_{H\phi}^{2}v_{H}^{2}x_{N}\lambda_{N}^{2}\cos^{2}\theta$
$\displaystyle+$ $\displaystyle
2\beta_{h\eta}\lambda_{H\phi}\lambda_{\phi}\lambda_{N}m_{N}v_{H}x_{N}\sin
2\theta\bigg{)}$

In the local $U(1)_{(B-L)}$ scenario, the reduced cross section of $NN\to\rho
Z^{\prime},Z^{\prime}Z^{\prime}$ are

$\displaystyle\hat{\sigma}(NN\to\rho Z^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{48\pi x_{N}}\bigg{(}\frac{\beta_{\rho
Z^{\prime}}x_{N}}{\lambda_{N}^{2}(x_{N}-r_{Z^{\prime}N})^{2}}\Big{(}384g_{B-L}^{6}(x_{N}-6)+3\lambda_{N}^{6}r_{\rho
N}^{2}x_{N}-12\lambda_{N}^{4}g_{B-L}^{2}$ (35) $\displaystyle\times$
$\displaystyle(2r_{\rho N}^{2}-2r_{\rho N}r_{Z^{\prime}N}+2r_{\rho
N}x_{N}-r_{Z^{\prime}N}x_{N}+2x_{N}^{2})+4\lambda_{N}^{2}g_{B-L}^{4}\big{(}3r_{\rho
N}^{2}+3r_{Z^{\prime}N}^{2}$ $\displaystyle-$ $\displaystyle
6r_{Z^{\prime}N}x_{N}-6r_{\rho
N}(x_{N}+r_{Z^{\prime}N}-8)+x_{N}(3x_{N}-\beta_{\rho
Z^{\prime}}^{2}x_{N}+60)\big{)}\Big{)}$ $\displaystyle-$
$\displaystyle\frac{64\lambda_{N}^{2}\beta_{\rho Z^{\prime}}x_{N}}{(r_{\rho
N}+r_{Z^{\prime}N}+x_{N}(\beta_{\rho Z^{\prime}}-1))(r_{\rho
N}+r_{Z^{\prime}N}-x_{N}(\beta_{\rho
Z^{\prime}}+1))}\Big{(}\lambda_{N}^{2}\big{(}r_{\rho
N}^{2}+r_{Z^{\prime}N}^{2}$ $\displaystyle+$ $\displaystyle r_{\rho
N}(6r_{Z^{\prime}N}-2x_{N})-x_{N}^{2}(\beta_{\rho
Z^{\prime}}^{2}-1)-2r_{Z^{\prime}N}(x_{N}+8)\big{)}+2g_{B-L}^{2}\big{(}96+r_{\rho
N}^{2}$ $\displaystyle+$ $\displaystyle 6r_{\rho
N}r_{Z^{\prime}N}+r_{Z^{\prime}N}^{2}+x_{N}^{2}-\beta_{\rho
Z^{\prime}}^{2}x_{N}^{2}-2r_{\rho
N}(12+x_{N})-2r_{Z^{\prime}N}(8+x_{N})\big{)}\Big{)}$ $\displaystyle+$
$\displaystyle\big{(}12\lambda_{N}^{2}g_{B-L}^{2}(r_{\rho
N}+r_{Z^{\prime}N}-x_{N}-12)+6\lambda_{N}^{4}(r_{\rho
N}+r_{Z^{\prime}N}-x_{N})\big{)}$ $\displaystyle\times$
$\displaystyle\log\left(\frac{(1-\beta_{\rho Z^{\prime}})x_{N}-r_{\rho
N}-r_{Z^{\prime}N}}{(1+\beta_{\rho Z^{\prime}})x_{N}-r_{\rho
N}-r_{Z^{\prime}N}}\right)+6\lambda_{N}^{2}\beta_{\rho
Z^{\prime}}x_{N}(\lambda_{N}^{2}-2g_{B-L}^{2})$ $\displaystyle+$
$\displaystyle\frac{12\lambda_{N}^{2}}{r_{\rho
N}+r_{Z^{\prime}N}-x_{N}}\Big{(}-\lambda_{N}^{2}r_{Z^{\prime}N}(r_{\rho
N}-4)+2g_{B-L}^{2}\big{(}r_{\rho N}(r_{Z^{\prime}N}-x_{N}+4)$ $\displaystyle-$
$\displaystyle
2(r_{Z^{\prime}N}-3x_{N}+12)\big{)}\Big{)}\log\left(\frac{(1-\beta_{\rho
Z^{\prime}})x_{N}-r_{\rho N}-r_{Z^{\prime}N}}{(1+\beta_{\rho
Z^{\prime}})x_{N}-r_{\rho
N}-r_{Z^{\prime}N}}\right)-\frac{12}{r_{Z^{\prime}N}-x_{N}}$
$\displaystyle\times$ $\displaystyle\bigg{(}-48\beta_{\rho
Z^{\prime}}g_{B-L}^{4}x_{N}+2\beta_{\rho
Z^{\prime}}\lambda_{N}^{2}g_{B-L}^{2}(r_{\rho
N}+r_{Z^{\prime}N})x_{N}+\Big{(}\lambda_{N}^{4}r_{\rho
N}(r_{Z^{\prime}N}-r_{\rho N}$ $\displaystyle+$ $\displaystyle
x_{N})+16g_{B-L}^{4}(2r_{\rho
N}+r_{Z^{\prime}N}+x_{N}-12)-4\lambda_{N}^{2}g_{B-L}^{2}\big{(}r_{\rho
N}(r_{Z^{\prime}N}+x_{N}-2)$ $\displaystyle-$ $\displaystyle
2(r_{Z^{\prime}N}+x_{N})\big{)}\Big{)}\log\left(\frac{(1-\beta_{\rho
Z^{\prime}})x_{N}-r_{\rho N}-r_{Z^{\prime}N}}{(1+\beta_{\rho
Z^{\prime}})x_{N}-r_{\rho N}-r_{Z^{\prime}N}}\right)\bigg{)}\bigg{)},$

$\displaystyle\hat{\sigma}(NN\to Z^{\prime}Z^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{32\pi
x_{N}}\bigg{(}-2\beta_{Z^{\prime}Z^{\prime}}\lambda_{N}^{2}(\lambda_{N}^{2}+4g_{B-L}^{2})x_{N}+\frac{\beta_{Z^{\prime}Z^{\prime}}(x_{N}-4)x_{N}}{(x_{N}-r_{\rho
N})^{2}}\Big{(}256g_{B-L}^{4}$ (36) $\displaystyle+$
$\displaystyle\lambda_{N}^{4}r_{\rho
N}^{2}-8\lambda_{N}^{2}g_{B-L}^{2}(r_{Z^{\prime}N}+2x_{N})\Big{)}-\frac{2\beta_{Z^{\prime}Z^{\prime}}x_{N}}{4r_{Z^{\prime}N}x_{N}-4r_{Z^{\prime}N}^{2}+(\beta_{Z^{\prime}Z^{\prime}}^{2}-1)x_{N}^{2}}$
$\displaystyle\times$
$\displaystyle\Big{(}\lambda_{N}^{4}\big{(}4r_{Z^{\prime}N}x_{N}-8r_{Z^{\prime}N}^{2}+(\beta_{Z^{\prime}Z^{\prime}}^{2}-1)x_{N}^{2}\big{)}-4\lambda_{N}^{2}g_{B-L}^{2}\big{(}8r_{Z^{\prime}N}^{2}-(\beta_{Z^{\prime}Z^{\prime}}^{2}-1)x_{N}^{2}$
$\displaystyle-$ $\displaystyle
4r_{Z^{\prime}N}(6+x_{N})\big{)}-4g_{B-L}^{4}\big{(}144+8r_{Z^{\prime}N}^{2}+x_{N}^{2}-\beta_{Z^{\prime}Z^{\prime}}^{2}x_{N}^{2}-4r_{Z^{\prime}N}(12+x_{N})\big{)}\Big{)}$
$\displaystyle+$
$\displaystyle\frac{4}{2r_{Z^{\prime}N}-x_{N}}\Big{(}\lambda_{N}^{4}r_{Z^{\prime}N}^{2}+4g_{B-L}^{4}\big{(}2r_{Z^{\prime}N}(x_{N}-1)-3(x_{N}-4)\big{)}+4\lambda_{N}^{2}g_{B-L}^{2}r_{Z^{\prime}N}$
$\displaystyle\times$
$\displaystyle\\!(r_{Z^{\prime}N}\\!-x_{N}\\!-2)\Big{)}\\!\log\left(\frac{(1-\beta_{Z^{\prime}Z^{\prime}})x_{N}\\!-r_{Z^{\prime}N}}{(1+\beta_{Z^{\prime}Z^{\prime}})x_{N}\\!-r_{Z^{\prime}N}}\right)\\!+2\Big{(}\lambda_{N}^{4}(2r_{Z^{\prime}N}\\!-x_{N})\\!+4\lambda_{N}^{2}g_{B-L}^{2}$
$\displaystyle\times$
$\displaystyle(2r_{Z^{\prime}N}-x_{N})+g_{B-L}^{4}\big{(}8r_{Z^{\prime}N}-4(x_{N}+18)\big{)}\Big{)}\log\left(\frac{(1-\beta_{Z^{\prime}Z^{\prime}})x_{N}\\!-r_{Z^{\prime}N}}{(1+\beta_{Z^{\prime}Z^{\prime}})x_{N}\\!-r_{Z^{\prime}N}}\right)$
$\displaystyle+$ $\displaystyle\frac{4}{r_{\rho
N}-x_{N}}\bigg{(}2\beta_{Z^{\prime}Z^{\prime}}x_{N}(16g_{B-L}^{4}+\lambda_{N}^{4}r_{\rho
N}-2\lambda_{N}^{2}g_{B-L}^{2}x_{N})+\Big{(}8\lambda_{N}^{2}g_{B-L}^{2}$
$\displaystyle\times$
$\displaystyle\big{(}r_{Z^{\prime}N}(x_{N}-2)-x_{N}\big{)}+\lambda_{N}^{4}r_{\rho
N}(x_{N}-2r_{Z^{\prime}N})-32g_{B-L}^{4}(r_{Z^{\prime}N}+x_{N}-6)\Big{)}$
$\displaystyle\times$
$\displaystyle\log\left(\frac{(1-\beta_{Z^{\prime}Z^{\prime}})x_{N}-r_{Z^{\prime}N}}{(1+\beta_{Z^{\prime}Z^{\prime}})x_{N}-r_{Z^{\prime}N}}\right)\bigg{)}\bigg{)}.$

As couplings of $\rho$ to SM fermions are suppressed by the small mixing
$\theta$, we only consider the $Z^{\prime}$ contribution in the process of
$NN\to f\bar{f}$. The result is

$\displaystyle\hat{\sigma}(NN\to f\bar{f})$ $\displaystyle=$
$\displaystyle\frac{N_{c}(f)Q_{f}^{2}\beta_{ff}g_{B-L}^{4}\Big{(}12r_{fN}(x_{N}-4)+x_{N}(3x_{N}+\beta_{ff}^{2}x_{N}-12)\Big{)}}{24\pi(x_{N}-r_{Z^{\prime}N})^{2}},$
(37)

where $N_{c}(f)$ and $Q_{f}$ represent the color multiplicity and the $B-L$
charge of the fermion $f$, respectively. Meanwhile, the s-wave $NN\to
hZ^{\prime}$ induced by the small mixing $\theta$ is calculated as

$\displaystyle\hat{\sigma}(NN\to hZ^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{8\pi
x_{N}m_{N}^{2}}\bigg{(}\frac{\beta_{hZ^{\prime}}x_{N}}{6\lambda_{N}^{2}(x_{N}-r_{Z^{\prime}N})^{2}}\bigg{(}4m_{N}^{2}\sin^{2}\theta\Big{(}96g_{B-L}^{6}(x_{N}-6)+3\lambda_{\phi}^{2}\lambda_{N}^{2}x_{N}$
(38) $\displaystyle+$
$\displaystyle\lambda_{N}^{2}g_{B-L}^{4}\big{(}3r_{hN}^{2}+3r_{Z^{\prime}N}^{2}-6r_{Z^{\prime}N}x_{N}-6r_{hN}(x_{N}+r_{Z^{\prime}N}-8)$
$\displaystyle+$ $\displaystyle
x_{N}(3x_{N}-\beta_{hZ^{\prime}}^{2}x_{N}+60)\big{)}+3\lambda_{N}^{2}g_{B-L}^{2}\big{(}4\lambda_{\phi}(r_{hN}-r_{Z^{\prime}N})$
$\displaystyle+$ $\displaystyle
x_{N}\lambda_{N}^{2}(r_{Z^{\prime}N}-r_{hN}-2x_{N})\big{)}\Big{)}+3\lambda_{H\phi}v_{H}\lambda_{N}^{3}\Big{(}\lambda_{H\phi}v_{H}\lambda_{N}x_{N}\cos^{2}\theta$
$\displaystyle+$ $\displaystyle 2m_{N}\sin
2\theta(2g_{B-L}^{2}(r_{hN}-r_{Z^{\prime}N})+\lambda_{\phi}x_{N})\Big{)}\bigg{)}$
$\displaystyle-$
$\displaystyle\frac{2m_{N}\sin\theta}{r_{Z^{\prime}N}-x_{N}}\bigg{(}-2\beta_{hZ^{\prime}}g_{B-L}^{2}m_{N}x_{N}\sin\theta\big{(}24g_{B-L}^{2}-\lambda_{N}^{2}(r_{hN}+r_{Z^{\prime}N})\big{)}$
$\displaystyle+$
$\displaystyle\\!\\!\Big{(}\lambda_{H\phi}v_{H}\lambda_{N}^{3}\\!\cos\theta(r_{hN}\\!-r_{Z^{\prime}N}\\!-x_{N})\\!+\\!2m_{N}\\!\sin\theta\big{(}8g_{B-L}^{4}\\!(2r_{hN}\\!+r_{Z^{\prime}N}\\!+x_{N}\\!-12)$
$\displaystyle+$
$\displaystyle\\!\lambda_{\phi}\lambda_{N}^{2}(r_{hN}\\!-r_{Z^{\prime}N}\\!-x_{N})-2g_{B-L}^{2}y_{N}^{2}(r_{hN}(r_{Z^{\prime}N}+x_{N}-2)-2r_{Z^{\prime}N}-2x_{N})\big{)}\Big{)}$
$\displaystyle\times$
$\displaystyle\log\left(\frac{(1-\beta_{hZ^{\prime}})x_{N}-r_{hN}-r_{Z^{\prime}N}}{(1+\beta_{hZ^{\prime}})x_{N}-r_{hN}-r_{Z^{\prime}N}}\right)\bigg{)}$
$\displaystyle-$
$\displaystyle\bigg{(}\frac{\lambda_{N}^{2}\beta_{hZ^{\prime}}x_{N}m_{N}^{2}\sin^{2}\theta}{(r_{hN}+r_{Z^{\prime}N}+x_{N}(\beta_{hZ^{\prime}}-1))(r_{hN}+r_{Z^{\prime}N}-x_{N}(\beta_{hZ^{\prime}}+1))}\Big{(}\lambda_{N}^{2}\big{(}r_{hN}^{2}+r_{Z^{\prime}N}^{2}$
$\displaystyle+$ $\displaystyle
r_{hN}(6r_{Z^{\prime}N}-2x_{N})-x_{N}^{2}(\beta_{hZ^{\prime}}^{2}-1)-2r_{Z^{\prime}N}(x_{N}+8)\big{)}+2g_{B-L}^{2}\big{(}96+r_{hN}^{2}$
$\displaystyle+$ $\displaystyle
6r_{hN}r_{Z^{\prime}N}+r_{Z^{\prime}N}^{2}+x_{N}^{2}-\beta_{hZ^{\prime}}^{2}x_{N}^{2}-2r_{hN}(12+x_{N})-2r_{Z^{\prime}N}(8+x_{N})\big{)}\Big{)}$
$\displaystyle-$ $\displaystyle
m_{N}^{2}\lambda_{N}^{2}\sin^{2}\theta\big{(}2g_{B-L}^{2}(r_{hN}+r_{Z^{\prime}N}-x_{N}-12)+\lambda_{N}^{2}(r_{hN}+r_{Z^{\prime}N}-x_{N})\big{)}$
$\displaystyle\times$
$\displaystyle\log\left(\frac{(1-\beta_{hZ^{\prime}})x_{N}-r_{hN}-r_{Z^{\prime}N}}{(1+\beta_{hZ^{\prime}})x_{N}-r_{hN}-r_{Z^{\prime}N}}\right)\bigg{)}$
$\displaystyle+$
$\displaystyle\bigg{(}\frac{2m_{N}^{2}\lambda_{N}^{2}\sin\theta^{2}}{r_{hN}+r_{Z^{\prime}N}-x_{N}}\Big{(}-\lambda_{N}^{2}r_{Z^{\prime}N}(r_{hN}-4)+2g_{B-L}^{2}\big{(}r_{hN}(r_{Z^{\prime}N}-x_{N}+4)$
$\displaystyle-$ $\displaystyle
2(r_{Z^{\prime}N}-3x_{N}+12)\big{)}\Big{)}\log\left(\frac{(1-\beta_{hZ^{\prime}})x_{N}-r_{hN}-r_{Z^{\prime}N}}{(1+\beta_{hZ^{\prime}})x_{N}-r_{hN}-r_{Z^{\prime}N}}\right)$
$\displaystyle+$ $\displaystyle
m_{N}^{2}\lambda_{N}^{2}\sin\theta^{2}\beta_{hZ^{\prime}}x_{N}(\lambda_{N}^{2}-2g_{B-L}^{2})\bigg{)}\bigg{)}$

In addition, for the $s$-channel processes with $\rho$ and $Z^{\prime}$
resonance, we have added the decay width $\Gamma_{\rho}$ and
$\Gamma_{Z}^{\prime}$ at the denominator of Breit-Wigner propagator in the
calculations.

As for the annihilation processes of $\chi\chi$, the reduced cross section can
refer to the corresponding one of $N$ pair annihilation, i.e., replacing
$\lambda_{N}$ and $m_{N}$ with $\lambda_{\chi}$ and $m_{\chi}$, respectively.
Meanwhile, the $\beta_{ij}$ function for $\chi$ pair annihilation is defined
as

$\displaystyle\beta_{ij}=x_{\chi}^{-1}\sqrt{(1-4x_{\chi}^{-1})(x_{\chi}^{2}+r_{i\chi}^{2}+r_{j\chi}^{2}-2x_{\chi}r_{i\chi}-2x_{\chi}r_{j\chi}-2r_{i\chi}r_{j\chi})}$
(39)

The reduced cross section of the conversion process $\chi\chi\to NN$ is

$\displaystyle\hat{\sigma}(\chi\chi\to NN)_{\rm global}$ $\displaystyle=$
$\displaystyle\frac{\lambda_{N}^{2}\lambda_{\chi}^{2}}{16\pi(x_{\chi}-r_{\rho\chi})^{2}}\beta_{NN}\Big{(}r_{\rho\chi}(r_{\rho\chi}-2x_{\chi})-4r_{N\chi}(x_{\chi}-4)+2x_{\chi}(x_{\chi}-2)\Big{)},$
(40) $\displaystyle\hat{\sigma}(\chi\chi\to NN)_{\rm local}$ $\displaystyle=$
$\displaystyle\frac{\beta_{NN}}{48\pi}\bigg{(}\frac{3\lambda_{\chi}^{2}\lambda_{N}^{2}(x_{\chi}-4)(x_{\chi}-4r_{N\chi})}{(x_{\chi}-r_{\rho\chi})^{2}}+\frac{1}{(x_{\chi}-r_{Z^{\prime}\chi})^{2}}\Big{(}\\!-24\lambda_{\chi}\lambda_{N}g_{B-L}^{2}r_{N\chi}^{1/2}x_{\chi}$
(41) $\displaystyle+$ $\displaystyle
3\lambda_{\chi}^{2}\lambda_{N}^{2}x_{\chi}^{2}+g_{B-L}^{4}\big{(}-12r_{N\chi}(x_{\chi}-8)+x_{\chi}(-12+3x_{\chi}+\beta_{NN}^{2}x_{\chi})\big{)}\Big{)}\bigg{)}.$

## References

  * (1) Y. Fukuda et al. [Super-Kamiokande], Phys. Rev. Lett. 81, 1562-1567 (1998) [arXiv:hep-ex/9807003 [hep-ex]].
  * (2) Q. R. Ahmad et al. [SNO], Phys. Rev. Lett. 89, 011301 (2002) [arXiv:nucl-ex/0204008 [nucl-ex]].
  * (3) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (4) C. S. Fong, E. Nardi and A. Riotto, Adv. High Energy Phys. 2012, 158303 (2012) [arXiv:1301.3062 [hep-ph]].
  * (5) G. Bertone, D. Hooper and J. Silk, Phys. Rept. 405, 279-390 (2005) [arXiv:hep-ph/0404175 [hep-ph]].
  * (6) M. Cirelli, A. Strumia and J. Zupan, [arXiv:2406.01705 [hep-ph]].
  * (7) T. Asaka and M. Shaposhnikov, Phys. Lett. B 620, 17-26 (2005) [arXiv:hep-ph/0505013 [hep-ph]].
  * (8) E. Ma, Mod. Phys. Lett. A 21, 1777-1782 (2006) [arXiv:hep-ph/0605180 [hep-ph]].
  * (9) A. Falkowski, J. T. Ruderman and T. Volansky, JHEP 05, 106 (2011) [arXiv:1101.4936 [hep-ph]].
  * (10) A. Falkowski, E. Kuflik, N. Levi and T. Volansky, Phys. Rev. D 99, no.1, 015022 (2019) [arXiv:1712.07652 [hep-ph]].
  * (11) T. Hugle, M. Platscher and K. Schmitz, Phys. Rev. D 98, no.2, 023020 (2018) [arXiv:1804.09660 [hep-ph]].
  * (12) M. Chianese, B. Fu and S. F. King, JCAP 03, 030 (2020) [arXiv:1910.12916 [hep-ph]].
  * (13) A. Liu, Z. L. Han, Y. Jin and F. X. Yang, Phys. Rev. D 101, no.9, 095005 (2020) [arXiv:2001.04085 [hep-ph]].
  * (14) P. Minkowski, Phys. Lett. B 67, 421-428 (1977)
  * (15) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980)
  * (16) M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45-47 (1986)
  * (17) S. Davidson, E. Nardi and Y. Nir, Phys. Rept. 466, 105-177 (2008) [arXiv:0802.2962 [hep-ph]].
  * (18) S. Blanchet and P. Di Bari, New J. Phys. 14, 125012 (2012) [arXiv:1211.0512 [hep-ph]].
  * (19) N. S. Manton, Phys. Rev. D 28, 2019 (1983)
  * (20) F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 2212 (1984)
  * (21) V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. B 155, 36 (1985)
  * (22) S. Y. Khlebnikov and M. E. Shaposhnikov, Nucl. Phys. B 308, 885-912 (1988)
  * (23) J. A. Harvey and M. S. Turner, Phys. Rev. D 42, 3344-3349 (1990)
  * (24) R. N. Mohapatra, S. Antusch, K. S. Babu, G. Barenboim, M. C. Chen, S. Davidson, A. de Gouvêa, P. de Holanda, B. Dutta and Y. Grossman, et al. Rept. Prog. Phys. 70, 1757-1867 (2007) [arXiv:hep-ph/0510213 [hep-ph]].
  * (25) M. Drewes, Int. J. Mod. Phys. E 22, 1330019 (2013) [arXiv:1303.6912 [hep-ph]].
  * (26) S. Davidson and A. Ibarra, Phys. Lett. B 535, 25-32 (2002) [arXiv:hep-ph/0202239 [hep-ph]].
  * (27) M. Flanz, E. A. Paschos, U. Sarkar and J. Weiss, Phys. Lett. B 389, 693-699 (1996) [arXiv:hep-ph/9607310 [hep-ph]].
  * (28) A. Pilaftsis, Nucl. Phys. B 504, 61-107 (1997) [arXiv:hep-ph/9702393 [hep-ph]].
  * (29) A. Pilaftsis, Phys. Rev. D 56, 5431-5451 (1997) [arXiv:hep-ph/9707235 [hep-ph]].
  * (30) A. Pilaftsis and T. E. J. Underwood, Nucl. Phys. B 692, 303-345 (2004) [arXiv:hep-ph/0309342 [hep-ph]].
  * (31) B. Dev, M. Garny, J. Klaric, P. Millington and D. Teresi, Int. J. Mod. Phys. A 33, 1842003 (2018) [arXiv:1711.02863 [hep-ph]].
  * (32) F. F. Deppisch, P. S. Bhupal Dev and A. Pilaftsis, New J. Phys. 17, no.7, 075019 (2015) [arXiv:1502.06541 [hep-ph]].
  * (33) Y. Cai, T. Han, T. Li and R. Ruiz, Front. in Phys. 6, 40 (2018) [arXiv:1711.02180 [hep-ph]].
  * (34) A. de Gouvea and P. Vogel, Prog. Part. Nucl. Phys. 71, 75-92 (2013) [arXiv:1303.4097 [hep-ph]].
  * (35) S. Alekhin, W. Altmannshofer, T. Asaka, B. Batell, F. Bezrukov, K. Bondarenko, A. Boyarsky, K. Y. Choi, C. Corral and N. Craig, et al. Rept. Prog. Phys. 79, no.12, 124201 (2016) [arXiv:1504.04855 [hep-ph]].
  * (36) S. Dodelson and L. M. Widrow, Phys. Rev. Lett. 72, 17-20 (1994) [arXiv:hep-ph/9303287 [hep-ph]].
  * (37) A. Datta, R. Roshan and A. Sil, Phys. Rev. Lett. 127, no.23, 231801 (2021) [arXiv:2104.02030 [hep-ph]].
  * (38) 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]].
  * (39) E. Ma, Phys. Rev. D 73, 077301 (2006) [arXiv:hep-ph/0601225 [hep-ph]].
  * (40) J. Kubo, E. Ma and D. Suematsu, Phys. Lett. B 642, 18-23 (2006) [arXiv:hep-ph/0604114 [hep-ph]].
  * (41) M. Pospelov, A. Ritz and M. B. Voloshin, Phys. Lett. B 662, 53-61 (2008) [arXiv:0711.4866 [hep-ph]].
  * (42) V. González-Macías, J. I. Illana and J. Wudka, JHEP 05, 171 (2016) [arXiv:1601.05051 [hep-ph]].
  * (43) M. Escudero, N. Rius and V. Sanz, Eur. Phys. J. C 77, no.6, 397 (2017) [arXiv:1607.02373 [hep-ph]].
  * (44) B. Batell, T. Han, D. McKeen and B. Shams Es Haghi, Phys. Rev. D 97, no.7, 075016 (2018) [arXiv:1709.07001 [hep-ph]].
  * (45) M. Blennow, E. Fernandez-Martinez, A. Olivares-Del Campo, S. Pascoli, S. Rosauro-Alcaraz and A. V. Titov, Eur. Phys. J. C 79, no.7, 555 (2019) [arXiv:1903.00006 [hep-ph]].
  * (46) P. Ballett, M. Hostert and S. Pascoli, Phys. Rev. D 101, no.11, 115025 (2020) [arXiv:1903.07589 [hep-ph]].
  * (47) A. Biswas, D. Borah and D. Nanda, JCAP 10, 002 (2021) [arXiv:2103.05648 [hep-ph]].
  * (48) R. Coy, A. Gupta and T. Hambye, Phys. Rev. D 104, no.8, 083024 (2021) [arXiv:2104.00042 [hep-ph]].
  * (49) D. Borah, M. Dutta, S. Mahapatra and N. Sahu, Phys. Rev. D 105, no.1, 015004 (2022) [arXiv:2110.00021 [hep-ph]].
  * (50) B. Fu and S. F. King, JHEP 12, 121 (2021) [arXiv:2110.00588 [hep-ph]].
  * (51) A. Das, S. Gola, S. Mandal and N. Sinha, Phys. Lett. B 829, 137117 (2022) [arXiv:2202.01443 [hep-ph]].
  * (52) L. Coito, C. Faubel, J. Herrero-García, A. Santamaria and A. Titov, JHEP 08, 085 (2022) [arXiv:2203.01946 [hep-ph]].
  * (53) A. Biswas, D. Borah, N. Das and D. Nanda, [arXiv:2205.01144 [hep-ph]].
  * (54) A. Liu, Z. L. Han, Y. Jin and H. Li, Phys. Rev. D 108, no.7, 075021 (2023) [arXiv:2306.14091 [hep-ph]].
  * (55) B. Batell, T. Han and B. Shams Es Haghi, Phys. Rev. D 97, no.9, 095020 (2018) [arXiv:1704.08708 [hep-ph]].
  * (56) E. Aprile et al. [XENON], Phys. Rev. Lett. 121, no.11, 111302 (2018) [arXiv:1805.12562 [astro-ph.CO]].
  * (57) E. Aprile et al. [XENON], Phys. Rev. Lett. 131, no.4, 041003 (2023) [arXiv:2303.14729 [hep-ex]].
  * (58) X. Cui et al. [PandaX-II], Phys. Rev. Lett. 119, no.18, 181302 (2017) [arXiv:1708.06917 [astro-ph.CO]].
  * (59) Y. Meng et al. [PandaX-4T], Phys. Rev. Lett. 127, no.26, 261802 (2021) [arXiv:2107.13438 [hep-ex]].
  * (60) A. Boveia and C. Doglioni, Ann. Rev. Nucl. Part. Sci. 68, 429-459 (2018) [arXiv:1810.12238 [hep-ex]].
  * (61) W. Adam and I. Vivarelli, Int. J. Mod. Phys. A 37, no.02, 2130022 (2022) [arXiv:2111.10180 [hep-ex]].
  * (62) M. Aguilar et al. [AMS], Phys. Rev. Lett. 110, 141102 (2013)
  * (63) M. Ackermann et al. [Fermi-LAT], Phys. Rev. Lett. 115, no.23, 231301 (2015) [arXiv:1503.02641 [astro-ph.HE]].
  * (64) R. E. Marshak and R. N. Mohapatra, Phys. Lett. B 91, 222-224 (1980)
  * (65) R. N. Mohapatra and R. E. Marshak, Phys. Rev. Lett. 44, 1316-1319 (1980) [erratum: Phys. Rev. Lett. 44, 1643 (1980)]
  * (66) N. Okada and O. Seto, Phys. Rev. D 82, 023507 (2010) [arXiv:1002.2525 [hep-ph]].
  * (67) M. Escudero, N. Rius and V. Sanz, JHEP 02, 045 (2017) [arXiv:1606.01258 [hep-ph]].
  * (68) S. Okada, Adv. High Energy Phys. 2018, 5340935 (2018) [arXiv:1803.06793 [hep-ph]].
  * (69) S. Iso, N. Okada and Y. Orikasa, Phys. Rev. D 83, 093011 (2011) [arXiv:1011.4769 [hep-ph]].
  * (70) J. Heeck and D. Teresi, Phys. Rev. D 94, no.9, 095024 (2016) [arXiv:1609.03594 [hep-ph]].
  * (71) P. S. B. Dev, R. N. Mohapatra and Y. Zhang, JHEP 03, 122 (2018) [arXiv:1711.07634 [hep-ph]].
  * (72) C. Garcia-Cely, A. Ibarra and E. Molinaro, JCAP 02, 032 (2014) [arXiv:1312.3578 [hep-ph]].
  * (73) C. Garcia-Cely, A. Ibarra and E. Molinaro, JCAP 11, 061 (2013) [arXiv:1310.6256 [hep-ph]].
  * (74) J. A. Casas and A. Ibarra, Nucl. Phys. B 618, 171-204 (2001) [arXiv:hep-ph/0103065 [hep-ph]].
  * (75) B. Barman, P. S. Bhupal Dev and A. Ghoshal, Phys. Rev. D 108, no.3, 035037 (2023) [arXiv:2210.07739 [hep-ph]].
  * (76) Y. Chikashige, R. N. Mohapatra and R. D. Peccei, Phys. Lett. B 98, 265-268 (1981)
  * (77) P. S. B. Dev, R. N. Mohapatra and Y. Zhang, JHEP 02, 186 (2016) [arXiv:1512.08507 [hep-ph]].
  * (78) F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz and S. Vogl, JHEP 02, 016 (2016) [arXiv:1510.02110 [hep-ph]].
  * (79) M. Duerr, F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz and S. Vogl, JHEP 09, 042 (2016) [arXiv:1606.07609 [hep-ph]].
  * (80) M. Plumacher, Z. Phys. C 74, 549-559 (1997) [arXiv:hep-ph/9604229 [hep-ph]].
  * (81) D. Borah, A. Dasgupta and D. Mahanta, Phys. Rev. D 104, no.7, 075006 (2021) [arXiv:2106.14410 [hep-ph]].
  * (82) J. Aalbers et al. [LZ], Phys. Rev. Lett. 131, no.4, 041002 (2023) [arXiv:2207.03764 [hep-ex]].
  * (83) T. Robens and T. Stefaniak, Eur. Phys. J. C 76, no.5, 268 (2016) [arXiv:1601.07880 [hep-ph]].
  * (84) T. Robens, Springer Proc. Phys. 292, 141-152 (2023) [arXiv:2209.15544 [hep-ph]].
  * (85) [ATLAS], ATLAS-CONF-2020-052.
  * (86) D. S. Akerib et al. [LZ], [arXiv:1509.02910 [physics.ins-det]].
  * (87) G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, Comput. Phys. Commun. 185, 960-985 (2014) [arXiv:1305.0237 [hep-ph]].
  * (88) J. Billard, L. Strigari and E. Figueroa-Feliciano, Phys. Rev. D 89, no.2, 023524 (2014) [arXiv:1307.5458 [hep-ph]].
  * (89) Y. Tan, X. Shi, R. Kiuchi, M. Ruan, M. Jing, D. Yu, K. Zhang, X. Lou, X. Mo and G. Li, et al. Chin. Phys. C 44, no.12, 123001 (2020) [arXiv:2001.05912 [hep-ex]].
  * (90) S. R. Gozzini et al. [ANTARES], PoS ICRC2023, 1375 (2023)
  * (91) M. D. Campos, F. S. Queiroz, C. E. Yaguna and C. Weniger, JCAP 07, 016 (2017) [arXiv:1702.06145 [hep-ph]].
  * (92) A. Acharyya et al. [CTA], JCAP 01, 057 (2021) [arXiv:2007.16129 [astro-ph.HE]].
  * (93) A. Mangipudi, E. Thrane and C. Balazs, JCAP 11, 010 (2022) [arXiv:2112.10371 [astro-ph.HE]].
  * (94) M. Lattanzi, R. A. Lineros and M. Taoso, New J. Phys. 16, no.12, 125012 (2014) [arXiv:1406.0004 [hep-ph]].
  * (95) G. Aad et al. [ATLAS], Phys. Lett. B 796, 68-87 (2019) [arXiv:1903.06248 [hep-ex]].
  * (96) A. M. Sirunyan et al. [CMS], JHEP 07, 208 (2021) [arXiv:2103.02708 [hep-ex]].
  * (97) G. Cacciapaglia, C. Csaki, G. Marandella and A. Strumia, Phys. Rev. D 74, 033011 (2006) [arXiv:hep-ph/0604111 [hep-ph]].
  * (98) S. Schael et al. [ALEPH, DELPHI, L3, OPAL and LEP Electroweak], Phys. Rept. 532, 119-244 (2013) [arXiv:1302.3415 [hep-ex]].
  * (99) C. Blanco, M. Escudero, D. Hooper and S. J. Witte, JCAP 11, 024 (2019) [arXiv:1907.05893 [hep-ph]].
  * (100) L. Basso, A. Belyaev, S. Moretti and C. H. Shepherd-Themistocleous, Phys. Rev. D 80, 055030 (2009) [arXiv:0812.4313 [hep-ph]].
  * (101) S. Iso, N. Okada and Y. Orikasa, Phys. Rev. D 80, 115007 (2009) [arXiv:0909.0128 [hep-ph]].
  * (102) P. Fileviez Perez, T. Han and T. Li, Phys. Rev. D 80, 073015 (2009) [arXiv:0907.4186 [hep-ph]].
  * (103) A. Tumasyan et al. [CMS], JHEP 11, 181 (2023) [arXiv:2307.06959 [hep-ex]].
  * (104) K. Fujii et al. [LCC Physics Working Group], [arXiv:1908.11299 [hep-ex]].
  * (105) S. K.A., A. Das, G. Lambiase, T. Nomura and Y. Orikasa, [arXiv:2308.14483 [hep-ph]].
  * (106) D. Hooper, Phys. Rev. D 91, 035025 (2015) [arXiv:1411.4079 [hep-ph]].
  * (107) J. Aalbers et al. [DARWIN], JCAP 11, 017 (2016) [arXiv:1606.07001 [astro-ph.IM]].
  * (108) A. Belyaev, N. D. Christensen and A. Pukhov, Comput. Phys. Commun. 184, 1729-1769 (2013) [arXiv:1207.6082 [hep-ph]].