#### The distribution of inelastic dark matter in the Sun

Eur. Phys. J. C
The distribution of inelastic dark matter in the Sun
Mattias Blennow 1 2
Stefan Clementz 2
Juan Herrero-Garcia 0
0 ARC Center of Excellence for Particle Physics at the Terascale (CoEPP), University of Adelaide , Adelaide, SA 5005 , Australia
1 Instituto de Física Teórica UAM/CSIC , Calle Nicolás Cabrera 13-15, Cantoblanco 28049, Madrid , Spain
2 Department of Physics, School of Engineering Sciences, KTH Royal Institute of Technology, AlbaNova University Center , 106 91 Stockholm , Sweden
If dark matter is composed of new particles, these may become captured after scattering with nuclei in the Sun, thermalize through additional scattering, and finally annihilate into neutrinos that can be detected on Earth. If dark matter scatters inelastically into a slightly heavier (O(10 − 100) keV) state it is unclear whether thermalization occurs. One issue is that up-scattering from the lower mass state may be kinematically forbidden, at which point the thermalization process effectively stops. A larger evaporation rate is also expected due to down-scattering. In this work, we perform a numerical simulation of the capture and thermalization process in order to study the evolution of the dark matter distribution. We then calculate and compare the annihilation rate with that of the often assumed Maxwell-Boltzmann distribution. We also check if equilibrium between capture and annihilation is reached. We find that, unless the mass splitting is very small ( 50 keV) and/or the dark matter has a sub-dominant elastic cross section, the dark matter distribution does not reach a stationary state, it is not isothermal, annihilation is severely suppressed, and equilibrium is generally not reached. We also find that evaporation induced by down-scattering is not effective in reducing the total dark matter abundance.
1 Introduction
A popular class of models for explaining a number of
observations in astrophysical systems is that of particle dark matter
(DM) [1–4]. In the event that DM interacts with the
particles of the standard model (SM), many different methods
for observing it have been proposed. Over the last decades
a number of experiments have been able to place impressive
bounds on the DM mass and its interaction cross sections.
One of the many ways of searching for DM is to look for
the effects that it may have on the Sun as it is captured by
scattering against solar material [5, 6]. If DM annihilates (e.g.,
thermal relics), SM particles can be produced, which in turn
decay or otherwise interact to give rise to a flux of high energy
neutrinos or, in more exotic scenarios, to other SM particles.
The neutrinos produced can be searched for in neutrino
telecopes on Earth [7–12], with various collaborations having
performed such searches with no positive detection [13–16].
The accumulation of large amounts of DM in the Sun may
also affect helioseismology and the solar temperature. These
modifications can potentially lead to observational effects
with the possibility to constrain DM properties or alleviate
the solar composition problem [17–32].
A collection of DM models that could possibly be probed
by the production of neutrinos from DM annihilations inside
the Sun is that of inelastic DM [33]. These models were
originally introduced to reconcile the annual modulation
observation of the DAMA/LIBRA experiment [34] with the
null results of the CDMS experiment [35], which ruled out
its explanation in terms of a standard elastically scattering
DM particle. Now DAMA is also incompatible with a
number of other experiments, including the large Xenon-based
experiments LUX [36], PandaX [37], and XENON1T [38].
It should be mentioned that it is also difficult to reconcile
DAMA with other direct detection experiments in the case
of inelastic scattering [39]. These models are also tightly
constrained by direct detection (DD) experiments [40–43].
The evasion of the bounds from CDMS in inelastic DM
models came from the introduction of a small mass splitting
δ that separated two different DM states which, upon
scattering, change from one to the other. In the scattering
process, some energy is converted from kinetic energy to mass,
which gives this type of model its name. The introduced mass
splitting has a large impact on the scattering kinematics of
DM and translates into altered solar capture rates of inelastic
DM. Capture of inelastic DM in the Sun has been discussed
for both the cases of endothermic and exothermic
scattering in, e.g., Refs. [44–49]. See also Refs. [50–52] for
studies of inelastic DM capture by compact stars such as white
dwarves and neutron stars. Inelastic DM has also been
proposed as a solution to the small scale structure problems in
models where a light mediator induces large self-scattering
cross sections [53–55], and in models where the more
massive state is unstable [56–58].
Upon being captured inside the Sun, DM is often assumed
to “instantaneously” thermalize with the surrounding plasma,
in which case its number density distribution is well
described by a Boltzmann distribution with a specific
temperature, i.e., f ∼ exp(−E / T ), at all times, see e.g. Refs. [12,
59]. In Ref. [60] it has also been shown using numerical
simulations that the Boltzmann distribution is a reasonable
assumption and that the thermalized DM is generally
concentrated in the core of the Sun. Therefore significant
annihilation can occur. However, in the case of inelastic DM with
sizable δ, scattering of a particle in the lower mass state is
only kinetically allowed if a large amount of kinetic energy
is supplied to the collision. This implies that the DM
particles scatter only off high velocity (and thus Boltzmann
suppressed) solar nuclei, so that the scattering probability is very
small, or that their orbit takes them from a radius with a large
gravitational potential into regions closer to the solar center,
which provides the necessary kinetic energy. Moreover, when
a DM particle subsequently scatters from the higher to the
lower mass state it may be boosted by the significant amount
of energy that is released to a velocity such that it is no longer
gravitationally bound. With this in mind, it is not obvious
that a thermalized distribution is to be expected in the case
of inelastic DM, nor is it clear if and when evaporation has
to be taken into account. Both of these phenomena can have
an impact on the annihilation rate of captured DM. A
nonthermal distribution could alter the annihilation rate so that a
larger population of DM must be present in order for
equilibrium between solar capture and annihilation to occur, while
evaporation would reduce the total number of particles that
can annihilate. The assumption of capture-annihilation
equilibrium allows one to bypass the annihilation rate in favour of
a direct link between neutrino rates and the solar capture rate.
Due to the effects of inelastic scattering on the annihilation
rate it is unclear if this is a justified assumption.
In this paper, we study the thermalization process of
inelastic DM using numerical simulations. We analyse the
impact that inelastic scattering has on the annihilation and
evaporation rates. In Sect. 2 we discuss the inelastic DM
framework and the relevant kinematic effects introduced by
a mass splitting. In Sect. 3, we describe our approach to the
problem and the numerical implementation of the simulation.
The results are presented and discussed in Sect. 4. Finally,
we summarize and give our conclusions in Sect. 5.
2 Inelastic dark matter
There are various scenarios in which inelastic DM appears
naturally [61–63]. In the simplest models inelastic DM
consists of two states χ and χ ∗, with masses mχ < mχ ∗ that
satisfy
(
1
)
mχ ∗ − mχ = δ
mχ .
Although the mass splitting is very small relative to the DM
masses, it has a significant impact on the resulting differential
rates in DD experiments, as was noted in Ref. [33].
According to our definition of the DM masses, χ ∗ is the slightly
more massive state, which indicates that the scattering
process is endothermic when the incoming DM particle is a χ ,
and exothermic when the incoming particle is a χ ∗. Below,
we discuss the endothermic case. The case of exothermic
scattering can be recovered by substituting δ → −δ. We
define up-scattering as the process in which a particle in the
lower mass state scatters with the solar nuclei labelled by A
into the higher mass state, i.e., χ A → χ ∗ A. Down-scatter
refers to the opposite reaction, χ ∗ A → χ A.
We primarily study inelastic DM that couples to protons
and neutrons through spin-independent interactions. We
disregard spin-dependent scattering as in the Sun it is mainly
hydrogen that carries spin. The reason is that capture of the
χ state requires δ to be so small that scattering is essentially
elastic for DM masses heavier than a few times that of the
proton, while χ ∗ can still be captured in significant
numbers [28]. However, for the latter, upon being captured and
as long as δ is non-negligible, these particles that are now
in the χ state would find themselves unable to up-scatter.
Over time, this would create a cloud of loosely bound DM
particles that is far too diffuse for annihilation to take place
efficiently. In other words, capture via spin-dependent
interactions of inelastic dark matter is only interesting in the limit
where δ is so tiny that the model is essentially elastic and
thermalization is expected.
In the following we focus on the case where the galactic
halo is composed of both χ and χ ∗ with equal abundances,
in which case both species are captured by the Sun. This is
plausibly the case, as the temperature at which the overall DM
abundance freeze-out occurs is of the order mχ /20, which
far exceeds the values of δ considered in this work.
2.1 Scattering kinematics
The scattering kinematics of DM is extensively covered in
Ref. [48]. Here we briefly review the scattering kinematics
that are important for the discussions that follow.
urel,lower =
2δ
μχ A
,
or there is simply not enough kinetic energy to produce
the heavier state. On the other hand, there is no such
constraint for exothermic scattering, which is always
kinematically allowed.
When a collision occurs, the solution to the energy and
momentum conservation equations yields the largest recoil
energy of the target solar nucleus, Emax, and the smallest one,
(
2
)
(
3
)
(
4
)
(
5
)
Emin:
Emax (min) = 2
2
μχ A
mχ m A
− μmχAA δ ,
Ei,kin
1 (+−)
mχ
1 − μχ A Ei,kin
δ
where m A is the mass of the target nucleus and Ei,kin =
Ei − φ (r ) is the kinetic energy of the incoming DM particle.
The allowed range of recoil energies E R for a given Ei,kin is
Emin ≤ E R ≤ Emax.
In the rest frame of the target nucleus, the angle θ between
the velocity vectors of the incoming and outgoing DM
particle satisfies
cos(θ ) =
Ei,kin + Ef,kin − m A ER/mχ ,
2 Ei,kin Ef,kin
where Ef,kin = Ef −φ (r ) is the kinetic energy of the outgoing
DM particle. This relation is needed to determine the change
of trajectory of a DM particle as it scatters.
When DM scatters inelastically, the only modification to
the cross section with respect to the elastic case is a
multiplying phase-space factor. If the elastic DM–nucleus scattering
cross section is σ0, then the inelastic scattering cross section
σinel would be
σinel =
1 −
where vrel is the relative speed between the DM particle and
its target and μχ A is the DM–target reduced mass. When
DM scatters endothermically, the relative velocity between
the DM particle and its target must exceed
f˙α =
αβ fβ + Cα − fα
αβ fβ .
β
β
3 Solar capture, thermalization, and annihilation
The method we employ to investigate the thermalization of
inelastic DM is similar to that of Refs. [64, 65]. This section
presents the implementation of the numerical simulation.
1 We assume spherical symmetry, i.e., orbits in different planes are
equivalent. We therefore only use L, the total magnitude of L, and not
its 3 components.
2 Note that in the literature the last term is often written as 2 ann N 2.
In this work we absorb the factor of 2 into ann and use a differently
normalized number density distribution function.
3.1 The phase space evolution of a captured population
The effective classical Hamiltonian H describing a particle
moving in a central potential φ (r ) is given by
H = mχ E = 21 mχ r˙2 +
≡ 21 mχ r˙2 + mχ Veff (L , r ).
It is very convenient, as we have done, to define the reduced
energy E as the energy divided by the DM mass, as well as
the reduced angular momentum L = r ×v, given the fact that
the orbit of a particle in a central potential is independent of
its mass. Since E and L = |L | are conserved quantities, it is
convenient to describe the DM particle orbit by these
quantities rather than position and velocity.1 For a given angular
momentum L , the smallest energy Emin(L ) of a particle with
a trajectory that intersects the Sun is given by
min Veff (L , r ) .
Emin(L ) = r ≤R
Taking the above into consideration, we find it convenient to
define the combination α = (E , L ) as a label for the phase
space position of a particular DM particle.
The time evolution of the total number of captured DM
particles N (t ) follows the differential equation
N˙ = C
− Eevap −
ann ,
where N˙ ≡ d N /dt , C is the solar capture rate, Eevap ∝ N is
the evaporation rate, and ann ∝ N 2 is the rate at which
particles are annihilated.2 If evaporation is neglected, the
equilibrium solution (N˙ = 0) of the evolution equation reads
ann ,
C
=
which implies that there is equilibrium between capture and
annihilations. In order to rigorously test this condition one
must calculate the annihilation rate. However, this requires
knowledge of how the DM is distributed in the Sun. In our
simulations, the distribution of particles is discretized in E
and L such that fα describes the number of particles in a
particular state α. The evolution of the distribution is then
governed by the equation
(
6
)
(
7
)
(
8
)
(
9
)
(
10
)
Each element in f contains the total number of particles in
state α, while each element in C gives the capture rate into the
corresponding state. The off-diagonal elements in αβ (α =
β) give the rate with which particles in state β scatter against
solar nuclei and end up in the state α. The diagonal entries
in are negative and correspond to the rate at which
particles scatter from the corresponding state to all other states,
including evaporation, i.e., positive energy states where the
DM particle escapes the Sun’s gravitational well. Finally,
αβ gives the rate at which a particle in state α annihilates
with a particle in state β.3
We also need to know the fraction of the time that a particle
in a given state α finds itself at a radius r as it travels between
the maximal radius, r+, and the minimal radius, r−, of the
complete orbit. These are found by solving Eq. (
6
) with the
substitution r˙ = 0. The time it takes the particle to move
between radius r1 to radius r2 can be found by isolating r˙ in
Eq. (
6
), which leads to
T (r1, r2) =
dt =
r2 dr
r1
r =
˙
r2
r1
dr
√2(E − Veff (L , r ))
(
11
)
.
Integrating from r1 = r− to r2 = r+ gives the time that a
particle needs to complete half an orbit.
3.2 Solar capture
The derivation of the solar capture rate of DM particles
originating from the DM halo dates back several decades [5, 6]. In
the standard calculation any information on the DM energy
and angular momentum post scattering is discarded, since
any particle with E < 0 is counted towards the total capture.
However, here we are not only interested in the total capture
rate but also in the distribution of the captured particles in
E -L space. We therefore give a short description of how we
compute Cα .
The solar capture rate of DM in differential form is given
by [6]
dC = π nχ
1 − (L /r w)2
f (u) dσ
u
d E R
2
(w, E R ) n A(r )
dr du d E R d L 2 .
(
12
)
In the above, we have assumed that the target nuclei are
stationary. Here r is the radius at which scattering occurs, n A(r )
is the local density of target particles of species A, u is the
DM velocity at a distance where the gravitational potential of
the Sun is negligible, and w = u2 + vesc(r )2 is the
velocity of the particle at radius r . The local halo number density
3 In the case of DM self-capture due to self-interactions, see e.g.,
Ref. [59], αβ would also incorporate that effect.
of DM and its speed distribution at the location of the Sun
enter explicitly through nχ and f (u). Note that the velocity
distribution f˜(u) is normalized such that
(
13
)
(
14
)
(
15
)
f˜(u) d3u =
f (u) du = 1,
where
f (u) =
f˜(u, θ , φ) u2 sin(θ ) dθ dφ .
The quantity L 2 is the square of the reduced angular
momentum of the incoming DM particle from the DM halo, and
dσ/d E R (w, E R ) is the differential cross section [33]
dσ
d E R
(w, E R ) =
m2μAA2χp2wσχ2p |F (E R )|2 .
Here it has been assumed that the coupling between DM
and nuclei is isospin conserving, leading to the A2
enhancement of the cross section, where A is the total number of
nucleons.4 Furthermore, σχ p is the DM–proton cross section
entering as σ0 in Eq. (
2
), μχ p is the DM–proton reduced
mass, and w is the relative velocity between the DM and
the nucleus. Interestingly, the phase-space factor relating the
inelastic and elastic scattering cross sections in Eq. (
2
) is
cancelled. The form factor F (E R ) accounts for the
decoherence in the DM–nucleus scattering process when the
momentum transfer q(E R ) is large. The latter is related to E R by
q = √2m A E R .
In order to calculate Cα , the integrand of Eq. (
12
) is
discretized over the region of relevant r , u, E R and L . The
integration range in r is 0 < r < R , where R is the solar
radius, and that for L 2 is 0 < L 2 < r 2w2. The limits in
u and E R are complicated and we refer to the discussion in
Ref. [48]. To ensure that the mesh is fine enough we calculate
the total capture rate as C = α Cα as well as by
integrating Eq. (
12
) as done in Ref. [48] (but using the Helm form
factor given in eq. (
42
) rather than the very frequently used
exponential form factor). We have verified that the two agree
to better than 1 % accuracy.
At each discretization point, the incoming DM velocity
vector wi can be reconstructed. When the DM particle
scatters, its energy post-collision is known from
Ef = Ei − E R − δ/mχ ,
(
16
)
where Ei (f) is the energy of the incoming (outgoing) DM
particle. The outgoing DM speed wf is known from the
equa4 For isospin violating DM, A2 is instead replaced by the factor (Z +
( A − Z ) fn / f p)2, where Z is the number of protons and f p ( fn ) is the
coupling of DM to protons (neutrons). Note that f p can be absorbed in
the definition of σχ p.
αβ =
Rβ→α (ri ).
i
αα = −
Rα (ri ),
i
We use Monte-Carlo methods to find the probability
distribution for a scattering DM particle at each discretization point
to end up in a state α. Finally, C is found by summing over
all discretized states weighted by their probability densities.
which also includes evaporation.
To find Pβ→α (r ), the DM velocity vector wi of a particle
in state β = (Ei, L i) before scattering is required. In the
Sun’s rest frame, its magnitude is found to be
in the matrix are given by the sum of contributions from
all shells that the particle passes through on its orbit
The diagonal elements of the matrix are given by the
negative of the total scattering rate from state α,
tion above. The angle θ between the incoming and outgoing
velocity vector is given by Eq. (
5
). There is also an azimuthal
angle ϕ around wi at which the outgoing velocity vector lies
that is randomly distributed in the interval 0 to 2π . In terms
of these angles, the angular momentum of the outgoing DM
particle L f is given by
⎡
3.3 Scattering among different states
When DM particles have been captured, occasional
scattering with solar nuclei takes a particle initially in the β = (Ei,
L i) state into the state α = (Ef , L f ). The differential
scattering rate at radius r of a DM particle with velocity w(r ),
travelling through a gas of nuclei of element A with
velocity v, number density n A(r ), and velocity distribution of the
nuclei f A(r, v), is given by [64]
d R(r ) = σ n A(r ) f A(r, v) |w(r ) − v| d3v.
The velocities of nuclei in the Sun follow the Boltzmann
distribution
f A(r, v) =
m A
2π T (r )
3/2
exp
m Av2
− 2T (r )
,
where T (r ) is the temperature of the solar plasma at radius
r . The cross section σ that enters is the integral over the
differential cross section in the frame in which the nucleon
is stationary.
In order to find αβ we discretize the Sun into thin
spherical shells with radii ri . Under the assumption that DM
particles complete many orbits between interactions, the rate at
which particles in state β scatter at radius ri and end up in
state α is given by
Rβ→α (ri ) = Rβ (ri ) Tβ (ri ) Pβ→α (ri ).
Here, Rβ (ri ) is the total scattering rate at radius ri and Tβ (ri )
is the fraction of the orbital time that the particle spends
inside the shell. The factor Pβ→α (ri ) is the probability that
the particle ended up in the particular state α after it scattered
at ri . Having calculated the above, the off-diagonal elements
(
18
)
(
19
)
(
20
)
wi(r ) =
2Ei − 2φ (r ).
We are free to choose a coordinate system in which r =
(r, 0, 0). The angle between r and wi(r ) is given by ξ =
sin−1(L i/r wi(r )). Thus the DM velocity vector can be
written as
wi(r ) = wi(r ) (cξ , sξ , 0),
where we use sx = sin(x ) and cx = cos(x ). The velocity
vector of the nucleus can be parametrized as
v = v cξ cη − sξ sηcϕ1 , sξ cη + cξ sηcϕ1 , sηsϕ1 ,
in terms of two other angles η and ϕ1 which are uniformly
distributed in the intervals 0 < η < π and 0 < ϕ1 < 2π ,
respectively. We have chosen the nuclei and DM velocity
vectors to be aligned if η = 0.
If scattering is kinematically allowed, a Galilean
transformation is made to the frame in which the nucleus is
stationary and the DM velocity is wi,sc = wi(r ) − v. It is in
this reference frame that the recoil energy E R is defined and
its allowed range is in the interval [Emin, Emax], given by
Eq. (
4
). For a given recoil energy, the angle θ between the
outgoing DM velocity wf,sc and wi,sc can be calculated using
Eq. (
5
). Transforming back to the Sun’s rest frame, the DM
velocity after scattering is
wf = R(ϕ2)wf,sc + v,
where the operator R(ϕ2) rotates wf,sc around wi,sc by the
angle ϕ2, which is uniformly distributed in the interval 0 <
ϕ2 < 2π . The state in which the DM particle ends up in is
given by
(
21
)
(
22
)
(
23
)
(
24
)
(
25
)
(
26
)
Ef = Ei + 21 wf2 − wi2(r ) ,
Lf2 = r 2wf2 − (r · wf )2.
number density distribution is often approximated as an
isothermal Maxwell–Boltzmann distribution that can be
written as [12,17,28,59,66,67]
The fractional time spent in the shell with inner radius rinner
and outer radius router is calculated as
fiso(r ) = n0 exp(−r 2/rχ2 ) N ,
Tβ (ri ) =
T (rinner, router) ,
T (r−, r+)
(
27
)
(
28
)
(
29
)
(
30
)
(
31
)
where T (r1, r2) is given in Eq. (
11
). The shell widths are
chosen such that
rinner =
ri−1 + ri ,
2
router =
αβ matrix is then found using Monte-Carlo methods.
3.4 The radial number distribution function and the
annihilation rate
In order to calculate the annihilation rate of DM, detailed
knowledge of the radial number density distribution function
f (r ) is necessary. The annihilation rate for self-annihilating
DM is given by
ann =
σann(vrel)|vrel| f (r , v1) f (r , v2) d3v1 d3v2 d3r,
where vrel = v1 − v2 is the relative velocity between the
two colliding particles, f (r , v) is the phase-space density
distribution for DM, which has been normalized such that
f (r , v) d3v d3r = N ,
where N is the total number of captured DM particles. We
assume that the spatial distribution is spherically
symmetric. However, a possible consequence of a very low number
of scattering events of a particle over a long time may be a
preference for some orbital planes over others due to a
directional dependence of the flux of incoming DM particles. Such
a directional dependence could be caused by, for example,
the solar motion through the DM halo or anisotropies in the
galactic distribution of DM. This would introduce an
angular dependence in the number density distribution so that the
local DM distribution is increased in some regions relative
to the spherically symmetric case, which would increase the
overall annihilation rate with respect to the latter case. In the
following we assume that spatial spherical symmetry holds
when evaluating Eq. (
30
), keeping the above caveat at the
back of the mind.
When the mean free path of DM inside the Sun is much
larger than the solar radius (as considered here), the radial
where n0 = π −3/2 rχ−3. Assuming constant solar density and
temperature, the length scale rχ of the distribution is given
by
r 2 3kB Tc
χ = 2π Gρcmχ
,
where kB is the Boltzmann constant and G is the
gravitational constant. The bulk of the DM distribution is generally
located in such a centralized region that the density ρc and
temperature Tc can be approximated by the corresponding
values at the Sun’s center. For elastic DM, this assumption
is generally valid for scattering cross sections that yield
significant capture rates, see e.g. Refs. [60,68]. In this case, the
annihilation rate of DM, written in terms of the thermally
averaged annihilation cross section σannv , is
ann,iso = σannv
fi2so(r ) d3r.
It is clear that altering f (r , v) could significantly modify the
expected annihilation rate and that such a change is expected
for inelastic DM if the sub-dominant elastic scattering is
negligible. As our simulated distributions are given in E -L space,
we must map them onto r -v space. With the assumption of
many orbits per scattering, a DM particle in state α spends
the fractional time Tα(r ) at radius r . The radial distribution
function can thus be calculated by distributing all particles
of each state into all possible radii, weighed by the fractional
time spent at that radii:
f˜num(r ) =
fαTα(r ).
α
This is just the angular averaged distribution of the full
threedimensional spatial distribution, i.e.,
f˜num(r ) = r 2
fnum(r, θ , ϕ) sin(θ ) dθ dφ.
Since spherical symmetry has been assumed, the relation
above informs us that the three-dimensional distribution
function f (r ) can be found from f˜(r ) as fnum(r ) =
(4πr 2)−1 f˜num(r ). Unfortunately, solving Eq. (
10
) to find
fα can be computationally infeasible. Neglecting the
annihilation rate, the analytic solution for f is
f (t ) =
0
t
e (t−t ) C (t ) dt ,
(
32
)
(
33
)
(
34
)
(
35
)
(
36
)
(
37
)
standard Maxwellian model for the galactic velocity
distribution, with a shift to the solar frame,
f (u) = √π v2
u
−exp
exp
.
The solar velocity through the Milky Way, v , is taken to be
220 km/s, and the velocity dispersion v¯ = 270 km/s.
Unless otherwise stated, we use the DM-proton cross
section σχ p = 10−42 cm2. Both in the case of capture and
subsequent scattering of captured particles, we use the Helm form
factor [73]
F (q) = 3 j1(q R) e−q2s2/2,
q R
where q = √2m A E R is the momentum transfer in the
scattering process, j1 is the spherical Bessel function of the first
kind, and R is given by
where a possible time-dependence in the solar capture rate
has been taken into account. This allows us to find out if
the total capture of DM is in equilibrium with the loss due
to evaporation after a solar lifetime. It also permits us to
calculate the annihilation rate, and compare it with that of an
isothermal distribution.
The annihilation cross section times relative velocity can
be expanded as
σ (vrel)vrel = a + bvr2el + . . . ,
where a is non-zero for s-wave annihilation. Below we make
a simple comparison of s-wave annihilating DM from a
Boltzmann distribution to the distribution that we extract
from our numerical data. Since the relative velocities of DM
particles in the Sun are small, if both a, b = 0, a dominates
and σ (vrel)vrel is constant. In this case, we can trivially
calculate the integrals over v1 and v2 in Eq. (
30
). We then obtain
the ratio between the s-wave annihilation rate for the derived
distribution
ann,num(t ) = σannv
fn2um(r , t ) d3r,
and the isothermal distribution in Eq. (
34
), i.e.,
ann,num(t )
ann,iso
=
fn2um(r , t ) d3r
fi2so(r ) d3r
.
This provides us with a quantitative measure of how much
the annihilation rate is affected due to the change in the DM
distribution relative to the Maxwell–Boltzmann distribution.
Note that if the annihilation cross section is velocity
dependent the full phase-space distribution is required. Under
the assumption of a spherical distribution, it can be found as
follows. Any state that contributes to f˜(r ) at some radius r
gives a contribution to the velocity distribution at this radius.
The magnitude of the velocity v can be found from Eq. (
23
),
while the angle between the radial coordinate and the
velocity vector is ψ1 = sin−1(L/ri v). These two relations can
be used to extract the two-dimensional velocity distribution
f (ri , v, ψ1). It should be recognized that, from the symmetry
of the problem, we have that f (ri , v, θ ) = f (ri , v, −θ ).
4 Numerical results
We must now make some assumptions in order to proceed.
Specifically, we must define the galactic velocity
distribution and the local background density of DM, as well as
the nuclear form factor. We use the value nχ = nχ∗ =
0.2 GeV/cm3, which is half of the local DM density [69–
72]. In any example where elastic scattering is considered,
we assign nχ the value 0.4 GeV/cm3. We also assume the
(
38
)
(
39
)
(
40
)
(
41
)
(
42
)
(
43
)
R =
We use a = 0.52 fm, s = 0.9 fm, b = (1.23 A1/3 −
0.6) fm [74]. In order to speed up the computation time of
Cα and αβ we only take into account scattering on the
elements hydrogen, helium, nitrogen, oxygen, neon, and iron,
with radial abundances provided by the AGSS09ph solar
model [75]. This is an excellent approximation as the
abundances of the other solar elements are negligible and
contribute very little to scattering rates.
We use 100 individual states in E that are uniformly
distributed over all possible bound state energies. For every
discretization point in E , L is uniformly discretized in 100 states
between 0 and Lmax, which is the largest allowed angular
momentum for the given energy and can be found by
inverting Eq. (
7
). Therefore, in total, we use 104 states.
The following plots are shown in units of energy (E )
of G M /R , and in units of angular momentum (L) of
√G M R , where M is the solar mass. One can easily
check that these quantities naturally correspond to the
typical energy and angular momentum of a DM particle orbiting
around the centre of the Sun:
Eχ = mχ E = mχ
Lχ = mχ L = mχ
G M
R
G M
R
20
1/2
R
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.2
0.4
0.6
0.8
1
1.2
1.4
0.2
0.4
0.6
0.8
1
1.2
1.4
Furthermore, notice that Eχ ∼ δ for typical WIMP masses,
and therefore it is expected that the excited state can be
created by endothermic scatterings (see also the discussions in
Refs. [33, 44]).
4.1 The distribution of captured particles
In Fig. 1, we show the density of capture in E -L space,
normalized by its maximum value, for the elastic case, taking the
DM mass mχ = 5 GeV in the left panel and mχ = 100 GeV
in the right. For mχ = 5 GeV, capture is dominated by
helium and oxygen, followed by a slightly lower capture
rate by hydrogen and nitrogen. The concentration of
capture in the region centred slightly above E = −G M / R
and L ∼ 0.3√G M R is due to scattering on hydrogen,
which absorbs little recoil energy due to its low mass
relative to the DM. For helium, oxygen and nitrogen, capture
tends to be concentrated towards more strongly bound orbits,
with a preference for more circular orbits, i.e., larger L . For
mχ = 100 GeV, capture is primarily due to scattering on
helium and oxygen in almost equal parts, at a rate that is a
few times larger than that for capture by iron and neon. As
can be seen, capture is now concentrated towards states that
are much less bound. This is expected since the ability to
lose energy in a collision for heavier DM is hampered by the
relatively low mass of hydrogen and oxygen.
Moving on to inelastic scattering, Fig. 2 shows an
example of the density of capture of DM, with capture of χ in the
left plot and χ ∗ in the right plot. We use mχ = 100 GeV and
δ = 100 keV. The capture rate of χ particles, Cχ , is roughly
half as large as the capture of χ ∗, Cχ ∗ , which is expected
from previous studies, see e.g., Ref. [48]. An interesting
difference between the capture of χ and χ ∗ particles is the fact
that the former are captured into more tightly bound orbits
than the latter. This is due to two reasons, the first of which
is that a significant amount of kinetic energy is lost in the
endothermic process to produce the χ ∗. This loss of energy
reduces the form-factor suppression as the momentum
trans0.2
0.4
0.6
0.8
1
1.2
1.4
0.2
0.4
0.6
0.8
1
1.2
1.4
fer is not as large. Scattering takes place primarily on iron,
which due to its large mass is also a superb target for
absorbing recoil energy relative to the other elements. On the other
hand, energy being released in the exothermic case translates
into a larger form factor suppression and thus a preference
for scattering events in which the DM particle loses as little
energy as possible, leaving it less tightly bound. In this case,
capture occurs primarily due to scattering with helium nuclei
and to a lesser degree with oxygen, both of which are not very
efficient at absorbing recoil energy. Overall, the shape of the
region into which capture proceeds through exothermic
scatterings is similar to the elastic case with mχ = 100 GeV.
4.2 Time evolution of the distribution in E -L space
Having calculated the αβ matrix, the total scattering rate
from each state can be found. The case with mχ = 100 GeV
and δ = 100 keV is shown in Fig. 3, where we plot the
base 10 logarithm of the total scattering rate times the solar
lifetime, for χ → χ ∗ in the left plot and for χ ∗ → χ in the
right one.
The largest rate for χ → χ ∗ scattering is found for
the states with low angular momentum and medium energy.
These are the particles that have enough energy to travel
fairly far out from the solar center. When they fall back into
the solar center, they regain a significant amount of kinetic
energy, which allows endothermic scattering to take place.
Particles with larger energies spend more of their time
outside the Sun, which decreases their scattering rate. There are
also two regions, one at very low energies and one at very
high energies and large angular momenta, where scattering
does not take place at all. This can be explained by the fact
that the total kinetic energy in collisions taking place in most
of the E -L plane is supplied almost entirely by the DM
particle. The only region where this is not the case is in the very
low E region in which DM particle orbits are confined to the
solar center. These particles have very low velocities and the
energy of nuclei, even though the temperature is high, is not
sufficient to provide conditions under which up-scattering
can occur. At large E and large L , the nuclei are essentially
stationary and the DM particles always have low velocities
due to their circular orbits, leading to the conclusion that
up-scattering is kinematically disallowed also in this region.
Even if scattering is allowed, the rates are suppressed due to
the DM particles travelling on orbits in which they spend the
vast majority of their time outside the Sun.
The scattering of χ ∗ → χ is never kinematically
suppressed since the process is exothermic. The rates are thus
largest for particles that are confined to the solar center, i.e.,
in the low E and L region. The only suppression in the
scattering rate occurs for states at large E , which spend more time
in less dense regions. The extreme case is thus for very large
E and L , with highly circular orbits in the outer regions of the
Sun, where the density of targets is the lowest, the DM
velocity is small, and most of the time is spent outside the Sun.
Interestingly, it can also be seen that the rate for exothermic
scattering is larger than the rate for endothermic scattering
in the entire E -L plane. This indicates that the form-factor
suppression, which is larger for exothermic scattering, is not
as strong as the kinematic suppression of endothermic
scattering.
Next, we can take a sample of freshly captured DM
particles in a time t which is small enough for no additional
scattering to have occurred post capture. The time evolution is
then found by solving Eq. (
10
) neglecting additional capture
and annihilation, with the initial distribution f (0) = C t .
The distribution then evolves in time as
f (t ) = e t C
t .
Figure 4 shows the base 10 logarithm of the distribution at
various times for mχ = 5 GeV. The distribution accumulates
(
46
)
0.2
0.4
0.6
0.8
1
1.2
1.4
0.2
0.4
0.6
0.8
1
1.2
1.4
20
15
10
5
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
0.2
0.4
0.6
0.8
1
1.2
1.4
Fig. 4 The base 10 logarithm of the distribution of dark matter at various times t /t
f (0) = C t . We use elastic scattering with mχ = 5 GeV
= 10−10, 10−7, 10−4 and 1 for an initial distribution
20
15
10
5
into the lower region of the E -L space very rapidly. Taking
the scale into account, the distribution comes close to
equilibrium at t ∼ 10−7 t , at which point most particles in the
Sun have gathered in orbits with very low energies. As time
evolves further, there is a constant flow of the few remaining
particles at larger E down towards the lower energy orbits. It
is also interesting to note that evaporation is negligible over
a solar lifetime, i.e., Nχ (t )/Nχ (0) = 1.5
We now use Eq. (
35
) to translate the distribution in E -L
space into a radial distribution. The results for elastic
scattering are shown in Fig. 5 for the times t = 10−10 t (left),
t = 10−8 t (middle) and t = 10−6 t (right). The
distribution is compared to the isothermal one of Eq. (
32
), with
5 The figure at t = t is in good agreement with Fig. 3.1 in Ref. [65].
the angular degrees of freedom integrated over. We see that
the distribution has essentially reached equilibrium already
at t = 10−8 t , changing only slightly at t = 10−6 t . The
Boltzmann distribution gives a fairly accurate description of
the distribution, although the numerically computed one is
slightly shifted towards larger radii, and its peak is not as
pronounced.
Moving on to the case of inelastic DM, we again focus
our discussion on the illustrative case of mχ = 100 GeV and
δ = 100 keV. Figure 6 shows the base 10 logarithm of the
χ (χ ∗) distribution in the E -L plane in the left (right) plot
at various times. The majority of particles in the distribution
of χ have concentrated in the low E region rather quickly as
particles at higher E tend to lose energy when scattering and
hence fall down the gravitational well. However, there is now
also a region at large E and large L that contains a
significant number of χ particles. Rather than particles scattered
into this region from other bound orbits, these particles have
been primarily captured directly into it, although this is not
apparent from Fig. 2 due to the scale used. Even though the
capture rate may be low, Fig. 3 explains the relatively large
concentration of particles as up-scattering in this region is
kinematically forbidden.
We next obtain the radial distributions for inelastic
scattering and show the results in Fig. 7. At very early times, the
distribution extends up to large radii. At t = 10−9 t , a large
concentration starts to form, shown below r/ R 0.3. It
very slowly moves towards smaller radii, forming a
distribution centred at r/ R 0.1 at t = 10−5 t . However, even
at t = t the distribution has yet to reach a stationary state.
Another important observation is that the Boltzmann
distribution is now a very poor description of the final distribution.
This is entirely due to particles being trapped with no
possibility of scattering further, in particular those in the region
with low E , which are the ones that contribute to fnum(r ) at
smaller radii. Due to the circular orbits of particles at large
E and L , their contribution to the radial distribution is at
significantly larger radii (r 0.6 R ) than that shown in Fig. 7.
This contribution is not significant due to their low abundance
relative to the distribution close to the solar center. Since the
number of χ ∗ particles is completely negligible, the total DM
distribution is practically identical to the χ distribution.
Finally, it is interesting to compare the DM distributions
at t = t between the elastic and the inelastic cases. In Fig. 8
we show, for mχ = 100 GeV, the elastic case (left panel)
and the inelastic one with δ = 100 keV (right panel). As can
be seen, the distribution is (as expected) extremely
concentrated towards the central region of the Sun. One can also
observe that no evaporation has taken place. The distribution
is in fact concentrated into so few states that a reliable radial
distribution cannot be derived unless a significant increase in
the number of low E states used in the simulation is made.
This is a problem that also appears for inelastic DM when δ
is small enough, which prevents us from calculating the
annihilation rate for arbitrary low values of δ using the method
described here.
4.3 Annihilation and evaporation
In order to investigate the effects that the altered distribution
has on the annihilation rate, we use Eq. (
37
) to calculate the
number density functions at t = t . This distribution
contains information on the total number of particles in the Sun
and therefore provides a more realistic distribution. In fact,
sets of particles that are captured at different times are
distributed differently in the Sun at t = t and thus contribute
differently to the overall distribution. Of course, the
number of particles that have evaporated is also affected by the
amount of time passed since they were originally captured.
We now calculate distributions of inelastic DM at t = t
for different masses and cross sections. In the left panel Fig. 9
we show the ratio of the numerically calculated annihilation
rate to the isothermal one, computed using Eq. (
40
), as a
function of δ. Again, the problem of deriving radial
distribution functions for low values of δ due to the E -L
discretization used is encountered, which is why in Fig. 9 we
only show annihilation rates for larger values of δ, where
the problem is avoided. The results are shown for two
scattering cross sections: σχ p = 10−42 cm2 (solid lines) and
σχ p = 10−45 cm2 (dashed lines), and for three different DM
masses: mχ = 20 GeV (black lines), mχ = 100 GeV (blue
lines), and mχ = 500 GeV (red lines). The annihilation rate
is severely suppressed for large values of δ. The reason is
that, as δ increases, so do the regions in which additional
scattering of χ is kinetically forbidden, which in turn leads
to a more diluted DM distribution.
One can also observe that the suppression for σχ p =
10−42 cm2 is not as severe as the one for σχ p = 10−45 cm2.
This implies that the distribution does not reach a steady
state within a solar lifetime, as considering different
scattering cross sections is equivalent to considering different
times. Note that the number of captured particles depends
on the scattering cross section. However, N drops out from
the ratio num(t )/ iso and therefore the larger ratio for the
larger scattering cross section is not due to the total number
of DM particles, but only to their different distributions.
Next, we use the number density distributions, computed
with the same parameters as before, to calculate the actual
annihilation rate using Eq. (
39
). We assume s-wave
annihilation and assign the thermal averaged cross section the value
σannv = 3×10−26 cm3/s. The results are shown in the right
panel of Fig. 9, where the annihilation rates (solid lines) are
compared to the solar capture rates (dashed lines). We only
show the results for σχ p = 10−42 cm2, keeping in mind that
the smaller the scattering cross section the lower the
annihilation rate. In order to understand the comparison between
-5
0
0
-1
-2
-3
-4
-5
0
0
-1
-2
-3
-4
-5
0
0.2
C and num one should now recall two assumptions that
have been made. First, we assumed that the distribution is
spherically symmetric, so that the overall annihilation rate
would decrease if some orbital plane was preferred. Second,
we have assumed that no annihilation has taken place, which
means that the actual annihilation rate is overestimated. Thus
num should be regarded as an upper bound on the
annihilation rate under the assumption of a spherically symmetric
distribution.
As can be seen in the right panel of Fig. 9, unless δ is small,
the capture rate exceeds the annihilation rate. The decrease of
the annihilation rate with increasing δ is not surprising, as the
region of no scattering in E -L space becomes larger, locking
some particles in larger orbits, which in turn flattens out the
number density over a larger radial range. The upper bound
on the annihilation rate is larger than the capture rate when
δ 25 − 35 keV for mχ = 20 GeV and mχ = 100 GeV,
which implies that equilibrium may have been reached. The
Fig. 7 The radial distribution of χ particles, fnum (r ), at various times t = 10−9 t (left panel), t = 10−5 t (middle panel) and t = t (right
panel). We use mχ = 100 GeV and δ = 100 keV
0.05
rates (dashed lines) for the numerically computed distributions using
σχ p = 10−42 cm2 and σannv = 3 × 10−26 cm3/s. In both panels we
show results for mχ = 20 GeV (black), mχ = 100 GeV (blue) and
mχ = 500 GeV (red)
0.99
0.98
Fig. 10 Ratio between the total number of particles, after the
distribution has evolved for a solar lifetime, and the initial number of particles.
Results for different dark matter masses are shown: mχ = 20 GeV
(black), mχ = 100 GeV (red) and mχ = 500 GeV (blue)
same trend is apparent for mχ = 500 GeV, but the lack of
resolution of the distribution below δ = 50 keV prohibits this
from being shown. When num exceeds C , one should make
sure that equilibrium has truly been achieved by including
annihilations in the simulation.
Finally, we are interested in knowing how much
evaporation affects the total number of DM particles. In the elastic
case, it is generally accepted that DM particles with masses
below mχ ∼ 3 − 4 GeV evaporate before being able to
annihilate after they are captured [64, 76]. The situation is not at
all as clear in the case of inelastic DM, due to the absorption
and release of energy as the DM particle scatters back and
forth between the heavier and lighter states. Evaporation is
thus a cause for concern, in particular for sizeable δ.
In Fig. 10 we show the total number of particles N (t ) in
the Sun at t = t divided by its initial value N (0), where N (t )
is calculated using Eq. (
46
) for different DM mass values:
mχ = 20 GeV (black line), mχ = 100 GeV (red line) and
mχ = 500 GeV (blue line). We see that there is a value of the
splitting δmax(mχ ) where the evaporation rate reaches a
maximum. This value increases with the DM mass. For a given
DM mass, at splittings much smaller or much larger than
δmax, the evaporation rate vanishes or becomes negligible.
The former case, δ δmax, corresponds to the well-known
elastic limit, where evaporation is important only for very low
DM masses (mχ ∼ 3 − 4 GeV). In the latter case, δ δmax,
evaporation becomes suppressed due to two reasons. First,
halo χ particles are captured into states with, on average,
lower E as δ is increased, which reduces the likelihood that
the particles evaporate as they subsequently transition into
the lower states. Halo χ ∗ are captured into high E states,
but as these particles scatter back into χ ∗ in the first
interaction after being captured, it is extremely likely that they
drop to a significantly lower E state. This inhibits their
evaporation. Second, the χ scatterings may not be kinematically
allowed, and if they are, the resulting χ ∗ end up with very
little energy, and therefore in tightly bound orbits. As can be
observed, the evaporation rate is extremely low over a solar
lifetime, with at most 1 (
2
) % percent of particles evaporated
for mχ 100 (500) GeV.
5 Summary and conclusions
In this paper we have studied the evolution of the distribution
of inelastic DM in the Sun. We have presented the results
of a numerical simulation of the process of DM capture and
further scattering with nuclei in the Sun. We were particularly
interested in the case of inelastic DM with mass splittings δ
ranging from tens to hundreds of keV.
Our goal was to quantitatively study the process of
thermalization. In order for our simulation to be
computationally feasible, we have neglected annihilations in the
evolution equation. For definiteness, we have assumed that the DM
halo consists of equal populations of the two different states.6
We have evolved some initially captured distributions of χ
and χ ∗ in order to study the final distributions at a time equal
to the solar lifetime. We have found that χ ∗ are absent in the
final distribution, see Fig. 6. We obtained a χ -distribution
that has not reached a stationary state at a solar lifetime, and
that is far from being isothermal (Maxwell–Boltzmann) with
a temperature equal to that of the solar core, see Fig. 7, unlike
in the case of elastic scattering (c.f. Fig. 4).
By assuming spherical symmetry, we have also computed
an upper bound on the annihilation rate and found that it is
quite suppressed for splittings larger than a few tens of keV.
The exact suppression factor depends on the DM mass and
the scattering rate, see Fig. 9. However, it is quite robust to
state that, for mass splittings above roughly 50 keV , unless
deviations from spherical symmetry in the distribution are
present, the annihilation rate may be too small to yield a
large enough neutrino rate to be observable at Earth based
neutrino telescopes.
Similarly, by comparing the numerically obtained
annihilation rates to the equilibrium case in which they are equal
to the capture rates (see dashed lines in Fig. 9), we have
found that equilibrium is not reached unless the splittings
are below roughly 50 keV. For the DM masses considered,
equilibrium is achieved at values of δ(keV)/mχ (GeV) equal
to or smaller than a few. Interestingly, the larger the DM
masses, the greater the slopes of the annihilation rates (solid
lines). Therefore, for equilibrium to occur, it is expected that
there is a maximum splitting δ that is independent of the DM
mass. We have also studied evaporation and found that it
plays a less important role than previously thought as it stays
6 We defer the case in which χ ∗ is unstable, and the case of light
mediators, for future work.
safely below a few percent for the splittings and DM masses
considered, c.f. Fig. 10.
The most phenomenologically relevant implications of
this work are regarding the detection of neutrinos from DM
annihilations in the Sun. The most promising cases to have a
large annihilation rate are those where a non-negligible
elastic cross section is also present or in the case of very small
mass splittings ( O(
10
) keV).
Finally we would like to point out that it would also be
interesting to pursue numerical studies of scenarios with large
DM self-interactions, which would contribute to both capture
and evaporation, in which case one could also incorporate
DM annihilations for the same price in terms of simulation
complexity.
Acknowledgements We would like to thank Sergio Palomares-Ruiz
for useful discussions, Simon Velander for participation during the early
stages of this project, and Simon Israelsson for checks of the numerical
code. We are grateful to Martin White for proof reading the revised
version of the manuscript. The authors acknowledge the support from
the Spanish MINECO through the “Ramón y Cajal” programme
(RYC2015-18132) and through the Centro de Excelencia Severo Ochoa
Program under grant SEV-2016-0597 [M.B.], the Göran Gustafsson
foundation [M.B., S.C.], and the Australian Research Council through the
ARC Centre of Excellence for Particle Physics at the Terascale (CoEPP)
(CE110001104) [J.H.-G.]. S.C. also acknowledges the hospitality of
Instituto de Física Teórica (IFT) and the support from the Roland
Gustafsson foundation for theoretical physics during his stay at IFT
as part of this work was carried out.
Open Access This article is distributed under the terms of the Creative
Commons Attribution 4.0 International License (http://creativecomm
ons.org/licenses/by/4.0/), which permits unrestricted use, distribution,
and reproduction in any medium, provided you give appropriate credit
to the original author(s) and the source, provide a link to the Creative
Commons license, and indicate if changes were made.
Funded by SCOAP3.
48. M. Blennow, S. Clementz, J. Herrero-Garcia, JCAP 1604, 004
(2016). arXiv:1512.03317
49. J. Smolinsky, P. Tanedo, Phys. Rev. D95, 075015 (2017).
arXiv:1701.03168
50. M. McCullough, M. Fairbairn, Phys. Rev. D81, 083520 (2010).
arXiv:1001.2737
51. D. Hooper, D. Spolyar, A. Vallinotto, N.Y. Gnedin, Phys. Rev. D81,
103531 (2010). arXiv:1002.0005
52. M. Baryakhtar, J. Bramante, S.W. Li, T. Linden, N. Raj, Phys. Rev.
Lett. 119, 131801 (2017). arXiv:1704.01577
53. K. Schutz, T.R. Slatyer, JCAP 1501, 021 (2015). arXiv:1409.2867
54. Y. Zhang, Phys. Dark Univ. 15, 82 (2017). arXiv:1611.03492
55. M. Blennow, S. Clementz, J. Herrero-Garcia, JCAP 1703, 048
(2017). arXiv:1612.06681
56. F.J. Sanchez-Salcedo, Astrophys. J. 591, L107 (2003).
arXiv:astro-ph/0305496
57. M. Abdelqader, F. Melia, Mon. Not. Roy. Astron. Soc. 388, 1869
(2008). arXiv:0806.0602
58. N.F. Bell, A.J. Galea, R.R. Volkas, Phys. Rev. D83, 063504 (2011).
arXiv:1012.0067
59. A.R. Zentner, Phys. Rev. D80, 063501 (2009). arXiv:0907.3448
60. A. Widmark, JCAP 1705, 046 (2017). arXiv:1703.06878
61. L.J. Hall, T. Moroi, H. Murayama, Phys. Lett. B424, 305 (1998).
arXiv:hep-ph/9712515
62. Y. Cui, D.E. Morrissey, D. Poland, L. Randall, JHEP 05, 076
(2009). arXiv:0901.0557
63. S. Chang, N. Weiner, I. Yavin, Phys. Rev. D82, 125011 (2010).
arXiv:1007.4200
64. A. Gould, Astrophys. J. 321, 560 (1987b)
65. Z.-L. Liang, Y.-L. Wu, Z.-Q. Yang, Y.-F. Zhou, JCAP 1609, 018
(2016). arXiv:1606.02157
66. J. Faulkner, R.L. Gilliland, Astrophys. J. 299, 994 (1985)
67. R. Garani, S. Palomares-Ruiz, JCAP 1705, 007 (2017).
arXiv:1702.02768
68. M. Nauenberg, Phys. Rev. D36, 1080 (1987)
69. R. Catena, P. Ullio, JCAP 1008, 004 (2010). arXiv:0907.0018
70. J.I. Read, J. Phys. G41, 063101 (2014). arXiv:1404.1938
71. M. Pato, F. Iocco, G. Bertone, JCAP 1512, 001 (2015).
arXiv:1504.06324
72. S. Sivertsson, H. Silverwood, J.I. Read, G. Bertone, P. Steger,
Submitted to: Mon. Not. Roy. Astron. Soc. (2017), arXiv:1708.07836
73. R.H. Helm, Phys. Rev. 104, 1466 (1956)
74. J.D. Lewin, P.F. Smith, Astropart. Phys. 6, 87 (1996)
75. A. Serenelli, S. Basu, J.W. Ferguson, M. Asplund, Astrophys. J.
705, L123 (2009). arXiv:0909.2668
76. G. Busoni, A. De Simone, W.-C. Huang, JCAP 1307, 010 (2013).
arXiv:1305.1817
1. L. Bergström , Rept. Prog. Phys . 63 , 793 ( 2000 ). arXiv:hep-ph/0002126
2. G. Bertone , D. Hooper , J. Silk , Phys. Rept . 405 , 279 ( 2005 ). arXiv:hep-ph/0404175
3. J.L. Feng , Ann. Rev. Astron. Astrophys. 48 , 495 ( 2010 ). arXiv: 1003 . 0904
4. G. Bertone , D. Hooper , Submitted to: Rev. Mod. Phys. ( 2016 ), arXiv: 1605 . 04909
5. W.H. Press, D.N. Spergel , Astrophys. J. 296 , 679 ( 1985 )
6. A. Gould , Astrophys. J. 321 , 571 (1987a)
7. J. Silk , K.A. Olive , M. Srednicki , Phys. Rev. Lett . 55 , 257 ( 1985 )
8. L.M. Krauss , M. Srednicki , F. Wilczek , Phys. Rev. D33 , 2079 ( 1986 )
9. J.S. Hagelin , K.W. Ng , K.A. Olive , Phys. Lett. B180 , 375 ( 1986 )
10. T.K. Gaisser , G. Steigman, S. Tilav, Phys. Rev. D34 , 2206 ( 1986 )
11. M. Srednicki , K.A. Olive , J. Silk , Nucl. Phys. B279 , 804 ( 1987 )
12. K. Griest , D. Seckel , Nucl. Phys. B283 , 681 ( 1987 )
13. A.D. Avrorin et al. (Baikal), Astropart. Phys . 62 , 12 ( 2015 ). arXiv: 1405 . 3551
14. K. Choi et al. ( Super-Kamiokande ) , Phys. Rev. Lett . 114 , 141301 ( 2015 ), arXiv: 1503 . 04858
15. S. Adrian-Martinez et al. (ANTARES) , Phys. Lett. B759 , 69 ( 2016 ), arXiv: 1603 . 02228
16. M. G. Aartsen et al. (IceCube) , Eur. Phys. J. C77 , 146 ( 2017 ), arXiv: 1612 . 05949
17. D.N. Spergel , W.H. Press, Astrophys. J. 294 , 663 ( 1985 )
18. I.P. Lopes , J. Silk , S.H. Hansen , Mon. Not. Roy. Astron. Soc . 331 , 361 ( 2002 ). arXiv:astro-ph/0111530
19. A. Bottino , G. Fiorentini, N. Fornengo , B. Ricci , S. Scopel , F.L. Villante , Phys. Rev. D66 , 053005 ( 2002 ). arXiv:hep-ph/0206211
20. M.T. Frandsen , S. Sarkar, Phys. Rev. Lett . 105 , 011301 ( 2010 ). arXiv: 1003 . 4505
21. D.T. Cumberbatch, J. Guzik , J. Silk , L.S. Watson , S.M. West , Phys. Rev. D82 , 103503 ( 2010 ). arXiv: 1005 . 5102
22. M. Taoso , F. Iocco , G. Meynet, G. Bertone, P. Eggenberger, Phys. Rev. D82 , 083509 ( 2010 ). arXiv: 1005 . 5711
23. I. Lopes, J. Silk , Astrophys. J. 757 , 130 ( 2012 ). arXiv: 1209 . 3631
24. I. Lopes, K. Kadota , J. Silk , Astrophys. J. Lett. 780 , L15 ( 2014a ). arXiv:1310.0673
25. I. Lopes, P. Panci , J. Silk , Astrophys. J. 795 , 162 ( 2014b ). arXiv:1402.0682
26. A.C. Vincent , P. Scott , A. Serenelli, Phys. Rev. Lett . 114 , 081302 ( 2015a ). arXiv:1411.6626
27. A.C. Vincent , A. Serenelli , P. Scott , JCAP 1508 , 040 ( 2015b ). arXiv:1504.04378
28. M. Blennow , S. Clementz, JCAP 1508 , 036 ( 2015 ). arXiv: 1504 . 05813
29. C.-S. Chen , G.-L. Lin , Y.-H. Lin , Phys. Dark Univ. 14 , 35 ( 2016 ). arXiv: 1508 . 05263
30. A.C. Vincent , P. Scott , A. Serenelli, JCAP 1611 , 007 ( 2016 ). arXiv: 1605 . 06502
31. B. Geytenbeek , S. Rao , P. Scott , A. Serenelli , A.C. Vincent , M. White , A.G. Williams , JCAP 1703 , 029 ( 2017 ). arXiv: 1610 . 06737
32. G. Busoni, A. De Simone , P. Scott , A.C. Vincent , JCAP 1710 , 037 ( 2017 ). arXiv: 1703 . 07784
33. D. Tucker-Smith , N. Weiner , Phys. Rev. D64 , 043502 ( 2001 ). arXiv:hep-ph/0101138
34. R. Bernabei et al. (DAMA , LIBRA), Eur. Phys. J. C67 , 39 ( 2010 ). arXiv: 1002 . 1028
35. R. Abusaidi et al. (CDMS) , Phys. Rev. Lett . 84 , 5699 ( 2000 ). arXiv:astro-ph/0002471
36. D. S. Akerib et al. (LUX) , Phys. Rev. Lett . 118 , 021303 ( 2017 ), arXiv: 1608 . 07648
37. X. Cui et al. (PandaX-II) , Phys. Rev. Lett . 119 , 181302 ( 2017 ), arXiv: 1708 . 06917
38. E. Aprile et al. (XENON) , Phys. Rev. Lett . 119 , 181301 ( 2017 ), arXiv: 1705 . 06655
39. N. Bozorgnia , J. Herrero-Garcia , T. Schwetz , J. Zupan, JCAP 1307 , 049 ( 2013 ). arXiv: 1305 . 3575
40. D. Yu . Akimov et al. (ZEPLIN-III), Phys . Lett. B692 , 180 ( 2010 ). arXiv: 1003 . 5626
41. Z. Ahmed et al. (CDMS, CDMS-II) , Phys. Rev. D83 , 112002 ( 2011 ). arXiv: 1012 . 5078
42. E. Aprile et al. ( XENON100 ), Phys. Rev. D84 , 061101 ( 2011 ). arXiv: 1104 . 3121
43. X. Chen et al. (PandaX-II) , Phys. Rev. D96 , 102007 ( 2017 ), arXiv: 1708 . 05825
44. S. Nussinov , L.-T. Wang, I. Yavin , JCAP 0908 , 037 ( 2009 ). arXiv: 0905 . 1333
45. A. Menon , R. Morris , A. Pierce , N. Weiner , Phys. Rev. D82 , 015011 ( 2010 ). arXiv:0905 .1847
46. J. Shu , P.-F. Yin , S.-H. Zhu , Phys. Rev. D81 , 123519 ( 2010 ). arXiv: 1001 . 1076
47. M. McCullough , L. Randall , JCAP 1310 , 058 ( 2013 ). arXiv: 1307 . 4095