# Evaporation of dark matter from celestial bodies

Raghuveer Garani<EMAIL_ADDRESS>INFN Sezione di Firenze, Via G. Sansone 1,
I-50019 Sesto Fiorentino, Italy  Sergio Palomares-Ruiz<EMAIL_ADDRESS>Instituto de F\́text{i}sica Corpuscular (IFIC), Universitat de València –
CSIC, Parc Cient\́text{i}fic, C/ Catedrático José Beltrán, 2, E-46980 Paterna,
Spain

###### Abstract

Scatterings of galactic dark matter (DM) particles with the constituents of
celestial bodies could result in their accumulation within these objects.
Nevertheless, the finite temperature of the medium sets a minimum mass, the
evaporation mass, that DM particles must have in order to remain trapped. DM
particles below this mass are very likely to scatter to speeds higher than the
escape velocity, so they would be kicked out of the capturing object and
escape. Here, we compute the DM evaporation mass for all spherical celestial
bodies in hydrostatic equilibrium, spanning the mass range
$[10^{-10}-10^{2}]\leavevmode\nobreak\ M_{\odot}$, for constant scattering
cross sections and $s$–wave annihilations. We illustrate the critical
importance of the exponential tail of the evaporation rate, which has not
always been appreciated in recent literature, and obtain a robust result: for
the geometric value of the scattering cross section and for interactions with
nucleons, at the local galactic position, the DM evaporation mass for all
spherical celestial bodies in hydrostatic equilibrium is approximately given
by $E_{c}/T_{\chi}\sim 30$, where $E_{c}$ is the escape energy of DM particles
at the core of the object and $T_{\chi}$ is their temperature. In that case,
the minimum value of the DM evaporation mass is obtained for super-Jupiters
and brown dwarfs, $m_{\rm evap}\simeq 0.7$ GeV. For other values of the
scattering cross section, the DM evaporation mass only varies by a factor
smaller than three within the range $10^{-41}\leavevmode\nobreak\
\textrm{cm}^{2}\leq\sigma_{p}\leq 10^{-31}\leavevmode\nobreak\
\textrm{cm}^{2}$, where $\sigma_{p}$ is the spin-independent DM-nucleon
scattering cross section. Its dependence on parameters such as the galactic DM
density and velocity, or the scattering and annihilation cross sections is
only logarithmic, and details on the density and temperature profiles of
celestial bodies have also a small impact.

##  I Introduction

The possibility of capture of dark matter (DM) particles by celestial bodies
has a long history dating back to the mid 1980s [1, 2, 3, 4]. Already in the
late 1970s, some of the potential effects of DM annihilations and scatterings
on energy transport within the Sun were realized [5] and then further proposed
to alleviate the so-called solar neutrino problem [6, 7, 8, 9, 10].
Nevertheless, these early works did not consider the process of DM capture,
but assumed the required amount of accumulated DM particles. The process of
capture of galactic DM particles by the Sun and the Earth was first studied in
Refs. [1, 2, 3, 4], which set the ground of further calculations.

The scattering of DM particles in galactic halos with the nuclei (or
electrons) of celestial objects could bring those particles into close orbits
and finally result in their gravitational capture within the objects. DM
particles would undergo further scatterings and, except for very small cross
sections, thermalize in a short period of time, so that this process can be
approximately considered as instantaneous [11, 12, 13, 14, 15, 16]. In
addition to the Sun [1, 17, 18, 19, 20, 21, 22, 2, 3], capture of DM particles
and its potential observational consequences have also been considered in
other celestial bodies, such as the Earth [18, 19, 21, 4, 23, 24], other solar
system planets and satellites [19, 25, 26, 27, 28, 29, 30, 31], exoplanets
[29, 32], brown dwarfs [33, 32], main-sequence [7, 34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46] and post-main sequence [47, 48, 49, 50] stars, and
compact objects such as white dwarfs [51, 52, 53, 54, 55, 56] and neutron
stars [57, 58, 59, 53, 60, 61, 62].

Below a minimum DM mass, scatterings of thermalized DM particles off ambient
targets (nuclei or electrons) would boost them very efficiently to speeds
above the escape velocity, so that they would evaporate from the capturing
celestial body. Therefore, only for DM masses above the evaporation mass,
effects of the capture of DM particles could have any impact on the evolution
of the objects or could produce any other observable signature (the exception
being the potential signal of evaporated DM particles at low-threshold direct
detection experiments [63]). The DM evaporation mass has been computed at
different levels of detail for the case of the Sun [5, 6, 7, 19, 2, 21, 3, 10,
64, 65, 66, 67], but has been less studied for other celestial bodies. Some
examples, however, are the calculations for the Earth [18, 19, 23, 31], the
Moon [31], Mars [30], giant planets [32, 68], brown dwarfs [33, 32], main-
sequence stars [33, 69], horizontal-branch stars [49, 70], and neutron stars
[71, 72].

The process of DM evaporation from celestial bodies depends on the temperature
of potential scatterers in the medium, which are assumed to be in thermal
equilibrium. Part of the finite thermal energy of nuclei (or electrons) would
be transferred to DM particles via elastic scatterings, which could then end
up with energies above the local escape energy. Unless further scatterings
take place, these DM particles would not be gravitationally bound any more and
thus would escape from the object. Obviously, in the limit of zero
temperature, the evaporation rate is zero.

Light thermalized DM particles typically have speeds higher than heavier
particles, so they are more likely to end up with speeds higher than the local
escape velocity after scattering off target particles in the medium. The
combination of the DM velocity distribution along with the probability to
escape after scattering results in the evaporation rate having an exponential
dependence [2, 3, 70]. The key quantity is the exponent. In the thin regime
(long mean free path), evaporation mostly occurs close to the center of the
object,111In the thick regime (short mean free path), DM particles approach
local thermodynamic equilibrium and evaporation occurs mostly in a shell
closer to the surface of the object [70], where the temperature, but also the
escape velocity, are lower. In any case, qualitatively, the discussion is
analogous. and the exponent is $E_{c}/T_{\chi}$, the ratio of the escape
energy of DM particles at the center of the object to their temperature
(assuming an isothermal distribution). The smaller the exponent, the larger
the number of DM particles close to the escape velocity and the higher the
probability to gain enough energy to escape, and hence, the higher the
evaporation rate. It is a well known fact since the early papers on the topic
that, in order to efficiently suppress DM evaporation in the Sun,
$E_{c}/T_{\chi}\simeq 30$ [2]. This implies that the DM evaporation mass is
set along the exponential tail of the evaporation rate [6, 21, 2, 3]. For the
case of the Sun, $m_{\rm evap}\simeq 3$ GeV [64, 66, 67], whereas for the
Earth, $m_{\rm evap}\simeq 12$ GeV [18, 19, 23, 31] and for the Moon, $m_{\rm
evap}\simeq 70$ GeV [31].222Note that Ref. [31] quotes $m_{\rm evap}\simeq 40$
GeV. Nevertheless, to obtain that DM evaporation mass for the Moon, the core
temperature was set to $T_{c}=700$ K, instead of $T_{c}=1700$ K, as written in
the text. For the latter core temperature, $m_{\rm evap}\simeq 70$ GeV. These
are the values we use as realistic model references. Thermal models of the
interior of the Moon predict a core temperature above 1000 K, but it could
reach a value of up to 2000 K [73].

In this work, we compute the DM evaporation mass for a wide range of celestial
bodies, from the smallest objects with spherical shape that can be in
hydrostatic equilibrium (small satellites and dwarf planets), $M\simeq
10^{-10}\leavevmode\nobreak\ M_{\odot}$ [74, 75], to the most massive main-
sequence stars, $M\simeq 100\leavevmode\nobreak\ M_{\odot}$. In addition, we
also discuss the DM evaporation mass for post-main-sequence stars, white
dwarfs and neutron stars. We consider a range of DM-nucleon (momentum and
velocity independent) scattering cross sections that covers ten decades,
$10^{-41}\leavevmode\nobreak\ \textrm{cm}^{2}\leq\sigma_{p}\leq
10^{-31}\leavevmode\nobreak\ \textrm{cm}^{2}$, and which runs over the thin
and thick regimes.

Recently, searches of effects of the accretion of sub-GeV DM particles in
giant planets and brown dwarfs have been proposed [32, 76, 68]. The ideas rely
on estimates of the DM evaporation mass for those objects as low as a few MeV,
based on their relatively large size and cool temperatures. Likewise, a low DM
evaporation mass for the Earth and Mars has also been claimed [30]. Here, we
show that these calculations, which assume constant scattering cross sections,
neglect the crucial exponential tail of the evaporation rate and underestimate
the DM evaporation mass by at least one order of magnitude, and we argue that
the suggested implications for DM masses below the correctly evaluated
(properly accounting for the exponential tail) DM evaporation mass do not
apply.

This paper is organized as follows. In Section II, we introduce the basic
ingredients required for the calculation of the DM evaporation mass: the
capture, annihilation and evaporation rates. Then, we describe in some detail
the process of DM evaporation in a generic celestial body and explain the
importance of the exponential tails in the determination of the DM evaporation
mass. In Section III, we describe the main average properties, relevant for
the calculation of the DM evaporation mass, of all celestial bodies we
consider in this work, spanning twelve orders of magnitude in mass,
$10^{-10}\leavevmode\nobreak\ M_{\odot}\leq M\leq 10^{2}\leavevmode\nobreak\
M_{\odot}$. The main results of this work are presented in Section IV, where
we show the value of the DM evaporation mass for all these objects and
describe its dependence on several parameters. Finally, in Section V, we
summarize our results and draw our conclusions.

##  II Basics of dark matter evaporation

DM particles in a galactic halo could scatter off the material of celestial
bodies and lose enough energy such that their velocity ends up smaller than
the escape velocity of the object, thus becoming gravitationally captured
within. The evolution of the total number of DM particles accumulated in a
celestial body, $N_{\chi}(t)$, is given by

$\frac{dN_{\chi}(t)}{dt}=\mathcal{C}-\mathcal{A}\,N^{2}_{\chi}(t)-\mathcal{E}\,N_{\chi}(t)\leavevmode\nobreak\
,$ (1)

where $\mathcal{C}$, $\mathcal{A}$ and $\mathcal{E}$ are the DM capture,
annihilation and evaporation rates, which throughout this work, are taken to
be constant in time and are defined in the next subsections. In the above
equation, canonical two-body annihilation processes are assumed and DM self-
interactions are not considered. The solution reads [21, 2]

$N_{\chi}(t)=\mathcal{C}\,\tau_{\rm eq}\,\frac{\tanh(\kappa\,t/\tau_{\rm
eq})}{\kappa+\frac{1}{2}\,\mathcal{E}\,\tau_{\rm
eq}\,\tanh(\kappa\,t/\tau_{\rm eq})}\leavevmode\nobreak\ ,$ (2)

where $\tau_{\rm eq}=1/\sqrt{\mathcal{A}\,\mathcal{C}}$ is the equilibration
time scale in the absence of evaporation and
$\kappa=\sqrt{1+(\mathcal{E}\,\tau_{\rm eq}/2)^{2}}$. Equilibrium is reached
when $\kappa\,t\gg\tau_{\rm eq}$. Before this occurs,
$N_{\chi}\simeq\mathcal{C}\,t$ and, as discussed below, the DM evaporation
mass grows with time. Once equilibrium is attained, if evaporation dominates
($\kappa\gg 1$), $N_{\chi}\simeq\mathcal{C}/\mathcal{E}$, and the number of
accumulated DM particles decreases exponentially with decreasing mass (see
below). If evaporation is not efficient, $\kappa\simeq 1$ (yet, $t\gg\tau_{\rm
eq}$), the number of DM particles is given by
$N_{\chi}\simeq\mathcal{C}\,\tau_{\rm eq}=\sqrt{\mathcal{C}/\mathcal{A}}$.
Bearing this in mind, the DM evaporation mass, $m_{\rm evap}$, can be defined
as [64] (see Refs. [5, 6, 21, 70, 66] for other definitions)

$\left|N_{\chi}(t;m_{\rm evap})-\frac{\mathcal{C}(m_{\rm
evap})}{\mathcal{E}(m_{\rm evap})}\right|\equiv 0.1\,N_{\chi}(t;m_{\rm
evap})\leavevmode\nobreak\ ,$ (3)

which is the mass for which the number of captured DM particles reduces to the
solution in the limit of evaporation dominance, $\mathcal{C}/\mathcal{E}$, at
the 10% level. When equilibrium is reached, this condition becomes time
independent and simplifies to $\mathcal{E}(m_{\rm evap})\,\tau_{\rm eq}(m_{\rm
evap})=1/\sqrt{0.11}$ [66]. The calculation of this minimum mass of DM
particles that can get efficiently trapped in celestial bodies depends on the
properties of the capturing object, and this is the main goal of this work.

###  II.1 Main inputs

In order to set the ground for the discussion about the evaporation process
and the minimum DM mass that can be efficiently captured by different
celestial bodies, here we define the three relevant rates introduced above in
terms of the DM properties (mass and scattering cross section off nuclei) in
the vicinity of the celestial objects (which depend on the local DM density
and velocity distribution), the properties of the capturing body (mass and
size, density and temperature profiles, and composition) and the kinematics of
elastic DM-nuclei scatterings. The general features of the DM capture process
by non-degenerate non-relativistic targets are well known in the literature.
Here, we only sketch the main ingredients and refer the reader to other works
for further details [1, 2, 3, 4, 77, 70, 24, 78, 79, 64, 80, 81, 65, 66, 82,
67, 83, 84].

We consider a celestial body of mass $M$ and radius $R$, composed of non-
degenerate non-relativistic elements with mass $m_{i}$, number density
$n_{i}(r)$ which are assumed to be in local thermodynamic equilibrium
described by a Maxwell-Boltzmann distribution of temperature $T(r)$. The total
number of targets $i$ is $N_{i}=(X_{i}/A_{i})\,M$, where $X_{i}$ and $A_{i}$
are the mass fraction and the mass of element $i$.

We consider the DM population in the vicinity of this object with mass
$m_{\chi}$, density $\rho_{\chi}$ and with velocity distribution $f_{v_{\rm
cb}}(u_{\chi})$, as seen by an observer moving at speed $v_{\rm cb}$ with
respect to the galactic rest frame,

$f_{v_{\rm cb}}(u_{\chi})=\frac{1}{2}\,\int_{-1}^{1}f_{\rm
gal}\left(\sqrt{u_{\chi}^{2}+v_{\rm cb}^{2}+2\,u_{\chi}\,v_{\rm
cb}\,\cos{\theta}}\,\right)\,\text{d}\cos{\theta}=\sqrt{\frac{3}{2\pi}}\frac{u_{\chi}}{v_{\rm
cb}\,v_{d}}\,\left(e^{-\frac{3\,(u_{\chi}-v_{\rm
cb})^{2}}{2\,v_{d}^{2}}}-e^{-\frac{3\,(u_{\chi}+v_{\rm
cb})^{2}}{2\,v_{d}^{2}}}\right)\leavevmode\nobreak\ ,$ (4)

where $u_{\chi}$ is the DM velocity at infinity in the celestial body’s rest
frame, $\cos{\theta}$ is the angle between the DM and the object velocities
and $f_{\rm gal}(u_{\rm gal})$ is the DM velocity distribution in the galactic
rest frame, which is usually assumed to be a Maxwell-Boltzmann distribution
(the so-called standard halo model) with dispersion velocity
$v_{d}=\sqrt{3/2}\,v_{\rm cb}$ at a given position in the halo. Throughout
this paper, we use $v_{d}=270\leavevmode\nobreak\ \textrm{km/s}$ and
$\rho_{\chi}=0.4\leavevmode\nobreak\ \textrm{GeV/cm}^{3}$.

These DM particles may eventually interact with the thermal distribution of
the object’s targets, with differential scattering cross section
$\text{d}\sigma_{i}/\text{d}v$. In this work, we focus on constant (i.e.,
momentum and velocity independent) total scattering cross sections, which can
be written in terms of the DM-proton scattering cross section (or the DM-
neutron scattering cross section, that we assume to be equal), for spin-
independent (SI) or spin-dependent (SD) interactions, at zero momentum
transfer, as

$\displaystyle\sigma_{i}^{\rm SI}$ $\displaystyle=$
$\displaystyle\left(\frac{\tilde{\mu}_{A_{i}}}{\tilde{\mu}_{p}}\right)^{2}\,A_{i}^{2}\,\sigma_{p}^{\rm
SI}\leavevmode\nobreak\ ,$ (5) $\displaystyle\sigma_{i}^{\rm SD}$
$\displaystyle=$
$\displaystyle\left(\frac{\tilde{\mu}_{A_{i}}}{\tilde{\mu}_{p}}\right)^{2}\,\frac{4\,(J_{i}+1)}{3\,J_{i}}\,\left|\langle
S_{p,i}\rangle+\langle S_{n,i}\rangle\right|^{2}\,\sigma_{p}^{\rm
SD}\leavevmode\nobreak\ ,$ (6)

