# The effectiveness of exoplanets and Brown Dwarfs as sub-GeV Dark Matter
detectors

C. Ilie,11footnotetext: Corresponding author. C. Levy  and J. Diks

###### Abstract

In this work we demonstrate that Dark Matter (DM) evaporation severely hinders
the effectiveness of exoplanets and Brown Dwarfs as sub-GeV DM probes.
Moreover, we find useful analytic closed form approximations for DM capture
rates for arbitrary astrophysical objects, valid in four distinct regions in
the $\sigma-m_{X}$ parameter space. As expected, in one of those regions the
Dark Matter capture saturates to its geometric limit, i.e. the entire flux
crossing an object. As a consequence of this region, which for many objects
falls within the parameter space not excluded by direct detection experiments,
we point out the existence of a DM parameter dependent critical temperature
($T_{crit}$), above which astrophysical objects lose any sensitivity as Dark
Matter probes. For instance, Jupiters at the Galactic Center have a $T_{crit}$
ranging from $700$ K (for a $3M_{J}$ Jupiter) to $950$ K (for $14M_{J}$). This
limitation is rarely (if ever) considered in the previous literature of
indirect Dark Matter detection based on observable signatures of captured Dark
Matter inside celestial bodies.

##  1 Introduction

One of the most profound problems in modern physics is the that of the nature
of Dark Matter (DM). The origin of the term can be traced back to the work of
Fritz Zwicky in 1933 [1] who predicted the existence of a non-luminous
component of matter in the Coma Cluster of galaxies to explain the measured
velocity dispersion. The idea of DM remained dormant in the literature until
the 1970s when Rubin and Ford inferred its existence on galactic scales by
measuring the rotational speed of H2 gas surrounding the Andromeda Galaxy [2].
Since then, a plethora of evidence has emerged supporting the DM hypothesis,
including gravitational lensing [3], the bullet cluster [4], and the Cosmic
Microwave Background [5, 6, 7, 8], which predicts that DM composes $\sim 85\%$
of the mass-budget of the universe. Despite the clear indicators of DM’s
existence through its gravitational pull at all scales, the exact nature of DM
has remained elusive despite the combined efforts from theorists and
experimentalists. Various candidates for the form of DM have emerged since its
existence was inferred. These include massive compact halo objects [9] and
various forms of particle DM, like Axions [10] and Weakly Interacting Massive
Particles (WIMPs) [11], to name a few. Massive compact halo objects (MACHOs)
have largely been ruled out as dark matter candidates from experimental data
and theoretical arguments222See Sec. 4.1 of [12] for a detailed description on
the constraints on the model of massive compact halo objects., while Axions,
WIMPs, and other particle DM candidates have been heavily constrained by
direct detection experiments (e.g. ADMX [13], XENON1T [14], PICO60 [15], etc.)
which are designed to measure the effect of DM particles on detectors deep
underground.

Despite our lack of understanding DM’s true form, a clear picture has emerged
for the role of DM in the structure of our universe in the standard
cosmological model, $\Lambda$CDM. Overdensities of DM form in the early
universe after a period of rapid expansion, known as cosmic inflation, which
amplify quantum fluctuations in the early universe [16]. These dense regions
of DM, known as halos, provide the gravitational impetus for baryonic matter
to collapse and form the first stars [17]. Through mergers, these DM halos
grow, providing the seeds for galaxy formation [18]. Simulations of DM in the
early universe confirm this picture and also posit the existence of a
filament-like structure connecting the DM halos [19]. This hierarchical
structure formation predicts that most, if not all, galaxies are formed in DM
halos that have a rich sub-structure.

Dark Matter capture is a mechanism for accumulating particle DM within a
compact object [20, 21, 22, 23]. DM particles incident on an object first
become accelerated by the object’s gravitational field, then, as they transit
the object, may collide with its constituents, losing kinetic energy. Provided
the particle loses enough energy, it will become gravitationally bound to the
object. From there, a number of processes may occur that leaves a DM footprint
on the object, including DM-DM annihilation, and kinetic heat transfer,
whereby DM collisions with the object’s baryons adds energy to the object
[e.g. 24]. Capture, and its observational effects, have long been studied and
are of serious interest in the search for DM since many physical situations
consist of astrophysical objects within dense DM environments, which
intuitively produces higher rates of capture. A non-exhaustive list of
previous literature on the observational effects of captured dark matter
include the study of neutron stars [21, 22, 24, 25], Population III stars [26,
27, 28], white dwarfs [23, 29], and, most relevant to this work, exoplanets
[30].

Soon after the formalism for capture was developed [31, 20], it was realized
that the reverse process, evaporation, could be relevant for light DM
candidates [32]. After sufficient time, captured DM thermalizes in the object
and through interactions with the stellar medium has a probability of gaining
enough energy to escape the gravitational well of the object. In [30], the
authors consider self-annihilating models of DM and use the potential
observable temperature of exoplanets to generate sensitivity estimates for
exoplanets as DM probes. To do so, the authors cut their sensitivity limits at
$\sim 10^{-2}$ GeV, claiming that below this mass evaporation becomes relevant
enough to erase sensitivity. They support this claim through a calculation of
the “evaporation mass” which is defined to identify the mass scale at which
evaporation becomes relevant. However, as pointed out in [33], the traditional
way to calculate the evaporation mass neglects an important exponential tail
in the calculation, which means the evaporation mass calculated in [30] should
likely be higher. In this work we explicitly confirm this expectation via
numerical calculations.

One of the main goals of this manuscript is to explicitly compute the
evaporation rate for Jupiter-like exoplanets (Jupiters) and show how this
affects their potential as DM probes for sub-GeV DM. As a result, in Sec. 3.3
we produce sensitivity limits that account directly for evaporation and
explicitly show how the evaporation mass approximation ($\sim 10^{-2}$ GeV)
used in [30] is an underestimate. We find instead, that for the objects
considered (Jupiters and Brown Dwarfs placed at the Galactic Center) the
evaporation mass for is around $\sim 1$ GeV (see Figs. 5 and 6). As such,
exoplanets and Brown Dwarfs are quite poor DM probes in the sub-GeV regime.
Additionally, we point out an inherent limitation of any astrophysical probe
that uses the potential observation of an enhancement from an expected
background surface temperature as a signature of Dark Matter. Namely, the
existence of a critical temperature 333For a definition of the critical
temperature ($T_{crit}$) see Eq. 3.19. Its implications with regards to
exoplanets as DM probes are discussed in Sec. 4 above which all DM
constraining power is lost. This effect was pointed out by us before in [22],
and it is a result of the saturation of the capture rates in a certain region
of the $\sigma$ vs $m_{X}$ parameter space. Additionally, in Sec. 2.1 we
present simple analytic approximations of the DM capture rates that are valid
in four distinct regions of the same parameter space (see Eqns. 2.12-2.15 and
Fig. 1). We thus extend our findings of [27], where we considered the case of
Pop III stars, when a definite hierarchy between the escape velocity at the
surface of the star and the DM dispersion velocity exists
($v_{\text{esc}}\gg\bar{v}$).

In this work, we use polytropes to approximate the structure of Jupiters and
Brown Dwarfs as this information is necessary to compute the rate of DM
evaporation. The polytropic model assumes a pressure-density relation of the
form $P=K\rho^{1+\frac{1}{n}}$, where $n$ is the polytropic index. Polytropes
have long been used to estimate the structure of self-gravitating gaseous
bodies [34] and is thus a good candidate model for gaseous exoplanets. For
Jupiters, we follow [33] and take the index as $n=1$ which is justified by
comparing to expected models of Uranus and Neptune [35], Saturn [36] and
Jupiter [37]. For Brown Dwarfs we take index $n=1.5$, a choice justified in
detail by [38] or [33]. In Appendix A we include useful approximations of
integrals involving DM profiles in the limit of high DM mass ($m_{X}\gg
1\,\mathrm{GeV}$) in objects that are modelled by polytropes. These
approximations are very useful for computing the rate of DM annihilation in
objects for Supermassive DM since numerical integration procedures become
increasingly difficult for the highly cored profiles associated with heavy DM.

Perhaps one of the most intriguing and counter-intuitive findings of our work
is the fact that under certain conditions, dark matter has a “floating
distribution” towards the surface of the object (see Sec. 3.2, and in
particular Fig. 4). This confirms similar results independently found in [39]
and [40].

##  2 Dark Matter processes in exoplanets

In this section, we briefly review three key Dark Matter processes: capture,
evaporation, and annihilation. The interplay between those will controll the
effectiveness of any given astrophysical object as a Dark Matter probe. If we
label $C$ as the rate at which DM particles accumulate in the object through
down-scattering, $\Gamma_{A}$ the rate at which they annihilate within the
object, and $E$ the per-particle rate at which they leave the object through
up-scattering, then the differential equation governing the total number of DM
particles $N_{X}$ in the object at time $t$ is:

$\frac{dN_{X}}{dt}=C-\Gamma_{A}(N_{X})-EN_{X}.$ (2.1)

The annihilation rate, $\Gamma_{A}(N_{X})$, is a model-dependent quantity. In
this work, we primarly consider the “Co-SIMP” [41] sub-GeV model, where
annihilation is given by: DM$+$DM$+$SM$\rightarrow$ DM$+$SM. In making this
choice we follow [30], the original paper that pointed out the usefulness of
exoplanets as sub-GeV DM probes. Co-SIMPs are a sub-GeV DM model that evade
naturally the Lee-Weinberg [42] lower bound on a thermal relic DM particle
mass. We also consider general models that annihilate like
DM$+$DM$\rightarrow$ SM$+$SM for high-mass DM when we discuss the
approximation of integrals involving highly-cored DM profiles. Both models
obey $\Gamma_{A}\sim N_{X}^{2}$ since only two DM particles are required. By
defining $C_{A}$ as an $N_{X}$-independent annihilation coefficient, we can
write the annihilation rate as:

$\Gamma_{A}(N_{X})=C_{A}N_{X}^{2},$ (2.2)

With this in hand, a closed-form solution for Eq. 2.1 can be found:

$N_{X}(t)=\sqrt{\frac{C}{C_{A}}}\frac{\tanh\left(\frac{\kappa
t}{\tau_{eq}}\right)}{\kappa+\frac{1}{2}E\tau_{eq}\tanh\left(\frac{\kappa
t}{\tau_{eq}}\right)},$ (2.3)

where $\tau_{eq}\equiv\left(CC_{A}\right)^{-1/2}$ and
$\kappa\equiv\sqrt{1+\left(\frac{1}{2}E\tau_{eq}\right)^{2}}$. One can verify
that when $t\gg\frac{\tau_{eq}}{\kappa}$, equilibrium is achieved and:

$N_{X}\simeq\sqrt{\frac{C}{C_{A}}}\frac{1}{\kappa+\frac{1}{2}E\tau_{eq}}.$

In the following subsections, we outline in more details the processes of
capture, annihilation, and evaporation and how to compute them in the context
of exoplanets.

###  2.1 Capture

Dark Matter capture is the process through which DM particles incident on an
astrophysical object become gravitationally bound as they lose energy
transiting this object. First studied in the context of the sun [31] and the
Earth [20] for DM that becomes captured after a single scattering event, and
the Earth for DM that requires multiple scattering events [43], this process
has since gained a lot of attention within the astroparticle physics community
for its ability to constrain DM in regions of parameter space inaccessible to
Earth-based detection experiments. An important extension of Gould’s original
DM capture at finite optical depth formalism [43] came recently with the work
of Bramante [21],444Also note the comment paper [22] in which several errors
in the original [21] are addressed, and most importantly, the existence of
critical temperature that limits the usefulness of astrophysical objects as DM
probes whenever $T_{eff}\gtrsim T_{crit}$. Regarding this point see Eq. 3.19
and its discussion, and Sec. 4 which we use throughout this work, and outline
below. This extension provides a convenient way to compute the capture rate in
the multiscatter limit.

We also present analytic approximations of the capture rate and discuss the
form of the capture rate in various regions of $\sigma-m_{X}$ parameter space
(where $\sigma$ is the DM-proton scattering cross-section), expressions which
are derived in Appendix B. A significant new finding of this work is an
extension of the analytic capture rates according to regions presented in
[27], where the objects of interest were Pop III stars. As such, the limit of
large escape velocity ($v_{\text{esc}}$) relative to to the DM velocity
dispersion ($\bar{v}$) in the surrounding halo, i.e. the limit of
$v_{\text{esc}}\gg\bar{v}$, was considered. In this work, we relax this
assumption and explore the parameter space without any apriori assumption on a
hierarchy between $v_{\text{esc}}$ and $\bar{v}$. Our results when the escape
velocity is high agree exactly with the findings of [27]. In the limit of
$\bar{v}\gg v_{\text{esc}}$, we find a slightly different picture: one of the
regions becomes erased (Region III), and another region is pushed upwards
(Region II). We also find that, unlike when $v_{\text{esc}}\gg\bar{v}$, the
boundary of Region I and Region II no longer coincide and a new region between
them develops. More details on these new findings can be found in [44], though
we summarize them below and outline the key steps in Appendix B.

Capture in its schematic form can be represented by the following differential
form:

$\frac{dC}{dVd^{3}u}=\Omega(n_{T}(r),w(r),\sigma,m,m_{X})dF(n_{X},u,v_{\text{obj}},v_{\text{esc}}^{\text{halo}}),$
(2.4)

where $dC$ is the differential capture rate per volume $dV$ of the object and
velocity $d^{3}u$ of the DM particle far away, $\Omega$ is the probability of
capture, and $dF$ is the incident flux of DM on the object. Further,
$n_{T}(r)$ is the number density of target nuclei of mass $m$ at radius
r,555From now on, $n_{T}$ will represent the average number density of the
most dominant nuclear species throughout the object. For Jupiters, the most
abundant species is hydrogen, making up about $84\%$ of the object [45]. In
this work, we assume the object is made entirely of hydrogen for simplicity,
though in this case it makes no significant difference. For a full treatement
of the capture rate for an arbitrary number of components, see [28]. $w(r)$ is
the speed of the accelerated DM particles at radius r, $\sigma$ is the DM-
nucleon scattering cross section, $m_{X}$ is the mass of the DM particle,
$n_{X}$ is the number density of DM particles in the surrounding medium,
$v_{\text{obj}}$ is the velocity of the object with respect to the DM halo,
and $v_{\text{esc}}^{\text{halo}}$ is the escape velocity of the DM halo.
Computing the capture rate for a spherical object then, in principle, amounts
to integrating over the capture rate per shell of the object, integrating over
the flux of dark matter particles, and integrating over all DM velocities far
from the object. In practice, a number of simplifying assumptions are made to
make the problem more digestible, all of which can be shown as valid in the
context of exoplanets, as outlined in [30].

Following [21, 22], we compute the capture rate, $C$ in the following way:

$C=\sum_{N=1}^{\infty}C_{N},$ (2.5)

where $C_{N}$ is the so-called “partial capture rate,” i.e. the capture rate
for exactly $N$ scattering events. With the assumptions of [21], the partial
capture rate is given by:

