††thanks<EMAIL_ADDRESS>0000-0002-4711-6516††thanks:
<EMAIL_ADDRESS>0000-0001-9344-736X††thanks<EMAIL_ADDRESS>0000-0002-1287-8780††thanks<EMAIL_ADDRESS>0000-0002-0001-7123††thanks:
<EMAIL_ADDRESS>0000-0002-3082-0929††thanks<EMAIL_ADDRESS>0000-0002-3415-0707

# Dark Matter Halo Parameters from Overheated Exoplanets
via Bayesian Hierarchical Inference

María Benito  Tartu Observatory, University of Tartu, Observatooriumi 1,
Tõravere 61602, Estonia  Konstantin Karchev  Theoretical and Scientific Data
Science, SISSA, Trieste, Italy  Rebecca K. Leane  SLAC National Accelerator
Laboratory, Stanford University, Stanford, CA 94039, USA  Kavli Institute for
Particle Astrophysics and Cosmology, Stanford University, Stanford, CA 94039,
USA  Sven Põder  National Institute of Chemical Physics and Biophysics, Rävala
10, Tallinn 10143, Estonia  Juri Smirnov  Department of Mathematical Sciences,
University of Liverpool, Liverpool, L69 7ZL, United Kingdom  Roberto Trotta
Theoretical and Scientific Data Science, SISSA, Trieste, Italy  Imperial
College London, Physics Department, Blackett Laboratory, Prince Consort Rd,
SW7 2AZ London, United Kingdom  Italian Research Center on High Performance
Computing, Big Data and Quantum Computing & National Institute for Nuclear
Physics, Sezione di Trieste, Italy

###### Abstract

Dark Matter (DM) can become captured, deposit annihilation energy, and hence
increase the heat flow in exoplanets and brown dwarfs. Detecting such a DM-
induced heating in a population of exoplanets in the inner kpc of the Milky
Way thus provides potential sensitivity to the galactic DM halo parameters. We
develop a Bayesian Hierarchical Model to investigate the feasibility of DM
discovery with exoplanets and examine future prospects to recover the spatial
distribution of DM in the Milky Way. We reconstruct from mock exoplanet
datasets observable parameters such as exoplanet age, temperature, mass, and
location, together with DM halo parameters, for representative choices of
measurement uncertainty and the number of exoplanets detected. We find that
detection of $\mathcal{O}(100)$ exoplanets in the inner Galaxy can yield
quantitative information on the galactic DM density profile, under the
assumption of 10% measurement uncertainty. Even as few as $\mathcal{O}(10)$
exoplanets can deliver meaningful sensitivities if the DM density and inner
slope are sufficiently large.

††preprint: SLAC-PUB-17772††preprint: LTH-1368

The milky way is estimated to be home to about 100 billion exoplanets and
brown dwarfs Mužić  _et al._ (2017). Each of them holds a wealth of data,
waiting to reveal their secrets and curiosities. One such secret is whether or
not exoplanets and brown dwarfs (hereafter collectively referred to as
“exoplanets”) harbor dark matter (DM) cores and consequently have temperatures
much higher than otherwise expected Leane and Smirnov (2021).

DM heating in exoplanets is expected to potentially arise in the following
way. Exoplanets sweep through the Galactic DM halo, collecting DM which has
scattered with the exoplanet particles and lost sufficient kinetic energy to
become bound to the exoplanet. The accumulated DM can then potentially
annihilate, depositing its annihilation energy into the exoplanet and
increasing its temperature beyond that expected from models of exoplanet
formation and evolution. Once sufficient DM has accumulated, the exoplanet can
reach equilibrium between DM scattering and annihilation, such that the
measured temperature provides a probe of the DM scattering rate and enables
inference of the DM density near the exoplanet, as the "excess" heating scales
with the DM density (see Figure 1). This signal can be potentially detected by
upcoming infrared telescopes, such as the James Webb Space Telescope (JWST),
or Roman Telescope Leane and Smirnov (2021).

Figure 1: Schematic of DM heated exoplanets. In the inner Galaxy, there is
more DM, leading to higher levels of exoplanet heating. Therefore, measuring
exoplanet temperatures as a function of Galactocentric radius provides a probe
of the DM density distribution in our Galaxy.

DM heating processes have also been considered in other celestial objects,
such as nuclear-burning stars, neutron stars and white dwarfs Goldman and
Nussinov (1989); Gould  _et al._ (1990); Kouvaris (2008); Bertone and
Fairbairn (2008); Spolyar  _et al._ (2008); Iocco (2008); Freese  _et al._
(2009); Taoso  _et al._ (2008); Sivertsson and Gondolo (2011); Freese  _et
al._ (2016); Salati and Silk (1989); Fairbairn  _et al._ (2008); Scott  _et
al._ (2009); Lopes and Lopes (2021); de Lavallaz and Fairbairn (2010);
Kouvaris and Tinyakov (2010); McDermott  _et al._ (2012); Kouvaris and
Tinyakov (2011a); Guver  _et al._ (2014); Bramante  _et al._ (2013); Bell  _et
al._ (2013); Bramante  _et al._ (2014); Bertoni  _et al._ (2013); Kouvaris and
Tinyakov (2011b); McCullough and Fairbairn (2010); Angeles Perez-Garcia and
Silk (2015); Bramante (2015); Graham  _et al._ (2015); Cermeno  _et al._
(2016); Krall and Reece (2018); Leane  _et al._ (2017); McKeen  _et al._
(2018); Baryakhtar  _et al._ (2017); Raj  _et al._ (2018); Bell  _et al._
(2018); Chen and Lin (2018); Hamaguchi  _et al._ (2019); Camargo  _et al._
(2019); Bell  _et al._ (2019); Garani and Heeck (2019); Acevedo  _et al._
(2020); Joglekar  _et al._ (2020a, b); Bell  _et al._ (2020); Garani  _et al._
(2021); Leane  _et al._ (2021); Acevedo  _et al._ (2024a); John  _et al._
(2023); Croon and Sakstein (2023); Acevedo  _et al._ (2024b), the Earth Mack
_et al._ (2007); Chauhan and Mohanty (2016); Bramante  _et al._ (2020), and
other planets and moons Adler (2009); Kawasaki  _et al._ (1992); Mitra (2004);
Bramante  _et al._ (2020); Garani and Tinyakov (2020); Chan and Lee (2020);
Leane and Linden (2023); Blanco and Leane (2023); Croon and Smirnov (2023);
Linden  _et al._ (2024). Exoplanets can be superior probes to these other
celestial objects in several ways Leane and Smirnov (2021); Leane and Tong
(2024). First, they have cooler cores while having relatively high escape
velocities, allowing lighter DM particles to not evaporate and be retained in
the system, providing a probe of sub-GeV mass DM. Second, they have large
radii, allowing them to be detectable far into the inner Galaxy. Third,
exoplanets are highly abundant, such that this signature can be probed
throughout the Galaxy with potentially high statistics. This allows for a
location-dependent probe of the DM density profile.

The exoplanet program is rapidly accelerating. Many new exoplanets are
expected to be identified soon, and many new measurements will be performed.
The Vera C. Rubin Observatory (Rubin), the Nancy Grace Roman Space Telescope
(Roman), and the Gaia Spacecraft have or will have targeted programs to
discover as many exoplanets as possible Perryman  _et al._ (2014); Green  _et
al._ (2012); Johnson  _et al._ (2020). While exoplanets have not yet had their
spectra measured at the low temperatures and large distances required for this
search, the technological breakthroughs of these instruments will allow
unprecedented sensitivities to distant exoplanets. There are also many surveys
such as the Optical Gravitational Lensing Experiment (OGLE), Two Micron All
Sky Survey (2MASS), and the Wide-field Infrared Survey Explorer (WISE), which
peer deep into our Galaxy. Detecting the temperature of either free-floating
exoplanets, or brown dwarfs, first requires target detection via gravitational
microlensing, which can be aided especially with simultaneous use of
telescopes. This has already been investigated for example with Roman and
Euclid Bachelet and Penny (2019). Microlensing has so far led to the
identification of candidates of free floating exoplanets and brown dwarfs at
large distances, for example OGLE-2015-BLG-1268, a brown dwarf candidate with
50 Jupiter masses and at $5.9\pm 1.0\text{ kpc}$ Zhu  _et al._ (2016).

Overall, the interplay of all these telescopes will only improve the prospects
for this search, making it timely to consider the precision required of
exoplanet-related measurements to be used in inference of DM distribution and
properties. Setting up a first statistical framework to attack this problem,
and quantifying the measurement uncertainties required alongside the number of
required target exoplanets to test the DM halo, is the goal of this work.