where $\tilde{\mu}_{A_{i}}$ and $\tilde{\mu}_{p}$ are the reduced masses of
the DM-nucleus $i$ and DM-proton systems, $J_{i}$ is total angular momentum of
nucleus $i$, $\sigma_{p}^{\rm SI,SD}$ is the SI/SD DM-proton scattering cross
section, and $\langle S_{p,i}\rangle$ and $\langle S_{n,i}\rangle$ are the
expectation values of the spins of protons and neutrons averaged over all
nucleons [85, 86, 87, 88, 89, 90]. Multiple mediators could exist in some
models that couple DM to the visible sector. As long as the total scattering
cross section remains (approximately) independent of Mandelstam variables in
the relevant momentum transfer regime, the discussion presented below applies.

The differential rate at which a DM particle with velocity $w$ scatters off a
target $i$ with velocity $u$ and relative angle $\theta_{\chi}$, in the
laboratory frame, to a final velocity $v$ is given by [3]

$\mathcal{R}_{i}^{\pm}(w\to v)=\int
n_{i}(r)\,\frac{\text{d}\sigma_{i}}{\text{d}v}\,|\bm{w}-\bm{u}|\,f_{i}(\bm{u},r)\,\text{d}^{3}\bm{u}=\frac{2}{\sqrt{\pi}}\,\frac{n_{i}(r)}{u_{i}^{3}(r)}\,\int_{0}^{\infty}\text{d}u\,u^{2}\,\int_{-1}^{1}\text{d}\cos\theta_{\chi}\,\frac{\text{d}\sigma_{i}}{\text{d}v}\,|\bm{w}-\bm{u}|\,e^{-u^{2}/u_{i}^{2}(r)}\leavevmode\nobreak\
.$ (7)

The superindex ± indicates $v\gtrless w$, $f_{i}(\bm{u},r)$ is the velocity
distribution of the target particles $i$ (assumed to be Maxwell-Boltzmann),
with density $n_{i}(r)$ and temperature $T(r)$, and
$u_{i}(r)\equiv\sqrt{2\,T(r)/m_{i}}$ is the most probable speed of the target
particles at position $r$. The relative velocity between the DM and the target
particles is given by
$|\bm{w}-\bm{u}|=\sqrt{w^{2}+u^{2}-2\,w\,u\,\cos{\theta_{\chi}}}$, where $w$
is related to $u_{\chi}$ by $w^{2}(r)=u_{\chi}^{2}+v_{e}^{2}(r)$, with
$v_{e}(r)$ the escape velocity at a radial distance $r$ from the center of the
object. Analytical expressions for these rates when the nuclear form factor is
not included or when the targets temperature is neglected (for the capture
process, i.e., for $\mathcal{R}^{-}(w\to v)$) were obtained long ago [3, 4].

###  II.2 DM capture by celestial bodies

With all the above ingredients, the capture rate of DM particles by celestial
bodies can be generically written as [4, 67]

$\mathcal{C}=\sum_{i}\int_{0}^{R}s_{\rm
cap}(r)\,4\pi\,r^{2}\,\text{d}r\int_{0}^{\infty}\text{d}u_{\chi}\,\left(\frac{\rho_{\chi}}{m_{\chi}}\right)\,\frac{f_{v_{\rm
cb}}(u_{\chi})}{u_{\chi}}\,w(r)\int_{0}^{v_{e}(r)}\mathcal{R}_{i}^{-}(w\to
v)\,\text{d}v\leavevmode\nobreak\ ,$ (8)

where the sum runs over all possible targets. We have also included the
suppression factor $s_{\rm cap}(r)$, which ignores multiple collisions, but
allows for a smooth transition between the optically thin (small cross
sections) and optically thick regimes (large cross sections) [67]. Not only
$\mathcal{R}^{-}_{i}(w\to v)$, but the total capture rate per unit volume can
also be analytically computed when either the targets temperature is neglected
or at finite temperature without including a non-trivial nuclear form factor
[4], assuming $s_{\rm cap}(r)=1$.

Thermal effects have a small impact on capture of DM particles by the nuclei
of celestial objects (i.e., $m_{\chi}v_{e}^{2}/2\gg T$), although they would
be more relevant in the interior of the smallest objects discussed in this
work. Moreover, effects from the lack of coherence are only important when the
inverse size of the target is smaller than the DM escape energy in the galaxy
and for heavy DM particles.333For instance, for the Earth they are negligible
for oxygen and only modest for iron, at mass resonance, whereas for the Sun
and $m_{\chi}\sim 15$ GeV, they suppress capture by iron, but not
significantly by oxygen or helium [4]. Larger celestial bodies have higher
escape velocities, so the suppression would affect more to lighter elements.
Nevertheless, in this work we are interested in the minimum mass of DM
particles that can be captured, for which these effects are less important
than for higher masses. In any case, they are included below in our numerical
computations and just neglected in this section for illustrative purposes.
Therefore, for the illustrative purposes of this section, we consider the
zero-temperature limit and do not include the nuclear form factor. In this
limit, the capture rate in the optically thin regime ($s_{\rm cap}(r)=1$)
reads [4]

$\mathcal{C}_{\rm weak}=\left(\frac{\rho_{\chi}}{m_{\chi}}\right)\langle
v\rangle_{0}\sum_{i}N_{i}\,\sigma_{i}\,\bigg{\langle}\frac{\hat{\phi}}{\langle\hat{\phi}\rangle_{i}}\,\left(1-\frac{1-e^{-\mathcal{B}_{i}^{2}}}{\mathcal{B}_{i}^{2}}\right)\,\xi_{\eta}(\mathcal{B}_{i})\bigg{\rangle}_{i}\,\left(\frac{3}{2}\frac{v_{e}^{2}(R)}{v_{d}^{2}}\langle\hat{\phi}\rangle_{i}\right)\leavevmode\nobreak\
,$ (9)

where $\langle v\rangle_{0}=\sqrt{8/(3\pi)}\,v_{d}$ is the average speed in
the galactic rest frame, which depends on the position in the halo, and

$\mathcal{B}_{i}^{2}(r)=\frac{3}{2}\frac{v_{e}^{2}(r)}{v_{d}^{2}}\frac{\mu_{i}}{\mu_{-,i}^{2}}\hskip
28.45274pt;\hskip
8.53581pt\hat{\phi}(r)=\frac{v_{e}^{2}(r)}{v_{e}^{2}(R)}\hskip
28.45274pt;\hskip
8.53581pt\langle\hat{\phi}\rangle_{i}=\frac{\int_{0}^{R}\hat{\phi}(r)\,n_{i}(r)\,4\pi\,r^{2}\,\text{d}r}{N_{i}}\leavevmode\nobreak\
,$ (10)

with $\mu_{i}=m_{\chi}/m_{i}$ and $\mu_{-,i}=(\mu_{i}-1)/2$. The term
$\xi_{\eta}(\mathcal{B}_{i})$ represents a suppression factor that accounts
for the motion of the celestial body with respect to the halo frame
($\xi_{0}(\mathcal{B}_{i})=1$), with $\eta^{2}=3\,v_{\rm
cb}^{2}/(2\,v_{d}^{2})$ [4]. This is a consequence of a higher DM kinetic
energy in the object’s frame. For $\eta=1$, $\xi_{1}(\mathcal{B}_{i})$
monotonically takes values in the interval $(0.37-0.75)$, with the extremes
reached for $\mathcal{B}_{i}^{2}(r)\ll 1$ and $\mathcal{B}_{i}^{2}(r)\gg 1$,
respectively. Note also that $\xi_{\eta}(\mathcal{B}_{i})$ includes a
dependence on $\mathcal{B}_{i}$, but it is a ratio of two sums over all
elements, so $\xi_{\eta}$ does not itself run over $i$, although it does
depend on the position within the celestial body.

In the thick regime (large cross section), the geometric limit for the capture
rate is usually considered [78, 91, 66]. This is an upper limit obtained from
purely geometrical arguments and assumes that both, the probability of
interaction and the probability of DM capture are one. This is, however, a
non-physical limit, as even if the interaction probability is one, the
probability of DM particles being captured (ending up with speeds below the
escape velocity) is dictated by kinematics and is always smaller than one.
This is corrected by the so-called saturation limit [67]. In the limit of high
escape velocities, $v_{e}>v_{d}$, the geometric limit only overestimates the
saturation value by ${\cal O}(10\%)$ (for masses around the DM evaporation
mass), in the case of the Sun, so it represents a reasonable
approximation.444In the limit $v_{e}\to\infty$, the geometric and saturation
limits coincide, as could be expected. Nevertheless, when $v_{e}<v_{d}$, the
geometric limit can grossly overestimate the maximum value for the capture
rate (except for DM masses closely matching the targets mass). The
overestimation scales as $v_{d}^{4}/v_{e}^{4}$ and therefore can be of several
orders of magnitude, being more important for smaller objects.555This factor
has been incorrectly neglected in some recent works that considered the
maximum value of the capture rate [92, 32, 93, 68]. Therefore, some of the
presented results must be re-scaled, changing significantly the conclusions in
some cases. Although, the impact on the DM evaporation mass is generically
mild, we use the more accurate saturation limit in what follows. In analogy to
the capture rate for weak cross sections, the saturation limit can be written
as

$\mathcal{C}_{\rm
sat}=\frac{\pi\,R^{2}}{\sum_{i}N_{i}\,\sigma_{i}}~\left(\frac{\rho_{\chi}}{m_{\chi}}\right)\langle
v\rangle_{0}\left[\sum_{i}N_{i}\,\sigma_{i}\,\left(1-\frac{1-e^{-\mathcal{B}_{i}^{2}(R)}}{\mathcal{B}_{i}^{2}(R)}\right)\right]\,\xi_{1}(\mathcal{B}_{i}(R))\,\left(\frac{3}{2}\frac{v_{e}^{2}(R)}{v_{d}^{2}}\right)\leavevmode\nobreak\
.$ (11)

In this section, for the sake of the discussion we consider this limit with a
single element.

###  II.3 DM annihilations in celestial bodies

After DM particles get captured, further scatterings with the target elements,
which are assumed to be in local thermodynamic equilibrium, would
approximately thermalize them at a temperature $T_{\chi}(r)$ and attain a
velocity distribution that can be approximated as Maxwell-Boltzmann.666For not
too small cross sections, the thermalization time is typically much shorter
than any relevant scale in this problem [11, 12, 13, 14, 15, 16] (otherwise,
the entire calculation would not be valid), such that the instantaneous
approximation is adequate. Note also that the assumptions of a uniform and
locally isotropic Maxwell-Boltzmann distribution for DM particles do not
exactly hold in a realistic situation [3, 94, 65, 12], although it is a
reasonable approximation. Moreover, we assume the distribution to be locally
truncated at the escape velocity, $v_{e}(r)$ [77]. In the case of weak cross
sections (optically thin regime), the DM radial distribution is approximately
isothermal [6, 7, 2],

$n_{\chi,{\rm
iso}}(r,t)=N_{\chi}(t)\,\frac{e^{-m_{\chi}\phi(r)/T_{\chi}}}{\int_{0}^{R}e^{-m_{\chi}\phi(r)/T_{\chi}}\,4\pi
r^{2}\,\text{d}r}\leavevmode\nobreak\ ,$ (12)

where
$\phi(r)=\int_{0}^{r}GM(r^{\prime})/{r^{\prime}}^{2}\,\text{d}r^{\prime}$,
with $G$ the gravitational constant. For all DM masses, the DM temperature is
a fraction of the central temperature of the object, being higher for heavier
DM particles, which are more centrally concentrated. We compute the DM
temperature, $T_{\chi}$, following Ref. [66]. In the case of large cross
sections (optically thick regime), DM particles would thermalize locally with
the medium and thus, $T_{\chi}(r)=T(r)$, with a radial distribution that can
be approximated as [10, 77]

$n_{\chi,{\rm LTE}}(r,t)=n_{\chi,{\rm
LTE,0}}(t)\,\left(\frac{T(r)}{T(0)}\right)^{3/2}\,{\rm
exp}\left(-\int_{0}^{r}\frac{\alpha(r^{\prime})\frac{\text{d}T(r^{\prime},t)}{\text{d}r^{\prime}}+m_{\chi}\frac{\text{d}\phi(r^{\prime})}{\text{d}r^{\prime}}}{T(r^{\prime})}\,\text{d}r^{\prime}\right)\leavevmode\nobreak\
.$ (13)

Here, $\alpha$ is the thermal diffusivity [77, 70] and $n_{\chi,{\rm
LTE},0}(t)$ is set by the normalization $\int_{0}^{R}n_{\chi,{\rm
LTE}}(r)\,4\pi r^{2}\,\text{d}r=N_{\chi}(t)$. The transition between the thin
and thick regimes can be described in terms of the Knudsen number,
$\mathcal{K}$, which is defined as the ratio of the mean free path to the
scale radius. Here, we follow Refs. [78, 46, 66] to interpolate between the
two regimes, using $\mathcal{K}_{0}=0.4$ as the pivot point [77] for all
cases. Note, however, that $\mathcal{K}$ depends on the properties of the
capturing body.

The radial distribution of thermalized DM particles clusters around the center
of the capturing body, being more centrally concentrated in the case of
heavier DM particles. Indeed, for all celestial bodies considered here, for
the DM evaporation mass (i.e., for the lightest DM particles that can get
efficiently captured), the scale radius of the distribution is $r_{s}\lesssim
0.2\,R$. For larger objects, like stars, $r_{s}\lesssim 0.1\,R$. Therefore,
for the case of $s$–wave DM annihilations, which is considered throughout this
work, the annihilation rate is given by

$\mathcal{A}=\langle\sigma_{A}v_{\chi\chi}\rangle\,\frac{\int_{0}^{R}n_{\chi}^{2}(r,t)\,4\pi\,r^{2}\,\text{d}r}{\left(\int_{0}^{R}\,n_{\chi}(r,t)\,4\pi\,r^{2}\,\text{d}r\right)^{2}}\simeq\frac{\langle\sigma_{A}v_{\chi\chi}\rangle}{V_{s}}\leavevmode\nobreak\
,$ (14)

where $\langle\sigma_{A}v_{\chi\chi}\rangle$ is the velocity-averaged DM
annihilation cross section times the relative velocity of two DM particles,
and throughout this paper we use the canonical value
$\langle\sigma_{A}v_{\chi\chi}\rangle=3\times 10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}$ as our default value. In the second equality we
have simply substituted the effective volume of integration by the volume at
the scale radius, $V_{s}=4/3\,\pi\,r_{s}^{3}$, where most of the DM particles
are concentrated. One can think of a region of approximately constant density
and temperature (those at the core), which is a reasonable approximation
within this small volume. In this section, we take $r_{s}=0.1\,R$ for
illustrative purposes, although a more precise value affects our estimates in
a negligible way.

###  II.4 DM evaporation off celestial bodies

In addition to the disappearance of captured DM particles due to
annihilations, scatterings off targets of the medium could boost these
particles to speeds above the local escape velocity, $v_{e}(r)$, such that
they become gravitationally unbound. This process is referred to as DM
evaporation, which is a finite-temperature process that depends on the
temperature of the thermal bath in the celestial body that can transfer
kinetic energy to DM particles. In the limit of zero temperature, there is no
DM evaporation, as target particles cannot impinge any extra energy to
thermalized DM particles. At finite temperature, energy conservation results
in a higher probability for light DM particles to end up with higher final
speeds, after scattering off thermal targets. Thus, for small DM masses, this
process is very efficient and sets a minimum mass of DM particles that can
remain trapped in a celestial body. Lighter DM particles would get kicked out
as they get captured. The evaporation rate is given by [3]

$\mathcal{E}=\sum_{i}\int_{0}^{R}s_{\rm
evap}(r)\,n_{\chi}(r,t)\,4\pi\,r^{2}\,\text{d}r\,\int_{0}^{v_{e}(r)}f_{\chi}(\bm{w},r)\,4\pi\,w^{2}\,\text{d}w\,\int_{v_{e}(r)}^{\infty}\mathcal{R}_{i}^{+}(w\rightarrow
v)\,\text{d}v\leavevmode\nobreak\ .$ (15)

where $f_{\chi}(\bm{w},r)$ is the thermal velocity distribution of DM
particles and $s_{\rm evap}(r)$ is a suppression factor that accounts for the
fraction of DM particles that, even with a speed higher than the escape
velocity, would actually escape due to further scatterings on their way out of
the celestial body [70]. We follow Refs. [91, 66] to implement this
suppression factor. As mentioned above, in the thin regime (long mean free
path, $s_{\rm evap}\sim 1$), the DM distribution can be approximated as
isothermal, with a temperature close to the central temperature of the
celestial body. Thus, evaporation takes places mainly in a small region around
the core. For large cross sections (thick regime, $s_{\rm evap}\ll 1$) the
shell that contributes most to the evaporation rate moves towards the surface
of the object, although it never reaches the last scattering surface (i.e.,
that for which the optical depth is equal to one) [70].

