11institutetext: Department of Physics and Institute of Theoretical Physics,
Nanjing Normal University, Nanjing, 210023, China22institutetext: Department
of Physics, Yantai University, Yantai 264005, China33institutetext: Department
of Physics, Chung-Ang University, Seoul 06974, Korea

# Axion Dark Radiation and Late Decaying Dark Matter in Neutrino Experiment

Yuchao Gu 1 Lei Wu 2,3 Bin Zhu
guyc@njnu.edu.cn;leiwu<EMAIL_ADDRESS>
###### Abstract

In this work, we investigate the dark sector of a supersymmetric axion model,
which consists of the late-decaying gravitino/axino dark matter and axion dark
radiation. In the early universe, the decay of the scalar superpartner of the
axion (saxion) will produce a large amount of entropy. This can not only
dilute the relic density of the gravitino/axino dark matter to avoid
overclosing the universe but also relax the constraint on the reheating
temperature $T_{R}$ after inflation. Meanwhile the axion dark radiation from
the saxion decay will increase the effective number of neutrino species
$N_{\rm eff}$, which can help to reduce the cosmological Hubble tension. In
the late universe, the decay of long-lived gravitino/axino dark matter will
produce the axions with MeV-scale kinetic energy, which can be tested through
the inverse Primakoff process $a+A\to\gamma+A$ in the neutrino experiments,
such as Hyper-Kamiokande.

##  1 Introduction

The Standard Model (SM) of particle physics is a successful theory that
describes the currently known elemental particles and their interactions.
Despite of its great success, there are still several problems in the Standard
Model, such as the absence of Dark Matter (DM) candidates and strong CP
problem Peccei and Quinn (1977). Supersymmetry is one of the most promising
extensions of the SM that can address these fundamental problems. When
$R$-parity is conserved, the lightest supersymmetric particle (LSP) can play
the role of DM particle naturally. In addition, the supersymmetric axion
models based on Peccei-Quinn mechanism also provide an elegant solution to the
strong CP problem Kim (1979); Shifman  _et al._ (1980); Zhitnitsky (1980);
Dine  _et al._ (1981); Kim and Nilles (1984); Kim and Carosi (2010).

As is known, the lightest neutralino as Weakly Interacting Massive Dark Matter
(WIMP) has been extensively studied in the past decades. However, the null
results of searching for WIMP dark matter have put stringent constraints on
its mass and interaction strength Cao  _et al._ (2014); Duan  _et al._
(2018a); Abdughani  _et al._ (2018); Duan  _et al._ (2018b); Abdughani and Wu
(2020); Abdughani  _et al._ (2019, 2021). This strongly motivates the
theoretical and experimental community to search for light or super-weakly
interacting DM particles. In the supersymmetric models, the gravitino (the
superpartner of graviton) or axino (the superpartner of axion) can serve as
such a kind of feeble DM particles (super-WIMP) Covi  _et al._ (1999, 2001);
Covi and Kim (2009); Bae  _et al._ (2015); Co  _et al._ (2017a); Hamaguchi
_et al._ (2017). Depending on the SUSY breaking scheme, the masses of these DM
particles can be very small. Their couplings with the SM particles are
extremely weak because of the large suppression of Planck scale $M_{P}$ or the
PQ symmetry breaking scale $V_{\rm PQ}$.

More interestingly, these super-WIMPs may have connection with some tentative
cosmological anomalies, such as the Hubble tension Hooper  _et al._ (2012); Gu
_et al._ (2020); Anchordoqui (2021), which refers to the tension that the
value of Hubble constant $H_{0}=(67.36\pm 0.54)$ km s-1Mpc-1 Verde  _et al._
(2019) inferred from the cosmic microwave background (CMB) in early universe
shows 4$\sigma$ to 6$\sigma$ discrepancy from that derived from the local
distance ladder in the late universe, $H_{0}=(74.03\pm 1.42)$ km s-1 Mpc-1
Riess  _et al._ (2019). This is challenging for the standard cosmology model
$\Lambda$CDM. In order to solve this problem, many new physics models have
been proposed in the past few years Di Valentino  _et al._ (2021). Among them,
the dark radiation Binder  _et al._ (2018); Pandey  _et al._ (2020); Vattis
_et al._ (2019); Xiao  _et al._ (2020); Blinov  _et al._ (2020) is one of
possible ways to reduce the tension, although a perfect solution is still
under investigation. The extra relativistic degrees of freedom can increase
the effective number of neutrino species $N_{\rm eff}$ and thus reduce the
sound horizon $r_{s}(z_{LS})$ during the time up to the last scattering. While
the angular scale $\theta_{\star}=r_{s}(z_{LS})/D_{A}(z_{LS})$ is fixed by the
measured acoustic peaks from the CMB power spectrum, where $D_{A}$ is the
angular diameter. Therefore, decreasing $D_{A}$ can keep $\theta_{\star}$
unchanged, which can be achieved by simultaneously increasing $H_{0}$.

In this paper, we investigate the phenomenology of the dark sector in a
supersymmetric axion model, where the gravitino/axino play the role of late-
decaying DM. In general, there will exist a serious cosmological
gravitino/axino problem of overclosure via thermal production Pagels and
Primack (1982); Moroi  _et al._ (1993); Asaka  _et al._ (2000); Bolz  _et al._
(2001); Roszkowski  _et al._ (2005); Pradler and Steffen (2007a, b); Cheung
_et al._ (2012); Co  _et al._ (2017b); Gu  _et al._ (2020). But thanks to the
saxion decay, a large amount of entropy can be produced to dilute the relic
density of the gravitino/axino dark matter in our model. Besides, the axions
from the saxion decay in the early universe can serve as dark radiation, which
can increase the effective neutrino species and thus relax the Hubble tension.
However, it should be mentioned too much dark radiation produced in the early
universe may affect predictions of the CMB power spectrum and big bang
nucleosynthesis(BBN) since extra dark radiation will accelerate the expansion
rate of the universe. Furthermore, if the late-decaying DM decays to the LSP
plus the SM particle directly, it may result in the non-thermal nuclear
reactions during or after BBN Kawasaki  _et al._ (2005), which can also
destroy the perfect predictions of light elements abundance. But in our
scenario, if the gravitino (axino) is the NLSP, it will decay into axino
(gravitino) LSP and axion so that the strong constraints from BBN can be
escaped. More interestingly, the axion from the decay of gravitino/axino in
the late universe can be tested through the absorptive effects in the new
direct detection experiments Hochberg  _et al._ (2017); Cui  _et al._ (2018);
Kuo  _et al._ (2021) or through the inverse Primakoff scattering Avignone  _et
al._ (1998); Bernabei  _et al._ (2001); Dent  _et al._ (2020) in neutrino
experiments, such as Hyper-Kamiokande Abe  _et al._ (2018).

This paper is laid out as follow. In Section 2, we briefly introduce our SUSY
DFSZ axion model and relevant interactions. In Section 3, we discuss the
saxion decay and implication for the Hubble tension. In Section 4, we will
consider the dilution effect originating from saxion decay and calculate the
relic density of gravitino/axino DM. In section 5, we will investigate the
possibility of detecting the axion with MeV kinetic energy from the late
decaying gravitino DM through the inverse Primakoff scattering in neutrino
experiments. Finally, we draw some conclusions in Section 6.