$C_{N}=\underbrace{\pi R^{2}}_{\textrm{capture
area}}\times\,\underbrace{n_{X}\int_{0}^{\infty}\dfrac{f(u)du}{u}\,(u^{2}+v_{\text{esc}}^{2})}_{\textrm{DM
flux}}\times\,\underbrace{p_{N}(\tau)}_{\textrm{prob. for $N$
collisions}}\times\,\underbrace{g_{N}(u)}_{\textrm{prob. of capture}},$ (2.6)

where $R$ is the radius of the object, $u$ is the speed of the DM particle far
from the object, $f(u)$ is the DM speed distribution in the medium,
$v_{\text{esc}}\equiv\sqrt{\frac{2GM}{R}}$ is the surface escape velocity of
the object, $\tau\equiv 2R\sigma n_{T}$ is the optical depth,666Physically,
this represents the average number of scatters a DM particle would undergo if
it traversed the diameter of the object. $p_{N}(\tau)$ is the probability a DM
particle undergoes $N$ scatters as it traverses the object, and $g_{N}(u)$ is
the probability a DM particle with initial speed $u$ becomes gravitationally
bound having scattered $N$ times. One can show that for $p_{N}(\tau)$ we have
[22]:

$p_{N}(\tau)=\frac{2}{\tau^{2}}\left(N+1-\frac{\Gamma(N+2,\tau)}{N!}\right).$
(2.7)

Further progress can be made on this expression if we assume a particular form
of $f(u)$. As is standard for thermal distributions, which is expected for DM
particles in the halos hosting exoplanets, $f(u)$ is taken to be a Maxwell-
Boltzmann distribution, and the partial capture rate is thus:

$C_{N}=\frac{1}{3}\pi
R^{2}p_{N}(\tau)\frac{\sqrt{6}n_{X}}{\sqrt{\pi}\bar{v}}\left(\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)-\left(2\bar{v}^{2}+3v_{N}^{2}\right)\exp\left(-\frac{3\left(v_{N}^{2}-v_{\text{esc}}^{2}\right)}{2\bar{v}^{2}}\right)\right),$
(2.8)

where $\bar{v}$ is the DM velocity dispersion in the capturing region of the
halo, and:

$v_{N}^{2}\equiv v_{\text{esc}}^{2}\left(1-\langle
z\rangle\beta_{+}\right)^{-N},$ (2.9)

where $\langle z\rangle\in[0,1]$ is the average of a kinematic parameter $z$,
defined as follows: a DM particle incident on the object will have an initial
kinetic energy $E_{0}=\frac{1}{2}m_{X}w^{2}$. After elastically colliding with
a nucleus, the particle will lose an amount of energy $\Delta
E=z\beta_{+}E_{0}$, where $z$ is related the scattering angle in the center of
mass frame by $z=\sin^{2}(\theta_{\text{CM}}/2)$ [23], and
$\beta_{+}\equiv\frac{4mm_{X}}{(m+m_{X})^{2}}$.

In [22, 27], analytic expressions for the total capture rate were found and
verified for the limiting regime of $v_{\text{esc}}\gg\bar{v}$ and $\langle
z\rangle\sim\frac{1}{2}$. The result, nicely summarized by Fig. 4 of [27], is
a set of analytic expressions for different regions of $\sigma-m_{X}$
parameter space, divided according to the optical depth, $\tau$, and the
dimensionless kinematic variable
$k\equiv\frac{3v_{\text{esc}}^{2}}{\bar{v}^{2}}\frac{\min(m_{X},m)}{\max(m_{X},m)}$.
In this work, we apply a more general analysis and acquire analytic capture
rates regardless of $\bar{v}-v_{\text{esc}}$ hierarchy, and for all possible
values of $\langle z\rangle$. This was motivated by the fact that many of the
exoplanet models considered stradle the boundary of $\bar{v}\sim
v_{\text{esc}}$, so it was necessary to have a clear picture of the analytic
capture rates in this more general scenario.

We summarize the results below, and derive them fully in Appendix B. First, it
is convenient to define a couple of parameters:

$\alpha\equiv\frac{3v_{\text{esc}}^{2}}{2\bar{v}^{2}}\langle z\rangle,$ (2.10)

and

$k\equiv\alpha\beta_{+}.$ (2.11)

Then, two distinct pictures emerge when $\alpha\geq 1$ and $\alpha<1$. The
case of $\alpha\gg 1$ was addressed in [27] where we found that the parameter
space naturally divided into 4 regions according to the values of $\tau$, $k$,
and $k\times\tau$, when compared to unity. These expressions are valid not
only in the $\alpha\geq 1$ regime, as a very similar, albeit slightly
different, picture emerges in the $\alpha<1$ regime. The four regions found in
[27] thus completely describe the analytic capture rates in most of the
$\sigma-m_{X}$ parameter space, though their distributions throughout the
space depend on $\alpha$. The analytic capture rates in the four regions
identified by us (see Fig. 1) are:

$C=2Ak\tau v_{\text{esc}}^{2},\quad\quad\text{(Region I)},$ (2.12)
$C=A\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right),\quad\text{(Region II)},$
(2.13) $C=2A\tau v_{\text{esc}}^{2},\quad\quad\text{(Region III)},$ (2.14)
$C=2Ak\tau v_{\text{esc}}^{2},\quad\quad\text{(Region IV)},$ (2.15)

where $A\equiv\frac{1}{3}\pi R^{2}\sqrt{\frac{6}{\pi}}\frac{n_{X}}{\bar{v}}$.
The simple capture rates in these regions stem from various approximations,
outlined in Appendix B. The key differences in the $\alpha<1$ regime, as one
can see from the schematic depiction in Fig. 1, is the erasure of Region III,
the decoupling of the Region I and Region II boundary, and the emergence of a
region between Region I and Region II where either approximation breaks down.
When $\alpha\geq 1$, the boundary of Region II, which coincides with the
boundary of Region I, is well-defined by the $k\tau=1$ line, a fact which no
longer holds in the $\alpha<1$ case. The new boundary for Region I is given by
the line $\frac{\langle z\rangle k\tau}{\alpha}=1$. However, as of now we
could not find an analytic form for the boundary of Region II in this
scenario, though we confirmed numerically that it lies between $k\tau=1$ and
$\frac{\langle z\rangle k\tau}{\alpha}=1$. Another very intriguing feature of
these equations, which is evidenced in Fig. 1 and Fig. 2 by the smooth
transition between Region I and IV, is the degeneracy of the capture rate in
Regions I and IV. This is a rather counter-intuitive result since, in deriving
each result, we assume opposite limits for the optical depth, $\tau$, i.e.
$\tau\gg 1$ for Region I , and $\tau\ll 1$ for Region IV.

In Fig. 1, we plot the $\sigma-m_{X}$ parameter space schematically,
indicating the regions identified in Eqs. (2.12)-(2.15). In this diagram, we
see the role that the $\bar{v}-v_{\text{esc}}$ hierarchy plays in determining
the location of the regions, and therefore where the different analytic
capture rates are valid. When $\alpha>1$, we see the already established
picture emerge as shown in [27]: four regions of parameter space, divided
according to the boundaries of $\tau=1$, $k=1$, and $k\tau=1$. When
$\alpha\leq 1$, we see the disappearance of Region III since $k\leq 1$ always
in this limit. We also begin to see the splitting of boundaries between Region
I and Region II, as the Region I approximations break down before the capture
rate is saturated. The new boundary for RI is determined by the line $\langle
z\rangle k\tau/\alpha=1$ as is shown in Appendix B. In between these
boundaries we have yet to produce a simplified analytic useful closed form for
the capture rate, though one can still compute it numerically. However, we
point out that in [44] we have obtained a full closed form (albeit very
cumbersome) of the total capture rates covering the entirety of the $\sigma-
m_{X}$ parameter space, for any $v_{\text{esc}}$ and $\bar{v}$. Lastly, we
point out that the boundary of Region II also appears lower than the $k\tau=1$
line, another difference when in the $\bar{v}>v_{\text{esc}}$ regime.

Figure 1: Schematic depiction of $\sigma-m_{X}$ parameter space indicating the
regions associated with the analytic capture rates presented in Eqs.
(2.12)-(2.15). The first panel describes the results found in [27] where
$v_{\text{esc}}>\bar{v}$, which corresponds to $\alpha>1$. When $\alpha\leq
1$, we see the erasure of Region III as seen in the central and right panels.
In the right panel, when $\alpha<1$, we see an interesting phenomenon occur:
the boundary between Region I (blue dashed line) and Region II (red dashed
line) become disassociated. The reason for this is that the approximation
leading to the emergence of Region I breaks down before the capture rate
becomes saturated in Region II. In the white region of the right panel,
between the boundaries of Region I and Region II, the capture rate can be
calculated numerically but we are unable to find a simple closed-form solution
(although in [44] we found an analytic–albeit very cumbersome–closed form for
the entire $\sigma-m_{X}$ parameter space). Intriguingly, the boundary of
Region II is also shifted below the $k\tau=1$ line, though finding its exact
location can only be done numerically at this point. However, we can say for
certain that it occurs somewhere between the $k\tau=1$ line and the Region I
boundary, given now by the following condition: $k\tau\frac{\langle
z\rangle}{\alpha}=1$.

In Fig. 2, we validate our analytic approximations for the DM capture rates in
the four distinct regions (RI-RIV) identified. Specifically, we plot the
relative error in the analytic capture rates of Eqs. (2.12)-(2.15) as compared
to the numerical capture rate (Eq. 2.5 combined with Eq. 2.8). Each panel
corresponds to a different value of $\alpha$, as labeled in the plot. In this
figure, we see the incredible accuracy in the approximations of the capture
rate, as well as the emergence of the regions depicted schematically in Fig.
1. Also, as mentioned above, we clearly see the degeneracy of our analytic
estimates in comparing Regions I and IV which, remarkably, produce the exact
same form for the analytic capture rate despite requiring different
assumptions in their derivations. Some key features that emerge when
$\alpha\leq 1$ are the erasure of a Region III (which is clearly present when
$\alpha>1$) and the failure to approximate the capture rate well between
Region I and Region II. As discussed above in reference to Fig. 1, we find
that the transition out of Region I and the transition into Region II no
longer coincide when $\alpha<1$. More details on this are found in Appendix B.

Figure 2: Relative error in the analytic approximations of the total capture
rate given by Eqs. (2.12)-(2.15) as compared to the numerical capture rate in
Eq. 2.5 and Eq. 2.8. Parameters are chosen conveniently to represent different
values of $\alpha$. In all panels, the solid red line represents $k\tau=1$,
while the dashed red lines represent $k\tau\frac{\langle z\rangle}{\alpha}=1$.
In the top panels, when $\alpha<1$, the gray region represents the decoupling
between regions I and II. We clearly see that the boundary of Region II, which
is the upper boundary of the grey region, lies below the $k\tau=1$ line and
above the Region I boundary, which is the lower boundary of the grey region.
This contrasts the lower panels, when $\alpha\geq 1$, where the boundaries of
Region I and II coincide at the solid red line, i.e. $k\tau=1$.

###  2.2 Evaporation

Evaporation is the reverse process of capture. DM particles that are
gravitationally bound will interact with the object’s constituents and,
through these interactions, it is possible for a particle to gain enough
energy to escape the object. Like capture, evaporation was first studied in
the context of the Sun, and has since been applied to a variety of
astrophysical objects. The rate of evaporation per particle in an object is
given by the following prescription [33]:

$E=\frac{1}{N_{X}}\int_{0}^{R}s_{\text{evap}}(r)n_{X}(r)~{}4\pi
r^{2}~{}dr\int_{0}^{v_{\text{esc}}(r)}f_{X}(\boldsymbol{w},r)~{}4\pi
w^{2}~{}dw\int_{v_{\text{esc}}(r)}^{\infty}\mathcal{R}^{+}(w\rightarrow
v)~{}dv,$ (2.16)

where $N_{X}$ is the total number of captured DM particles, $v$ represents the
up-scattered speeds sufficient for escape, $v_{\text{esc}}(r)$ is the escape
velocity at radius $r\leq R$, $\boldsymbol{w}$ represents the velocity of DM
particles with local distribution $f_{X}(\boldsymbol{w},r)$ (taken to be a
Maxwell-Boltzmann distribution in this paper), $\mathcal{R}^{+}(w\rightarrow
v)$ is the rate of up-scattering from speed $w$ to $v$, $n_{X}(r)$ is the DM
number density distribution within the object, and $s_{\text{evap}}(r)$ is a
suppression factor which accounts for two distinct possibilities: 1. DM
leaving in non-radial trajectories, and 2. DM down-scattering while on the way
out. The suppression factor is given by [33]:

$s(r)=e^{-\tau(r)}\underbrace{\frac{7}{10}\frac{1-e^{-\frac{10}{7}\tau(r)}}{\tau(r)}}_{\text{Angular
Factor}}\underbrace{{}_{0}F_{1}\left(;1+\frac{2}{3}\hat{\phi}(r);\tau(r)\right)}_{\text{Down-
scatter Factor}},$ (2.17)

where ${}_{0}F_{1}\left(;b;z\right)$ is the confluent hypergeometric limit
function and:

$\tau(r)\equiv\int_{r}^{R}\sigma
n(r^{\prime})dr^{\prime},~{}\hat{\phi}(r)\equiv\frac{\frac{1}{2}m_{X}v_{\text{esc}}^{2}(r)}{T(r)}$

Assuming a velocity-independent scattering cross section, the rate of
upscattering can be shown to carry the following functional dependence:

$\mathcal{R}^{+}(w\rightarrow
v)=\mathcal{R}^{+}(w,~{}v,~{}m,~{}m_{X},~{}T(r),~{}n(r),~{}\sigma),$ (2.18)

where $n(r)$ is the hydrogen number density distribution and $T(r)$ is its
temperature distribution. More details can be found in Appendix A of [46] for
the full form as well as details on how to arrive there.

To compute the evaporation rate, then, a number of features of the object and
how DM distributes itself in the object are required. In the case where DM
interacts very little with the object’s constituents, the isothermal limit,
its distribution is given by a rather straightforward isothermal distribution.
However, if there are a large number of interactions between DM and the
object’s constituents, the local thermal equilibrium (LTE) limit,
understanding the object’s distribution properties is essential to
understanding the DM distribution. In Section 3.2, we discuss how to approach
this problem for Jupiters and Brown Dwarfs, outlining how polytropes with
index $n=1-1.5$ can be used to model these objects and how that information is
further employed to determine the DM distribution in the high-interaction
limit. Typically, one finds themselves in a regime of either very high or very
low interaction and in these cases, computing the evaporation rate is
relatively simple, since finding the profile is straightforward. However, when
straddling the border between isothermal and LTE, some care must be taken. To
tackle this problem, we follow [46] and interpolate between the two regimes
via a parameter called the Knudsen Number. It is defined as:

$\text{Kn}=\frac{l_{\text{mfp}}}{L},$ (2.19)

