# A halo of trapped interstellar matter surrounding the solar system

Jorge Peñarrubia1,2
1Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
2Centre for Statistics, University of Edinburgh, School of Mathematics,
Edinburgh EH9 3FD, UK<EMAIL_ADDRESS>
###### Abstract

This paper shows that gravitating bodies travelling through the Galaxy can
trap lighter interstellar particles that pass nearby with small relative
velocities onto temporarily-bound orbits. The capture mechanism is driven by
the Galactic tidal field, which can decelerate infalling objects to a degree
where their binding energy becomes negative. Over time, trapped particles
build a local overdensity – or ‘halo’– that reaches a steady state as the
number of particles being captured equals that being tidally stripped. This
paper uses classical stochastic techniques to calculate the capture rate and
the phase-space distribution of particles trapped by a point-mass. In a steady
state, bound particles generate a density enhancement that scales as
$\delta(r)\sim r^{-3/2}$ (a.k.a ‘density spike’) and follow a velocity
dispersion profile $\sigma_{h}(r)\sim r^{-1/2}$. Collisionless $N$-body
experiments show excellent agreement with these theoretical predictions within
a distance range $r\gtrsim r_{\epsilon}$, where $r_{\epsilon}\simeq
0.8\,\exp[-V_{\star}^{2}/(2\sigma^{2})]\,Gm_{\star}/\sigma^{2}$ is the thermal
critical radius of a point-mass $m_{\star}$ moving with a speed $V_{\star}$
through a sea of particles with a velocity dispersion $\sigma$. Preliminary
estimates that ignore collisions with planets and Galactic substructures
suggest that the solar system may be surrounded by a halo that contains the
order of $N^{\rm ISO}(<0.1\,{\rm pc})\sim 10^{7}$ energetically-bound
’Oumuamua-like objects, and a dark matter mass of $M^{\rm DM}(<0.1\,{\rm
pc})\sim 10^{-13}M_{\odot}$. The presence of trapped interstellar matter in
the solar system can affect current estimates on the size of the Oort Cloud,
and leave a distinct signal in direct dark matter detection experiments.

###### keywords:

Galaxy: kinematics and dynamics – Cosmology: dark matter– Galaxy: local
interstellar matter–comets: general – minor planets— minor planets, asteroids:
individual: 1I/’Oumuamua – Oort Cloud.

##  1 Introduction

The solar system consists of the Sun and the population of planets and minor
bodies gravitationally bound to it. Although this definition is unambiguous,
it ignores the Milky Way. To date, the dynamical interplay between the Sun and
the interstellar surroundings is poorly understood, specially in the outer-
most regions of the solar system, where external forces become comparable to
the Sun’s gravitational attraction. In this region, the boundary with the
Galaxy is blurred, as there is a continuous exchange of material that goes in
both ways: from the interstellar medium into the solar system via
gravitational capture, and from the solar system out to the Milky Way via
tidal heating and stripping111In this paper, the word ‘capture’ is broadly
used to describe any dynamical process wherein a body undergoes a transition
from Milky Way orbit to a heliocentric one..

Most planet formation models assume that the solar system formed about 4.5
billion years ago in a dense molecular cloud. Planets, moons, asteroids, and
other small bodies formed in a protoplanetary disc that became progressively
eroded due to collisions with the planets and the action of external tidal
forces (e.g. Kokubo & Ida 1998; Portegies Zwart 2021a and references therein).
In this scenario, the material currently orbiting the Sun represents the
surviving remnants of a nearly obliterated disc (e.g. Charnoz & Morbidelli
2003; Johansen et al. 2021, Portegies Zwart 2021; Portegies Zwart et al.
2021).

This picture has been recently challenged by the discovery of two visitors
from the interstellar space, the ’Oumuamua object (Meech et al. 2017) and the
comet Borisov (Jewitt & Luu 2019), which provided undeniable evidence that the
solar system may also host a sizeable population of minor bodies of extra-
solar origin (e.g., Siraj & Loeb 2019; Namouni & Morais 2020). Like ’Oumuamua
and Borisov, the vast majority of interstellar objects are expected to move
across the solar system on unbound trajectories. Although these visitors from
outer space represent transient occurrences, they point to the possibility
that the Sun may have captured an unknown number of these objects in the past.

It is important to bear in mind that gravitational capture cannot happen in
the classical two-body problem. Consider for example an interstellar body
approaching the Sun on a hyperbolic, unbound trajectory with a specific energy
$E_{\infty}=v^{2}_{\infty}/2>0$, where $v_{\infty}$ is the asymptotic relative
velocity in the limit $t\to-\infty$. If there is no internal or external force
acting on this object other than the gravitational attraction of the Sun, then
energy conservation demands $E=v^{2}/2-Gm_{\star}/r=E_{\infty}$. This in turn
means that the infalling speed progressively accelerates as it approaches the
Sun as $v^{2}=v^{2}_{\infty}+2Gm_{\star}/r$, which is typically known as
‘gravitational focusing’, and experiences a mirrored deceleration on its way
out of the solar system. Therefore, gravitational capture is not possible
unless one or several mechanisms ‘slow down’ the infalling body and dissipate
its kinetic energy to a degree wherein the binding energy becomes negative,
$E<0$.

The best studied dissipative mechanism for gravitational capture is 3-body
interactions with planets (e.g., Valtonen & Innanen 1982; Valtonen 1983; Siraj
& Loeb 2019, 2021; Lingam & Loeb 2018; Hands & Dehnen 2020; Napier et al.
2021a,b; Dehnen & Hands 2022; Dehnen et al. 2021). Dark Matter (DM) particles
–if they exist– can also be captured out of the Milky Way DM halo via
collisions with planets. This scenario has been extensively explored for the
case of the solar system (Gould 1987, 1988; Lundberg & Edsjö 2004; Xu & Siegel
2008; Peter 2009a,b,c; Edsjö & Peter 2010; Iorio 2010; Anderson et al. 2020;
Lehmann et al. 2021). More recently, Moro-Martín & Norman (2022) have also
studied the effect of gas drag on the capture of interstellar objects during
the early stages of star and planet formation.

In past studies, capture rates by the solar system are typically calculated
under the assumption of isolation (i.e. neglecting external forces), and
sampling a distribution of encounter speeds that follows either the velocity
distribution of Milky Way stars in the Solar neighbourhood, or that
appropriate for the solar birth environment (e.g. Portegies Zwart 2009; Adams
2010; Pfalzner 2013; Parker 2020).

To date, the role of the Milky Way potential in capturing interstellar matter
remains poorly understood. To address this problem, it becomes necessary to
drop the isolation assumption and account for the motion the Sun and that of
the interstellar visitors around the Milky Way centre. To this aim, two sets
of coupled differential equations must be solved simultaneously. The first one
describes the orbit of the Sun around the Milky Way

$\displaystyle\ddot{\textbf{{R}}_{\star}}=-\nabla\Phi_{g}({\textbf{{R}}_{\star}}).$
(1)

where ${\textbf{{R}}_{\star}}$ is the galactocentric location of the Sun and
$\Phi_{g}$ is the Galactic potential. The second set describes the
trajectories of interstellar objects in the vicinity of the solar system

$\displaystyle\ddot{\textbf{{R}}}=-\nabla\Phi_{g}({\textbf{{R}}})+{\textbf{{F}}}_{\star}+{\textbf{{F}}}_{p}.$
(2)

where
${\textbf{{F}}}_{\star}=-Gm_{\star}\,({\textbf{{R}}}-{\textbf{{R}}}_{\star})//|{\textbf{{R}}}-{\textbf{{R}}}_{\star}|^{3}$
is the specific force induced by the Sun, and
${\textbf{{F}}}_{p}=-\sum_{i=1}^{N_{p}}Gm_{p,i}({\textbf{{R}}}-{\textbf{{R}}}_{i})/|{\textbf{{R}}}-{\textbf{{R}}}_{i}|^{3}$
is the total specific force generated by $N_{p}$ planets (or any other massive
bodies) moving in or around the solar system. For simplicity, the remainder of
this work purposely sets $F_{p}=0$. We come back to this issue in Section 4.
Note also that the ‘collisionless’ approximation corresponds to the case
$F_{\star}=F_{p}=0$, wherein the orbits of interstellar objects are solely
governed by the smooth galactic potential $\Phi_{g}$.

Figure 1: Trajectories of particles trapped around a point mass
$m_{\star}/M_{g}=5.3\times 10^{-4}$ moving on a circular orbit with
galactocentric radius $R=8\,{\rm kpc}$ (dotted-green line), where
$M_{g}=1.84\times 10^{12}M_{\odot}$ is the mass of the host dark matter halo
(go to §3 for details on the numerical experiments). Left and middle columns
show orbits in a reference frame centred at the host Galactic and Keplerian
potentials, respectively. Black/red solid lines show the orbits of trapped
particles with positive/negative orbital energies $E=v^{2}/2-Gm_{\star}/r$.
The separation between the point-mass and the dark matter particle is plotted
in the right column as a function of time in units of the scale-length
$d\equiv|\nabla\rho/\rho|^{-1}$. Trajectories can be broadly divided in three
categories: a) Transient captures have a duration comparable to the local
dynamical time, $T(d)$, given by (15). (b) Semi-stable captures perform
several revolutions in the Keplerian potential before being tidally stripped.
(c) Repeated captures move on intricate trajectories dictated by the combined
Sun+Galactic potentials.

Defining the relative separation vector
${\textbf{{r}}}={\textbf{{R}}}-{\textbf{{R}}}_{\star}$, and Taylor expanding
$\nabla\Phi_{g}$ around the solar location at leading order $\mathcal{O}(r/R)$
yields a single differential equation for the relative separation between the
pair, which is known as the tidal approximation

$\displaystyle\ddot{\textbf{{r}}}\equiv\ddot{\textbf{{R}}}-\ddot{\textbf{{R}}}_{\star}\approx-\frac{Gm_{\star}}{r^{3}}{\textbf{{r}}}+{\bf
T}\cdot{\textbf{{r}}};$ (3)

where

$\displaystyle T^{ij}=-\frac{\partial^{2}\Phi_{g}}{\partial x_{i}\partial
x_{j}}\bigg{|}_{{\textbf{{R}}}_{\star}},$ (4)

is a 3x3 tidal tensor. The presence of a tidal field induces a velocity
variation

$\displaystyle\Delta{\textbf{{v}}}=\int_{S}{\rm d}t\,{\textbf{{F}}}_{\rm
tide}(t),$ (5)

where ${\textbf{{F}}}_{\rm tide}={\bf T}\cdot{\textbf{{r}}}$ is the tidal
force, and $S$ denotes the trajectory of an infalling object from
$t\to-\infty$ to the present. Hence, in general the orbital energy of
infalling Milky Way particles is not a conserved quantity

$\displaystyle\Delta
E=E-E_{\infty}=\Delta{\textbf{{v}}}\cdot{\textbf{{v}}}+\frac{1}{2}|\Delta{\textbf{{v}}}|^{2}.$
(6)

here, $E=v^{2}/2-\Phi_{\star}$ is the binding energy in a Keplerian potential
$\Phi_{\star}=-Gm_{\star}/r$, and $r$ and $v$ are the relative distance and
velocity to the point-mass $m_{\star}$. Crucially, Equations (5) and (6) show
that, depending on their individual trajectories, Galactic tidal forces can
either accelerate ($\Delta E>0$) or decelerate ($\Delta E<0$) infalling bodies
approaching the point-mass $m_{\star}$. This mechanism is therefore akin to
slingshot manoeuvres, which can be designed to increase/decrease the speed of
a spacecraft or redirect its path. In this paper, we are interested in objects
that lose enough kinetic energy as to become bound to the solar system, such
that $E<0$. For clarity, we shall refer to this particular type of capture
process as tidal trapping. The questions we aim to answer are the following:
how many interstellar bodies are expected to be tidally trapped on bound
orbits around the Sun at present? what is their baryon-to-dark matter ratio?
how are they distributed across the Solar system? how long do they remain
energetically bound?

Unfortunately, tidally-trapped particles move on very intricate trajectories
in the proximity of the Sun, hindering an analytical approach to these
questions. For illustration, Fig. 1 shows the motion of three particles that
become temporarily bound to a massive point-mass moving on a circular orbit
around the host galaxy potential (see §3 for details on the numerical
calculations). All trapped particles found in our numerical analysis remain
bound to the Keplerian potential for a limited amount of time and therefore
represent transient captures. We also found that the trajectories of these
objects can be broadly separated in three categories: a) Transient capture. In
the majority of capture events, interstellar particles only remain bound to
the point-mass $m_{\star}$ for the duration of a single encounter. This is
illustrated in the top panels of Fig. 1, which colours points along the
trajectory with negative and positive binding energies with red and black
lines, respectively. Top-left panel shows a particle with an original
galactocentric apocentre of $\sim 9\,{\rm kpc}$ being decelerated during a
single interaction with $m_{\star}$. As result, orbital energy is lost and the
particle is transferred to a new galactocentric orbit with apocentre $\sim
7\,{\rm kpc}$. The complexity of the fly-by approach is highlighted in the
middle panel, which shows the trajectory in the reference frame of the point-
mass $m_{\star}$. Notice that the binding energy only turns negative close to
perihelion ($t=0$), suggesting that the same dynamical mechanism that allows
this particle to lose kinetic energy in the first place is also responsible
for its later unbinding. b) Semi-stable capture (middle row). In this case,
tidally-trapped particles undergo several revolutions around the point-mass
$m_{\star}$ before being released back to the Galactic potential. These
objects are particularly interesting given that the probability of collisions
with planets is proportional to the time that an interstellar particle spends
in the solar system. c) Repeated capture (bottom row). The most convoluted
trajectories corresponds to particles that move in the Galactic potential and
experience periodic interactions with the point-mass $m_{\star}$. The sign of
the binding energy can be seen to alternate between positive and negative
values within short time intervals. A crucial result highlighted in Fig. 1 is
that tidal trapping always lead to temporary captures. We study to this
important aspect below with the aid of numerical experiments.

Figure 2: Left panel: Specific binding energy, $E=v^{2}/2-Gm_{\star}/r$, as a
function of time for the particles shown in Fig. 1. Energy is measured in
units of the escape velocity $v_{e}(r)=\sqrt{2Gm_{\star}/r}$ measured at the
tidal radius $r=r_{t}$, Equation (43), whereas time is normalized by the
fluctuation mean-life, $T(d)$, Equation (15) (see text). Right panel: Relative
distance and velocity of the three orbits plotted in the left panel. The
dashed line marks the value of the escape velocity as a function of distance,
and the vertical arrow shows the location of the tidal radius. Points
above/below the dashed line have positive/negative binding energies. In this
diagram, particles approach the point-mass from the upper-right corner, and
leave the system following a similar trajectory. Note that particles with a
constant energy move parallel to the dashed line (a.k.a ‘gravitational
focusing’ effect).

For a better understanding of the mechanics of tidal injection it is useful to
examine how orbital energy varies during the capture process. This is plotted
in the left panel of Fig. 2 for the three trajectories shown in Fig. 1. The
first noteworthy result is the confirmation that binding energy $E$ can not be
assumed to remain constant during capture. Indeed, all trapped particles show
very strong variations of energy as they approach the point-mass ($E$ sign
becomes negative), and after they are released back to the Galactic potential
($E$ sign becomes positive). The variation of energy is much smaller after the
injection onto the Keplerian system has taken place, but even then it does not
remain constant. As expected, the amount of time that a particle remains bound
to the point mass $m_{\star}$ scales directly with the binding energy attained
during injection. In particular, the trajectories of weakly-bound particles a)
(‘transient capture’) and c) (’repeated capture’) remain bound for a shorter
time span than b) (’semi-stable capture’), which reaches much higher binding
energies.

Where does capture take place? To answer this question, we plot in the right
panel of Fig. 2 the phase-space coordinates of trajectories shown in Fig. 1,
where $r$ and $v$ respectively denote distances and velocities relative to the
point mass $m_{\star}$. In this diagram, interstellar particles approach
$m_{\star}$ from the upper-right corner, and leave the Keplerian system
approximately following the inverse trajectory. At small separations, trapped
particles display an anticorrelation between distance and velocity expected
from gravitational focusing. However, at large separations the relative
distance and the relative velocity decrease as particles move towards the
point-mass, which is not what one would naively expect from gravitational
focusing. Indeed, if $E$ were a conserved quantity the relative velocity
should converge towards an asymptotic limit $v_{\infty}$ at $r\to\infty$ as
$v^{2}=v^{2}_{\infty}+2Gm_{\star}/r$. In contrast, Fig. 2 shows shows no
asymptotic behaviour, which casts doubts on the physical meaning of
$v_{\infty}$ for systems that are not in isolation. The transition between
these two regimes roughly happens at the tidal radius of the point-mass,
$r_{t}$ (marked with a vertical arrow, see §3 for details on its derivation).