##  2 SUSY DFSZ axion model

In the supersymmetric axion model, the axion superfield at the energy scale
below the PQ symmetry breaking scale $V_{PQ}$ is given by

$\hat{A}=\frac{s+ia}{\sqrt{2}}+\sqrt{2}\theta\tilde{a}+\theta^{2}F,$ (1)

where $s$ is the saxion field, $a$ is the axion field and $\tilde{a}$ is the
fermionic superpartner of axion i.e. axino. The advantage of SUSY DFSZ model
is that the $\mu$ problem and the strong CP problem can be solved
simultaneously. Since two Higgs doublets carry the PQ charge, the bare $\mu$
term is forbidden. Only the higher dimensional operator is allowed,

$W_{\rm
DFSZ}\ni\lambda\frac{{\hat{S}}^{2}}{M_{P}}\hat{H}_{u}\hat{H}_{d}=\mu\hat{H}_{u}\hat{H}_{d}.$
(2)

When the PQ symmetry is broken, the singlet superfield acquires the vacuum
expectation value $\langle S\rangle\sim V_{PQ}$. Then the higgsino can obtain
the mass $\mu\sim\lambda V_{PQ}^{2}/M_{pl}$. It is interesting that a natural
$\mu$ term with $\mathcal{O}(10^{2})$ GeV can be obtained by assuming
$V_{PQ}\sim\mathcal{O}(10^{10})$ GeV and $\lambda\sim\mathcal{O}(1)$.
Therefore, the interaction of the axion superfield and the Higgs fields can be
expressed by the spurion analysis from the following superpotential,

$\mathcal{L}_{\rm DFSZ}=\int
d^{2}\theta\,\mu\exp\bigg{[}c_{\mu}\frac{\hat{A}}{V_{PQ}}\bigg{]}\hat{H}_{u}\hat{H}_{d}\sim\int
d^{2}\theta\mu\hat{H}_{u}\hat{H}_{d}+\mu\frac{c_{\mu}\hat{A}}{V_{PQ}}\hat{H}_{u}\hat{H}_{d}+\cdots.$
(3)

Here $c_{\mu}$ is the PQ charge of Higgs doublets bi-linear operator, which is
usually taken as $c_{\mu}=2$ because the Higgs field is charged by $+1$. In
addition, the interactions of saxion-axion-axino and the saxion self-
interaction come from the Kähler potential

$K=\sum_{i}v_{i}^{2}{\rm
exp}\bigg{[}c_{i}\left(\frac{\hat{A}+\hat{A}^{\dagger}}{V_{PQ}}\right)\bigg{]}$
(4)

By expanding Eq. 4 up to the cubic term, we can obtain the relevant the saxion
interactions in our study,

$\mathcal{L}_{saa,s\tilde{a}\tilde{a}}=-\frac{\kappa}{\sqrt{2}V_{PQ}}s\partial_{\mu}a\partial^{\mu}a-\frac{\kappa}{\sqrt{2}V_{PQ}}s\partial_{\mu}s\partial^{\mu}s+\frac{\kappa}{\sqrt{2}V_{PQ}}sm_{\tilde{a}}\left(\tilde{a}\tilde{a}+\tilde{a}^{\dagger}\tilde{a}^{\dagger}\right),$
(5)

where $\kappa=\sum_{i}c_{i}^{3}v_{i}^{2}/V_{PQ}^{2}$ is dimensionless
coupling. We assume a single PQ charge and take $\kappa=1$ in our calculation.

Figure 1: The lifetimes of the gravitino/axino on the plane of mass ratio $R$
versus the gravitino mass $m_{3/2}$ and the axino mass $m_{\tilde{a}}$,
respectively.

Besides, the gravitino exists naturally in SUSY axion model, which is the
prediction of the local supersymmetry. Its interactions with the SM particles
is given by,

$\mathcal{L}_{\tilde{G}}=-\frac{i}{2M_{P}}\mathcal{J}^{\mu}\tilde{G}+h.c.$ (6)

with

$\mathcal{J}^{\mu}=\sqrt{2}\sigma^{\nu}\bar{\sigma}^{\mu}\chi\partial_{\nu}\phi^{\dagger}-\sigma^{\alpha\beta}\sigma^{\mu}\lambda^{a\dagger}F_{\alpha\beta}^{a}$
(7)

where the supercurrent $\mathcal{J^{\mu}}$ includes the gauge supermultiplets
$(V,\lambda)$ and the chiral supermultiplets $(\phi,\chi)$. This will induce
the interaction of the gravitino with the axion supermultiplet. Then, we can
derive the interaction of the gravitino, the axino and the axion,

$\mathcal{L}=-\frac{1}{2M_{P}}\partial_{\nu}a\bar{\psi}_{\mu}\gamma^{\nu}\gamma^{\mu}i\gamma_{5}\tilde{a}$
(8)

where the $\psi_{\mu}$ is the gravitino field. If the gravitino is heavier
than the axino, the decay channel $\tilde{G}\rightarrow a\tilde{a}$ is allowed
and the corresponding decay width is given by

$\Gamma_{3/2}=\frac{m_{3/2}^{3}}{192\pi M_{P}^{2}}(1-R)^{2}(1-R^{2})^{3}$ (9)

where $R\equiv m_{\tilde{a}}/m_{3/2}$ is the mass ratio. On the contrary, the
axino will decay to the gravitino via the process $\tilde{a}\to\tilde{G}a$
provided that $m_{\tilde{a}}$ is larger than $m_{3/2}$. In this case, the
decay width of the axino is given by

$\Gamma_{\tilde{a}}=\frac{m_{\tilde{a}}^{5}}{96\pi
m_{3/2}^{2}M_{P}^{2}}\left(1-R^{-1}\right)^{2}\left(1-R^{-2}\right)^{3}$ (10)

In Figure 1, we show the numerical results of the lifetimes of the gravitino
($\tau_{3/2}$) and the axino ($\tau_{\tilde{a}}$) for above two cases on the
plane of the mass ratio $R$ versus the gravitino mass $m_{3/2}$ and axino mass
$m_{\tilde{a}}$ respectively. The green, red and cyan lines illustrate the
lifetime $\tau=10^{12}$, $4.3\times 10^{17}$ and $10^{22}$s, respectively. It
can be seen that the lifetime of the gravitino/axino will increase with the
decrease of its mass for a given value of the mass ratio $R$. When the
gravitino/axino mass approaches to the LSP mass, i.e. $R\rightarrow 1$, their
lifetimes will be independent of the their masses. On the other hand, when
$R\rightarrow\infty$, the lifetime of the axino not only depends on its mass
$m_{\tilde{a}}$ but also the mass ratio $R$. While if $R\rightarrow 0$, the
lifetime of the gravitino only depends on its mass $m_{3/2}$. These features
will lead to the restriction on our axion from the decay of these late-
decaying DMs.

##  3 Saxion Decay and Hubble Tension

As discussed in Sec. 2, the self-interaction of the axion supermultiplet is
induced by the Kähler potential term. Besides, the axion supermultiplet can
interact with the SM particles through the Higgs field becasue the two Higgs
doublets $H_{u}$ and $H_{d}$ carry the PQ charge. Therefore, the saxion can
decay into axion pairs, axino pairs and visible particles, whose decay widths
are given by

