On the evolution process of two-component dark matter in the Sun

We introduce dark matter (DM) evolution process in the Sun under a two-component DM (2DM) scenario. Both DM species $\chi$ and $\xi$ with masses heavier than 1 GeV are considered. In this picture, both species could be captured by the Sun through DM-nucleus scattering and DM self-scatterings, e.g. $\chi\chi$ and $\xi\xi$ collisions. In addition, the heterogeneous self-scattering due to $\chi$ and $\xi$ collision is essentially possible in any 2DM models. This new introduced scattering naturally weaves the evolution processes of the two DM species that was assumed to evolve independently. Moreover, the heterogeneous self-scattering enhances the number of DM being captured in the Sun mutually. This effect significantly exists in a broad range of DM mass spectrum. We have studied this phenomena and its implication for the solar-captured DM annihilation rate. It would be crucial to the DM indirect detection when the two masses are close. General formalism of the 2DM evolution in the Sun as well as its kinematics are studied.


Introduction
Dark matter (DM) composes five times as prevalent as ordinary matter, yet its particle nature is still elusive. The essence of DM is often portrayed as Weakly Interacting Massive Particle (WIMP) with one specie. Experiments built to probe the interaction between the Standard Model (SM) particles and DM are running in progress [1][2][3][4][5][6]. Besides, terrestrial neutrino detectors [7][8][9] and satellite detectors [10][11][12] are designed to detect the SM particle fluxes from the annihilation of DM. Though much more stringent constraints on the DM properties have been set, primary implication for DM is generally from its gravitational influence. The understanding of DM is still in its budding stage.
Nevertheless, no strong evidence indicates that there exists only one-component DM in the dark sector (DS). DM with n-component (nDM) is also a plausible option. Each one contributes relic abundance Ω α h 2 to the the total relic abundance Ω DM h 2 where Ω DM h 2 = This paper is structured as follows. In section 2, we briefly introduce the 2DM scenario including notations and general assumptions in this work. In section 3, the evolution equations for χ and ξ are given. Coupling terms from heterogeneous self-interaction are introduced in the equations. In section 4, we present all the rates of interaction. Physical implications are shown. In section 5, numerical results of the evolution equations are calculated along with the DM total annihilation rates for both species. If DM can annihilate to final state with SM particles, the DM total annihilation rate characterizes the intensity of such SM particle flux. Finally, we summarize in section 6. Mathematical derivations of the rates relative to the heterogeneous self-interaction are given in the appendix.
2 Remarks on the 2DM scenario