It is important to stress that interstellar particles cannot visit the area
under the escape velocity curve, $v<v_{e}$ (black-dashed line) if the solar
system is assumed to be in isolation. The key difference in our models is
therefore the galactic potential, which sources a tidal field that actively
injects & removes interstellar particles from this region of phase-
space222From the point of view of statistical mechanics, tidal forces act as a
‘cooling’ mechanism that lowers the kinetic energy of infalling interstellar
objects. Notice, however, that heat transport goes in the opposite direction
to what is expected in fluids and gases, i.e. from the relatively cold
Keplerian system towards a relatively hot interstellar environment. This is
due to the negative specific heat of gravitating systems (Antonov 1961; see
Lynden-Bell 1999 for a review)., as illustrated in Fig. 2.

Equations (1) and (2) correspond to a special case of the restricted three-
body problem, where the most massive body (the host Galaxy) has an extended
mass distribution. Capture processes in the ‘classical’ three-body problem,
where the three bodies are represented as point-masses, have been studied in
the past. For example, the pioneer work of Szebehely (1967) showed that a
finite number of solutions exists where the lightest object is transferred
from one distinct mode of motion around the most massive body to another
distinct mode around the intermediate-mass one. This mechanism was used by
Hunter (1967), Heppenheimer (1975) and Heppenheimer & Porco (1976) to study
the formation of Jupiter’s outer satellites, and by Singer (1968) to analyze
the origin of the Moon as a planetoid captured by the Earth. More recently,
Suetsugu & Ohtsuki (2013) investigate the general case of temporary capture of
planetesimals by a giant planet, while Jílková et al. (2015) explore the
scenario where the inner Oort Cloud was captured from another star during a
close encounter in their birth cluster. Following up on the above results,
Higuchi & Ida (2017) analyze capture of asteroids by a planet moving on an
eccentric orbit, finding that temporary capture becomes more difficult as the
planet’s eccentricity increases. Daniel et al. (2017) analyzed the opposite
process, namely the escape of energetically-unbound particles from their
progenitor stellar cluster. It was found that a number of “potential escapers”
remained within two Jacobi radii for several dynamical periods, even though
their binding energy remained positive during this time. Recently, three-body
captures in accretion discs have also gained attention as a possible source of
Black-Hole binaries. E.g Li et a. (2022) and Boekholt et al. (2022) show that
close encounters between Black Hole pairs moving on circular disc orbits
around a supermassive Black Hole can form bound pairs with the help of
gravitational wave dissipation.

Here, it is worth highlighting the work of Petit & Henon (1986), who used
numerical method to follow the interaction of two small satellites, both
initially moving on circular orbits around a Keplerian potential. Two
important results in that paper are relevant in our analysis: first, their
work confirmed earlier claims that captures in a three-body system are
temporary events that ultimately end up the dissolution of the bound pair.
Furthermore, it was found that capture only happens for very precise
combinations of impact parameters and relative velocities which exhibit a
self-similar, Cantor-like structure. More recent work of Boekholt et al.
(2022) confirmed these results and found that the phase space structure that
leads to capture resembles a Cantor set with a fractal dimension of 0.4.

Taking these results at face value portrays the outer solar system as a
dynamically active region full of interstellar objects trapped on intricate,
unstable orbits that cross the tidal radius repeatedly and blur the boundary
between the solar system and the Milky Way. Alas, following the capture and
subsequent escape of tidally-trapped particles via numerical solutions to the
equations of motion (1) and (2) quickly becomes a CPU-demanding task when one
considers the astronomically-large number of minor bodies populating the
interstellar space.

As an alternative, this paper shows that a deeper theoretical understanding of
ensemble properties of tidally-trapped bodies can be obtained by abandoning
the classical approach of solving deterministic equations of motion for
individual objects, and focusing instead on a statistical description of the
phase-space distribution of interstellar particles in the vicinity of a point-
mass $m_{\star}$. Section 2 introduces a statistical theory that allow us to
compute the number of tracer particles trapped by a point-mass $m_{\star}$ as
it travels through an extended galaxy, and derive their steady-state
distribution in the Keplerian potential $\Phi_{\star}(r)=-Gm_{\star}/r$.
Section 3 shows a number of numerical experiments devised to test the theory.
In these experiments, a massive point-mass $m_{\star}$ captures tracer
particles from a spherical, non-rotating halo initially in dynamical
equilibrium. Section 4 discusses the implication of our findings in the
context of the solar system. In particular, we estimate the number of
Interstellar Objects (ISOs) and DM particles trapped on bound orbits around
the Sun.

##  2 Statistical theory

This Section presents a statistical theory that allows us to (i) estimate the
umber of particles that become gravitationally bound to a point-mass
$m_{\star}$ (the “Sun”) as it travels through a galaxy made of lighter
particles with individual masses $m\ll m_{\star}$ (the “galaxy”), as well as
(ii) derive their steady-state distribution in the Keplerian potential
$\Phi_{\star}(r)=-Gm_{\star}/r$.

Analytical solutions to this problem can be found using classical stochastic
techniques insofar as field particles in the vicinity of $m_{\star}$ can be
assumed to follow random trajectories. As shown below, this condition is met
in the collisionless limit, in which field particles move on unperturbed
orbits around a smooth Galactic potential $\Phi_{g}$. From a statistical
framework, the collisionless approximation is accurate insofar as the number
of particles affected by $m_{\star}$ represents a very small fraction of the
total number of field particles located within a volume element $V$ centred
around the point-mass at any given time. We will see that this condition
requires the local velocity dispersion of the host galaxy ($\sigma$) to be
much higher than the escape speed of the Keplerian potential ($v_{e}$), i.e.
$v_{e}\ll\sigma$, such that adding or removing the point-mass $m_{\star}$
leads to negligible effects on the local distribution function of the host
galaxy.

For clarity, the problem is broken in three parts. First, in §2.1 derive an
analytical equation for the number of field particles with relative velocities
$v<v_{e}$ within a volume $V$ around point-mass ($N_{b}$).

Second, in §2.2 we assume that particles move on statistically uncorrelated
(random) trajectories, which demands $N_{b}\ll N=n\,V$, where $n$ is the local
number density of field particles. This condition allows us to use
Smoluchowski (1916) statistical theory to derive an analytical expression for
the rate at which field particles enter $V$ with negative binding energies,
$E<0$.

The third and most difficult aspect of the analysis is to account for the
impact of $m_{\star}$ on the trajectories of infalling objects. To this aim,
§2.3 develops a statistical description of the limited life time of tidally-
trapped particles in the Keplerian potential $\Phi_{\star}$, i.e. the time
over which a trapped particle continuously has $E<0$ before being lost to
galactic tides. This in turn can be used to estimate the number of bound
interstellar particles in steady-state ($N_{\rm ss}$). We will see that that
$N_{\rm ss}=\alpha\,N_{b}$, where $\alpha$ is a dimension-less, order-unity
parameter that is empirically computed from the distribution of survival times
in Section 3. Lastly, §2.4 derives the equilibrium distribution of field
particles with $E<0$ orbiting around $\Phi_{\star}$.

We will also study the limitations of the theory in detail. Of particular
relevance is the singular behaviour of the Keplerian forces at $r=0$, which
translates in to a escape velocity that diverges
$v_{e}=2Gm_{\star}/r\to\infty$ in the limit $r\to 0$. This shortcoming is
analyzed in Section 3 with the aid of $N$-body experiments.

###  2.1 Bound particles in the vicinity of the point-mass

Let us first estimate the number of background particles expected to be found
in a volume element centred at the point-mass ($V$) with relative velocities
below the escape speed, $v<v_{e}$ This is done under two assumptions: (i)
field particles move on uncorrelated trajectories, and (ii) their distribution
in phase-space follows the local distribution function of the host galaxy.
Both conditions are valid in the collisionless limit, where the presence of
$m_{\star}$ induces negligible effects on the local galactic
background333Under this condition, the boundary $E=0$ loses any dynamical
meaning and the problem becomes purely statistical, i.e. it reduces to
counting particles with relative distances ($r$) and velocities ($v$) such
that $v<v_{e}(r)$ (or $E<0$), which we use as the main criterium for capture..
We will re-visit this approximation in Section 2.3, which accounts for the
effect of the gravitational attraction of $m_{\star}$ on the motion of bound
particles.

For simplicity, in what follows the host galaxy is assumed to be made of light
particles ($m\ll m_{\star}$) in a state of dynamical equilibrium, such that
the number density can be derived from the constant local matter density as
$n(\textbf{{R}})=\rho(\textbf{{R}})/m$. In our notation, phase-space
quantities measured in the galaxy frame are shown with capital letters, and
those measured with respect to the point-mass with lower letters.

The size of the volume element centred at $\textbf{{R}}_{\star}$ is to be
sufficiently small to guarantee that $r\ll
d(\textbf{{R}}_{\star})\equiv|\nabla\rho/\rho|^{-1}_{\textbf{{R}}_{\star}}$,
such that the number density can be assumed to be approximately constant, i.e.
$n({\textbf{{R}}_{\star}}+{\textbf{{r}}})\approx
n({\textbf{{R}}_{\star}})\equiv n$ for $r\ll d$. Hence, the distribution of
background particles becomes homogeneous across $V$, i.e. $p(\textbf{{r}}){\rm
d}^{3}r\approx 4\pi r^{2}\,n\,{\rm d}r$ at $r\ll d$, which is known as the
local approximation.

The relative velocity distribution of particles within the volume element $V$
is assumed to follow a Maxwellian distribution displaced by the reflex
velocity of the point-mass

$\displaystyle
p(\textbf{{v}})=\frac{1}{(2\pi\sigma^{2})^{3/2}}\exp\bigg{[}-\frac{(\textbf{{v}}+\textbf{{V}}_{\star})^{2}}{2\sigma^{2}}\bigg{]},$
(7)

where $\sigma=\sigma(\textbf{{R}}_{\star})$ is the local, one-dimensional
velocity dispersion of the host galaxy. The mean squared (relative) velocity
between the background particles and the point-mass $m_{\star}$ can be
straightforwardly derived from (7) as

$\displaystyle\langle v^{2}\rangle$ $\displaystyle=\int{\rm
d}^{3}v\,p(\textbf{{v}})\,v^{2}=3\sigma^{2}+V_{\star}^{2}.$ (8)

Following Peñarrubia (2021, hereafter P21), the number of background particles
within the volume $V$ with specific energy $E=v^{2}/2-Gm_{\star}/r<0$ can be
can be calculated under the local approximation and Maxwellian approximations
as

$\displaystyle N_{b}(r)$ $\displaystyle=\int_{V}{\rm
d}^{3}r\,n(\textbf{{r}})\int_{E<0}{\rm d}^{3}v\,p(\textbf{{v}})$ (9)
$\displaystyle=\frac{8\pi^{2}n}{(2\pi\sigma^{2})^{3/2}}\int_{0}^{r}{\rm
d}r^{\prime}\,r^{\prime 2}\int_{0}^{v_{e}(r^{\prime})}{\rm
d}v\,v^{2}\int_{-1}^{+1}{\rm
d}x\,e^{-\frac{v^{2}+V_{\star}^{2}+vV_{\star}x}{2\sigma^{2}}}$
$\displaystyle=2\sqrt{2\pi}\frac{n}{\sigma^{3}}\int_{0}^{r}{\rm
d}r^{\prime}\,r^{\prime 2}\int_{0}^{v_{e}(r^{\prime})}{\rm
d}v\,v^{2}\frac{\sigma^{2}}{vV_{\star}}e^{-(v+V_{\star})^{2}/(2\sigma^{2})}\bigg{(}e^{2vV_{\star}/\sigma^{2}}-1\bigg{)}$
$\displaystyle=2\sqrt{2\pi}\frac{n}{\sigma^{3}}\int_{0}^{r}{\rm
d}r^{\prime}\,r^{\prime
2}\bigg{\\{}-\frac{\sigma^{4}}{V_{\star}}e^{-(v_{e}+V_{\star})^{2}/(2\sigma^{2})}\bigg{(}e^{2v_{e}V_{\star}/\sigma^{2}}-1\bigg{)}+\sqrt{\frac{\pi}{2}}\sigma^{3}\bigg{[}{\rm
erf}\bigg{(}\frac{V_{\star}+v_{e}}{\sqrt{2}\sigma}\bigg{)}-{\rm
erf}\bigg{(}\frac{V_{\star}-v_{e}}{\sqrt{2}\sigma}\bigg{)}\bigg{]}\bigg{\\}}$

here ${\rm erf}(x)$ is the error function, and $v_{e}(r)=\sqrt{2Gm_{\star}/r}$
is the escape speed. Using the condition $v_{e}\ll\sigma$, the integrand term
within brackets can be approximated as

$\displaystyle-\frac{\sigma^{4}}{V_{\star}}e^{-(v_{e}+V_{\star})^{2}/(2\sigma^{2})}\bigg{(}e^{2v_{e}V_{\star}/\sigma^{2}}-1\bigg{)}+\sqrt{\frac{\pi}{2}}\sigma^{3}\bigg{[}{\rm
erf}\bigg{(}\frac{V_{\star}+v_{e}}{\sqrt{2}\sigma}\bigg{)}-{\rm
erf}\bigg{(}\frac{V_{\star}-v_{e}}{\sqrt{2}\sigma}\bigg{)}\bigg{]}=\frac{2}{3}e^{-V_{\star}^{2}/(2\sigma^{2})}v_{e}^{3}+\mathcal{O}\bigg{(}\frac{v_{e}}{\sigma}\bigg{)}^{5}.$
(10)

Thus, at leading order Equation (9) becomes

$\displaystyle N_{b}(r)$
$\displaystyle\approx\frac{4\sqrt{2\pi}}{3}e^{-V_{\star}^{2}/(2\sigma^{2})}\frac{n}{\sigma^{3}}\int_{0}^{r}{\rm
d}r^{\prime}\,r^{\prime 2}v_{e}^{3}(r^{\prime})$ (11)
$\displaystyle=\frac{32\sqrt{\pi}}{9}(Gm_{\star})^{3/2}\frac{n}{\sigma^{3}}e^{-V_{\star}^{2}/(2\sigma^{2})}r^{3/2}.$

which recovers Equation (4) of P21 for the case of a point-mass at rest
($V_{\star}=0$). The number of bound particles within the volume $V$ is
therefore proportional to the phase-space density of light background
particles at the location of the point-mass, $Q_{g}\equiv
Q_{g}(\textbf{{R}}_{\star})=n/\sigma^{3}$. The value of $N_{b}$ drops
exponentially for point-masses travelling across the background
($V_{\star}\neq 0$).

The above formula plays an important role in the remainder of this work, and
is accurate insofar as the local approximation holds, $r\ll d$, and the point-
mass is sufficiently light, such that the escape velocity
$v_{e}=(2Gm_{\star}/r)^{1/2}\ll\sigma$. This latter condition introduces a
critical radius ($r_{0}$) at small distances where the escape speed becomes
comparable to the velocity dispersion of the background particles. Setting
$v_{e}(r_{0})=\sigma$ yields $r_{0}=2Gm_{\star}/\sigma^{2}$. Therefore, the
statistical theory is valid within the range $r_{0}\lesssim r\ll d$, with the
lower limit approaching $r_{0}\to 0$ as $m_{\star}\to 0$.

Recall that the theory is built under the assumption that bound particles
represent a negligible fraction of the total number of field particles
enclosed within the same volume, $N=n\,V$. Using (11) it is straightforward to
show that the condition $v_{e}\ll\sigma$ meets this demand

$\displaystyle\frac{N_{b}}{N}\bigg{|}_{r}=\frac{2\sqrt{2}}{3\sqrt{\pi}}\bigg{(}\frac{v_{e}}{\sigma}\bigg{)}^{3}e^{-V_{\star}^{2}/(2\sigma^{2})}\lll
1~{}~{}~{}~{}~{}~{}{\rm for}~{}~{}~{}~{}\frac{v_{e}}{\sigma}\ll 1.$ (12)

###  2.2 Entrapment rate

Next goal is to compute the number of particles that enter the volume $V$ with
energy $E<0$ within a time interval, $t_{0},t_{0}+\Delta t$, where $\Delta t$
can be arbitrarily small. To this aim, it is useful to define the entrapment
rate as

$\displaystyle C_{\rm trap}\equiv\lim_{\Delta t\to 0}\frac{\Delta
N_{b}}{\Delta t}.$ (13)

The chief difficulty in deriving $\Delta N_{b}=N_{b}(t)-N_{b}(t_{0})$ from
Equation (11) lies in the statistical correlation between the location of
particles at two consecutive times, $t=t_{0}$ and $t_{0}+\Delta t$. Indeed, it
is easy to see that if the time interval becomes small, $\Delta t\to 0$, one
should expect a strong correlation between the values of $N_{b}$ measured at
the two instants of time because many particles found within the volume $V$ at
$t_{0}$ still remain in this region at $t_{0}+\Delta t$. Such correlations
vanish when the length of $\Delta t$ is larger than the average time that a
particle spends in the vicinity of $m_{\star}$. In this case, particles
identified at $t_{0}$ leave the region under observation and are replaced by a
new set of particles at $t_{0}+\Delta t$. Hence, consecutive measurements of
$N_{b}$ become statistically independent, and the problem can be treated
within a Markovian framework.