$\displaystyle\Gamma_{s\rightarrow aa}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}m_{s}^{3}}{64\pi V_{PQ}^{2}},$ (11)
$\displaystyle\Gamma_{s\rightarrow\tilde{a}\tilde{a}}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}m_{\tilde{a}}^{2}m_{s}}{8\pi V_{PQ}^{2}},$ (12)
$\displaystyle\Gamma_{s\rightarrow XX}$ $\displaystyle=$
$\displaystyle\mathcal{D}\frac{c_{\mu}^{2}\mu^{4}}{16\pi m_{s}V_{PQ}^{2}}.$
(13)

where $\mathcal{D}$ is the number of final state particles, such as $D=4$ in
the SM ($X=$ $h$, $Z$ and $W^{\pm}$) and $D=8$ in the MSSM ($X=$ $h$, $H$,
$A$, $H^{\pm}$, $Z$ and $W^{\pm}$). $c_{\mu}$ represents the PQ charge of
$\mu$ term. In Fig 2, we present the branching ratios of the saxion decay. It
can be seen that the branching ratio of $s\to aa$ is close to one when
$m_{s}\leq 10^{3}$ GeV, while that of $s\to XX$ will be dominant when
$m_{s}\geq 3\times 10^{3}$ GeV. This is because that the decay width of $s\to
aa$ is proportional to $m_{s}^{3}$ while that of $s\to XX$ is proportional to
$m^{-1}_{s}$. In the region $10^{3}$ GeV $<m_{s}<$ $3\times 10^{3}$ GeV, these
two branching ratios are comparable, but still larger than the branching ratio
of the saxion decay into axino pairs, $s\to\tilde{a}\tilde{a}$.

Figure 2: The branching ratios of the saxion as the function of saxion mass.
Here we consider the decay processes: $s\to aa$ (cyan line),
$s\to\tilde{a}\tilde{a}$ (magenta line) and $s\to XX$ (green line). Other
parameters are taken as $\kappa=1$, $c_{\mu}=2$, $m_{\tilde{a}}=10$ GeV,
$\mu=700$ GeV and $\mathcal{D}=4$ (The additional non-SM Higgs bosons are
decoupled.)

Due to the saxion decay, it should be mentioned that there will be extra
radiation energy from the axion, which can be injected into the early universe
in the form of dark radiation. During the radiation dominant epoch, the dark
radiation can alter the expansion rate of the universe via increasing the
radiation energy in the Friedmann equation in flat universe, which is related
with the Hubble constant by

$H^{2}(t)\simeq\frac{8\pi G}{3}(\rho_{\gamma}+\rho_{\nu}+\rho_{a})$ (14)

where $\rho_{\gamma}$, $\rho_{\nu}$ and $\rho_{a}$ is the energy density of
photon, neutrino and axion, respectively. The extra axion radiation will
contribute to the effective number of neutrino species $N_{\rm eff}$, which
has a positive relation to the Hubble constant. The additional $\Delta N_{\rm
eff}$ is parameterized by

$\Delta N_{\rm eff}=N_{\rm eff}\frac{\rho_{a}}{\rho_{\nu}},$ (15)

where $N_{\rm eff}=3.046$ in the SM. When the neutrinos were in the thermal
bath, the neutrino energy density $\rho_{\nu}$ can be related to the Standard
Model energy density $\rho_{\rm SM}$

$\rho_{\nu}(T)=\frac{g_{\nu}(T)}{g_{\rm SM}(T)}\rho_{\rm SM}(T)$ (16)

where $g_{\nu}(T)$ and $g_{\rm SM}(T)$ are the number of relativistic degree
of freedom of neutrinos and Standard Model particles at the temperature $T$,
respectively. Since the decay width $\Gamma_{s\rightarrow\tilde{a}\tilde{a}}$
is much smaller than $\Gamma_{s\rightarrow aa}$ and $\Gamma_{s\rightarrow{\rm
visible}}$, the energy density of saxion are mostly injected in the form of
axion and visible particles via saxion decay. Assuming that the additional
energy density comes from the axion and the SM particles in the saxion decay,
the ratio of the axion energy density to the SM energy density equals to the
ratio of their decay widths Co  _et al._ (2016),

$\frac{\rho_{a}(T)}{\rho_{SM}(T)}=\frac{\Gamma_{s\to aa}}{\Gamma_{s\to XX}}$
(17)

Therefore, the $\Delta N_{\rm eff}$ arising from the saxion decay is given by

$\Delta N_{\rm eff}=\frac{4}{7}\frac{\Gamma_{s\rightarrow
aa}}{\Gamma_{s\rightarrow{\rm XX}}}g_{\ast}\left(T_{\nu}^{\rm dec}\right),$
(18)

where the $g_{\ast}\left(T_{\nu}^{\rm dec}\right)$ is the number of
relativistic degree of freedom at the neutrino decoupling temperature. Since
both decay widths are proportional to the PQ symmetry breaking scale
$1/V_{PQ}^{2}$, the value of $\Delta N_{\rm eff}$ is independent of $V_{PQ}$,
but relies on $m_{s}$. In addition, the relation between the current Hubble
constant $H_{0}$ and $H_{\rm CMB}$ inferred by CMB data is Vagnozzi (2020);
Anchordoqui (2021)

$H_{0}=H_{\rm CMB}+6.2\Delta N_{\rm eff}$ (19)

with $H_{\rm CMB}=67.9$km/s/Mpc.

Figure 3: The current Hubble constant $H_{0}$ as the function of the saxion
mass $m_{s}$. The gray, light-gray and blue bands are 1$\sigma$, 2$\sigma$ and
3$\sigma$ ranges of $H_{0}$ for solving the Hubble tension, respectively. The
magenta and cyan dashed lines are the bounds on the $\Delta N_{\rm eff}$ from
CMB and BBN.

In Fig 3, we show the current Hubble constant $H_{0}$ as the function of the
saxion mass $m_{s}$. Since the heavier saxion can produce the more axion dark
radiation, the value of $H_{0}$ increases with the saxion mass $m_{s}$. We can
see that the Hubble tension will be settled within $1\sigma$ range when 1.2
TeV $<m_{s}<$ 1.3 TeV. However, it should be noted that the extra dark
radiation can accelerate the expansion rate of the universe, which will
results in earlier freeze-out with a larger neutron fraction and earlier onset
of nucleosynthesis. Furthermore, this eventually leads to a larger 4He
abundance, destroying the perfect predictions of standard BBN theory. In
addition, the extra radiation also delays the matter-radiation equality,
altering the observed CMB spectrum. Therefore, the extra effective number of
neutrino species $\Delta N_{\rm eff}$ arising from the extra dark radiation is
strongly constrained by both of the CMB and the BBN, which further produce a
bound on the saxion mass via Eq. 18. It can be seen that the current BBN bound
$\Delta N_{\rm eff}<0.5$ Blinov  _et al._ (2019) is weaker than the CMB and
large scale structure bound $\Delta N_{\rm eff}<0.28$ Aghanim  _et al._
(2020). But the BBN limit is not as sensitive to the choice of cosmological
model. Under this constraint, the Hubble tension is only marginally reduced to
$2\sigma$ via simply increasing $\Delta N_{\rm eff}$. While the currently
allowed range $10^{-2}<\Delta N_{\rm eff}<0.5$ will be within the reach of
CMB-S4 predictions Abazajian  _et al._ (2016).