Brief review of the 2DM models
The 2DM models typically require two different discrete symmetries assigned to each DM to sustain their stability. However, simple extension of the SM group by a global discrete symmetry can be violated by gravity [58] or induce the cosmic defects that are not compatible with cosmological observations [59]. These problems can be evaded if one retains the wanted discrete symmetries by breaking a gauge group.
One of the simplest 2DM models is based on an Abelian U (1) d , where d refers to hidden charge, and then assign some integer quantum numbers n 1 and n 2 to the hidden scalar fields χ and ξ [60]. After χ and ξ fields develop vacuum expectation values, the U (1) d will break with the residual discrete symmetries Z n 1 ⊗ Z n 2 . Examples are Z 2 ⊗ Z 2 [16] or Z 2 ⊗ Z 4 [19]. As a result χ and ξ can be both stable and will be the DM candidates. Other interesting models are the DM can be the multiplet of vector bosons of some hidden non-Abelian gauge groups [14,15,17,18]. For example, a hidden SU (2) d gauge group with a hidden fundamental representation of scalar field φ, where a = 1, 2, 3, F aµν is the hidden field strength, D µ φ = ∂ µ φ − i g d 2 τ a · A aµ φ, g d is the hidden coupling constant and A µ is the hidden gauge field. After φ developing a vacuum expectation value, SU (2) d breaks and the corresponding hidden gauge bosons will become degenerate massive particles. In such case, a residual custodial SO(3) symmetry remains due to the fact of scalar field φ being the fundamental representation. A Z 2 ⊗ Z 2 discrete symmetry, will apply to the hidden gauge bosons. Therefore, A a µ can be stable and DM candidates. One can extend the non-Abelian SU (2) d to larger groups [17]. It is worth of mentioning that in such models the DM can interact with the SM sectors via Higgs portal, or gauge boson kinetic mixing, where B µν and X µν are the field strength of U (1) Y and U (1) d respectively. Therefore the DM annihilation final states are SM particles or new scalars if it is kinematically allowed. These DM particles behave as thermal WIMPs and can retain the observed DM relic abundance. Interesting phenomena provided by these DM particles, e.g. excess cosmic rays, direct search, cosmology and collider physics, can be found in refs. [14,15,17,18,[61][62][63][64][65][66][67] and references therein.

Notation conventions and general assumptions
Suppose the two DM species χ and ξ only differ in mass such that m χ = m ξ in general.
The corresponding relic abundances are Ω α h 2 where α = χ and ξ. Assuming DM as the thermal relic and its total abundance Ω DM h 2 is made up of Ω χ h 2 and Ω ξ h 2 , e.g. Ω DM h 2 = Ω χ h 2 + Ω ξ h 2 . Ratio of the two relic abundances can be defined by If the annihilation is dominated by s-wave process at freeze-out epoch, we have the relic abundance inversely proportional to its thermal relic annihilation cross section σv 0 that is given by [46] which is mass-independent up to logarithmic corrections. Hence, we can express the annihilation cross sections by r ρ : Therefore, with eqs. (2.6) and (2.7), we can further define an effective annihilation cross section σ eff v 0 for Ω DM h 2 by To produce thermal relic abundance Ω DM h 2 ≈ 0.12 [13], it is reasonable to assume In addition, the number of DM particles captured by the Sun is relevant to the local DM density ρ DM around our solar neighborhood. Thus, the relation that should hold. Without loss of generality, we let ρ α ∝ Ω α h 2 , thus r ρ = ρ ξ /ρ χ . We can rewrite the above identity as ] m χ = 100 GeV, m ξ = 10 GeV Figure 1. Local DM number density n α vs. r ρ . ρ DM = 0.3 GeV cm −3 is the local DM density near the solar neighborhood. A similar result can be also found in the figure 1 of ref. [22]. and the local DM number density n α = ρ α /m α is plotted in figure 1 versus r ρ . A similar plot can be found in the figure 1 of ref. [22] as well. When r ρ = 1, the DM with lighter mass has the higher number density. On the other hand, the number densities happen to be equal when r ρ = m ξ /m χ . In the later analysis, we assume the validations of eqs. (2.7), (2.8) and (2.10) pass to the present day. Once r ρ is assigned, σ α v and ρ α are specified consequently.
3 General formalism of dark matter evolution in the Sun

The 1DM evolution equation
When the Sun sweeps the DM halo, DM particles are attracted by the solar gravity. The subsequent scatterings with the solar nuclei and other DM particles already trapped in the Sun could happen. DM particles can be captured by the Sun when its final velocity is smaller than the escape velocity of the Sun after scattering. Alternatively, DM particles trapped inside the Sun will be kicked out if its final velocity after the scattering with the nuclei is larger than the escape velocity. The inclusion of DM self-interaction will also have effects on the capture and evaporation of DM particles in the Sun. Incorporating all these effects, the general equation describes the DM evolution process is given by where N DM is the DM number in the Sun, C c the rate at which DM is captured by the solar nuclei [47][48][49]52], C s the self-capture rate at which DM is captured due to scattering with other trapped DM inside the Sun [68], C e the evaporation rate due to DM-nucleus scattering [69], C a the annihilation and C se the self-evaporation rate that caused by DM-DM scattering [70]. However, recent study shows unless DM mass m DM 4 GeV, the evaporation effect is much inefficient even with the inclusion of C se [70]. Thus, in the absence of evaporations, eq. (3.1) reads along with an analytical solution where τ = 1/ C c C a + C 2 s /4 is the equilibrium timescale. When t τ , dN DM /dt = 0.