Here, the theory of Smoluchowski (1916) is used to compute the entrapment rate
(13) between arbitrarily-close time intervals. This theory was originally
devised to model the fluctuations of colloid concentrations in a liquid, and
later used by Chandrasekhar (1941, 1943) in the context of self-gravitating
systems. The general idea rests on the notion of “probability after-effect”,
or Wahrscheinlichkeitsnachwirkung, which describes the statistics of
correlated observations at different instants of time. In this paper, we use
the results of Smoluchowski (1916) without giving formal proof (interested
readers are directed to §III of Chandrasekhar 1943 for a detailed account of
the theory). The main assumptions are (1) that the motions of individual
particles are independent of each other, and (2) that all positions within the
volume $V$ have the same a priori probability to be sampled from an
homogeneous distribution. Under these conditions, Smoluchowski (1916) shows
that the probability $P_{N}(t)$ that at some later time there are still $N$
particles inside $V$ follows a law of decay that is analogous to the law of
decay of radioactive substances

$\displaystyle P_{N}(t){\rm d}t=e^{-t/T}{\rm d}t/T,$ (14)

with $T$ being the mean-life of a state in which the number of particles
within the volume $V$ remains constant. Under the local approximation, $r\ll
d$, the fluctuation mean-life roughly corresponds to the time that a particle
moving on a straight-line trajectory takes to cross the size of the sphere

$\displaystyle T(r)=\sqrt{\frac{2\pi}{3}}\frac{r}{\langle
v^{2}\rangle^{1/2}},$ (15)

where $\langle v^{2}\rangle^{1/2}$ is the average relative velocity between
$m_{\star}$ and the background particles, which can be estimated from (8)
under the Maxwellian approximation. Note that the time-scale (15) roughly
corresponds to the time that a field particle moving on a straight line with a
velocity $\langle v^{2}\rangle^{1/2}$ takes to cross the volume $V$.

Smoluchowski (1916) then shows that the probability that $N_{e}$ particles
enter the volume $V$ at any arbitrary time between $t_{0},t_{0}+\Delta t$
follows a Poisson distribution

$\displaystyle W(N_{e})=\frac{e^{-NP}(NP)^{N_{e}}}{N_{e}!},$ (16)

where $N=nV$ is the average number of particles in $V$, and $P$ is the so-
called probability after-effect factor

$\displaystyle P=\frac{\Delta t}{T},$ (17)

which represents the probability that a particle initially inside a given
volume emerges from $V$ within the time-interval $\Delta t$ (see also
Chandrasekhar 1943, p.53). As expected, when the time-interval is much shorter
than the fluctuation mean-life particles do not have sufficient time to exit
this region, hence the probability after-effect approaches $P\to 0$ in the
limit $\Delta t/T\to 0$. Direct inspection of the Poisson distribution (16)
shows that the average number of particles leaving the volume is
$N_{l}=NP=N\Delta t/T$. In dynamical equilibrium, this is equal to the number
of entering particles, hence $N_{e}=N_{l}=N\Delta t/T$.

Under the assumption that particles within the volume $V=4\pi r^{3}/3$ are
statistically uncorrelated, the above arguments must hold regardless their
binding energy to the point-mass $m_{\star}$.

The entrapment rate (13) at which particles with negative binding energies
enter the volume $V$ then follows from Equations (16) and (17), where $N$ is
replaced by $N_{b}$. Inserting (11) and (15) yields

$\displaystyle C_{\rm trap}(r)$ $\displaystyle=\lim_{\Delta t\to
0}\,N_{b}\frac{P}{\Delta t}=\frac{N_{b}}{T}$ (18)
$\displaystyle=\frac{16}{3}\bigg{(}\frac{2}{3}\bigg{)}^{1/2}\sqrt{\langle
v^{2}\rangle}(Gm_{\star})^{3/2}\frac{n}{\sigma^{3}}e^{-V_{\star}^{2}/(2\sigma^{2})}r^{1/2}.$

Time-integrating (18) returns a number of entrapment events that grows
linearly with time.

$\displaystyle N_{\rm trap}(r,t)$ $\displaystyle=\int_{0}^{t}{\rm d}t\,C_{\rm
trap}=N_{b}\,\frac{t}{T}$ (19)
$\displaystyle=\frac{16}{3}\bigg{(}\frac{2}{3}\bigg{)}^{1/2}\sqrt{\langle
v^{2}\rangle}(Gm_{\star})^{3/2}\frac{n}{\sigma^{3}}e^{-V_{\star}^{2}/(2\sigma^{2})}r^{1/2}\,t.$

Similarly, it is useful to compute the number of field particles leaving or
entering the volume $V$ in the time interval $t$

$\displaystyle N_{e}(r,t)$ $\displaystyle=N_{l}(r,t)=N\frac{t}{T}$ (20)
$\displaystyle=\bigg{(}\frac{8\pi}{3}\bigg{)}^{1/2}\sqrt{\langle
v^{2}\rangle}\,n\,r^{2}\,t,$

with $N=nV$.

Notice that (19) and (20) cross each other at a radius comparable to the
critical radius, $r_{0}=2Gm_{\star}/\sigma^{2}$. More precisely, solving
$N_{\rm trap}(r_{\epsilon},t)=N_{e}(r_{\epsilon},t)$ yields

$\displaystyle
r_{\epsilon}=\bigg{(}\frac{16}{9\pi}\bigg{)}^{1/3}e^{-V_{\star}^{2}/(3\sigma^{2})}\frac{Gm_{\star}}{\sigma^{2}}\simeq
0.4\,e^{-V_{\star}^{2}/(3\sigma^{2})}r_{0}.$ (21)

Below this “thermal” critical radius444Here, the word “thermal” highlights the
Maxwellian factor $\exp[-V_{\star}^{2}/(3\sigma^{2}_{\star})]$ multiplying the
critical radius. the theory is not valid, as it predicts that the number of
newly trapped particles exceeds the total number of particles entering the
volume, $N_{\rm trap}>N_{e}$ at $r<r_{\epsilon}$, a result with no physical
meaning within our idealized theoretical framework. This scale becomes
arbitrarily small as the attractor mass vanishes ($m_{\star}\to 0$), or its
speed increases ($V_{\star}\to\infty$).

The thermal critical radius $r_{\epsilon}$ defines the sphere of influence of
the point-mass $m_{\star}$, such that the dynamics of field particles found
within the volume $V_{\epsilon}=4\pi r^{3}_{\epsilon}/3$ are completely
dominated by the Keplerian potential, $\Phi_{\star}$. Indeed, comparison of
the mean kinetic energy of field particles $\mathcal{T}=3\sigma^{2}/2$ and the
potential energy at $r=r_{\epsilon}$ derived from Equation (21) yields

$\displaystyle\Phi_{\star}(r_{\epsilon})$
$\displaystyle=-\frac{Gm_{\star}}{r_{\epsilon}}$ (22)
$\displaystyle=-\frac{1}{2}\bigg{(}\frac{9\pi}{2}\bigg{)}^{1/3}e^{V_{\star}/(3\sigma^{2})}\,\sigma^{2}\approx-0.81\,e^{V_{\star}/(3\sigma^{2})}\,\mathcal{T}.$
(23)

Hence, particles located within a distance $r\lesssim r_{\epsilon}$ have
binding energies that exceed the mean kinetic energy of the field,
$|\Phi_{\star}|\gtrsim{\mathcal{T}}$. In general, the probability to find
unbound interlopers crossing the volume $V_{\epsilon}$ is small, and it
decreases further for point-masses that are not at rest with the galactic
background ($V_{\star}\neq 0$).

###  2.3 Survival

As outlined in the Introduction, tidal forces can both decelerate and
accelerate interstellar particles in the proximity of the point-mass
$m_{\star}$. This means that the same dynamical mechanism that allows
interstellar bodies to lose energy will also allow them to escape from the
Sun’s gravitational field. Following Napier et al. (2021b), it is useful to
define a dynamical lifetime function, $f_{\rm surv}(t)$, which determines the
fraction of captured objects that remain in bound orbits as a function of time
since entrapment occurred. Time-integrating (18) yields

$\displaystyle N_{\rm surv}(r,t)=\int_{0}^{t}{\rm d}t\,f_{\rm surv}(t)\,C_{\rm
trap}=N_{b}\,\frac{T_{s}}{T}\alpha(t),$ (24)

where $T_{s}$ is a characteristic survival time, and $\alpha(t)$ is a
dimensionless abundance parameter defined as

$\displaystyle\alpha(t)\equiv\frac{1}{T_{s}}\int_{0}^{t}{\rm d}t\,f_{\rm
surv}(t).$ (25)

Numerical tests in Section 3 show that Equation (25) converges asymptotically
to a value of order unity, $\alpha(t)\to\alpha$ at $t\to\infty$, such that the
number of surviving particles reaches a steady-state value

$\displaystyle N_{\rm surv}(r,t)\to N_{\rm
ss}(r)=N_{b}\,\frac{T_{s}}{T}\alpha~{}~{}~{}~{}~{}~{}~{}{\rm
for}~{}~{}~{}t/T_{s}\to\infty.$ (26)

Equation (26) exhibits three different limiting cases depending on the type of
trajectories on which tidally-trapped objects move:

  1. 1.

Unstable orbits. Binding energies only become negative during the duration of
a single encounter. The survival time in this case is comparable to the
average time that a random body spends within the volume $V$, i.e.
$T_{s}\approx T\approx r/\sqrt{\langle v^{2}\rangle}$. Setting $T_{s}=T$ in
Equation (26) returns an asymptotic, steady-state value $N_{\rm
ss}=N_{b}\,\alpha$ for $t\gg T_{s}$.

  2. 2.

Semi-stable orbits. Particles remain bound during several orbital revolutions
before being released back to the galactic potential. The survival time is
long but finite, $T_{s}\gg T$. Hence, the number of bound objects converges to
$N_{\rm ss}\approx N_{b}(T_{s}/T)\alpha\gg N_{b}$ on a time-scale $t\gg
T_{s}$.

  3. 3.

Stable orbits. Trapped objects remain bound for an arbitrarily long interval
of time. The survival time diverges $T_{s}/T\to\infty$, and the survival
fraction approaches $f_{\rm surv}\to 1$. Comparison of (19) and (24) yields
$N_{\rm surv}(t)=N_{\rm trap}(t)=N_{b}\,(t/T_{s})$, which grows linearly with
time. As a result, no steady-state exists in this limit.

In the collisionless approximation, interstellar particles cross the volume
$V$ on statistically uncorrelated trajectories. Under this condition, the
amount of time that a random field particle spends in the volume $V$ is
governed by Smoluchowski (1916) probability function (14). Accordingly, the
dynamical lifetime is expected to decay exponentially as $f_{\rm
surv}(t)=\exp(-t/T)$, where $T$ is the fluctuation mean-life (15). The
characteristic survival time that a random bound particle spends in the volume
$V$ is therefore comparable to the duration a fly-by encounter, $T_{s}=T$,
which means that tidally-trapped objects move on unstable orbits in the
Keplerian potential $\Phi_{\star}$. Inserting an exponential-decaying lifetime
function in Equation (25) and setting $T_{s}=T$ yields an abundance parameter

$\displaystyle\alpha=\lim_{t\to\infty}\frac{1}{T}\int_{0}^{t}{\rm
d}t\,\exp(-t/T)=1.$ (27)

Hence, the steady-state number of bound particles around the point-mass
$m_{\star}$ can be derived from (11) and (26) by setting $T_{s}=T$ and
$\alpha=1$, which yields

$\displaystyle N_{\rm ss}(r)$ $\displaystyle=N_{b}(r)$ (28)
$\displaystyle=\frac{32\sqrt{\pi}}{9}(Gm_{\star})^{3/2}\frac{n}{\sigma^{3}}e^{-V_{\star}^{2}/(2\sigma^{2})}r^{3/2}~{}~{}~{}~{}~{}~{}~{}{\rm
for}~{}~{}~{}~{}t\gg T.$

As shown in §2.2, these predictions fail on scales comparable or smaller than
the thermal critical radius, where Equation (28) predicts a number of bound
particles that exceeds the number of particles enclosed in the volume element
$V$. More precisely, solving $N_{\rm
ss}(r^{\prime}_{\epsilon})=N(r^{\prime}_{\epsilon})$ and inserting (21) yields
$r^{\prime}_{\epsilon}=2^{2/3}r_{\epsilon}$. Therefore, Equation (28) should
not be applied on scales below the thermal critical radius $r\lesssim
r_{\epsilon}$. Section 3 tests the above theoretical expectations with live
$N$-body models.

###  2.4 Steady-state profiles

The presence of trapped interstellar bodies around the central star leads to
the formation of an extended ‘halo’ that approaches an equilibrium state as
the number of particles being trapped equals the number being tidally
stripped.

In the collisionless limit, the steady-state density enhancement due to the
presence of bound interstellar material can be directly derived from (28) as

$\displaystyle\delta(r)$
$\displaystyle\equiv\frac{1}{4\pi\,n\,r^{2}}\frac{{\rm d}N_{\rm ss}}{{\rm
d}r}$ (29)
$\displaystyle=\frac{4}{3\sqrt{\pi}}\frac{(Gm_{\star})^{3/2}}{\sigma^{3}}e^{-V_{\star}^{2}/(2\sigma^{2})}\frac{1}{r^{3/2}},$

which scales as $\delta\sim(v_{c}/\sigma)^{3}\sim r^{-3/2}$, where
$v_{c}(r)=\sqrt{Gm_{\star}/r}$ is the circular velocity at a distance $r$ from
the Keplerian potential. This profile – generally known as a ‘density spike’–
was originally found by Gondolo & Silk (1999) by analyzing the spatial
distribution of particles with a constant phase-space density around a Black
Hole in isolation.

Of particular relevance is the distance at which the density of trapped
particles equals the background interstellar density. Using (29), it is
straightforward to show that this happens at the thermal critical radius,
$\delta(r_{\epsilon})=1$, which calls for caution when extrapolating Equation
(29) at $r\lesssim r_{\epsilon}$ (or $\delta\gtrsim 1$). We will come back to
this issue in Section 3.

The velocity dispersion of the bound halo is found from (29) by solving the
isotropic Jeans equations

$\displaystyle\sigma_{h}^{2}(r)$
$\displaystyle=\frac{1}{\delta(r)}\int_{r}^{\infty}{\rm
d}r^{\prime}\delta(r^{\prime})\bigg{|}\frac{{\rm d}\Phi_{\star}}{{\rm
d}r}\bigg{|}$ (30) $\displaystyle=\frac{2}{5}\frac{Gm_{\star}}{r},$

which corresponds to a constant fraction of the circular velocity,
$\sigma_{h}/v_{c}=(2/5)^{1/2}\approx 0.632$ at all radii.

It is straightforward to show that the mean phase-space density of halo
particles is constant across the volume $V$. Combination of (29) and (30)
yields

$\displaystyle Q_{h}$
$\displaystyle\equiv\frac{n\,\delta(r)}{\sigma^{3}_{h}(r)}$ (31)
$\displaystyle=\frac{5}{3}\bigg{(}\frac{10}{\pi}\bigg{)}^{1/2}e^{-V_{\star}^{2}/(2\sigma^{2})}\,Q_{g},$

where $Q_{g}=n/\sigma^{3}$ is the field phase-space density. Interestingly,
the ratio $Q_{h}/Q_{g}$ solely depends on the speed of the attractor relative
to the host galaxy, $V_{\star}/\sigma$. Interstellar particles bound to a
point-mass at rest have higher phase-space densities than those in the
Galactic background, $Q_{h}(V_{\star}=0)\simeq 2.97\,Q_{g}$. The halo mean
phase-space density drops for moving objects, such that $Q_{h}<Q_{g}$ for
objects travelling with a speed $V_{\star}>1.476\,\sigma$.

###  2.5 Orbits in isolation

Up to this point, we have not made any explicit assumption on the motion of
particles trapped in the halo. This Section derives the equilibrium
distribution of a large ensemble of these objects orbiting around a point-mass
in isolation, such that the only force acting on bound particles corresponds
to the gravitational attraction of $m_{\star}$. Invoking isolation is a useful
theoretical fudge that allows us to compute ensemble-averaged properties of
trapped particles as a superposition of orbits that conserve specific energy
$E=v^{2}/2-Gm_{\star}/r$ and angular momentum
$\textbf{{L}}=\textbf{{r}}\times\textbf{{v}}$, making the problem
mathematically tractable. However, one should bear in mind that this
assumption is particularly poor for weakly-bound particles, as illustrated in
Fig. 2.