##  4 Late Decaying DM: Gravitino/Axino

In supersymmetric axion models, the gravitino/axino can play the role of dark
matter particle. The gravitino mass $m_{3/2}$ depends on the scheme of
supersymmetry breaking Pagels and Primack (1982); Weinberg (1982); Khlopov and
Linde (1984); Dine  _et al._ (1995); Giudice and Rattazzi (1999); Randall and
Sundrum (1999); Ellis  _et al._ (2004); Buchmuller  _et al._ (2006); Dudas
_et al._ (2017), which can be at the same order of the axino mass
$m_{\tilde{a}}$ Goto and Yamaguchi (1992); Chun  _et al._ (1992); Chun and
Lukas (1995); Kim and Seo (2012). Therefore, they can consist of a two-
component dark matter sector, in which the heavier one is considered as late
decaying DM and the lighter one is the LSP dark matter. However, if the mass
of the gravitino is heavier than keV, the abundance of the gravitino will
easily overclose the universe. In order to give the present observed dark
matter relic density, the reheating temperature $T_{R}$ after inflation is
usually much lower than that for the thermal leptogenesis. Fortunately, in our
model, the saxion decay can help to relax the above cosmological constraint on
the gravitino mass and the reheating temperature $T_{R}$. After inflation, the
saxion field obtains an enormous potential for $s_{I}$, which results in an
initial displacement of the saxion field for saxion condensate oscillation.

In the early universe, the energy density of saxion $\rho_{s}$ begin to
dominate the universe at the temperature $T_{M}$ when $\rho_{s}=\rho_{r}$.
During the dominant epoch, the saxion started to decay and ended at the
temperature $T_{s}$. Due to saxion decay, a large amount of the entropy was
injected into the universe and thus the dark matter ($\tilde{G}$ and
$\tilde{a}$) relic density produced before can be diluted Banks  _et al._
(2003); Kawasaki and Nakayama (2008). The effect of dilution can be
parameterized by the factor $D_{s}$, which is given by

$D_{s}=\left(\frac{9g_{\ast}(T_{s})}{g_{\ast}(T_{M})}\right)^{1/4}\frac{m_{s}^{1/2}s_{I}^{2}}{\sqrt{\Gamma_{s}}M_{pl}^{2}},$
(20)

where $g_{\ast}(T_{s})=10.75$ and $g_{\ast}(T_{M})=228.75$ are the number of
relativistic degree of freedom at the temperature $T_{s}$ and $T_{M}$
respectively.

The relic density of the gravitino depends on whether it is in the thermal
bath Fujii and Yanagida (2002). The freeze-out temperature of the gravitino is
given by

$T^{f}_{3/2}\approx 10^{13}{\rm
GeV}\left(\frac{g_{\ast}}{230}\right)^{\frac{1}{2}}\left(\frac{m_{3/2}}{1{\rm
GeV}}\right)^{2}\left(\frac{1{\rm TeV}}{m_{\tilde{g}}}\right)^{2},$ (21)

where $m_{\tilde{g}}$ is the mass of gluino. $g_{\ast}$ is the effective
degrees of freedom of relativistic particles when the gravitinos are out of
equilibrium. If the reheating temperature $T_{R}>T^{f}_{3/2}$, the pre-
existing gravitinos are in the thermal equilibrium due to the high reheating
temperature $T_{R}$ and then freeze out as the expansion of universe. The
undiluted freeze-out gravitino yield at the equilibrium

$Y^{FO}_{3/2}=\frac{135\zeta\left(3\right)g_{3/2}}{8\pi^{4}g_{\ast}},$ (22)

where $g_{3/2}=4$ is the internal degrees of freedom of the gravitino, $\zeta$
is the zeta function. On the contrary, when $T_{R}<T_{3/2}^{f}$, the pre-
existing gravitinos may be out of the thermal equilibrium. In this case, the
gravitino can also be produced by the scattering processes, such as
$\tilde{g}\tilde{g}\to\tilde{G}\tilde{G}$, which is so-called freeze-in
mechanism. The resulting undiluted freeze-in gravitino yield can be obtained
by

$Y^{FI}_{3/2}=\sum^{3}_{i=1}y_{i}g^{2}_{i}\left(T_{R}\right)\left(1+\frac{M^{2}_{i}\left(T_{R}\right)}{3m^{2}_{3/2}}\right)\times\ln\left(\frac{k_{i}}{g_{i}\left(T_{R}\right)}\right)\left(\frac{T_{R}}{10^{10}{\rm
GeV}}\right)$ (23)

where $g_{i}$ and $M_{i}$ are the gauge coupling constants and the gaugino
masses for $M_{1}=M_{2}=500$GeV, $M_{3}=1$TeV. The constants $c_{i}$, $k_{i}$
and $y_{i}$ are associated with the Standard Model gauge group
$SU(3)_{c}\times SU(2)_{L}\times U(1)_{Y}$ Pradler and Steffen (2007a). Thus,
the total undiluted gravitino abundance depends on the relation between
reheating temperature $T_{R}$ and freeze-out temperature $T_{3/2}^{f}$, which
can be written as

$\Omega_{3/2}h^{2}=m_{3/2}s\left(T_{0}\right)h^{2}/\rho_{c}\left[Y^{FO}_{3/2}\Theta(T_{R}-T_{3/2}^{f})+Y^{FI}_{3/2}\Theta(T_{3/2}^{f}-T_{R}))\right].$
(24)

where $h$ is the Hubble constant in the units of 100km/s/Mpc.
$s\left(T_{0}\right)$ is the current entropy density and $\rho_{c}$ is the
critical density with
$\rho_{c}/\left[s\left(T_{0}\right)h^{2}\right]=3.6\times 10^{-9}{\rm GeV}$.
Similarly, if the axino is the LSP, the corresponding freeze-out temperature
of the axino $T_{\tilde{a}}^{f}$ is calculated by,

$T^{f}_{\tilde{a}}=10^{11}{\rm GeV}\left(\frac{\sqrt{2}V_{PQ}}{N_{\rm
DW}10^{12}{\rm GeV}}\right)^{2}\left(\frac{0.1}{\alpha_{s}}\right)^{3},$ (25)

where $N_{\rm DW}$ is the domain wall number and $\alpha_{s}$ is strong
interaction coupling constant. Provided that the reheating temperature $T_{R}$
is higher than the freeze-out temperature $T_{\tilde{a}}^{f}$, the effect of
$T_{R}$ will compensate the feeble interaction between the axino and other
particles and axino is in thermal equilibrium. Then the undiluted yield of the
axino in the thermal bath is given by

$Y^{FO}_{\tilde{a}}=\frac{135\zeta\left(3\right)g_{\tilde{a}}}{8\pi^{4}g_{\ast}},$
(26)