The 2DM evolution equations
On the other hand, eq. (3.1) only characterizes 1DM scenario. Once the second DM specie is included, additional evolution equation should be added. In this scenario, the selfinteractions are not only due to χχ and ξξ scatterings as well as χξ scattering. Thus, the evolution process is determined by The above equations are modified from eq. (3.1). Four additional coefficients are introduced. C χ(ξ)→ξ(χ) s denote the heterogeneous self-capture rates due to halo χ(ξ) scatters with trapped ξ(χ) in the Sun. C χ(ξ)→ξ(χ) se denote the heterogeneous self-evaporation rates due to the χ(ξ) scatters with ξ(χ) in the Sun. The rest are C α c the solar captures, C α s the self-capture rates, C α e the evaporation rates, C α se the self-evaporation rates and C α a the annihilation rates as those in the 1DM case.
Both eqs. (3.4a) and (3.4b) correlate together through the terms subject to N ξ in dN χ /dt and N χ in dN ξ /dt. The DM numbers of χ and of ξ in the Sun are mutually dependent. Without correlation terms, the evolution processes for both DM species are decoupled. Generally, evaporation is inefficient unless the DM mass is light enough, typically when m ev α 4 GeV. Even including the extra contribution C χ(ξ)→ξ(χ) se , we have numerically justified that this effect does not change the m ev α dramatically. Thus, when m α > 4 GeV, we can safely ignore the evaporation from eqs. (3.4a) and (3.4b). Therefore, But in the later numerical calculations, we will always use the general expressions eqs. (3.4a) and (3.4b). Note that the evolution equations have no analytical expressions on N α and τ in the 2DM scenario. However, approximated expressions can be obtained in certain situations. It will be discussed in section 5.

The solar capture rate
The solar capture rate due to DM-nucleus scattering can be numerically approximated as [55,71] C SI for spin-independent (SI) case and for spin-dependent (SD) case. ρ α is the DM local density andv = 270 km s −1 the DM velocity dispersion. σ SI,SD H,He is the DM-nucleus scattering cross section for hydrogen or helium. Taking proton mass m p is close to neutron mass m n . The DM-nucleus cross section σ A at which interaction is undergoing that is related to DM-nucleon cross section σ αp by for SI interaction interaction and for SD interaction, where A is the atomic number, m A the corresponding nucleus mass, J the total angular momentum of the nucleus and S p and S n the spin expectation values of proton and of neutron averaged over the entire nucleus [72][73][74][75][76][77]. To apply the above results, we have assumed χ and ξ obey the same Maxwell-Boltzman velocity distribution. The effect of uncertainties in velocity distributions to the capture rate is minor [78]. We note that σ αp is a model-dependent parameter in general. In addition, following earlier work [71], refined calculation on solar capture rate including the contributions from elements beyond hydrogen and helium can be found in refs. [49,52]. In ref. [52], constant scattering cross section as well as velocity-dependent and transfer-momentum-dependent cases are fully considered. We adopt the numerical procedure in ref. [52] to calculate C c in the 2DM scenario.