Following Equation (31), we assume that energetically-bound particles are
distributed homogeneously in phase-space on distance scales $r\ll d$, with a
distribution function that is approximately constant across the volume $V=4\pi
r^{3}/3$, i.e. $f({\textbf{{r}}},{\textbf{{v}}})=f_{0}$. The spectrum of semi-
major axes and eccentricities is found by mapping points in integral-of-motion
space onto the phase-space volume ${\rm d}^{6}\Omega={\rm d}^{3}r\,{\rm
d}^{3}v$ by means of the Jacobian $\omega(E,L)=|\partial^{6}\Omega/\partial
E\partial L|=8\pi^{2}L\,P$, also known as density of states (e.g. see Appendix
A of P15), where $P=2\pi a^{3/2}/\sqrt{Gm_{\star}}$ is the orbital period and
$L=[Gm_{\star}\,a(1-e^{2})]^{1/2}$ is the angular momentum of particles with a
given energy. Writing ${\rm d}E=Gm_{\star}{\rm d}a/(2a^{2})$ and $L\,{\rm
d}L=-Gm_{\star}\,a\,e{\rm d}e$ yields ${\rm d}^{6}\Omega=8\pi^{2}L\,P\,{\rm
d}E{\rm d}L=8\pi^{3}(Gm_{\star})^{3/2}\,e\,a^{1/2}\,{\rm d}(-e){\rm d}a$.
Hence,

$\displaystyle\omega(a,e)=8\pi^{3}(Gm_{\star})^{3/2}\,e\,a^{1/2},$ (32)

which recovers Equation (19) of Dehnen et al. (2021). The probability to find
particles with semi-major axes and eccentricities in the intervals $(a,a+{\rm
d}a)$ and $(e,e+{\rm d}e)$ is

$\displaystyle
p(a,e)=\omega(a,e)\,f(a,e)=8\pi^{3}(Gm_{\star})^{3/2}\,e\,a^{1/2}\,f_{0}.$
(33)

Equation (33) reveals a few points of interest. Notice first that the
distribution of integrals has a separable form, $p(a,e)=p(a)p(e)$, with a
semi-major axis distribution that follows a power-law

$\displaystyle p(a){\rm d}a\sim a^{1/2}{\rm d}a.$ (34)

The fact that $p(a)\to 0$ in the limit $a\to 0$ highlights the low probability
to trap particles onto tightly-bound orbits by chance. The eccentricity
distribution in Equation (33) is ‘thermal’

$\displaystyle p(e){\rm d}e=2e\,{\rm d}e,$ (35)

which indicates that trapped particles tend to fall in on eccentric orbits. As
a sanity check, Appendix A shows that the number density profile associated
with the distribution function (33) recovers Equation (29) derived from
Smoluchowski’s (1916) statistical theory. Appendix B proves that Equations
(33) also describes the orbital distribution of wide binaries that form via
chance entrapment of uncorrelated stars with Maxwellian velocities.

The closest approach of a trapped object to the sun (i.e. the orbital
pericentre) is related to the semi-major axis and eccentricity as
$r_{p}=a(1-e)$. The fraction of objects within a pericentre interval
$r_{p},r_{p}+{\rm d}r_{p}$ can be straightforwardly calculated from (35) as
(c.f. eqs 3 and 4 of Hills 1981)

$\displaystyle p(r_{p})\,{\rm
d}r_{p}=\frac{2}{a}\bigg{(}1-\frac{r_{p}}{a}\bigg{)}\,{\rm d}r_{p}.$ (36)

Notice that for $r_{p}\ll a$ the distribution $p(r_{p})$ becomes independent
of $r_{p}$, which therefore implies that trapped halo particles reaching the
inner solar system have a homogeneous distribution of perihelia. This has
important implications for the detection of trapped interstellar objects from
Earth, as discussed in Section 4.

##  3 Numerical tests

Given the heuristic nature of the theory presented in §2, it is useful to run
controlled $N$-body tests that falsify two key assumptions on which the theory
rests. Namely, we use galaxy models where (i) the background density is not
homogeneous, and (ii) the velocity function is not Maxwellian. In these
experiments, the motion a point-mass $m_{\star}$ (the “Sun”) is followed as it
travels on a circular orbit through a self-gravitating Dehnen (1993) sphere in
dynamical equilibrium (the “DM halo”). At each time-step, we identify
particles that become gravitationally bound to the point mass (i.e. binding
energy flips from positive to negative), compute their orbital elements in the
Keplerian potential $\Phi_{\star}$, and count the number of particles stripped
by Galactic tides (i.e. binding energy flips from negative to positive).

To compute the orbits of the point-mass $m_{\star}$ and halo particles in the
host galaxy, we solve the equations of motion (1) and (2), respectively, using
a leap-frog technique with varying time-steps (e.g. Press 1992). Details on
the numerical integration are given in Appendix C.

In addition, for a better understanding of the effect of $m_{\star}$ on the
local galactic background we run $N$-body experiments where the gravitational
force induced by the point-mass is removed from the equations of motion
($F_{\star}=0$). Hence, in these so-called “collisionless” experiments field
particles follow orbits in a smooth galaxy potential $\Phi_{g}$ and do not
feel the gravitational attraction of the point-mass, $m_{\star}$. Their phase-
space distribution obeys the collisionless Boltzmann equation, which provides
a closer match to the idealized conditions on which the statistical theory
presented in §2 rests.

###  3.1 Set up

The first step of this Section is to generate $N$-body realizations of a
Dehnen (1993) sphere in dynamical equilibrium. The gravitational potential of
these objects has an analytical form

$\displaystyle\Phi_{g}(R)=\frac{4\pi
G\rho_{0}}{3-\gamma}\times\begin{cases}-\frac{1}{2-\gamma}\big{[}1-\big{(}\frac{R}{R+R_{0}}\big{)}^{2-\gamma}\big{]}&,\gamma\neq
2\\\ \ln\big{(}\frac{R}{R+R_{0}}\big{)}&,\gamma=2.\end{cases}$ (37)

where $R$ is the galactocentric radius. The total mass is
$M_{g}=4\pi\rho_{0}R_{0}^{3}/(3-\gamma)$, and the number density profile is

$\displaystyle n(R)=\frac{n_{0}}{(R/R_{0})^{\gamma}[1+(R/R_{0})^{4-\gamma}]},$
(38)

with $n_{0}$ chosen such that $4\pi\int_{0}^{\infty}{\rm
d}r\,r^{2}\,n(r)=N_{g}$. The density scale-length can be directly derived from
(38)

$\displaystyle d(R)=\bigg{|}\frac{\nabla
n}{n}\bigg{|}^{-1}_{R}=R\frac{R+R_{0}}{4R+\gamma R_{0}}.$ (39)

Recall that at a given galactocentric radius, the local approximation is valid
on distance scales $r\ll d$, where $r$ is measured relative to the point-mass
location.

To assign orbital velocities in a way that guarantees dynamical equilibrium,
the distribution function is calculated using Eddington (1916) inversion

$\displaystyle f(E)=\frac{1}{\sqrt{8}\pi^{2}}\bigg{[}\int_{E}^{0}\frac{{\rm
d}\Phi}{\sqrt{\Phi-E}}\frac{{\rm d}^{2}n}{{\rm
d}\Phi^{2}}+\frac{1}{\sqrt{-E}}\bigg{(}\frac{{\rm d}n}{{\rm
d}\Phi}\bigg{)}_{\Phi=0}\bigg{]},$ (40)

where $n[R(\Phi_{g})]$ corresponds to the profile (38) expressed as a function
of the potential (37). Note that the velocity distribution resulting from
Eddington’s equation (40) is not Maxwellian. To generate halo models in
equilibrium, positions and velocities for $N_{g}=5\times 10^{7}$ particles are
drawn randomly from the distribution function (40), with unity vectors
isotropically distributed over the surface of a sphere.

To guide the physical interpretation of our results, we choose a dark matter
halo with a mass and scale radius that match those of the Aquarius simulation,
$M_{g}=1.84\times 10^{12}M_{\odot}$, $R_{0}=15.3\,{\rm kpc}$ and $\gamma=1$
(Springel et al. 2008). Subsequently, we place point-mass attractors at a
fixed distance from the halo centre, $R_{\star}=0.065\,R_{s}$, with a circular
velocity $V_{\star}=V_{c}(R_{\star})=172.4\,{\rm km\,s^{-1}}$. The density
scale-length (39) at this galactocentric radius is $d(R_{\star})=0.84\,{\rm
kpc}$, and the local velocity dispersion is $\sigma(R_{\star})=196\,{\rm
km\,s^{-1}}$. Under the Maxwellian approximation, the average relative
velocity between the point-mass and the surrounding DM particles is therefore
$\sqrt{\langle v^{2}\rangle}=\sqrt{V_{\star}^{2}+3\sigma^{2}}\simeq 381\,{\rm
km\,s^{-1}}$. The fluctuation mean-life (15) at a relative distance $r=d$ is
$T(d)=3.12\,{\rm Myr}$, whereas the mean phase-space density of trapped halo
particles predicted by Equation (31) is $Q_{h}/Q_{g}\approx 2$.

For illustration, we consider three values for $m_{\star}$ that are
specifically chosen to test various aspects of the theory presented in §2.
Recall that the theoretical model works on distance scales where the local
approximation is valid, $r\ll d$, and the escape velocity is much smaller than
the velocity dispersion of the field, $v_{\rm esc}(r)/\sigma\ll 1$. With this
in mind, we choose $m_{\star}/\,M_{\odot}=10^{8},3\times 10^{8}$ and $10^{9}$,
which have characteristic escape velocities $v_{\rm
esc}(d)=(2Gm_{\star}/d)^{1/2}=\\{0.16,0.28,0.51\\}\,\sigma$, respectively. We
shall see below that the lightest model mostly evolves in a collisionless
regime, while the most massive one shows noticeable deviations from a
collisionless behaviour.

Finally, the integration time ($t_{f}$) is chosen to obtain a statistically-
meaningful number of trapping events. In particular, setting $t_{f}/T(d)=100$
in Equation (19) yields $N_{\rm trap}(d,t_{f})\sim\\{1.1,57.3,348.6\\}\times
10^{4}$ trapped particles for $m_{\star}/\,M_{\odot}=10^{8},3\times 10^{8}$
and $10^{9}$, respectively.

Notice that the $N$-body set-up is still far from reaching the resolution
required to model stellar-size objects. Unfortunately, given that the number
of bound particles (28) drops as $N_{\rm trap}\sim m_{\star}^{3/2}$, it
becomes numerically demanding to run experiments with lighter point-mass
attractors. For example, in order to find $N_{\rm trap}\sim 10^{2}$ objects
bound to a stellar-size object ($m_{\star}/M_{g}\sim 10^{-12}$) one would need
to increase the number of halo particles by $\sim 12$ orders of magnitude
(i.e. $N_{g}\sim 10^{20}$), which lies far beyond current computational
capabilities.

###  3.2 Results

Figure 3: Distribution of survival times, $t_{\rm surv}$, for collisionless
and collisional models (top and bottom panels, respectively) measured in units
of the fluctuation mean-life, $T$, Equation (15). Here, $t_{\rm surv}$ is
defined as the time interval over which a particle continuously has $E<0$
within the volume element $V=4\pi r^{3}/3$ centred at a point-mass
$m_{\star}$. Lines are colour-coded according to volume radius ($r$) given in
units of the local density scale-length $(d)$. Black-dashed lines show an
exponentially-decaying dynamical survival fraction (14) with a characteristic
survival time set equal to the fluctuation mean-life, $f_{\rm
surv}(t)=\exp(-t/T)$. Right panels plot the abundance parameter $\alpha$
derived from Equation (25) as a function of the volume size ($r$) measured in
units of the thermal critical radius ($r_{\epsilon}$), Equation (21).

####  3.2.1 Survival time

The theory presented in §2.2 assumes that the gravitational force of
$m_{\star}$ induces negligible perturbations on the trajectories of nearby
particles, which is the so-called “collisionless” approximation. Under this
condition, the amount of time that an interstellar object spends in the volume
$V$ follows Smoluchowski (1916)’s exponentially-decaying function (14),
$\exp(-t/T)$, where the time-scale $T$ is defined by Equation (15) and roughly
corresponds to the average time that a particle moving on a straight line with
a relative velocity $\langle v^{2}\rangle^{1/2}$ takes to cross the volume
$V$. Thus, in a collisionless framework the typical time-span that a particle
remains energetically bound to $m_{\star}$ is comparable to the duration of a
flyby encounter, whereas the probability to find bound particles with survival
times $t_{\rm surv}\gg T$ drops exponentially.

To test these expectations, Fig. 3 plots the dynamical lifetime function
($f_{\rm surv}$) of particles bound to three different point-masses
$m_{\star}$. Top and bottom panels respectively correspond to “collisionless”
models, in which the force induced by the point-mass is removed from the
equations of motion (2) ($F_{\star}=0$), and to “collisional” models, which do
account for the gravitational force of the point-mass ($F_{\star}\neq 0$).
Recall that $t_{\rm surv}$ is defined as the time over which a particle
located within a volume element $V=4\pi r^{3}/3$ continuously has $E<0$.

As expected, models that neglect the point-mass self-gravity lead to a
distribution of survival times that converges towards Smoluchowski (1916)
exponentially-decaying function, $\exp(-t/T)$ (black-dashed lines), in the
limit $r/d\to 0$. Consequently, the abundance parameter (25) plotted in the
right-most panel asymptotically approach unity, i.e. $\alpha\to 1$ as $r/d\to
0$. Notice that varying the attractor $m_{\star}$ has a visible impact on how
fast these models converge towards the collisionless limit, with more massive
models converging faster than low mass ones. Measuring the volume size $r$ in
units of the thermal critical radius $r_{\epsilon}$ defined in Equation (21)
reveals the self-similar behaviour of $\alpha$, which is solely controlled by
the ratio $r/r_{\epsilon}$, such that $\alpha\approx 1$ for
$r/r_{\epsilon}\lesssim 1$. Following Equation (28), this implies that the
steady-state number of bound particles in the volume element $V$ converges to
the value derived from the unperturbed distribution function of the host, i.e.
$N_{\rm ss}=N_{b}$ for $r/r_{\epsilon}\lesssim 1$.

Comparison with the corresponding collisional models in the bottom panels
reveals a number of interesting features. The first noticeable result is the
emergence of a sizeable population of particles that remain bound to the
Keplerian potential $\Phi_{\star}$ for time intervals that can be as long as
the integration time of the models, $t_{\rm surv}\sim t_{f}=100\,T$. This
translates into a distribution of survival times that departs from the
exponential shape predicted by Smoluchowski (1916) and becomes more akin to a
step function, in which most trapped particles move on unstable orbits that
are lost to galactic tides soon after being entrapped, $t_{\rm surv}\lesssim
T$, and a sizeable population of objects that move on semi-stable orbits with
long survival times $t_{\rm surv}\gg T$ (c.f. panel c) in Fig. 1). As
discussed in §4, long-live objects are particularly interesting for detection
experiments in the solar system.

Interestingly, the abundance parameters shown in the bottom right-most panel
of Fig. 3 also exhibit a self-similar behaviour when the distance is plotted
in units of the thermal critical radius, but the radial dependence has
noticeable differences with the values derived from collisionless simulations.
In particular, we find that the gravitational attraction of the point-mass
reduces the average survival time of trapped particles at small distances,
such that $\alpha\ll 1$ at $r\ll d$. This in turn implies a steady-state
number of bound particles $N_{\rm ss}\lesssim N_{b}$ at $r\lesssim
r_{\epsilon}$. On the other hand, $\alpha$ exhibits systematically larger
values than in collisionless models at distances $r\gtrsim r_{\epsilon}$,
suggesting that the gravitational attraction of the point-mass may affect the
overall spatial distribution of tidally-trapped particles. We inspect this
issue in §3.2.3.

Figure 4:  Time-evolution of the number of particles energetically-bound to
three different point-masses $m_{\star}$ within the volume element $V=4\pi
r^{3}/3$ for a fixed volume size $r=0.02\,d$. Time is given in units of the
fluctuation mean-life, $T(r)$, Equation (15). Red and blue symbols denote
collisional and collisionless models, respectively. In collisionless models,
the cumulative number of particles entering the volume $V$ (light-coloured
symbols) grows linearly with time and is accurately described by Equation (19)
on time scales $t\gtrsim T$ (black-dashed lines) . The number of bound
particles within $V$ (dark-coloured symbols) matches the steady-state value
predicted by Equation (28) (dotted lines). In collisional models, the number
of particle entering the volume $V$ also shows a linear growth with time, but
with values that lie systematically above the relation (19). The discrepancy
increases systematically for more massive point-masses $m_{\star}$. Similarly,
the population of bound particles is larger than expected from Equation (28).
Note that in the middle and right panels, the number of bound particles that
exceeds the total number of random field particles expected to be found within
$V$, i.e. $N_{b}\gtrsim N=n\,V$. This effect is caused by the self-gravity of
the point mass $m_{\star}$ and occurs in models with a volume size smaller
than the thermal critical radius, i.e. $r\lesssim r_{\epsilon}$ (see text).
Figure 5: Left panels: Time-evolution of the cumulative number of particles
that enter the volume element $V=4\pi r^{3}/3$ at the time $t_{f}=100\,T$ as a
function of volume size ($r$) normalized by the local density scale ($d$).
Here, time is given in units of the fluctuation mean-life, $T$, Equation (15).
Upper and bottom panels correspond to collisionless ($F_{\star}=0$) and
collisional ($F_{\star}\neq 0$) experiments, respectively. Right panels:
Number of particles with negative binding energies at $t=t_{f}$ for the models
shown in the left panels. Black-dotted lines correspond to the steady-state
population of bound particles derived in a collisionless regime ($N_{\rm
ss}$), Equation (28), whereas black-solid lines mark the relation $N=4\pi
r^{3}/3$ expected in a constant-density background of field particles. Dotted
and solid lines cross each other at at a distance $r=2^{2/3}\,r_{\epsilon}$,
where $r_{\epsilon}$ is the thermal critical radius, Equation (21) (marked
with vertical arrows for reference). In collisionless models (upper panel),
the number of bound particles cannot exceed the total number of particles in
the volume element, $N_{\rm ss}\leq N$. Most field particles entering the
radius $r\lesssim r_{\rm eps}$ are energetically bound. Bottom panels show
that the point-mass attraction removes this upper threshold in collisional
models due to the presence of particles with long survival times $t_{\rm
surv}\gg T$ (see Fig. 3).