In the limit of the thin regime, $s_{\rm evap}(r)=1$, an analytical solution
for the evaporation rate per unit volume exists [3]. For $m_{\chi}=m_{i}$, in
the limit $E_{e}=m_{\chi}\,v_{e}^{2}/2\gg T$, the solution is rather simple.
Given that the DM evaporation mass is typically of the same order of the
targets mass and the DM escape energy is larger than the thermal energy, we
can consider this solution to illustrate the main features of the evaporation
rate, which can be approximated as [2, 3]

$\mathcal{E}\simeq\sum_{i}\left[\frac{1}{V_{s}}\frac{2}{\sqrt{\pi}}\left(\frac{2\,T_{\chi}}{m_{\chi}}\right)^{1/2}\,\left(\frac{E_{c}}{T_{\chi}}\right)\,e^{-E_{c}/T_{\chi}}\,\right]\,N_{i}(r_{0.95})\,\sigma_{i}\leavevmode\nobreak\
,$ (16)

where $E_{c}=m_{\chi}v_{e}^{2}(r=0)/2$ is the escape energy of DM particles at
the core of the celestial body and $N_{i}(r_{0.95})$ is the number of targets
$i$ within a radius $r_{0.95}$ such that $T(r_{0.95})=0.95\,T_{\chi}$, which
typically represents a small fraction of the total mass of the object. For the
sake of illustration, in the next subsection we simply consider
$N_{i}(r_{0.95})=0.1\,M/m_{i}$, which is approximately correct for the Sun
[3]. The chosen value affects very little the calculation of the DM
evaporation mass, because the evaporation rate depends exponentially on the DM
mass. Moreover, note that Eq. (16) is obtained for $\mu_{i}=1$
($m_{\chi}=m_{i}$), but the evaporation rate depends very weakly on this
ratio, except if $m_{i}\gg m_{\chi}$ [3]. Yet, in general, evaporation becomes
efficient for $m_{\chi}$ not very different from $m_{i}$.

For large cross sections (thick regime), an analogous approximate expression
for the evaporation rate can be written [70], including a suppression factor
which accounts for the short mean free path of kicked DM particles. The
expression, though, depends more critically on the temperature and density
profiles of the celestial body than in the thin regime, and we do not
explicitly consider it in this section, although we do compute it in detail
for our numerical results. In any case, the discussion is qualitatively
analogous.

###  II.5 DM evaporation mass: a tale of two tails

With all the relevant quantities briefly discussed and the simple
approximations established, we turn to the estimation of the DM evaporation
mass for a generic object. This is well known for the Sun and we simply
restate and explain the results in the context of a generic celestial body
with mass $M$, radius $R$ and core temperature $T_{c}$, whose properties
remain approximately constant in time. We stress that the inputs used in this
section are just presented for illustrative purposes, but we use the complete
expressions to obtain our numerical results in Section IV.

We consider the saturation value for the capture rate, Eq. (11), and the
annihilation and evaporation rates given in Eqs. (14) and (16), respectively.
Before equilibration ($t\ll\tau_{\rm eq}$), the DM evaporation mass grows with
time ($\mathcal{E}\,t\simeq\ln(11)$) up to its maximum value at equilibrium.
At this point, the DM evaporation mass is independent of time (assuming all
the properties of the capturing object remain the same) and, from Eq. (3), can
be defined as [66]

$\mathcal{E}(m_{\rm evap})\,\tau_{\rm eq}(m_{\rm
evap})\simeq\frac{1}{\sqrt{0.11}}\leavevmode\nobreak\ .$ (17)

For $3\,v_{e}^{2}\,\mu_{i}\gg 2\,v_{d}^{2}\,\mu_{-,i}^{2}$
($\mathcal{B}_{i}^{2}\gg 1$), representative of stars like the Sun (except for
heavy DM, which is not the focus of this paper), and for the geometric cross
section, $\sum N_{i}\,\sigma_{i}^{\rm geom}=\pi\,R^{2}$, which approximately
sets the maximum capture rate (the saturation value), the equilibration time
approximately scales as (using $r_{s}=0.1\,R$)

$\tau_{\rm eq}\simeq 7\times 10^{12}\leavevmode\nobreak\ {\rm
s}\,\left(\frac{M_{\odot}}{M}\right)^{1/2}\,\left(\frac{R}{R_{\odot}}\right)\,\left(\frac{0.4\leavevmode\nobreak\
\textrm{GeV/cm}^{3}}{\rho_{\chi}}\right)^{1/2}\,\left(\frac{m_{\chi}}{\textrm{GeV}}\right)^{1/2}\,\left(\frac{v_{d}}{270\leavevmode\nobreak\
\textrm{km/s}}\right)^{1/2}\,\left(\frac{3\times 10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}}{\langle\sigma_{A}v_{\chi\chi}\rangle}\right)^{1/2}\leavevmode\nobreak\
.$ (18)

The evaporation rate approximately scales as (using $\sigma_{\rm
geom}\,N(r_{0.95})=0.1\,\pi\,R^{2}$)

$\mathcal{E}\simeq 0.06\leavevmode\nobreak\
\textrm{s}^{-1}\,\left(\frac{R_{\odot}}{R}\right)\,\left(\frac{T_{\chi}}{1.5\times
10^{7}\leavevmode\nobreak\
\textrm{K}}\right)^{1/2}\,\left(\frac{\textrm{GeV}}{m_{\chi}}\right)^{1/2}\,\left(\frac{E_{c}}{T_{\chi}}\right)\,e^{-E_{c}/T_{\chi}}\leavevmode\nobreak\
,$ (19)

where we have kept the dependence on $E_{c}/T_{\chi}$ explicit. For celestial
bodies like the Sun (and in general for other objects), the prefactor is very
large, so it is obvious that $E_{c}/T_{\chi}\gg 1$ in order to suppress the
evaporation rate at the level of $\tau_{\rm eq}^{-1}$. Putting these two
quantities together, the equation for the DM evaporation mass reads

$\left(\frac{E_{c}}{T_{\chi}}\right)e^{-E_{c}/T_{\chi}}\simeq 7\times
10^{-12}\,\left(\frac{M}{M_{\odot}}\right)^{1/2}\left(\frac{1.5\times
10^{7}\leavevmode\nobreak\
\textrm{K}}{T_{\chi}}\right)^{1/2}\left(\frac{\rho_{\chi}}{0.4\leavevmode\nobreak\
\textrm{GeV/cm}^{3}}\right)^{1/2}\left(\frac{270\leavevmode\nobreak\
\textrm{km/s}}{v_{d}}\right)^{1/2}\left(\frac{\langle\sigma_{A}v_{\chi\chi}\rangle}{3\times
10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}}\right)^{1/2}\leavevmode\nobreak\ .$ (20)

For the representative (solar) values used in the above equation, the solution
is $E_{c}/T_{\chi}\simeq 29$, which results in a DM evaporation mass for the
Sun, $m_{\rm evap}\simeq 3.2$ GeV, using $v_{e}^{2}(r=0)=5\,v_{e}^{2}(r=R)$
and $T_{\chi}=0.9\,T_{c}$. This is known in the literature for over three
decades [6, 21, 2, 3]. Already the authors of Ref. [2] explicitly wrote that
the relevant value to compute the DM evaporation mass for the Sun is
$E_{c}/T_{\chi}\simeq 30$. Furthermore, note that for smaller scattering cross
sections, the right-hand side scales as $\sigma^{-1/2}$, so the DM evaporation
mass is slightly smaller. For larger cross sections, there is an extra
(exponential) suppression term in the evaporation rate (left-hand side in the
above equation) and the DM evaporation mass is also smaller, and more
pronouncedly than for smaller cross sections. Notice also that the details of
the density and temperature profiles are embedded in $v_{e}^{2}(0)$ and
$T_{\chi}$, which cannot vary much given the mass, radius and core temperature
of the celestial body, and thus, affect little the value of the DM evaporation
mass. In particular, variations on $v_{e}^{2}(0)$ are expected to be $\lesssim
10\%$, which also applies to the smaller objects we discuss next.

In the opposite limit, representative of the Earth (and other solar system
planets), $\mathcal{B}_{i}^{2}\ll 1$. Therefore, in addition to the
$v_{e}^{2}/v_{d}^{2}$ suppression, there is an extra suppression factor in the
capture rate,
$\mathcal{B}_{i}^{2}/2=(3\,v_{e}^{2}\,\mu_{i})/(4\,v_{d}^{2}\,\mu_{-,i}^{2})$,
except in a very narrow mass range, where the mass of the targets and of the
DM particles closely match (in that case, the results for $\mathcal{B}_{i}\gg
1$ are reproduced). This implies a longer equilibration time and therefore, a
higher DM evaporation mass. In this limit, the equations for the equilibration
time, the evaporation rate and the DM evaporation mass read

$\displaystyle\tau_{\rm eq}$ $\displaystyle\simeq$ $\displaystyle 1.5\times
10^{15}\leavevmode\nobreak\ {\rm
s}\,\left(\frac{\mu_{-}^{2}}{3\,\mu}\right)^{1/2}\left(\frac{M_{\oplus}}{M}\right)\left(\frac{R}{R_{\oplus}}\right)^{3/2}$
$\displaystyle\hskip 42.67912pt\times\left(\frac{0.4\leavevmode\nobreak\
\textrm{GeV/cm}^{3}}{\rho_{\chi}}\right)^{1/2}\left(\frac{m_{\chi}}{\textrm{GeV}}\right)^{1/2}\left(\frac{v_{d}}{270\leavevmode\nobreak\
\textrm{km/s}}\right)^{3/2}\left(\frac{3\times 10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}}{\langle\sigma_{A}v_{\chi\chi}\rangle}\right)^{1/2}\leavevmode\nobreak\
,$ $\displaystyle\mathcal{E}$ $\displaystyle\simeq$ $\displaystyle
0.1\leavevmode\nobreak\
\textrm{s}^{-1}\,\left(\frac{R_{\oplus}}{R}\right)\left(\frac{T_{\chi}}{6000\leavevmode\nobreak\
\textrm{K}}\right)^{1/2}\left(\frac{m_{\chi}}{\textrm{GeV}}\right)^{-1/2}\left(\frac{E_{c}}{T_{\chi}}\right)\,e^{-E_{c}/T_{\chi}}\leavevmode\nobreak\
,$ (22)
$\displaystyle\left(\frac{E_{c}}{T_{\chi}}\right)\,e^{-E_{c}/T_{\chi}}$
$\displaystyle\simeq$ $\displaystyle 2\times
10^{-14}\,\left(\frac{3\,\mu}{\mu_{-}^{2}}\right)^{1/2}\left(\frac{M}{M_{\oplus}}\right)\left(\frac{R_{\oplus}}{R}\right)^{1/2}\left(\frac{6000\leavevmode\nobreak\
\textrm{K}}{T_{\chi}}\right)^{1/2}$ $\displaystyle\hskip
34.14322pt\times\left(\frac{\rho_{\chi}}{0.4\leavevmode\nobreak\
\textrm{GeV/cm}^{3}}\right)^{1/2}\left(\frac{270\leavevmode\nobreak\
\textrm{km/s}}{v_{d}}\right)^{3/2}\left(\frac{\langle\sigma_{A}v_{\chi\chi}\rangle}{3\times
10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}}\right)^{1/2}\leavevmode\nobreak\ ,$

which results in $E_{c}/T_{\chi}\simeq 34$ for the Earth, and in a DM
evaporation mass $m_{\rm evap}\simeq 13$ GeV, where we have used
$v_{e}^{2}(r=0)=1.9\,v_{e}^{2}(r=R)$, $T_{\chi}=T_{c}$ and $\mu=1/3$. Note
that $E_{c}/T_{\chi}$ is similar to the value obtained in the limit
$\mathcal{B}_{i}^{2}\gg 1$ and $m_{\rm evap}$ is very close to the value
quoted in Ref. [23]. Indeed, the fact that the DM evaporation mass is
approximately given by $E_{c}/T_{\chi}\sim 30$ is a robust result that applies
to all round objects in hydrostatic equilibrium. The same considerations
regarding the dependence on the scattering cross section can be drawn for
$\mathcal{B}_{i}^{2}\ll 1$. In this limit, however, the equilibration time is
longer. Yet, this is approximately the shortest possible equilibration time
(we are using a capture rate close to maximum), so cross sections smaller than
the geometric value could result in equilibration times longer than the age of
the object. In those cases, the approximate equation to be used is
$\mathcal{E}\,t\simeq\ln(11)$ and the DM evaporation mass grows with time.
Indeed, it is unlikely the smallest objects we consider reach equilibrium,
unless the DM mass closely matches the mass of some of the targets or they are
close to the center of their host halo (high $\rho_{\chi}$ and low $v_{d}$).

All in all, the exponential dependence of the evaporation rate, which sets the
DM evaporation mass along the exponential tail, is really a tale of two tails
[21]. On one hand, a DM particle in the high-velocity tail of its distribution
could scatter off a target with typical thermal speed, and be promoted to a
speed higher than the escape velocity. On the other hand, a DM particle with
typical thermal speed may be kicked off the celestial object due to the
scattering with a target in the high-velocity tail of its distribution. The
first process is the most important one [21].

Figure 1: The two tails of the evaporation rate. The probability for a DM
particle with speed $w$ to scatter off a target in the medium, with
temperature $T$, and gain enough energy to escape from the gravitational
potential of the capturing object, $\Omega_{v_{e}}^{+}(w)$ (solid lines), and
the velocity distribution of DM particles, with temperature $T_{\chi}=0.9\,T$,
$f_{\chi}(w)$ (dashed lines), both in arbitrary units. We show these
distributions for $E_{c}/T_{\chi}=10$ (blue lines) and $E_{c}/T_{\chi}=30$
(green lines), using $\mu=1$.

One way to visualize the exponential suppression of the evaporation rate is by
plotting the thermal velocity distribution of DM particles (with temperature
$T_{\chi}$), which scales as
$\sim\left(E_{e}/T_{\chi}\right)^{3/2}(w/v_{e})^{2}\,\exp[-(E_{e}/T_{\chi})\,(w/v_{e})^{2}]$,
and the probability for a DM particle with speed $w$ to interact with a target
particle (with a thermal distribution of temperature $T(r)$) and gain enough
energy to escape. The latter probability is proportional to
$\Omega_{v_{e}}^{+}(w)=\int_{v_{e}}^{\infty}\mathcal{R}^{+}(w\to
v)\,\text{d}v$ [3], and approximately scales as
$\sim\exp[-(E_{e}/T)(1-(w/v_{e})^{2})]$. We evaluate $E_{e}$ at the core,
$E_{c}$, because in the thin regime DM evaporation mostly occurs close to the
center of the celestial body. In Fig. 1 we show both distributions as a
function of $w/v_{e}$, for two values of $E_{c}/T_{\chi}$. As can be clearly
seen, the smaller $E_{c}/T_{\chi}$, the larger the overlap of the two
distributions, or in other words, the higher the evaporation rate. For small
$E_{c}/T_{\chi}$, more DM particles have speeds close to the escape velocity
(which after evaporation get repopulated by the thermalization process), so
they need less energy to evaporate. Moreover, the smaller $E_{c}/T_{\chi}$,
the higher the probability for DM particles to end up with speeds higher than
the escape velocity. Therefore, the overall probability for this to happen is
higher, and in relative terms, the ratio of probabilities for two values of
$E_{c}/T_{\chi}$ approximately scales as the exponential of the difference
between the two values of $E_{c}/T_{\chi}$. These two effects go in the same
direction and result in a huge evaporation rate when a significant overlap
occurs. As a consequence, in order to suppress sufficiently the evaporation
rate to maintain an equilibrium population of DM particles, the two
distributions should only overlap far in their exponential tails and
$E_{c}/T_{\chi}\sim 30$ is generically required. We stress that this is a
robust result, which not only applies to the Sun, but to all round celestial
bodies in hydrostatic equilibrium, and provides the correct result, for
constant scattering cross sections at the geometric value, within $\lesssim
30\%$ accuracy.