In this work, we develop a Bayesian hierarchical model (BHM) that takes as
input measured properties of exoplanets and infers simultaneously: parameters
describing their population, such as the shape of their initial mass function
and Galactic radial number density profile; the individual properties of all
objects; and parameters of the Galactic DM density profile. We use this model
because our problem contains many correlated uncertainties, and the model is
particularly powerful at accounting for the variability within the exoplanet
population and measurement uncertainties. We use this Bayesian hierarchical
model to investigate the sample size and measurement precision required for
inference of DM halo parameters.

This paper is organized as follows. In Section I, we describe our statistical
model and provide details of the luminosities of exoplanets, including DM
capture and heating. We discuss our results, including the sensitivity to DM-
heating detectability, in Section II. Finally, we conclude and discuss the
outlook for detecting DM-heated exoplanets in Section III.

##  I The Bayesian Hierarchical Model

The interplay between the latent ($i.e.$, unobserved) distribution of
exoplanet parameters and noise arising from observational error leads to a
well-known statistical effect that astronomers call “Eddington bias” Eddington
(1913) – the fact that the parameters of the underlying population will be
incorrectly estimated from the observed sample (for an illustration of how and
why this arises in Gaussian linear models, see Refs. Kelly (2007); March  _et
al._ (2011)). To counter the bias, it is necessary to model the process
leading from population parameters to observed (noisy) samples. This can be
achieved by means of a hierarchical structure, which models separately the
population distribution and the noise. Conditioning on the observed data then
leads to posterior distributions for all – global and latent – parameters in
the model that correctly capture the different properties at each level where
uncertainties are introduced (either through intrinsic variability at the
population level or through noise). In other words, this requires a Bayesian
hierarchical model of the forward generation of the data.

In this spirit, we present in this section a Bayesian hierarchical model for
DM-heated exoplanets, which allows one to infer the normalisation and slope of
the Galactic DM halo density profile jointly with parameters describing the
population of exoplanets and the a-priori uncertain properties of individual
objects from noisy measurements of their temperatures, masses, ages, and
positions within the Galaxy.

Figure 2 shows a graphical network which illustrates the structure of this
model, and is summarised in Table LABEL:tab:BHM. We proceed to elaborate it
below.

Parameter |  | Prior | Mock value
---|---|---|---
Measured age, mass, | $\hat{X}$ | $\mathcal{N}(X,(\sigma\cdot X)^{2})$ |
Gal. distance & temp. |  |  |
Mean | $X$ |  |
Uncertainty | $\sigma\cdot X$ |  | $\sigma\in[1,20]\%$
Latent age | $A$ | $\mathcal{U}(a,b)$ |
Min. age | a |  | 8 Gyr
Max. age | b |  | 10 Gyr
Latent mass | $M$ | $\mathcal{P}(\gamma_{M},M_{\rm min})$ | $\in[14,55]\,\rm M_{\rm jup}$
Slope | $\gamma_{M}$ | $\mathcal{N}(0.6,0.1^{2})$ | 0.6
Min. mass | $M_{\rm min}$ |  | $14\,\rm M_{\rm jup}$
Latent Gal. dist. | $R$ | $\mathcal{E}(\gamma_{R})$ | $\in[0.1,1]\,\rm kpc$
Slope | $\gamma_{R}$ | $\mathcal{U}(1,2)$ | 1.46
Latent temp. | $T$ |  |
DM dens. slope | $\alpha$ | $\mathcal{U}(0,3)$ | $\in[0,3]$
DM dens. norm. | $C$ | $\mathcal{U}(0.5,40)$ | $6,20\,\rm GeV/cm^{3}$
DM vel. dispersion | $\sigma_{\rm DM}$ |  | $100\,\rm km/s$
Table 1: Exoplanet parameters, priors and fiducial values used to generate
mock observations (see Subsection I.5).  Figure 2: Directed acyclic graph
depicting our Bayesian hierarchical model used to analyse overheated
exoplanets and properties of their population and of the DM density profile.
Deterministic and probabilistic connections between variables are indicated by
solid and dotted arrows, respectively. Within the $N$ plate, representing
enumeration over N objects, data are shown in green, while independent and
dependent (unobserved, latent) parameters are denoted by purple and grey
nodes, respectively. Global parameters are depicted in blue, with those held
fixed in this study shown in pink.

###  I.1 Observables

The observed data $\bm{d}$ consists of noisy measurements $\hat{A}_{i}$,
$\hat{M}_{i}$, $\hat{R}_{i}$, and $\hat{T}_{i}$ ($i=1,\dots,N$) of the ages,
masses, Galactocentric distances, and temperatures, respectively, of a sample
of $N$ exoplanets; the data are shown in green in the Directed Acyclic Graph
(DAG) of Fig. 2. Under the assumption of normally distributed errors, the
sampling distribution of $\bm{d}\in\mathbb{R}^{4N}$ is Gaussian:

$\bm{d}\sim\mathcal{N}([\bm{A},\bm{M},\bm{R},\bm{T}],\hat{\bm{\Sigma}}),$ (1)

where $\bm{A},\bm{M}$ and $\bm{R}$ are $\mathbb{R}^{N}$ vectors containing the
latent (unobserved, true) age $A_{i}$, mass $M_{i}$, and Galactocentric
distance $R_{i}$, respectively, of individual exoplanets, while
$\bm{T}\in\mathbb{R}^{N}$ is the latent temperature of each exoplanet, which
is a deterministic function of the other exoplanet parameters and of
properties of the DM halo (see Eq. 4). Latent variables are depicted in the
DAG as purple and grey nodes, the latter indicating that they are obtained
from other variables via a deterministic relationship (dotted arrows). The
observational covariance matrix $\hat{\bm{\Sigma}}\in\mathbb{R}^{4N\times 4N}$
accounts for correlations between measurements of individual exoplanets and
among different exoplanets. We will assume for simplicity that all
measurements are independent, which means that $\hat{\bm{\Sigma}}$ reduces to
a diagonal matrix. In order to explore the effect of measurement uncertainty
on our results, we will adopt different noise levels, expressed through the
relative uncertainty $\sigma$ in measurements of each quantity
$X\in\\{A,M,R,T\\}$, as detailed in Table LABEL:tab:BHM and Subsection I.5.
Thus, the covariance matrix takes the form:

$\hat{\bm{\Sigma}}=\text{diag}(\sigma^{2}\cdot A^{2}_{i},\sigma^{2}\cdot
M^{2}_{i},\sigma^{2}\cdot R^{2}_{i},\sigma^{2}\cdot T^{2}_{i}),\,i=1,\dots,N.$

###  I.2 Heat Model for Exoplanets

We now discuss the heat contributions to exoplanets, both from DM and internal
background processes. Without DM, exoplanets have an effective temperature as
a function of their age $A_{i}$ and mass $M_{i}$. This heat occurs from $e.g.$
burning processes in the object, as well as heat left over from formation. For
the internal heat without DM, $\Gamma^{\rm int}_{{\rm heat},i}(A_{i},M_{i})$,
we use the theoretical evolution model Atmo2020 Phillips  _et al._ (2020).

The amount of DM heating in an exoplanet depends in part on the flux of DM
particles available to be captured. This is characterized by the phase-space
distribution of DM particles, which provides a particular DM flux at a given
location in the Galaxy. The most widely adopted generalized Navarro-Frenk-
White (gNFW) density profile (see Eq. 13) exhibits a 2D degeneracy between the
scale radius $r_{s}$ and the inner slope $\gamma$ Karukes  _et al._ (2019).
When the Galactocentric position is smaller than $r_{s}$, this profile
approximates a power-law form (see App. B for details). Since we are
interested in the inner 1 kpc of the Galaxy and values of $r_{s}$ less than 10
kpc are disfavoured by observations of the Milky Way’s rotation curve Benito
_et al._ (2019, 2021), we adopt a power-law density profile given by

$\rho_{\rm DM}(R;\alpha,C)=C\left(\frac{R}{1\,\rm kpc}\right)^{-\alpha},$ (2)

where $C$ denotes the DM density at a Galactocentric radius of 1 kpc.

For the DM velocities, we assume that the DM particles follow an isotropic
velocity profile described by a Maxwell-Boltzmann distribution, with an
averaged velocity $\bar{v}_{\rm DM}(R)$ related to the velocity dispersion of
DM particles $\sigma_{\rm DM}(R)$ as $\bar{v}_{\rm
DM}(R)=\sqrt{8/(3\pi)}\,\sigma_{\rm DM}(R)$. The mean velocity and velocity
dispersion of DM particles in the inner part of the Galaxy are not well
constrained by observations and must be estimated from simulations (see $e.g.$
Ref. Board  _et al._ (2021)). Under the assumption of an isothermal sphere,
the DM velocity dispersion can be related to the circular velocity curve by
$\sigma_{\rm DM}(r)=\sqrt{3/2}v_{c}(r)$. However, for simplicity we assume a
constant value for the velocity dispersion. Specifically, we adopt
$\sigma_{\rm DM}=100\,\rm km/s$, which is compatible with the circular
velocities at a radius of $0.1\,{\rm kpc}-1\,\rm kpc$, obtained by combining
the contribution of the stellar bulge, the stellar disk and the DM halo (see
e.g. Refs. Benito  _et al._ (2021); Põder  _et al._ (2023)). This assumption
can be refined in future investigations.