where the internal degrees of freedom $g_{\tilde{a}}=2$ is for the axino. On
the other hand, when the axino is out of the equilibrium, its undiluted
abundance arising from the scattering of gluon and gluino is estimated by

$Y^{FI}_{\tilde{a}}\simeq 2\times 10^{-6}{\rm
ln}\left(\frac{3}{g_{3}}\right)\left(\frac{N_{DW}}{6}\right)^{2}\left(\frac{T_{R}}{10^{10}{\rm
GeV}}\right)\left(\frac{10^{14}{\rm GeV}}{V_{PQ}}\right)^{2},$ (27)

where $g_{3}$ is the strong gauge coupling and we take $N_{\rm DW}=6$.
Therefore, the undiluted relic abundance of the axino can be written as

$\Omega_{\tilde{a}}h^{2}=m_{\tilde{a}}s\left(T_{0}\right)h^{2}/\rho_{c}\left[Y^{FI}_{\tilde{a}}\Theta(T_{\tilde{a}}^{f}-T_{R})+Y^{FO}_{\tilde{a}}\Theta(T_{R}-T_{\tilde{a}}^{f})\right].$
(28)

Figure 4: The constraint of DM relic density $0.075<\Omega_{\rm
DM}h^{2}<0.126$ (green band) for the gravitino LSP (upper panel) and the axino
LSP (lower panel) on the plane of the reheating temperature $T_{R}$ versus
their masses $m_{3/2}$ and $m_{\tilde{a}}$.
$f_{1}=\Omega_{3/2}h^{2}/\Omega_{\rm DM}h^{2}$ and
$f_{2}=\Omega_{\tilde{a}}h^{2}/\Omega_{\rm DM}h^{2}$ are the fraction of the
energy density of the gravitino and axino NLSP in the total DM energy density,
respectively. We assume the saxion mass $m_{s}=10^{3}$GeV. The right side of
the gray dash line corresponds to the lifetime of the NLSP $\tau_{\rm
NLSP}<t_{0}$, where $t_{0}=13.7$Gyr is the lifetime of universe. While the
region above the red line denotes the reheating temperature $T_{R}$ is higher
than freeze-out temperature of the gravitino $T_{3/2}^{f}$.

Assuming that the axino and gravitino were produced before $\rho_{s}$
dominated the universe, both will be diluted by the saxion decay. Moreover, we
require that the lifetime of the NLSP is longer than the age of universe
($t_{0}$). As a result, the total DM relic density is

$\Omega_{DM}h^{2}=\frac{1}{D_{s}}\left(\Omega_{3/2}h^{2}+\Omega_{\tilde{a}}h^{2}\right).$
(29)

In Fig 4, we show the constraint of DM relic density $0.075<\Omega_{\rm
DM}h^{2}<0.126$ for the gravitino LSP and the axino LSP, respectively. In the
region above the red line, the gravitino is in the thermal equilibrium and is
produced by the freeze-out mechanism. While in the region below the red line,
it never reaches the thermal equilibrium and is produced by the freeze-in
mechanism. It should be mentioned that the axino never reaches the thermal
equilibrium in our plots because the axino freeze-out temperature $T_{a}^{f}$
is less than the reheating temperature $T_{R}$. With the increase of the
reheating temperature $T_{R}$, the large dilution effect is required to avoid
overclosing the universe. In other words, the strong constraint on $T_{R}$
from the observed DM relic density can be relaxed by dilution. Since the high
$V_{PQ}$ can produce the large dilution factor, the reheating temperature
$T_{R}$ in the right panel is much larger than that in the left panel. In
addition, for a given mass ratio $R$, when the gravitino mass $m_{3/2}$
becomes small, the freeze-out temperature of the gravitino $T_{3/2}^{f}$ is
lower than that of the axino $T_{\tilde{a}}^{f}$. On the other hand, when
$m_{3/2}$ becomes heavy, $T_{3/2}^{f}$ will be larger than
$T_{\tilde{a}}^{f}$. Therefore, the gravitino will dominate the DM relic
density in the small $m_{3/2}$ region because of small freeze-out temperature.
While the axino will be the main component of the DM abundance in large
$m_{3/2}$ region. In the upper panel, we also find that the gravitino DM yield
$Y_{3/2}^{FO}$ will not change with the variation of the reheating temperature
$T_{R}$ when the gravitino is in thermal bath. This makes the fraction $f_{2}$
become larger with the increase of high $T_{R}$ because more axinos are
produced. While for the axino DM in the lower panel, it can be also seen that
there is a tension between the life time of the decaying gravitino and the
fraction of it in the total DM energy density that is required to be less then
10% by the mesurements of the CMB and matter power spectrum Poulin  _et al._
(2016). Therefore we will focus on the the gravitino DM plus the axino
decaying DM in the following study.

##  5 Axion from Late Decaying DM in Neutrino Experiments

Owing to the mass splitting between the NLSP and LSP, the NLSP can decay into
the LSP plus an axion in our model. As known, the differential events caused
by the absorption of the axion always depends on the energy differential flux
of the axion arriving at earth, $d\Phi_{a}/dE_{a}$, which include the
contributions in the galaxy and outside the galaxy. The energy spectrum of the
axion produced from the late decaying DM is given by,

$\frac{dN}{dE_{a}}=N_{a}\delta(E_{a}-E_{\rm em})$ (30)

where $E_{\rm em}$ is the energy of axion at emission and $N_{a}$ is the
number of axion in final state. Assuming that the parent particle $X$ is at
rest, we can calculate $E_{\rm em}$ by using four momentum conservation in the
massless limit $m_{a}=0$,

$E_{\rm em}=\frac{m_{X}}{2}(1-y^{2})$ (31)

where $y=R$ for the gravitino NLSP and $y=R^{-1}$ for the axino NLSP.

The differential flux of axion from DM decay within the galaxy arriving at
earth is obtained by the line of sight integral,

$\frac{d\Phi_{gal}}{dE_{a}}=\frac{fBr_{a}e^{-t_{0}/\tau_{X}}}{\tau_{X}m_{X}}\frac{dN}{dE_{a}}\times
R_{\rm sol}\rho_{\rm sol}\mathcal{J}$ (32)

where $R_{\rm sol}=8.33$ kpc is the distance to the galaxy center and
$\rho_{\rm sol}=0.3$ GeV${}^{3}/$cm is the DM energy density at the position
of Earth; $t_{0}=13.7$ Gyr and $\tau_{X}$ is the lifetime of universe and late
decaying DM respectively. $Br_{a}$ is the branching ratio of $X$ decay to
axions. Same as defined in Fig. 4, $f$ is the fraction of the decaying DM to
the total DM energy density , which has to be less than about 10% as the
constraints of CMB and matter power spectrum Poulin  _et al._ (2016). In
addition, assuming that $\mathcal{J}$ are independent of the angle and
direction, we average the $\mathcal{J}$ over all directions and take the value
of $\mathcal{J}=2.1$ derived from the NFW density profile Navarro  _et al._
(1996). Note that we take 5% Gaussian distribution to display the differential
galactic flux.

Figure 5: The dependence of the flux of supersymmetric axion dark radiation on
the axion kinetic energy in Milky Way galaxy and outside galaxy.