####  3.2.2 Entrapment rates

Next step is to check the analytical capture rates derived in §2. For
illustration, we choose a relatively small radius ($r=0.02\,d$) centred around
$m_{\star}$, and record the number of particles that enter & leave the volume
$V=4\pi r^{3}/3$ centred around $m_{\star}$, and monitor their binding
energies to the Keplerian potential $\Phi_{\star}(r)=-Gm_{\star}/r$ at each
time step.

Fig. 4 plots the time-evolution of the cumulative number of particles that
enter the volume $V$ (light-coloured symbols), and the size of the population
of particles with negative binding energies $E<0$ (dark-coloured symbols).
Here, time measured in units of the fluctuation mean-life, $T(r=0.02d)$,
Equation (15). Red and blue symbols denote collisional and collisionless
models, respectively.

This Figure highlights a few important aspects of the problem. First, on
sufficiently long time-scales $t\gtrsim T$ the cumulative number of particles
crossing the volume $V$ grows linearly with time, which is to be expected from
a constant entrapment rate. However, while collisionless models follow very
closely the analytical expression (19), $N_{e}=N\,t/T$ (black-dashed lines),
in collisional simulations the number of field particles entering $V$ lies
systematically above the theoretical expectations. The fact that the
discrepancy grows in proportion to $m_{\star}$ indicates that the
gravitational attraction of the point mass causes a systematic increase the
local density of field particles, which in turn boosts the number of objects
moving through $V$.

Crucially, in all our models the number of bound particles within the volume
$V$ converges to a steady-state value $N_{\rm surv}(t)\to N_{\rm ss}$ on time-
scales $t\gtrsim T$. Here, the size of the bound population is measured as
$N_{\rm surv}(t)=N_{\rm trap}(t)-N_{\rm unb}(t)$ (dark-coloured symbols),
where $N_{\rm trap}(t)$ and $N_{\rm unb}(t)$ respectively are the cumulative
number of particles trapped by and stripped from the Keplerian potential
$\Phi_{\star}$ within a time interval $t$. For reference, solid lines mark the
number of random particles enclosed with the volume element $V$ in a constant-
density field, $N=n\,V$.

For light point-mass models ($m_{\star}=10^{8}\,M_{\odot}$, left panel), the
steady-state population of bound objects accurately matches the analytical
expression (28), $N_{\rm surv}(t)\approx N_{b}$ (dotted-black lines),
independently of whether or not the gravitational attraction of the point-mass
is included in the equations of motion. However, as $m_{\star}$ increases, we
start to observe noticeable differences between collisionless and collisional
models. Namely, whereas in collisionless experiments the number of trapped
particles is limited by the number of particles within a given volume element,
such that $N_{\rm surv}(t)\leq N$, collisional models do not obey this
threshold. Instead, the steady-state number of bound particles lies between
the number expected from a constant-density field and that predicted by
Equation (28), i.e. $N\lesssim N_{\rm surv}(t)\lesssim N_{b}$. This result can
be explained by the boosted rate of particles entering the volume $V$ in
collisional models, and the fact that the abundance parameter falls below
unity at small distances from the point-mass (see bottom-right panel in Fig.
3).

For a better understanding of the above results, Fig. 5 shows how the total
number of particles that enter the volume (left panels) and the size of the
bound population (right panels) vary as a function of volume radius ($r$). As
above, the we run models in both collisionless (upper panels) and collisional
modes (bottom panels).

The first noteworthy result in the upper-left panel is the excellent match
between the number of particles entering the volume measured in the
simulations and those predicted by Equation (20) (black-dashed line),
$N_{e}(t_{f})=N\,t_{f}/T$, over 9 orders of magnitude, which gives evidence
that Smoluchowski’s (1916) probability after-effect factor $T$ defined in
Equation (17) successfully accounts for the correlation in the number of
particles measured at consecutive time intervals. Note also that collisionless
models are insensitive to the value of $m_{\star}$ by construction. However,
this symmetry is broken in collisional models plotted in the bottom-left
panel, which show that the attractive force of the point-mass tends to
increase the number of particles at small distances with respect to the values
predicted by Equation (20), a discrepancy that grows as the value of
$m_{\star}$ increases.

Right panels show the number of bound particles found in the volume $V$ as a
function of volume size for the three point-masses considered above. In
collisionless models (upper-right panel), Equation (28) (black-dashed lines)
slightly underestimates the numerical values at large radii from the point-
mass. The source of the discrepancy can be traced back to the raising
abundance parameter away from the point-mass, $\alpha\gtrsim 1$ at $r\gtrsim
r_{\epsilon}$, shown in Fig. 3, whereas the analytical expression (28) adopts
$\alpha=1$. However, the main mismatch is found at small volume sizes, where
Equation (28) grossly overestimates the number of particles in the volume $V$
with $E<0$. This can be easily understood by noting that while the number of
bound particles scales as $N_{\rm ss}\sim r^{3/2}$, the total number of random
particles within a given volume element has a steeper radial dependence,
$N\sim r^{3}$. The two curves approximately cross at the thermal critical
radius, $r\approx r_{\epsilon}$, with $r_{\epsilon}$ given by Equation (21).
As discussed in §2.2, Equation (28) should not be extrapolated down to
arbitrarily small volume sizes around $m_{\star}$. A simple way to account for
this shortcoming is by modifying Equation (28) as follows

$\displaystyle N_{\rm ss}(r)=\begin{cases}N&,r\lesssim r_{\epsilon}\\\
N_{b}&,r_{\epsilon}\lesssim r\ll d,\end{cases}$ (41)

which can be now applied at arbitrarily small distances from the point-mass.
In what follows, it is worth bearing in mind that all field particles found
within the thermal critical radius $r\ll r_{\epsilon}$ have negative binding
energies. In other words, one should expect no interstellar visitors with
$E>0$ at radii $r\ll r_{\epsilon}$.

Adding the self-gravity of $m_{\star}$ to the equations of motion
predominantly affects the distribution of bound particles in the vicinity of
the point-mass. In particular, bottom-right panel shown that in collisional
models the number of bound particles within a volume element $V$ can now
exceed the number of random particles expected in a constant-density
background, such that $N_{\rm ss}\gtrsim N=n\,V$ at $r\lesssim r_{\epsilon}$.
This enhancement is due to the systematic increase in the number of particles
entering the volume $V$, as shown in Fig. 4. As a result, in collisional
models the steady-state number of bound particles lies systematically above
the value found in collisionless experiments and below a naive extrapolation
of Equation (28) below the critical thermal radius, such that $N\lesssim
N_{\rm ss}\lesssim N_{b}$ at $r\lesssim r_{\epsilon}$, with $N_{b}$ given by
Equation (11) (shown with black-dotted lines).

Figure 6: Left panels: Spatial locations of halo particles bound to a point-
mass $m_{\star}=3\times 10^{8}\,M_{\odot}$. Colours denote mean density
enhancements measured at the location of each individual particle ($\delta$).
For reference, the thermal critical radius ($r_{\epsilon}$) is marked with a
black circle. Middle-left panels: density enhancement profile of halo
particles, $\delta(r)$, as a function of distance from the central point-mass.
Theoretical curves (29) are shown with black-dashed lines. Black-dotted
particles denote the mean density of field particles ($\delta=1$). Middle-
right panels: One-dimensional velocity dispersion profiles of bound halo
particles as a function of distance from the point-mass. The theoretical
profile (30) is shown with a black-dashed line. Right panels: Mean phase-space
density of bound halo particles ($Q_{h}=n\,\delta/\sigma^{3}_{h}$) as a
function of distance. Black-dotted particles denote the phase-space density of
the field ($Q_{h}/Q_{g}=1$).

####  3.2.3 Profiles

Over time, tidally-trapped particles build up a‘halo’ around the point-mass
$m_{\star}$. Fig. 6 shows the spatial and kinematic distribution of particles
bound to the Keplerian potential $\Phi_{\star}$, together with the analytical
profiles derived in §2.4. Left panels plot the positions of bound particles on
the orbital plane of $m_{\star}$ measured from the intermediate-mass model
with $m_{\star}=3\times 10^{8}\,M_{\odot}$. Particles are colour-coded
according to the density of bound particles at their individual locations
($\delta$). Comparison between collisionless (top panel) and collisional
models (bottom panel) indicates that the differences introduced by the
gravitational attraction of $m_{\star}$ are mostly confined within the thermal
critical radius, $r_{\epsilon}$, which is marked with a solid-black circle for
reference. In particular, we find that collisional haloes exhibit much higher
densities within the sphere of influence of the point-mass, $r\lesssim
r_{\epsilon}$, than the collisionless counterparts.

This can be better seen in the middle-left panels, which show the density
enhancement of field particles with negative binding energies,
$\delta(r)\equiv n_{h}(r)/n$, in collisionless experiments. For ease of
reference, dotted-black lines mark the density of the local galactic
background ($\delta=1$). The first noteworthy result is that by plotting
distances in units of the thermal critical radius, $r_{\epsilon}$, the
profiles associated with different point-masses collapse into a single curve,
thus revealing a self-similar behaviour of tidally-trapped haloes. Comparison
with the theoretical expectation (black-dashed lines) shows that Equation (29)
is accurate at large distances $r\gtrsim r_{\epsilon}$, but systematically
overestimates the density of bound particles within the thermal critical
radius, $r\lesssim r_{\epsilon}$. In particular, in their central regions halo
profiles progressively depart from the power-law $\delta\sim r^{-3/2}$ and
become “cored”. The presence of a central core in the trapped halo is indeed
expected from Equation (41)

$\displaystyle\delta(r)=\frac{1}{4\pi\,n\,r^{2}}\frac{{\rm d}N_{\rm ss}}{{\rm
d}r}=\begin{cases}1&,r\lesssim r_{\epsilon}\\\
4/(3\sqrt{\pi})(Gm_{\star})^{3/2}\sigma^{-3}e^{-V_{\star}^{2}/(2\sigma^{2})}r^{-3/2}&,r\gtrsim
r_{\epsilon}.\end{cases}$ (42)

In contrast, bottom middle-left panel indicates that collisional models do not
follow a constant-density core at small distances. Instead, the power-law
‘spike’ $\delta\sim r^{-3/2}$ arises below the thermal critical radius, albeit
with a lower normalization than outside the critical radius. As shown in the
left panels of Fig. 5, this can be traced back to the enhanced number of field
particles crossing the volume $V$ in collisional models.

Top middle-right panels of Fig. 6 show the one-dimensional velocity dispersion
profiles of trapped particles in a steady-state. Although not shown here, we
find that the trapped halo has an isotropic velocity distribution at all
radii. As in previous plots, measuring distances in units of the thermal
critical radius ($r_{\epsilon}$) removes the point-mass dependence from the
profiles and yields a single, self-similar function, which follows very
closely a Keplerian profile (30), $\sigma(r)\sim r^{-1/2}$ (black-dashed
lines) at large distances $r\gtrsim r_{\epsilon}$. Again, the main differences
between collisionless and collisional models are restricted to distances below
the thermal critical radius. Namely, in collisionless experiments the velocity
dispersion profile departs from a Keplerian power-law, converging instead
towards the local velocity dispersion of the host galaxy, i.e.
$\sigma_{h}\approx\sigma$. In contrast, collisional models do not show a
convergent velocity dispersion at small distances, but rather a continuation
of the Keplerian profile towards arbitrarily small distances from the point-
mass.

We can now combine the density and velocity dispersion profiles measured above
in order to test whether the mean phase-space density of halo particles
($Q_{h}=n\,\delta/\sigma_{h}^{3}$) varies with radius. Recall that the
statistical theory outlined in §2 predicts a constant phase-space density
across small volume elements with a size $r\ll d$. The halo-to-field phase-
space density ratio $Q_{h}/Q_{g}$ given by Equation (31) solely depends on the
speed at which the point-mass travels across the Galaxy relative to the local
velocity dispersion. In our experiments, $V_{\star}/\sigma\simeq 0.88$, which
yields $Q_{h}/Q_{g}\simeq 2$. In contrast, right panel of Fig. 6 that the
phase-space density of bound halo particles varies with radius, and only
becomes constant at large and small distances from the point-mass. The value
predicted by Equation (31) (marked with a horizontal black-dashed line)
matches the empirical results at the thermal critical radius, $r\approx
r_{\epsilon}$. Below this distance scale, $Q_{h}$ approaches a constant value
that depends on whether or not the point-mass self-gravity is incorporated
into the equations of motion. In collisionless models, the halo phase-space
density of the halo converges towards the field value, $Q_{h}/Q_{g}\approx 1$
at $r\lesssim r_{\epsilon}$, whereas in the collisional counterparts the
phase-space density ratio falls below unity, $Q_{h}/Q_{g}\lesssim 1$, which
reflects the shorter dynamical lifetime of trapped particles in collisional
simulations (see right panels of Fig. 3).

It is also interesting to highlight the effect of the point-mass force on the
spatial and kinematic distribution of halo particles. Namely, while
collisionless simulations exhibit constant density profiles and flat velocity
dispersion profiles below the thermal critical radius, collisional models
follow power-law profiles $\delta\propto r^{-3/2}$ and $\sigma_{h}\propto
r^{-1/2}$ in regions dominated by the point-mass attraction. Remarkably, both
collisionless and collisional models exhibit a constant phase-space density
below the thermal critical radius, $Q_{h}\approx{\rm const}^{\prime}$ at $r\ll
r_{\epsilon}$.

Figure 7: Left panels: Semimajor axis distribution of particles trapped by a
point-mass, $m_{\star}/M_{g}=1.3\times 10^{-4}$, within a volume $V=4\pi
r^{3}/3$. For reference, volume sizes ($r$) are marked with vertical dotted
lines. The theoretical relation (53) is shown with a blacked-dashed line. Note
the sharp drop in the probability to find particles with a semi-major axis
$a\gtrsim r/2$, which correspond to orbits with apocentres larger than the
volume size, $r_{a}\geq 2a$. Middle panels: Eccentricity distribution of the
models shown in the left panel. The ‘thermal’ function (35) is shown with a
blacked-dashed line. Regardless of volume size, all models exhibit “super-
thermal” eccentricity distributions due to the presence of particles with
large apocentres $r_{a}\gtrsim r$ penetrating the volume $V$ on highly-
eccentric orbits. Right panels: Distribution of orbital pericentres. By
construction, no particle found within the volume $V$ has orbital pericentres
larger than the radius of the sphere, which translates into a truncation in
the distribution $p(r_{p})$ at $r_{p}=r$. Notice that at close distances from
the point-mass, $r_{p}\ll r$, the probability to detect an object within a
pericentre interval $r_{p},r_{p}+{\rm d}r_{p}$ becomes independent of $r_{p}$,
as expected from Equation (36).

####  3.2.4 Orbits

The final task of this Section is to inspect the orbits of particles trapped
in a Keplerian potential $\Phi(r)=-Gm_{\star}/r$. For simplicity, we compute
the orbital elements of individual particles using the relative position and
velocity measured at the time when the binding energy flips from a positive to
a negative value, and work under the assumption that the integrals of motion
are conserved quantities. With this simplification in mind, Fig. 7 plots the
semi-major axis (left panel) and eccentricity (right panel) distributions of
particles accreted by a point-mass $m_{\star}=3\times 10^{8}\,M_{\odot}$
within a volume element with different sizes (marked with vertical dotted
lines for reference). For a better understanding of the effect of the point-
mass self-gravity, we show collissionless and collisional modes in the upper
and bottom panels, respectively.