where $l_{\text{mfp}}\equiv(\sigma n_{T})^{-1}$ is the mean free path, $n_{T}$
is the average number density of hydrogen, and $L$ is the length scale of the
system. Thus, when $\text{Kn}\gg 1$, the so-called “Knudsen limit,” DM
interacts very little with the surrounding medium since $l_{\text{mfp}}\gg L$.
Conversely, when $\text{Kn}\ll 1$, there is a high level of interaction.
Determining the length scale of the system can be done by finding the radius
$r_{X}$ at which the DM particles are in thermodynamic equilibrium with the
core [31]. In other words, solving the following:

$\frac{3}{2}T_{c}=m_{X}\Phi(r_{X}),$ (2.20)

where $\Phi(r)$ is the gravitational potential of the object. The solution is
given by [31]:

$L=r_{X}=\left(\frac{9}{4\pi}\frac{T_{c}}{G\rho_{c}m}\right)^{1/2}\sqrt{\frac{m}{m_{X}}}.$
(2.21)

Following [46], we define the Knudsen transition function by:

$f(\text{Kn})=\frac{1}{1+\left(\text{Kn}/0.4\right)^{2}},$ (2.22)

which is the interpolation function between the isothermal and the LTE limits.
Note the following, which is consistent with this interpretation of $f(K)$ as
an interpolation:

$\lim_{\text{Kn}\to\infty}f(K)=0,$ $\lim_{\text{Kn}\to 0}f(K)=1.$

With this, [46] defines the product $n_{X}(r)f_{X}(\mathbf{w},r)$ in the
intermediary regime ($\text{Kn}\sim 1$) to be:777If one is interested in a
more accurate approach to the Knudsen transition regime, see [47], where Monte
Carlo simulations were carried out to model energy transport of DM in
astrophysical bodies.

$n_{X}(r)f_{X}(\mathbf{w},r)=f(\text{Kn})n_{X}^{\text{LTE}}f_{X}(\mathbf{w},r)^{\text{LTE}}+\left[1-f(\text{Kn})\right]n_{X}^{\text{ISO}}f_{X}(\mathbf{w},r)^{\text{ISO}}.$
(2.23)

By integrating over the DM velocity space ($\mathbf{w}$), one gets:

$n_{X}(r)=f(\text{Kn})n_{X}^{\text{LTE}}+\left[1-f(\text{Kn})\right]n_{X}^{\text{ISO}}.$
(2.24)

Combining Eqs. 2.16 and 2.23 and noting that $f(K)$ has no dependence on the
variables being integrated, we see that we can generally write the evaporation
rate as:

$E=f(\text{Kn})E_{\text{LTE}}+\left[1-f(\text{Kn})\right]E_{\text{ISO}},$
(2.25)

where $E_{\text{LTE}}$ and $E_{\text{ISO}}$ are given by Eq. 2.16 where the
distributions $n_{X}$ and $f_{X}$ are taken in the LTE/isothermal limits. Note
here that, in both the LTE and isothermal case, we do not need to actually
know the total number of DM particles in the object, $N_{X}$, despite its
appearance in Eq. 2.16. This is because we have:

$E\sim\frac{\int
n_{X}(r)g(r)dV}{N_{X}}\sim\frac{n_{X}(r=0)\int\frac{n_{X}(r)}{n_{X}(r=0)}g(r)dV}{n_{X}(r=0)\int\frac{n_{X}}{n_{X}(r=0)}dV}\sim\frac{\int\frac{n_{X}}{n_{X}(r=0)}g(r)dV}{\int\frac{n_{X}}{n_{X}(r=0)}dV},$

where $g(r)$ represents the rest of whats being integrated in Eq. 2.16, and
$n_{X}/n_{X}(r=0)$ is the dimensionless DM distribution, whose form we know
exactly in the isothermal and LTE regimes, as discussed in Section 3.2. One
must take care when applying the interpolating prescription in Eq. 2.25,
however, since you may enter certain regimes where it is no longer valid. This
is as a result of the exponential cutoff of the evaporation rate above the so-
called “evaporation mass.” For example, if $\text{Kn}\ll 1$ so that we should
have $E=E_{\text{LTE}}$ but whenever, for a given set of parameters,
$E_{\text{LTE}}$ is exponentially suppressed more than $E_{\text{ISO}}$, one
will find $E\approx\left[1-f(\text{Kn})\right]E_{\text{ISO}}$. This is because
the first term would be more exponentially suppressed than the second term and
the $\left[1-f(\text{Kn})\right]$ factor does not send $E_{\text{ISO}}$ to $0$
quickly enough. To avoid this issue, we stipulate a condition when computing
the evaporation rate that when one is well into the $\text{Kn}\ll 1$ or
$\text{Kn}\gg 1$ regime, we switch to $E=E_{\text{LTE}}$ or $E_{\text{ISO}}$
respectively.

###  2.3 Annihilation

Once DM is captured, if it is self-annihilating it could contribute a heat
source to the object through model-dependent annihilation channels. As
mentioned in the beginning of the section, we consider the CoSIMP DM model for
the purpose of obtaining sensitivity limits in the sub-GeV mass region, where,
as suggested by [30], exoplanets should excel as DM detectors. For future
reference we consider here other $2\to 2$ Dark Matter models, including WIMPs
and Superheavy DM models, in approximating annihilation coefficients for heavy
($m_{X}\gg 1\,\mathrm{GeV}$) DM. The CoSIMP model assumes a strong coupling
between DM and SM particles, and the thermally averaged annihilation cross
sections, denoted by $\langle\sigma v^{2}\rangle$, are constrained to produce
the correct relic abundance of DM in the universe [41]. We focus on CoSIMP
annihilation over SIMP annihilation since, while these number changing ($3\to
2$) process can happen simultaneously [48], in exoplanets the CoSIMP model
dominates the annihilation process due to the high density of SM particles in
exoplanets that aid the efficiency of annihilation [30]. In natural units, the
annihilation interaction cross-section that gives the correct relic abundance
is [30]:

$\langle\sigma
v^{2}\rangle=10^{3}\left(\frac{1\,\mathrm{GeV}}{m_{X}}\right)^{3}\,\mathrm{GeV}^{-5}.$
(2.26)

For $2\to 2$ models outside of the WIMP window, one could consider thermally
averaged cross sections $\langle\sigma v\rangle$ bounded by the Griest-
Kamionkowski bound [49]:

$\langle\sigma v\rangle=\frac{4\pi}{m_{X}^{2}v_{X}},$ (2.27)

or, in the WIMP regime, the “WIMP miracle” thermally averaged cross section:

$\langle\sigma v\rangle_{\text{WIMP}}\simeq 10^{-23}\,\mathrm{cm^{3}s^{-1}}.$
(2.28)

With this in mind, the $N_{X}$-independent annihilation coefficient for the
CoSIMP model is [30]:

$C_{A}=\frac{\int dV~{}n_{X}^{2}n\langle\sigma v^{2}\rangle}{\left(\int
dVn_{X}\right)^{2}},$ (2.29)

while for the $2\to 2$ models it is:

$C_{A}=\frac{\int dV~{}n_{X}^{2}\langle\sigma v\rangle}{\left(\int
dVn_{X}\right)^{2}}.$ (2.30)

If we take Eq. 2.29, for example, and rewrite it in the following way:

$\Gamma=C_{A}N_{X}^{2}=\int dV~{}n_{X}^{2}n\langle\sigma v^{2}\rangle,$ (2.31)

we see that it is just an explicit way to calculate the total annihilation
rate by summing the local annihilation rates $d\Gamma=\langle\sigma
v^{2}\rangle n_{X}^{2}n~{}dV$ over the entire object. Note that, given a DM
distribution of the form $n_{X}(r)=n_{X}(0)\eta(r)$, for some $\eta(r)$ such
that $\eta(0)=1$, that the dimensional component of the distribution cancels
between numerator and denominator of Eq. 2.29, so only functional form
$\eta(r)$ of the distribution is necessary. This is the benefit of considering
an annihilation coefficient $C_{A}$: it gives an insight into the annihilation
efficiency of DM in the object without needing to know the number of DM
particles present. Furthermore, it is the annihilation coefficient, along with
the capture rate and the per-particle evaporation rate, that determine the
time-scale at which equilibrium between these processes occur, as per Eq. 2.1.
In section 3.2 we find the functional form of $\eta(r)$ and in Appendix A we
use this result to analytically approximate $C_{A}$ in the high DM mass limit
for $2\to 2$ DM models.

##  3 Constraining Dark Matter with exoplanets

###  3.1 Background on exoplanets and exoplanet searches

When considering the promise of astrophysical objects as probes of DM through
the processes of capture and subsequent annihilation, there are a number of
desirable features that will tend to make a DM signal stronger. Ideally, one
would like an object with a high baryonic density, large surface area, small
internal heat sources, and in a dense DM environment. However, many times only
one or a few of the conditions listed above are met. For instance, as the
universe expands, the density of DM in the universe dilutes, so objects in
early DM halos, such as Population III stars, have the advantage of being in
very dense DM environments. From an observational point of view, however, such
objects pose a challenge due to their large distance. On the other hand,
exoplanets in the Milky Way can provide nearby astrophysical objects with
little to no internal heat sources, and as such they can be used as probes of
DM, as shown initially by [30].

Exoplanets are planetary bodies lying outside of our solar system. These
bodies can take a wide variety of forms, and may be broadly categorized into
“Earth-like” planets, Jupiters, and Brown Dwarfs based on their size and make-
up. Earth-like planets are characterized by their relatively small size and
rocky interiors. Jupiters are composed of gas, predominantly hydrogen, and
have radii of approximately Jupiter’s radius and range in mass from the mass
of Jupiter to about 14 times the mass of Jupiter. Brown Dwarfs are the largest
of the three ranging in mass from about $14M_{J}$ to $75M_{J}$ and straddle
the planet-star boundary. Remarkably, Brown Dwarfs have approximately the same
radii as Jupiters, making them much denser objects and ideal as DM probes.
However, with this increased density comes additional internal heating and
longer cooling timescales that may pose a challenge in identifying a DM
signal. Because of their larger size and density, we focus on Jupiters and
Brown Dwarfs as DM probes in this work. To give a sense of exoplanet abundance
in our own Milky Way, note that microlensing observations have been used to
estimate that every star in our galaxy has, on average, at least one bound
exoplanet [50], and the number of stars in the Milky Way is of order $\sim
100$ billion. Furthermore, [50] find that about $20\%$ of stars in the Milky
Way have a Jupiter to Brown Dwarf sized companion planet. Exoplanets are ideal
probes of DM primarily due to their proximity. Searches targeting exoplanets
within our own galaxy of course have the advantage of location. To date, there
have been a number of exoplanet detections within our own galaxy [51], with a
few recent hot gas giants being detected in one single month alone [52, 53],
making these promising candidates as probes of DM. However, the candidate HIP
99770 b is at most 414Myr in age [52], which accounts for its relatively high
temperature. Although the candidate TOI 4727 b in [53] is more than 1Gyr in
age, its temperature is typical for exoplanets at its orbital distance
assuming radiative properties similar to Jupiter. While those two examples
(TOI 4727 b HIP 99770 b) of hot gas giants did not prove to be particularly
useful as DM probes, the promise of future detections of rogue overheated
exoplanets with dedicated exoplanet surveys is high. Moreover, with its
excellent field of view, the upcoming Roman Space Telescope [54] adds to the
prospects of discovering numerous objects that can serve as Dark Matter
probes, such as those considered in this paper: rogue overheated Jupiters and
Brown Dwarfs.

###  3.2 Particle distributions in exoplanets

In this section, we outline how we go about approximating the hydrogen and DM
distribution in Jupiters and Brown Dwarfs. We start by describing
distributions of hydrogen before discussing how we find the DM distributions
in these objects. In Appendix A, we use the results of this section to
approximate the annihilation coefficient for highly-cored profiles. It is
worth emphasizing here that, as done in [30] and the majority of works
studying the astrophysical observables of annihilating DM, we are making the
simplifying assumption that the effects of DM heating on the object does not
change its structure. Going beyond this commonly used approximation would
require incorporating the effects of DM on the heating and energy transport in
stellar evolution codes, such as MESA [55].

Jupiters and Brown Dwarfs are made predominantly of hydrogen and, to a good
approximation, can be modelled by polytropes of with an index ranging from
$n=1$ to $n=\frac{3}{2}$. We follow [33] in our choices of polytrope index,
who chose $n=1$ for Jupiters by comparing to the expected models of gas giants
in our solar system [35, 36, 37] and $n=3/2$ for Brown Dwarfs, as they are
dominated by convective energy transport [38]. The $n=1$ polytrope has a
closed-form solution while the $n=1.5$ must be solved numerically. Polytropic
distributions of index $n$ are of the form
$\rho(\xi)=\rho_{c}\theta^{n}(\xi)$, where $\xi$ is a dimensionless radial
variable and $\theta(\xi)$ is a solution to the Lane-Emden equation:

$\frac{1}{\xi^{2}}\frac{d}{d\xi}\left(\xi^{2}\frac{d\theta}{d\xi}\right)=-\theta^{n}.$
(3.1)

The above equation is a restatement of the condition of hydrostatic
equilibrium of a self-gravitating, spherically symmetric fluid. The boundary
conditions are given by: $\theta(0)=1$ and $\theta^{\prime}(0)=0$. The
dimensionless radial parameter is defined by $\xi=r\frac{\xi_{1}}{R}$, where
$\xi_{1}$ is the first positive solution of $\theta(\xi)=0$. Thus,
$\xi=\xi_{1}\iff r=R$. If we assume the hydrogen behaves like an ideal gas, as
is the case for Jupiters, the temperature profile of the object is given by:
$T(\xi)=T_{c}\theta(\xi)$. Remarkably, the same relationship is also valid the
interior of Brown Dwarfs [56].

Furthermore, one can show that for a general polytrope, the mass enclosed by a
radius $\xi$ is:

$M(\xi)=-4\pi\rho_{c}\left(\frac{R}{\xi_{1}}\right)^{3}\xi^{2}\frac{d\theta}{d\xi}.$
(3.2)

If the minus sign is concerning, note that $\theta^{\prime}(\xi)\leq 0$ for
all $\xi\in\left[0,\xi_{1}\right]$. The gravitational potential is defined in
the usual way:

$\Phi(r)\equiv\int_{0}^{r}dr^{\prime}\frac{GM(r^{\prime})}{r^{\prime 2}}.$
(3.3)

Using the dimensionless variable $\xi$, as appropriate for polytropes, this
becomes:

$\Phi(\xi)\equiv\frac{\xi_{1}}{R}\int_{0}^{\xi}d\xi^{\prime}\frac{GM(\xi^{\prime})}{\xi^{\prime
2}}.$ (3.4)

Then, we can compute this in terms of the solution to the Lane-Emden equation
as:

$\Phi(\xi)=\frac{1}{2}\frac{\xi_{1}}{\alpha_{n}}v_{esc}^{2}[1-\theta(\xi)],$
(3.5)