The self-capture rate
The self-capture happens when the halo DM particles scatter off the DM particles that are already trapped inside the Sun. Starting with halo χ particle captured by the trapped ξ particle in the Sun. Therefore, the coefficient of heterogeneous self-capture rate can be expressed as where dC χ→ξ s /dV is the heterogeneous self-capture rate in the Sun within a given shell. It is determined by m χ , m ξ and σ χξ where σ χξ is the heterogeneous self-scattering cross section. The analytical form of dC χ→ξ s /dV is given in eq. (A.11). Assuming the heat exchanges among DM particles are very efficient after capture. They will quickly reach the thermal equilibrium temperature T α = T . Therefore, n α (r) = n 0 α e −mαφ(r)/T where n 0 α is the DM number density in the solar core, φ(r) = r 0 GM (r )/r 2 dr and M (r ) the solar mass enclosed by radius r . The case for halo ξ particle captured by trapped χ particle is essentially the same. Simply swaps χ and ξ and replace all χ's parameters by ξ's.
The expression given in eq. (4.5) is generally for m χ = m ξ . Nevertheless, if the capture is due to the same DM specie, e.g. χχ or ξξ scatterings, it is done by letting m χ = m ξ and denotes as C α s . Such that we have a rather simple analytical expression [68]: where σ α is the self-scattering cross section, φ α 5.1 [47] a dimensionless average solar potential experienced by the captured DM within the Sun and v esc (R ) ≈ 632 km s −1 the Sun's escape velocity at surface. We will characterize σ α by in the later numerical analysis. Such relation appears to alleviate the diversity of galactic rotation curve in the presence of baryonic effect as well as small-scale structure problems [42]. Plot for C χ→ξ s in the m χ − m ξ plane is shown in figure 2.

The annihilation rate
When more and more DM particles accumulate in the Sun, the rate of annihilation becomes stronger. The coefficient of annihilation rate is expressed as where σv is the DM annihilation cross section. An approximation for C α a is given by [46] C α a = σ α v is the DM effective volume and T = 1.54 × 10 7 K the solar core temperature. Essentially, the DM temperature T is not necessary the same as the solar temperature T [52,79] and depends on m α . But the temperature deviation from T is generally small and has little impact on the final number of DM particles in the Sun. It is reasonable to impose T = T in our later discussion.

The evaporation and self-ejection rates
In the Sun, after the collision happens between two particles, if one gets velocity larger than the escape velocity v esc , it won't be captured. This effect is called evaporation. When the evaporation is caused by DM-nucleus scattering, an approximation is given in ref. [80] C α where v esc (0) = 1366 km s −1 is the escape velocity at solar core,r the mean DM distance from the solar center. The quantity Σ evap is a factor that relates to the DM-nucleus scattering cross section [80]. The expression C α e is valid when m α /m A > 1. In our calculation, we modified the numerical procedure done in ref. [52] for the 2DM scenario.
Likewise, evaporation can happen as a result of χ particle scattering off ξ particle in the Sun and vice versa. Hence we obtain coefficient of heterogeneous self-evaporation is given by 4πr 2 (dC χ→ξ se /dV )dr ( 4πr 2 n χ (r)dr)( 4πr 2 n ξ (r)dr) (4.11) where dC χ→ξ se /dV is the heterogeneous self-evaporation in the Sun within a given shell. It is a function of m χ , m ξ and σ χξ and its analytical form is expressed in eq. (A.18). If the self-evaporation is due to either χ or ξ itself, simply let m χ = m ξ with additional factor of 1/2 to avoid over counting. The plot of C χ→ξ se against m χ is given in figure 3. In addition, if a DM particle in the halo transports enough kinetic energy to the trapped DM particle in the Sun and leads to the trapped DM particle being ejected to the interstellar space. This is called the (heterogeneous) self-ejection. Numerical calculations σ χξ = 10 -24 cm 2 Figure 3. The self-evaporation coefficient C χ→ξ se . It is easily seen from the plot that when m χ > 4 GeV, the effect is small enough to ignore it.
show the self-ejection effect is always small, compared to the self-capture one. Thus, we can ignore this effect safely from the calculation. Discussion on self-ejection rate is presented in the appendix A.3.