The total DM heat power $\Gamma_{{\rm heat},i}^{\rm DM}$ in the exoplanet is
given (in natural units) by Leane and Smirnov (2021)

$\Gamma_{{\rm heat},i}^{\rm DM}=f\pi R_{E,i}^{2}\rho_{\rm
DM}(R_{i})\bar{v}_{\rm
DM}\left(1+\frac{3}{2}\frac{v_{esc}^{2}(R_{E,i},M_{i})}{\sigma_{\rm
DM}^{2}}\right),$ (3)

where $f$ is the fraction of DM particles captured (and will depend on the
particle physics model), $R_{E,i}$ is the radius of the exoplanet, $v_{{\rm
esc},i}$ is the escape velocity at its surface. We neglect the relative motion
of the exoplanet, which is only a minor correction, and assume that the object
is at rest with respect to the DM halo. For the radius $R_{E,i}$, we adopt
$R_{\rm jup}$ for all exoplanets since the radius of exoplanets older than 500
Myr is roughly equal to that of Jupiter according to current evolutionary
models (see $e.g.$ Saumon and Marley (2008); Phillips  _et al._ (2020); Marley
_et al._ (2021)) and leave for future work the relaxation of this assumption.

Under the assumption that there is no nearby host star ($i.e.$ we consider
free-floating exoplanets or brown dwarfs), the total heat or luminosity of the
exoplanet is then the sum of the contributions from DM heating and internal
heat, given by

$\Gamma_{{\rm heat},i}^{\rm int}+\Gamma_{{\rm heat},i}^{\rm
DM}=4\pi{R_{E,i}}^{2}\sigma_{\rm SB}T_{i}^{4}\epsilon_{i},$ (4)

where $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant, and $\epsilon_{i}$
is the dimensionless emissivity of the given exoplanet. The emissivity is a
measure of heat retention of the object; $\epsilon_{i}=1$ corresponds to a
blackbody, while values $\epsilon_{i}<1$ are a greybody and can produce
spectra peaked in more favorable wavelength regions for detection (see Ref.
Leane and Smirnov (2021) for an extended discussion of this point). As such,
we take the simplifying yet conservative choice of a blackbody, i.e.
$\epsilon_{i}=1$. Note that here we have assumed only DM heating from contact
interactions. Dark kinetic heating of exoplanets due to long-range forces can
increase the signal sensitivity substantially Acevedo  _et al._ (2024c), but
we do not consider this here.

Figure 3 shows expected exoplanet temperatures as a function of age, with and
without DM, for a range of planetary masses. We have used an example NFW DM
profile for the DM heating rates; even larger amounts of DM heating occur for
greater inner slopes as per gNFW profiles. We see that DM heating is most
striking when objects are old, or when objects are smaller such that
backgrounds are comparably lower. Note however that this does not take into
account telescope detectability, which can be easier to achieve for hotter
objects which have higher luminosities. We therefore will focus on candidates
with masses greater than 14 $M_{\rm jup}$ (sufficiently hot with DM heating to
detect in the inner Galaxy), but not larger than 55 $M_{\rm jup}$ (backgrounds
become too high relative to DM heating), which is the optimal detection range
Leane and Smirnov (2021).

Figure 4 depicts the expected excess in exoplanet effective temperature, as a
function of age and mass for exoplanets across three Galactocentric distance
ranges: $[0.1,0.4)\,\rm kpc$, $[0.4,0.7)\,\rm kpc$ and $[0.7,1]\,\rm kpc$. The
decrease in exoplanet abundance moving away from the Galactic center is
attributed to the postulation of an exponential radial number density profile
for the exoplanet population (see Subsections I.3 and I.5 for details), and
the scatter at a given age and mass is due to the distance mixing.

Figure 3: Theoretical effective temperature for the evolution of Jupiters,
Super-Jupiters, and Brown Dwarfs, as per Ref. Phillips  _et al._ (2020), with
(solid) and without (dashed) DM heating. For the DM-heating lines, a
Galactocentric distance of 0.1 kpc and a NFW DM profile with scale radius
$r_{s}=20\,\text{kpc}$, and DM velocity dispersion $\sigma_{\rm DM}=100\text{
km/s}$, are assumed for illustration. Larger DM heating rates occur for
greater inner slope values.  Figure 4: Theoretical excess in the exoplanet
effective temperature due to DM-heating for exoplanets located between
Galactocentric distances of 0.1–0.4 kpc (left panel), 0.4–0.7 kpc (middle
panel) and 0.7–1 kpc (right panel), as a function of exoplanet age and mass.
We consider a sample of 1000 exoplanets under the assumption of a power-law DM
density profile characterized by a normalization of $C=10\,\rm GeV/cm^{3}$ and
an inner slope of $\alpha=1$.

###  I.3 Population Distribution of Exoplanets

For the population of exoplanets we assume independent distributions of
$R_{i}$, $A_{i}$ and $M_{i}$ controlled by global parameters as explained
below.

  * •

Exoplanet Age: We consider a uniform ($i.e.$ constant birthrate) age
distribution from 8–10 Gyr,

$A_{i}\sim\mathcal{U}(8\,\text{Gyr};10\,\text{Gyr}),$ (5)

compatible with observations of an old population of Galactic bulge stars (see
Ref. Barbuy  _et al._ (2018); Joyce  _et al._ (2023) and references therein).

  * •

Exoplanet Mass: We use a Pareto distribution, namely

$M_{i}\sim\frac{\gamma_{M}-1}{M_{\rm min}}\left(\frac{M_{i}}{M_{\rm
min}}\right)^{-\gamma_{M}},$ (6)

which describes a power-law initial mass function. As noted above, we set
$M_{\rm min}=14\,\rm M_{jup}$ and truncate the distribution at $M_{\rm
max}=55\,\rm M_{jup}$ to restrict the simulation within the optimal DM
detection range Leane and Smirnov (2021).

  * •

Exoplanet Spatial Distribution: We assume that exoplanets follow the
distribution of bulge stars, hence we employ an exponential distribution,
given by

$R_{i}\sim\gamma_{R}e^{-\gamma_{R}R_{i}}.$ (7)

Note that for $\gamma_{R}=1.46\,\rm kpc^{-1}$ this is equivalent to the
exoplanets being spatially distributed according to the E2 bulge profile from
Ref. Stanek  _et al._ (1997) (see App. A for details). Since the closer to the
centre of the Galaxy, the higher the expected temperature excess, we limit the
simulation to $R_{i}\in[0.1,1]$, with the lower boundary representing JWST’s
approximate minimum sensitivity threshold for exoplanet temperature
measurement Leane and Smirnov (2021).

###  I.4 Priors

The priors on the population variables are chosen to encompass observational
constraints. Specifically, the population parameter $\gamma_{M}$ is inferred
to be $0.6\pm 0.1$ Kirkpatrick  _et al._ (2021) based on a sample of
exoplanets located within 20 pc from the Sun. This value is consistent with
other estimates in the literature, e.g. Mužić  _et al._ (2017). Conversely,
$\gamma_{R}$ is assigned a uniform distribution that accounts for the
uncertain shape of the distribution of stars and sub-stellar objects in the
central region of our Galaxy. Mathematically, this is represented as:

$\displaystyle\begin{split}\gamma_{M}&\sim\mathcal{N}(0.6,0.1^{2})\textrm{
and}\\\ \gamma_{R}&\sim\mathcal{U}(1,2).\end{split}$ (8)

Finally, we adopt uniform priors for the DM halo parameters within the ranges
$0<\alpha<3$ and $0.5<C[{\rm GeV/cm^{3}}]<40$.

###  I.5 Simulator Configuration

In the absence of real data and in order to verify the inference procedure, we
analyze mock data. Herein, we justify the selected configuration of global
parameters to generate the simulations and outline potential improvements for
future investigations.

  * •

Regarding exoplanet ages, we opted for an old population with uniform ages
spanning from 8 to 10 Gyr. Given the ongoing debate surrounding the age
distribution of stars in the Galactic bulge (see Refs. Barbuy  _et al._
(2018); Joyce  _et al._ (2023) and references therein), we defer accounting
for intermediate-aged stars in the Galactic bulge ($e.g.$, Bensby  _et al._
(2017)) and possible contamination from disc stars to future work.

  * •