where
$\alpha_{n}\equiv\xi_{1}^{2}\left(\frac{d\theta}{d\xi}\Bigr{|}_{\xi_{1}}\right)^{2}$.
We note here that for the case of n=1, $\xi_{1}=\alpha_{n}=\pi$. As we shall
see later on, having an expression for the gravitational potential in terms of
the Lane-Emden function will help us tremendously in computing $n_{X}(r)$ in
the limit of high interactions.

For index $n=1$, Eq. 3.1, along with the boundary conditions, has the
following closed form solution:

$\theta(\xi)=\frac{\sin(\xi)}{\xi}$ (3.6)

Clearly, then, $\xi_{1}=\pi$ for $n=1$. One can approximate $\theta(\xi)$
quite straightforwardly using Taylor expansion around $\xi=0$:

$\theta(\xi)\simeq
1-\frac{1}{6}\xi^{2}+\frac{1}{120}\xi^{4}-\frac{1}{5040}\xi^{6}+\mathcal{O}(\xi^{8}).$
(3.7)

Note, however, that in this approximation, the first zero of $\theta(\xi)$ is
at $\xi_{1}\approx 3.08$, which one must keep in mind if integrating radially.
In Appendix A, which focuses on approximating integrals of the form:

$\int dVn_{X}(\xi),$

for very heavy DM, we will utilize this expansion but keep only the first two
terms, which is actually a valid expansion for all index $n$ polytropes [34].
In addition to being general (i.e. $n$ independent), the reason for truncating
at the second order term is that the DM distribution profile $n_{X}$ becomes
highly cored, whenever $m_{X}\gg 1$ GeV. Later in this section, we show how
the Lane-Emden function appears in integrals over the DM distribution (See
Eqs. 3.11, 3.14). Equation 3.7 is a sufficient, in most cases. In fact,
without this approximation, computing such integrals for high DM mass
numerically is challenging because of how cored the profiles are. For sub-GeV
DM, however, we do not make any approximations on the Lane-Emden function
since these integrals are straightforward to compute numerically.

With the hydrogen distribution in hand, we turn to approximating the
distribution of captured DM in these objects. We must first determine how
frequently the DM particles interact with the hydrogen distribution. To do so,
we use the Knudsen number, as defined in Eq. 2.19. In Fig. 3, we scan over a
part of the $\sigma-m_{X}$ parameter space for $M=1M_{J}$ and $M=14M_{J}$
Jupiters and compute the Knudsen number. In this figure, it is clear that,
depending on the object, there is a region where we are not entirely in the
$\text{Kn}\gg 1$ or $\text{Kn}\ll 1$ limit and so we must interpolate between
the two regimes as we traverse the $\text{Kn}\sim 1$ boundary using the
interpolation function defined in Eq. 2.22. Moreover, we note that this
transition happens within the region of the $\sigma$ vs. $m_{X}$ parameter
space probed by Jupiters (see Fig. 3 of [30] and Fig. 5 in this paper).

Figure 3: Knudsen number for a $M=1M_{J}$ (left panel) and $M=14M_{J}$ (right
panel) jupiter-like exoplanet as a function of the DM-hydrogen scattering
cross section ($\sigma$) and DM mass ($m_{X}$). The solid black line
represents the $\text{Kn}=1$ contour, and so we see that the region of
$\sigma-m_{X}$ parameter space in which we are interested straddles the
$\text{Kn}\sim 1$ boundary.

Let us first suppose that we are in the Knudsen limit ($\text{Kn}\gg 1$).
Then, as shown in [31], the DM distributes itself isothermally at a
temperature $T_{X}$. The resulting profile is given by:

$n_{X}(\xi)=n_{X}(0)e^{-\frac{m_{X}\Phi(\xi)}{T_{X}}}.$ (3.8)

The main question to address in an isothermal profile is the temperature of
the captured DM distribution. Following [31], we find the distribution’s
temperature by solving the first energy moment of the collisional Boltzmann
equation. In doing so, we are designating by $T_{X}$ the average DM
temperature over a single orbit since, in reality, as DM particles orbit
within the object, their temperature will undergo slight changes as they
interact with the local medium. For time-independent distributions, such as
the one we are considering, requiring that the first moment is satisfied is
equivalent to there being no net flow of energy into the dark matter
distribution from the surrounding hydrogen in the planet. The equation
necessary to solve the DM temperature is [31]:

$\int_{0}^{\xi_{1}}\theta(\xi)\exp\left(-\frac{\mu}{\Theta}\Phi(\xi)\right)\left(\frac{\Theta+\mu\theta(\xi)}{\mu}\right)^{1/2}[\Theta-\theta(\xi)]\xi^{2}d\xi=0,$
(3.9)

where we define the following dimensionless variables for computational
convenience:

$\mu\equiv\frac{m_{X}}{m}\text{ ,
}\tilde{\Phi}(r)\equiv\frac{m\Phi(r)}{T_{c}}\text{ and
}\Theta\equiv\frac{T_{X}}{T_{c}}.$ (3.10)

Combining Eq. 3.8 and Eq. 3.5, we see that the dependence of $n_{X}(\xi)$ in
the isothermal limit is given by:

$n_{X}\sim\exp\left\\{-\gamma\left(1-\theta(\xi)\right)\right\\},$ (3.11)

where $\gamma\sim m_{X}/T_{X}$. Additionally, one can show that if $m_{X}\gg
m$, $T_{X}\approx T_{c}$ since in Eq. 3.9, if $\mu\to\infty$, the exponential
term cuts off the integral to the area around $\xi\to 0$. Thus, the only way
the LHS of Eq. 3.9 is zero is if $\Theta\to\theta(1)=1$, i.e. $T_{X}\to
T_{c}$. In Appendix A, we will explore the $m_{X}\to\infty$ limit to find
approximations for various integrals involving $n_{X}$. This was done
primarily by the recognition that numerical integration on exponential
functions with large decay constants is challenging. In the opposite limit of
$m_{X}\to 0$, one can show that $T_{X}\approx\text{const.}$, and so $\gamma\to
0$. In this case, the isothermal profile simply becomes constant, $n_{X}\simeq
n_{X}(0)$. This is a fact that is unique to the isothermal distribution. As we
will show below, the DM profile in the opposite ($\text{Kn}\ll 1$) limit is
not flat if $m_{X}\to 0$, but rather obtains a “floating” profile with higher
densities toward the edge of the object.

Now, we shall explore the limit where $\text{Kn}\ll 1$, also known as the
“Local Thermal Equilibrium” or LTE limit. In [57], it was shown that the DM
profile in this limit is given by:

$n_{X}(r)=n_{X}(0)\left(\frac{T(r)}{T_{c}}\right)^{3/2}\exp\left(-\int_{0}^{r}\frac{\alpha(r^{\prime})\frac{dT}{dr^{\prime}}+m_{\chi}\frac{d\Phi}{dr^{\prime}}}{T(r^{\prime})}dr^{\prime}\right),$
(3.12)

where $\alpha(r)$ is a factor related to the thermal diffusivity of DM. As
explained in [57], $\alpha(r)$ is actually a constant, $\alpha_{0}$, if we
consider one nuclear species and can be approximated analytically as shown in
Appendix A of [57]. With this in mind, and defining
$E_{e}\equiv\frac{1}{2}m_{X}v_{\text{esc}}^{2}$, we start by making progress
on the integral in the exponent:

$\begin{split}\int_{0}^{r}\frac{\alpha(r^{\prime})\frac{dT}{dr^{\prime}}+m_{X}\frac{d\Phi}{dr^{\prime}}}{T(r^{\prime})}dr^{\prime}&=\int_{0}^{\xi}\frac{\alpha_{0}T_{c}\frac{d\theta}{d\xi^{\prime}}-\frac{\xi_{1}}{2\alpha_{n}}m_{X}v_{\text{esc}}^{2}\frac{d\theta}{d\xi^{\prime}}}{T_{c}\theta(\xi^{\prime})}d\xi^{\prime},\\\
&=\left(\alpha_{0}-\frac{\xi_{1}}{\alpha_{n}}\frac{E_{e}}{T_{c}}\right)\int_{0}^{\xi}\frac{1}{\theta(\xi^{\prime})}\frac{d\theta}{d\xi^{\prime}}d\xi^{\prime},\\\
&=\left(\alpha_{0}-\frac{\xi_{1}}{\alpha_{n}}\frac{E_{e}}{T_{c}}\right)\int_{0}^{\xi}\frac{d}{d\xi^{\prime}}\left[\log\left(\theta(\xi^{\prime})\right)\right]d\xi^{\prime},\\\
&=\left(\alpha_{0}-\frac{\xi_{1}}{\alpha_{n}}\frac{E_{e}}{T_{c}}\right)\log\left(\theta(\xi)\right),\end{split}$
(3.13)

We have used $\Phi(\xi)$ for polytropes, as given by Eq. 3.5 and the fact that
$T\sim\theta$, as appropriate for gas pressure dominated objects, such as
exoplanets.888As it turns out, $T\sim\theta$ even for Brown Dwarfs [56], the
other class of objects we consider in this paper.Plugging in to Eq. 3.12, we
get:

$\begin{split}n_{X}(\xi)&=n_{X}(0)\theta^{3/2}(\xi)\exp\left(-\left(\alpha_{0}-\frac{\xi_{1}}{\alpha_{n}}\frac{E_{e}}{T_{c}}\right)\log\left(\theta(\xi)\right)\right)\\\
&=n_{X}(0)\theta^{3/2}(\xi)\theta^{\frac{\xi_{1}}{\alpha_{n}}\frac{E_{e}}{T_{c}}-\alpha_{0}}(\xi)\\\
&=n_{X}(0)\theta^{p}(\xi),\end{split}$ (3.14)

where
$p\equiv\frac{\xi_{1}}{\alpha_{n}}\frac{E_{e}}{T_{c}}+\frac{3}{2}-\alpha_{0}$.
Thus, we see that, in the LTE limit, the DM profile behaves as $n_{X}\sim
n^{p}$ where $n(\xi)=n_{c}\theta(\xi)$ is the hydrogen number density profile
and $p(m_{X})$ is a mass-dependent quantity. In fact, one can show that as
$m_{X}\to 0$, $\alpha_{0}\approx 2.1$. Since $E_{e}\sim m_{X}$ and
$\alpha_{0}>3/2$, it follows that for sufficiently small $m_{X}$ the exponent
$p$ becomes negative. For example, assuming $n=1$, as appropriate for
Jupiters, get $p\approx-0.6$, in the limit of $m_{X}\to 0$. This result is
intriguing and somewhat counter-intuitive since as $\xi\to\xi_{1}$,
$\theta(\xi)\to 0$ and so $n_{X}\to\infty$. In other words, for small enough
DM mass, i.e. where $p(m_{X})<0$, the profile forms a ring-like structure with
higher densities toward the edge of the object. In Fig. 4 (bottom row) this
phenomenon becomes clear in the plots for the density profiles for DM inside a
$1M_{J}$ Jupiter. A similar result was found recently for the Earth [39] and
general celestial bodies [40].

Thus far, we have examined the density profile of DM particles in exoplanets
in the extremes of high and low interaction. However, as is evidenced by Fig.
3, much of the $\sigma-m_{X}$ parameter space in which we are interested
exists in areas where $\text{Kn}\sim 1$. Thus, following [46], we approximate
the true DM profile at a point in $\sigma-m_{X}$ parameter space via Eq. 2.24.
Note that, as expected:

$\lim_{\text{Kn}\to\infty}f(\text{Kn})=0\implies\lim_{\text{Kn}\to\infty}n_{X}(r)=n_{X}^{\text{ISO}}(r),$

and so the Isothermal profile is recovered in the Knudsen limit. Similarly:

$\lim_{\text{Kn}\to 0}f(\text{Kn})=1\implies\lim_{\text{Kn}\to
0}n_{X}(r)=n_{X}^{\text{LTE}}(r),$

and the local thermal equilibrium profile is recovered in the appropriate
limit.

Figure 4: Dark Matter density profiles in a $1M_{J}$ Jupiter for
$m_{X}=0.1\,\mathrm{GeV}$ (top panels) and $m_{X}=3.4\times
10^{-2}\,\mathrm{GeV}$ (bottom panels) particles in three distinct regimes:
local thermal equilibrium (left panels), Knudsen transition (central panels),
and isothermal (right panels). The $m_{X}=3.4\times 10^{-2}\,\mathrm{GeV}$ DM
particle was chosen since at this mass, $p(m_{X})\sim-0.6$, making it nice to
visualize the floating profile, while for $m_{X}=0.1\,\mathrm{GeV}$, the
profile is closer to the hydrogen profile. In each panel, the red line
indicates the true DM profile given by Eq. 2.24, while the dashed green line
is the LTE profile $n_{X}^{\text{LTE}}$ as per Eq. 3.14 and the dashed blue
line is the isothermal profile $n_{X}^{\text{ISO}}$ as per Eq. 3.8.

In Fig. 4, we plot DM density profiles for a DM mass of
$m_{X}=0.034\,\mathrm{GeV}$ and $m_{X}=0.1\,\mathrm{GeV}$ with varying values
of $\sigma$ to demonstrate the three possible regimes: local thermal
equilibrium, Knudsen transition, and isothermal. The range of masses were
chosen to highlight the transition between a cored and floating profile in the
LTE limit. In the left and right panels of this plot, we see the red line,
which is the true DM profile, converge on the green and blue lines
respectively, demonstrating the reduction of the full profile to the two
limiting regimes of LTE and isothermal distributions. The central panels
represents a scenario when $\text{Kn}\sim 1$, and we clearly see that the true
profile is a combination of the LTE and isothermal profiles. Another
intriguing feature in Fig. 4 is the “floating” nature of the LTE profile in
the bottom panel. This occurs when the exponent $p$ in Eq. 3.14 becomes
negative which, for a $1M_{J}$ Jupiter, occurs for DM masses below $m_{X}\sim
0.04\,\mathrm{GeV}$. As mentioned in the derivation of Eq. 3.14, a similar
result of higher DM densities towards the edge of the planet was recently
shown to apply to the Earth [39], though not through the use of polytropes. In
[39], the authors apply the same formalism for a distribution of DM that
interacts frequently with the object’s distribution, established in [57], and
also find a floating DM profile in the Earth.

###  3.3 Dark Matter sensitivity estimates with exoplanets

For self-annihilating DM models, the annihilation of captured DM particles in
exoplanets could heat the object and through observations of the object’s
temperature, it is possible to constrain this DM model in $\sigma-m_{X}$
parameter space. As suggested by [30], exoplanets are excellent sub-GeV DM
probes, and as such, we focus on such DM models here. In order to bypass the
Lee-Weinberg limit [42] for WIMPs ($m_{X}\gtrsim 1$ GeV), and still have a
thermal relic for which the annihilation cross section is fixed by the relic
abundance, one needs to allow for DM interacting stronger that the weak scale
[e.g. see 58]. Here we follow [30] and restrict our attention to the Co-SIMP
DM model of [41].