In all our models, the semi-major axis distribution exhibits at prominent peak
at $a=r/2$, which for radial orbits corresponds to apocentre distance that
equals the volume size, $r_{a}=2a=r$. We find that the numerical curves
approach the scale-free behaviour $p(a)\sim a^{1/2}$ predicted by (34) at
semi-major axes that are comparable or smaller than half of the volume size,
i.e. $a\ll r/2$. The probability to find particles with larger semi-major axes
drops significantly, but it is far from zero, which means that some particles
found at this instant of time within the volume $V$ move on orbits that take
them in and out of the region under observation, i.e their orbital aponcentres
are larger than the volume size, $r_{a}>r$.

These particles move on highly eccentric orbits, which leads to a “super-
thermal” eccentricity distribution in the right panel of Fig. 7. In
particular, there is a significant excess of halo particles with $e\gtrsim
0.7$ with respect to the thermal distribution (35), $p(e)=2e$. As expected,
the eccentricity distribution progressively thermalizes as the volume size
shrinks and the population of particles with large apocentres $r_{a}>r$
decreases.

Comparison between collisionless and collisional models highlights the effect
of incorporating the point-mass self-gravity into the equations of motion (2).
Notice first that the fraction of particles with apocentres $r_{a}=2a\gtrsim
r$ shows a systematic increase in collisional models with respect to the
collisionless counterparts. This is in line with the results plotted in Fig.
5, which show that the gravitational attraction of $m_{\star}$ boosts the
number of objects entering the volume $V$ from distant regions on very
eccentric orbits. As a result, the eccentricity distribution becomes less
thermal as the volume size shrinks and the gravitational attraction of the
point-mass becomes dominant, thus reversing the trend found in collisionless
simulations.

Right panel of Fig. 7 shows the distribution of orbital pericentres in these
models, $r_{p}=a(1-e)$. As expected from Equation (36), we find that the
probability to find trapped particles in the range $r_{p},r_{p}+{\rm
d}\,r_{p}$ becomes approximately constant for orbital pericentres much smaller
than the volume size, $r_{p}\ll r$. This result holds independently of whether
the attraction of the point-mass is included in the equations of motion, which
has important implications for the detection of trapped interstellar objects
in the solar system, as discussed in §4.

####  3.2.5 Galactocentric distance

The above experiments place point-masses on circular orbits at a fixed
galactocentric distance from the galactic potential, $R=0.065~{}R_{0}$, where
$R_{0}$ is the host scale radius (see §3.1). In this Section, we vary the
orbital radii of the point-mass in order to evaluate how the halo size changes
across the host galaxy. To characterize the halo size, we measure the steady-
state number of bound particles enclosed within tidal radius of the point-
mass, which is computed as (e.g. Renaud et al. 2011; Peñarrubia et al. 2016)

$\displaystyle
r_{t}(R)=\bigg{[}\frac{Gm_{\star}}{\gamma(R)\,\Omega^{2}(R)}\bigg{]}^{1/3}=\bigg{[}\frac{m_{\star}}{\gamma(R)\,M_{g}(<R)}\bigg{]}^{1/3}R,$
(43)

where $\gamma(R)=|{\rm d}\log n/{\rm d}\log R|$ is the power-law slope of the
host’s density profile, and $\Omega^{2}(R)=GM_{g}(<R)/R^{3}$ is the circular
frequency about a host galaxy with a mass profile
$M_{g}(<R)=4\pi\int_{0}^{R}{\rm d}r\,r^{2}\,n(r)$. §3 uses Dehnen spheres to
model the galactic potential. These objects have logarithmic density slopes
that approach $\gamma(R)\to\gamma$ in the limit $R\to 0$, and $\gamma(R)\to 4$
in the limit $R\to\infty$. Note that in a Keplerian potential the power-law
index and the enclosed mass are radially-independent quantities, $\gamma=3$
and $\Omega^{2}=GM_{g}/R^{3}$, such that Equation (43) reduces to the well-
known Jacobi radius $r_{t}=[m_{\star}/(3M_{g})]^{1/3}R$.

As discussed in Section 2, the statistical theory is only accurate insofar as
the local approximation can be applied, which demands the tidal radius to be
much smaller than the scale length associated with the radial density profile
of the host, i.e. $r_{t}\ll d$. It is straightforward to show that this
condition fails within small distances from the galaxy centre, $R\leq R_{d}$.
For a Dehnen (1993) galaxy model, the condition $r_{t}(R_{d})=d(R_{d})$ can be
solved analytically using (39) and (43), which yields

$\displaystyle
R_{d}=\bigg{(}\gamma^{2}\frac{m_{\star}}{M_{g}}\bigg{)}^{1/(3-\gamma)}R_{0},$
(44)

where $M_{g}=4\pi\,\rho_{0}\,R_{0}^{3}/(3-\gamma)$ is the total halo mass. For
low-mass point-masses, $m_{\star}\ll M_{g}$, Equation (44) returns $R_{d}\ll
R_{0}$. It is thus clear that the local approximation fails at small distances
from the centre of the host galaxy555Notice the special case of host galaxies
with a constant density profile $\gamma=0$. In these systems the problem is
not well defined due to centrally-divergent behaviour of the tidal radius,
$r_{t}\to\infty$ at $r\to 0$, see Renaud et al. (2011) for a detailed
discussion..

Fig. 8 plots the number of energetically-bound particles in steady-state for
point-masses moving on circular orbits at different galactocentric radii. For
reference, the numbers predicted by Equation (28) are shown with black-dashed
lines. The number of trapped objects declines with galactocentric radius,
which simply reflects the decreasing phase-space density of galactic particles
as the distance increases. As expected, we find a mismatch between the
theoretical curves and the numerical values at small galactocentric radii. In
particular, Equation (28) overpredicts the number of trapped halo particles at
$R\lesssim 3R_{d}$, where $R_{d}$ is given by Equation (44) and marked with
vertical arrows for reference. Note that since $R_{d}\to 0$ as $m_{\star}\to
0$, the range of galactocentric distances where the theory is accurate
increases as the perturber mass decreases.

Figure 8: Steady-state number of energetically-bound particles located within
the tidal radius of point-masses moving on circular orbits at different
galactocentric distances ($R$). Here, $R$ is measured in units of the host
scale-length $R_{0}$ (see §3.1). Black-dashed lines show the values predicted
by Equation (28). Note that the theoretical curves over-predict the number of
trapped particles at distances $R\lesssim 3R_{d}$, where $R_{d}$ corresponds
to the galactocentric distance at which the point-mass tidal radius $r_{t}$ is
equal to the local scale-length $d$ (values derived from Equation (44) are
marked with arrows).

##  4 Discussion: Trapped matter in the solar system

According to the statistical theory presented above, the solar system must be
embedded in an extended ‘halo” of interstellar objects temporarily trapped on
energetically-bound orbits around the Sun. This Section provides a rough
estimate of the steady-state number of interstellar objects and Dark Matter
(DM) particles in the halo, and how they are distributed across the solar
system. It is worth noting that these estimates contain a significant number
of uncertainties and must be therefore viewed as preliminary steps in an
ongoing effort to solve a difficult problem. For example, the statistical
analysis neglects the response of trapped halo particles to interactions with
planets as well as Galactic substructures in the vicinity of the Sun, which
could potentially affect the numerical values derived below (see Section 4.4).

For the estimates below, the Sun is placed at
${\textbf{{R}}}_{\odot}=(-8.17,0.0,0.02)$ kpc (Gravity Collaboration et al.
2019). At this location, the circular velocity is $V_{c}(R_{\odot})=237\,{\rm
km\,s^{-1}}$ (McMillan et al. 2017) and the Local Standard of Rest (LSR) is
${\textbf{{V}}}_{\rm LSR}=(11.1,12.2,7.3)\,{\rm km\,s^{-1}}$ (Schöenrich et
al. 2010). The Galactic potential is modelled as a disc, bulge and dark matter
halo with the following parameters. The disc is a Miyamoto-Nagai (1975) model
with a mass $M_{d}=6.6\times 10^{10}M_{\odot}$, and radial and vertical scale
lengths $a=8$ kpc and $b=0.3$ kpc. The MW bulge follows a Hernquist (1990)
profile with a mass $M_{b}=2.3\times 10^{10}M_{\odot}$ and a scale radius
$c=1.2$ kpc. The MW dark matter halo is modelled as a Navarro, Frenk & White
(1997) profile with a virial mass $M_{\rm vir}=1.2\times 10^{12}M_{\odot}$,
virial radius $r_{\rm vir}=222$ kpc and scale-radius $r_{s}=14.1\,{\rm kpc}$
(Bland-Hawthorn & Gerhard 2016). In this model, the velocity dispersion of the
DM halo at the solar radius is $\sigma_{\rm DM}\simeq 135\,{\rm km\,s^{-1}}$.
For the stellar velocity dispersion in the Solar vicinity we adopt the values
measured by Anguiano et al. (2020) using Gaia DR2 and APOGEE data, who find
$\sigma_{\rm thin}\simeq 28\,{\rm km\,s^{-1}}$, $\sigma_{\rm thick}\simeq
54\,{\rm km\,s^{-1}}$ and $\sigma_{\rm halo}\simeq 120\,{\rm km\,s^{-1}}$, for
the thin and thick discs and the stellar halo, respectively. For simplicity,
we will assume that interstellar particles move on isotropic trajectories.

Figure 9: Left panel: Density enhancement, $\delta(r)\equiv n_{h}(r)/n$, as a
function of distance from the Sun for baryonic ISOs and DM particles. Dotted
lines correspond to the analytical formula (29), whereas solid lines show the
empirical curves found in collisionless models, Equation (42). The Solar
radius is marked with a vertical-dotted line for reference. Right panel:
Velocity dispersion profile of ISOs and DM particles bound to the Solar
potential (black dashed lines). For reference, coloured-horizontal lines show
the velocity dispersion of each Milky Way component considered in the
analysis. Arrows mark the radii where $\delta=1$, Equation (21).

###  4.1 Halo composition

As the Sun moves through the Milky Way, it traps light bodies that populate
the interstellar background (comets, rocks, asteroids, DM particles, etc). The
frequency of the entrapment events strongly depends on the kinematics of these
objects in the vicinity of the solar system.

Left panel of Fig. 9 plots the density enhancement of halo objects
($\delta\equiv n_{h}/n$) as a function of distanced from the Sun. Lines are
colour-coded according to whether ISOs are tidally injected from the thin and
thick discs, and the stellar and Dark Matter (DM) haloes. Dotted and solid
lines correspond to Equation (29) and the empirically-corrected expression
(42), respectively. Recall that these estimates adopt $\alpha=1$, which is a
conservative choice at distances $r\gtrsim r_{\epsilon}$ (see numerical
experiments plotted in the right panel of Fig. 3). As expected, tidal trapping
is most efficient for bodies that belong to the thin disc (red lines), which
co-rotate with the Sun and have a low velocity dispersion. Interestingly, the
radius where the population of trapped ISOs dominates over the number of
(unbound) Milky Way interlopers roughly coincides with the radius of the
Earth’s orbit, $r_{\epsilon}^{\rm ISO}\approx 0.8\,{\rm AU}$. Therefore, the
majority of ISOs reaching the inner-most regions of the solar system will be
energetically bound to the Sun, and by definition members of the trapped halo.

In addition, the halo of trapped matter also contains a dark matter component.
The DM density enhancement is relatively low on account on the large velocity
dispersion of DM particles ($\sigma_{\rm DM}\simeq 135\,{\rm km\,s^{-1}}$) and
the fact that the Sun is speeding through the Galactic halo at $\sim 245\,{\rm
km\,s^{-1}}$. According to the estimates plotted in Fig. 9, the density of
bound DM particles exceeds the local DM background in the proximity of the
Solar surface, $r_{\epsilon}^{\rm DM}\approx 0.015\,{\rm AU}\simeq
3\,r_{\odot}$, where $r_{\odot}$ is the nominal Solar radius, suggesting that
the vast majority of DM particles cross the solar system on unbound
trajectories.

Right panel of Fig. 9 plots the velocity dispersion of particles trapped in
the halo (black-dashed line) as a function of distance from the Sun, Equation
(30). In our statistical model, all bound particles follow the same dispersion
profile regardless of their Galactic origin. In this sense, the process of
tidal trapping erases any dynamic ‘memory’ of their past orbits in the
Galactic potential. Comparison with the velocity dispersion of the Milky Way
components (plotted with horizontal lines for reference) shows that the halo
of trapped objects is much ‘colder’ than the Galactic environment at $r\gtrsim
1\,{\rm AU}$. At smaller distances, the halo velocity dispersion raises
towards the Sun as $\sigma_{h}\sim r^{-1/2}$, exceeding the background
velocity dispersion at distances comparable to the thermal critical radius,
$r\lesssim r_{\epsilon}$, where $r_{\epsilon}$ is given by Equation (21)
(marked with vertical arrows for ease of comparison). Recall that in
collisionless approximation, the radius at which the density enhancement
reaches unity coincides with the thermal critical radius, i.e.
$\delta(r_{\epsilon})=1$ (see §3.2.3).

###  4.2 Interstellar Objects

In our theory, the number of tidally-trapped ISOs is determined by the local
background density and the velocity dispersion of these objects in the solar
vicinity. The discoveries of the irregular body ’Oumuamua (Meech et al. 2017)
and the comet Borisov (Jewitt & Luu 2019), both travelling on a unbound
orbits, suggests that the local number density of ISOs is of the order of
$n^{\rm ISO}\sim 2\times 10^{15}\,{\rm pc}^{-3}$ (Do et al. 2018; Jewitt et
al. 2020). The local velocity distribution of ISOs is still unknown. Here, we
use the estimates plotted in Fig. 9, which suggest that the contribution from
the thick disc and stellar halo to the population of trapped ISOs can be
neglected at leading order, such that $\sigma^{\rm ISO}\approx\sigma_{\rm
thin}$.

Left panel of Fig. 10 shows the expected number of ISOs as a function of
distance to the Sun. For reference, we plot the cumulative number of visitors
expected from a constant-density background, $N^{\rm ISO}(r)=4\pi n^{\rm
ISO}r^{3}/3$, is plotted with a blue line, while the steady-state number of
energetically-bound ISOs (orange line) is estimated from Equation (41) as

$\displaystyle N_{\rm ss}^{\rm ISO}(r)\approx 4\times
10^{6}\,\bigg{(}\frac{m_{\star}}{M_{\odot}}\bigg{)}^{3/2}\bigg{(}\frac{n^{\rm
ISO}}{2\times 10^{15}\,{\rm pc}^{-3}}\bigg{)}\bigg{(}\frac{\sigma^{\rm
ISO}}{28\,{\rm km\,s^{-1}}}\bigg{)}^{-3}\bigg{(}\frac{r}{2\times 10^{4}\,{\rm
AU}}\bigg{)}^{3/2}$ (45)

at $r\gg r_{\epsilon}$. We stress that the estimate (45) becomes uncertain
beyond the region populated by Oort comets, $r\gtrsim 2\times 10^{4}\,{\rm
AU}$, which is poorly constrained observationally and may be subject to
Galactic perturbations (see §4.4).

According to (45), there are of the order of $N_{\rm ss}^{\rm ISO}(30\,{\rm
AU})\sim 226$ ‘Oumuamua-like objects trapped in the solar system at distances
comparable to the orbit of Neptune. Interestingly, Namouni & Morais (2020)
found several Centaurs and a few Trans Neptunian objects that may have been
captured from the interstellar background, which would go in line with our
theoretical predictions (but see Morbidelli et al. 2020).

The population of bound ISOs is dwarfed by the number of Milky Way interlopers
currently passing through the solar system, which exceeds $N^{\rm ISO}\sim
10^{13}$ at $r\sim 2\times 10^{4}\,{\rm AU}$. In contrast, the Oort Cloud is
thought to host some $\sim 10^{11}$ energetically-bound comets (Oort 1950;
Heisler 1990; Brasser & Morbidelli 2013; Feng & Bailer-Jones 2014; Rickman et
al. 2017). In agreement with Siraj & Loeb (2021), our estimates suggest that
unbound Milky Way interlopers may outnumber bound Oort cloud members by
roughly $\sim 2$ orders of magnitude. Note, however, that the size and spatial
distribution of the Oort Cloud is a long-standing matter of debate (for
reviews see Dones et al. 2004; Rickman 2014; Dones et al. 2015; Jewitt &
Seligman 2022). Direct observations in this region are difficult due to the
vast distances involved. Current estimates can only (indirectly) inferred by
comparing the predictions from planet formation models against the
(infrequent) visits of long-period comets to the inner solar system, or from
the modelling of the cratered surfaces of the Moon and terrestrial planets.
Existing estimates on the number of objects in the Oort cloud are model
dependent and still uncertain because the flux of long-period comets that
reach the inner solar system with pericentres $r_{p}<4\,{\rm AU}$ is only 2–3
per year, so an error of only 1 comet per year can significantly affect the
existing constraints (e.g. Moro-Martín 2019).