Exoplanets are set to have a power-law initial mass function of
$\gamma_{M}=0.6$, consistent with findings in the literature (e.g. Kirkpatrick
_et al._ (2021); Mužić  _et al._ (2017)).

  * •

Furthermore, exoplanets are spatially distributed according to the E2 bulge
profile from Ref. Stanek  _et al._ (1997) (see App. A for details), setting
$\gamma_{R}=1.46\,\rm kpc^{-1}$.

  * •

For global DM parameters, we simulate datasets from true underlying values of
$\alpha$ drawn from its prior range and two settings for the normalization
constant, specifically $C=6\,\rm GeV/cm^{3}$ and $C=20\,\rm GeV/cm^{3}$. The
former value is adopted to reproduce the normalization of the DM density as in
the NFW profile with scale radius of $20~{}\rm kpc$ and local DM density of
$\rho_{0}=0.42\,\rm GeV/cm^{3}$ Karukes  _et al._ (2019); Benito  _et al._
(2021), while the latter value could account for an increase of the DM density
due to, for example, adiabatic contraction Blumenthal  _et al._ (1986).

Furthermore, we adopt relative measurement uncertainties, denoted by $\sigma$,
across a wide range to investigate the impact of noise levels. In particular,
$\sigma$ is fixed to values within $[1,20]\,\%$, where the optimistic (but
unrealistic in the near future) value of 1% allows verification of the
inference procedure, while a more realistic albeit still optimistic scenario
is 20% uncertainty Zhu  _et al._ (2016). Additionally, we vary the number of
observed exoplanets $N$ from 10 to 200. In this first study, we do not
consider selection effects and focus on the construction and validation of an
idealised simulator. Modelling selection effects in a Bayesian hierarchical
model framework is computationally demanding, as it usually involves a
selection probability correction factor that depends in a non-factorizable way
on the parameters of interest, and can only be computed analytically if a
number of assumptions about selection probability are met Kelly (2007). The
next step will be to investigate selection effects with, for example,
Simulation-Based Inference (SBI) methods, which offer more flexibility and
ease of implementation in modelling selection processes.

###  I.6 Inference

To sample the high-dimensional posterior over all 3N latent and 4 global
parameters, we employ the No-U-turn sampler (NUTS) Hoffman and Gelman (2011),
a self-tuning variant of Hamiltonian Monte Carlo (HMC) Neal (2011), as
implemented in the Pyro Bingham  _et al._ (2018) probabilistic programming
language. After 1400 adaptation and warmup steps, we run a chain of 3200
samples (with near-unity autocorrelation length due to the efficiency of
NUTS), which are enough to adequately explore the posterior distribution.
Inference took around 4 minutes (25 minutes) per data set with 10 (100)
objects on a single CPU.

##  II Results

###  II.1 Global Parameter Inference

Figure 5: Posterior marginalized 2D and 1D distributions for the exoplanet
population parameters and DM halo density profile from analyses of mock data
(true parameter values indicated by crosses / vertical lines in the 2D / 1D
plots) with different noise levels and number of detected objects as indicated
in the legends. The 2D contours depict the 2$\sigma$ highest posterior density
(HPD) regions (i.e. with 86% credibility).  Figure 6: Marginal central 68%
credible intervals for the DM density slope $\alpha$, inferred from mock data
with true simulated $\alpha$ values drawn from the prior distribution and
plotted on the abscissa. The legends indicate the number of exoplanets
analysed. As in Fig. 5, in the left and right panels the noise levels are
fixed to $1\%$ and $10\%$, respectively. In the simulation, $C$, $\gamma_{R}$
and $\gamma_{M}$ were fixed to $20\,\rm GeV/cm^{3}$, $1.43$ and $0.6$
respectively and then inferred jointly with $\alpha$ and the latent variables
(exactly as in Fig. 5).

Figure 5 shows example joint posteriors for the global (DM halo and exoplanet
population) parameters obtained from mock data sets comprising 10 or 100
exoplanets, with noise levels set to 1% or 10%. This figure illustrates the
increase in precision obtained with noise level reduced to 1% or, similarly,
with sample size increased to 100 objects. Regarding the population
parameters, the inference is dominated by the initial prior on the slope of
the initial mass function (IMF) power law denoted as $\gamma_{M}$, since it is
small compared to the expected constraints of 100 objects or less. On the
contrary, for the exponential slope of the number density of exoplanets,
$\gamma_{R}$, the prior distribution is broad compared to the expected
constraints of 100 or fewer objects (even with zero uncertainty)111The trivial
estimator $\hat{\gamma}_{R}=\sum\hat{R}_{i}/N$ has variance
$1/(N\gamma_{R}^{2})$ (if the observational noise is negligible), which
amounts to about 10% of the prior range for $\gamma_{R}$ and is an indicator
for the optimal precision in inferring it. However, the fact that we are
restricting $R\in[0.1,1]\,\rm kpc$ further broadens the $\gamma_{R}$, and so
the $\pm 0.1$ precision is not achieved even for 1% observational
uncertainty..

Figure 6 focuses on the marginal results (under the configurations of $N$ and
$\sigma$ indicated in the legend) for the inner DM profile slope $\alpha$.
Here we plot central 68% credible intervals, defined as the 16th and 84th
percentiles of the 1D marginal posterior, versus the true222Here and
elsewhere, we use the word ‘true’ to indicate the parameter values adopted in
generating the data. $\alpha$ values used to simulate the analysed data sets.
From this figure, the increase in precision with increasing number of observed
exoplanets and decreasing noise level is also observed. Additionally, we see
that larger values of the slope can be recovered with greater precision. This
is expected, since for a fixed normalisation $C$, the DM density increases
with increasing $\alpha$ value.

Figure 7 displays the normalized bias of inferring $\alpha$ and $C$ as
measured by the offset of the posterior mean from the true value divided by
the standard deviation of the marginal posterior, i.e.

$\textrm{Normalized bias}=\frac{\langle X\rangle-X_{\rm true}}{\sqrt{\langle
X^{2}\rangle-\langle X\rangle^{2}}},$ (9)

where $X\in\\{\alpha,C\\}$, and inferred from mock data with true parameter
values drawn from their respective priors. We observe a systematic shift of
the posteriors towards the middle of the prior range when the true values of
$\alpha$ are at the extremes of this range, particularly for larger noise
level and smaller number of exoplanets. This is a reflection of a statistical
effect called ‘shrinkage’, which in the context of the Bayesian hierarchical
model employed here ‘shrinks’ posterior estimates towards the prior mean when
the data are less constraining. These figures also reveal a subtle pattern in
which the mean estimate of $\alpha$ tends to be biased low, while that of $C$
tends to be biased high, explained by the anticorrelation between these
parameters (see Fig. 5). This is, in turn, driven by the intrinsic scatter in
temperatures at a given Galactocentric distance due to the scatter in
exoplanet masses and ages which produces an irreducible variance in the
distribution of their values.

Figure 7: Normalized bias for the DM density profile parameters $\alpha$ (top)
and $C$ (bottom) as inferred for a given $(N,\sigma)$ configuration using mock
data with true values drawn from their priors. The left panels display the
true parameter value on the abscissa, while the right panels depict histograms
of the bias across the different $(N,\sigma)$ configurations.  Figure 8: 1D
marginalized posterior for the inner slope of a gNFW density profile, derived
from the analyses of the mock data in Figure 5. The colour scheme, denoting
varying noise levels and number of detected objects, matches that of the
reference figure. The grey, dashed line depicts the prior, while the vertical
black line is the ground truth.

Figure 8 illustrates how the inferred posterior distributions, depicted in
Fig. 5, translate into the marginal posterior for the inner slope $\gamma$ of
a gNFW profile (via Eq. (13)). The transformation into a gNFW assumes equal
inner steepness in both profiles ($\alpha=\gamma$), a local DM density of
$\rho_{0}=0.43\pm 0.02\,\rm GeV/cm^{3}$ Karukes  _et al._ (2019); Benito  _et
al._ (2021), and equivalence in density at 1 kpc from the Galactic center. As
expected, with increasing numbers of objects and smaller observational error
the posterior distribution for $\gamma$ becomes increasingly concentrated
around the true value. As for the scale radius $r_{s}$ of the gNFW profile,
our choice of priors translates into an informative prior on $r_{s}$, which
cannot anyway be meaningfully constrained using only data in the inner kpc and
a single value for the DM density at the solar location. As a consequence of
the informative prior and the weak constraining power of the data, $r_{s}$
cannot be strongly constrained in this setup.

Figure 9: Inferred (yellow, with $1\sigma$ posterior ranges) and true (green)
temperature excesses as a function of inferred (with $1\sigma$ posterior
ranges) and true Galactocentric distance, respectively, for mock observations
of $N=100$ exoplanets generated assuming different DM halo shapes and a
relative uncertainty in exoplanet parameters of $1\%$. The bottom panels show
the normalized bias. Note that errorbars in the inferred Galactocentric
distance are tiny due to the low noise level of 1% and are consequently not
discernible.