An exoplanet in thermal equilibrium will emit a luminosity $L$ at a
temperature $T$ according to the Stefan-Boltzmann law:

$L_{\text{obs}}=4\pi\epsilon R^{2}\sigma_{\text{SB}}T_{\text{obs}}^{4},$
(3.15)

where $\sigma_{\text{SB}}$ is the Stefan-Boltzmann constant, and $\epsilon$ is
the emisivity of the atmosphere. For simplicity, we shall assume the
emissivity of the object is $\epsilon\simeq 1$. If $\epsilon<1$, the object
will become hotter by trapping outgoing heat-flow. The total flux from the
exoplanet will stay constant on the basis of energy conservation, but the
spectral flux density curve becomes shifted to shorter wavelengths, in view of
the Wien’s displacement law, and the constancy between the $\epsilon
T_{obs}^{4}$ product. In the context of JWST searches, this is beneficial
since short-wavelength bands generally have greater flux sensitivity [30].

In principle, the total luminosity comes from a combination of sources
including internal heating, DM heating, and any other external heat sources.
Assuming thermal equilibrium, we write the total luminosity as a sum of three
components:

$L_{\text{obs}}=L_{\text{int}}+L_{\text{DM}}+L_{\text{ext}},$

where $L_{\text{int}}\equiv 4\pi R^{2}\sigma_{\text{SB}}T_{\text{int}}^{4}$
represents the internal heating rate (from all sources except DM),
$L_{\text{DM}}$ stands for the heating rate due to captured DM annihilations,
and $L_{\text{ext}}$ is the rate of heating from external sources, such as the
parent star. If, instead
$L_{\text{obs}}>L_{\text{int}}+L_{\text{DM}}+L_{\text{ext}}$, then the object
would be cooling and our sensitivity estimates become bounds. Following [30],
we set $L_{\text{ext}}=0$ since we are mostly interested in old, isolated
exoplanets for cleaner measurements of their temperature. We note here that
determining if a given exoplanet is in thermal equilibrium or cooling is
difficult, particularly for isolated objects. Provided one can find a way to
estimate their age, we can determine whether or not the object’s temperature
is above the expected value which would indicate a DM signal. What follows is
valid regardless of $L_{\text{ext}}$, except that the sensitivity limits are
instead transformed to bounds. For the internal heating, we follow [59] where
evolutionary models are used to obtain the temperature of old exoplanets of
various masses and we interpolate to find $T_{\text{int}}$. This quantity
($T_{int}$) represents the exoplanet (mass dependent) effective temperature
contribution from internal heat sources (excluding Dark Matter annihilations)
to the object’s surface temperature. One can also view $T_{\text{int}}$ as the
minimum surface temperature the object will have, i.e. its temperature in the
absence of DM heating and external heating, and we emphasize that this is not
the internal temperature of the object. In appendix C, we outline how we
estimate $T_{\text{int}}$ for Jupiters and Brown Dwarfs.

The final component of the exoplanet’s heating comes from DM annihilation.
Then, given that equilibrium between capture, annihilation, end evaporation is
achieved in relatively short time-scales for the Co-SIMP models considered
here [30] we can express the total, constant, number of DM particles in the
object as:

$N_{X}=\sqrt{\frac{C}{C_{A}}}\frac{1}{\kappa+\frac{1}{2}E\tau_{\text{eq}}}.$

To arrive at this, consider the limit of $t\gg\tau_{\text{eq}}/\kappa$ in Eq.
2.3, which gives the total number of DM particles in the object at
equilibrium. Then, since $L_{\text{DM}}=m_{X}\Gamma_{A}=m_{X}C_{A}N_{X}^{2}$,
we arrive at the following equation for the DM luminosity:

$L_{\text{DM}}=\frac{Cm_{X}}{\left(\kappa+\frac{1}{2}E\tau_{\text{eq}}\right)^{2}},$
(3.16)

where we remind the reader that
$\kappa\equiv\sqrt{1+\left(\frac{1}{2}E\tau_{eq}\right)^{2}}$. Note that, if
$E\approx 0$, then $L_{\text{DM}}=m_{X}C$, so the DM luminosity with
evaporation is the DM luminosity without evaporation but suppressed by the
factor $\left(\kappa+1/2E\tau_{\text{eq}}\right)^{-2}$.

Ultimately, we are interested in solving
$L_{\text{obs}}=L_{\text{int}}+L_{\text{DM}}$ in $\sigma-m_{X}$ parameter
space. This equation can be recast in the following way:

$\frac{Cm_{X}}{\left(\kappa+\frac{1}{2}E\tau_{\text{eq}}\right)^{2}}=4\pi
R^{2}\sigma_{\text{SB}}\left(T_{\text{obs}}^{4}-T_{\text{int}}^{4}\right).$
(3.17)

We again stress that $T_{\text{int}}$ is not the internal temperature of the
object, but rather the contribution to the object’s surface temperature from
any sources excluding: Dark Matter or external heating.

At this point, we draw attention to an important and often overlooked
consequence of the existence of Region II in the $\sigma-m_{X}$ parameter
space (see Fig. 1), where the total capture rate saturates. Namely, for all
$m_{X}$ for which DM evaporation is insignificant (compared to capture and
annihilation) the constraining power of exoplanets as DM probes is completely
erased if the object is hotter than a critical temperature, as we will explain
below. Simply put, the lhs of Eq. 3.17 becomes $\sigma$ independent, and, as
such, this becomes a completely degenerate parameter. The statement above can
be recast in terms of the existence of a DM and astrophysical parameters
dependent critical temperature, $T_{\text{crit}}$, above which any
astrophysical object cannot serve as a DM probe. It is straightforward to
verify that the product $Cm_{X}$ is maximized in Region II of $\sigma-m_{X}$
parameter space, which corresponds to a cross section so high that all DM
particles crossing an object are captured, i.e. the geometric limit. Thus, let
us consider Eq. 3.17 in the limit that $E\to 0$ and assume we are in Region
II, where the capture rate is maximized:

$C_{\text{max}}m_{X}=4\pi
R^{2}\sigma_{\text{SB}}\left(T_{\text{obs;max}}^{4}-T_{\text{int}}^{4}\right).$
(3.18)

Now, since we are assuming the only heat sources are from DM annihilation and
internal heat, this equation actually implies an upper bound on
$T_{\text{obs}}$ which we label $T_{\text{crit}}$. For any
$T_{\text{obs}}>T_{\text{crit}}$, it is impossible for the planet to be heated
to the observed temperature purely by DM and internal heat sources because we
have found maximal heating at $T_{\text{crit}}$ (by definition of $T_{crit}$).
Solving for $T_{\text{crit}}$ by substituting Eq. 2.13 for $C_{\text{max}}$ we
find:

$T_{\text{crit}}=\left(\frac{1}{\sqrt{24\pi}}\frac{\rho_{X}}{\sigma_{\text{SB}}\bar{v}}\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)+T_{\text{int}}^{4}\right)^{\frac{1}{4}}.$
(3.19)

Note that the critical temperature implicitly depends on the location of the
exoplanet since the DM density and velocity dispersion play a key role in
maximizing annihilation and are themselves dependent on location. To sum up,
for every object for which we know $T_{\text{int}}$, its escape velocity
($v_{\text{esc}}$), the DM density at its location ($\rho_{X}$), and the DM
dispersion velocity ($\bar{v}$) we have identified a critical temperature
($T_{\text{crit}}$) according to Eq. 3.19. If the object is observed with an
inferred $T_{\text{obs}}>T_{\text{crit}}$, then we necessarily know that DM
heating cannot be, alone, responsible for this overheating when compared to
the expected $T_{\text{int}}$. This, in turn, implies the existence of
external heat sources. Unless measurements are accurate enough to determine
$L_{\text{ext}}$, one cannot, under those circumstances, make any predictions
about $L_{\text{DM}}$, as those two are now degenerate. In turn, this implies
that the object lost all of its constraining power, whenever
$T_{\text{obs}}>T_{\text{crit}}$ and there is no accurate determination of
$L_{\text{ext}}$ possible. If instead we are dealing with an isolated (rogue)
exoplanet, for which $L_{\text{ext}}\equiv 0$ by definition, and for which one
measures $T_{\text{obs}}>T_{\text{crit}}$, then this implies the existence of
other unaccounted for internal heating mechanisms, excluding DM. Regardless,
this would imply a complete loss of DM constraining power for this object.

With the recent launch of the James Webb Space Telescope (JWST), there has
been a lot of excitement about the possibility for exoplanet detection. An
extensive discussion on exoplanets and their usefulness as DM probes is very
nicely summarized in [30], along with useful discussions on the possibilities
for detection with JWST. Included in this analysis is a determination of
JWST’s lowest possible temperature reach for objects located near the center
of the galaxy. This limit is nicely summarized in Figure S2 of [30].
Specifically, they find that a likely JWST detection near the galactic center
implies a minimum observable temperature of $T_{\text{min}}=650\,\mathrm{K}$
given an exposure time of $10^{5}\,\mathrm{s}$. In other words, it is not
possible for JWST to detect an exoplanet with an observed temperature below
$T_{\text{min}}=650\,\mathrm{K}$ given the telescopes sensitivity, and
assuming a reasonable exposure time. We stress that this is a mere limitation
of the instrument and not representative of any fundamental limit. In fact,
one can show that, under the commonly made assumption of Poisson “noise,” we
have $T_{\text{min}}\sim t_{\text{exp}}^{-1/8}$. This leads to the following
exposure time dependence of $T_{\text{min}}$ at for Jupiters and BDs near the
galactic center:

$T_{\text{{min}}}(t_{\text{{exp}}})=650\times\left(\frac{t_{\text{{exp}}}}{10^{5}\mathrm{s}}\right)^{-1/8}\mathrm{K}$
(3.20)

Also note that Jupiters and Brown Dwarfs have similar radii at $\sim R_{J}$,
so $T_{\text{min}}$ is object-independent. As we are focusing on planets found
toward the galactic center, we thus adopt this value (i.e. 650 K) as the
smallest observed temperature considered, thus assuming a hypothetical
$10^{5}$ s exposure. For an actual observation, one can simply find the
rescaled $T_{\text{min}}$ using Eq. 3.20. Thus, when computing sensitivity
limits to DM through exoplanet detection, we do no consider objects with
observed temperature below $T_{\text{min}}$, since these are out of reach of
JWST. Because of this, it does not make sense to consider objects whose
critical temperature is below $T_{\text{min}}$ (i.e.
$T_{\text{crit}}<T_{\text{min}}$) because if they are observed, their
temperature $T_{\text{obs}}$ is necessarily above $T_{\text{min}}$, which
means $T_{\text{obs}}>T_{\text{crit}}$, and, as discussed above, no
constraints can be placed on the DM parameter space if this is the case. Thus,
we can only consider those objects whose critical temperature surpasses the
minimum observational temperature, i.e. $T_{\text{crit}}>650\,\mathrm{K}$. For
Jupiters around $1M_{J}$, the internal heat gives rise to an expected
temperature (in absence of any other heating sources) $T_{\text{int}}\sim
80\,\mathrm{K}$. Thus, one can verify by varying Jupiter mass in Eq. 3.19 and
taking galactocentric DM parameters (see paragraph below), that the constraint
of $T_{\text{crit}}>T_{\text{min}}$ restricts us to exoplanets with mass
$M\gtrsim 3M_{J}$, as shown in Appendix C.

After this important discussion regarding the possible limitations of the
usefulness of exoplanets (and astrophysical probes in general) whenever their
surface temperature is larger than the critical temperature
($T_{\text{crit}}$) we identified, we go back to the scenarios where
$T_{\text{obs}}<T_{crit}$. To solve Eq. 3.17, we need to have a grasp on
further parameters like the DM dispersion velocity and DM density around the
object. These parameters depend on the location of the exoplanet and will
affect the sensitivity limits. For the purpose of being concrete, we shall
consider Jupiter-like exoplanets of masses $3M_{J}$ and $14M_{J}$ near the
galactic center and a Brown Dwarf of mass $70M_{J}$ near the galactic center.
The choice of the location (GC) is motivated by the high DM density, and
therefore it offers the deepest reach (sensitivity) to DM any probe can have.
The lower mass limit is chosen in consideration of the requirement that
$T_{\text{crit}}>650\,\mathrm{K}$. For each object we obtain the central
temperature ($T_{c}$) using the numerical simulations of [38].999The values of
$T_{c}$ estimated this way are consistent, to within a few percent, with those
presented in Eqns. 3.13 and 3.16 of [33]. Focusing on exoplanets found near
the galactic center means that for the DM dispersion velocity we have
$\bar{v}\approx 50\,\mathrm{km/s}$ and for the DM density, $\rho_{X}\approx
10^{3}\,\mathrm{GeV/cm^{3}}$ [30]. The density does depend on the profile
chosen, though they are all of order $10^{3}\,\mathrm{GeV/cm^{3}}$ [60].

Figure 5: Sensitivity limits in $\sigma-m_{X}$ space for the detection of
$M=3M_{J}$ (left panels) and $M=14M_{J}$ (right panels) Jupiters at
temperatures $650\,\mathrm{K}$ (top panels) and $T_{\text{crit}}$ (bottom
panels, values of $T_{\text{crit}}$ labeled in the title of each panel). We
use a GC DM velocity and density. The dashed lines indicate the true
sensitivity limits which account for evaporation, including the ping pong
suppression effect of Eq. 2.17, while the solid lines are sensitivity limits
if evaporation is neglected. For consitency with previous literature, we also
include sensitivity with evaporation where the suppression factor is set to
$1$, as seen by the Dash-dot lines. The green region is excluded by spin-
dependent (XENON1T and CDEX-1B) direct detection experiments [14, 61, 62].
Spin independent experiments are even more sensitive, and would have already
ruled out the entire region of parameter space probed by Jupiters. and This
plot demonstrates that DM evaporation severely limits the usefulness of
Jupiters as sub-GeV DM probes. Figure 6: Sensitivity limits in $\sigma-m_{X}$
space for the detection of $M=70M_{J}$ Brown Dwarf at temperatures
$T_{obs}=1100\,\mathrm{K}$ (left panel) and $T_{obs}\simeq T_{crit}\simeq
1500\,\mathrm{K}$ (right panel). We use a GC DM velocity and density. The
dashed lines indicate the true sensitivity limits which account for
evaporation, including the ping-pong effect, the dash-dot lines account for
evaporation without the ping-pong effect, and the solid colored lines are
sensitivity limits if evaporation is neglected. The green region is excluded
by spin dependent direct detection experiments. This plot demonstrates that DM
evaporation severely limits the usefulness of Brown Dwarfs as sub-GeV DM
probes.