Number of dark matter particles in the Sun
In the following analysis, when r ρ is assigned, the annihilation cross section σ α v and ρ α can be specified through eqs. (2.8) and (2.10). Thus, thermal relic abundance Ω DM h 2 ≈ 0.12 and ρ DM = 0.3 GeV cm −3 would be satisfied automatically. In the later numerical analysis, we have taken that σ SI αp = 10 −46 cm 2 as a benchmark value for SI case. It is slightly smaller than the most stringent value of LUX when the DM mass is roughly around 30 GeV [3]. For SD case, σ SD αp = 10 −42 cm 2 that is chosen not to violate the results from Super-K [8], PICO-60 [4] and IceCube [9]. The self-scattering cross section σ α is indicated from eq. (4.7) and we set σ χξ = 10 −23 cm 2 and 10 −24 cm 2 as the benchmark values. These two values are within 0.1 cm 2 g −1 ≤ σ DM /m DM ≤ 10 cm 2 g −1 in the whole interested mass range. The value of σ χξ tells us how χ and ξ intertwine during the evolution. The number of DM particles in the Sun, N α , is plotted in figure 4 versus time t. We have fixed m χ at 1000 GeV and calculated with m ξ = 100 GeV and 10 GeV. The case for σ χξ = 0 is labeled as decoupled for comparison.
Number of DM in the Sun, N α , is proportional to its local number density n α = ρ α /m α when reaches equilibrium stage. Suppose r ρ = 1 and m χ m ξ . We have n χ n ξ . Thus, χ affects little on the evolution of ξ. Hence, in the equilibrium stage, we could drop C ξ→χ s N eq in eq. (3.5b). In this way, simple expressions for N eq χ and N eq ξ can be given by are the correction factors due to the (heterogeneous) self-captures. We have verified eqs. (5.1a) and (5.1b) and they agree with numerical solutions of eqs. (3.5a) and (3.5b) very well after reaching the equilibrium state. See dot-dashed lines in figure 4. When n ξ n χ , ξ evolves solely in the Sun. But N eq χ is subject to a correction that is proportional to σ χξ N eq ξ . In addition, we take eq. (4.7) as the benchmark value of σ α . It results in a very strong self-interacting effect. Therefore, the DM numbers in the equilibrium state, N eq α ,  is mostly determined by the interactions in the DS. This agrees with the conclusion in ref. [79] for 1DM case.
On the other hand, equilibrium must achieves simultaneously for both DM species. When n ξ n χ , the equilibrium timescale can be determined by ξ solely. It is given by τ ξ = 1/ C ξ c C ξ a + (C ξ s ) 2 /4. Therefore, the role plays by σ χξ that is insignificant and can be omitted. In figure 5, equilibrium region for a given m ξ is indicated by its corresponding color contour. The place enclosed by the contour indicates t /τ < 1 as well as the nonequilibrium region.
Note that when m χ = m ξ , C χ→ξ s = C χ s and C ξ→χ s = C ξ s . The evolution equations eqs. (3.4a) and (3.4b) are degenerate. It can be considered as an 1DM scenario.