Crucially, existing estimates on the number of objects in the Oort Cloud do
not account for the presence of tidally-trapped ISOs in the solar system,
which are expected to contribute to the population of minor bodies with large
apocentres and negative binding energies. For example, Figs. 9 and 10 indicate
that most ISOs detected at $r\lesssim 1\,{\rm AU}$ will be bound to the Sun’s
potential and may not be easily distinguishable from genuine members of the
Oort Cloud. This issue has been recently discussed by Dehnen & Hands (2022),
who also point out that captured ISOs follow orbits akin to those of known
long-period comets, and that some of these comets could be of extra-solar
origin. These results open up interesting questions regarding the composition
of the outer solar system that may be worth investigating in a separate
contribution.

Figure 10: Left panel: Number of (baryonic) Interstellar Objects (ISOs) as a
function of distance from the Sun estimated from the number density and
velocity dispersion of ISOs in the Solar vicinity (see text). Orange line
plots the number of bound objects derived from Equation (41). For reference,
dotted-blue line shows the number of interlopers that one would expect from
the (unperturbed) field, $N=4\pi nr^{3}/3$. Horizontal and vertical dotted
lines mark $N=1$ and the canonical distance to the Oort Cloud, $2\times
10^{4}\,{\rm AU}$, respectively. Notice that most ISOs crossing the radius
$r\lesssim 1\,{\rm AU}$ are expected to be energetically bound to the sun.
Middle panel: Cumulative Dark Matter mass in the solar system as a function of
distance from the Sun. Right panel: Maxwellian velocity distribution of Dark
Matter particles expected at the Earth location. The presence of a bound DM
halo surrounding the solar system manifests as a narrow velocity peak centred
at the Earth-Sun relative velocity, $v_{\oplus}\simeq 29.7\,{\rm km\,s^{-1}}$,
which is marked with an orange arrow for reference (see text). In contrast,
unbound Milky Way dark matter particles crossing the solar system have a much
wider velocity distribution that peaks at the relative velocity between the
Sun and the Milky Way centre, $V_{\odot}\simeq 245\,{\rm km\,s^{-1}}$.

###  4.3 Dark Matter

The Sun can also trap DM particles from the Milky Way DM halo. For our
estimates, we adopt a background DM density of $\rho_{\rm
DM}=0.012\,M_{\odot}/\,{\rm pc}^{3}$ (e.g. Read 2014, and references therein),
and a velocity dispersion $\sigma_{\rm DM}=135\,{\rm km\,s^{-1}}$. Middle
panel of Fig. 10 plots the DM mass profile as a function of heliocentric
distance, with orange/blue lines denoting bound/unbound particles,
respectively. This plot highlights a number of interesting results. First,
notice first that close to the Solar surface, $r\lesssim r_{\epsilon}^{\rm
DM}\approx 0.015\,{\rm AU}$, all DM particles are energetically-bound to the
Sun, and are therefore members of the bound halo. This result may be relevant
in the context of capture of DM by the Sun, which is typically estimated by
adopting a velocity distribution that describes the (unbound) Milky Way DM
component at the solar location (for a review, see Nuñez-Castiñeyra et al.
2019), thus ignoring the effect of the Sun’s gravity on the trajectories of DM
particles as they approach the solar system.

At larger radii, $r\gg r_{\epsilon}^{\rm DM}$, the halo mass in bound DM
particles can be derived from (41) as

$\displaystyle M_{\rm h}^{\rm DM}(r)\approx 6\times
10^{-14}\,M_{\odot}\,\bigg{(}\frac{m_{\star}}{M_{\odot}}\bigg{)}^{3/2}\bigg{(}\frac{\rho_{\rm
DM}}{0.012\,M_{\odot}\,{\rm pc}^{-3}}\bigg{)}\bigg{(}\frac{\sigma_{\rm
DM}}{135\,{\rm km\,s^{-1}}}\bigg{)}^{-3}\bigg{(}\frac{r}{2\times 10^{4}\,{\rm
AU}}\bigg{)}^{3/2},$ (46)

which indicates that the amount of DM trapped in the solar system is small,
$M_{\rm h}^{\rm DM}(0.1\,{\rm pc})\sim 10^{-13}M_{\odot}$. This value sits
well below current upper limits derived from the motion of planets, which
constrain the halo mass within Saturn’s orbit below $<1.7\times
10^{-10}M_{\odot}$ (Pitjev & Pitjeva 2013). The presence of trapped DM
particles induces an extra gravitational attraction of $\sim 2\times
10^{-23}{\rm m}/{\rm s}^{2}$ at $r\sim 70\,{\rm AU}$, which cannot explain the
so-called Pioneer anomaly, a uniform acceleration of $\sim 8\times
10^{-10}{\rm m}/{\rm s}^{2}$ towards the Sun detected by spacecrafts Pioneer
10 and 11 at a distance between $20\,{\rm AU}$ and $70\,{\rm AU}$ (Turyshev &
Toth 2010 and references therein).

On Earth, the density of bound DM particles is subdominant with respect the
Galactic background. According to Equation (42) the density enhancement
induced by the bound DM halo is $\delta^{\rm DM}(r_{\oplus})\approx 1.8\times
10^{-3}$ at $r=r_{\oplus}=1\,{\rm AU}$. However, despite their relatively
small contribution to the DM density on Earth, trapped DM particles have a
relatively high phase-space density. This property is key for direct detection
attempts using elastic atomic-recoil experiments (Lewin & Smith 1996), as well
as to estimate capture rates of DM particles by the Earth, which are expected
to self-annihilate in the core and produce neutrinos that can be detected
(e.g. Gould 1987, 1988).

This issue is illustrated in the right panel of Fig. 10, which shows the
relative velocity distribution of DM particles expected at $r=1\,{\rm AU}$,
which under the Maxwellian approximation (7) can be written as

$\displaystyle p({\textbf{{v}}|r_{\oplus}})=p_{\rm
h}({\textbf{{v}}|r_{\oplus}})+p_{\rm
MW}({\textbf{{v}}|r_{\oplus}})=\frac{\delta(r_{\oplus})}{[2\pi\sigma_{h}^{2}(r_{\oplus})]^{3/2}}\exp\bigg{[}-\frac{(v-v_{\oplus})^{2}}{2\sigma_{h}^{2}(r_{\oplus})}\bigg{]}+\frac{1}{(2\pi\sigma_{\rm
DM}^{2})^{3/2}}\exp\bigg{[}-\frac{(v-V_{\odot})^{2}}{2\sigma_{\rm
DM}^{2}}\bigg{]},$ (47)

where the velocity vectors have been aligned with the Solar motion for ease of
comparison. In Equation (47), the left-hand term is due to the presence of
tidally-trapped DM particles on Earth, which manifest in a cold velocity peak
centred at the Earth-Sun velocity, $v_{\oplus}=29.7\,{\rm km\,s^{-1}}$ (marked
with an orange arrow for reference) with a dispersion
$\sigma_{h}(r_{\oplus})=(2/5)^{1/2}v_{\oplus}\approx 18.8\,{\rm km\,s^{-1}}$
given by (30). The right-hand term accounts for Milky Way DM particles
crossing the solar system on unbound trajectories, and exhibits a much wider
velocity distribution centred at the Galactocentric velocity,
$V_{\odot}=245\,{\rm km\,s^{-1}}$ (blue arrow) with a dispersion $\sigma_{\rm
DM}=135\,{\rm km\,s^{-1}}$. The relative amplitude between the two velocity
peaks is set by the mean phase-space density of the halo relative to the
background

$\frac{p_{\rm h}({\textbf{{v}}_{\oplus}|r_{\oplus}})}{p_{\rm
MW}({\textbf{{V}}_{\odot}|r_{\oplus}})}=\delta(r_{\oplus})\frac{\sigma^{3}_{\rm
DM}(R_{\odot})}{\sigma^{3}_{\rm
h}(r_{\oplus})}=\frac{Q_{h}(r_{\oplus})}{Q_{g}}.$

From Equation (31) we find $Q_{h}/Q_{g}\simeq 0.6$ at the Earth location.
Furthermore, since the Earth moves on a nearly circular orbit around the Sun,
experiments that aim to detect bound DM particles should expect no annual
modulation. Crucially, current exclusion limits on the DM particle mass &
cross section derived from atomic-recoil experiments & neutrino telescopes
neglect the cold velocity peak generated by tidally-trapped particles shown in
Fig. 10, and may therefore be more sensitive to higher particle masses than
previously thought.

###  4.4 Collisional effects

For simplicity, our analysis ignores the collisional element of the problem.
E.g., the above estimates do not account for collisions between infalling
stellar particles and the planets (or any other massive object) orbiting
around the Sun. It also ignores gravitational perturbations from Galactic
substructures in the vicinity of the solar system, either visible (e.g. stars
& gas clouds) or invisible (black holes, free-floating planets, dark matter
subhaloes, etc).

Galactic substructures generate random fluctuations of the local tidal field
that inject kinetic energy onto bound orbits (e.g. Chandrasekhar 1941;
Peñarrubia 2018). Over time, the cumulative effect of tidal fluctuations
causes tidal heating and a progressive erosion of the Oort cloud
(Chandrasekhar 1941; Hut & Tremaine 1985; Peñarrubia 2019; Torres et al. 2019;
Portegies Zwart 2021; Portegies Zwart et al. 2021). For example, Monte-Carlo
experiments carried by P21666Note that P21 experiments follow the evolution of
wide binaries in a clumpy potential. However, the results are applicable here
on account that the relative motion between binary stars is governed by the
same differential equations that determine the orbits of tracer particles
orbiting around a Keplerian potential.  show that tidal evaporation truncates
the semi-major axis distribution of halo particles, and that the truncation is
smaller than the tidal radius of the Keplerian potential (c.f. Fig. 8 in P21).
At present, it is unclear whether the solar system extends beyond the Oort
Cloud, which calls for caution when extrapolating the results of §4.2 and §4.3
at radii $r\gtrsim 2\times 10^{4}\,{\rm AU}$.

The role of planets in shaping the population of tidally-trapped objects is
more difficult to forecast, as their gravitational attraction induces
competing effects. On the one hand, interactions with planets may inject
energy onto the orbits of minor bodies orbiting around the Sun, thereby
shortening their dynamical lifetimes (e.g. Yabushita 1980; Napier et al.
2021b). Yet, in the case of tidally-trapped objects the number of hard
encounters is limited by the short time span that these objects spend in the
solar system, which is comparable to the duration of a fly-by encounter,
$t_{\rm surv}\sim T$ (see Fig. 3). This picture changes for semi-stable
captures with long survival times, $t_{\rm surv}\gg T$, which may experience a
significant amount of heating if their orbits bring them into the regions of
the solar system populated by planets.

On the other hand, 3-body collisions with planets provides an independent
mechanism for capturing interstellar visitors onto bound orbits, thus
contributing to the bound halo population. For example, Hands & Dehnen (2020)
simulate $\sim 4\times 10^{8}$ trajectories of interstellar objects that
approach a solar system in isolation (i.e. external forces are ignored), and
finds a steady-state population of up to $\sim 10^{4}$–$10^{5}$ ‘Oumuamua-
style objects, most of them located at $r\lesssim 100\,{\rm AU}$ from the Sun
(e.g. Fig. 5 of Dehnen et al. 2021).

The complex interplay between the above dynamical processes calls for follow-
up experiments that include collisional terms in the equations of motion (2).

##  5 Summary

This paper shows that Galactic tidal forces can inject ISOs and DM particles
onto temporarily-bound orbits around the Sun. The mechanics of tidal injection
are akin to slingshot manoeuvres, where the motion of the infalling object in
the Galactic potential combines with the gravity of the Sun to decelerate the
relative speed an infalling body to a degree where the binding energy becomes
negative. The same mechanism responsible for capturing interstellar particles
also leads to their later unbinding. As a result, tidally-trapped objects have
a limited life span in the solar system before being lost back to the Galactic
potential.

The presence of tidally-trapped objects creates an overdensity of interstellar
material surrounding the point-mass, or “halo”, which reaches a steady state
as the number of objects tidally injected from the field becomes comparable to
those being tidally stripped. It should be stressed that the halo does not
consist of the same bodies which are permanently captured, but rather of
temporarily-trapped bodies moving through the overdensity with negative
binding energies.

The stochastic model presented in this paper applies Smoluchowski (1916)’s
kinetic theory to derive the entrapment rate and the steady-state phase-space
distribution of tidally-trapped particles in a Keplerian potential,
$\Phi_{\star}(r)=-Gm_{\star}/r$. In a steady state, bound particles generate a
local density enhancement $\delta(r)\sim r^{-3/2}$ (also known as a ‘density
spike’) and follow an isotropic velocity dispersion profile that scales as
$\sigma_{h}(r)\sim r^{-1/2}$. Tests against live $N$-body models that allow
the Milky Way to respond to the gravitational attraction of a point-mass
$m_{\star}$ show that the analytical derivations are accurate on distance
range $r\gtrsim r_{\epsilon}$. Here, $r_{\epsilon}$ is the thermal critical
radius, Equation (21), which marks the distance from a point-mass at rest
where the Keplerian potential roughly equals the mean kinetic energy of field
particles, $|\Phi_{\star}(r_{\epsilon})|\approx 3\sigma^{2}/2$. Numerical
experiments show that the stochastic theory tends to over-estimate the density
of trapped particles below this scale (see §3).

It is worth noting that although tidal trapping does not require additional
mechanism(s) to remove kinetic energy from infalling interstellar material,
3-body interactions with planets, resonant capture and friction in the natal
gas cloud will also contribute to the population of bound interstellar
particles surrounding the solar system. Alas, these dynamical processes are
strongly intertwined, which adds significant complexity when making
predictions on the amount and composition of the interstellar matter currently
bound to the Sun.

Our theory predicts that other massive bodies travelling through the Milky Way
must also be embedded in tenuous envelopes of trapped interstellar particles.
Here we focus on the solar system, but similar arguments can be extended to
individual planets (e.g. Adler 2009), as well as Black Holes moving through
the Galaxy (e.g. Gondolo & Silk 1999; Boudaud et al. 2021). This possibility
will be exploring in separate contributions.

## Acknowledgements

This work has greatly benefited from the insightful comments of Malcolm
Fairbairn, Henrique Araujo, Michael Petersen and Raphael Errani. The author
would also like to thank the anonymous referee for very helpful suggestions.

## Data Availability

No data were generated for this study.