In Fig. 5, we plot the sensitivity limits of two Jupiters near the galactic
center, $M=3M_{J}$ (left panels) and $M=14M_{J}$ (right panels), at both the
minimum temperature for JWST detection, $T_{\text{min}}=650\,\mathrm{K}$ (top
row), and a temperature very close to the critical temperature,
$T_{\text{crit}}$ (bottom row). In Fig. 6, we plot sensitivity limits for a
$M=70M_{J}$ Brown Dwarf at temperatures $1100K$ and $1500K$, where the first
temperature is chosen to be above the temperature due the internal heat source
(i.e. $T_{int}$) while the second is just below $T_{\text{crit}}$. Near the
galactic center, a $3M_{J}$ Jupiter has a critical temperature
$T_{\text{crit}}\simeq 700\,\mathrm{K}$, a $14M_{J}$ Jupiter has a critical
temperature $T_{\text{crit}}\simeq 1000\,\mathrm{K}$, and a $70M_{J}$ Brown
Dwarf has a critical temperature $T_{\text{crit}}\simeq 1550\,\mathrm{K}$. In
both Figs. 5 and 6 we consider three cases for the sensitivity limits: i.
including evaporation suppressed by the $s(r)$ factor given by Eq. 2.17
(dashed lines); ii. maximum evaporation, i.e. $s(r)$ is set to one (dash-
dotted lines); and iii. no evaporation considered at all (solid lines). Since
evaporation decreases the number of DM particles available for annihilation in
an object, we expect that the inclusion of evaporation lowers the effect of DM
heating and thus leads to weaker constraints in $\sigma-m_{X}$ space. This is
confirmed by Figs. 5 and 6, where we see the dashed and dash-dot lines
(Evaporation) diverge upwards from the solid line (No Evaporation), as
evaporation begins to kick in. The fact that the constraints coincide for
masses above some object-dependent cutoff is expected since evaporation is
exponentially suppressed for DM masses above the “evaporation mass.” This
quantity (the evaporation mass) is estimated for arbitrary objects in [33], by
carefully including the effects of the tail of the Maxwell-Boltzmann
distribution. This important effect was previously neglected when naive
estimates for the evaporation mass were made in the literature. For $m_{X}$
above this mass, $E\approx 0$, and so the factor suppressing the sensitivity
limits, $\left(\kappa+\frac{1}{2}E\tau_{\text{eq}}\right)^{2}$, is
approximately unity, and the sensitivity with or without evaporation is
identical. One important takeaway message is that DM evaporation (even when
one includes its suppression) severely restricts the usefulness of Jupiters
and Brown Dwarfs as sub-GeV DM probes, as can be explicitly seen from Figs. 5
and 6.

In what follows we discuss the expected shape of the sensitivity limits,
whenever evaporation is not included (i.e. the solid lines). First, remember
that those are generated as solutions to Eq. 3.17, which, in turn, depend on
the region of $\sigma-m_{X}$ parameter space. For example, if we consider the
$M=14M_{J}$ sensitivity limits for $T_{\text{obs}}=650\,\mathrm{K}$ as seen in
the upper right panel of Fig. 5, we see that the solution below $m_{X}\sim
0.1\,\mathrm{GeV}$ sits at $\sigma\sim 10^{-36}\,\mathrm{cm^{2}}$, which one
can confirm sits in Region III by computing $k$ and $\tau$ as defined in
Section 2.1. Also, since we restrict our attention on the shape of the solid
lines (no evaporation), the problem of finding the sensitivity limits is
reduced to:

$\left(2A\tau v_{\text{esc}}^{2}\right)m_{X}=4\pi
R^{2}\left(T_{\text{obs}}^{4}-T_{\text{int}}^{4}\right),$ (3.21)

where we have substituted $C$ with the analytic version in Region III given in
Eq. 2.14. Then, recalling that $A\sim m_{X}^{-1}$ and $\tau\sim\sigma$, we see
that the sensitivity limits obey $\sigma\sim m_{X}^{0}$. This expectation is
confirmed in all panels of Figs. 5 and 6, where in each case one can find a
flat region of the solid lines (sensitivity limits without evaporation), which
corresponds to the $m_{X}$ range of Region III of the $\sigma$ vs $m_{X}$
parameter space. We note here that the general shape of those limits is a
“bowl shaped” curve, symmetric with respect to $m_{X}\sim 1$ GeV, which
corresponds to the mass of the target particle inside the object, i.e. in this
case protons. This symmetric shape can be traced back to the symmetry in the
regions identified by us in Fig. 1 with respect to the same vertical axis
corresponding to $m_{X}\sim 1$ GeV.

From Figs. 5 and 6 we also see that at either very large or very small DM mass
(i.e. outside of the shallow part of the “bowl" described above) there is some
sensitivity loss, even if DM evaporation is completely negligible (see solid
lines). This effect can be traced back to the difference in the capture rates
in Regions III and IV, i.e. to the appearance of a factor of $k$ in Region IV
when compared to Region III. As $k$ is either linear with $m_{X}$ (whenever
$m_{X}<m$) or inversely proportional with $m_{X}$ (whenever $m_{X}>m$). This
leads to sensitivity limits that are slanted as $\sigma\sim m_{X}$ or
$\sigma\sim m_{X}^{-1}$ once we cross into Region IV or parameter space, at
either end of Region III.

We point out that, everything else being the same, cooler objects probe deeper
in to the $\sigma$ vs $m_{X}$ parameter space. This can be explicitly seen
from contrasting for instance the top and bottom panels of Fig. 5 or the left
and right panel of Fig. 6 and is a direct consequence of Eq. 3.17. As the
observed temperature increases, the sensitivity limits shift upwards, until
$T_{obs}\simeq T_{crit}$, when the limits depicted by the solid lines (no
evaporation) in both Fig. 5 and 6 perfectly trace the boundary of Region II of
parameter space, as expected. As previously explained, once $T_{obs}>T_{crit}$
the DM constraining of those objects power is lost.

##  4 Summary and discussion

In this work, we (re)examined the effectiveness of exoplanets and Brown Dwarfs
as probes of Dark Matter, by including the effects of DM evaporation. Our key
result is that DM evaporation, even when suppressed, almost completely erases
all sensitivity at $m_{X}\lesssim 1\,\mathrm{GeV}$ (see Figs. 5 and 6). While
in contrast with the findings of [30] (who neglect the role of DM
evaporation), our results are in agreement with the findings of [33], who
estimate the evaporation mass for arbitrary celestial bodies by including the
often neglected role of the tail of the DM distribution.

Furthermore, we stress the importance of a new parameter, the critical
temperature $T_{\text{crit}}$, that is essential to the viability a given
object has a probe of annihilating dark matter via potential detection of an
extra heating signal. Note that the critical temperature, as defined in Eq.
3.19, is not only a function of the exoplanet mass, but implicitly its
location, since the DM heating signal depends on the local DM density and
velocity dispersion. This upper limit on observed exoplanet temperatures
arises from considering the maximum heating signal that DM can produce,
corresponding to the geometric capture limit. Simply put, exoplanets hotter
than their corresponding critical temperature are not sensitive to DM since no
regions in $\sigma-m_{X}$ space could produce a signal where
$T_{\text{obs}}>T_{\text{crit}}$.101010This assumes there is no external
heating source, though introducing one would simply mean adding a term in the
parentheses of Eq. 3.19 and thus increasing the critical temperature.

The existence of the critical temperature ($T_{\text{crit}}$) also implies a
minimum mass for a Jupiter such that it can be considered a useful DM probe.
Namely, in order for an object to be detected in a reasonable exposure time,
it has to have a minimum surface temperature ($T_{\text{min}}$). However,
objects for which $T_{\text{min}}>T_{\text{cirt}}$ would not be sensitive to
DM, as explained above. Conversely, one should only consider objects for which
$T_{\text{crit}}>T_{\text{min}}$. As explained in Appendix C, this inequality
can be converted into a lower bound on the mass of the DM probe. For Jupiters
at the galactic center this becomes $M\gtrsim 3M_{J}$.

Another striking result worthy of discussion is the floating profile for light
DM particles in the Local Thermal Equilibrium limit. As discussed in Section
3.2, this feature arises when DM interacts frequently with the Exoplanet’s
hydrogen distribution and is a highly counter-intuitive result. Typically, we
expect DM to develop a cored profile, with higher densities toward the center
of the object. However, when the DM mass falls below a point where
$p\equiv\frac{E_{e}}{T_{c}}+\frac{3}{2}-\alpha_{0}$ becomes negative, the DM
will tend to accumulate towards the surface of the object. Similar results in
recent work [39, 40] are exciting confirmations of this effect.

We end by summarizing our main finding. In principle, exoplanets and Brown
Dwarfs could provide rich laboratories for the exploration of Dark Matter, as
pointed out by [30]. However, DM evaporation severely hinders their usefulness
as DM probes in the poorly constrained sub-GeV regime, leaving only small
swaths of the $\sigma-m_{X}$ parameter space that are not already ruled out by
direct detection experiments.

##  Appendix A Analytic approximations of annihilation coefficients

In this section, we use the results of Section 3.2 to derive analytic
approximations for the annihilation and evaporation coefficients of Eq. 2.1.
For the high-mass limit, we will compute the annihilation coefficient since
the DM profiles become extremely cored in this limit, making numerical
integration very difficult. We will also compute an analytic evaporation rate,
though it is worth pointing out that in the high-mass limit this is not very
useful since evaporation is exponentially suppressed with mass making it
irrelevant above some cutoff.

We will start by looking at integrals of the form

$\int dVn_{X},~{}\text{and}~{}\int dVn_{X}^{2},$

since these are the primary integrals we wish to approximate. We are
interested in these integrals because they appear in the computation of
annihilation coefficients as per Eqs. 2.29 and 2.30, and they allow one to
compute the total number of DM particles in the object. First, we remind the
reader of the Taylor expansion of the Lane-Emden function $\theta(\xi)$ per
Eq. 3.7, though we only keep up to 2nd order:

$\theta(\xi)\simeq 1-\frac{1}{6}\xi^{2}.$

Then, we can analytically compute the following DM distribution integral in
the LTE limit:

$\begin{split}\int
dVn_{X}^{\text{LTE}}(\xi)&=4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}(0)\int_{0}^{\xi_{1}}d\xi~{}\xi^{2}\theta^{p}(\xi),\\\
&\simeq
4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}(0)\int_{0}^{\sqrt{6}}d\xi~{}\xi^{2}\left(1-\frac{1}{6}\xi^{2}\right)^{p},\\\
&=4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}(0)\sqrt{\frac{27\pi}{2}}p^{-3/2}.\end{split}$
(A.1)

We remind the reader that $p\equiv\frac{E_{e}}{T_{c}}+\frac{3}{2}-\alpha_{0}$.
Note how in the second line the bounds of the integral change from
$0\to\xi_{1}$ to $0\to\sqrt{6}$, since $\sqrt{6}$ is the first zero of the
approximated function and if this is not changed then, depending on $p$, we
may get non-physical values for $n_{X}^{\text{LTE}}(\xi)$. In the high DM mass
limit, this choice is inconsequential because the profiles become increasingly
cored and so the integral converges long before $\xi=\sqrt{6}$. This will not
be the case for the purely isothermal result. We can then proceed to do the
same integral for the isothermal distribution:

$\begin{split}\int
dVn_{X}^{\text{ISO}}(\xi)&=4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}(0)\int_{0}^{\xi_{1}}d\xi~{}\xi^{2}\exp\left\\{-\frac{m_{X}v_{\text{esc}}^{2}}{2T_{X}}\left(1-\theta(\xi)\right)\right\\},\\\
&\simeq
4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}(0)\int_{0}^{\xi_{1}}d\xi~{}\xi^{2}\exp\left\\{-\frac{m_{X}v_{\text{esc}}^{2}}{12T_{X}}\xi^{2}\right\\},\\\
&=4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}(0)\left[\frac{\sqrt{\pi}}{4A^{3/2}}\text{Erf}\left(\sqrt{A}\xi_{1}\right)-\frac{\xi_{1}}{2A}e^{-A\xi_{1}^{2}}\right],\end{split}$
(A.2)

where:

$A\equiv\frac{m_{X}v_{\text{esc}}^{2}}{12T_{X}}.$ (A.3)

With these in hand, we can utilize the general form of $n_{X}(\xi)$ in Eq.
LABEL:eq:nX-xi_FullProfile and simply apply linearity of the integral since
the Knudsen transition function, $f(\text{Kn})$, is radially independent:

$\int dVn_{X}(\xi)=f(\text{Kn})\int
dVn_{X}^{\text{LTE}}(\xi)+\left[1-f(\text{Kn})\right]\int
dVn_{X}^{\text{ISO}}(\xi).$ (A.4)

Now, the more tricky approximation occurs when we look to compute the integral
of $n_{X}^{2}(\xi)$. This is because of the cross term
$n_{X}^{\text{LTE}}n_{X}^{\text{ISO}}.$ Let us write the expansion of
$n_{X}^{2}(\xi)$:

$n_{X}^{2}(\xi)=\left[f(\text{Kn})n_{X}^{\text{LTE}}(\xi)\right]^{2}+\left[\left(1-f(\text{Kn})\right)n_{X}^{\text{ISO}}(\xi)\right]^{2}+2f(\text{Kn})\left(1-f(\text{Kn})\right)n_{X}^{\text{LTE}}(\xi)n_{X}^{\text{ISO}}(\xi).$
(A.5)

Thus, we are further interested in the following integrals:

$\int dV\left[n_{X}^{\text{LTE}}(\xi)\right]^{2},~{}\int
dV\left[n_{X}^{\text{ISO}}(\xi)\right]^{2},~{}\text{and}~{}\int
dVn_{X}^{\text{LTE}}(\xi)n_{X}^{\text{ISO}}(\xi).$

The first two are straightforward. In the LTE case, we have:

$\begin{split}\int dV\left[n_{X}^{\text{LTE}}(\xi)\right]^{2}&\simeq
4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}^{2}(0)\sqrt{\frac{27\pi}{2}}p^{-3/2},\\\
&=\frac{n_{X}(0)}{\sqrt{8}}\int dVn_{X}^{\text{LTE}}(\xi).\end{split}$ (A.6)

In the isothermal limit, the integral is:

$\int dV\left[n_{X}^{\text{ISO}}(\xi)\right]^{2}\simeq
4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}^{2}(0)\left[\frac{\sqrt{2\pi}}{16A^{3/2}}\text{Erf}\left(\sqrt{2A}\xi_{1}\right)-\frac{\xi_{1}}{4A}e^{-2A\xi_{1}^{2}}\right].$
(A.7)

The cross-term integral can, remarkably, be solved analytically using the same
approximation techniques as above. The result is as follows:

$\begin{split}\int
dVn_{X}^{\text{LTE}}(\xi)n_{X}^{\text{ISO}}(\xi)&=4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}^{2}(0)\int_{0}^{\xi_{1}}d\xi~{}\xi^{2}\theta^{p}(\xi)\exp\left\\{-\frac{m_{X}v_{\text{esc}}^{2}}{2T_{X}}\left(1-\theta(\xi)\right)\right\\},\\\
&\simeq
4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}^{2}(0)\int_{0}^{\sqrt{6}}d\xi~{}\xi^{2}\left(1-\frac{1}{6}\xi^{2}\right)^{p}e^{-A\xi^{2}},\\\
&=4\pi\left(\frac{R}{\xi_{1}}\right)^{2}n_{X}^{2}(0)\sqrt{\frac{27\pi}{2}}\frac{\Gamma(p+1)}{\Gamma\left(p+\frac{5}{2}\right)}{}_{1}F_{1}\left(\frac{3}{2},~{}p+\frac{5}{2},-6A\right),\end{split}$
(A.8)