Finally, note that some recent works, which also assumed constant scattering
cross sections, did not follow these arguments and incorrectly estimated the
DM evaporation mass [30, 95, 32, 76, 68]. In Ref. [30], the DM evaporation
mass was defined as the DM mass for which the thermal radius (obtained from
the virial theorem) is equal to the radius of the capturing body, which is
similar to setting $E_{e}/T_{\chi}\sim 1$. This implies a DM evaporation mass
for the Earth of $\sim 100$ MeV, which is not correct. For larger masses with
smaller thermal radius than the radius of the object, the evaporation rate is
already very high and, as explained above, those DM particles would also
evaporate. Under the very same assumptions, the correct value of the DM
evaporation mass for the Earth is $m_{\rm evap}\sim 12$ GeV [18, 19, 23, 31].
On the other hand, Ref. [32] proposes the observation of exoplanets and brown
dwarfs to search for effects of captured DM particles with sub-GeV mass.
Nevertheless, these authors set the condition for DM particles to remain
trapped to $E_{e}/T_{\chi}>1$,777Actually, the authors of Ref. [32] use the
condition $3\,T(r)/2<GM(r)\,m_{\chi}/(2\,r)$, which, as argued above, is not
correct. Additionally, the right-hand side is claimed to be the gravitational
potential at position $r$, but notice that this not correct. Moreover, it is
not clear how the minimum bounded DM masses which are quoted in that paper
were obtained. neglecting the crucial exponential tail of the evaporation
rate. Again, this implies a gross underestimation of the DM evaporation mass
and explains why Ref. [32] (and also Refs. [76, 68]) incorrectly claimed DM
evaporation masses well in the sub-GeV regime. A DM evaporation mass for the
Earth of $\sim 100$ MeV was also incorrectly estimated, as well as for other
planets and brown dwarfs, whose DM evaporation mass was also greatly
underestimated to be as low as $m_{\rm evap}\sim 4.5$ MeV.888In addition to
neglecting the exponential tail of the evaporation rate, Refs. [32, 76] also
consider core temperatures for brown dwarfs which are about a factor of six
smaller than what we do (see below). Overall, this amounts to about a factor
of $\sim 200$ underestimation of the DM evaporation mass with respect to our
result. A more accurate estimation for the most massive brown dwarfs was
$m_{\rm evap}\sim{\cal O}(1)$ GeV [33]. Below, this is explicitly discussed in
more detail and contrasted with our results.

###  II.6 DM evaporation off compact bodies: white dwarfs and neutron stars

As we illustrate below, for white dwarfs and neutron stars, the DM evaporation
mass is much smaller than the targets masses, $\mu_{i}\ll 1$. Therefore, the
approximation we considered above to estimate the evaporation rate (i.e.,
$E_{e}/T\gg 1$) is not appropriate in these cases, as for these objects,
$\mu_{i}\,E_{e}/T\ll 1$. In addition, so far in this section, we have only
considered the case of non-degenerate targets. This is valid for all cases
discussed in this paper except for neutron stars (and for brown and white
dwarfs if the targets are electrons, which we do not consider in this paper).
In the case of neutron stars, the targets we study are degenerate neutrons and
the calculation of the capture and evaporation rates of DM particles (and also
of the DM thermalization time [11, 16]) needs to be modified to properly
include the effect of Pauli blocking. This is approximately accounted for by
adding a correction factor, which in the case of the evaporation rate is
$\zeta_{i}\simeq{\rm min}\\{2\,T_{c}/p_{{\rm F},i}\,,1\\}$ [71], where
$p_{{\rm F},i}$ is the Fermi momentum of target particles $i$. As will become
obvious in the next paragraphs, for the case of interactions with nucleons,
the correction factor in the capture rate is not relevant in the discussion of
the DM evaporation mass.

Including the impact of Pauli blocking and taking the limit
$\mu_{i}\,E_{e}/T\ll 1$ (with $E_{e}/T\gg 1$), the evaporation rate for equal
temperatures ($T_{c}=T_{\chi}$), using Ref. [3], can be written as

$\mathcal{E}\simeq\sum_{i}\,\frac{4}{3}\left(\frac{\mu_{i}}{\pi}\right)^{1/2}\,\left(\frac{E_{c}}{T_{\chi}}\right)^{1/2}\,\zeta_{i}\,\left[\frac{1}{V_{s}}\frac{2}{\sqrt{\pi}}\left(\frac{2\,T_{\chi}}{m_{\chi}}\right)^{1/2}\,\left(\frac{E_{c}}{T_{\chi}}\right)\,e^{-E_{c}/T_{\chi}}\,\right]\,N_{i}(r_{x})\,\sigma_{i}\leavevmode\nobreak\
.$ (24)

For non-relativistic and degenerate target particles, such as neutrons in
neutron stars, the Fermi momentum is $p_{\rm F}\simeq\sqrt{2\,m_{i}\,\mu_{\rm
F}}$, with $\mu_{\rm F}\simeq 350$ MeV [71]. For relativistic and fully
degenerate target particles, such as electrons in white dwarfs and in neutron
stars, $p_{\rm F}\simeq\mu_{\rm F}$. As done above for other celestial bodies,
for the following estimates we take
$\sum_{i}N_{i}(r_{x})\,\sigma_{i}=0.1\,\pi\,R^{2}$, although the DM
evaporation mass is very little sensitive to its precise value.

Now, we establish a distinction between white dwarfs and neutron stars. For
the former, we just proceed as in the preceding subsection, but using Eq. (24)
for the DM evaporation rate. For the latter, however, $m_{\rm evap}\lesssim
100$ keV, so the annihilation cross section cannot be taken as the canonical
value (otherwise, the DM contribution would overclose the Universe!). Thus,
for neutron stars, we assume a very small annihilation cross section, so that
equilibration between capture and annihilation does not take place. Using Eq.
(3) with $\mathcal{A}=0$, the (time-dependent) DM evaporation mass is defined
by $\mathcal{E}(m_{\rm evap})\,t_{\rm NS}=\ln(11)$, where $t_{\rm NS}$ is the
age of the neutron star. This is effectively equivalent to asymmetric DM
scenarios.

Therefore, for white dwarfs, and considering only interactions with nuclei, as
done throughout this work (i.e., $\zeta_{i}=1$), the DM evaporation mass can
be obtained by solving

$\left(\frac{E_{c}}{T_{\chi}}\right)^{3/2}e^{-E_{c}/T_{\chi}}\simeq 6\times
10^{-11}\,\left(\frac{M}{M_{\odot}}\right)^{1/2}\left(\frac{4\times
10^{5}\leavevmode\nobreak\
\textrm{K}}{T_{\chi}}\right)^{1/2}\left(\frac{\rho_{\chi}}{0.4\leavevmode\nobreak\
\textrm{GeV/cm}^{3}}\right)^{1/2}\left(\frac{270\leavevmode\nobreak\
\textrm{km/s}}{v_{d}}\right)^{1/2}\left(\frac{\langle\sigma_{A}v_{\chi\chi}\rangle}{3\times
10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}}\right)^{1/2}\leavevmode\nobreak\ .$ (25)

For a white dwarf with $M=M_{\odot}$ ($R\simeq 0.9\leavevmode\nobreak\
R_{\oplus}$) and $T_{c}=T_{\chi}=4\times 10^{5}$ K, this results in
$E_{c}/T_{\chi}\simeq 29$, which is in the same ballpark as the values for
other objects. In this case, the DM evaporation mass is $m_{\rm evap}\simeq
1.0$ MeV, where we have used $v_{e}^{2}(r=0)=3.9\,v_{e}^{2}(r=R)$.

In the case of neutron stars, the equation for the DM evaporation mass,
$\mathcal{E}(m_{\rm evap})\,t_{\rm NS}=\ln(11)$, can be written as

$\left(\frac{E_{c}}{T_{\chi}}\right)^{3/2}e^{-E_{c}/T_{\chi}}\simeq 3\times
10^{-12}\,\left(\frac{p_{\textrm{F}}}{0.8\leavevmode\nobreak\
\textrm{GeV}}\right)\,\left(\frac{10^{5}\leavevmode\nobreak\
\textrm{K}}{T_{\chi}}\right)^{3/2}\left(\frac{R}{11.5\leavevmode\nobreak\
\textrm{km}}\right)\left(\frac{4.5\leavevmode\nobreak\ \textrm{Gyr}}{t_{\rm
NS}}\right)\leavevmode\nobreak\ .$ (26)

For a neutron star with $M=1\leavevmode\nobreak\ M_{\odot}$, $R=11.5$ km,
$T_{c}=T_{\chi}=10^{5}$ K and $t_{\rm NS}=4.5$ Gyr, this results in
$E_{c}/T_{\chi}\simeq 32$, again very similar to the values for other objects.
Using $v_{e}^{2}(r=0)=1.5\,v_{e}^{2}(r=R)$, the DM evaporation mass is $m_{\rm
evap}\simeq 1.4$ keV.

Thus, the DM evaporation mass for white dwarfs and neutron stars is also
approximately given by $E_{c}/T_{\chi}\sim 30$. We stress that this is a
general and robust result which applies to all the objects we consider in this
work, that is, to all spherical celestial bodies in hydrostatic equilibrium.
In any case, these are approximate estimates, which allow us, nonetheless, to
obtain the DM evaporation mass with a precision of a few tens of percent. We
obtain our results in Section IV following Ref. [71] for the calculation of
the evaporation rate in neutron stars for DM scattering off non-relativistic
degenerate neutrons, and we follow the discussion above [66] for all other
objects. Although for the calculation of DM evaporation mass in neutron stars
a more accurate treatment must use relativistic kinematics [96, 72], the
required corrections do not significantly change the results obtained here.

##  III Main properties of celestial bodies

Figure 2: Radius of planetary bodies, brown dwarfs and main-sequence stars, as
a function of the mass of the object and using the mass–radius relation
described in the text. We also show a compilation of data corresponding to
planetary bodies [97] (solar system planets and satellites in blue), brown
dwarfs with measured radius [98, 99, 100, 101, 102, 103, 104, 105], low-mass
stars [106], intermediate-mass and massive stars [107]. Figure 3: Core
temperature of planetary bodies, brown dwarfs and main-sequence stars, as a
function of the mass of the object and using the mass–core temperature
relation described in the text. From detailed models and data, the estimated
core temperature for the Moon, Earth and Sun is also indicated. Figure 4:
Escape velocity at the surface of planetary bodies, brown dwarfs and main-
sequence stars, as a function of the mass of the object and using the
mass–radius relation reported in the text. The values for the Moon, the Earth
and the Sun are shown. We also indicate two values of the galactic dispersion
velocity $v_{d}$ (dashed lines), which are references for the local
neighborhood and near the galactic center. The DM capture rate by objects with
$v_{e}<v_{d}$ is suppressed by a factor proportional to $(v_{e}/v_{d})^{4}$,
but with resonance-like features for DM masses matching targets masses.

In this section, we describe the average properties of celestial round bodies
in hydrostatic equilibrium, spanning a wide mass range,
$[10^{-10}-10^{2}]\leavevmode\nobreak\ M_{\odot}$. In order to determine the
DM evaporation mass, the required main characteristics of the capturing
objects are the mass $M$, radius $R$, and their density $\rho(r)$ and
temperature $T(r)$ profiles. We stress that, for a given object’s mass, radius
and core temperature, the DM evaporation mass depends little on the shape of
the density (mainly via the ratio of the gravitational potential at the center
and at surface) or temperature profiles. The kinematics of elastic scattering
depends on the mass of the DM particles as well as that of the targets, so
another important factor is the composition of the material. For the sake of
simplifying the discussion, we include hydrogen (in terms of the $X_{\rm
H}\equiv X$ mass fraction), helium ($X_{\rm He}\equiv Y$ mass fraction) and as
representative of heavier elements we consider carbon, oxygen, water, silicate
perovskite (MgSiO3) and iron (the mass fraction of these heavier elements is
generically denoted by $Z$).

Given that some generic features of celestial bodies can be approximately
described in terms of polytropes, we first briefly introduce the properties of
objects with this kind of equation of state. We shall later use them as
ballpark models for some cases. Next, we provide an overview of the general
properties of planetary objects, brown dwarfs, main-sequence stars, post-main-
sequence evolutionary phases of stars, white dwarfs and neutron stars. The
process of DM capture is assumed not to modify their properties in a
significant way, so that we can still use results from standard modeling
without including DM effects.

Throughout this work, we consider the mass as the single variable that
determines the rest of the properties of celestial bodies, in an average way.
We provide parameterizations for the radius, core temperature, density and
temperature profiles and composition, as a function of the mass of the object.
All of them are based on actual data and modeling. We just impose continuity
at the transitions from one mass range to another.

The mass–radius and mass–core temperature relations reported in this section
are shown in Figs. 2 and 3, respectively. The escape velocity at surface is
shown as a function of the mass of celestial bodies in Fig. 4. For those
objects located far enough from the center of the host halo such that their
local dispersion velocity is higher than the escape velocity of the object,
the capture rate is suppressed (except at mass-matching, where resonance-like
features appear) and hence, the equilibration time is longer (see below) and
the DM evaporation mass is larger (assuming equilibrium is reached).

###  III.1 Polytropic models

The interiors of some celestial bodies are reasonably well described by gases
with polytropic equations of state, such that $P(r)=K\,\rho(r)^{1+1/n}$, where
$P(r)$ is the pressure, $\rho(r)$ is the density, $K$ is a proportionality
constant and $n$ is the polytropic index. This index evolves in mass from
$n\gtrsim 0$, corresponding to rocky planets with Earth-like masses, to $n=3$,
corresponding to massive stars with a radiative core [108]. For intermediate
masses, from Jovian planets up to low-mass brown dwarfs, objects are well
approximated by $n\simeq 1$, and brown dwarfs by $n\simeq 3/2$. This variation
covers a range of about nine orders of magnitude in mass. The main advantage
of this type of models is that pressure only depends on density, so only the
hydrostatic and Poisson equations are needed, with no reference to heat
transfer or thermal balance. Although this might seem an oversimplification,
these models have proven to be remarkably useful in the interpretation of many
features of the structure of celestial bodies and have already been used in
the context of DM capture and evaporation in stars [6, 49, 50, 43, 109, 95],
so we consider them to obtain a generic description. Therefore, we first
describe the distribution of density, pressure and temperature of polytropic
models.

We consider a celestial body with mass $M$ and radius $R$ constituted of a
material with an equation of state of a polytrope of index $n$. The Poisson
and hydrostatic equations (assuming spherical symmetry) can be written as

$\displaystyle\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(r^{2}\,\frac{\partial\Phi}{\partial r}\right)$ $\displaystyle=$
$\displaystyle 4\pi G\,\rho(r)\leavevmode\nobreak\ ,$ (27)
$\displaystyle\frac{\partial P}{\partial r}$ $\displaystyle=$
$\displaystyle-\frac{\partial\Phi}{\partial r}\,\rho(r)\leavevmode\nobreak\ .$
(28)

By substituting the expression of the pressure in terms of the density for a
polytrope, the hydrostatic equation can be integrated. Setting the
gravitational potential $\Phi(R)=0$ at the surface $\rho(R)=0$, one gets

$\rho(r)=\left(\frac{-\Phi(r)}{(n+1)\,K}\right)^{n}\leavevmode\nobreak\ .$
(29)

Now, we introduce the dimensionless quantities $\theta_{n}$ and $\xi$,

$\theta_{n}=\frac{\Phi}{\Phi_{c}}=\left(\frac{\rho}{\rho_{c}}\right)^{1/n}\hskip
28.45274pt;\hskip 28.45274pt\xi=\frac{r}{r_{n}}\hskip 14.22636pt,\hskip
14.22636ptr_{n}^{2}=\frac{(n+1)^{n}\,K^{n}}{4\pi
G\,(-\Phi_{c})^{n-1}}=\frac{(n+1)\,K}{4\pi
G\,\rho_{c}^{1-1/n}}\leavevmode\nobreak\ ,$ (30)

where the subscript c indicates quantities evaluated at the center of the
object. Note that $P(r)=P_{c}\,\theta^{1+n}_{n}(r)$. Substituting Eq. (29)
into Eq. (27), in terms of the variables defined in Eq. (30), the Lane-Emden
equation is obtained,

$\frac{1}{\xi^{2}}\,\frac{d}{d\xi}\left(\xi^{2}\,\frac{d\theta_{n}}{d\xi}\right)+\theta_{n}^{n}=0\leavevmode\nobreak\
.$ (31)

Finite solutions require $d\theta_{n}(0)/d\xi=0$ and in order for $\rho_{c}$
to represent the density at the center, $\theta_{n}(0)=1$. With these boundary
conditions, the equation can be solved, although in most cases only
numerically. These are the so-called E-solutions. Given the polytropic index
$n$, the density profile is fully determined in terms of the total mass and
radius. The mass is given by

$M=4\pi\,\rho_{c}\,R^{3}\,\left.\left(-\frac{1}{\xi}\frac{d\theta_{n}}{d\xi}\right)\right|_{\xi=\xi_{n}}\equiv
4\pi\,\rho_{c}\,R^{3}\,\frac{|\theta_{n}^{\prime}(\xi_{n})|}{\xi_{n}}\leavevmode\nobreak\
,$ (32)

where $\xi_{n}=R/r_{n}$, so $\theta_{n}(\xi_{n})=0$. This results in