In the calculation of the differential flux from extra-galactic distances, one
should take the red-shift effect into account. Since the momentum of the axion
is inversely proportional to the scale factor $a$, the relation between the
momentum at emission $p_{\rm em}$ and that absorbed by detector $p_{a}$ is
related by the red-shift factor $z$,

$p_{\rm em}(z)=(1+z)\times p_{a}.$ (33)

Besides, the cosmic time at redshift $z$ for a flat universe is obtained

$t(z)=\frac{1}{3H_{0}\sqrt{\Omega_{\Lambda}}}{\rm
ln}\bigg{[}\frac{\sqrt{1+(\Omega_{m}/\Omega_{\Lambda})(1+z)^{3}}+1}{\sqrt{1+(\Omega_{m}/\Omega_{\Lambda})(1+z)^{3}}-1}\bigg{]},$
(34)

where $H_{0}$ is the current Hubble constant; $\Omega_{m}=0.315$ is the ratio
of matter energy density $\rho_{m}$ to the critical density
$\rho_{c}=3M_{pl}^{2}H_{0}^{2}$, while $\Omega_{\Lambda}=1-\Omega_{m}$ is the
ratio of cosmological constant $\rho_{\Lambda}$ to $\rho_{c}$. As a result,
the extra-galactic flux for a 2-body decay is given by

$\frac{d\Phi_{eg}}{dE_{a}}=\frac{N_{a}}{p_{a}}\frac{fBr_{a}\Omega_{\rm
dm}\rho_{c}}{\tau_{X}m_{X}H_{0}}\frac{e^{-t(z)/\tau_{X}}}{\sqrt{(1+z)^{3}\Omega_{m}+\Omega_{\Lambda}}}\Theta(z),$
(35)

where $\Omega_{\rm dm}=0.2607$ is the current DM density Aghanim  _et al._
(2020). The total differential flux $d\Phi_{a}/dE_{a}$ consists of the
galactic and the extra-galactic flux,

$\frac{d\Phi_{a}}{dE_{a}}=\left(\frac{d\Phi_{gal}}{dE_{a}}+\frac{d\Phi_{eg}}{dE_{a}}\right)$
(36)

From Eq. 32 and 35, it can be seen that the flux of the axion radiation
depends on the life-time $\tau_{NLSP}$ of the NLSP, the mass of NLSP $m_{X}$
and the fraction $f$. The light NLSP mass and large fraction can lead to a
large flux. However, for a given mass ratio $R$, a lighter NLSP mass, a longer
NLSP life time, which will suppress the flux greatly. Besides, if the NLSP is
the gravitino, the late decaying DM will dominate the DM relic density, which
is excluded by the CMB constraint. Hence, we only consider the axino NLSP as
late decaying DM.

It should be noted that the kinetic energy of the axion from the axino decay
depends on our model parameters. When the kinetic energy lies in the keV
range, the flux of our axion will be much smaller than that of the solar axion
Derbin  _et al._ (2011); Arisaka  _et al._ (2013) so that we cannot obtain the
sensitivity of our axion radiation in the solar axion experiments. Therefore,
we focus on the NLSP mass $m_{NLSP}$ is larger than keV scale, such as the MeV
axino, which will produce the MeV axion from the axino decay. On the other
hand, the MeV axion is weakly limited by cosmological and astrophysical
observations.

In Figure 5, we show the differential flux of the axion radiation. The
relevant parameters, such as the life-time $\tau_{\tilde{a}}$, the mass
$\tilde{a}$, the mass ratio $R$, the reheating temperature $T_{R}$ and the PQ
symmetry breaking scale $V_{PQ}$, are chosen to satisfy the aforementioned
constraints. As above discussed, the flux strongly relies on the fraction of
late-decaying DM, like in left-bottom panel, which reach about 10% and has the
much larger flux other three cases. Besides, we find that the galactic
contribution is lightly smaller than the extra-galactic contribution. Due to
the red-shift effect, the extra-galactic differential flux spreads and ends at
the corresponding the axion kinetic energy $E_{\rm em}$ in each panel.

Figure 6: The number of events of the photon from the inverse Primarkoff
process $a+A\to\gamma+A$ as function of the axion-photon coupling
$g_{a\gamma\gamma}$ in the Hyper-Kamiokande experiment. The four benchmark
points are the same as the Fig. 5

In order to detect such MeV axions, we consider the inverse Primakoff
scattering process $a+A\to\gamma+A$, in which the axion scatters off with the
atom $A$ induced by the axion-photon coupling. Since the inelastic scattering
cross section is suppressed by the atomic number $Z$, we study the elastic
scattering process. The corresponding elastic cross section is given by Abe
_et al._ (2021)

$\sigma_{\rm el}\left(E_{a}\right)=\frac{\alpha
Z^{2}g_{a\gamma\gamma}^{2}}{2}\bigg{[}{\rm
ln\left(2E_{a}r_{A}\right)-\frac{1}{2}}\bigg{]}$ (37)

where $\alpha$ is the fine-structure constant, $Z$ is the number of atom $A$,
$g_{a\gamma\gamma}$ is the axion-photon coupling and $r_{A}$ is the radius of
the atom. With this, we can estimate the sensitivity of our MeV axion from the
axino decay in the future neutrino experiments such as Hyper-Kamiokande. The
target atom for the Hyper-Kamiokande is the oxygen, whose atomic number is
$Z=8$. There are $N_{O}\simeq 6.26\times 10^{33}$ oxygen atoms for a fiducial
mass of 187 kt water. By using the relation of the axion production rate
$\dot{N_{a}}$ and the differential axion flux $d\Phi_{a}/dE_{a}$ Payez  _et
al._ (2015), we can obtain the expected number of events as

$N_{\rm event}=\int\frac{d\Phi_{a}}{dE_{a}}\sigma_{\rm
el}\left(E_{a}\right)dE_{a}\times\Delta t\times N_{O}$ (38)

where the $\Delta t$ is the exposure time of Hyper-Kamiokande experiment and
$r_{A}=2\times 10^{-10}$m is the radius of oxygen atom. In Figure 6, we show
the number of events of the inverse Primakoff scattering of our axion with the
atoms, where we take the exposure time $\Delta t=1$ year as an example. The
four benchmark points are the same as the Fig. 5. Assuming a possion
distribution, we can see that axion-photon coupling $g_{a\gamma\gamma}$ for
$m_{\tilde{a}}=30$ MeV and $R=10$ can be excluded to $4\times 10^{-9}$ GeV-1
by the Hyper-Kamiokande experiment. The bounds for other cases are weaker
because of the small flux.

##  6 Conclusion

In this paper, We have investigated the dark sector of supersymmetric DFSZ
model, where the gravitino/axino is the late decaying DM. Owing to the saxion
decay, the gravitino/axino problem that its relic density easily overclose the
universe can be elegantly solved by the dilution effect. On the other hand,
the early axion produced by the saxion decay as dark radiation contributes the
extra effective number of neutrino species $\Delta N_{\rm eff}$. However, due
to the cosmological constraints on $\Delta N_{\rm eff}$ by CMB observations,
increasing simply $\Delta N_{\rm eff}$ only alleviate the Hubble tension to
$3\sigma$. On the other, we find that a MeV axion is emitted from the late
decaying DM decay in the late universe can be tested in the neutrino
experiment like Hyper-Kamiokande via the inverse Primakoff scattering. With an
exposure time $\Delta t=1$ year of Hyper-Kamiokande, the axion-photon coupling
$g_{a\gamma\gamma}$ can be excluded to $4\times 10^{-9}$ GeV-1.