## References

  * Adams (2010) Adams F. C., 2010, ARA&A, 48, 47. doi:10.1146/annurev-astro-081309-130830
  * Adler (2009) Adler S. L., 2009, PhLB, 671, 203. doi:10.1016/j.physletb.2008.12.023
  * Anderson, Partenheimer, & Wiser (2020) Anderson N. B., Partenheimer A., Wiser T. D., 2020, arXiv, arXiv:2007.11016
  * Antonov (1961) Antonov, V. A. 1961, Soviet Ast., 4, 859
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn J., Gerhard O., 2016, ARA&A, 54, 529. doi:10.1146/annurev-astro-081915-023441
  * Boekholt, Rowan, & Kocsis (2022) Boekholt T. C. N., Rowan C., Kocsis B., 2022, arXiv, arXiv:2203.09646
  * Boudaud et al. (2021) Boudaud M., Lacroix T., Stref M., Lavalle J., Salati P., 2021, JCAP, 2021, 053. doi:10.1088/1475-7516/2021/08/053
  * Chandrasekhar (1941) Chandrasekhar, S. 1941, ApJ, 94, 511
  * Chandrasekhar (1943) Chandrasekhar, S. 1943, Reviews of Modern Physics, 15, 1
  * Charnoz & Morbidelli (2003) Charnoz S., Morbidelli A., 2003, Icar, 166, 141. doi:10.1016/S0019-1035(03)00213-6
  * Daniel, Heggie, & Varri (2017) Daniel K. J., Heggie D. C., Varri A. L., 2017, MNRAS, 468, 1453. doi:10.1093/mnras/stx571
  * Dehnen (1993) Dehnen, W. 1993, MNRAS, 265, 250
  * Dehnen & Hands (2022) Dehnen W., Hands T. O., 2022, MNRAS, 512, 4062. doi:10.1093/mnras/stab3670
  * Dehnen, Hands, & Schönrich (2021) Dehnen W., Hands T. O., Schönrich R., 2021, MNRAS.tmp. doi:10.1093/mnras/stab3666
  * Do, Tucker, & Tonry (2018) Do A., Tucker M. A., Tonry J., 2018, ApJL, 855, L10. doi:10.3847/2041-8213/aaae67
  * Dones et al. (2004) Dones L., Weissman P. R., Levison H. F., Duncan M. J., 2004, ASPC, 323, 371
  * Dones et al. (2015) Dones L., Brasser R., Kaib N., Rickman H., 2015, SSRv, 197, 191. doi:10.1007/s11214-015-0223-2
  * Drimmel & Poggio (2018) Drimmel R., Poggio E., 2018, RNAAS, 2, 210. doi:10.3847/2515-5172/aaef8b
  * Eddington (1916) Eddington, A. S. 1916, MNRAS, 76, 572
  * Edsjo & Peter (2010) Edsjo J., Peter A. H. G., 2010, arXiv, arXiv:1004.5258
  * Feng & Bailer-Jones (2014) Feng F., Bailer-Jones C. A. L., 2014, MNRAS, 442, 3653. doi:10.1093/mnras/stu1128
  * Gondolo & Silk (1999) Gondolo P., Silk J., 1999, PhRvL, 83, 1719. doi:10.1103/PhysRevLett.83.1719
  * Gould (1987) Gould A., 1987, ApJ, 321, 571. doi:10.1086/165653
  * Gould (1988) Gould A., 1988, ApJ, 328, 919. doi:10.1086/166347
  * Gravity Collaboration et al. (2019) Gravity Collaboration, Abuter R., Amorim A., Bauböck M., Berger J. P., Bonnet H., Brandner W., et al., 2019, A&A, 625, L10. doi:10.1051/0004-6361/201935656
  * Hands & Dehnen (2020) Hands T. O., Dehnen W., 2020, MNRAS, 493, L59. doi:10.1093/mnrasl/slz186
  * Heggie (1975) Heggie D. C., 1975, MNRAS, 173, 729
  * Heisler (1990) Heisler J., 1990, Icar, 88, 104. doi:10.1016/0019-1035(90)90180-H
  * Heisler & Tremaine (1986) Heisler J., Tremaine S., 1986, Icar, 65, 13. doi:10.1016/0019-1035(86)90060-6
  * Heppenheimer (1975) Heppenheimer T. A., 1975, Icar, 24, 172. doi:10.1016/0019-1035(75)90094-9
  * Heppenheimer & Porco (1977) Heppenheimer T. A., Porco C., 1977, Icar, 30, 385. doi:10.1016/0019-1035(77)90173-7
  * Higuchi & Ida (2017) Higuchi A., Ida S., 2017, AJ, 153, 155. doi:10.3847/1538-3881/aa5daa
  * Hills (1981) Hills J. G., 1981, AJ, 86, 1730. doi:10.1086/113058
  * Hunter (1967) Hunter R. B., 1967, MNRAS, 136, 245. doi:10.1093/mnras/136.3.245
  * Hut & Tremaine (1985) Hut P., Tremaine S., 1985, AJ, 90, 1548. doi:10.1086/113868
  * Iorio (2010) Iorio L., 2010, JCAP, 2010, 018. doi:10.1088/1475-7516/2010/05/018
  * Jeans (1928) Jeans J. H., 1928, asco.book
  * Jewitt et al. (2017) Jewitt D., Luu J., Rajagopal J., Kotulla R., Ridgway S., Liu W., Augusteijn T., 2017, ApJL, 850, L36. doi:10.3847/2041-8213/aa9b2f
  * Jewitt & Luu (2019) Jewitt D., Luu J., 2019, ApJL, 886, L29. doi:10.3847/2041-8213/ab530b
  * Jewitt et al. (2020) Jewitt D., Hui M.-T., Kim Y., Mutchler M., Weaver H., Agarwal J., 2020, ApJL, 888, L23. doi:10.3847/2041-8213/ab621b
  * Jewitt & Seligman (2022) Jewitt D., Seligman D. Z., 2022, arXiv, arXiv:2209.08182
  * Jílková et al. (2015) Jílková L., Portegies Zwart S., Pijloo T., Hammer M., 2015, MNRAS, 453, 3157. doi:10.1093/mnras/stv1803
  * Johansen et al. (2021) Johansen A., Ronnet T., Bizzarro M., Schiller M., Lambrechts M., Nordlund A., Lammer H., 2021, SciA, 7, eabc0444. doi:10.1126/sciadv.abc0444
  * Kokubo & Ida (1998) Kokubo E., Ida S., 1998, Icar, 131, 171. doi:10.1006/icar.1997.5840
  * Lehmann et al. (2021) Lehmann B. V., Ross O. G., Webber A., Profumo S., 2021, MNRAS, 505, 1017. doi:10.1093/mnras/stab1121
  * Lewin & Smith (1996) Lewin J. D., Smith P. F., 1996, APh, 6, 87. doi:10.1016/S0927-6505(96)00047-3
  * Li, Lai, & Rodet (2022) Li J., Lai D., Rodet L., 2022, ApJ, 934, 154. doi:10.3847/1538-4357/ac7c0d
  * Lingam & Loeb (2018) Lingam M., Loeb A., 2018, AJ, 156, 193. doi:10.3847/1538-3881/aae09a
  * Lundberg & Edsjö (2004) Lundberg J., Edsjö J., 2004, PhRvD, 69, 123505. doi:10.1103/PhysRevD.69.123505
  * Lynden-Bell (1999) Lynden-Bell, D. 1999, Physica A Statistical Mechanics and its Applications, 263, 293
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76. doi:10.1093/mnras/stw2759
  * Meech et al. (2017) Meech K. J., Weryk R., Micheli M., Kleyna J. T., Hainaut O. R., Jedicke R., Wainscoat R. J., et al., 2017, Natur, 552, 378. doi:10.1038/nature25020
  * Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
  * Morbidelli et al. (2020) Morbidelli A., Batygin K., Brasser R., Raymond S. N., 2020, MNRAS, 497, L46. doi:10.1093/mnrasl/slaa111
  * Moro-Martín (2019) Moro-Martín A., 2019, AJ, 157, 86. doi:10.3847/1538-3881/aafda6
  * Namouni & Morais (2020) Namouni F., Morais M. H. M., 2020, MNRAS, 494, 2191. doi:10.1093/mnras/staa712
  * Napier, Adams, & Batygin (2021) Napier K. J., Adams F. C., Batygin K., 2021a, PSJ, 2, 53. doi:10.3847/PSJ/abe76e/53
  * Napier, Adams, & Batygin (2021) Napier K. J., Adams F. C., Batygin K., 2021b, PSJ, 2, 217. doi:10.3847/PSJ/ac29bb
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Nuñez-Castiñeyra, Nezri, & Bertin (2019) Nuñez-Castiñeyra A., Nezri E., Bertin V., 2019, JCAP, 2019, 043. doi:10.1088/1475-7516/2019/12/043
  * Oort (1950) Oort J. H., 1950, BAN, 11, 91
  * Oncins et al. (2022) Oncins M., Miralda-Escudé J., Gutiérrez J. L., Gil-Pons P., 2022, arXiv, arXiv:2205.13003
  * Parker (2020) Parker R. J., 2020, RSOS, 7, 201271. doi:10.1098/rsos.201271
  * Peñarrubia (2015) Peñarrubia J., 2015, MNRAS, 451, 3537. doi:10.1093/mnras/stv1146
  * Peñarrubia et al. (2016) Peñarrubia J., Ludlow A. D., Chanamé J., Walker M. G., 2016, MNRAS, 461, L72. doi:10.1093/mnrasl/slw090
  * Peñarrubia (2018) Peñarrubia J., 2018, MNRAS, 474, 1482. doi:10.1093/mnras/stx2773
  * Peñarrubia (2019) Peñarrubia J., 2019, MNRAS, 484, 5409. doi:10.1093/mnras/stz338
  * Peñarrubia (2021) Peñarrubia J., 2021, MNRAS, 501, 3670. doi:10.1093/mnras/staa3700 (P21)
  * Peter (2009) Peter A. H. G., 2009a, PhRvD, 79, 103531. doi:10.1103/PhysRevD.79.103531
  * Peter (2009) Peter A. H. G., 2009b, PhRvD, 79, 103532. doi:10.1103/PhysRevD.79.103532
  * Peter (2009) Peter A. H. G., 2009c, PhRvD, 79, 103533. doi:10.1103/PhysRevD.79.103533
  * Petit & Henon (1986) Petit J.-M., Henon M., 1986, Icar, 66, 536. doi:10.1016/0019-1035(86)90089-8
  * Pfalzner (2013) Pfalzner S., 2013, A&A, 549, A82. doi:10.1051/0004-6361/201218792
  * Pitjev & Pitjeva (2013) Pitjev N. P., Pitjeva E. V., 2013, AstL, 39, 141. doi:10.1134/S1063773713020060
  * Portegies Zwart (2009) Portegies Zwart S. F., 2009, ApJL, 696, L13. doi:10.1088/0004-637X/696/1/L13
  * Portegies Zwart (2021) Portegies Zwart S., 2021, A&A, 647, A136. doi:10.1051/0004-6361/202038888
  * Portegies Zwart et al. (2021) Portegies Zwart S., Torres S., Cai M. X., Brown A. G. A., 2021, A&A, 652, A144. doi:10.1051/0004-6361/202040096
  * Press et al. (1992) Press, W. H., et al. 1992, Cambridge: University Press, 2nd ed.
  * Read (2014) Read J. I., 2014, JPhG, 41, 063101. doi:10.1088/0954-3899/41/6/063101
  * Renaud et al. (2011) Renaud, F., Gieles, M., & Boily, C. M. 2011, MNRAS, 418, 759
  * Rickman et al. (2014) Rickman H., Wiśniowski T., Wajer P., Gabryszewski R., Valsecchi G. B., 2014, A&A, 569, A47. doi:10.1051/0004-6361/201423966
  * Rickman et al. (2017) Rickman H., Wiśniowski T., Gabryszewski R., Wajer P., Wójcikowski K., Szutowicz S., Valsecchi G. B., et al., 2017, A&A, 598, A67. doi:10.1051/0004-6361/201629376
  * Schönrich, Binney, & Dehnen (2010) Schönrich R., Binney J., Dehnen W., 2010, MNRAS, 403, 1829. doi:10.1111/j.1365-2966.2010.16253.x
  * Singer (1968) Singer S. F., 1968, GeoJI, 15, 191. doi:10.1111/j.1365-246X.1968.tb05758.x
  * Siraj & Loeb (2019) Siraj A., Loeb A., 2019, ApJL, 872, L10. doi:10.3847/2041-8213/ab042a
  * Siraj & Loeb (2021) Siraj A., Loeb A., 2021, MNRAS, 507, L16. doi:10.1093/mnrasl/slab084
  * Smoluchowski (1916) Smoluchowski, M. V. 1916, Zeitschrift fur Physik, 17, 557
  * Suetsugu & Ohtsuki (2013) Suetsugu R., Ohtsuki K., 2013, MNRAS, 431, 1709. doi:10.1093/mnras/stt290
  * Torres et al. (2019) Torres S., Cai M. X., Brown A. G. A., Portegies Zwart S., 2019, A&A, 629, A139. doi:10.1051/0004-6361/201935330
  * Szebehely (1967) Szebehely V., 1967, torp.book
  * Turyshev & Toth (2010) Turyshev S. G., Toth V. T., 2010, LRR, 13, 4. doi:10.12942/lrr-2010-4
  * Valtonen & Innanen (1982) Valtonen M. J., Innanen K. A., 1982, ApJ, 255, 307. doi:10.1086/159830
  * Yabushita (1980) Yabushita S., 1980, A&A, 85, 77

##  Appendix A Density profile implied by a constant Distribution Function

In order to derive density and velocity dispersion profiles from a
distribution function, one needs to compute the velocity element

$\displaystyle{\rm d}^{3}v=2\pi v_{t}{\rm d}v_{t}{\rm d}v_{r}=2\pi\frac{L{\rm
d}L}{r^{2}}\frac{{\rm d}E}{v_{r}},$ (48)

where $v_{r}=\\{2[E+Gm_{\star}/r-L^{2}/(2r^{2})]\\}^{1/2}$ and $v_{t}=L/r$ are
the radial and tangential components of the velocity vector, respectively. For
particles moving on a Keplerian orbit the condition $v_{r}=0$ is met at
$r_{\rm a}=a(1+e)$ (apocentre) and $r_{\rm p}=a(1-e)$ (pericentre). Hence,

$\displaystyle
v_{r}^{2}=2\bigg{[}-\frac{Gm_{\star}}{2a}+\frac{Gm_{\star}}{r}-\frac{Gm_{\star}a(1-e^{2})}{2r^{2}}\bigg{]}=\frac{Gm_{\star}}{a}\bigg{(}\frac{a}{r}\bigg{)}^{2}\bigg{[}\frac{r}{a}-(1-e)\bigg{]}\bigg{[}(1+e)-\frac{r}{a}\bigg{]}.$
(49)

Inserting (49) into (48) and expressing energy and angular momentum as a
function of semimajor axis and eccentricity yields

$\displaystyle{\rm d}^{3}v=\frac{\pi(Gm_{\star})^{3/2}}{r}e{\rm
d}(-e)\frac{{\rm d}a}{a^{3/2}}\frac{1}{\\{[r/a-(1-e)][(1+e)-r/a]\\}^{1/2}}.$
(50)

The number density of halo particles is found by integrating (33) over (50)

$\displaystyle n_{h}(r)$ $\displaystyle=\int{\rm d}^{3}v\,f$ (51)
$\displaystyle=\pi
f_{0}\frac{(Gm_{\star})^{3/2}}{r}\int_{0}^{\infty}\frac{{\rm
d}a}{a^{3/2}}\int_{0}^{1}{\rm d}e\frac{e}{\\{[r/a-(1-e)][(1+e)-r/a]\\}^{1/2}}$
$\displaystyle=\pi
f_{0}\frac{(Gm_{\star})^{3/2}}{r}\int_{0}^{\infty}\frac{{\rm
d}a}{a^{3/2}}\bigg{\\{}\bigg{[}\frac{(2a-r)r}{a^{2}}\bigg{]}^{1/2}-\bigg{[}\frac{(2a-r)r-a^{2}}{a^{2}}\bigg{]}^{1/2}\bigg{\\}}$
$\displaystyle=\frac{4\pi\sqrt{2}}{3}f_{0}\bigg{(}\frac{Gm_{\star}}{r}\bigg{)}^{3/2},$

which recovers the profile (29) for

$\displaystyle
f_{0}=\frac{1}{2\sqrt{2\pi}}e^{-V_{\star}^{2}/(2\sigma^{2})}\,Q_{g},$ (52)

with $Q_{g}=n/\sigma^{3}$ being the mean phase-space density of the
interstellar background.

It is worth noting that the power-law profile (51) arises from particle
ensembles orbiting in a Keplerian potential with constant phase-space density.
This profile has been studied by other authors in several contexts (e.g.
Gondolo & Silk 1999; Oncins et al. 2022)

##  Appendix B Wide binaries

It is straightforward to show that the distribution function (33) also
describes the orbits of wide binaries formed in an homogeneous sea of
particles with Maxwellian velocities. Expressing (34) as a function of orbital
energy via the transformation $p(E)\,{\rm d}E=p(a)\,{\rm d}a$ with
$a=Gm_{\star}/(-2E)$ leads to an energy distribution

$\displaystyle p(E)\,{\rm d}E\sim\frac{{\rm d}E}{(-E)^{5/2}}.$ (53)

Hence, re-writing (33) in terms of energy yields

$\displaystyle p(E,e)\,{\rm d}E{\rm d}e\sim\frac{1}{(-E)^{5/2}}\,(2e)\,{\rm
d}E\,{\rm d}e,$ (54)

which matches the distribution of integrals found by Jeans (1928) (see also
Heggie 1975)

$\displaystyle p_{\rm pairs}(E,e){\rm d}E{\rm
d}e\sim\frac{\exp(-E/\mathcal{T})}{(-E)^{5/2}}(2e)\,{\rm d}E\,{\rm d}e,$ (55)

in the limit where stellar pairs move on orbits with small relative velocities
$v<v_{e}\ll\sigma$, and low binding energies,
$|E|\ll\mathcal{T}=3\sigma^{2}/2$, such that $\exp(-E/\mathcal{T})\approx 1$.

##  Appendix C Numerical integration

This Section briefly explains how the coupled differential equations (1) and
(2) are solved in this work. The main difficulty lies in the vastly different
dynamical ranges of particles bound to the point-mass and those moving the in
Galactic potential. For illustration, consider the models shown in Section 3,
which resolve distances below the critical radius of the point-mass,
$r_{0}=2Gm_{\star}/\sigma^{2}$. At this distance, the circular frequency is
$\Omega_{\star}(r_{0})=(Gm_{\star}/r_{0}^{3})^{1/2}=2^{-3/2}\sigma^{3}/(Gm_{\star})$.
This number can be directly compared with the circular frequency of the point-
mass around the host galaxy centre,
$\Omega_{g}(R_{\star})=[GM_{g}(<R_{\star})/R^{3}_{\star}]^{1/2}$. Estimating
the local velocity dispersion as $\sigma^{2}\approx
V_{c}^{2}/2=GM_{g}(<R_{\star})/(2R_{\star})$ yields a frequency ratio
$\Omega_{g}/\Omega_{\star}\sim 8m_{\star}/M_{g}(<R_{\star})$. Given that the
numerical models in Section 3 have mass ratios $m_{\star}/M_{g}\ll 1$, it
seems crucial to solve equations (1) and (2) taking into account the disparate
dynamical range between particles moving in the host galaxy and those
temporarily bound to the Keplerian potential.

To this aim, we apply a leap-frog scheme with a varying time-step. In