where ${}_{1}F_{1}\left(a,b,z\right)$ is the confluent hypergeometric function
of the first kind (also known as Kummer’s function), $\Gamma(z)$ is the gamma
function and we remind the reader that
$p\equiv\frac{E_{e}}{T_{c}}+\frac{3}{2}-\alpha_{0}$ and $A$ is defined in Eq.
A.3. Note that, as with the integrals of pure LTE distributions, we cut the
integral at $\xi=\sqrt{6}$ to avoid non-physical results. Again, we emphasize
the fact that, in the high-mass limit, the cored nature of the profiles
guarantees integral convergence before $\xi=\sqrt{6}$.

We now have all we need to compute the annihilation coefficient $C_{A}$ for
high-mass $2\to 2$ Superheavy DM models. Recall that the annihilation
coefficient is given by:

$C_{A}=\langle\sigma v\rangle\frac{\int dVn_{X}^{2}}{\left(\int
dVn_{X}\right)^{2}}.$ (A.9)

Thus, with the above integrals completed, one can write the annihilation
coefficient analytically, though it would certainly be very messy. One useful
way to proceed is to understand the behavior of the coefficient in the limit
of asymptotic DM mass. To do this, we shall start by looking at the asymptotic
behavior of the integrals computed above. By doing so, we will actually be
able to show that the asymptotic behavior is the same regardless of the
Knudsen number. Starting with the LTE profile, first recall that for large
$m_{X}$, $p\sim m_{X}$. Then, if we look at Eq. A.1 and divide out $n_{X}(0)$,
our only mass-dependence comes from $p$, so we have:

$\int dV\frac{n_{X}^{\text{LTE}}(\xi)}{n_{X}(0)}\sim p^{-3/2}\sim
m_{X}^{-3/2}.$

Furthermore, since the integral of $\left[n_{X}^{\text{LTE}}\right]^{2}$
scales like the integral of $n_{X}^{\text{LTE}}$ as per Eq. A.6, if we divide
Eq. A.6 by $n_{X}^{2}(0)$, we get the same scaling relation as above:

$\int dV\left[\frac{n_{X}^{\text{LTE}}(\xi)}{n_{X}(0)}\right]^{2}\sim
p^{-3/2}\sim m_{X}^{-3/2}.$

In the isothermal case, it is first worth recalling that as $m_{X}\to\infty$,
$T_{X}\to T_{c}$, in other words the temperature loses its mass dependence.
With this, we see that in the limit of $m_{X}\to\infty$, $A\sim
m_{X}\to\infty$. Recalling the following limit:

$\lim_{x\to\infty}\text{Erf}(x)=1,$

we get the following from Eq. A.2:

$\lim_{m_{X}\to\infty}\int dV\frac{n_{X}^{\text{ISO}}(\xi)}{n_{X}(0)}\sim
A^{-1}\left[\frac{\sqrt{\pi}}{4A^{1/2}}-\frac{\xi_{1}}{2}e^{-A\xi_{1}^{2}}\right].$

It is then clear that only the first term survives since the exponential term
will decay far more quickly than the $A^{-1/2}$ term in the limit
$m_{X}\to\infty$ (or equivalently $A\to\infty$). Thus, we obtain the following
scaling relation for the integral of the isothermal profile:

$\int dV\frac{n_{X}^{\text{ISO}}(\xi)}{n_{X}(0)}\sim A^{-3/2}\sim
m_{X}^{-3/2},$

Given the similarities of Eq. A.2 and Eq. A.7, we can use an identical
argument to show:

$\int dV\left[\frac{n_{X}^{\text{ISO}}(\xi)}{n_{X}(0)}\right]^{2}\sim
m_{X}^{-3/2}.$

Thus, both limits have the same scaling with mass in the high mass limit.
Intriguingly, the integral with the cross term also exhibits this behavior.
First, we can expand Kummer’s function, ${}_{1}F_{1}\left(a,b,z\right)$, in
the limit that $z\to-\infty$. This is given by [63]:

$\lim_{z\to-\infty}{}_{1}F_{1}\left(a,b,z\right)=\frac{\Gamma(b)}{\Gamma(b-a)}(-z)^{-a}.$
(A.10)

Applying this to our cross term integral, we find:

$\int
dV\frac{n_{X}^{\text{LTE}}(\xi)n_{X}^{\text{ISO}}(\xi)}{n_{X}^{2}(0)}\sim\frac{\Gamma(p+1)}{\Gamma\left(p+\frac{5}{2}\right)}\frac{\Gamma\left(p+\frac{5}{2}\right)}{\Gamma(p+1)}(6A)^{-3/2}\sim
m_{X}^{-3/2}.$ (A.11)

With these scalings in mind, it is straightforward to show that the asymptotic
behavior of $C_{A}$ is given by:

$C_{A}\sim\langle\sigma v\rangle m_{X}^{3/2}.$

With this, it is then possible to, for example, determine how
$N_{X}(t\gg\tau_{\text{eq}})$ behaves in different regions of $\sigma-m_{X}$
parameter space, or determine how $n_{X}(0)$ behaves asymptotically with mass.
These integrals are thus very useful in understanding the high-mass behavior
of captured DM distributions and also allow for much easier computations in
this limit as it avoids the difficulties of numerical integration for these
very cored profiles.

##  Appendix B Deriving analytic capture rates

In this appendix, we outline how we go about arriving at closed-form
expressions for the DM capture rate in regions of $\sigma-m_{X}$ parameter
space. Note that most of the difficulty in approximating the capture rate
occurs in the multiscatter regime, when $\tau\gg 1$. The reason for this is
that when $\tau\ll 1$, the partial capture rate for $N>1$ is suppressed
significantly, so $C_{tot}\simeq C_{1}$. Thus, we start with the more
difficult task of approximating the multiscatter capture rate. With this in
mind, we start with the general multiscatter capture rate, which we reproduce
here for convenience:

$C_{N}=\frac{1}{3}\pi
R^{2}p_{N}(\tau)\frac{\sqrt{6}n_{X}}{\sqrt{\pi}\bar{v}}\left(\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)-\left(2\bar{v}^{2}+3v_{N}^{2}\right)\exp\left(-\frac{3\left(v_{N}^{2}-v_{\text{esc}}^{2}\right)}{2\bar{v}^{2}}\right)\right),$
(B.1)

where we remind the reader of the following definitions:

$p_{N}(\tau)=\frac{2}{\tau^{2}}\left(N+1-\frac{\Gamma(N+2,\tau)}{N!}\right),$
$v_{N}=v_{\text{esc}}\left(1-\langle z\rangle\beta_{+}\right)^{-N/2},$
$\beta_{+}=\frac{4mm_{X}}{\left(m+m_{X}\right)^{2}},$
$k\equiv\frac{3v_{\text{esc}}^{2}}{2\bar{v}^{2}}\langle
z\rangle\beta_{+}\equiv\alpha\beta_{+}.$

When considering DM-nucleon scattering cross sections that are independent of
the scattering angle, we can typically make the approximation $\langle
z\rangle\approx 1/2$. We also have, for all DM masses, $\beta_{+}\leq 1$,
which means that we can approximate $v_{N}$ in the following way for $\langle
z\rangle\beta_{+}N_{\text{max}}\ll 1$:

$v_{N}\approx v_{\text{esc}}\left(1+\frac{1}{2}\langle
z\rangle\beta_{+}N\right).$ (B.2)

The derivations for the capture rates in most of $\sigma-m_{X}$ space relies
on this approximation. However, it is very crucial to note that this
approximation will not always be valid when $\alpha\ll 1$, and hence the
emergence of a region without a closed form solution in Fig. 1. We will
discuss this more below when we derive the capture rates in different regions.
Next, we can approximate the probability function $p_{N}(\tau)$ in the
$\tau\gg 1$ limit:

$p_{N}(\tau)\approx\frac{2}{\tau^{2}}\left(N+1\right)\Theta\left(\tau-N\right).$
(B.3)

With these approximations, we can rewrite the partial capture rate as:

$C_{N}\simeq\frac{2A_{X}}{\tau^{2}}\Theta\left(\tau-N\right)\left[\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)-\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)e^{-kN}-2\bar{v}^{2}kNe^{-kN}\right],$
(B.4)

where

$A_{X}\equiv\frac{1}{3}\pi R^{2}\sqrt{\frac{6}{\pi}}\frac{n_{X}}{\bar{v}}.$

The presence of the Heaviside step function means that we need only sum
$C_{N}$ up to $\tau$. Thus, we define the following sums for notational
convenience:

$B_{1}\equiv\sum_{N=1}^{\tau}(N+1),$ (B.5)
$B_{2}\equiv\sum_{N=1}^{\tau}(N+1)e^{-kN},$ (B.6)
$B_{3}\equiv\sum_{N=1}^{\tau}N(N+1)e^{-kN}.$ (B.7)

The total capture rate can then be expressed succinctly:

$C_{tot}\approx\frac{2A_{X}}{\tau^{2}}\left[\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)\left(B_{1}-B_{2}\right)-2\bar{v}^{2}kB_{3}\right].$
(B.8)

These sums can be closed and are given by:

$\displaystyle B_{1}$ $\displaystyle=\frac{1}{2}\tau(\tau+3),$ (B.9)
$\displaystyle B_{2}$
$\displaystyle=\frac{e^{-k\tau}\left(1+e^{k\tau}\left(2e^{k}-1\right)-e^{k}\left(\tau+2\right)+\tau\right)}{\left(e^{k}-1\right)^{2}},$
(B.10) $\displaystyle B_{3}$
$\displaystyle=\frac{e^{-k\tau}\left(-e^{2k}(\tau+2)(\tau+1)+2e^{k(\tau+2)}+2e^{k}\tau(\tau+2)-\tau(\tau+1)\right)}{\left(e^{k}-1\right)^{3}}.$
(B.11)

We thus have completely determined the total capture rate in the multiscatter
regime with two assumptions: $\langle z\rangle\beta_{+}N_{\text{max}}<1$ and
$\tau\gg 1$. Note that, since the sums are being truncated at $\tau$, the
first condition is the same as $\langle z\rangle\beta_{+}\tau<1$. However, in
the above form, the capture rate is difficult to compute, largely because it
involves exponentiation of very large numbers. Thus, we can make further
assumptions, effectively dividing the $\sigma-m_{X}$ DM parameter space into
different regions, which is what is presented in Section 2.1.

There are two distinct cases we focus on in the multiscattering regime:
$\alpha>1$ and $\alpha<1$. In both cases we see the emergence of Regions I and
II, however there are important differences between the two and, in fact,
there emerges a region in the $\alpha<1$ regime where both approximations
break down, as mentioned above when discussing the approximation of $v_{N}$ in
Eq. B.2. We will start with the $\alpha>1$ case. The transition here surrounds
the condition $k\tau\sim 1$. We define Region I as when $k\tau\ll 1$ and
$\tau\gg 1$ in this case. Before proceeding further, we will show that the
approximation of $v_{N}$ will always be true in this regime. That
approximation hinged on the fact that $\langle z\rangle\beta_{+}\tau\ll 1$.
Note that this constraint can be rewritten as the following:

$\frac{\langle z\rangle k\tau}{\alpha}\ll 1.$ (B.12)

In this form, it is clear that when $k\tau\ll 1$ and $\alpha\geq 1$, this
condition is satisfied. We also see how we might run into trouble when
$\alpha<1$. Now that we have verified the validity of the approximation we
started with, we can proceed to expand around $k\tau\ll 1$ and $k\ll 1$. The
reason we know $k\ll 1$ is because $\tau\gg 1$ and $k\tau\ll 1$. One can
verify that by expanding they get the Region I capture rate in Eq. 2.12. To
find the Region II capture rate, one can expand around $k\tau\gg 1$ and
recover Eq. 2.13. A keen reader may object since we are no longer guaranteed
that the approximation of $v_{N}$ holds. While this may be the case, we find
numerically that a transition to the Region II capture rate, which is, after
all, the geometric capture rate, does occur at $k\tau=1$ (See Fig. 2).

Now we tackle the case when $\alpha<1$. The key difference here is that there
is no longer a transition line at $k\tau=1$ between Regions I and II, but
rather a transition region. To see why, let’s re-examine Eq. B.12. Here, we
must re-define the border of Region I to be when this approximation breaks
down. However, this is not the same as when $k\tau=1$, since, atleast for
$\alpha<1/2$, when $k\tau=1$, the condition is violated. So the region of
transition occurs between the lines $k\tau=1$ and $\frac{\langle z\rangle
k\tau}{\alpha}=1$. We do admit here that there is still a puzzle as to where
exactly the transition into Region II occurs. We know that we can recover the
Region II capture rate for when $k\tau\gg 1$, and this gives us an upper bound
on the Region II boundary. However, numerical explorations, such as in Fig. 2,
seem to suggest that the Region II boundary occurs somewhere below the
$k\tau=1$ line but above $\frac{\langle z\rangle k\tau}{\alpha}=1$. We thus
actually have a lower and upper bound on the Region II boundary. We do also
want to point out that, even though we have yet to determine that boundary
analytically, it is straightforward to use numerical techniques to estimate
it.

Now we discuss the analytic approximations in the single-scatter regime, i.e.
when $\tau\ll 1$. As with the multiscatter case, there are differences between
the $\alpha>1$ and $\alpha<1$ scenarios. Specifically, Region III only emerges
in the case where $\alpha>1$. This is because when $\alpha>1$, there are two
regions to consider: $k>1$ and $k<1$. However, if $\alpha<1$, then $k<1$
always and there is only need for Region IV. To arrive at the capture rate in
these regions, we make the approximation $C_{\text{tot}}\approx C_{1}$. This
is valid because of the following expansion for $p_{N}(\tau)$ when $\tau\ll
1$:

$p_{N}(\tau)\approx\frac{2\tau^{N}}{N!(N+2)},$ (B.13)

which significantly suppresses $C_{N}$ for any $N>1$. Note that the
approximation of $v_{N}$ in Eq. B.2 is always valid here since:

$\langle z\rangle\beta_{+}N_{\text{max}}=\langle z\rangle\beta_{+}<1.$

We can thus write the total capture rate in the single-scatter regime as:

$C_{tot}\simeq\frac{2}{3}A_{X}\tau\left[\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)-\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)e^{-k}-2\bar{v}^{2}ke^{-k}\right].$
(B.14)

Then, by expanding for $k\gg 1$ and $k\ll 1$, we get the capture rates Region
III and IV respectively (Eqs. 2.14 and 2.15).