$R=\left(\frac{(n+1)\,K}{\left(4\pi\right)^{1/n}G}\,\xi_{n}^{1+1/n}\,|\theta_{n}^{\prime}(\xi_{n})|^{1-1/n}\right)^{-\frac{n}{n-3}}\,M^{\frac{1-n}{3-n}}\leavevmode\nobreak\
,$ (33)

which is the well-known polytropic mass–radius relation. Similarly, the
central density can be written in terms of the total mass,

$\rho_{c}=\left(\frac{(n+1)\,K}{(4\pi)^{1/3}\,G}\,\left(\xi_{n}^{2}\,|\theta_{n}^{\prime}(\xi_{n})|\right)^{2/3}\right)^{\frac{3n}{n-3}}\,M^{\frac{2n}{3-n}}\leavevmode\nobreak\
,$ (34)

and $\rho(r)=\rho_{c}\,\theta_{n}^{n}(r)$. The escape velocity at the core can
be written in terms of the polytropic index and the escape velocity at the
surface,
$v_{e}^{2}(r=0)=\left(1+\left(\xi_{n}\,|\theta_{n}^{\prime}(\xi_{n})|\right)^{-1}\right)\,v_{e}^{2}(r=R)$.

These results correspond to zero temperature (full degeneracy). In the case of
partial degeneracy, temperature is proportional to the degeneracy parameter
and the Fermi energy. For a non-relativistic partially degenerate gas, the
Fermi energy is proportional to $\rho^{2/3}$ and pressure to $\rho^{5/3}$, as
in the case of a polytrope with index $n=3/2$. Therefore, $P\sim\rho\,T$,
which coincides with the equation of state of an ideal monoatomic gas. In
general, for such an equation of state, the temperature profile goes as

$T(r)=T_{c}\,\theta_{n}(r)\sim\rho(r)^{1/n}\hskip
14.22636pt\longrightarrow\hskip 14.22636ptT_{c}\sim
M^{\frac{2}{3-n}}\sim\frac{M}{R}\leavevmode\nobreak\ .$ (35)

Note that this is just the virial theorem for a body of gas in hydrostatic
equilibrium. With these tools at hand, along with empirical mass–radius and
mass–core temperature relations, we can now describe the relevant properties
of different celestial objects. We comment on each type of them below.

###  III.2 Planetary bodies

Planetary bodies are the smallest celestial objects heavy enough so that self-
gravity can force them into spherical shape. Among them, those orbiting a star
and having cleared the neighborhood around their orbit are called planets,
whereas if the last condition is not met they are dwarf planets or satellites.
Their composition is classified in three main types: rock made up of silicates
(compounds of Mg, Si, O) and iron; ice (water and other nolecules as CH4, NH3,
CO, N2 and CO2); and gas (hydrogen and helium).

Both the mechanical and thermal profiles in the interior of planets are fully
determined by the equation of state of the material. For spherical bodies of
homogeneous composition and at zero temperature, in hydrostatic equilibrium,
their radius increases with mass up to a maximum value [110, 111, 112, 113],
which depends on the planet composition. Below this critical mass,
$\sim\textrm{few}\leavevmode\nobreak\ M_{J}\sim 10^{3}\leavevmode\nobreak\
M_{\oplus}\sim\textrm{few}\leavevmode\nobreak\ 10^{-3}\leavevmode\nobreak\
M_{\odot}$, where $M_{J}$ is Jupiter’s mass,999Note that a significant H/He
envelope would result in a larger radius [114, 115]. Another effect that could
enlarge the radius of a planet is irradiation due to the proximity to the
parent star [116, 117], which slows down contraction. These effects could
modify the estimate of the maximum radius at zero temperature. Moreover, for
large masses, the hydrogen in the core can fuse and thermal effects become
important, so the zero-temperature estimate does not apply. Indeed, young, and
still hot, giant planets can be larger than the quoted maximum size.  Coulomb
forces balance gravity and the density varies little and roughly corresponds
to the value at zero pressure, so that $M\sim R^{1/3}$ (or $n\sim
0$).101010Solid and liquid materials are not infinitely incompressible, so
$n>0$ and thus, the exponent in the mass–radius scaling relation is smaller
than 1/3 [118, 119]. Around this critical mass, pressure is large enough to
ionize the material and electron degeneracy pressure starts playing a role.
Thus, for $M>\textrm{few}\leavevmode\nobreak\ M_{J}$, the mass–radius relation
tends to $R\sim M^{-1/3}$ (or $n=3/2$), as will be described below.

The minimum mass required for a planetary body to achieve a nearly spherical
shape in hydrostatic equilibrium is estimated to be $\sim 3\times
10^{-5}\leavevmode\nobreak\ M_{\oplus}\sim 10^{-10}\leavevmode\nobreak\
M_{\odot}$ [74, 75], which roughly corresponds to the smallest round
satellites and dwarf planets in the solar system [120, 121, 122, 123].
Although there exist spherical satellites as massive as $3\times
10^{-4}\leavevmode\nobreak\ M_{\oplus}$ [120] which are not in hydrostatic
equilibrium, in this section we consider the range of masses that represents
planetary bodies to be $(3\times 10^{-5}\leavevmode\nobreak\
M_{\oplus}\simeq)\,10^{-10}\leavevmode\nobreak\ M_{\odot}\lesssim M\lesssim
13\leavevmode\nobreak\ M_{J}\,(\simeq 0.012\leavevmode\nobreak\ M_{\odot})$.

Regarding the structural properties of planets, although the polytropic model
was considered as a first attempt to describe them, that equation of state
does not incorporate the approximate incompressibility of solids and liquids
at low pressures. To account for it, many different equations of state have
been considered [124, 125, 126, 127, 128, 129, 130, 131, 115, 132, 133, 134].
The resulting shape of the mass–radius relation, in the case of a homogeneous
composition, is very similar for different materials, all cases presenting a
maximum radius, with planets with heavier elements resulting in smaller radii
[114, 115]. In fact, this can be used to infer the bulk composition of
planets. But not only does the size of planets depend on mass, but their
composition also does. As a consequence of low-mass cores not being able to
efficiently accrete large amounts of gas during the early formation stages,
low-mass planets are composed mainly of heavy elements and high-mass planets
of gas (see, e.g., Refs. [135, 136, 137]). We account for this in a simplified
way, by using a variable composition along the following empirical piece-wise
mass–radius relation [97] (see also Refs. [138, 119, 139, 140, 141, 142,
143]),

$\left(\frac{R}{R_{\oplus}}\right)=\begin{cases}\,1.0\,\left(\frac{M}{M_{\oplus}}\right)^{0.28}\hskip
32.72049pt;\hskip 17.07164pt10^{-10}\,M_{\odot}<M\leq 2\,M_{\oplus}\hskip
24.18483pt\textrm{(rocky planets and satellites)}\\\\[8.61108pt]
\,0.81\,\left(\frac{M}{M_{\oplus}}\right)^{0.59}\hskip 28.45274pt;\hskip
36.98857pt2\,M_{\oplus}<M\leq 130\,M_{\oplus}\hskip 14.22636pt\textrm{(icy and
small gaseous planets)}\\\\[8.61108pt]
\,20.6\,\left(\frac{M}{M_{\oplus}}\right)^{-0.076}\hskip 19.91692pt;\hskip
22.76219pt130\leavevmode\nobreak\ M_{\oplus}<M\leq 13\,M_{J}\hskip
22.76219pt\textrm{(gas giants)}\end{cases}\leavevmode\nobreak\ ,$ (36)

and we have extrapolated the fit down to $M=10^{-10}\leavevmode\nobreak\
M_{\odot}$.111111This relation approximately applies to Earth-like satellites
and dwarf planets. Note that icy satellites, which are a large fraction of the
satellites in our solar system, are less dense and cooler. In any case, these
bodies are too small to be of much interest in the context of the capture of
DM and the estimate of the DM evaporation mass. We keep them for completeness,
but we only consider rock structures. We only consider the fit in Ref. [97] up
to the icy–gas giants limiting mass, $M=130\leavevmode\nobreak\ M_{\oplus}$.
In the gas giants mass range we ensure continuity with the brown dwarfs regime
(see next subsection) with a power-law relation. For brown dwarfs, we do not
use the results of Ref. [97] as only a few objects are included, so the fit is
likely dominated by giant planets. The parameterization we use results in
slightly smaller ($\lesssim 10\%$) radii, which implies smaller DM evaporation
masses.

These three regimes correspond to rocky planets (Earth-like planets, but also
satellites and dwarf planets), ice and small gas giant planets (Neptune- and
Saturn-like) and gas giant planets (Jupiter-like and super-Jupiters). The
first transition is expected to be caused by the accretion of substantial
volatile gas envelopes [144], but also by the range of possible compositions
of super-Earths, which can have substantial amounts of water and hence, can
have a larger radii. Note that the purely statistical value for this
transition obtained with the fit is slightly lower than the theoretically
expected minimum mass for planets to retain a gaseous atmosphere by
gravitational instability, $\sim 10\leavevmode\nobreak\ M_{\oplus}$ [145,
146]. Additionally, planets more massive than $\sim(30-60)\leavevmode\nobreak\
M_{\oplus}$ are likely to be predominantly composed of gas (hydrogen and
helium) [147], which represents the transition between icy and small gaseous
giant planets, although this is not visible in the mass–radius relation above.
The second transition is likely a consequence of gravitational self-
compression as mass increases, such that the growth of planets with mass
stops. It is also a consequence of a possible bias in observations, as
transiting exoplanets are easier to observe if they are close to their parent
star, where irradiation effects are important and could result in larger
volumes [116, 117].

For rocky planets, the inner temperature is too low to have a significant
impact on their size [118, 148, 114, 115]. Nevertheless, it is high enough to
(totally or partially) melt the interior materials, causing their separation
according to density. This results in differentiated planets, with denser
material lying beneath less dense material, as is the case, with different
degrees, of all planets in our solar system. And the radius of realistic
planets with mixed composition lies between that of homogeneous planets
composed of the denser and less dense materials [114, 115].

Unlike for the Earth, the interior of Earth-like planets is poorly known.
Nevertheless, by solving planetary structures with realistic equations of
state across the mass range $(0.1-30)\leavevmode\nobreak\ M_{\oplus}$ for two-
layer rocky planets, a core radius fraction $\textrm{CRF}=\sqrt{1/3}$
(appropriate for the Earth) is found to represent a reasonable assumption. In
this way, we consider as a default model for rocky planets, the following two-
layer profile [149],

$\rho(r)=\begin{cases}\,\frac{1}{\textrm{CRF}}\,\bar{\rho}\hskip
42.67912pt;\hskip 39.83368pt0\leq r\leq(\textrm{CRF})\,R\\\\[8.61108pt]
\,\frac{2}{3}\frac{R}{r}\,\bar{\rho}\hskip 46.94687pt;\hskip
5.69054pt(\textrm{CRF})\,R<r\leq R\end{cases}\hskip 14.22636pt;\hskip
5.69054pt\textrm{(rocky planets and satellites)}\leavevmode\nobreak\ ,$ (37)

where $\bar{\rho}=3\,M/(4\,\pi\,R^{3})$ is the average density of the planet.
This profile results in a core mass fraction
$\textrm{CMF}=\textrm{CRF}^{2}=1/3$. Regarding the constituents, we consider a
mass composition with a 2:1 ratio of silicates (MgSiO3) to iron [150], with an
iron core and a rocky mantle, which is representative of the Earth. Note,
however, that the composition can be substantially different, including iron
planets and ice satellites. This has little impact on the DM evaporation mass,
though. Similarly to other quantities, the DM evaporation mass depends only
logarithmically (except at mass-matching values for $v_{e}<v_{d}$) on the
atomic number of the target particles. We have checked some extreme
(unrealistic) composition models and the DM evaporation mass gets modified by
$\lesssim 10\%$. Being beyond the scope of this work, we do not pursue a more
detailed account of variations on planetary composition, but realistic
uncertainties on the DM evaporation mass caused by this are expected to be
small.

As the energy transport in the interior of planets is likely dominated by
convection within layers [138], the temperature profile is quasi-adiabatic
within layers. Similar to the simplified density profile, we consider a two-
layer constant-temperature profile, with the temperature scaling as
$T(r)\propto\rho(r)$ (i.e., assuming a Grüneisen parameter $\gamma=1$).
Estimating the interior temperatures is a rather non-trivial issue, even for
the Earth, and depends on the poorly known form of the equation of state at
high pressures and its temperature dependence. In this work, we consider a
scaling of the core temperature with the mass of Earth-like and dwarf planets
and satellites as

$T_{c}=4500\leavevmode\nobreak\
\textrm{K}\,\left(\frac{M}{M_{\oplus}}\right)^{0.3}\hskip 28.45274pt;\hskip
14.22636pt10^{-10}\leavevmode\nobreak\ M_{\odot}\leq M\leq
2\leavevmode\nobreak\ M_{\oplus}\leavevmode\nobreak\ ,$ (38)

which slightly underestimates the expectations for planets and satellites for
$M\gtrsim 10^{-8}\leavevmode\nobreak\ M_{\odot}$ and slightly overestimates
them for smaller masses [151, 152, 153, 154, 73, 155]. For instance, it
underestimates Earth’s expectations by $\sim 30\%$ [155].

For icy and gaseous giant planets, we consider the density and temperature
profiles corresponding to a polytrope of index $n=1$, although density
profiles get steeper with mass as a consequence of the increasing electron
degeneracy. This approximation is in reasonable agreement with expectations
for continuous radial density profiles of planets like Neptune or Uranus [156]
and Saturn or Jupiter [157, 158]. Discontinuous density profiles have also
been derived, with a core extending up to 70% of the planet radius [159, 160,
161, 162]; in this case with an equivalent polytropic index $n\sim 0.6-1$
[163]. These models cannot be distinguished, however, as they both fit well
the measured gravitational potential.

As for the mass scaling of the core temperature, we use estimates for gas
giant planets in the solar system and from Ref. [164] (for solar metallicity
at 5 Gyr) to perform a fit. This facilitates continuity at
$M=13\leavevmode\nobreak\ M_{J}$ with the relation for brown dwarfs (see
below). Note, however, that a significant range of temperatures is possible,
depending on a number of factors (atmospheric conditions, composition, age,
proximity to the parent star, etc.). The following approximation is intended
to describe relatively cool bodies in this mass range,

$\log\left(\frac{T_{c}}{10^{4}\leavevmode\nobreak\
\textrm{K}}\right)=0.28+0.56\,\log{\left(\frac{M}{M_{J}}\right)}+0.22\,\log^{2}{\left(\frac{M}{M_{J}}\right)}+0.035\,\log^{3}{\left(\frac{M}{M_{J}}\right)}\hskip
14.22636pt;\hskip 14.22636pt2\leavevmode\nobreak\ M_{\oplus}<M\leq
13\leavevmode\nobreak\ M_{J}\leavevmode\nobreak\ ,$ (39)

and reproduces reasonable well the results in Refs. [138, 136, 165, 166, 167,
168, 169] for ice and gas giant planets, in particular, for the ones in our
solar system. The increase of the slope with mass is mainly due to the phase
transition from molecular to metallic (liquid) hydrogen at around $10^{11}$ Pa
and $T\sim(10^{3}-10^{4})$ K, when the average energy of electrons becomes
higher than the ionization potential of hydrogen (see, e.g., Ref. [163] and
references therein). This results in a faster increase of temperature with
pressure, and thus, with mass. Note also that super-Earths are usually warmer
(and denser), so Eq. (39) also slightly underestimates the core temperature of
the most massive rocky planets [118, 138]. This results in a slight
underestimation of the DM evaporation mass, as we discuss below.

Concerning the element composition, different materials would likely result in
different radius for a given mass. In general, modeling the overall structure
and composition is fraught with degeneracies, which we cannot carefully
account for in our simplified treatment. Here, we consider a continuous
transition from rocky/icy planets to gaseous giant planets based on pebble
accretion [170, 171, 172, 173, 174], which we parameterize by performing a fit
to the results in Ref. [175] corresponding to planet formation at 10 au (about
the average distance of Saturn to the Sun). In this way, the mass fraction in
elements heavier than hydrogen or helium is given by

$\log
Z=\begin{cases}\,-0.015-0.042\,\log\left(\frac{M}{M_{\oplus}}\right)+0.027\,\log^{2}\left(\frac{M}{M_{\oplus}}\right)-0.020\,\log^{3}\left(\frac{M}{M_{\oplus}}\right)\hskip
14.22636pt;\hskip 25.60747pt2\leavevmode\nobreak\ M_{\oplus}\leq M\leq
33.7\leavevmode\nobreak\ M_{\oplus}\\\\[8.61108pt]
\,1.44-\log\left(\frac{M}{M_{\oplus}}\right)\hskip 223.3543pt;\hskip
14.22636pt33.7\leavevmode\nobreak\ M_{\oplus}<M\leq 13\leavevmode\nobreak\
M_{J}\end{cases}\leavevmode\nobreak\ .$ (40)