##  7 Acknowledgments

We thank Shan-Shan Weng for their helpful suggestions and discussions. This
work is supported by the National Natural Science Foundation of China (NNSFC)
under grant No. 11805161 and 11705093, by Jiangsu Specially Appointed
Professor Program, and Natural Science Foundation of Shandong Province under
the grants ZR2018QA007. BZ is also supported by the Basic Science Research
Program through the National Research Foundation of Korea (NRF) funded by the
Ministry of Education, Science and Technology (NRF-2019R1A2C2003738), and by
the Korea Research Fellowship Program through the NRF funded by the Ministry
of Science and ICT (2019H1D3A1A01070937).

## References

  * Peccei and Quinn (1977) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977).
  * Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ (1980) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 166, 493 (1980).
  * Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980).
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Kim and Nilles (1984) J. E. Kim and H. P. Nilles, Phys. Lett. B 138, 150 (1984).
  * Kim and Carosi (2010) J. E. Kim and G. Carosi, Rev. Mod. Phys. 82, 557 (2010), [Erratum: Rev.Mod.Phys. 91, 049902 (2019)], arXiv:0807.3125 [hep-ph] .
  * Cao  _et al._ (2014) J. Cao, C. Han, L. Wu, P. Wu, and J. M. Yang, JHEP 05, 056 (2014), arXiv:1311.0678 [hep-ph] .
  * Duan  _et al._ (2018a) G. H. Duan, W. Wang, L. Wu, J. M. Yang, and J. Zhao, Phys. Lett. B 778, 296 (2018a), arXiv:1711.03893 [hep-ph] .
  * Abdughani  _et al._ (2018) M. Abdughani, L. Wu, and J. M. Yang, Eur. Phys. J. C 78, 4 (2018), arXiv:1705.09164 [hep-ph] .
  * Duan  _et al._ (2018b) G. H. Duan, K.-I. Hikasa, J. Ren, L. Wu, and J. M. Yang, Phys. Rev. D 98, 015010 (2018b), arXiv:1804.05238 [hep-ph] .
  * Abdughani and Wu (2020) M. Abdughani and L. Wu, Eur. Phys. J. C 80, 233 (2020), arXiv:1908.11350 [hep-ph] .
  * Abdughani  _et al._ (2019) M. Abdughani, K.-I. Hikasa, L. Wu, J. M. Yang, and J. Zhao, JHEP 11, 095 (2019), arXiv:1909.07792 [hep-ph] .
  * Abdughani  _et al._ (2021) M. Abdughani, Y.-Z. Fan, L. Feng, Y.-L. Sming Tsai, L. Wu, and Q. Yuan, (2021), arXiv:2104.03274 [hep-ph] .
  * Covi  _et al._ (1999) L. Covi, J. E. Kim, and L. Roszkowski, Phys. Rev. Lett. 82, 4180 (1999), arXiv:hep-ph/9905212 .
  * Covi  _et al._ (2001) L. Covi, H.-B. Kim, J. E. Kim, and L. Roszkowski, JHEP 05, 033 (2001), arXiv:hep-ph/0101009 .
  * Covi and Kim (2009) L. Covi and J. E. Kim, New J. Phys. 11, 105003 (2009), arXiv:0902.0769 [astro-ph.CO] .
  * Bae  _et al._ (2015) K. J. Bae, H. Baer, E. J. Chun, and C. S. Shin, Phys. Rev. D 91, 075011 (2015), arXiv:1410.3857 [hep-ph] .
  * Co  _et al._ (2017a) R. T. Co, F. D’Eramo, and L. J. Hall, JHEP 03, 005 (2017a), arXiv:1611.05028 [hep-ph] .
  * Hamaguchi  _et al._ (2017) K. Hamaguchi, K. Nakayama, and Y. Tang, Phys. Lett. B 772, 415 (2017), arXiv:1705.04521 [hep-ph] .
  * Hooper  _et al._ (2012) D. Hooper, F. S. Queiroz, and N. Y. Gnedin, Phys. Rev. D 85, 063513 (2012), arXiv:1111.6599 [astro-ph.CO] .
  * Gu  _et al._ (2020) Y. Gu, M. Khlopov, L. Wu, J. M. Yang, and B. Zhu, Phys. Rev. D 102, 115005 (2020), arXiv:2006.09906 [hep-ph] .
  * Anchordoqui (2021) L. A. Anchordoqui, Phys. Rev. D 103, 035025 (2021), arXiv:2010.09715 [hep-ph] .
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A. G. Riess, Nature Astron. 3, 891 (2019), arXiv:1907.10625 [astro-ph.CO] .
  * Riess  _et al._ (2019) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Astrophys. J. 876, 85 (2019), arXiv:1903.07603 [astro-ph.CO] .
  * Di Valentino  _et al._ (2021) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, (2021), arXiv:2103.01183 [astro-ph.CO] .
  * Binder  _et al._ (2018) T. Binder, M. Gustafsson, A. Kamada, S. M. R. Sandner, and M. Wiesner, Phys. Rev. D 97, 123004 (2018), arXiv:1712.01246 [astro-ph.CO] .
  * Pandey  _et al._ (2020) K. L. Pandey, T. Karwal, and S. Das, JCAP 07, 026 (2020), arXiv:1902.10636 [astro-ph.CO] .
  * Vattis  _et al._ (2019) K. Vattis, S. M. Koushiappas, and A. Loeb, Phys. Rev. D 99, 121302 (2019), arXiv:1903.06220 [astro-ph.CO] .
  * Xiao  _et al._ (2020) L. Xiao, L. Zhang, R. An, C. Feng, and B. Wang, JCAP 01, 045 (2020), arXiv:1908.02668 [astro-ph.CO] .
  * Blinov  _et al._ (2020) N. Blinov, C. Keith, and D. Hooper, JCAP 06, 005 (2020), arXiv:2004.06114 [astro-ph.CO] .
  * Pagels and Primack (1982) H. Pagels and J. R. Primack, Phys. Rev. Lett. 48, 223 (1982).
  * Moroi  _et al._ (1993) T. Moroi, H. Murayama, and M. Yamaguchi, Phys. Lett. B 303, 289 (1993).
  * Asaka  _et al._ (2000) T. Asaka, K. Hamaguchi, and K. Suzuki, Phys. Lett. B 490, 136 (2000), arXiv:hep-ph/0005136 .
  * Bolz  _et al._ (2001) M. Bolz, A. Brandenburg, and W. Buchmuller, Nucl. Phys. B 606, 518 (2001), [Erratum: Nucl.Phys.B 790, 336–337 (2008)], arXiv:hep-ph/0012052 .
  * Roszkowski  _et al._ (2005) L. Roszkowski, R. Ruiz de Austri, and K.-Y. Choi, JHEP 08, 080 (2005), arXiv:hep-ph/0408227 .
  * Pradler and Steffen (2007a) J. Pradler and F. D. Steffen, Phys. Lett. B 648, 224 (2007a), arXiv:hep-ph/0612291 .
  * Pradler and Steffen (2007b) J. Pradler and F. D. Steffen, Phys. Rev. D 75, 023509 (2007b), arXiv:hep-ph/0608344 .
  * Cheung  _et al._ (2012) C. Cheung, G. Elor, and L. J. Hall, Phys. Rev. D 85, 015008 (2012), arXiv:1104.0692 [hep-ph] .
  * Co  _et al._ (2017b) R. T. Co, F. D’Eramo, L. J. Hall, and K. Harigaya, JHEP 07, 125 (2017b), arXiv:1703.09796 [hep-ph] .
  * Kawasaki  _et al._ (2005) M. Kawasaki, K. Kohri, and T. Moroi, Phys. Rev. D 71, 083502 (2005), arXiv:astro-ph/0408426 .
  * Hochberg  _et al._ (2017) Y. Hochberg, T. Lin, and K. M. Zurek, Phys. Rev. D 95, 023013 (2017), arXiv:1608.01994 [hep-ph] .
  * Cui  _et al._ (2018) Y. Cui, M. Pospelov, and J. Pradler, Phys. Rev. D 97, 103004 (2018), arXiv:1711.04531 [hep-ph] .
  * Kuo  _et al._ (2021) J.-L. Kuo, M. Pospelov, and J. Pradler, (2021), arXiv:2102.08409 [hep-ph] .
  * Avignone  _et al._ (1998) F. T. Avignone, III _et al._ (SOLAX), Phys. Rev. Lett. 81, 5068 (1998), arXiv:astro-ph/9708008 .
  * Bernabei  _et al._ (2001) R. Bernabei _et al._ , Phys. Lett. B 515, 6 (2001).
  * Dent  _et al._ (2020) J. B. Dent, B. Dutta, J. L. Newstead, and A. Thompson, Phys. Rev. Lett. 125, 131805 (2020), arXiv:2006.15118 [hep-ph] .
  * Abe  _et al._ (2018) K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 [physics.ins-det] .
  * Acciarri  _et al._ (2015) R. Acciarri _et al._ (DUNE), (2015), arXiv:1512.06148 [physics.ins-det] .
  * Abi  _et al._ (2018) B. Abi _et al._ (DUNE), (2018), arXiv:1807.10334 [physics.ins-det] .
  * Abi  _et al._ (2020) B. Abi _et al._ (DUNE), (2020), arXiv:2008.12769 [hep-ex] .
  * An  _et al._ (2016) F. An _et al._ (JUNO), J. Phys. G 43, 030401 (2016), arXiv:1507.05613 [physics.ins-det] .
  * Co  _et al._ (2016) R. T. Co, F. D’Eramo, and L. J. Hall, Phys. Rev. D 94, 075001 (2016), arXiv:1603.04439 [hep-ph] .
  * Vagnozzi (2020) S. Vagnozzi, Phys. Rev. D 102, 023518 (2020), arXiv:1907.07569 [astro-ph.CO] .
  * Blinov  _et al._ (2019) N. Blinov, K. J. Kelly, G. Z. Krnjaic, and S. D. McDermott, Phys. Rev. Lett. 123, 191102 (2019), arXiv:1905.02727 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Abazajian  _et al._ (2016) K. N. Abazajian _et al._ (CMB-S4), (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Weinberg (1982) S. Weinberg, Phys. Rev. Lett. 48, 1303 (1982).
  * Khlopov and Linde (1984) M. Y. Khlopov and A. D. Linde, Phys. Lett. B 138, 265 (1984).
  * Dine  _et al._ (1995) M. Dine, A. E. Nelson, and Y. Shirman, Phys. Rev. D 51, 1362 (1995), arXiv:hep-ph/9408384 .
  * Giudice and Rattazzi (1999) G. F. Giudice and R. Rattazzi, Phys. Rept. 322, 419 (1999), arXiv:hep-ph/9801271 .
  * Randall and Sundrum (1999) L. Randall and R. Sundrum, Nucl. Phys. B 557, 79 (1999), arXiv:hep-th/9810155 .
  * Ellis  _et al._ (2004) J. R. Ellis, K. A. Olive, Y. Santoso, and V. C. Spanos, Phys. Lett. B 588, 7 (2004), arXiv:hep-ph/0312262 .
  * Buchmuller  _et al._ (2006) W. Buchmuller, K. Hamaguchi, and J. Kersten, Phys. Lett. B 632, 366 (2006), arXiv:hep-ph/0506105 .
  * Dudas  _et al._ (2017) E. Dudas, Y. Mambrini, and K. Olive, Phys. Rev. Lett. 119, 051801 (2017), arXiv:1704.03008 [hep-ph] .
  * Goto and Yamaguchi (1992) T. Goto and M. Yamaguchi, Phys. Lett. B 276, 103 (1992).
  * Chun  _et al._ (1992) E. J. Chun, J. E. Kim, and H. P. Nilles, Phys. Lett. B 287, 123 (1992), arXiv:hep-ph/9205229 .
  * Chun and Lukas (1995) E. J. Chun and A. Lukas, Phys. Lett. B 357, 43 (1995), arXiv:hep-ph/9503233 .
  * Kim and Seo (2012) J. E. Kim and M.-S. Seo, Nucl. Phys. B 864, 296 (2012), arXiv:1204.5495 [hep-ph] .
  * Banks  _et al._ (2003) T. Banks, M. Dine, and M. Graesser, Phys. Rev. D 68, 075011 (2003), arXiv:hep-ph/0210256 .
  * Kawasaki and Nakayama (2008) M. Kawasaki and K. Nakayama, Phys. Rev. D 77, 123524 (2008), arXiv:0802.2487 [hep-ph] .
  * Fujii and Yanagida (2002) M. Fujii and T. Yanagida, Phys. Lett. B 549, 273 (2002), arXiv:hep-ph/0208191 .
  * Poulin  _et al._ (2016) V. Poulin, P. D. Serpico, and J. Lesgourgues, JCAP 08, 036 (2016), arXiv:1606.02073 [astro-ph.CO] .
  * Navarro  _et al._ (1996) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996), arXiv:astro-ph/9508025 .
  * Derbin  _et al._ (2011) A. V. Derbin, A. S. Kayunov, V. N. Muratova, D. A. Semenov, and E. V. Unzhakov, Phys. Rev. D 83, 023505 (2011), arXiv:1101.2290 [hep-ex] .
  * Arisaka  _et al._ (2013) K. Arisaka, P. Beltrame, C. Ghag, J. Kaidi, K. Lung, A. Lyashenko, R. D. Peccei, P. Smith, and K. Ye, Astropart. Phys. 44, 59 (2013), arXiv:1209.3810 [astro-ph.CO] .
  * Abe  _et al._ (2021) T. Abe, K. Hamaguchi, and N. Nagata, Phys. Lett. B 815, 136174 (2021), arXiv:2012.02508 [hep-ph] .
  * Payez  _et al._ (2015) A. Payez, C. Evoli, T. Fischer, M. Giannotti, A. Mirizzi, and A. Ringwald, JCAP 02, 006 (2015), arXiv:1410.3747 [astro-ph.HE] .