Implication for the dark matter total annihilation rate in the Sun
When an appreciated amount of DM particles accumulate in the solar core, the total annihilation rate 1 as a result of these particles is given by  for a given DM specie α. If it is in the equilibrium state, we can apply eqs. (5.1a) and (5.1b) and obtain where R χ,ξ is given in eq. (5.1c). The above equations assume ξ dominates the DM population over χ. Counter case is vice versa. For a more general discussion, unless specified, we will not assume which specie is dominant over the other. The plot of Γ χ versus m χ is shown in figure 6 with r ρ = 1. In this figure, we fixed m ξ = 10 GeV and 100 GeV while m χ runs from 1 GeV to 1000 GeV. In the above choice of parameters, both DM species are all in equilibrium state today. As a consequence of large interactions in the DS, the number of DM in the the equilibrium state is affected little from the DM-nucleus interaction. Results from SI and SD cases are both similar. Therefore, we focus on the SI case only in the following discussion. In figure 6, the lower panel shows the ratio between coupled and decoupled cases. Such ratio indicates how strong is the correction from σ χξ .
On the left panel of figure 6, we fixed m ξ = 10 GeV and Γ ξ is indicated by the blue line. When m χ > m ξ , it is true that n χ < n ξ . Hence ξ is the dominant specie in the  Sun and Γ ξ can be considered as independent of χ particles. But Γ χ is usually subject to a correction from ξ when m χ > m ξ . However, when m χ is close to m ξ , both numbers N χ,ξ are nearly equivalent. Mutual influence is strong in this region. Not only Γ χ is enhanced by ξ particles, as well as Γ ξ is increased by χ particles in the Sun. The ratio of correction is shown in the lower panel. A quick drop of Γ χ when m χ 4 GeV is due to the evaporation effect. The discussion is similar to the right figure of figure 6, instead of raising m ξ to 100 GeV. Again, the correction from ξ to Γ χ is not significant when m χ > m ξ here. Nonetheless, in the range m χ is smaller than m ξ (n χ > n ξ ), Γ ξ is subject to a correction from χ. Note that m χ = 100 GeV (m ξ = 10 GeV) of the left figure is the same as the right figure of m χ = 10 GeV (m ξ = 100 GeV). It can be realized from the symmetry of the evolution equations given in eqs. (3.5a) and (3.5b).
The case for r ρ = 1 is shown in figure 7. Parameters are the same as in the case of r ρ = 1 and the DM particles are also in the equilibrium state. We know that Γ α ∝ C 2 s /C a = n 2 / σv and the ratio Γ χ /Γ ξ ∝ (m ξ /m χ ) 2 /r 3 ρ in terms of eq. (2.7). Hence we can deduce that Γ χ /Γ ξ ∼ 1 when m χ ≈ 32m ξ for r ρ = 0.1. This statement is partially correct since the exact Γ α is subject to an extra correction factor from R α in eqs. (5.3a) and (5.3b). We have numerically verified the correction factor is roughly 3. Precisely speaking, when m χ ≈ 100m ξ , Γ χ /Γ ξ ∼ 1. From the upper panel left in figure 7, it is clearly seen that Γ χ /Γ ξ ∼ 1 when m χ ∼ 1000 GeV. This argument agrees with our numerical result well. Similarly, it applies to the case of m ξ = 100 GeV. In this case, Γ χ /Γ ξ ∼ 1 when m χ ≈ 10 4 GeV. For r ρ = 10, we can use the approximation above and obtain that Γ χ /Γ ξ ∼ 1 when m χ ≈ 0.01m ξ . Therefore, m χ = 0.1(1) GeV when m ξ = 10(100) GeV that we would have Γ χ /Γ ξ ∼ 1. However, one should bear in mind that for m χ 4 GeV all χ particles have evaporated already.

Summary
In this paper, we address the issue of 2DM evolution in the Sun. We consider a scenario, where the heterogeneous χξ self-scattering happens. Such interaction weaves the evolution processes for both DM species that was assumed to evolve independently. We found that when one DM specie is sub-dominant, its number of particles in the Sun is subject to a correction from the dominant specie. This correction always enhances the number of DM particles being captured. When the masses of the two DM species are close, the enhancement is mutual and has the largest impact. However, the sub-dominant specie in general has smaller total annihilation rate, the effect of heterogeneous self-capture would be tiny to the detection unless its annihilation final state is distinct from the dominant one.
Though the heterogeneous self-interaction causes extra self-evaporation and self-ejection, we have demonstrated that these negative effects are either small or it must happens when the DM mass is sufficient light, m ev α 4 GeV. Therefore, in most of the interested mass range that relates to our study, they can be safely ignored.
To summarize, we would like to point out that the heterogeneous self-interaction is a natural consequence of any 2DM or nDM models. This effect will eventually reflect in the DM annihilation rates. Potentially, if the DM annihilates to the SM particles in the final state, such signal could be detected in the terrestrial detectors. Therefore, the strength of the heterogeneous self-interaction could be probed. Moreover, any sign of such interaction could be considered as a possible existence of DM beyond one-component.