###  II.2 Latent Parameter Inference

Our Bayesian hierarchical model allows us to further infer the age, mass, and
Galactocentric distance of each exoplanet in a given simulated data set, along
with the population and DM halo parameters. We leave for App. C a more
thorough examination of the marginal posteriors of the latent parameters.

Figure 9 compares the true and inferred temperature excesses for simulated
datasets with $100$ exoplanets, noise level $\sigma=1\%$ and different shapes
of the DM halo as indicated in the text insets. The normalized residuals,
given by Eq. (9), for each DM halo configuration is shown in the bottom panels
of this figure. The patterns observed in the residuals in the bottom panels of
Fig. 9 can be explained in the context of hierarchical Bayesian modelling by
the trends in the inference of global parameters, as described in the previous
section, namely the shrinkage of posterior estimates and the 2D degeneracy
between $\alpha$ and $C$.

###  II.3 Detection

Our goal is to detect annihilating DM particles through their impact on the
temperature rise of a _population_ of exoplanets. A global indication of this
signal would be an increase in the temperature of exoplanets with decreasing
Galactocentric distance, as the expected excess in effective temperature
scales with the DM density in the exoplanet’s surrounding. Note however that
there is a limitation to the lowest exoplanet temperature that can be measured
at large distances. In Ref. Leane and Smirnov (2021), it was estimated that
JWST could observe exoplanets with Jupiter’s radius down to about 650 K.
However, a number of simplifying assumptions were made in obtaining this
number, including observation time, and that the thermal emission is a
blackbody; considering more accurate spectra can produce higher signals in
telescope bins which are more detectable Leane and Smirnov (2021).
Furthermore, other telescopes such as Roman may be more optimal for this
search. We therefore do not implement this temperature cut in this work, and
while including a cut can reduce signal sensitivity, we defer to future work
for a more detailed analysis of this issue.

In the hierarchical Bayesian framework we adopt throughout, the “detection
significance” can be quantified by the Bayes factor, i.e. the ratio of
Bayesian evidences333Adopting equal _model_ priors (i.e. prior beliefs in the
existence/non-existence of DM) makes the Bayes factor equal to the ratio of
the two models’ posterior probabilities.:

$B=\frac{p(\bm{d}|\text{{DM}})}{p(\bm{d}|\text{{no-DM}})}.$ (10)

Here no-DM indicates a model that only considers the standard exoplanet
evolution, while DM includes the contribution of DM heating. It should be
noted that no-DM is _nested_ within DM: it corresponds to setting $C=0$. In
turn, the Bayesian evidence for each model is the average likelihood of all
its global and latent parameters, distributed according to the hierarchical
prior. Since the number of parameters scales with $N$, the high-dimensional
integral this requires prohibits an exact computation even for as few as 10
exoplanets.

Figure 10: Data realizations (dots) are colour-coded according to their
approximate Bayes factor, equal to the posterior odds ratio between the DM and
no-DM models, calculated for simulated observations from the DM model with
varying numbers of exoplanets and noise levels, $\alpha$ drawn from the prior,
and fiducial $\gamma_{R}=1.43$ and $\gamma_{M}=0.6$. A log Bayes factor larger
than 0 indicates (correctly) a preference for the DM model. The left (right)
panel assumes a lower (larger) DM density. The background shading illustrates
the median trend of Bayes factors, with the black line indicating the locus of
equal posterior odds.

Instead, we employ a simple approximation scheme for the required evidences,
which makes use of the Hamiltonian Monte Carlo posterior inference described
above. For a given data set, we first run the Hamiltonian Monte Carlo with
each model (labelled $M$). From each set of posterior samples, we then derive
a Gaussian posterior _density_ estimate centred on the maximum-a-posteriori
(MAP) value, $\theta_{0}$, and with a covariance matrix matching the empirical
covariance of the samples. This allows us to approximate
$p(\theta_{0}|\mathbf{d},M)$; dividing it by the joint probability at
$\theta_{0}$, given model $M$,
$p(\theta_{0},\mathbf{d}|M)=p(\mathbf{d}|\theta_{0},M)p(\theta_{0}|M)$ (which
is already calculated and recorded in the Hamiltonian Monte Carlo run), gives
(an approximation for) the evidence $p(\mathbf{d}|M)$ by virtue of Bayes’
theorem. This approximation allows us to explore via analyses of mock data the
expected “detectability” of DM under different data set configurations: number
of observed objects, $N$, and noise level, $\sigma$, and for different halo
parameters: slope, $\alpha$, and normalisation, $C$.

Figure 10 illustrates the effect of these four parameters. Both panels of the
figure show the mean of the Bayes factors within 2D bins of ($N$, $\sigma$),
calculated from data set realizations with different $N\in[5,200]$ and
$\sigma\in[1/\%,20/\%]$, drawn uniformly from their range. Furthermore, the
true simulated values of $\alpha$ are drawn from its prior distribution, while
$\gamma_{R}$ and $\gamma_{M}$ are fixed to $1.43$ and $0.6$, respectively.
Finally, the two panels present the results with $C$ fixed at $6\,\rm
GeV/cm^{3}$ (left) and $20\,\rm GeV/cm^{3}$ (right). As expected, the figure
shows that the feasibility of detecting DM-overheated exoplanets – $i.e.$,
values of the log Bayes factor larger than 0 – increases with increasing
number of observed objects and DM density, and decreasing uncertainty, $i.e.$,
towards the top left of each panel. It is also shown that for small $N$ and
large noise level, there is a large scatter of Bayes factors at fixed
$(\alpha,C)$ that can be explained by whether we observe the less massive and
older exoplanets close to the Galactic centre. For example, for a sample
dataset of about 10 exoplanets, there is the potential to detect the signal
when the density slope $\alpha$ exceeds about 1, along with $C=6\,\rm
GeV/cm^{2}$, a normalization constant in agreement with rotation curve
constraints.

##  III Conclusions and Outlook

We have demonstrated the application of Bayesian hierarchical modelling of
overheated exoplanets with the goal of detecting and characterising dark
matter and applied it to simulated data sets representing realisitic exoplanet
populations within the inner kpc of the Galaxy. Considering different
configurations of sample size and uncertainty level, our study has
demonstrated that it is possible to achieve high accuracy in inferring the
shape of the DM halo, both in the inner slope and in the normalisation. In the
case of 1% measurement uncertainties, observations of a sample of about 10
exoplanets can be sufficient for this task, depending on the underlying DM
density at a Galactocentric distance of 1 kpc and the inner DM slope. Despite
the potential variability, the detection of DM-overheated exoplanets within 1
kpc of our Galaxy remains feasible for data sets containing about 100
exoplanets and a DM density profile in agreement with extrapolations of the
inference of the DM density profile in our Galaxy via the stellar disc
rotation curve. However, more realistic but still very optimistic measurement
uncertainties of about 10% require around 100 exoplanets for DM parameter
inference. To optimise detection strategies, follow-up observations with
multiple telescopes should prioritise accurate measurements of the mass, age
and Galactocentric distance of these exoplanets

While our model correctly infers both global and per-object latent parameters,
an application to real data requires the inclusion of observational selection
effects. A viable solution for this is simulation-based inference (SBI), which
offers a straightforward approach to handling arbitrarily complex
probabilistic models both for parameter inference and model selection.
Therefore, integrating the Bayesian hierarchical model developed here into an
SBI framework represents a logical progression in advancing this novel and
promising approach to indirectly search for DM particles through overheated
exoplanets.

## Acknowledgments

We thank Bruce Macintosh, Rain Kipper, Elmo Tempel and Peeter Tenjes for
helpful discussions and comments. RKL was supported in part by the U.S.
Department of Energy under Contract DE-AC02-76SF00515. This work was supported
by the Estonian Research Council grants PRG1006, PRG803 and PSG938, and
partially supported by the ETAg CoE grant “Foundations of the Universe"
(TK202) and University of Tartu ASTRA project 2014-2020.4.01.16-0029 KOMEET.
RT acknowledges co-funding from Next Generation EU, in the context of the
National Recovery and Resilience Plan, Investment PE1 – Project FAIR “Future
Artificial Intelligence Research”. This resource was co-financed by the Next
Generation EU [DM 1555 del 11.10.22]. RT is partially supported by the
Fondazione ICSC, Spoke 3 “Astrophysics and Cosmos Observations”, Piano
Nazionale di Ripresa e Resilienza Project ID CN00000013 “Italian Research
Center on High-Performance Computing, Big Data and Quantum Computing” funded
by MUR Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di
ricerca e creazione di “campioni nazionali di R&S (M4C2-19 )” - Next
Generation EU (NGEU). This article/publication is based upon work from COST
Action COSMIC WISPers CA21106, supported by COST (European Cooperation in
Science and Technology)