The hydrogen plus helium mass fraction is given by $X+Y=1-Z$, and we assume
solar relative composition (i.e., X/Y = 3). For $M>2\leavevmode\nobreak\
M_{\oplus}$, as a benchmark for heavy elements, we consider water. Note that
the so-called pebble isolation mass ($M_{\rm PIM}\simeq 34\leavevmode\nobreak\
M_{\oplus}$ in this case) is smaller the closer the orbit is, which results in
a larger fraction of heavy elements for $M\lesssim M_{\rm PIM}$ and in a
smaller fraction for $M\gtrsim M_{\rm PIM}$ [175].

###  III.3 Brown dwarfs

Now we consider more massive celestial bodies up to masses for which hydrogen
fusion takes place efficiently, which defines stars. As mentioned above, at
zero temperature there is a critical mass that sets the maximum size of
planets and the transition to brown dwarfs, which is caused by a critical
pressure, $\sim(10-10^{3})$ GPa (smaller values for lighter elements), above
which chemical bonds of the material get broken. Thus, after this pressure is
reached, the core of a planet shrinks and its density grows. At this point,
ionized electrons start becoming partially degenerate and their pressure
increasingly significant. As mass grows, the description of planetary bodies
in terms of polytropes experiences a transition from an index $n\gtrsim 0$
(constant density) to $n=3/2$ (non-relativistic degenerate gas). In the limit
of a non-relativistic fully degenerate gas, from Eq. (33) we would expect
$R\sim M^{-1/3}$. However, Coulomb pressure by ions (which tends to set a
constant density with a fixed inter-particle distance scale, $R\sim M^{1/3}$)
almost cancels the electron degeneracy pressure, which results in $R\sim
M^{-1/8}$ (or $n=11/9$ for a polytrope) [164]. The accurate description of
these effects is governed by the equation of state around these pressures,
which requires accounting for strongly correlated, polarizable, partially
degenerate quantum and classical plasmas, in a medium where partial ionization
by pressure becomes important. Moreover, the equation of state is not the only
factor affecting this relation, which also changes with time as the object
cools down, resulting in a smaller radius for older objects and with a
mass–radius relation which is even flatter for younger bodies [176]. Thus, in
realistic models, the transition between planets and brown dwarfs takes place
at higher masses than the prediction at zero temperature due to different
effects.

The transition between giant planets and brown dwarfs is a matter of a long-
standing debate [177, 178, 164, 179, 180, 181]. In general, it is believed
that planets are richer in heavy elements than their parent stars, whereas
brown dwarfs share a similar composition, which would hint at a distinction
based on the formation mechanism [177]. Nevertheless, the current distinction
is set in terms of mass limits based on nuclear fusion processes. For
$M\gtrsim 13\leavevmode\nobreak\ M_{J}\simeq 0.0124\,M_{\odot}$ (for solar
metallicity) [182, 183, 184], deuterium can start to fuse ($T_{c}\gtrsim
5\times 10^{5}$ K), whereas effective hydrogen burning does not take place for
masses below $\sim 0.07\,M_{\odot}\simeq 73\leavevmode\nobreak\ M_{J}$ (which
defines stars, $T_{c}\gtrsim 3\times 10^{6}$ K) [185, 186]. Within this range
of masses, celestial bodies are usually defined as brown dwarfs.

The empirical mass–radius relation for brown dwarfs can be approximately given
by a fit to the results of Ref. [164] (for solar metallicity and an age of 5
Gyr),121212Gas giants continue to cool down and contract after $\sim$1 Myr,
whereas brown dwarfs heat up before $\sim$1 Gyr, but when electron degeneracy
dominates, they cool down and contract [187, 164, 179]. In any case, we
consider relations relevant for most part of the life of celestial objects,
corresponding to late evolutionary stages.

$\log\left(\frac{R}{R_{J}}\right)=0.063+0.0036\,\log\left(\frac{M}{M_{J}}\right)-0.055\,\log\left(\frac{M}{M_{J}}\right)^{2}\hskip
28.45274pt;\hskip 14.22636pt13\leavevmode\nobreak\ M_{J}<M\leq
0.07\,M_{\odot}\leavevmode\nobreak\ .$ (41)

As already mentioned, brown dwarfs are well described by non-relativistic
partially degenerate polytropes ($n\simeq 3/2$), which are fully determined in
terms of their mass and radius. Thus, the mass–radius relation reduces the
number of free parameters needed to determine the shape of the density and
temperature profiles to one, say the mass. In the case of the density, the
constraint on the volume integral (determined by the mass), completely fixes
its central value. Thus, our default model follows a polytropic density
profile with the temperature profile given by $T(r)\sim\rho(r)^{2/3}$, with
the mass–core temperature relation obtained from a fit to the different brown
dwarfs in Ref. [164] (for solar metallicity and an age of 5 Gyr),

$\log\left(\frac{T_{c}}{{\rm
K}}\right)=0.34+9.36\,\log\left(\frac{M}{M_{J}}\right)-6.20\,\log^{2}\left(\frac{M}{M_{J}}\right)+1.56\,\log^{3}\left(\frac{M}{M_{J}}\right)\hskip
28.45274pt;\hskip 14.22636pt13\leavevmode\nobreak\ M_{J}<M\leq
0.07\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ .$ (42)

These results are in agreement with others in the literature (see, e.g., Refs.
[187, 179]) and were obtained with the equation of state for hydrogen and
helium from Ref. [188], which covers the relevant range of densities and
temperatures (see also Refs. [165, 166, 189, 167, 168, 169] for recent
calculations with ab initio equations of state and Ref. [190] for a comparison
of them). Note that updated equations of state predict not only higher central
temperatures, but also larger radii, and different atmospheric boundary
conditions could also result in higher central temperatures [189]. On another
hand, metallicities smaller than solar values result in slightly smaller radii
and higher core temperatures, although these variations only amount to a few
percent (see, e.g., Ref. [164]) and have a small impact on the calculation of
the DM evaporation mass. At late evolutionary stages, brown dwarfs cool down.
We consider $\sim 5$ Gyr as a representative value, but the core temperature
of the oldest ($\gtrsim 10$ Gyr) and most massive ($M\gtrsim
40\leavevmode\nobreak\ M_{J}$) brown dwarfs is expected to be smaller by
$\lesssim 20\%$ (see, e.g., Refs. [164, 179]). Note, however, that accounting
for the formation of dust in the equation of state results in slightly slowing
down the cooling process for the most massive brown dwarfs [191], so the
actual difference in core temperature between 5 Gyr and 10 Gyr is expected to
be even smaller.

Finally, for the composition of brown dwarfs, we take 98% hydrogen and helium
(with solar 3:1 proportion) and 2% heavy elements (water), $X+Y=0.98$ and
$Z=0.02$ [164, 179].

###  III.4 Main-sequence stars

For masses around $\sim 0.07\leavevmode\nobreak\ M_{\odot}$, celestial bodies
contract as they radiate from the surface and the core heats up. This
contraction eventually stops either by the appearance of electron degeneracy
or by thermal pressure from hydrogen burning. The most dominant process is
mainly determined by the attained temperatures. For low enough temperatures,
electrons become partially degenerate and hydrogen cannot fuse. Given the weak
dependence of pressure on temperature, further radiation from the surface does
not result in contraction, but in the cooling of the core. These objects are
called brown dwarfs and have been discussed above. On the other hand, if
central temperatures are high enough for hydrogen to ignite before electrons
become degenerate, energy losses from the surface are balanced by the thermal
energy from hydrogen burning and the object inflates. These celestial bodies
are called stars. The transition mass between these two types of objects is
$\sim 0.07\leavevmode\nobreak\ M_{\odot}$ [185, 186, 178]. Therefore,
thermonuclear reactions involving hydrogen create a thermal pressure that can
sustain the gravitational pull of the gas and make stars bigger. This results
in the non-degeneracy of the gas, which approximately behaves as a classical
ideal gas, $P\propto\rho\,T$.

In this section we only consider stars during main-sequence. This is the
evolutionary phase since ignition after gravitational collapse until the
hydrogen core is consumed and converted into helium. This stage is the longest
one and represents about 90% of the stars life, when they can be approximately
considered to be in hydrostatic equilibrium, with properties varying
relatively little.

For masses above the hydrogen-burning limit, $M\gtrsim
0.07\leavevmode\nobreak\ M_{\odot}$, electron degeneracy decreases quickly
with mass, but up to $M\sim 0.4\leavevmode\nobreak\ M_{\odot}$, it still
impacts star evolution and the inner structure can be described in terms of a
polytrope of index $n=3/2$, as it happens for brown dwarfs. For $M\lesssim
0.2\leavevmode\nobreak\ M_{\odot}$, stars are fully convective and their
internal structure strongly depends on the boundary conditions. As in the case
of brown dwarfs, convection deep into the optically thin layers is favored by
the formation of molecular hydrogen in the outer parts. The presence of
molecules enhances collision-induced absorption, so radiative opacity
increases, which in turn, reduces the radiative transport efficiency.
Additionally, the adiabatic gradient is reduced, which favors convection
[185]. For slightly heavier stars, convection penetrates less efficiently
their inner parts and a radiative core forms. The transition from a convective
star to the formation of a significant radiative core occurs at $M\sim
0.4\leavevmode\nobreak\ M_{\odot}$. However, hydrogen burning is not very
efficient up to $M\sim 0.7\leavevmode\nobreak\ M_{\odot}$ and stars can be
approximately described by a polytrope of index $n=2$ in that mass range
(i.e., $0.4\leavevmode\nobreak\ M_{\odot}\lesssim M\lesssim
0.7\leavevmode\nobreak\ M_{\odot}$). For $0.7\leavevmode\nobreak\
M_{\odot}\lesssim M\lesssim 1.5\leavevmode\nobreak\ M_{\odot}$, the decrease
in the central abundance of hydrogen results in an increase of the molecular
weight, which heats up further and inflates the star more efficiently. For
$M\gtrsim 1.5\leavevmode\nobreak\ M_{\odot}$, stars reach high enough
temperatures ($T_{c}\gtrsim 2\times 10^{7}$ K) so that the CNO burning cycle
ignites in most of their core. The extreme temperature dependence of these
reactions results in convective instability and the formation of a convective
core, which is larger the heavier the star is. This is also favored by the
increasing importance of radiation pressure with mass, which reduces the
adiabatic gradient.131313Note that the condition for convection instability is
that the temperature gradient is steeper than the adiabatic one. This can be
achieved either by reducing the adiabatic gradient (as is the case in the
outer layers in low-mass stars) or by increasing the temperature gradient (as
happens when the CNO cycle is efficient in high-mass stars). Thus, for masses
$M\gtrsim 0.7\leavevmode\nobreak\ M_{\odot}$, stars can be approximately
described by a polytrope of index $n=3$.

As a consequence of all the above features, from the hydrogen-burning limiting
mass to very massive stars, the radius and core temperature continuously
increase with mass. At the very low-mass end, $0.07\leavevmode\nobreak\
M_{\odot}\lesssim M\lesssim 0.4\leavevmode\nobreak\ M_{\odot}$, the density of
stars decreases with mass. Nevertheless, within the mass range
$0.4\leavevmode\nobreak\ M_{\odot}\lesssim M\lesssim 1.5\leavevmode\nobreak\
M_{\odot}$, the gas is non-degenerate and the nuclear rate of energy
production via hydrogen fusion grows in efficiency, resulting in a warmer
($T_{c}$ with a stronger dependence on the mass than in less massive stars)
and a denser core, so the luminosity grows steeply with mass. For $M\gtrsim
1.5\leavevmode\nobreak\ M_{\odot}$, radiation pressure reduces the mass
dependence of the luminosity, and the very steep temperature dependence of the
CNO rate results in hotter and less dense stars the more massive they are.

Although a proper treatment of the evolution equations should include the full
system of hydrostatic equilibrium, mass conservation, energy transfer and
energy conservation equations (along with the chemical equation that governs
the changes in composition), the approximation of the density and temperature
profiles in terms of polytropes works reasonable well and it is good enough
within our simplistic modeling. In particular, as mentioned above, we consider
$n=3/2$ for $0.07\leavevmode\nobreak\ M_{\odot}\leq M\leq
0.4\leavevmode\nobreak\ M_{\odot}$, $n=2$ for $0.4\leavevmode\nobreak\
M_{\odot}<M\leq 0.7\leavevmode\nobreak\ M_{\odot}$ and $n=3$ for
$M>0.7\leavevmode\nobreak\ M_{\odot}$ [192] to obtain the density and
temperature profiles.

As for the other less massive celestial bodies discussed in previous
subsections, empirical mass–radius relations for stars have also been proposed
for a long time (see, e.g., Refs. [193, 194, 195, 107] for some recent
relations). Here, we consider the following empirical mass–radius relation,
inspired by and in agreement with the data and models of Refs. [164, 107, 106]

$\log\left(\frac{R}{R_{\odot}}\right)=\begin{cases}\,1.17\,\log\left(\frac{M}{M_{\odot}}\right)+0.67\,\log^{2}\left(\frac{M}{M_{\odot}}\right)+0.43\,\log^{3}\left(\frac{M}{M_{\odot}}\right)\hskip
14.22636pt;\hskip 14.22636pt0.07\leavevmode\nobreak\ M_{\odot}<M\leq
1.5\leavevmode\nobreak\ M_{\odot}\\\\[8.61108pt]
\,0.13+0.56\,\log\left(\frac{M}{M_{\odot}}\right)\hskip 145.10922pt;\hskip
17.07164pt1.5\leavevmode\nobreak\ M_{\odot}<M\leq 100\leavevmode\nobreak\
M_{\odot}\end{cases}\leavevmode\nobreak\ .$ (43)

For low-mass stars, we perform a fit using the data compiled in Ref. [106],
which is in good agreement with Refs. [164, 107] for $M\leq
1.5\leavevmode\nobreak\ M_{\odot}$. The (log) linear parameterization in the
high-mass range is a simplification of the piece-wise fit in Ref. [107].

As for any polytropic model, this relation reduces the number of parameters
required to describe the internal stellar structure to one and, as done
throughout this paper, we use the mass of the star. Note that using the
mass–radius relation above, the core density, $\rho_{c}\sim M/R^{3}$, always
decreases with the stellar mass. Nevertheless, we have mentioned above that in
the mass interval $0.4\leavevmode\nobreak\ M_{\odot}\lesssim M\lesssim
1.5\leavevmode\nobreak\ M_{\odot}$, detailed stellar models predict the core
density to slightly grow with mass. We cannot correct for this trend using the
polytrope assumption, but given that it is known to provide an approximately
correct overall description, we use it in what follows. This approximation has
a negligible impact on the calculation of the DM evaporation mass.

The situation is different for the core temperature, which is not fixed by the
mass–radius relation. We obtain a monotonously increasing mass–core
temperature relation using the results corresponding to the models for low-
mass stars from Ref. [164] and for high-mass stars from Ref. [196] (at zero-
age main sequence),

$\log\left(\frac{T_{c}}{10^{4}\leavevmode\nobreak\
\textrm{K}}\right)=\begin{cases}3.17+0.82\,\log\left(\frac{M}{M_{\odot}}\right)-0.25\,\log^{2}\left(\frac{M}{M_{\odot}}\right)-2.15\,\log^{3}\left(\frac{M}{M_{\odot}}\right)-1.60\,\log^{4}\left(\frac{M}{M_{\odot}}\right)\\\
\hskip 298.75394pt;\hskip 14.22636pt0.07\leavevmode\nobreak\ M_{\odot}<M\leq
1.5\leavevmode\nobreak\ M_{\odot}\\\
3.24+0.30\,\log\left(\frac{M}{M_{\odot}}\right)-0.054\,\log^{2}\left(\frac{M}{M_{\odot}}\right)\hskip
14.22636pt;\hskip 14.22636pt1.5\leavevmode\nobreak\ M_{\odot}<M\leq
100\leavevmode\nobreak\ M_{\odot}\end{cases}\leavevmode\nobreak\ .$ (44)

We have also confirmed that these are reasonable parameterizations by
comparing them with the evolutionary tracks obtained with the MIST code [197,
198], that in turn uses the stellar evolution code MESA [199, 200, 201].

Note that for $M\sim(0.2-0.4)\leavevmode\nobreak\ M_{\odot}$, the core
temperature is roughly given by $T_{c}\simeq 1.6\,\mu_{\rm mol}\,G\,M/(3\,R)$,
where $\mu_{\rm mol}$ is the mean molecular weight of the material (in units
of the atomic mass unit), which is $\mu_{\rm mol}\simeq 0.6$ for solar
abundances. This is the expected result for a polytrope of index $n=3/2$, for
a monoatomic and fully ionized ideal gas, Eq. (35). This is a consequence of
degeneracy becoming less important with mass, while the star being largely
convective. For $0.4\leavevmode\nobreak\ M_{\odot}\lesssim M\lesssim
1.5\leavevmode\nobreak\ M_{\odot}$, had we used the mass–radius relation
above, assuming the equation of state of an ideal gas would not result in
$T_{c}$ increasing with mass, though. For higher masses, $M\gtrsim
1.5\leavevmode\nobreak\ M_{\odot}$, radiation pressure grows in importance,
which results in a dependence of the core temperature on the mass less steep
than at lower masses, as mentioned above.