respectively.
A Derivations of the 2DM heterogeneous self-scattering rates A.1 The self-capture rate To the capture rate of different classes of particles has been fully discussed in refs. [47,48,68]. In this appendix, we only present the mathematical key point to derive the heterogeneous self-capture rate.
Following earlier works [48,68], the problem begins by considering capture in a spherical shell of material (solar interior) on which capture is happening of radius r and local escape velocity v esc (r). Now at an imaginary surface bounding a region of radius R, which the solar gravity is negligible at R. The DM flux goes inward across the surface is [81] where f (u) is the DM velocity distribution at infinity, J = Ru sin θ the angular momentum per unit mass and θ the angle relative to the radial direction. Taking Ω(w) is the rate at which a DM particle enters the shell r with velocity w = v 2 esc (r) + u 2 and scatters to velocity less than v esc (r). The probability of such a DM to be captured is [68] dP = Ω(w) w where Θ is the Heaviside step function. The differential rate of capture can be easily obtained by multiplying eqs. (A.1) and (A.2) then integrate over all angular momentum J 2 . Replacing dV = 4πr 2 dr we have, Thus, the total DM capture rate per unit shell volume is given by In the above equation, w depends on u explicitly. The remaining task is to determine Ω(w). The scattering in the shell is simply nσw, with the the scattering cross section σ and the target number density n. Practically we assume nearly isotropic and velocity-independent σ. The incoming particle with m χ and scatters off bounded particle with m ξ . In order to be captured, χ particle must loses a fractional of kinetic energy over the interval where µ and µ ± are expressed as and η 2 = 3(v /v) 2 /2, v = 220 km s −1 the solar moving velocity andv = 270 km s −1 the DM velocity dispersion. Therefore, the capture probability in each scattering is The rate of capture is simply the scattering rate n ξ σw times the capture probability p cap . Hence, With respect to the solar moving frame, we can expressed f (u) as in terms of the dimensionless variables x 2 = 3(u/v) 2 /2 and η 2 = 3(v /v) 2 /2. Integrating eq. (A.9) over u, we have where we have replaced σ by σ χξ to indicate the heterogeneous self-scattering cross section and Thus, the coefficient of heterogeneous self-capture rate is evaluated as where n ξ (r) is the number distribution of ξ particles in the Sun. The case for halo ξ particle scatters with solar trapped χ particle is essentially identical.

A.2 The self-evaporation rate
Self-evaporation happens when two DM particles collide, one gets velocity larger than the escape velocity v esc . Such calculation is similar to the evaporation between DM and nucleus presented in ref. [47]. Here we show the key to obtain the heterogeneous self-evaporation rate. To scatter a ξ particle from velocity w to v esc > w by χ particle, the rate is In order to evaluate eq. (A.17), we assumed the that T χ = T ξ = T and w c = v esc . 2 Therefore, dC χ→ξ se dV = 2 π 2T m ξ n ξ (r)n χ (r)σ χξ e −Ee/T −β + β − − 1 2µ X(β − , β + ) where E e = m ξ v 2 esc (r)/2. Therefore, we have the coefficient of the heterogeneous selfevaporation rate, C χ→ξ se = 4πr 2 (dC χ→ξ se /dV )dr ( 4πr 2 n χ (r)dr)( 4πr 2 n ξ (r)dr) .
(A. 19) When m χ = m ξ , it reduces to the 1DM case and a symmetric factor 1/2 should be introduced to avoid over counting.

A.3 The self-ejection rate
Once the incoming χ particle loses a fraction of energy ∆E/E > v 2 esc (r)/w 2 to a trapped particle ξ. The ξ particle will be ejected from the Sun. Following the derivation in the appendix A.1 but replacing p cap by the ejection probability [68] p ejec = µ 2 Thus, the rate of ejection, Integrating over the χ number distribution in the halo f (u) given in eq. (A.10) we have By changing of variable we get Therefore, dC χ→ξ ej dV = 4 √ π σ χξ n χ n ξ (r)v 3η e −K 2 + (e 4Kη K + − K − ) Our final result of the heterogeneous self-ejection rate is evaluated as C χ→ξ ej = 4πr 2 (C χ→ξ ej /dV )dr 4πr 2 n ξ (r)dr . (A.24) However, due to the large escape velocity in the Sun, such self-ejection effect is always insignificant comparing to other effects. Thus, we can safely ignore this correction in the DM evolution.