##  Appendix A Spatial Distribution

The spatial distribution of exoplanets is assumed to follow the E2 bulge
profile from Stanek  _et al._ (1997). Thus, the number density profile, given
in arbitrary units, takes the form444We adopt Galactocentric spherical
coordinates with $\phi$ increasing in the direction of Galactic rotation, i.e.
clock-wise as seen from the North Galactic Pole (NGP), and where $\theta$ is
measured from the NGP.

$n(r,\phi,\theta)=n_{0}\exp{\left[-\gamma_{R}(\theta,\phi)r\right]}$ (11)

with

$\gamma_{R}(\theta,\phi)=\sqrt{\sin^{2}{\theta}\left[\left(\frac{\cos{(\phi+\alpha)}}{x_{0}}\right)^{2}+\left(\frac{\sin{(\phi+\alpha)}}{y_{0}}\right)^{2}\right]+\left(\frac{\cos{\theta}}{z_{0}}\right)^{2}},$

(12)

$x_{0}=$0.899\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$,
$y_{0}=$0.386\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$,
$z_{0}=$0.250\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$, and $\alpha=23.8^{\rm
o}$ is the angle between the bulges’s major axis and the Sun-Galactic center
line. The scale radius and scale length of the disc are
$2.15\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$ and
$0.40\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$, respectively555Scale parameters
in the planed defined by $\theta=\pi/2$ (i.e. $xy$-plane), which are taken
from observational studies, are re-scaled as $R_{0}$..

We are only interested in the distances of exoplanets to the Galactic center,
rather than in their three-dimensional position. Is it for this reason that we
get rid of $\theta$ and $\phi$ by integrating $n(r,\phi,\theta)$ over these
coordinates. The integration must be done over the region covered by the
observations. However, for simplicity, we have assumed $\phi=0$ and
$\theta=\pi/2$, which gives a value of $\gamma_{R}=1.46\,\rm kpc^{-1}$.

##  Appendix B Power-law Bias

As discussed in the main text for radii that are significantly smaller than
the scale radius of the profile $r_{s}$, we can expand the gNFW equation

$\rho_{\mathrm{DM}}(R;\gamma,r_{s},\rho_{0})=\rho_{0}\left(\frac{R_{0}}{R}\right)^{\gamma}\left(\frac{r_{s}+R_{0}}{r_{s}+R}\right)^{3-\gamma},$
(13)

which yields a power law expression

$\rho_{\rm DM}(R;C,\gamma)\approx C\,R^{-\gamma}\quad\text{for }R\ll r_{s}.$
(14)

As discussed, this helps to simplify the fitting procedure and breaks the
degeneracy between the scale radius and the inner steepness parameter
$\gamma$. However, this procedure introduces a systematic bias to slightly
larger $\gamma$ values. A simple shape analysis shows that comparing the shape
of the full expression in Eq. 13 and the expanded expression in Eq. 14 leads
to very similar graphs for slightly different scale parameters $\gamma$,
modulo a total normalization correction.

Table 2 shows several values for the $\gamma$ parameters that lead to
identical inner slopes for the DM profile. We can make two observations. One,
is that the bias when using the power law leads to systematically larger
values for $\gamma$. The second is, that the bias can be somewhat substantial
for $\gamma$ values smaller than one, and vanishes entirely for steeper
profiles, when $\gamma$ approaches two. In the intermediate regime, that
corresponds to the standard NFW profile, the bias is about $10\%$. We conclude
that the induced bias is well below the experimental accuracy that we can hope
to archive in the near and intermediate future, and thus the expansion
prescription is a valid method for the $\gamma$ factor reconstruction.

$\gamma$ from full gNFW |  $\gamma$ from Power Law | Induced bias
---|---|---
0.5 | 0.57 | $+14\%$
1.0 | 1.1 | $+10\%$
1.2 | 1.3 | $+8\%$
1.5 | 1.6 | $+6\%$
1.8 | 1.8 | $0\%$
Table 2: The induced bias for several benchmark values of the DM profile
parameter $\gamma$.

##  Appendix C Closer Examination of Latent Parameter Inference

For a given mock dataset, the accuracy in the inferred age, mass and
Galactocentric distance, as measured by the normalized bias, is plotted in the
top, middle and bottom panels, respectively, of Figure 11, as a function of
the true values of these latent parameters. The true latent parameter space
has been divided into three equal-width intervals, with color bands within
each interval representing the bias region encompassing 1$\sigma$ of the
points is shown to ease visualization. These bands highlight the effect of
‘shrinkage’ discussed in the text, particularly notable at high noise levels,
on latent parameter inference (most notably for $A$), demonstrating a tendency
for the marginal posteriors to be shrunk, in virtue of the Bayesian
hierarchical model, towards the mean of the prior range.

Figure 11: Normalized bias for the latent ages (top), masses (middle) and
Galactocentric distance (bottom) inferred, along with the global parameters,
from a mock observation with 100 exoplanets and $\sigma$ as indicated in the
legend. In the simulation, $\alpha$, $C$, $\gamma_{R}$ and $\gamma_{M}$ were
fixed to $1.5$, $20\,\rm GeV/cm^{3}$, $1.43$ and $0.6$, respectively. The
shaded bands indicate the standard deviation computed within the top-hat bin
shown.