Furthermore, unlike what happens for less massive celestial objects (planetary
bodies and brown dwarfs), stars remain in main sequence for a period of time
which is shorter the heavier the star is. It is, thus, relevant to discuss how
long this stage lasts, such that their overall properties remain approximately
constant, without much change. For the DM capture process to be efficient and
the discussion about the DM evaporation mass to be of most interest, the
equilibration time $\tau_{\rm eq}$ (actually, $\tau_{\rm eq}/\kappa$) must be
shorter than the stars age. Otherwise, the DM evaporation mass is time
dependent.

As the total energy release in the $pp$ chain and the CNO cycle is the same, a
simple estimate for the lifetime in main sequence is given by

$t_{\rm life,MS}=t_{\rm
life,\odot}\,\left(\frac{M/M_{\odot}}{L/L_{\odot}}\right)\leavevmode\nobreak\
,$ (45)

where $t_{\rm life,\odot}=10$ Gyr and $L_{\odot}=3.828\times 10^{26}$ W are
the lifetime and current luminosity of the Sun. To obtain the lifetime of
other stars we consider the empirical mass–luminosity relation from Ref.
[107], which we extrapolate, in the low- and high-mass extremes, to cover the
entire stellar mass range discussed in this work. Moreover, we simplify the
multiple piece-wise form by performing a smooth fit in the mass interval
$0.1\leavevmode\nobreak\ M_{\odot}\lesssim M\leq 100\leavevmode\nobreak\
M_{\odot}$,

$\log\left(\frac{L}{L_{\odot}}\right)=4.24\,\log\left(\frac{M}{M_{\odot}}\right)+0.239\,\log^{2}\left(\frac{M}{M_{\odot}}\right)-0.759\,\log^{3}\left(\frac{M}{M_{\odot}}\right)+0.216\,\log^{4}\left(\frac{M}{M_{\odot}}\right)\leavevmode\nobreak\
.$ (46)

Notice that only stars with masses $M\lesssim 2\leavevmode\nobreak\ M_{\odot}$
have lifetimes in the main sequence longer than $\sim 1$ Gyr. Very massive
stars, $M\gtrsim 10\leavevmode\nobreak\ M_{\odot}$ leave the main sequence
after less than $\sim 10$ Myr.

###  III.5 Post-main-sequence stars

The internal properties of stars during post-main-sequence evolution undergo
significant changes in relatively short periods of time, which crucially
depend on the stellar mass they had in main sequence and on environmental
conditions [196]. Given that evaluating the implications of non-static
configurations on the DM evaporation mass is beyond the scope of this paper,
we briefly and qualitatively describe some of the stellar properties after
stars leave the main sequence in this subsection, and leave the description of
compact stellar remnants to the next subsections.

The main-sequence phase ends when hydrogen is exhausted in the stellar core.
At this point, an approximately isothermal helium core develops and hydrogen
burns in a shell surrounding it [202]. This results in the production of more
helium, which increases the mass of the core, bounded by the Schönberg-
Chandrasekhar limit [203], and reduces its size. While the core temperature
slightly increases, the envelope of the star expands and cools, becoming
mostly convective, to favor the transport of additional energy. This is the
so-called subgiant branch for low- and intermediate-mass stars and lasts for
$\sim\textrm{few Gyr}$ for $M\lesssim 1\leavevmode\nobreak\ M_{\odot}$ and
$\sim 10^{7}$ yr for $M\gtrsim 2\leavevmode\nobreak\ M_{\odot}$. During this
phase, although the core temperature does not change much (it is just slightly
higher than the terminal age main-sequence value), the central density can
increase by more than an order of magnitude with respect to the main-sequence
value [204].

Eventually, stars enter the red giant branch, delimited by their Hayashi line,
and their luminosity increases. Very low-mass stars ($M\lesssim
0.3\leavevmode\nobreak\ M_{\odot}$), however, are fully convective and do not
reach temperatures high enough to burn helium, so they do not become red
giants and end up as helium white dwarfs [205], when all hydrogen has been
fused to helium. For low-mass stars, $0.3\leavevmode\nobreak\
M_{\odot}\lesssim M\lesssim 2.3\leavevmode\nobreak\ M_{\odot}$, the density of
the core can get so high that electrons become degenerate and for $M\gtrsim
0.8\leavevmode\nobreak\ M_{\odot}$, the temperature increases up to a point
when helium starts to burn ($T_{c}\sim 10^{8}$ K) [202] and the core mass is
$\sim 0.5\leavevmode\nobreak\ M_{\odot}$ [206]. At this point, as temperature
rises (in a run-away process, as pressure does not depend on temperature), the
core turns non-degenerate again, which finally results in a helium
flash.141414Helium burning does not start at the center of the core. Large
amounts of energy lost in the form of neutrinos during the final stages in the
red giant branch produce a temperature inversion, so that helium burning
starts in a shell around the helium-rich core [207, 208, 209]. This phase
lasts for about $10^{9}$ yr for $M\sim 1\leavevmode\nobreak\ M_{\odot}$ and
$10^{7}$ yr for $M\sim 2\leavevmode\nobreak\ M_{\odot}$. During this phase,
the core temperature increases by a factor of a few, and the central density
is several orders of magnitude larger than during main sequence, although it
decreases much faster towards the edge of the core. On the other hand, the
stellar radius can increase by two orders of magnitude for the lightest stars
and by a factor of a few for the most massive ones in this interval, and a
significant mass loss is experienced.

The end of this phase for low-mass stars, which results in stellar contraction
and core cooling, marks the entry into the horizontal branch stage. This phase
lasts for up to $\sim 10^{8}$ years [206], while helium burns and a carbon-
oxygen core develops, and it is similar, but shorter, than the hydrogen-
burning phase that defines main sequence. The position on the horizontal
branch depends on the mass lost during the red giant phase. The larger the
mass loss, the bluer and slightly smaller the star becomes. In particular, a
$\sim 1\leavevmode\nobreak\ M_{\odot}$ star on this branch is an order of
magnitude larger than during main sequence.

After the helium core is depleted, the immediate stellar evolution strongly
depends on the mass [210, 211, 212]. For low masses ($0.3\leavevmode\nobreak\
M_{\odot}\lesssim M\lesssim 2.3\leavevmode\nobreak\ M_{\odot}$), stars enter
the asymptotic giant branch, which is similar to the red giant branch, but
with slightly higher temperatures. The non-degeneracy of the core, however,
implies that the increase of pressure leads to expansion and cooling of the
region where helium burns.

In heavier stars, $2.3\leavevmode\nobreak\ M_{\odot}\lesssim
M\lesssim(8-10)\leavevmode\nobreak\ M_{\odot}$ [210, 211, 212, 213], helium
ignites before the core becomes degenerate, so there is no helium flash. Stars
can move away from the red giant branch getting bluer and then becoming redder
again (blue loops), and a degenerate carbon-oxygen core develops. The most
massive stars in this mass range might not even reach the red giant branch
before turning into red supergiants (if they retain their envelope) or getting
bluer and becoming Wolf-Rayet stars (if they do not retain their envelope).
This phase lasts for a short period of time, the shorter the more massive
stars are. All these stars, $M\lesssim(8-10)\leavevmode\nobreak\ M_{\odot}$,
end up their lives as white dwarfs. For more massive stars,
$(8-10)\leavevmode\nobreak\ M_{\odot}\lesssim M\lesssim 25\leavevmode\nobreak\
M_{\odot}$, the helium-burning phase is also the red supergiant stage. These
stars, though, proceed through a series of nuclear burning phases and end up
as core-collapse supernovae, leaving neutron stars as remnants. The remnants
of even more massive stars are black holes.

###  III.6 White dwarfs

Once stars with initial masses $M\lesssim(8-10)\leavevmode\nobreak\ M_{\odot}$
reach the tip of the asymptotic red giant branch, they attain their maximum
size and continue losing mass at a fast pace, ejecting whole shells of
material (creating planetary nebulae). At this point, stars get hotter and
become blue supergiants in a very short period of time, $\sim 10^{4}$ yr,
until only the hot core remains. This results in a narrow final range of
masses. However, temperatures are not high enough for nuclear fusion to
proceed and, lacking any other energy sources, stars contract. This increase
in density finally results in electron degeneracy throughout the star, which
creates a new source of pressure that halts gravitational contraction. From
then on, they become white dwarfs (see, e.g., Refs. [214, 215, 216, 196, 217])
and this is the final evolutionary stage for stars with masses
$M\lesssim(8-10)\leavevmode\nobreak\ M_{\odot}$ [218], which represent the
fate of about 97% of all stars. Due to electron degeneracy, energy losses from
the surface cannot be balanced any more by reducing the stellar size, so white
dwarfs keep on cooling for the rest of their lives, unless they accrete matter
from a nearby star.

Analogously to the description of polytropes, the mechanical properties of
white dwarfs can be decoupled from the thermal ones. Nevertheless, the
equation of state is not that of a polytrope. In this case, pressure is
generated by a gas of fully degenerate electrons, but with a varying degree of
relativistic Fermi momentum. The equation that describes the internal
structure of white dwarfs, which is similar to the polytrope model equation,
is the Chandrasekhar equation [219],

$\frac{\text{d}^{2}\varphi}{\text{d}\zeta^{2}}+\frac{2}{\zeta}\,\frac{\text{d}\varphi}{\text{d}\zeta}+\left(\varphi^{2}-\frac{1}{z_{c}^{2}}\right)^{3/2}=0~\leavevmode\nobreak\
,$ (47)

with boundary conditions $\varphi(0)=1$ and $\varphi^{\prime}(0)=0$, and

$\varphi\equiv\frac{z}{z_{c}}\hskip 28.45274pt;\hskip
28.45274pt\zeta\equiv\left(\mu_{e}\,m_{u}\,m_{e}\,\sqrt{\frac{4\,G}{3\,\pi}}\,z_{c}\right)\,r\leavevmode\nobreak\
.$ (48)

Here $z^{2}=\left(p_{\rm F}/m_{e}\right)^{2}+1$, with $p_{\rm F}$ the Fermi
momentum, and $z_{c}=z(r=0)$, with $z\in[1,\infty)$. Note that this equation
reduces to the Lane-Emden equation for polytropes in the limits of
$z\to\infty$ ($n=3$) and $z\to 1$ ($n=3/2$).

The density profile, radius and mass of white dwarfs are given by

$\displaystyle\rho(r)$ $\displaystyle=$
$\displaystyle\left(\frac{\mu_{e}\,m_{u}\,m_{e}^{3}}{3\,\pi^{2}}\right)\left(z_{c}^{2}\,\varphi^{2}-1\right)^{3/2}\simeq
1.95\times 10^{6}\leavevmode\nobreak\
\textrm{g/cm}^{3}\,\left(z_{c}^{2}\,\varphi^{2}-1\right)^{3/2}\leavevmode\nobreak\
,$ (49) $\displaystyle R$ $\displaystyle=$
$\displaystyle\sqrt{\frac{3\,\pi}{4\,G}}\,\frac{\zeta_{1}}{\left(\mu_{e}\,m_{u}\,m_{e}\,z_{c}\right)}\simeq
0.61\,\frac{\zeta_{1}}{z_{c}}\,R_{\oplus}\leavevmode\nobreak\ ,$ (50)
$\displaystyle M$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3\,\pi}}{2}\,\frac{1}{(\mu_{e}\,m_{u})^{2}\,G^{3/2}}\,\left|\zeta_{1}^{2}\,\varphi^{\prime}(\zeta_{1})\right|\simeq
0.72\,\left|\zeta_{1}^{2}\,\varphi^{\prime}(\zeta_{1})\right|\,M_{\odot}\leavevmode\nobreak\
,$ (51)

where $\mu_{e}\simeq 2/(1+X)$ is the mean molecular weight per free electron
and $\zeta_{1}\equiv\zeta(z=1)$. The Chandrasekhar solution does not imply a
lower bound on the radius ($R\to 0$ for $z_{c}\to\infty$), which decreases as
the central density increases, keeping the total mass at a fixed value,
$M_{\rm Ch}\simeq(2/\mu_{e})^{2}\,1.456\leavevmode\nobreak\ M_{\odot}$; this
is the so-called Chandrasekhar mass. For $z_{c}\to 1$, $\zeta_{1}\to\infty$
and $\left|\zeta_{1}^{2}\,\varphi^{\prime}(\zeta_{1})\right|\to 0$ and thus,
$R\to\infty$ and $M~\to 0$. These limiting conditions are modified once
several corrections (general relativity, equation of state, Coulomb
interactions) are incorporated. The resulting maximum values, which slightly
depend on composition, are $M_{\rm max}\sim 1.3\leavevmode\nobreak\ M_{\odot}$
and $R(M_{\rm max})\sim 0.02\leavevmode\nobreak\ M_{\odot}$.

White dwarfs have been observed with masses
$M\sim(0.2-1.3)\leavevmode\nobreak\ M_{\odot}$, with a distribution peaked at
$M\simeq 0.6\leavevmode\nobreak\ M_{\odot}$ [220, 221]. Here, we consider this
range of masses and solve the Chandrasekhar equation to obtain the mass–radius
relation and the density profile. The former is shown in the left panel of
Fig. 5, along with a selection of observed cool white dwarfs with $T_{\rm
eff}<10^{4}$ K [221], which illustrates that Chandrasekhar solution represents
a reasonable approximation, and in general, a conservative one for the
calculation of the DM evaporation mass.

Figure 5: Mass–radius relation for white dwarfs (left panel) and neutron stars
(right panel). We also show a selection of observed cool white dwarfs with
$T_{\rm eff}<10^{4}$ K and $\textrm{S/N}\geq 10$ [221]. For neutron stars, we
use the equation of state from Ref. [222], based on the nuclear energy-density
functional BSk20 [223].

Most white dwarfs have a carbon-oxygen core surrounded by a thin helium
envelope ($Y\lesssim 0.01$), which is in turn surrounded by an even thinner
hydrogen envelope ($X\lesssim 10^{-4}$), although low-mass white dwarfs can
have a helium core (those stars not massive enough, $M\lesssim
0.3\leavevmode\nobreak\ M_{\odot}$, to burn helium) and the most massive white
dwarfs ($M\gtrsim 1\leavevmode\nobreak\ M_{\odot}$) can develope an oxygen-
neon core. As representative of the core composition, we consider white dwarfs
with $Z_{\rm C}=0.4$ and $Z_{\rm O}=0.6$ (see, e.g., Refs. [224, 225]),
although differences on the composition do not affect our results of the DM
evaporation mass.

In order to approximately describe the thermal properties of white dwarfs (at
least of relatively cool and evolved ones), we point out that the core
constitutes more than 99% of their mass and that the main contribution to the
heat capacity comes from the non-degenerate gas of ions. Degenerate electrons
are very efficient in transporting energy outwards, so the core can be
approximately described to be isothermal. Finally, energy is radiated away
through the non-degenerate envelope, which cools down the white dwarf. A
simple description in terms of a two-layer model [226] provides a very good
agreement with more refined predictions of the cooling evolution of white
dwarfs [227, 228]. Thus, in this work we consider an isothermal core with a
temperature in the range $4\times 10^{5}\leavevmode\nobreak\ \textrm{K}\leq
T_{c}\leq 4\times 10^{6}\leavevmode\nobreak\ \textrm{K}$, which is appropriate
for white dwarfs older than $\sim 3$ Gyr [228].

###  III.7 Neutron stars

In even heavier stars than those discussed above, $(8-10)\leavevmode\nobreak\
M_{\odot}\lesssim M\lesssim 25\leavevmode\nobreak\ M_{\odot}$, the carbon-
oxygen core reaches high enough temperatures ($T_{c}\gtrsim 10^{9}$ K) to
continue nuclear fusion processes before becoming degenerate, resulting in a
shell-like structure, which ends up in an iron core for $M\gtrsim
10\leavevmode\nobreak\ M_{\odot}$. Once nuclear fusion cannot proceed further,
stars approach the end of their lives as core-collapsed supernovae explosions
(or via electron-capture supernovae in a narrow range of masses). The remnants
of these explosions are neutron stars (see, e.g., Refs. [229, 214, 196]),
which are born very hot ($T\gtrsim 10^{10}\leavevmode\nobreak\ $ K), but cool
down very fast by neutrino emission (which lasts for $\lesssim 10^{5}$ yr),
down to $T\sim 10^{8}$ K after $\sim$100 yr. As density increases, neutron-
rich nuclei start releasing free neutrons, which are increasingly degenerate
and would eventually be the main source of pressure. Further increase of the
density would result in the formation of a degenerate bath of neutrons plus a
small admixture of electrons, muons and protons.

