The distribution of inelastic dark matter in the Sun

The European Physical Journal C, May 2018

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 (\({\mathcal {O}} (10-100)\,\hbox {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 (\(\lesssim 50\,\hbox { 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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

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, 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

This is a preview of a remote PDF:

Mattias Blennow, Stefan Clementz, Juan Herrero-Garcia. The distribution of inelastic dark matter in the Sun, The European Physical Journal C, 2018, 386, DOI: 10.1140/epjc/s10052-018-5863-4