The distribution of inelastic dark matter in the Sun

If dark matter is composed of new particles, these may become captured after scattering with nuclei in the Sun, thermalise 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)$ keV) state it is unclear whether thermalisation occurs. One issue is that up-scattering from the lower mass state may be kinematically forbidden, at which point the thermalisation 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 thermalisation 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$ 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.


Introduction
A popular class of models for explaining a number of observations in astrophysical systems is that of particle dark matter (DM) [1][2][3][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 e-mail: m.blennow@csic.es b e-mail: scl@kth.se c e-mail: juan.herrero-garcia@adelaide.edu.au 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][8][9][10][11][12], with various collaborations having performed such searches with no positive detection [13][14][15][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][18][19][20][21][22][23][24][25][26][27][28][29][30][31][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][41][42][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][45][46][47][48][49]. See also Refs. [50][51][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][54][55], and in models where the more massive state is unstable [56][57][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.

Inelastic dark matter
There are various scenarios in which inelastic DM appears naturally [61][62][63]. In the simplest models inelastic DM consists of two states χ and χ * , with masses m χ < m χ * that satisfy 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, 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.

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.
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 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 E min (L) of a particle with a trajectory that intersects the Sun is given by 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 equatioṅ whereṄ ≡ d N/dt, C is the solar capture rate, E evap ∝ 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 (Ṅ = 0) of the evolution equation reads 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 equatioṅ 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ṙ = 0. The time it takes the particle to move between radius r 1 to radius r 2 can be found by isolatingṙ in Eq. (6), which leads to Integrating from r 1 = r − to r 2 = r + gives the time that a particle needs to complete half an orbit.

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] 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 = u 2 + v esc (r ) 2 is the velocity of the particle at radius r . The local halo number density of DM and its speed distribution at the location of the Sun enter explicitly through n χ and f (u). Note that the velocity distributionf ( u) is normalized such that where 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σ Here it has been assumed that the coupling between DM and nuclei is isospin conserving, leading to the A 2 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 2 w 2 . 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 w i can be reconstructed. When the DM particle scatters, its energy post-collision is known from where E i (f) is the energy of the incoming (outgoing) DM particle. The outgoing DM speed w f is known from the equa-tion above. The angle θ between the incoming and outgoing velocity vector is given by Eq. (5). There is also an azimuthal angle ϕ around w i 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 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.

Scattering among different states
When DM particles have been captured, occasional scattering with solar nuclei takes a particle initially in the β = (E i , 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] The velocities of nuclei in the Sun follow the Boltzmann distribution 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 r i . Under the assumption that DM particles complete many orbits between interactions, the rate at which particles in state β scatter at radius r i and end up in state α is given by Here, R β (r i ) is the total scattering rate at radius r i and T β (r i ) is the fraction of the orbital time that the particle spends inside the shell. The factor P β→α (r i ) is the probability that the particle ended up in the particular state α after it scattered at r i . Having calculated the above, the off-diagonal elements 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 α, which also includes evaporation. To find P β→α (r ), the DM velocity vector w i of a particle in state β = (E i , L i ) before scattering is required. In the Sun's rest frame, its magnitude is found to be We are free to choose a coordinate system in which r = (r, 0, 0). The angle between r and w i (r ) is given by ξ = sin −1 (L i /r w i (r )). Thus the DM velocity vector can be written as where we use s x = sin(x) and c x = cos(x). The velocity vector of the nucleus can be parametrized as 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 w i,sc = w i (r ) − v. It is in this reference frame that the recoil energy E R is defined and its allowed range is in the interval [E min , E max ], given by Eq. (4). For a given recoil energy, the angle θ between the outgoing DM velocity w f,sc and w i,sc can be calculated using Eq. (5). Transforming back to the Sun's rest frame, the DM velocity after scattering is where the operator R(ϕ 2 ) rotates w f,sc around w i,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 The fractional time spent in the shell with inner radius r inner and outer radius r outer is calculated as where T (r 1 , r 2 ) is given in Eq. (11). The shell widths are chosen such that The αβ matrix is then found using Monte-Carlo methods.

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 where v rel = v 1 − v 2 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 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 number density distribution is often approximated as an isothermal Maxwell-Boltzmann distribution that can be written as [12,17,28,59,66,67] where n 0 = π −3/2 r −3 χ . Assuming constant solar density and temperature, the length scale r χ of the distribution is given by where k B 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 T c 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 σ ann v , is 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: This is just the angular averaged distribution of the full threedimensional spatial distribution, i.e., Since spherical symmetry has been assumed, the relation above informs us that the three-dimensional distribution function f ( r ) can be found fromf (r ) as f num ( r ) = (4πr 2 ) −1f num (r ). Unfortunately, solving Eq. (10) to find f α can be computationally infeasible. Neglecting the annihilation rate, the analytic solution for f is 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 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 σ (v rel )v rel is constant. In this case, we can trivially calculate the integrals over v 1 and v 2 in Eq. (30). We then obtain the ratio between the s-wave annihilation rate for the derived distribution and the isothermal distribution in Eq. (34), i.e., ann,num (t) 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 tof (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/r i v). These two relations can be used to extract the two-dimensional velocity distribution f (r i , v, ψ 1 ). It should be recognized that, from the symmetry of the problem, we have that

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/cm 3 , which is half of the local DM density [69][70][71][72]. In any example where elastic scattering is considered, we assign n χ the value 0.4 GeV/cm 3 . We also assume the standard Maxwellian model for the galactic velocity distribution, with a shift to the solar frame, The solar velocity through the Milky Way, v , is taken to be 220 km/s, and the velocity dispersionv = 270 km/s. Unless otherwise stated, we use the DM-proton cross section σ χ p = 10 −42 cm 2 . Both in the case of capture and subsequent scattering of captured particles, we use the Helm form factor [73] where q = √ 2m A E R is the momentum transfer in the scattering process, j 1 is the spherical Bessel function of the first kind, and R is given by We use a = 0.52 fm, s = 0.9 fm, b = (1.23A 1/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 L max , which is the largest allowed angular momentum for the given energy and can be found by inverting Eq. (7). Therefore, in total, we use 10 4 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: 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]).

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 trans- 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.

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 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), 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 f num (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.

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 cm 2 (solid lines) and σ χ p = 10 −45 cm 2 (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 cm 2 is not as severe as the one for σ χ p = 10 −45 cm 2 . 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 σ ann v = 3×10 −26 cm 3 /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 cm 2 , keeping in mind that the smaller the scattering cross section the lower the annihilation rate. In order to understand the comparison between 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 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.

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