The description of the interior of neutron stars is particularly challenging.
Unlike main-sequence stars or other post-main-sequence phases, the huge
gravitational field of neutron stars requires the use of general relativity.
This modifies the Newtonian hydrostatic equilibrium equation, which is
replaced by the Tolman-Oppenheimer-Volkoff equation [230, 231]. This just
represents a calculational complication, but the most important issue is that
the correct equation of state to be used is not yet known. Nevertheless, the
determination of relatively large masses for neutron stars from gravitational
wave observations of merging systems of binary neutron stars and of a black
hole and a neutron star, favor stiff equations of state, which predict
relatively large maximum masses, $\sim(2-3)\leavevmode\nobreak\ M_{\odot}$
[232, 233, 234, 235, 236, 237, 238, 239, 240]. In this work, we consider the
unified equation of state from Ref. [222], based on the nuclear energy-density
functional BSk20 [223], which covers the mass range $0.09\leavevmode\nobreak\
M_{\odot}<M<2.3\leavevmode\nobreak\ M_{\odot}$. The resulting mass–radius
relation is shown in the right panel of Fig. 5. Notice that differences with
respect to more recent equations of state [241] have a negligible impact on
the results presented in this work.

For old neutron stars ($t_{\rm NS}\gtrsim 10^{5}$ yr), the main cooling
mechanism is electromagnetic cooling. At this evolutionary stage the internal
temperature is expected to be similar to the surface one and, in idealized
scenarios, temperatures as low as $\sim 10^{3}$ K after $\sim 10^{7}$ yr, and
even lower for older neutron stars, are expected [242]. Nevertheless, for
theroretical models to be consistent with observations, a heating mechanism is
required. Different possibilities have been suggested, which could result in
temperatures as high as $10^{6}$ K [243]. Therefore, we consider an isothermal
profile within the interval $10^{5}\leavevmode\nobreak\ \textrm{K}\leq
T_{c}\leq 10^{6}\leavevmode\nobreak\ \textrm{K}$, as representative of
relatively old neutron stars, $t_{\rm NS}>10^{7}$ yr. Neutron stars younger
than $\sim 10^{6}$ yr are expected to have typical central temperatures of
$\sim 10^{8}$ K [244].

##  IV DM evaporation mass in celestial bodies

After the overall description of celestial bodies in the previous section, we
now compute the minimum mass DM particles must have, such that evaporation
from the capturing object is not efficient. The calculation of all the
elements required for the computation of the DM evaporation mass follows Ref.
[66], including the correction to the capture rate from Ref. [67], and Ref.
[71] for neutron stars, as described in Section II. We first evaluate the DM
evaporation mass for the geometric cross section,
$\sum_{i}N_{i}\,\sigma_{i}^{\rm geom}=\pi\,R^{2}$, which results in a capture
rate close to its maximum value (the saturation value), for planetary objects,
brown dwarfs and main-sequence stars. This cross section is different for each
object and depends on the DM mass. Next, we compute the DM evaporation mass
for all those objects and within a wide range of cross sections, in the thin
and thick regimes. For post-main-sequence stellar phases we qualitatively
discuss the evolution of the DM evaporation mass, and for white dwarfs and
neutron stars we compute the DM evaporation mass for the geometric cross
section. All results are obtained for SI interactions with constant scattering
cross section, but we also comment on the DM evaporation mass in the SD case.
Moreover, we study its dependence with several factors.

###  IV.1 DM evaporation mass in planetary bodies, brown dwarfs and main-
sequence stars

Figure 6: DM equilibration time, for the DM evaporation mass, as a function of
the mass of the capturing object, for planetary bodies, brown dwarfs and main-
sequence stars. We take the geometric SI cross section,
$\sum_{i}N_{i}\,\sigma_{i}^{\rm geom}=\pi\,R^{2}$, and
$\langle\sigma_{A}v_{\chi\chi}\rangle=3\times 10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}$. Also depicted are the current solar age,
$t_{\odot}=4.5$ Gyr (black dashed line) and the stellar lifetime, $t_{\rm
life,MS}$, when shorter than $t_{\odot}$ (blue dashed line). The jump at
$M=2\leavevmode\nobreak\ M_{\oplus}$ is mainly due to the non-smooth
transition in composition and the fact that $v_{e}<v_{d}$.

Before discussing the main results, it is important to evaluate whether
equilibrium between DM capture and annihilation is reached for each particular
celestial body. Otherwise, the DM evaporation mass is not constant but grows
with time, assuming the properties of the capturing object remain
approximately unchanged. For very low-mass objects, the equilibration time can
be much longer than their age, and even than the age of the Universe. For
intermediate-mass and the most massive stars, the equilibration time (for the
DM evaporation mass) decreases with stellar mass, although with a much weaker
dependence than the lifetime.

In Fig. 6 we compare the equilibration time, $\tau_{\rm eq}$, corresponding to
the geometric cross section and to the DM evaporation mass shown in Fig. 7,
with other characteristic time scales (solar system current age,
$t_{\odot}\simeq 4.5$ Gyr, and stars lifetime). This comparison is shown as a
function of the mass of celestial bodies, which spans the range
$10^{-10}\leavevmode\nobreak\ M_{\odot}\leq M\leq 10^{2}\leavevmode\nobreak\
M_{\odot}$, from small satellites to massive stars. The general trend,
although it does not apply to the entire mass range, can be understood from
the discussion in Section II: the equilibration time for the DM evaporation
mass decreases with the mass of the capturing body. The discontinuity at
$M=2\leavevmode\nobreak\ M_{\oplus}$ can be understood from our modeling of a
non-smooth transition in composition and the fact that $v_{e}<v_{d}$, so the
$1/\mu_{-}^{2}$ factor plays a key role (see Eq. (II.5)). As already
mentioned, if the equilibration time is longer than the age of the object, the
DM evaporation mass grows with time, as the number of capture DM particles
does (as long as the properties of the object do not change). This occurs
until equilibration is reached or the system is destroyed or it evolves in a
significant manner.

For the local DM density, the canonical thermal annihilation cross section and
the geometric cross section, DM capture and annihilation do not reach
equilibrium for the smallest objects, $M\lesssim 3\times
10^{-8}\leavevmode\nobreak\ M_{\odot}$, during the solar system age
($t_{\odot}=4.5$ Gyr), so the DM evaporation mass would increase with time.
Likewise, equilibration does not take place during the lifetime of the most
massive stars, $M\gtrsim 60\leavevmode\nobreak\ M_{\odot}$. Furthermore,
planetary systems are usually bounded to stellar objects, so their lifetimes
could be linked to stellar lifetimes. Objects in the range $2\times
10^{-8}\leavevmode\nobreak\ M_{\odot}\lesssim M\lesssim 60\leavevmode\nobreak\
M_{\odot}$ are in the regime in which the DM capture and annihilation rates
could be in equilibrium, so the DM evaporation mass would remain approximately
constant, as long as the properties of the objects can be regarded as
constant. This is so for the DM annihilation and scattering cross sections and
the DM density and velocity dispersion we have considered in this figure.
Larger annihilation cross sections, higher DM densities or lower velocity
dispersion, would result in shorter equilibration times. For instance, for the
smallest objects or the most massive stars, if close to the galactic center,
equilibrium is more likely. Note that capture is already assumed to be close
to maximum, so the equilibration time could only increase by considering a
smaller scattering cross section and remains almost unchanged for larger cross
sections.

Figure 7: DM evaporation mass as a function of the mass of the capturing
object, for planetary bodies, brown dwarfs and main-sequence stars. We take
the geometric SI scattering cross section, $\sum_{i}N_{i}\,\sigma_{i}^{\rm
geom}=\pi\,R^{2}$, the canonical value of the DM annihilation cross section
$\langle\sigma_{A}v_{\chi\chi}\rangle=3\times 10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}$, and assume a position within the local
neighborhood, $\rho_{\chi}=0.4\leavevmode\nobreak\
\textrm{GeV}/\textrm{cm}^{3}$ and $v_{d}=270\leavevmode\nobreak\
\textrm{km/s}$, although with other values similar results are obtained. We
also indicate the DM evaporation mass using detailed models and data for the
Moon, Earth and Sun. The shaded band depicts the range $E_{c}/T_{c}=(20-40)$.

In Fig. 7, we show the DM evaporation mass as a function of the mass of the
capturing object, for the geometric SI cross section and assuming
$t=\textrm{min}\\{t_{\odot},\tau_{\rm life,MS}\\}$. We see that it decreases
from values $m_{\rm evap}\sim 400$ GeV for the smallest objects with spherical
shape that can attain hydrostatic equilibrium, $M\sim
10^{-10}\leavevmode\nobreak\ M_{\odot}$, to $m_{\rm evap}\simeq 0.7$ GeV for
super-Jupiters and small brown dwarfs, $M\sim 10^{-2}\leavevmode\nobreak\
M_{\odot}$. For more massive brown dwarfs, the DM evaporation mass is slightly
larger and grows when entering the stellar regime. Note, however, that for
$M\gtrsim M_{\oplus}$, the DM evaporation mass only varies within an order of
magnitude, $m_{\rm evap}\sim(1-10)$ GeV, and grows for smaller objects (rocky
planets and satellites) due to their small size. This behavior follows the
scaling $m_{\rm evap}\propto T_{\chi}R/(M\hat{\phi_{c}})$, with
$\hat{\phi}_{c}\equiv v_{e}^{2}(r=0)/v_{e}^{2}(r=R)$, which can be understood
from the fact that $E_{c}/T_{\chi}\sim 30$, as discussed in Section II. To
illustrate the robustness of this result, we show the range
$E_{c}/T_{c}=(20-40)$ with a band, which fully embeds the values of the DM
evaporation mass for all objects and roughly accounts for systematics in
modeling of celestial bodies properties. Furthermore, the small variation of
the DM evaporation mass can be understood by considering the virial theorem,
which implies that the factor $T_{c}\,R/M$ varies little for a given class of
objects. Similarly to the jump in the equilibration time at
$M=2\leavevmode\nobreak\ M_{\oplus}$, the discontinuity at that value on the
DM evaporation mass is caused by the abrupt transition in the composition and
density profile, from rocky planets to icy planets and to the fact that
$v_{e}<v_{d}$. For the considered parameters, equilibration is not reached for
$M\lesssim 3\times 10^{-8}\leavevmode\nobreak\ M_{\odot}$ (see Fig. 6) and
this explains the slight bending of the curve towards smaller DM evaporation
masses, as in those cases, the DM evaporation mass grows with time until
reaching equilibrium. The same occurs for $M\gtrsim 60\,M_{\odot}$, but in
those cases equilibrium cannot be reached, as it would require a time longer
than the age of those stars. Super-Earths, $M\lesssim 10\leavevmode\nobreak\
M_{\oplus}$, with a larger fraction of metals than what is assumed here, would
have a slightly larger DM evaporation mass. Note, however, that for a given
mass, heavier compositions generically imply smaller sizes [114, 115].

Remarkably, as evident from Eqs. (20) and (II.5), changes in parameters as
cross sections, DM density or velocity dispersion, affect the DM evaporation
mass only logarithmically. Therefore, given that $E_{c}/T_{\chi}\sim 30$, in
order to obtain a value of the DM evaporation mass smaller by a factor of two,
the term
$\log\left(\rho_{\chi}\langle\sigma_{A}v_{\chi\chi}\rangle/v_{d}\right)$ for
$v_{e}\gg v_{d}$, or
$\log\left(\rho_{\chi}\langle\sigma_{A}v_{\chi\chi}\rangle/v_{d}^{3}\right)$
for $v_{e}\ll v_{d}$, must be larger by a factor of the order of
$\sim(12-15)$. This implies that the sensitivity of the DM evaporation mass to
changes on these parameters is relatively weak and thus, its value is rather
stable against different particle physics models (with constant scattering
cross sections) or for different locations of celestial bodies within the host
galactic halo.

Figure 8: 2-D contour of the DM evaporation mass as a function of the SI
scattering cross section and the mass of the capturing objects, for planetary
bodies, brown dwarfs and main-sequence stars. We take the canonical value of
the DM annihilation cross section
$\langle\sigma_{A}v_{\chi\chi}\rangle=3\times 10^{-26}\leavevmode\nobreak\
\textrm{cm}^{3}/\textrm{s}$ and assume a position within the local
neighborhood, $\rho_{\chi}=0.4\leavevmode\nobreak\
\textrm{GeV}/\textrm{cm}^{3}$ and $v_{d}=270\leavevmode\nobreak\
\textrm{km/s}$, although with other values similar results are obtained. The
value of the SI geometric scattering cross section (dashed line),
$\sum_{i}N_{i}\,\sigma_{i}^{\rm geom}=\pi\,R^{2}$, is approximately
proportional to $R^{2}/M$, as expected.

Nevertheless, the results in Fig. 7 correspond to approximately the largest
possible value of the DM evaporation mass for each celestial body (aside from
variations on the average properties here considered), that is, the one
obtained for a capture rate close to the saturation value. Next, we also study
its variation with the value of the SI scattering cross section. This is
depicted in Fig. 8. The maximum value for all celestial bodies is visible in
the figure, which is indeed close to that obtained for the geometric cross
section (dashed line). For $M\gtrsim 2\times 10^{-7}\leavevmode\nobreak\
M_{\odot}$, due to the exponential suppression of the evaporation rate in the
thick regime, the smallest DM evaporation mass is achieved for the largest
value of the scattering cross section we consider,151515For larger cross
sections, the assumed scaling with the atomic mass is not entirely reliable
for contact interactions [245]. $\sigma_{p}=10^{-31}\,\textrm{cm}^{2}$. For
$M\lesssim 2\times 10^{-7}\leavevmode\nobreak\ M_{\odot}$, the transition
between the thin and thick regimes takes place at values of the cross section
closer to $\sigma_{p}=10^{-31}\,\textrm{cm}^{2}$, so the smallest DM
evaporation mass is found in the thin regime, for the smallest cross section
we consider, $\sigma_{p}=10^{-41}\,\textrm{cm}^{2}$. Note, however, that, for
the average properties of celestial bodies considered in this work, the DM
evaporation mass never reaches values below $\sim 250$ MeV in the parameter
space shown in Fig. 8. Therefore, even in extreme situations with
$\rho_{\chi}\langle\sigma_{A}v_{\chi\chi}\rangle/v_{d}$ for $v_{e}\gg v_{d}$
(or $\rho_{\chi}\langle\sigma_{A}v_{\chi\chi}\rangle/v_{d}^{3}$ for $v_{e}\ll
v_{d}$) being many orders of magnitude larger than what we have assumed for
this figure, DM evaporation masses much below that mass are very unlikely.

Along with the most massive stars, the objects for which the smallest values
are obtained are super-Jupiters and small brown dwarfs, as correctly pointed
out in Refs. [32, 76, 68]. Nevertheless, in those papers, the estimated values
of the DM evaporation mass down to a few MeV were obtained neglecting the
critical exponential tail of the evaporation rate and thus, are incorrect and
the conclusions reached for those very low masses are not valid; under the
assumptions in this paper, the DM evaporation mass is rather $m_{\rm
evap}>250$ MeV in the entire parameter space shown in Fig. 8. This can be
understood from two facts: the critical exponential tail of the DM evaporation
rate was incorrectly neglected, as explained in Section II, and the core
temperature of brown dwarfs was underestimated by a factor $\gtrsim 4$ (cf.
Refs. [187, 164, 179, 189, 199]). The importance of this tail, for the case of
the Sun, is not a new finding, but has been known for over three decades [21,
2, 3, 70]. We do stress that the fact that the DM evaporation mass is
approximately given by $E_{c}/T_{\chi}\sim 30$ can be generalized to all
spherical celestial bodies in hydrostatic equilibrium. Indeed, note that Ref.
[33] did correctly estimate the DM evaporation mass for the most massive brown
dwarfs and low-mass stars in the context of asymmetric DM scenarios.
Additionally, notice that we have considered temperatures of brown dwarfs
corresponding to late evolutionary stages. From this point of view, our
assumption is conservative, as younger brown dwarfs are warmer, which results
in a higher DM evaporation mass. Note that even older ($\gtrsim 10$ Gyr) brown
dwarfs would have a slightly lower core temperature ($\lesssim 20\%$ cooler
than at $5$ Gyr for the most massive ones and even more similar for the least
massive ones), which implies a correspondingly lower DM evaporation mass. On
another hand, let us stress again that uncertainties on the density or
temperature profiles have a small impact on the DM evaporation mass. All in
all, these uncertainties are approximately accounted for by the gray band in
Fig. 7.

As we have discussed above and can be seen from Fig. 8, the DM evaporation
mass for a given object mass is rather stable, within about a factor of two at
most, against variations by ten orders of magnitude in the scattering cross