On the evaporation of solar dark matter: spin-independent effective operators

As a part of the effort to investigate the implications of dark matter (DM)-nucleon effective interactions on the solar DM detection, in this paper we focus on the evaporation of the solar DM for a set of the DM-nucleon spin-independent (SI) effective operators. In order to put the evaluation of the evaporation rate on a more reliable ground, we calculate the non-thermal distribution of the solar DM using the Monte Carlo methods, rather than adopting the Maxwellian approximation. We then specify relevant signal parameter spaces for the solar DM detection for various SI effective operators. Based on the analysis, we determine the minimum DM masses for which the DM-nucleon coupling strengths can be probed from the solar neutrino observations. As an interesting application, our investigation also shows that evaporation effect can not be neglectd in a recent proposal aiming to solve the solar abundance problem by invoking the momentum-dependent asymmetric DM in the Sun.


INTRODUCTION
As the nearest celestial body that is well understood and is capable of stimulating and responding to the phenomena associated with the Dark Matter (DM), the Sun is presumed to be an ideal host for the DM detection. For one thing its deep gravitational well attracts and traps the Galactic DM particles through the scatter off solar elements, if there exists a DM-nucleon interaction at the weak scale. For another thing these captured DM particles may accumulate in the solar core and subsequently annihilate to primary and secondary high energy neutrino flux that escape from the dense solar plasma, leaving a smoking-gun for their presence in the Sun. At present, a number of terrestrial neutrino detection projects such as IceCube [1,2], Super-Kamiokande [3], Baikal Neutrino Project [4] and ANTARES [5] are dedicated to such observational mission.
In general, the neutrino flux at the detector location is related to the solar DM annihilation through the following schematic relation: where d is the Sun-Earth distance, dΦ ν /dE ν and dN ν /dE ν represent the neutrino differential flux at the Earth and the neutrino energy spectrum per DM annihilation event in the Sun, respectively. The total annihilation rate Γ A can be expressed in terms of the number of the trapped DM particles N χ : where A denotes twice the annihilation rate of a pair of DM particles. The evolution of the solar DM number N χ is depicted with the following equation: which involves the DM capture (evaporation) rate C (E ) by scattering off atomic nuclei in the Sun, as well as the annihilation rate A . Eq. (1.3) has an analytic solution with τ e = C A + E 2 /4 −1/2 (1.5) the time scale for the capture, evaporation and annihilation processes to equilibrate. Once the equilibrium is reached at the present day, i.e., tanh (t /τ e ) 1, with t = 4.5 × 10 9 yr being the solar age, the annihilation output Γ A also reaches its maximum value. As will be shown in Sec. 3 3.2, a GeV increment in the DM mass parameter results in 1 ∼ 2 orders of magnitude reduction in the evaporation rate E in the few-GeV region. Thus depending on the ratio E 2 / (C A ), such equilibrium can be categorized into two different scenarios: (1) 1, that's when the evaporation effect can be neglected and the equilibrium is between annihilation and solar capture, i.e., Γ A C /2, so we can either determine or constrain the strength of the DM-nucleon interaction from solar neutrino observation; 1, under this circumstance evaporation overwhelms annihilation for the DM depletion, and the balance between evaporation and solar capture yields Γ A A C 2 / 2 E 2 , which not only implies a heavy suppression of the neutrino flux, but also prevents us from drawing the coupling strength of the DM-nucleon interaction from the possible observed signals.
Therefore, from the theoretical point of view it is interesting to pin down the parameter space where the neutrino observation is relevant for the DM detection. Conventionally, such purpose is fulfilled with a characteristic quantity, the evaporation mass m evap , which is defined with equation E (m evap ) = t −1 for the given DM-nucleon coupling. Above the evaporation mass one can safely assume that the capture-annihilation equilibrium is reached.
The key point of the problem is to calculate the distribution of the solar DM. While in Ref. [6,7] authors adopts a Maxwellian distribution to describe the non-thermal equilibrium between the solar DM particles and solar nuclei, the studies in Refs. [8,9] indicate a deviation from the Maxwellian form, in a manner that the actual velocity distribution is suppressed at the tail and tends to be anisotropic at large radius. Such deviation can be attributed to the fact that the energetic collisions that send the DM particles into high orbits occur predominantly near the hot core of the Sun, so as a result one expects a lower angular momentum distribution for the high-energy orbits. In order to well describe the physical processes such as evaporation and energy transfer of the solar DM, an accurate description of the tail of the velocity distribution is necessary.
In addition, since the evaporation mass has been studied thoroughly in the literature under the assumption of a constant DM-nucleon cross section [6,[8][9][10], the quest to the extend the discussion to a broader set of DM-nucleon effective interaction operators naturally arises. For instance, it is tempting to evaluate the evaporation rate for the light asymmetric DM particle with a DM-nucleon scattering amplitude linearly proportional to the square of the transferred momentum q 2 , because while the authors of Refs. [11,12] manage to resolve the disagreement between the solar model and helioseismological data with preferred DM mass of 3 GeV and coupling strength of 10 −37 cm 2 , the evaporation effect is not included in their discussion. Given small DM masses as such, evaporation may no longer be neglected in the buildup of the solar DM, and a quantitative analysis is needed on this issue.
Thus, as a tentative study we investigate the implications of the non-relativistic spinindependent (SI) effective operators on the solar DM distribution and evaporation mass. The set of 15 Galilean invariant operators is introduced in Ref. [13] * as a comprehensive and convenient treatment for the DM-nucleus interaction in the DM direct detection. Following Ref. [9] we calculate the non-thermal distribution of the solar DM by Monte Carlo methods, and numerically compute evaporation rates for different SI DM-nucleus effective operators. Moreover, based on the calculated capture and evaporation rates, we also discuss the parameter space relevant for the DM detection. This paper is organized as follows. In Sec. 2  [13]. m N is the mass of the nucleon.
we take a brief review on the effective interaction between the DM particle and nucleus. In Sec. 3 we calculate the solar DM distribution and evaporation rate for various SI DM-nucleus interaction operators, and discuss relevant implications for the high-energy solar neutrino signals. Some interesting discussions are arranged in Sec. 4.

EFFECTIVE INTERACTION BETWEEN DM AND NUCLEUS
We discuss the DM-nucleus scattering at low-energy scale in the context of the nonrelativistic (NR) effective interaction theory [13,[15][16][17][18], in which a set of linearly independent operators listed in Tab. 2.1 can be generated from the following five Hermitian operators: q is the transferred momentum from nucleon to the DM particle in a collision, and the transverse velocity is defined asv ⊥ = v + q/ (2µ N ), which satisfies q ·v ⊥ = 0 for the on-shell process, where v = v χ,i − v N,i is the relative initial velocity between the DM particle and nucleon, and µ N = m χ m N / (m χ + m N ) is the reduced mass of the system.Ŝ χ andŜ N are the spins of the DM particle and the nucleon, respectively. While the operators presented in Tab. 2.1 exhaust all the possible NR reduction of the Lorentz invariant spin-1/2 DM-nucleon interaction, up to corresponding coefficients dependent on the Galilean invariant scalar q 2 , in this study we shall investigate the implication of all the SI operatorsÔ 1 ,Ô 5 ,Ô 8 andÔ 11 † for the DM evaporation mass. Since the atomic nucleus is a composite of bound nucleons, its structural effect has to †Ô 2 is out of consideration because it will not be induced as the leading order term in non-relativistic expansion from the relativistic operators, unless there exists significant fine tuning that leads to a delicate cancellation among the leading pieces [13]. be taken into consideration in the analysis of the DM-nucleus interaction. Interestingly, in addition to the conventional nuclear form factor that describes the mass distribution within a nucleus, other types of DM and nuclear response functions arise from various underlying DM-nucleon interactions. For example, the operatorv ⊥ can be divided into the centre-ofmass and the relative motion components aŝ where µ A is reduced mass of the DM-nucleus system, v N,i (v N,f ) and v A,i (v A,f ) denote the initial (final) velocities of the constituent nucleon and the whole nucleus, respectively. Whilev ⊥ A ≡ v χ,i − v A,i + q/ (2µ A ) represents the nucleus transverse velocity, the latter term 1 2 x N in coordinate space, and gives rise to a nuclear response function (∆ response in Ref. [13]) associated with the nuclear orbital angular momentum in the long-wavelength limit. Nevertheless, compared with the conventional form factor that corresponds to W M in Refs. [15][16][17][18], response functions coming from the nuclear intrinsic motion (W ∆ in Refs. [15][16][17][18]) can be safely neglected if the isospin symmetry is respected. This is a direct observation from the nuclear response functions provided in Ref. [19]: isoscalar response functions (µ A /m N ) 2 W 00 ∆ are much smaller than W 00 M for the unpaired solar elements (e.g., 14 N, 23 Na and 27 Al). Not even to mention that these W ∆ responses associated with the unpaired elements suffer significant abundance suppression in the Sun.
Therefore, assuming the DM particle couples to the proton and neutron with equal strengths, the effects of response ∆ can be neglected for operatorsÔ 5 andÔ 8 , and hence we simply utilize the conventional Helm form factor to account for the nuclear internal structure, when investigating the implication of various SI interactions on the DM evaporation on a case-by-case basis. As a result, the DM-nucleus differential cross section for operators i = 1, 5, 8, 11 can be expressed in terms of the transferred momentum q as follows where c i carrying a dimension of mass −2 is the nucleon coupling constant for operatorÔ i , A is the atomic number of the target nucleus A, v rel = v χ,i − v A,i is the relative incoming velocity of the DM-nucleus system, and P i (v 2 rel , q 2 ) is the corresponding DM response function listed explicitly in Tab. 2.2 . In Tab. 2.2 , j χ represents the spin of the DM particle, and v ⊥2 the Helm form factor, with j 1 (x) = sin (x) /x 2 −cos (x) /x being the Bessel spherical function of the first kind, R 1 = R 2 0 − 5s 2 with R 0 1.23 A 1/3 fm, and s 1 fm [20].

DISTRIBUTION AND EVAPORATION OF SOLAR DM
In this section we will discuss the distribution and evaporation of the solar DM. Since the evaporation occurs predominantly at the high end of the velocity distribution, its evaluation relies on an accurate description thereof. We determine the solar DM distribution by solving the Boltzmann equation in a numerical way, and then separately calculate the evaporation rate for various effective SI DM-nucleon interaction operators. Now we delve into the details.

high end of the velocity distribution in the Sun
To date, there are two effective strategies in literature for determining the solar DM distribution. In the "Brownian motion" method that is pioneered by the author of Ref. [8], the distribution sample is obtained by simulating the motion of a single DM particle wandering in the Sun ‡ . While the "Brownian motion" method is efficient in describing the bulk of the velocity distribution, it turns impractical in computing the tail of the distribution for which a huge and uneconomical base of event samples is required to generate sufficient statistics. Therefore in order to determine the distribution of the solar DM, we resort to essentially the same method as the one outlined in Ref. [9].
Here we take a brief introduction to the methodology. Our discussion begins with the assumption that the presence of the solar DM does not bring any significant impact on the solar structure, i.e., the feedback from the accumulating DM particles is assumed to be negligible. The Boltzmann equation is linear due to the absence of the DM self-interaction, and can be further simplified as the following master equation if expressed with a convenient choice of parameters E (total energy per unit mass) and L (angular momentum per unit mass) [9]: where f (E, L) is the distribution function of the solar DM, and S (E, L; E , L ) represents the scattering matrix element for transition process (E, L) → (E , L ). In fact, to fully describe the physical state of the bound DM particle we still need an extra parameter, say, ‡ See Appendix A in Ref. [21] for an example. a temporal parameter τ , to label the position in the periodic orbit defined by energy and angular momentum. However, we approximate both the distribution function and scattering matrix elements as independent of parameter τ in Eq. (3.1). The reason is because a small DM-nucleus cross section, or equivalently, a large mean free path leads to a slowly increasing probability for a renewal collision, which implies an insensitive reliance of the distribution and scattering matrix on parameter τ .
The scattering matrix S (E, L; E , L ) is determined with simulation approach and the weighting method is adopted to facilitate the computation. Specifically speaking, we first calculate the probability for a trapped DM particle to collide with the solar elements on its trajectory at a fixed time interval ∆t, and then as a weight this probability is multiplied with the tally of the simulating transition events, so as to evaluate the scattering matrix in a more efficient manner. The numerical integration of the bound DM orbits is based on the Standard Sun Model (SSM) GS98 [22] and 5 solar elements H, 4 He, 14 N, 16 O and 56 Fe are included in the simulation of the DM-nucleus scattering. With random numbers that help pick out both the colliding solar element and its velocity, as well as the scattering angle in the centre-of-mass (CM) frame, we determine the outgoing state of the scattered DM particle after a coordinate transformation back to the solar reference. Further details of the discussion on the thermal collision are arranged in Appendix A.
It is also worth mentioning that in principle all kinetically allowed states of (E, L), including both the bound and unbound states that are connected to each other through capture and evaporation, should be involved in Eq. (3.1) for a realistic description of the solar DM. In practice, however, we model the captured DM particles as a closed system; that is to say, the number of the solar DM particles is assumed to be conserved within a timescale comparable to the relaxation time of the system, and the transitions are confined to only the Finally, by convoluting f χ (E, L) with φ EL ( r, v χ ), the distribution function of radius r and velocity v χ for orbit (E, L) , we obtain the DM distribution function For illustration, we present the distribution function of radius r after integrating out velocity v χ and vice versa for the orbit E = −1.225, L = 0.124 in Fig. 3.2.
Although the Maxwellian form of DM velocity distribution fails to describe the tail of the actual velocity distribution, as mentioned in Sec. 1, it suffices to approximate the bulk of the non-thermal distribution, on which physical processes such as DM annihilation can be evaluated easily and accurately. The approximate thermal distribution is expressed as f th ∝ exp (−m χ E/T χ ), with the effective temperature parameter T χ . T χ is determined by the demand that there be no net energy transfer from the solar nuclei to the shuttling DM particles once the steady state has been achieved, a requirement corresponds to the following energy-moment equation [6]:   where m A and n A (r) are the mass and the local number density of element A, T (r) is the temperature within the Sun, and V (r) is the gravitational potential as the function of radius r. In Tab. 3.3 shown is the effective temperature T χ for some benchmark DM masses from 1 GeV to 100 GeV. For a DM particle weighing tens of GeV, the effective temperature T χ can be approximated as the solar centre temperature T (0). For contrast, we compare the simulated velocity distribution f χ to the approximate thermal one f th in Fig. 3.3 and Fig. 3.4 for effective operatorsÔ 1 ,Ô 5 ,Ô 8 andÔ 11 in terms of the ratio f χ /f th . To estimate the errors that propagate from the simulated scattering matrices, we also present the standard deviations of the discrete limiting distributions for each set of parameters (Ô i , m χ ) in Fig. 3.3 and Fig. 3.4. Since in simulation the transitions are restricted to only the bound states, the DM velocity v χ stretches to no further than the escape velocity at the solar core v esc (0) ≈ 3.17. Echoing the studies in Refs. [8,9], while the ratio f χ /f th turns out to be suppressed at the high end of the velocity distribution, such suppression tends to be more significant for larger DM masses.

evaporation, capture and the minimum testable mass of the solar DM
In Ref. [9], the author provided a thorough discussion on the DM evaporation, under the assumption of a constant DM-nucleon cross section, which corresponds to the operatorÔ 1 in the context of the effective operators. Now we extends the discussion to include other SI effective operatorsÔ 5 ,Ô 8 andÔ 11 . Our interest are focused on the scenario in which the Sun is optically thin to the DM particles, so an evaporation event is counted once the speed of scattered DM particle exceeds the local escape velocity. For large DM-nucleus cross section, the blocking effect due to multiple collisions has to be taken into consideration, which turns out to heavily suppress the evaporation [23]. However, as will be shown later, the DM direct detections disfavour the coupling parameters relevant for the optically thick regime for these SI effective operators. As a consequence, a large optical depth for the solar DM particles amounts to a satisfactory approximation within the scope of this work.
Following Ref. [9], we start with the quantity R A (w → v) which represents the possibility of a DM particle with initial velocity w scattered to final velocity v by nucleus A in a unit volume, where dσ χA (|w − u A |) /dv is the differential cross section for the DM-nucleus system, which depends on their relative velocity w − u A , and · · · denotes the average over the thermal velocity distribution of element A. The Maxwellian distribution f A (u A ) is written as where u 0 = 2 T /m A . For the purpose of concision, we postpone the explicit expression of Eq. (3.4) to Appendix B. Next, given DM velocity w and the escape velocity v esc , the evaporation rate in the unit volume can be written as where the summation is taken over all solar elements. Finally, by convoluting Ω + (w |v esc ) with DM distribution f χ (r, w) determined from simulation, we express the DM evaporation rate as follows where Ω + (w |v esc ) depends on the radial coordinate r through the distributions of solar nuclei and the escape velocity v esc (r), which are both described with the SSM GS98 [22]. Given j χ = 1/2, the evaporation rate for various SI effective operators are expressed with the following fitting functions: which approximate the numerical results with an accuracy better than 10% in the DM mass range 2 ≤ m χ ≤ 5 GeV. For the sake of convenience, we invoke the DM-nucleon cross section σ p = c 2 1 µ 2 N /π instead of coupling parameter c 1 in Eq. (3.8a). Here we take a short review of the solar capture rate C and the annihilation coefficient A . The standard procedure for evaluating the DM capture rate C is developed in the literature [24][25][26]. Given the Galactic DM distribution unperturbed by solar influence, we first derive the collision event rate using the Liouville theorem and angular momentum conservation in the solar central force field, and by demanding the momentum transfer be large enough for the capture, we then extract the capture rate out of the total collision event rate. While discussions on capture rates for various DM-nucleon effective operators can be found in Refs. [19,27], here we present the numerical results for j χ = 1/2 in the DM mass range 2 GeV ≤ m χ ≤ 5 GeV as the following fitting functions dependent on the DM mass x = (m χ /1 GeV): (3.9d) In above evaluation of the capture rates, we adopt the isothermal DM halo model with a local density ρ χ = 0.3 GeV · cm −3 and a Maxwellilan velocity distribution with the dispersion v 0 = 220 km · s −1 , truncated at the Galactic escape velocity of 544 km · s −1 .
The annihilation coefficient A can be expressed in terms of the thermal cross section σv and the effective occupied volume of the solar DM V eff as the following: Now we are ready to explore the parameter space where the solar neutrino observational approach is effective for the DM detection, putting our intuitive discussion in Sec. 1 onto concrete computation. On one hand, as mentioned in Sec. 1, to ensure the full strength of the neutrino flux it is required that tanh (t /τ e ) 1, for which we adopt the criterion t /τ e 3.0. On the other hand, to specify the parameter region for the annihilation-and evaporationdominated scenarios, we set the criteria as E 2 / (4 C A ) ≤ 0.1 and E 2 / (4 C A ) ≥ 10, respectively. For concreteness, in Fig. 3.5 we show the relevant parameter regions for the annihilation-and evaporation-dominated regimes for SI effective operatorsÔ 1 ,Ô 5 ,Ô 8 and O 11 , by assuming the canonical s-wave thermal annihilation cross section σv = 3 × 10 −26 cm 2 , although the p-wave annihilation is also possible. Also shown in Fig. 3.5 (in yellow dashed lines) are the 90% C. L. upper limits on the DM-nucleon couping strengths imposed by the second run of the CDMSlite [28], which are derived using the Poisson statistics based on the event spectrum, signal efficiency, and detector resolution presented in Ref. [28], along with the astrophysical parameters consistent with the calculation of the capture rate. The new CDMSlite constraints are strong enough for narrowing our investigation to the optically thin regime. To illustrate this, takingÔ 1 for example, we note that the upper bound of σ p ≈ 10 −39 cm 2 corresponds to a mean free path at the solar centre l χ (0) = A n A (0) σ χA −1 ≈ 10 R , with σ χA the DM-nucleus cross section. So the assumption of a large optical depth is justified.
Given above quantitative analysis, we are able to draw clear boundaries among different signal topologies. For instance, for the effective interactionÔ 1 with a DM-nucleon cross section σ p = 10 −40 cm 2 , the assumption of an equilibrium between capture and annihilation is only valid for a DM particle heavier than 2.96 GeV, while for a DM mass smaller than 2.67 GeV, one can no longer extract the coupling strength of the DM-nucleon interaction from the observed neutrino flux, because the number of DM particles N χ = C /E becomes independent of cross section σ p [10] § . In addition, if the DM-nucleon cross section σ p is smaller roughly than 10 −44 cm 2 , the equilibrium among capture, evaporation and annihilation has not yet been achieved at the present day. As a consequence, the signal flux is suppressed and the unsaturated number of the solar DM (Eq. (1.4)) needs to be specified for neutrino telescopes to determine or constrain the coupling strength (see, e.g., Ref. [30]).

DISCUSSIONS
As mentioned in Sec. 1, authors of Refs. [11,12] introduce the weakly interacting asymmetric dark matter (ADM) with generalised form factors in an attempt to solve the solar abundance problem. Without annihilation the ADM may accumulate to such amount that their presence can slightly affect the solar structure. Assuming the evaporation rate is zero, it is found that the following SI interaction between a 3 GeV ADM and nucleon gives the best result: where the coupling σ 0 = 10 −37 cm 2 , and the reference momentum q 0 = 40 MeV. The translation between the contexts of the generalised form factor and the effective operatorÔ 11 is realised through the relation which gives c 11 = 1.87 × 10 −3 GeV −2 for j χ = 1/2. For the best-fit parameters given above, we calculate the evolution of the solar DM with and without evaporation in Fig. 4.1. It is evident that the presence of evaporation significantly constrain the increment of the DM number N χ and freezes it at a number O (10 4 ) smaller than the value without evaporation, which indicates an inconsistency for the model in Eq. (4.1) to alleviate the discrepancies between the SSM and helioseismological observables. Note that although we evaluate the evaporation rate by neglecting the interplay between the accumulated DM population and solar nuclei background, our calculation still holds in the ADM scenario because the relevant effects only result in minor changes in the solar structure. It should be also note that such inconsistency has been confirmed by the DM direct detection from the experimental aspect: CRESST-II ruled out this particular model at 90% C.L. [31]. In order to evade the constraints from the direct detection, the same authors of Refs. [11,12] recently propose a spin-dependent (SD) v 2 interaction as an alternative solution in Ref. [32]. We leave the discussion on the relevant evaporation effect in the SD scenario for future work.
Finally, we discuss a subtlety underlying the methodology applied to calculate the steady distribution f χ (E, L) in Sec. 3, i.e., to what extent the Markov chain approach describes the realistic evolution of the solar DM distribution, considering that both the replenishment and the leakage of DM particles are not reflected in the master equation Eq. (3.1). To this end, we explicitly write down the differential increment of the solar DM number in a time step δt, where vector ξ T = (ξ 1 , ξ 2 , · · · , ξ n ) and ξ T = (ξ 1 , ξ 2 , · · · , ξ n ) denotes the normalised probability for the n states at time t and t + δt, respectively, and η T = (η 1 , η 2 , · · · , η n ) represents the distribution for the newly captured DM particles in time interval δt. The Markov transition matrix S is expressed as describes the leakage due to evaporation, with e i being the evaporation rate for the i-th state. It is evident from Eq. (4.3) that the equilibrium distribution of the Markov chain ξ eq which satisfies the equation S · ξ eq = ξ eq well approximates the realistic distribution so long as the fractional change of DM number is negligible in the relaxation time δt = t relax , i.e., Therefore, for a time step δt t relax , it is reasonable to assume that solar DM equilibrates to its limit distribution instantaneously, and the descriptions of the distribution and the total number of the solar DM decouple and thus can be treated separately. Under such circumstance, one determines the evaporation rate using the steady distribution function and in turn integrates Eq. (1.3) to obtain the number of the solar DM in a self-consistent way. Note that for simplicity the annihilation is not included in our discussion, which however, will not cause any loss of generality of our conclusion. where is implicitly dependent on time once the DM trajectory is determined. The Galilean invariant σ (|w − u A |) can be obtained by integrating the differential cross section in Eq. (2.3). However, it should be noted that for DM mass around a few GeV, the typical momentum transfer in the thermal collision is of order of MeV, so we can neglect the Helm form factor for the bound DM scattering process. Here we take operatorÔ 11 as a specific example to illustrate how to calculate λ. First it is not difficult to obtain the cross section and then we input the v 2 rel reliance into integration in Eq. (A.2) as follows For simplicity we omit the summation notation over various solar elements A.
The analytic integration is performed using Mathematica. So once the DM particle motion is specified, the collision probability can be evaluated explicitly with Eq. (A.1). As an illustration, a segment of the solar DM trajectory is shown in Fig. A.1. Similar depiction is presented in Ref. [33], where the bound orbit is calculated using an analytic approximation for the solar potential.
Appendix B: calculation of the scattering event rate In this appendix we provide a detailed discussion on the scattering event rate R A (w → v) at which a DM particle scatters from initial velocity w to final one v, off a thermal bath composed of element A per unit volume. Except for a few notations, our discussion follows closely the original calculation in Refs. [9,34]. In short, after a coordinate transformation from the solar system to the CM system, Eq. (3.4) is expressed as an integration over the transformed coordinates (s, t) as the following: where η + A ≡ 1 + η A ≡ 1 + m χ /m A , s = (m χ w + m A u A ) / (m A + m χ ) and t = m A (w − u A ) / (m A + m χ ) are the CM velocity and the DM incoming velocity in the CM frame, respectively. u 2 A = η + A s 2 + η A η + A t 2 − η A w 2 , M is the relevant scattering amplitude dependent on (s, t) through the transferred momentum q = m χ (t − t), with t the DM outgoing velocity in the CM frame, and Θ is the Heaviside step function. By illustrating the relevant kinetic relation in Fig. B.1, we express the term |M| which also leads to the substitution for variables (s, t) in the expression of u 2 A : Applying these substitutions to Eq. (B.1) and assuming v > w for evaporation, we have (B.11) In practice, we simply numerically calculate Eq. (B.10) for various DM-nucleon effective interactions, rather than finding an analytic expression as has been done for the simplest caseÔ 1 in Ref. [9].