## References

  * Mužić  _et al._ (2017) K. Mužić, R. Schödel, A. Scholz, V. C. Geers, R. Jayawardhana, J. Ascenso, and L. A. Cieza, MNRAS 471, 3699 (2017), arXiv:1707.00277 [astro-ph.SR] .
  * Leane and Smirnov (2021) R. K. Leane and J. Smirnov, Phys. Rev. Lett. 126, 161101 (2021), arXiv:2010.00015 [hep-ph] .
  * Goldman and Nussinov (1989) I. Goldman and S. Nussinov, Phys. Rev. D40, 3221 (1989).
  * Gould  _et al._ (1990) A. Gould, B. T. Draine, R. W. Romani, and S. Nussinov, Phys. Lett. B238, 337 (1990).
  * Kouvaris (2008) C. Kouvaris, Phys. Rev. D77, 023006 (2008), arXiv:0708.2362 [astro-ph] .
  * Bertone and Fairbairn (2008) G. Bertone and M. Fairbairn, Phys. Rev. D77, 043515 (2008), arXiv:0709.1485 [astro-ph] .
  * Spolyar  _et al._ (2008) D. Spolyar, K. Freese, and P. Gondolo, Phys. Rev. Lett. 100, 051101 (2008), arXiv:0705.0521 [astro-ph] .
  * Iocco (2008) F. Iocco, Astrophys. J. Lett. 677, L1 (2008), arXiv:0802.0941 [astro-ph] .
  * Freese  _et al._ (2009) K. Freese, P. Gondolo, J. Sellwood, and D. Spolyar, Astrophys. J. 693, 1563 (2009), arXiv:0805.3540 [astro-ph] .
  * Taoso  _et al._ (2008) M. Taoso, G. Bertone, G. Meynet, and S. Ekstrom, Phys. Rev. D 78, 123510 (2008), arXiv:0806.2681 [astro-ph] .
  * Sivertsson and Gondolo (2011) S. Sivertsson and P. Gondolo, Astrophys. J. 729, 51 (2011), arXiv:1006.0025 [astro-ph.CO] .
  * Freese  _et al._ (2016) K. Freese, T. Rindler-Daller, D. Spolyar, and M. Valluri, Rept. Prog. Phys. 79, 066902 (2016), arXiv:1501.02394 [astro-ph.CO] .
  * Salati and Silk (1989) P. Salati and J. Silk, Astrophys. J. 338, 24 (1989).
  * Fairbairn  _et al._ (2008) M. Fairbairn, P. Scott, and J. Edsjo, Phys. Rev. D 77, 047301 (2008), arXiv:0710.3396 [astro-ph] .
  * Scott  _et al._ (2009) P. Scott, M. Fairbairn, and J. Edsjo, Mon. Not. Roy. Astron. Soc. 394, 82 (2009), arXiv:0809.1871 [astro-ph] .
  * Lopes and Lopes (2021) J. Lopes and I. Lopes, Astron. Astrophys. 651, A101 (2021), arXiv:2107.13885 [astro-ph.SR] .
  * de Lavallaz and Fairbairn (2010) A. de Lavallaz and M. Fairbairn, Phys. Rev. D81, 123521 (2010), arXiv:1004.0629 [astro-ph.GA] .
  * Kouvaris and Tinyakov (2010) C. Kouvaris and P. Tinyakov, Phys. Rev. D82, 063531 (2010), arXiv:1004.0586 [astro-ph.GA] .
  * McDermott  _et al._ (2012) S. D. McDermott, H.-B. Yu, and K. M. Zurek, Phys. Rev. D85, 023519 (2012), arXiv:1103.5472 [hep-ph] .
  * Kouvaris and Tinyakov (2011a) C. Kouvaris and P. Tinyakov, Phys. Rev. Lett. 107, 091301 (2011a), arXiv:1104.0382 [astro-ph.CO] .
  * Guver  _et al._ (2014) T. Guver, A. E. Erkoca, M. Hall Reno, and I. Sarcevic, JCAP 1405, 013 (2014), arXiv:1201.2400 [hep-ph] .
  * Bramante  _et al._ (2013) J. Bramante, K. Fukushima, and J. Kumar, Phys. Rev. D87, 055012 (2013), arXiv:1301.0036 [hep-ph] .
  * Bell  _et al._ (2013) N. F. Bell, A. Melatos, and K. Petraki, Phys. Rev. D87, 123507 (2013), arXiv:1301.6811 [hep-ph] .
  * Bramante  _et al._ (2014) J. Bramante, K. Fukushima, J. Kumar, and E. Stopnitzky, Phys. Rev. D89, 015010 (2014), arXiv:1310.3509 [hep-ph] .
  * Bertoni  _et al._ (2013) B. Bertoni, A. E. Nelson, and S. Reddy, Phys. Rev. D88, 123505 (2013), arXiv:1309.1721 [hep-ph] .
  * Kouvaris and Tinyakov (2011b) C. Kouvaris and P. Tinyakov, Phys. Rev. D83, 083512 (2011b), arXiv:1012.2039 [astro-ph.HE] .
  * McCullough and Fairbairn (2010) M. McCullough and M. Fairbairn, Phys. Rev. D81, 083520 (2010), arXiv:1001.2737 [hep-ph] .
  * Angeles Perez-Garcia and Silk (2015) M. Angeles Perez-Garcia and J. Silk, Phys. Lett. B744, 13 (2015), arXiv:1403.6111 [astro-ph.SR] .
  * Bramante (2015) J. Bramante, Phys. Rev. Lett. 115, 141301 (2015), arXiv:1505.07464 [hep-ph] .
  * Graham  _et al._ (2015) P. W. Graham, S. Rajendran, and J. Varela, Phys. Rev. D92, 063007 (2015), arXiv:1505.04444 [hep-ph] .
  * Cermeno  _et al._ (2016) M. Cermeno, M. Perez-Garcia, and J. Silk, Phys. Rev. D94, 063001 (2016), arXiv:1607.06815 [astro-ph.HE] .
  * Krall and Reece (2018) R. Krall and M. Reece, Chin. Phys. C42, 043105 (2018), arXiv:1705.04843 [hep-ph] .
  * Leane  _et al._ (2017) R. K. Leane, K. C. Y. Ng, and J. F. Beacom, Phys. Rev. D95, 123016 (2017), arXiv:1703.04629 [astro-ph.HE] .
  * McKeen  _et al._ (2018) D. McKeen, A. E. Nelson, S. Reddy, and D. Zhou, Phys. Rev. Lett. 121, 061802 (2018), arXiv:1802.08244 [hep-ph] .
  * Baryakhtar  _et al._ (2017) M. Baryakhtar, J. Bramante, S. W. Li, T. Linden, and N. Raj, Phys. Rev. Lett. 119, 131801 (2017), arXiv:1704.01577 [hep-ph] .
  * Raj  _et al._ (2018) N. Raj, P. Tanedo, and H.-B. Yu, Phys. Rev. D97, 043006 (2018), arXiv:1707.09442 [hep-ph] .
  * Bell  _et al._ (2018) N. F. Bell, G. Busoni, and S. Robles, JCAP 1809, 018 (2018), arXiv:1807.02840 [hep-ph] .
  * Chen and Lin (2018) C.-S. Chen and Y.-H. Lin, JHEP 08, 069 (2018), arXiv:1804.03409 [hep-ph] .
  * Hamaguchi  _et al._ (2019) K. Hamaguchi, N. Nagata, and K. Yanagi, Phys. Lett. B795, 484 (2019), arXiv:1905.02991 [hep-ph] .
  * Camargo  _et al._ (2019) D. A. Camargo, F. S. Queiroz, and R. Sturani, JCAP 1909, 051 (2019), arXiv:1901.05474 [hep-ph] .
  * Bell  _et al._ (2019) N. F. Bell, G. Busoni, and S. Robles, JCAP 1906, 054 (2019), arXiv:1904.09803 [hep-ph] .
  * Garani and Heeck (2019) R. Garani and J. Heeck, Phys. Rev. D100, 035039 (2019), arXiv:1906.10145 [hep-ph] .
  * Acevedo  _et al._ (2020) J. F. Acevedo, J. Bramante, R. K. Leane, and N. Raj, JCAP 03, 038 (2020), arXiv:1911.06334 [hep-ph] .
  * Joglekar  _et al._ (2020a) A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Physics Letters B 809, 135767 (2020a), arXiv:1911.13293 [hep-ph] .
  * Joglekar  _et al._ (2020b) A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Phys. Rev. D 102, 123002 (2020b), arXiv:2004.09539 [hep-ph] .
  * Bell  _et al._ (2020) N. F. Bell, G. Busoni, S. Robles, and M. Virgato, JCAP 2020, 028 (2020), arXiv:2004.14888 [hep-ph] .
  * Garani  _et al._ (2021) R. Garani, A. Gupta, and N. Raj, Phys. Rev. D 103, 043019 (2021), arXiv:2009.10728 [hep-ph] .
  * Leane  _et al._ (2021) R. K. Leane, T. Linden, P. Mukhopadhyay, and N. Toro, (2021), arXiv:2101.12213 [astro-ph.HE] .
  * Acevedo  _et al._ (2024a) J. F. Acevedo, R. K. Leane, and L. Santos-Olmsted, JCAP 03, 042 (2024a), arXiv:2309.10843 [hep-ph] .
  * John  _et al._ (2023) I. John, R. K. Leane, and T. Linden, (2023), arXiv:2311.16228 [astro-ph.HE] .
  * Croon and Sakstein (2023) D. Croon and J. Sakstein, (2023), arXiv:2310.20044 [astro-ph.HE] .
  * Acevedo  _et al._ (2024b) J. F. Acevedo, J. Bramante, Q. Liu, and N. Tyagi, (2024b), arXiv:2404.10039 [hep-ph] .
  * Mack  _et al._ (2007) G. D. Mack, J. F. Beacom, and G. Bertone, Phys. Rev. D 76, 043523 (2007), arXiv:0705.4298 [astro-ph] .
  * Chauhan and Mohanty (2016) B. Chauhan and S. Mohanty, Phys. Rev. D 94, 035024 (2016), arXiv:1603.06350 [hep-ph] .
  * Bramante  _et al._ (2020) J. Bramante, A. Buchanan, A. Goodman, and E. Lodhi, Phys. Rev. D 101, 043001 (2020), arXiv:1909.11683 [hep-ph] .
  * Adler (2009) S. L. Adler, Phys. Lett. B 671, 203 (2009), arXiv:0808.2823 [astro-ph] .
  * Kawasaki  _et al._ (1992) M. Kawasaki, H. Murayama, and T. Yanagida, Prog. Theor. Phys. 87, 685 (1992).
  * Mitra (2004) S. Mitra, Phys. Rev. D 70, 103517 (2004), arXiv:astro-ph/0408341 .
  * Garani and Tinyakov (2020) R. Garani and P. Tinyakov, Phys. Lett. B 804, 135403 (2020), arXiv:1912.00443 [hep-ph] .
  * Chan and Lee (2020) M. H. Chan and C. M. Lee, Phys. Rev. D 102, 023024 (2020), arXiv:2007.01589 [astro-ph.HE] .
  * Leane and Linden (2023) R. K. Leane and T. Linden, Phys. Rev. Lett. 131, 071001 (2023), arXiv:2104.02068 [astro-ph.HE] .
  * Blanco and Leane (2023) C. Blanco and R. K. Leane, (2023), arXiv:2312.06758 [hep-ph] .
  * Croon and Smirnov (2023) D. Croon and J. Smirnov, (2023), arXiv:2309.02495 [hep-ph] .
  * Linden  _et al._ (2024) T. Linden, T. T. Q. Nguyen, and T. M. P. Tait, (2024), arXiv:2402.01839 [hep-ph] .
  * Leane and Tong (2024) R. K. Leane and J. Tong, (2024), arXiv:2405.05312 [hep-ph] .
  * Perryman  _et al._ (2014) M. Perryman, J. Hartman, G. Á. Bakos, and L. Lindegren, Astrophys. J. 797, 14 (2014), arXiv:1411.1173 [astro-ph.EP] .
  * Green  _et al._ (2012) J. Green _et al._ , “Wide-field infrared survey telescope (wfirst) final report,” (2012), arXiv:1208.4012 [astro-ph.IM] .
  * Johnson  _et al._ (2020) S. A. Johnson, M. Penny, B. S. Gaudi, E. Kerins, N. J. Rattenbury, A. C. Robin, S. C. Novati, and C. B. Henderson, The Astronomical Journal 160, 123 (2020).
  * Bachelet and Penny (2019) E. Bachelet and M. Penny, The Astrophysical Journal 880, L32 (2019).
  * Zhu  _et al._ (2016) W. Zhu, S. Calchi Novati, A. Gould, A. Udalski, C. Han, Y. Shvartzvald, C. Ranc, U. G. Jørgensen, R. Poleski, V. Bozza, C. Beichman, G. Bryden, S. Carey, B. S. Gaudi, C. B. Henderson, R. W. Pogge, I. Porritt, B. Wibking, J. C. Yee, SPITZER Team, M. Pawlak, M. K. Szymański, J. Skowron, P. Mróz, S. Kozłowski, Ł. Wyrzykowski, P. Pietrukowicz, G. Pietrzyński, I. Soszyński, K. Ulaczyk, OGLE Group, J. Y. Choi, H. Park, Y. K. Jung, I. G. Shin, M. D. Albrow, B. G. Park, S. L. Kim, C. U. Lee, S. M. Cha, D. J. Kim, Y. Lee, KMTNET Group, M. Friedmann, S. Kaspi, D. Maoz, WISE Group, M. Hundertmark, R. A. Street, Y. Tsapras, D. M. Bramich, A. Cassan, M. Dominik, E. Bachelet, S. Dong, R. Figuera Jaimes, K. Horne, S. Mao, J. Menzies, R. Schmidt, C. Snodgrass, I. A. Steele, J. Wambsganss, RoboNeT Team, J. Skottfelt, M. I. Andersen, M. J. Burgdorf, S. Ciceri, G. D’Ago, D. F. Evans, S. H. Gu, T. C. Hinse, E. Kerins, H. Korhonen, M. Kuffmeier, L. Mancini, N. Peixinho, A. Popovas, M. Rabus, S. Rahvar, R. Tronsgaard, G. Scarpetta, J. Southworth, J. Surdej, C. von Essen, Y. B. Wang, O. Wertz, and MiNDSTEP Group, Astrophys. J. 825, 60 (2016), arXiv:1510.02097 [astro-ph.SR] .
  * Eddington (1913) A. S. Eddington, MNRAS 73, 359 (1913).
  * Kelly (2007) B. C. Kelly, Astrophys. J. 665, 1489 (2007), arXiv:0705.2774 [astro-ph] .
  * March  _et al._ (2011) M. C. March, R. Trotta, P. Berkes, G. D. Starkman, and P. M. Vaudrevange, MNRAS 418, 2308 (2011), arXiv:1102.3237 [astro-ph.CO] .
  * Phillips  _et al._ (2020) M. W. Phillips, P. Tremblin, I. Baraffe, G. Chabrier, N. F. Allard, F. Spiegelman, J. M. Goyal, B. Drummond, and E. Hébrard, Astron. Astrophys. 637, A38 (2020), arXiv:2003.13717 [astro-ph.SR] .
  * Karukes  _et al._ (2019) E. V. Karukes, M. Benito, F. Iocco, R. Trotta, and A. Geringer-Sameth, JCAP 09, 046 (2019), arXiv:1901.02463 [astro-ph.GA] .
  * Benito  _et al._ (2019) M. Benito, A. Cuoco, and F. Iocco, JCAP 2019, 033 (2019), arXiv:1901.02460 [astro-ph.GA] .
  * Benito  _et al._ (2021) M. Benito, F. Iocco, and A. Cuoco, Physics of the Dark Universe 32, 100826 (2021), arXiv:2009.13523 [astro-ph.GA] .
  * Board  _et al._ (2021) E. Board, N. Bozorgnia, L. E. Strigari, R. J. J. Grand, A. Fattahi, C. S. Frenk, F. Marinacci, J. F. Navarro, and K. A. Oman, JCAP 2021, 070 (2021), arXiv:2101.06284 [astro-ph.CO] .
  * Põder  _et al._ (2023) S. Põder, M. Benito, J. Pata, R. Kipper, H. Ramler, G. Hütsi, I. Kolka, and G. F. Thomas, Astron. Astrophys. 676, A134 (2023), arXiv:2309.02895 [astro-ph.GA] .
  * Saumon and Marley (2008) D. Saumon and M. S. Marley, Astrophys. J. 689, 1327 (2008), arXiv:0808.2611 [astro-ph] .
  * Marley  _et al._ (2021) M. S. Marley, D. Saumon, C. Visscher, R. Lupu, R. Freedman, C. Morley, J. J. Fortney, C. Seay, A. J. R. W. Smith, D. J. Teal, and R. Wang, arXiv e-prints , arXiv:2107.07434 (2021), arXiv:2107.07434 [astro-ph.SR] .
  * Acevedo  _et al._ (2024c) J. F. Acevedo, R. K. Leane, and A. J. Reilly, (2024c), arXiv:2405.02393 [astro-ph.EP] .
  * Barbuy  _et al._ (2018) B. Barbuy, C. Chiappini, and O. Gerhard, Annu. Rev. Astron. Astrophys. 56, 223 (2018), arXiv:1805.01142 [astro-ph.GA] .
  * Joyce  _et al._ (2023) M. Joyce, C. I. Johnson, T. Marchetti, R. M. Rich, I. Simion, and J. Bourke, Astrophys. J. 946, 28 (2023), arXiv:2205.07964 [astro-ph.SR] .
  * Stanek  _et al._ (1997) K. Z. Stanek, A. Udalski, M. SzymaŃski, J. KaŁuŻny, Z. M. Kubiak, M. Mateo, and W. KrzemiŃski, Astrophys. J. 477, 163 (1997), arXiv:astro-ph/9605162 [astro-ph] .
  * Kirkpatrick  _et al._ (2021) J. D. Kirkpatrick, C. R. Gelino, J. K. Faherty, A. M. Meisner, D. Caselden, A. C. Schneider, F. Marocco, A. J. Cayago, R. L. Smart, P. R. Eisenhardt, M. J. Kuchner, E. L. Wright, M. C. Cushing, K. N. Allers, D. C. Bardalez Gagliuffi, A. J. Burgasser, J. Gagné, S. E. Logsdon, E. C. Martin, J. G. Ingalls, P. J. Lowrance, E. S. Abrahams, C. Aganze, R. Gerasimov, E. C. Gonzales, C.-C. Hsu, N. Kamraj, R. Kiman, J. Rees, C. Theissen, K. Ammar, N. S. Andersen, P. Beaulieu, G. Colin, C. A. Elachi, S. J. Goodman, L. Gramaize, L. K. Hamlet, J. Hong, A. Jonkeren, M. Khalil, D. W. Martin, W. Pendrill, B. Pumphrey, A. Rothermich, A. Sainio, A. Stenner, C. Tanner, M. Thévenot, N. V. Voloshin, J. Walla, Z. Wędracki, and Backyard Worlds: Planet 9 Collaboration, The Astrophysical Journal Supp. Series 253, 7 (2021), arXiv:2011.11616 [astro-ph.SR] .
  * Bensby  _et al._ (2017) T. Bensby, S. Feltzing, A. Gould, J. C. Yee, J. A. Johnson, M. Asplund, J. Meléndez, S. Lucatello, L. M. Howes, A. McWilliam, A. Udalski, M. K. Szymański, I. Soszyński, R. Poleski, Ł. Wyrzykowski, K. Ulaczyk, S. Kozłowski, P. Pietrukowicz, J. Skowron, P. Mróz, M. Pawlak, F. Abe, Y. Asakura, A. Bhattacharya, I. A. Bond, D. P. Bennett, Y. Hirao, M. Nagakane, N. Koshimoto, T. Sumi, D. Suzuki, and P. J. Tristram, Astron. Astrophys. 605, A89 (2017), arXiv:1702.02971 [astro-ph.GA] .
  * Blumenthal  _et al._ (1986) G. R. Blumenthal, S. M. Faber, R. Flores, and J. R. Primack, Astrophys. J. 301, 27 (1986).
  * Hoffman and Gelman (2011) M. D. Hoffman and A. Gelman, arXiv e-prints , arXiv:1111.4246 (2011), arXiv:1111.4246 [stat.CO] .
  * Neal (2011) R. Neal, in  _Handbook of Markov Chain Monte Carlo_ (2011) pp. 113–162.
  * Bingham  _et al._ (2018) E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. Szerlip, P. Horsfall, and N. D. Goodman, arXiv e-prints , arXiv:1810.09538 (2018), arXiv:1810.09538 [cs.LG] .