##  Appendix C Estimating $T_{\text{int}}$, i.e. the contribution of internal
heating to Jupiters’ and Brown Dwarfs’ surface temperatures

In this section, we outline how we estimate the temperature of the internal
heating component ($T_{\text{int}}$ in Eq. 3.17) of Jupiters and Brown Dwarfs.
We also derive the minimum Jupiter mass for exoplanet detection in
consideration of the critical temperature $T_{\text{crit}}$ introduced in
Section 3.3. In Figure 5 of [59], evolutionary models are used to find the
cooling curves of Jupiters and Brown Dwarfs. Since we are interested in
exoplanets heated by DM in the late stages of their life, after they have
cooled sufficiently from internal heat sources, we estimate the internal
heating component of our objects of interest towards the end of their cooling
curves, at an age of $\sim 1\,\mathrm{Gyr}$. To find the temperature
contribution from internal heat sources, we fit the points at the end of the
cooling curve in Figure 5, where each point corresponds to a different object
mass. The luminosity at these points corresponds to a time $t\sim
1\,\mathrm{Gyr}$.

Figure 7: Interpolation of the temperature contribution to internal heating
($T_{\text{int}}$ in Eq. 3.17) for Jupiters and Brown Dwarfs as a function of
their mass. The “Data” represents the temperature from the cooling curves of
[59] at time $\sim 1\,\mathrm{Gyr}$.

In Figure 7, we show the data points from Figure 5 of [59], as well as an
interpolation of this data for all Jupiter and Brown Dwarf Masses.

We now go on to compute the minimum Jupiter mass to be considered as useful
probes of DM for JWST searches given the critical temperature discussed in
Section 3.3. As discussed there, we only consider exoplanets that can be
feasibly detected by JWST in a reasonable exposure time. This implies a
minimum temperature $T_{\text{min}}=650\,\mathrm{K}$. We then only consider
those objects with critical temperature $T_{\text{crit}}>T_{\text{min}}$, else
we are not sensitive to DM signals. We thus have:

$\begin{split}T_{\text{crit}}&>T_{\text{min}},\\\
\left(\frac{1}{\sqrt{24\pi}}\frac{\rho_{X}}{\sigma_{\text{SB}}\bar{v}}\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)+T_{\text{int}}^{4}\right)^{\frac{1}{4}}&>T_{\text{min}},\\\
\frac{1}{\sqrt{24\pi}}\frac{\rho_{X}}{\sigma_{\text{SB}}\bar{v}}\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)+T_{\text{int}}^{4}&>T_{\text{min}}^{4},\\\
T_{\text{int}}&>\left(T_{\text{min}}^{4}-\frac{1}{\sqrt{24\pi}}\frac{\rho_{X}}{\sigma_{\text{SB}}\bar{v}}\left(2\bar{v}^{2}+3v_{\text{esc}}^{2}\right)\right)^{\frac{1}{4}}.\end{split}$
(C.1)

Note that the object-dependence enters in both $T_{\text{int}}$ and
$v_{\text{esc}}$. One can verify that the smallest mass satisfying this
equation is $M\approx 3M_{J}$.

## References

  * [1] F. Zwicky, _Die Rotverschiebung von extragalaktischen Nebeln_ , _Helvetica Physica Acta_ 6 (1933) 110.
  * [2] V. C. Rubin and J. Ford, W. Kent, _Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions_ , _ApJ_ 159 (1970) 379.
  * [3] J. A. Tyson, G. P. Kochanski and I. P. Dell'Antonio, _Detailed mass map of CL 0024 $+$1654 from strong lensing_, _The Astrophysical Journal_ 498 (1998) L107.
  * [4] D. Clowe, M. Bradač, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones et al., _A direct empirical proof of the existence of dark matter_ , _The Astrophysical Journal_ 648 (2006) L109.
  * [5] WMAP collaboration, _Five-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation_ , _Astrophys. J. Suppl._ 180 (2009) 330 [0803.0547].
  * [6] E. Komatsu, K. M. Smith, J. Dunkley, C. L. Bennett, B. Gold, G. Hinshaw et al., _Seven-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Interpretation_ , _ApJS_ 192 (2011) 18 [1001.4538].
  * [7] Planck collaboration, _Planck 2015 results. XIII. Cosmological parameters_ , _Astron. Astrophys._ 594 (2016) A13 [1502.01589].
  * [8] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _A &A_ (2018) [1807.06209].
  * [9] R. N. Mohapatra and V. L. Teplitz, _Mirror dark matter_ , in _Sources and Detection of Dark Matter and Dark Energy in the Universe_ , D. B. Cline, ed., (Berlin, Heidelberg), pp. 211–221, Springer Berlin Heidelberg, 2001.
  * [10] S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * [11] G. Servant and T. M. Tait, _Is the lightest kaluza–klein particle a viable dark matter candidate?_ , _Nuclear Physics B_ 650 (2003) 391.
  * [12] K. Freese, _Status of Dark Matter in the Universe_ , _Int. J. Mod. Phys._ 1 (2017) 325 [1701.01840].
  * [13] T. Braine, R. Cervantes, N. Crisosto, N. Du, S. Kimes, L. Rosenberg et al., _Extended search for the invisible axion with the axion dark matter experiment_ , _Physical Review Letters_ 124 (2020) .
  * [14] XENON Collaboration 4 collaboration, _Constraining the spin-dependent wimp-nucleon cross sections with xenon1t_ , _Phys. Rev. Lett._ 122 (2019) 141301.
  * [15] PICO collaboration, _Dark Matter Search Results from the Complete Exposure of the PICO-60 C 3F8 Bubble Chamber_, _Phys. Rev. D_ 100 (2019) 022001 [1902.04031].
  * [16] A. R. Liddle and D. H. Lyth, _Cosmological Inflation and Large-Scale Structure_. Cambridge University Press, 2000, 10.1017/CBO9781139175180.
  * [17] V. Bromm and R. B. Larson, _The First stars_ , _Ann. Rev. Astron. Astrophys._ 42 (2004) 79 [astro-ph/0311019].
  * [18] V. Bromm, N. Yoshida, L. Hernquist and C. F. McKee, _The formation of the first stars and galaxies_ , _Nature_ 459 (2009) 49 [0905.0929].
  * [19] J. Wang, S. Bose, C. S. Frenk, L. Gao, A. Jenkins, V. Springel et al., _Universal structure of dark matter haloes over a mass range of 20 orders of magnitude_ , _Nature_ 585 (2020) 39.
  * [20] A. Gould, _Direct and Indirect Capture of Wimps by the Earth_ , _Astrophys. J._ 328 (1988) 919.
  * [21] J. Bramante, A. Delgado and A. Martin, _Multiscatter stellar capture of dark matter_ , _Phys. Rev. D_ 96 (2017) 063002 [1703.04043].
  * [22] C. Ilie, J. Pilawa and S. Zhang, _Comment on “multiscatter stellar capture of dark matter”_ , _Phys. Rev. D_ 102 (2020) 048301.
  * [23] B. Dasgupta, A. Gupta and A. Ray, _Dark matter capture in celestial objects: Improved treatment of multiple scattering and updated constraints from white dwarfs_ , _JCAP_ 08 (2019) 018 [1906.04204].
  * [24] M. Baryakhtar, J. Bramante, S. W. Li, T. Linden and N. Raj, _Dark Kinetic Heating of Neutron Stars and an Infrared Window on WIMPs, SIMPs, and Pure Higgsinos_ , _Phys. Rev. Lett._ 119 (2017) 131801 [1704.01577].
  * [25] R. Garani, A. Gupta and N. Raj, _Observing the thermalization of dark matter in neutron stars_ , _Phys. Rev. D_ 103 (2021) 043019 [2009.10728].
  * [26] K. Freese, D. Spolyar and A. Aguirre, _Dark Matter Capture in the first star: a Power source and a limit on Stellar Mass_ , _JCAP_ 0811 (2008) 014 [0802.1724].
  * [27] C. Ilie, C. Levy, J. Pilawa and S. Zhang, _Constraining dark matter properties with the first generation of stars_ , _Phys. Rev. D_ 104 (2021) 123031.
  * [28] C. Ilie and C. Levy, _Multicomponent multiscatter capture of dark matter_ , _Phys. Rev. D_ 104 (2021) 083033.
  * [29] C. J. Horowitz, _Nuclear and dark matter heating in massive white dwarf stars_ , _Phys. Rev. D_ 102 (2020) 083031 [2008.03291].
  * [30] R. K. Leane and J. Smirnov, _Exoplanets as sub-gev dark matter detectors_ , _Phys. Rev. Lett._ 126 (2021) 161101.
  * [31] D. N. Spergel and W. H. Press, _Effect of hypothetical, weakly interacting, massive particles on energy transport in the solar interior_ , _Astrophysical Journal_ 294 (1985) 663.
  * [32] A. Gould, _Weakly Interacting Massive Particle Distribution in and Evaporation from the Sun_ , _Astrophysical Journal_ 321 (1987) 560.
  * [33] R. Garani and S. Palomares-Ruiz, _Evaporation of dark matter from celestial bodies_ , _Journal of Cosmology and Astroparticle Physics_ 2022 (2022) 042.
  * [34] S. Chandrasekhar, _An introduction to the study of stellar structure_. 1939.
  * [35] R. Helled, J. D. Anderson, M. Podolak and G. Schubert, _INTERIOR MODELS OF URANUS AND NEPTUNE_ , _The Astrophysical Journal_ 726 (2010) 15.
  * [36] M. S. Marley and J. J. Fortney, _Chapter 33 - interiors of the giant planets_ , in _Encyclopedia of the Solar System (Third Edition)_ , T. Spohn, D. Breuer and T. V. Johnson, eds., (Boston), pp. 743–758, Elsevier, (2014), DOI.
  * [37] D. J. Stevenson, _Jupiter’s interior as revealed by juno_ , _Annual Review of Earth and Planetary Sciences_ 48 (2020) 465.
  * [38] G. Chabrier and I. Baraffe, _Theory of low-mass stars and substellar objects_ , _Annual Review of Astronomy and Astrophysics_ 38 (2000) 337–377.
  * [39] J. Bramante, J. Kumar, G. Mohlabeng, N. Raj and N. Song, _Light dark matter accumulating in terrestrial planets: Nuclear scattering_ , .
  * [40] R. K. Leane and J. Smirnov, _Floating dark matter in celestial bodies_ , .
  * [41] J. Smirnov and J. F. Beacom, _New freezeout mechanism for strongly interacting dark matter_ , _Phys. Rev. Lett._ 125 (2020) 131301.
  * [42] B. W. Lee and S. Weinberg, _Cosmological Lower Bound on Heavy Neutrino Masses_ , _Phys. Rev. Lett._ 39 (1977) 165.
  * [43] A. Gould, _Big bang archeology - WIMP capture by the earth at finite optical depth_ , _Astrophysical Journal_ 387 (1992) 21.
  * [44] C. Ilie, _Closed-form expressions for multiscatter Dark Matter capture rates (In Preparation)_ , _JCAP (Intended)_ (2024) .
  * [45] T. Guillot and D. Gautier, _Giant Planets_.  12, 2009.  0912.2019.
  * [46] R. Garani and S. Palomares-Ruiz, _Dark matter in the sun: scattering off electrons vs nucleons_ , _Journal of Cosmology and Astroparticle Physics_ 2017 (2017) 007.
  * [47] H. Banks, S. Ansari, A. C. Vincent and P. Scott, _Simulation of energy transport by dark matter scattering in stars_ , _Journal of Cosmology and Astroparticle Physics_ 2022 (2022) 002.
  * [48] Y. Hochberg, E. Kuflik, T. Volansky and J. G. Wacker, _Mechanism for thermal relic dark matter of strongly interacting massive particles_ , _Phys. Rev. Lett._ 113 (2014) 171301.
  * [49] K. Griest and M. Kamionkowski, _Unitarity limits on the mass and radius of dark-matter particles_ , _Phys. Rev. Lett._ 64 (1990) 615.
  * [50] A. Cassan, D. Kubas, J.-P. Beaulieu, M. Dominik, K. Horne, J. Greenhill et al., _One or more bound planets per milky way star from microlensing observations_ , _Nature_ 481 (2012) 167.
  * [51] “Nasa exoplanets.” https://exoplanets.nasa.gov/.
  * [52] T. Currie, G. M. Brandt, T. D. Brandt, B. Lacy, A. Burrows, O. Guyon et al., _Direct imaging and astrometric detection of a gas giant planet orbiting an accelerating star_ , _Science_ 380 (2023) 198 [https://www.science.org/doi/pdf/10.1126/science.abo6192].
  * [53] A. F. Gupta, J. M. Jackson, G. Hebrard, A. S. Lin, K. G. Stassun, J. Dong et al., _A high-eccentricity warm jupiter orbiting toi-4127_ , 2023.
  * [54] S. A. Johnson, M. Penny, B. S. Gaudi, E. Kerins, N. J. Rattenbury, A. C. Robin et al., _Predictions of the nancy grace roman space telescope galactic exoplanet survey. II. free-floating planet detection rates_ , _The Astronomical Journal_ 160 (2020) 123.
  * [55] B. Paxton, L. Bildsten, A. Dotter, F. Herwig, P. Lesaffre and F. Timmes, _Modules for Experiments in Stellar Astrophysics (MESA)_ , _ApJS_ 192 (2011) 3 [1009.1622].
  * [56] S. Auddy, S. Basu and S. R. Valluri, _Analytic models of brown dwarfs and the substellar mass limit_ , _Advances in Astronomy_ 2016 (2016) 1–15.
  * [57] A. Gould and G. Raffelt, _Thermal Conduction by Massive Particles_ , _ApJ_ 352 (1990) 654.
  * [58] J. L. Feng and J. Kumar, _The WIMPless Miracle: Dark-Matter Particles without Weak-Scale Masses or Weak Interactions_ , _Phys. Rev. Lett._ 101 (2008) 231301 [0803.4196].
  * [59] J. A. Caballero, _A review on substellar objects below the deuterium burning mass limit: Planets, brown dwarfs or what?_ , _Geosciences_ 8 (2018) .
  * [60] T. Cohen, M. Lisanti, A. Pierce and T. R. Slatyer, _Wino dark matter under siege_ , _Journal of Cosmology and Astroparticle Physics_ 2013 (2013) 061.
  * [61] XENON collaboration, _Search for Light Dark Matter Interactions Enhanced by the Migdal Effect or Bremsstrahlung in XENON1T_ , _Phys. Rev. Lett._ 123 (2019) 241803 [1907.12771].
  * [62] Z. Liu, Q. Yue, L. Yang, K. Kang, Y. Li, H. Wong et al., _Constraints on spin-independent nucleus scattering with sub-GeV weakly interacting massive particle dark matter from the CDEX-1b experiment at the china jinping underground laboratory_ , _Physical Review Letters_ 123 (2019) .
  * [63] M. Abramowitz and I. Stegun, _Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables_ , Applied mathematics series. Dover Publications, 1965.