# The Population of Interstellar Objects Detectable with the LSST and
Accessible for In Situ Rendezvous with Various Mission Designs

Devin J. Hoover Dept. of Astronomy and Astrophysics, University of Chicago,
Chicago, IL 60637  Darryl Z. Seligman Dept. of the Geophysical Sciences,
University of Chicago, Chicago, IL 60637  Matthew J. Payne Harvard-Smithsonian
Center for Astrophysics, 60 Garden St., MS 51, Cambridge, MA 02138, USA

###### Abstract

The recently discovered population of interstellar objects presents us with
the opportunity to characterize material from extrasolar planetary and stellar
systems up close. The forthcoming Vera C. Rubin Observatory Legacy Survey of
Space and Time (LSST) will provide an unprecedented increase in sensitivity to
these objects compared to the capabilities of currently operational
observational facilities. We generate a synthetic population of ‘Oumuamua-like
objects drawn from their galactic kinematics, and identify the distribution of
impact parameters, eccentricities, hyperbolic velocities and sky locations of
objects detectable with the LSST, assuming no cometary activity. This
population is characterized by a clustering of trajectories in the direction
of the solar apex and anti-apex, centered at orbital inclinations of $\sim
90^{\circ}$. We identify the ecliptic or solar apex as the optimal sky
locations to search for future interstellar objects as a function of survey
limiting magnitude. Moreover, we identify the trajectories of detectable
objects that will be reachable for in-situ rendezvous with a dedicated mission
with the capabilities of the forthcoming Comet Interceptor or proposed Bridge
concept. By scaling our fractional population statistics with the inferred
spatial number density, we estimate that the LSST will detect of order $\sim
15$ interstellar objects over the course of its $\sim 10$ year observational
campaign. Furthermore, we find that there should be $\sim 1-3$ and $\sim
0.0007-0.001$ reachable targets for missions with propulsion capabilities
comparable to Bridge and Comet Interceptor, respectively. These numbers are
lower limits, and will be readily updateable when the number density and size
frequency distribution of interstellar objects are better constrained.

Interstellar Objects

##  1 Introduction

The recent discoveries of the first two interstellar objects (ISOs), 1I/2017
U1 (‘Oumuamua) and C/2019 2I (Borisov), imply that a galactic-wide population
of roughly $10^{26}$ similar bodies exist, with spatial number densities of
order $n_{o}\sim 1-2\times 10^{-1}\,$au-3 (Trilling et al., 2017; Laughlin &
Batygin, 2017; Do et al., 2018; Moro-Martín, 2019, 2018; Levine et al., 2021).
These number densities have been estimated by incorporating the detailed
detection capabilities of current observational facilities. Forthcoming
surveys, such as the Vera C. Rubin Observatory Legacy Survey of Space and Time
(LSST), will provide unprecedented completeness for populations of minor
bodies in the Solar System (Jones et al., 2009; Ivezić et al., 2019). It has
been estimated that the LSST will detect between $1-10$ interstellar objects
every year (Engelhardt et al., 2017; Cook et al., 2016; Trilling et al., 2017;
Seligman & Laughlin, 2018). It is important to note that the detection
frequency depends sensitively on the spatial distribution, size-frequency
distribution and levels of cometary activity of these objects. The future
detections and characterizations of these objects will yield insights into the
physical mechanisms at work in extra-solar planetary systems and star
formation regions throughout the galaxy.

The first interstellar object, ‘Oumuamua, defied all expectations for an
interstellar comet. It exhibited a distinct lack of dust and typical cometary
volatiles (Meech et al., 2017; Jewitt et al., 2017; Trilling et al., 2018), an
extreme $6:6:1$ aspect ratio (Knight et al., 2017; Drahus et al., 2018; Bolin
et al., 2018; McNeill et al., 2018; Belton et al., 2018; Mashchenko, 2019), a
reddened reflection spectrum (Masiero, 2017; Fitzsimmons et al., 2018;
Bannister et al., 2017; Ye et al., 2017), and a non-gravitational acceleration
(Micheli et al., 2018). Curiously, it also exhibited a surprisingly low
incoming velocity with respect to the Local Standard of Rest (LSR) in the
direction of the solar apex (Mamajek, 2017), implying a young $<40$Myr age
(Almeida-Fernandes & Rocha-Pinto, 2018; Hallatt & Wiegert, 2020; Hsieh et al.,
2021).

This peculiar combination of unique physical properties has led to a variety
of theories describing the provenance of ‘Oumuamua. The non-gravitational
acceleration has been postulated to be driven by radiation pressure (Micheli
et al., 2018), which could be explained if ‘Oumuamua was an ultra low-density
fractal aggregate (Moro-Martín, 2019; Luu et al., 2020; Sekanina, 2019), or a
millimeter-thin artificial light sail (Bialy & Loeb, 2018). If the non-
gravitational acceleration was powered by cometary outgassing (Micheli et al.,
2018; Seligman et al., 2019), it has been demonstrated that the energetics are
consistent with the sublimation of H2 (Füglistaler & Pfenniger, 2015; Seligman
& Laughlin, 2020; Levine & Laughlin, 2021), N2 (Jackson & Desch, 2021; Desch &
Jackson, 2021), or CO (Seligman et al., 2021). Recently, Grude Flekkøy &
Brodin (2022) calculated infrared and optical spectral signatures that would
differentiate between the various proposed scenarios.

In stark contrast to ‘Oumuamua, the interstellar comet Borisov displayed
physical characteristics that are broadly consistent with the census of
cometary bodies in the Solar System. It had dusty cometary activity (Jewitt &
Luu, 2019; Bolin et al., 2020b; Fitzsimmons et al., 2019; Ye et al., 2020;
McKay et al., 2020; Guzik et al., 2020; Hui et al., 2020; Kim et al., 2020;
Cremonese et al., 2020; Yang et al., 2021), and early observations confirmed
the presence of carbon and nitrogen based species in the outflow (Opitom et
al., 2019; Kareta et al., 2020; Lin et al., 2020; Bannister et al., 2020; Xing
et al., 2020; Aravind et al., 2021). It was particularly enriched in CO
(Bodewits et al., 2020; Cordiner et al., 2020), which could be explained if it
formed in the outer regions of a protoplanetary disk (Price et al., 2021). It
exhibited a brightening event (Drahus et al., 2020) and subsequent breakup in
the spring of 2020 (Jewitt et al., 2020a, b; Bolin et al., 2020a; Zhang et
al., 2020), which was likely due to seasonal effects (Kim et al., 2020).
Recently, Guzik & Drahus (2021) reported a spectroscopic detection of atomic
nickel vapor in the coma. Bagnulo et al. (2021) presented polarimetric
observations that demonstrated that the dust coma had abnormally high
polarization compared to Solar System comets. Astrometric analyses of the non-
gravitational acceleration are consistent with observed levels of outgassing
(Hui et al., 2020; de la Fuente Marcos & de la Fuente Marcos, 2020; Manzini et
al., 2020).

Future detections and subsequent characterization of interstellar objects will
be a welcome development. They should contextualize and/or explain the
discrepancy in the observed characteristics of the first two discovered
members of the population. In addition to ground and space-based observations
of these future objects, in situ measurements would provide an opportunity to
characterize their composition and physical properties in detail (Seligman &
Laughlin, 2018). ESA’s Comet Interceptor mission (Jones & ESA Comet
Interceptor Team, 2019; Pau Sánchez et al., 2021) or the NASA concept study
Bridge (Moore et al., 2021) may be well-positioned to provide these
observations.

In this paper, we investigate the population of interstellar objects that will
be detected by the LSST, with a focus on objects that could serve as realistic
targets for interception missions. In §2, we describe the galactic
distribution of interstellar objects from which we draw our simulated targets.
In §3, we describe the numerical simulations that we perform in order to
evaluate the trajectories and distances from the Earth that each of these
objects attain. In §4 and §5, we present the distribution of trajectories and
orbital elements for objects that are (i) detectable by the LSST and (ii)
realistic targets for an interception mission. We also provide an up-to-date
estimation for the number of objects that will be detected for both of these
populations. In §6, we summarize our calculations and conclude.

##  2 The Galactic Population of Interstellar Objects

In this section, we describe the methodology with which we simulate a galactic
population of interstellar objects. We assume that the population is
dynamically relaxed into the same kinematic distribution as local stars via
scattering events off of giant molecular clouds and dark matter substructure
(Seligman & Laughlin, 2018). Since this population consists of $\sim 10^{26}$
objects, it represents an excellent realization of the fine-grain assumption
in the collisionless Boltzmann equation.

We assume that the interstellar objects within the solar neighborhood follow a
Schwarzschild velocity distribution with respect to the LSR, with velocity
dispersions $\sigma_{R}$, $\sigma_{\phi}$, $\sigma_{z}$. The right-handed
galactic coordinate system is oriented by the vertex deviation, $v_{d}$, for
the primary eigenvector of the ellipsoidal distribution function. The spatial
number density of ISOs, $n({\bf\Delta v})$, in a specific volume of phase
space ${\bf\Delta v}=(\Delta v_{1},\Delta v_{2},\Delta v_{3})$, is given by,

$n({\bf\Delta
v})=\frac{n}{8}\,\bigg{(}\,\prod_{i=1}^{3}\textrm{erf}\big{(}\frac{v_{i}}{\sqrt{2}\sigma_{i}}\big{)}|_{v_{i}-\Delta
v_{i}}^{v_{i}+\Delta v_{i}}\bigg{)}\,.$ (1)

Note that in this equation the numbered subscript $i$ is used to indicate the
orthogonal components of $\sigma$ in the order previously indicated.

We adopt the values of $\sigma_{R}$, $\sigma_{\phi}$, $\sigma_{z}$, and
$v_{d}$ for G stars (Binney & Merrifield, 1998; Seligman & Laughlin, 2018). It
is important to note that different stellar populations correspond to
different typical ages and velocity dispersions, as summarized in Table 1. The
typical ages of interstellar objects is somewhat unconstrained. We verified
that the results of our simulations are not sensitively dependant on the
spectral type that we use. It would be informative to perform the calculations
presented in §4 and §5 for spectral classes representative of older and
younger populations, but this is outside of the scope of this paper.

Gravitational focusing preferentially increases the flux of interstellar
objects with low heliocentric velocities by a factor $\xi$, which we
approximate as,

$\xi\approx 1+(\frac{42km/s}{v_{\infty}+\beta})^{2}\,.$ (2)

Here, $v_{\infty}$ is the hyperbolic excess velocity with respect to the Sun
and $\beta$ is an artificial smoothing factor that we set to $\beta=$2 km/s to
avoid the hyperbolic case where $v_{\infty}=0$ (Seligman & Laughlin, 2018). We
compute the probability that an interstellar object enters the Solar System
with a given velocity using Equations 1 and 2.

We construct a three-dimensional galactic velocity grid by generating 100
evenly-spaced values for each component of heliocentric velocity between -100
and 100 km/s. We calculate the probability that an interstellar object will
enter the Solar System for each grid cell and corresponding velocity vector.

An incoming interstellar object will have a velocity at $r=\infty$ drawn from
this distribution, and will approach the Sun at a given impact parameter, $b$.
If we assume that the impact parameter is $b=0$, then as the object enters the
vicinity of the Sun, its position vector is parallel and opposite to its
velocity vector. Therefore, every velocity vector isomorphically maps to a
unit vector pointing away from the Sun towards a location on the unit sphere
(which we arbitrarily define as the 1 au unit sphere). We convert each
incoming heliocentric velocity to a corresponding declination and right
ascension. We show the resulting sky position probability distributions in
Figure 1 for kinematics representative of three different stellar populations,
described in Table 1. The majority of interstellar objects enter the Solar
System in the direction of the solar apex, regardless of the assumed stellar
population. This is because the apex is the direction with which the Sun is
moving with respect to the LSR.

Figure 1: Sky map illustrating the probability that a future interstellar object will enter the Solar System in a given direction, assuming $b=0$. The probabilities use the Schwarzschild distributions for M dwarfs (upper panel), G dwarfs (middle panel), and O/B giants (lower panel). Warmer colors indicate a higher probability. In all panels, the initial position of ‘Oumuamua is indicated with a red star. The majority of interstellar objects enter the Solar System from the direction of the solar apex. Stellar Type | $\sigma_{R}$ | $\sigma_{\phi}$ | $\sigma_{Z}$ | $v_{d}(^{\circ})$
---|---|---|---|---
M | 31 | 23 | 16 | 7
G | 26 | 18 | 15 | 12
O/B | 12 | 11 | 9 | 36
Table 1: Parameters in the Schwarzschild distributions for different stellar
types. All velocities are expressed in kilometers per second. Since different
stellar populations correspond to different velocity distributions, the
initial locations of our simulated objects depend on the distribution from
which we draw samples.

To generate samples from our grid, we define a multivariate spline that maps
velocity vectors to their associated probabilities. For each sample, we draw a
probability from a uniform distribution with bounds marked by the minimum and
maximum values of probabilities associated with each point on the grid. Then,
we use differential evolution to find the velocity vector that minimizes the
difference between the spline interpolation and probability. We repeat this
process $10^{6}$ times, generating a sample of 1,000,000 objects. We verified
that the distribution of our simulated objects was converged for this sample
size.

Figure 2: A schematic diagram showing the geometry of an incoming interstellar
object. Each object has an initial velocity vector drawn from the galactic
distribution, which maps to a location on the unit sphere in the heliocentric
reference frame. We then place the interstellar object at a position given by
an impact parameter, $b$, that we draw.

Using the distribution of incoming velocity trajectories, we construct a
realistic distribution of incoming interstellar object trajectories. We denote
our initial velocity vector as $\vec{v}_{\infty}$, where the magnitude,
$|\vec{v}_{\infty}|$, is the hyperbolic excess velocity. For the simplest case
where $b=0$, $\vec{v}_{\infty}$ is oriented in the direction of the Sun. The
incoming position vector, $\vec{R}$, is given by,

$\vec{R}=-k\vec{v}_{\infty}\,,$ (3)

for some constant, $k\in\mathbb{R}$. In reality, each incoming interstellar
object will have a nonzero impact parameter. For every initial velocity vector
in our synthetic population, we draw $b$ from a parabolic probability
distribution, $p(b)\sim b^{2}$ with lower and upper bounds at 0 and 5 au,
respectively. The parabolic distribution is a product of the cylindrical
radial cross section of the Solar System. We also draw an angle, $\psi$, from
a uniform distribution from $(0,2\pi)$ to represent the direction of $b$ in
the plane normal to the velocity vector. In other words, $\psi$ represents a
rotation within the tangent bundle of the 5 au sphere. We require that each
interstellar object be placed at an initial position exactly $R=5$ au away
from the Sun. Since we assume that the synthetic population consists of
‘Oumuamua-like objects without bright cometary activity, it is unlikely that
they will be detected by the LSST at large distances. For the purposes of this
paper, it is unnecessary to simulate any objects outside 5au. Although outside
the scope of this paper, it would be straightforward and worthwhile to create
a synthetic population of interstellar comets that included the effects of
cometary brightening, similar to the analysis presented in Cook et al. (2016).
In Figure 2, we show a schematic diagram of the geometry. We can substitute
$\vec{R}_{sky}$ for $\vec{R}$ in Equation 3 and solve for $k$ using,

$\vec{R}_{sky}=-\,\bigg{(}\sqrt{R^{2}-b^{2}}\bigg{)}\,\,\hat{v}_{\infty}\,,$
(4)

where $\hat{v}_{\infty}$ is the unit vector pointing in the direction of
$\vec{v}_{\infty}$. We define $\vec{b}$ as the vector with magnitude $b$ that
lies in the tangent plane rotated by $\psi$ and has the property
$\vec{b}\cdot\vec{v}_{\infty}=0$. We can write $\vec{b}$ as

$\vec{b}=\begin{bmatrix}b\cos{\phi}\cos{\psi}\\\ b\cos{\phi}\sin{\psi}\\\
b\sin{\phi}\\\ \end{bmatrix}\,,$ (5)

where the angle $\phi$ is given by,

$\tan\phi=-\,\bigg{(}\frac{\cos{\psi}\,v_{\infty,x}+\sin{\psi}\,v_{\infty,y}}{v_{\infty,z}}\bigg{)}\,.$
(6)

Here the subscripts on $v_{\infty}$ indicate the component of the vector. This
allows us to construct $\vec{b}$ for each initial velocity vector. The initial
position vector, $\vec{R},$ is given by

$\vec{R}=\vec{R}_{sky}+\vec{b}\,.$ (7)

We place each interstellar object at an initial distance of 5 au away from the
Sun in the direction of $\hat{R}$. In order to account for the Solar
acceleration, we multiply each velocity vector by a factor
$\sqrt{v_{\infty}^{2}+v_{esc}^{2}}/v_{\infty}$, where $v_{esc}$ is the solar
escape velocity at 5 au.

##  3 Numerical Simulations

For a given interstellar object initial condition, we use the N-body code
REBOUND (Rein & Liu, 2012) to integrate the trajectories of the Earth and the
interstellar object for a 3-year time period. This timespan is sufficient to
capture the interstellar object’s trajectory through the inner solar system,
which we define as the 5 au sphere centered at the Sun. By comparison, an
object on a parabolic orbit with perihelion at 1 au crosses the 5 au sphere in
$\sim 2$ years. Given the variability of the hyperbolic excess velocities
within our sample, a higher integration time is needed to accommodate ISOs
with low initial velocities and higher perihelia. For computational
efficiency, we neglect the contributions from the other planets in the Solar
System. In reality, the giant planets will provide gravitational perturbations
to a small subset of these interstellar objects. In the extreme case, Siraj &
Loeb (2019) investigated the efficiency with which Jupiter captures
interstellar objects. However, it has been shown that there is no evidence for
bound interstellar objects currently in the Solar System (Morbidelli et al.,
2020). Hands & Dehnen (2020) demonstrated that the volume capture rate of
interstellar objects was 0.051 ${\rm au^{3}\,yr^{-1}}$, of which only
$0.033\%$ are within 6 au at any time. The initial condition that we use could
be modified to investigate the possibility of capture for objects that start
at further distances, as studied by Napier et al. (2021), which is outside of
the scope of this paper. In any case, the majority of the interstellar objects
in our simulated population ($>99\%$) will be unaffected by perturbations from
the giant planets, which have small interaction cross sections in comparison
to the 5 au sphere.

An object with absolute magnitude, $H$, has an apparent magnitude, $m$, which
may be calculated using,

$m=H+2.5\log_{10}{\bigg{(}\frac{d_{BS}^{2}d_{BO}^{2}}{p(\theta)d_{OS}^{4}}\bigg{)}}\,.$
(8)

The parameters $d_{BO}$, $d_{OS}$, and $d_{BS}$ represent the distances
between the body and observer, Sun and observer, and body and Sun
respectively. $\theta$, the phase angle, and $p(\theta)$, the phase integral,
are defined by

$\cos{\theta}=\frac{d_{B0}^{2}+d_{BS}^{2}-d_{OS}^{2}}{2d_{BO}d_{BS}}\,,$ (9)

and

$p(\theta)=\frac{2}{3}((1-\frac{\theta}{\pi})\cos{\theta}+\frac{1}{\pi}\sin{\theta})\,.$
(10)

Since only two interstellar objects have been detected, the absolute magnitude
and size-frequency distribution of the population is unconstrained. Therefore,
we assume that all interstellar objects in our sample have the same value of
$H$ as ‘Oumuamua, $H=22.4$. In reality, the absolute magnitude will depend on
the size of the body. Moreover, if an interstellar object exhibits cometary
activity like Borisov, then it will be much brighter and the absolute
magnitude should change as a function of heliocentric distances. For the
purposes of this paper, we conservatively evaluate the detection capabilities
for the LSST for ‘Oumuamua-like objects. However, these distributions can be
updated to account for differences in intrinsic brightness, once the size-
frequency distribution of the population is better constrained.

We initialize the position of the Earth in our simulations from a distribution
of 24 start dates evenly distributed between January 1st and December 31st of
2021, although the year is arbitrary. For each simulated trajectory, we
randomly select one of these dates, which defines the initial condition of the
Earth. At 100 evenly-spaced times within the 3 year integration period, we
evaluate the position, velocity, apparent magnitude as seen from Earth,
declination/right ascension, altitude/azimuth, orbital elements, as well as
other parameters, for each interstellar object. Note that we implemented a
cubic interpolation to reconstruct aspects of the trajectories for some of the
results presented in §4. We perform these calculations for a synthetic
population of $\sim 10^{6}$ interstellar objects. We verified that the shape
of the distributions presented in the following two sections and fractional
percentages presented in Table 2 were numerically converged.

Figure 3: The initial positions of interstellar objects that are detectable by
the LSST (i.e. their positions when they enter the 5 au sphere prior to
detection). We list our detectability criteria in Section 4. While the
majority of simulated interstellar objects enter the inner Solar System in the
direction of the solar apex, a significant yet smaller fraction enter from the
anti-apex. We include the ecliptic, the LSST’s nominal field-of-view, and the
positions of ‘Oumuamua and Borisov at the moment that they entered the Solar
System for reference. While the LSST’s observational strategy is complex, with
the Northern Ecliptic Spur extending to +$30^{\circ}$ declination only for
some right ascensions (Bianco et al., 2021), the Wide-Fast-Deep, which
includes all right ascensions, extends to +$12^{\circ}$ in most regions. Thus,
we adopt +$12^{\circ}$ as the LSST’s nominal field-of-view.

Figure 4: Sky positions of interstellar objects when they are detectable for
all sky surveys with limiting magnitudes of 24 and 22. The top panel shows
histograms of the sky positions for the simulated objects that satisfy our
detectability criteria for the LSST the moment that they reach $m\leq 24$. The
nominal field of view of the LSST is over-plotted in red lines. Compared to
Figure 3, the detectable objects are less clustered at the apex due to their
subsequent motion along their trajectories prior to detection. The bottom
panel displays the positions of objects that reach $m\leq 22$ and satisfy our
detectability criteria, when they reach $m\sim 22$. The brightest objects
($m\leq 22$) are preferentially found along/near the ecliptic. The positions
of ‘Oumuamua and Borisov when they were detected are shown, which is distinct
from their positions in Figure 3.

##  4 The Population Detectable with the LSST

In this section, we identify which members of the simulation will be
detectable by the LSST, given its detailed observational capabilities. We
accomplish this using the following procedure. For each simulated object, we
calculate the altitude/azimuth of the object and the Sun at each point along
its 3-year trajectory. We consider an interstellar object detectable by the
LSST if it meets all of the following 3 criteria at any time along its
trajectory through the 5 au sphere:

  1. 1.

Its apparent magnitude is less than 24 and higher than 16. We consider $m=24$
to be the apparent magnitude of the dimmest objects that the LSST is capable
of observing. Because of the LSST’s 30-second exposures, sources brighter than
the 16th magnitude would saturate, meaning they would not be detected by the
LSST.

  2. 2.

Its altitude is higher than $30^{\circ}$ relative to the latitude/longitude of
Rubin Observatory, and its declination is less than $+12^{\circ}$. This
ensures that the object is within the nominal field of view of the LSST, as
the Wide-Fast-Deep (WFD) extends to $+12^{\circ}$ declination for most right
ascensions, according to Figure 1 of Bianco et al. (2021).

  3. 3.

The Sun’s altitude is lower than $-18^{\circ}$ relative to the
latitude/longitude of the LSST. This ensures that the object is observed
before/after astronomical twilight, allowing for detailed observations of
faint point sources.

These criteria ensure that when each object reaches $m\leq 24$, the angular
distance between the object and the Sun is large. We find that out of the
$\sim 10^{6}$ simulated interstellar objects, $6.956\%$ of them reach a
minimum apparent magnitude brighter than 24, and $2.828\%$ of them are
detectable by the LSST. We define a “reachable” object as one that can be
intercepted with a rocket sent from Earth with an impulsive change of
velocity, $\Delta$v$<20$ km/s. We further explore the reachable population in
Section 5.

In Figure 3, we show the position on the sky of the interstellar objects that
are detectable with the LSST when they are entering the Solar System, while in
the top panel of Figure 4, we show the sky positions of the same population of
objects when they reach $m\leq 24$, which corresponds to the LSST’s limiting
magnitude. We include the ecliptic and galactic planes for reference, since
precise observations near the latter are difficult. It is clear that there is
a strong clustering in the sky location of these objects in the vicinity of
the solar apex as they enter the 5 au sphere. At the point within the 5 au
sphere that the objects first become detectable, there is still a clustering
of objects towards the solar apex and anti-apex, but the detection locations
are more isotropic than the locations at which they enter the 5 au sphere.
This is because the object has moved significantly through its trajectory by
the time it is detectable.

In the bottom panel of Figure 4, we show the subset of this population that
reach an apparent magnitude of $m\leq 22$. This subset would be representative
of the detected population in a less sensitive observational search. There is
a strong clustering of detections in the vicinity of the ecliptic as well as
minor groupings of objects around the solar apex and anti-apex, which are
remnants of the incoming distribution. The distribution of locations in
Figures 3 and 4 reveal the likelihood of detecting interstellar objects with
similar trajectories as ‘Oumuamua or Borisov. Because of ‘Oumuamua’s initial
position near the solar apex, one would expect to detect more objects from the
same region in the sky. Borisov’s initial position, far from the solar apex
and anti-apex, was unusual. Because of Borisov’s high incoming velocity, its
age was high. According to Figure 1, older objects have more isotropic
incoming trajectories, explaining Borisov’s behavior. However, its relatively
close approach to Earth and, by extension, the ecliptic, rendered it
detectable.

To explain the clustering of bright objects near the ecliptic, we evaluate the
phase angle, $\theta$, that minimizes the apparent magnitude of an object for
different values of their minimum distance $d_{BO}$. We calculate the apparent
magnitude as a function of phase angle using Equations 8-10, assuming that the
Earth has a circular orbit ($d_{OS}=1$au). In Figure 5 we show the apparent
magnitude as a function of phase angle for three values of $d_{BO}$ which are
representative of the simulated objects that attain $m\sim 22$. In our
simulated population, the vast majority of detected objects do not pass close
to the Earth, and hence $d_{BO}\gg R_{\oplus}$ (Figure 9). A phase angle of
$\theta=0^{\circ}$ minimizes the apparent magnitude for a constant value of
$d_{BO}$. This corresponds to the configuration where the object, Earth, and
Sun are aligned, with the Earth in the center. Therefore, an interstellar
object that achieves a distance from the Earth of $R_{\oplus}\ll d_{BO}<0.5$au
is brightest in the ecliptic plane. The clustering around the ecliptic in the
bottom panel of Figure 4 extends to $\sim\pm 20^{\circ}$. For an object with
$d_{BO}=0.45$au, this corresponds to a phase angle of $\sim 50^{\circ}$,
which, according to Figure 5, produces an apparent magnitude fainter than 22.
Therefore, the clustering of the brightest interstellar objects around the
ecliptic is due the limiting magnitude of 22. This limits the detected objects
to the ones that come within $<0.5$au of the Earth, whose brightness depends
most sensitively on the phase angle.

Figure 5: Apparent magnitude, $m$, as a function of phase angle, $\theta$, for
an ‘Oumuamua like object with $d_{BO}=0.25$au (red), $d_{BO}=0.35$au (green),
and $d_{BO}=0.45$au (blue). The apparent magnitude is minimized for
$\theta\approx 0^{\circ}$, when the object is in the ecliptic plane. The
closer an object comes to the Earth in this range, the more sensitively its
brightness depends on phase angle. Figure 6: Orbital inclination of the
initial population of interstellar objects, and all of the objects that reach
apparent magnitudes of 24 and 21.5. As the magnitude cutoff decreases, the
distribution of inclinations trends towards uniformity. Our distribution for
all objects is similar to that displayed in Figure 2c of Engelhardt et al.
(2017). A typical interstellar object enters the Solar System on a polar
trajectory.

In Figure 6, we show the distribution of orbital inclination for the entire
population of interstellar objects and the objects that reach apparent
magnitudes of 24, the limiting magnitude for the LSST. Since other surveys,
such as Pan-STAARS and the Zwicky Transient Facility (Chambers et al., 2016;
Bellm, 2014), have fainter limiting magnitudes than the LSST, we also show the
distribution for objects that reach apparent magnitudes of 21.5. The initial
and detectable populations are centered at inclinations of $\sim 90^{\circ}$.
This is due to the orientation and velocity of the Solar System with respect
to the galactic mid-plane. There are two compounding properties that produce
inclined orbits for interstellar objects. The first is that the ecliptic plane
is inclined with respect to the galactic mid-plane by $\sim 60^{\circ}$, and
second is that the apex is at $\sim 30^{\circ}$ in declination. Since the
shape of the inclination distributions of the initial conditions and $m\leq
24,21.5$ populations are similar, we conclude that the detectability criteria
is independent of the inclination. Engelhardt et al. (2017) found a $\sin{i}$
inclination distribution that was also maximized at $i=90^{\circ}$. The
differences in the shape of the distributions presented here can be attributed
to the methodology for initializing trajectories.

Figure 7: The seasonal variability with which interstellar objects are
detected by the LSST (representing $m\leq 24$) and Pan-STARRs (representing
$m\leq 19$). We show the fraction of objects that meet observability criteria
and reach $m\leq 24$ and $m\leq 19$ as a function of month detected in blue
and red histograms. The frequency of ISO detections with the LSST is slightly
higher than average close to December/January. Figure 8: Minimum apparent
magnitudes of each interstellar object belonging to the detectable and
reachable ($\Delta v\leq 20$ km/s) populations. The black line indicates the
minimum apparent magnitude that ‘Oumuamua reached. The majority of detectable
ISOs do not get brighter than $m\sim 19$. Figure 9: The minimum distances to
the Earth of each interstellar object in the detectable and reachable
populations. Detectable and reachable objects pass within $\sim 1.5$ and $\sim
1$ au of the Earth, respectively. Figure 10: The wait time, or time between
detection and closest approach to the Earth, of each interstellar object in
the detectable and reachable populations. Most ISOs have wait times of
$\sim$50 days or less, while the upper limit is $\sim 150$ days for both
populations.

In Figure 7, we show the seasonal variability with which the the LSST detects
objects in our simulated population. We also show the seasonal variability for
the subset of these objects that reach $m<19$, corresponding to roughly the
observational capabilities of currently operational all-sky surveys such as
Pan-STARRs. While there appears to be a slight enhancement in the detection
rate in December and January, there is no statistically significant seasonal
variability in detection rate. Our results do not account for the LSST’s
seasonal observational strategy.

In Figure 8, we show histograms of the minimum apparent magnitudes of
simulated interstellar objects. Both the detectable and reachable ($\Delta
v\leq 20$ km/s) populations have distributions that are skewed to fainter
apparent magnitudes. Almost all of the simulated objects ($>99\%$) within both
subgroups do not get brighter than $m\sim 19$. If ‘Oumuamua is a
representative member of its population, then it is not surprising that we
have not detected closer, brighter objects. With the increase in sensitivity
that the LSST represents, we will probe a vastly larger population of
interstellar objects than ever before.

In Figure 9, we show the distributions of the minimum distances to the Earth
for the detectable and reachable interstellar objects. The minimum distance is
calculated by extremizing a cubic spline interpolation on the relative
distance as a function of time. Detectable and reachable objects pass within
$\sim$1.5 and $\sim 1$ au of the Earth, respectively. Although the galactic
population of interstellar objects should produce a $b^{2}$ probability
distribution, the detectable and reachable criteria preferentially selects
objects that come close to the Earth.

In Figure 10, we show the distribution of the wait time for every object, or
time between detection by the LSST and closest approach to the Earth. This
distribution is roughly log normal. The LSST should detect interstellar
objects before they reach their perihelia, with sufficient time to send a
rendezvous mission. It is possible that we will detect ISOs with wait times up
to $\sim$150 days.

##  5 Potential Rendezvous Targets

Figure 11: The $\Delta$v required to guarantee an interception for each
interstellar object belonging to the detectable population. We consider any
object that satisfies $\Delta v<20$ km/s as reachable, which corresponds to
$\sim 20\%$ of detectable objects. Figure 12: Distributions of dynamical
impact parameters for the detectable and reachable populations. Although the
population is drawn from a $b^{2}$ distribution, the LSST will detect objects
out to 3 au. Note that $b$ here represents the orbital impact parameter, and
is not related to a hypothetical impactor mission. Figure 13: Distributions of
perihelia for the detectable and reachable populations. The LSST will detect
objects with perihelia out to 2 au.

In this section we further refine our simulated trajectories to examine the
distribution of objects that are not only detected by the LSST but will also
be feasible targets for an in-situ interception mission. Such a mission could
include imaging, spectroscopy and, if $\Delta$v capabilities are sufficient,
an impactor. An impactor collision of this type would excavate a substantial
amount of sub-surface material from an incoming interstellar object. This
material could be closely examined by a companion flyby spacecraft, with a
suite of instruments similar to those proposed for Bridge (Moore et al.,
2021), allowing for detailed characterization of the composition. Such a
strategy was employed by the Deep Impact mission (A’Hearn et al., 2005).

Trajectories that will be “reachable” from Earth can be categorized by the
magnitude of the impulsive $\Delta$v required to guarantee an interception for
a spacecraft sent from the Earth or the Earth’s Lagrange points L1 or L2. Such
a spacecraft would execute a maneuver that would put it on an intercept
trajectory before its target would reach its closest approach to Earth. While
an ISO intercept mission in this manner would be feasible only for objects
with sufficiently large wait times, Figure 10 demonstrates the existence of a
minimum wait time at $\sim 25$ days for the reachable population. This would
provide a mission with the time needed to design and execute a transfer
maneuver. We note that this is a highly simplified estimation, and that, while
outside of the scope of this paper, it would be worthwhile to investigate the
dynamics of optimal trajectories for such a spacecraft. Seligman & Laughlin
(2018) provided an order-of-magnitude estimate for an attainable $\Delta$v,
using the quoted payload capability of a SpaceX Falcon Heavy to Mars. They
calculated that a maximum $\Delta$v$\sim 15$km/s would be attainable for such
a mission concept to match the impact kinetic energy achieved by the Deep
Impact Tempel I mission.

For our simulated trajectories, we use a fiducial $\Delta v<20$ km/s for
reachable trajectories. This is meant to include the distribution of objects
that may be reachable with additional $\Delta$v from deep space maneuvers. We
estimate the required $\Delta$v for each of our simulated objects by dividing
the minimum distance to the Earth in Figure 9 by the wait time in Figure 10.
The resulting distribution of $\Delta$v is shown in Figure 11. In Figure 8, we
show the minimum apparent magnitude for the subset of objects that reach
$\Delta v<20$ km/s. We find that $\sim 20\%$ of detectable objects are
reachable for such a mission. This implies that roughly 1 in 5 interstellar
objects detected by the LSST will be viable rendezvous targets. This is
significantly better than the estimates in Seligman & Laughlin (2018).

Figure 14: Distributions of eccentricities for the detectable and reachable
populations. The eccentricities of ‘Oumuamua and Borisov are shown in black
and grey dashed lines, respectively. The typical interstellar object has
eccentricity close to unity. Figure 15: The distribution of hyperbolic excess
velocity (incoming velocity) for the detectable and reachable populations. The
incoming velocities of ‘Oumuamua and Borisov are shown in black and grey
dashed lines, respectively.

In Figure 12, we show the distribution of impact parameters for detectable and
reachable objects. The limiting magnitude of the LSST creates a sharp cutoff
at $b\sim 3$ au for objects that will be detectable. If future interstellar
objects have brighter absolute magnitudes than ‘Oumuamua, this cutoff will
increase. These simulations can easily be scaled to probe the impact parameter
cutoff as a function of intrinsic magnitude. In Figure 13, we show the same
distributions, but for perihelia, $r_{min}$, instead of impact parameter. In
terms of perihelia, the LSST cutoff is at $r_{min}\sim 2$au, with this limit
being imposed by the detectability criteria of $m\leq 24$ for an $H=22.4$
object (as per Figure 9, all detectable objects must pass within $\sim 1.5$ au
of the Earth).

In Figures 14 and 15, we show the distributions of eccentricity and hyperbolic
velocity, as in Figure 12. The distribution of eccentricity is roughly log-
normal, while the distribution of hyperbolic velocity has two peaks and
significant substructure. The distributions for objects that are reachable for
an interception mission roughly mirrors those of the detected population. This
implies that the objects that we will be able to probe with an in-situ
measurement roughly represents the detectable sample for these parameters.

H=22.4

Criterion | Percent | Conservative Rate per Year | Optimistic Rate per Year
---|---|---|---
$m\leq 24$ | 6.956% | $\sim 2.3$ | $\sim 4.6$
Detectable by the LSST | 2.828% | $\sim 0.9$ | $\sim 1.9$
Detectable, $\Delta v<30$ km/s | 1.055% | $\sim 0.35$ | $\sim 0.7$
Detectable, $\Delta v<15$ km/s | 0.424% | $\sim 0.1$ | $\sim 0.3$
Detectable, $\Delta v<2$ km/s | 0.002% | $\sim 0.0007$ | $\sim 0.001$

H=23.4

Criterion | Percent | Conservative Rate per Year | Optimistic Rate per Year
---|---|---|---
$m\leq 24$ | 4.241% | $\sim 1.4$ | $\sim 2.8$
Detectable by the LSST | 1.457% | $\sim 0.5$ | $\sim 1$
Detectable, $\Delta v<30$ km/s | 0.793% | $\sim 0.3$ | $\sim 0.5$
Detectable, $\Delta v<15$ km/s | 0.347% | $\sim 0.1$ | $\sim 0.2$
Detectable, $\Delta v<2$ km/s | 0.002% | $\sim 0.0007$ | $\sim 0.001$

H=24.4

Criterion | Percent | Conservative Rate per Year | Optimistic Rate per Year
---|---|---|---
$m\leq 24$ | 1.933% | $\sim 0.6$ | $\sim 1.3$
Detectable by the LSST | 0.378% | $\sim 0.1$ | $\sim 0.25$
Detectable, $\Delta v<30$ km/s | 0.284% | $\sim 0.09$ | $\sim 0.2$
Detectable, $\Delta v<15$ km/s | 0.182% | $\sim 0.06$ | $\sim 0.1$
Detectable, $\Delta v<2$ km/s | 0.002% | $\sim 0.0007$ | $\sim 0.001$

H=25.4

Criterion | Percent | Conservative Rate per Year | Optimistic Rate per Year
---|---|---|---
$m\leq 24$ | 0.739% | $\sim 0.2$ | $\sim 0.5$
Detectable by the LSST | 0.099% | $\sim 0.03$ | $\sim 0.07$
Detectable, $\Delta v<30$ km/s | 0.091% | $\sim 0.03$ | $\sim 0.06$
Detectable, $\Delta v<15$ km/s | 0.073% | $\sim 0.02$ | $\sim 0.05$
Detectable, $\Delta v<2$ km/s | 0.001% | $\sim 0.0003$ | $\sim 0.0007$

Table 2: The annual frequency that interstellar objects are detectable by the
LSST and reachable with various $\Delta$v criteria. Top panel: Statistics for
our sample of interstellar objects, which assumes that each object’s absolute
magnitude is equal to that of ‘Oumuamua. Second panel: Statistics for the same
sample after increasing each object’s absolute magnitude by 1. Third panel:
Statistics for the original sample after increasing each object’s absolute
magnitude by 2. Bottom panel: Statistics for the original sample after
increasing each object’s absolute magnitude by 3. Our percentages are
significantly different even if we assume that each object’s absolute
magnitude is greater than that of ‘Oumuamua by 1. The frequency of detectable
objects is therefore sensitive to their intrinsic absolute magnitude.

In Table 2, we show the fractional percentages of objects that will be
detected by the LSST and attainable targets for a range of limiting $\Delta$v.
To convert these fractions to detection rates, we consider a typical ISO
travelling within the 5 au sphere centered at the Sun. We assume a spatial
number density of $n_{o}\sim 1\times 10^{-1}\,$au-3, lower than that
calculated by Do et al. (2018) given the non-detection of an additional
‘Oumuamua like object since 2017 (Levine et al., 2021). With this assumed
spatial density, there are $\sim 50$ interstellar objects within the 5au
sphere at any given time.

The crossing time of the 5au sphere varies based on the trajectory and
incoming velocity vector of a given interstellar object. We use incoming
velocities as our values for mean crossing velocities, implicitly averaging
over the Solar acceleration. Although the distribution of hyperbolic excess
velocity does not decrease monotonically, its frequency is highest at $\sim
20$ km/s and begins to decrease at $\sim 40$ km/s, as demonstrated in Figure
15. Therefore, we use mean crossing velocities of 20 and 40 km/s to represent
our conservative and optimistic predictions for the detection frequency.

The typical distance traveled can be approximated as the average distance,
$||d||$, of a chord between two points on a sphere, as described by Solomon
(1978). For a 2-dimensional circle, this may be calculated using the following
equation,

$\begin{split}||d||=\frac{1}{(2\pi)^{2}}\int_{0}^{2\pi}\int_{0}^{2\pi}\,\bigg{(}2-2\cos{\big{(}\theta_{1}-\theta_{2}\big{)}}\bigg{)}^{1/2}r\,d\theta_{1}d\theta_{2}\,,\end{split}$
(11)

which yields $||d||=4/\pi r$. For a sphere,

$\begin{split}||d||=\frac{1}{(4\pi)^{2}}\int_{0}^{2\pi}\int_{0}^{\pi}\int_{0}^{2\pi}\int_{0}^{\pi}\,\bigg{(}2-2\cos{\phi_{1}}\cos{\phi_{2}}\\\
-2\cos{\theta_{1}}\cos{\theta_{2}}\sin{\phi_{1}}\sin{\phi_{2}}-2\sin{\theta_{1}}\sin{\theta_{2}}\sin{\phi_{1}}\sin{\phi_{2}}\bigg{)}^{1/2}\\\
r\,\sin{\phi_{1}}\sin{\phi_{2}}d\phi_{1}d\theta_{1}d\phi_{2}d\theta_{2}\,,\end{split}$
(12)

which gives $||d||=4/3r$. We verified these values numerically by calculating
the average of sets of points drawn uniformly from the boundary of the unit
circle and sphere.

It is noteworthy that the average chord length depends on the method by which
each chord is constructed (see Table 5 in Solomon (1978)). For the case
considered here, where two random points on the surface of the sphere are
connected, $||d||=4/3r$, which is also 4 times the volume divided by the
surface area. If the chord is generated from a normal to the plane of a random
great circle, through a random point in the circle of intersection,
$||d||=4/3r$ still. However, if the chords are formed by choosing a single
random point on the sphere, and then a random direction vector, $||d||=r$. If
the chords are generated by selecting two random points within the sphere,
then $||d||=12/7r$. For our problem, the method of selecting two random points
on the surface is most appropriate. This is a simplified version of the set of
problems considered by Alagar (1976); Dunbar (1997); Burgstaller &
Pillichshammer (2009), who calculated similar quantities for higher
dimensional cases, such as for a hypersphere $X\subseteq\mathbb{R}^{s}$.
Interestingly, this geometric probability problem regarding the distribution
of chords within a sphere was developed as an application to cellular biology.
Specifically, this formalism was required to quantify the extent with which
members of pairs of chromosomes were positioned randomly with respect to other
members during mitosis (Barton et al., 1963; David & Fix, 1964).

These values correspond to crossing times of $\sim 0.8$ and $\sim 1.6$ years.
To convert the fractions to detection rates, we use the following equation,

$\begin{split}{\rm rate}=33.11\,{\rm yr}^{-1}\bigg{(}\frac{n_{o}}{1\times
10^{-1}\,\rm{au}^{-3}}\bigg{)}\\\
\bigg{(}\frac{V}{4/3\pi(5\rm{au})^{3}}\bigg{)}\bigg{(}\frac{v}{20\,{\rm
km/s}}\bigg{)}\bigg{(}\frac{(4/3)\,5{\rm au}}{||d||}\bigg{)}\,,\end{split}$
(13)

where $V$ is the simulated volume and v is the assumed typical velocity.
Therefore, we multiply the percentages in Table 2 by $\sim 66$ and $\sim 33$
objects to represent the yearly steady state population of the 5au sphere.

We estimate that the LSST should detect between $0.9-1.9$ interstellar objects
like ‘Oumuamua ($H=22.4$) every year. Therefore, by the end of its 10 year
observational campaign, the census of detected interstellar objects should be
roughly $9-19$ in magnitude, which is a significant increase from previous
estimates. Moreover, with a dedicated mission such as Bridge, which should be
able to generate a $\Delta$v$\sim 15$km/s, there should be of order $\sim 1-3$
objects detected in the lifetime of the LSST that would be attainable targets.
If the Comet Interceptor is capable of generating a $\Delta$v$\sim 2$km/s,
then it appears that there is a $\sim 0.07\%$ chance that a target is
detected.

It is important to note that these estimates do not take into account the
size-frequency distribution of interstellar objects and the possibility of
cometary activity, as discussed previously in this paper. This could provide a
significant increase to the number of detected interstellar objects and
feasible targets for rendezvous missions.

However, we consider the H-frequency distribution of interstellar objects by
running our simulations with the assumption that the absolute magnitudes of
each object are higher (fainter) than that of ‘Oumuamua by 1, 2, and 3. Table
2 also displays the fractional percentages of objects that meet various
criteria under these assumptions. We find that the percentages are sensitive
to our initial choice of absolute magnitude. For example, while the percentage
of objects that reach $m\leq 24$ is 6.956% if we assume an absolute magnitude
of 22.4, the percentage of objects that meet the same criterion decreases by a
factor of $\sim 1.6$ if we assume an absolute magnitude of 23.4. Until an
intrinsic absolute magnitude distribution can be determined for ISOs,
population estimates will vary significantly depending on the assumptions
made.

These detection rates are significantly higher than the ones calculated in
Trilling et al. (2017) and Cook et al. (2016). The rates are still very
uncertain, since they are based on the detection of only two objects. Future
detections of interstellar objects will refine our estimates of spatial number
density in the inner Solar System. These simulations and the resulting
percentages in Table 2 can be readily updated when the intrinsic spatial
number density is better constrained.

Because the LSST is capable of surveying the entire sky in a few days, we
assumed an ISO detection efficiency of $100\%$. However, it is likely that the
detection efficiency is lower, which would decrease all of our estimates
presented in Table 2. However, the ability of the LSST to detect transient
objects has been demonstrated for the Near Earth Objects, JFCs and LPCs
(Solontoi et al., 2011; Vereš & Chesley, 2017a, b; Jones et al., 2018).
Engelhardt et al. (2017) found that there is a $\sim 35\%$ chance that a
detectable ISO will not be recognized as such due to lack of follow-up
observations. However, their investigation involved surveys that have lower
ISO detection efficiencies than those of the LSST. If we adopt $65\%$ as a
minimum value for the LSST’s ISO detection efficiency, then we still expect
that the LSST will detect at least $\sim 6$ ISOs. The detection of these
objects will drastically expand our knowledge of this novel field of
astrophysics and planetary science.

##  6 Conclusions

In this paper, we simulated a population of interstellar objects drawn from
their galactic kinematic distribution. Using this synthetic population, we
evaluated the distribution of interstellar objects that should be detectable
with the forthcoming LSST. We showed that the incoming interstellar objects
are strongly clustered in the direction of the solar apex. The distribution of
detectable locations on the sky is more isotropic than the distribution of
initial locations with slight enhancements in the vicinity of the apex and
anti-apex. This is due to the fact that by the time that the objects are
detectable by the LSST, they have moved significantly through their
trajectories. We showed that less sensitive all sky surveys are more biased to
detecting objects close to the ecliptic.

We calculated the distribution of orbital trajectories for the objects that
are detected by the LSST. Moreover, we evaluated a subset of these objects
that will be reachable for in-situ interception missions. We estimated that
the LSST should detect of order $\sim 9-19$ interstellar objects every year,
and that $\sim 3-7$ reachable targets with a dedicated mission, such as Bridge
will be detected over the survey’s lifetime. There is a $\sim 0.07\%$ chance
that a target for the Comet Interceptor will be detected.

For a larger portion of interstellar objects to be reachable, more efficient
technologies, such as solar sails, are required. Linares et al. (2020)
proposed a “statite” (static-satellite) concept in which solar sails are
utilized to maintain a probe’s orbit and to power an intercept trajectory.
Such a mission could hypothetically achieve $\Delta v\approx 100$ km/s, far
higher than the $\Delta$v requirements of Bridge or Comet Interceptor.
However, these missions require that the orbits of the potential target be
known 4-16 months in advance. Figure 10 establishes a tentative upper limit on
the wait times of reachable objects of $\sim 5$ months. Therefore, a statite
must be in the optimal position prior to its intercept maneuver. This type of
mission would most likely require a pre-existing cluster of statites in Earth-
like heliocentric orbits to be feasible.

Since a dedicated ISO survey would have virtually the same methodology as a
dedicated NEO survey, existing NEO surveys may increase the number of detected
ISOs beyond the number of ISOs that we expect the LSST to detect. Since
‘Oumuamua and Borisov were detected using Pan-STARRS and the Crimean
Astrophysical Observatory (Meech et al., 2017; de León et al., 2019), an
international observational campaign may be a worthwhile use of resources.
Dedicated surveys with limiting magnitudes around or fainter than 22, such as
Pan-STAARS or ATLAS (Chambers et al., 2016; Tonry et al., 2018), may observe
an ISO while searching in the vicinity of the ecliptic. Targeted searches that
are equally or more sensitive than the LSST could focus their searches on the
vicinity of the solar apex.

##  7 Acknowledgements

This work utilized the computational resources of the NIH HPC Biowulf cluster
(http://hpc.nih.gov). We thank David Hoover for access to the Biowulf
computing cluster. We thank Adina Feinstein, Amir Siraj, Jacob Bean, Darin
Ragozzine, and Ivan Levcovitz for useful conversations and suggestions. We
acknowledge support from the University of Chicago Quad Faculty Research Grant
program and Dean’s Fund for Undergraduate Research. We also thank the two
anonymous reviewers, whose reports led to modifications that improved the
quality of this paper.

## References

  * A’Hearn et al. (2005) A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2005, Science, 310, 258, doi: 10.1126/science.1118923
  * Alagar (1976) Alagar, V. S. 1976, Journal of Applied Probability, 13, 558–566, doi: 10.2307/3212475
  * Almeida-Fernandes & Rocha-Pinto (2018) Almeida-Fernandes, F., & Rocha-Pinto, H. J. 2018, MNRAS, 480, 4903, doi: 10.1093/mnras/sty2202
  * Aravind et al. (2021) Aravind, K., Ganesh, S., Venkataramani, K., et al. 2021, MNRAS, 502, 3491, doi: 10.1093/mnras/stab084
  * Bagnulo et al. (2021) Bagnulo, S., Cellino, A., Kolokolova, L., et al. 2021, Nature Communications, 12, 1797, doi: 10.1038/s41467-021-22000-x
  * Bannister et al. (2017) Bannister, M. T., Schwamb, M. E., Fraser, W. C., et al. 2017, ApJL, 851, L38, doi: 10.3847/2041-8213/aaa07c
  * Bannister et al. (2020) Bannister, M. T., Opitom, C., Fitzsimmons, A., et al. 2020, arXiv e-prints, arXiv:2001.11605.  https://arxiv.org/abs/2001.11605
  * Barton et al. (1963) Barton, D. E., David, F. N., & Fix, E. 1963, Biometrika, 50, 23.  http://www.jstor.org/stable/2333743
  * Bellm (2014) Bellm, E. 2014, in The Third Hot-wiring the Transient Universe Workshop, ed. P. R. Wozniak, M. J. Graham, A. A. Mahabal, & R. Seaman, 27–33.  https://arxiv.org/abs/1410.8185
  * Belton et al. (2018) Belton, M. J. S., Hainaut, O. R., Meech, K. J., et al. 2018, ApJL, 856, L21, doi: 10.3847/2041-8213/aab370
  * Bialy & Loeb (2018) Bialy, S., & Loeb, A. 2018, ApJ, 868, L1
  * Bianco et al. (2021) Bianco, F. B., Ivezić, Ž., Jones, R. L., et al. 2021, arXiv e-prints, arXiv:2108.01683.  https://arxiv.org/abs/2108.01683
  * Binney & Merrifield (1998) Binney, J., & Merrifield, M. 1998, Galactic Astronomy
  * Bodewits et al. (2020) Bodewits, D., Noonan, J. W., Feldman, P. D., et al. 2020, Nature Astronomy, 4, 867, doi: 10.1038/s41550-020-1095-2
  * Bolin et al. (2020a) Bolin, B. T., Bodewits, D., Lisse, C. M., et al. 2020a, The Astronomer’s Telegram, 13613, 1
  * Bolin et al. (2018) Bolin, B. T., Weaver, H. A., Fernandez, Y. R., et al. 2018, ApJL, 852, L2, doi: 10.3847/2041-8213/aaa0c9
  * Bolin et al. (2020b) Bolin, B. T., Lisse, C. M., Kasliwal, M. M., et al. 2020b, AJ, 160, 26, doi: 10.3847/1538-3881/ab9305
  * Burgstaller & Pillichshammer (2009) Burgstaller, B., & Pillichshammer, F. 2009, Bulletin of the Australian Mathematical Society, 80, 353–359, doi: 10.1017/S0004972709000707
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints.  https://arxiv.org/abs/1612.05560
  * Cook et al. (2016) Cook, N. V., Ragozzine, D., Granvik, M., & Stephens, D. C. 2016, ApJ, 825, 51, doi: 10.3847/0004-637X/825/1/51
  * Cordiner et al. (2020) Cordiner, M. A., Milam, S. N., Biver, N., et al. 2020, Nature Astronomy, 4, 861, doi: 10.1038/s41550-020-1087-2
  * Cremonese et al. (2020) Cremonese, G., Fulle, M., Cambianica, P., et al. 2020, ApJ, 893, L12, doi: 10.3847/2041-8213/ab8455
  * David & Fix (1964) David, F. N., & Fix, E. 1964, Biometrika, 51, 373.  http://www.jstor.org/stable/2334143
  * de la Fuente Marcos & de la Fuente Marcos (2020) de la Fuente Marcos, C., & de la Fuente Marcos, R. 2020, A&A, 643, A18, doi: 10.1051/0004-6361/202037447
  * de León et al. (2019) de León, J., Licandro, J., Serra-Ricart, M., et al. 2019, Research Notes of the American Astronomical Society, 3, 131, doi: 10.3847/2515-5172/ab449c
  * Desch & Jackson (2021) Desch, S. J., & Jackson, A. P. 2021, Journal of Geophysical Research: Planets, e2020JE006807
  * Do et al. (2018) Do, A., Tucker, M. A., & Tonry, J. 2018, ApJ, 855, L10, doi: 10.3847/2041-8213/aaae67
  * Drahus et al. (2018) Drahus, M., Guzik, P., Waniak, W., et al. 2018, Nature Astronomy, 2, 407, doi: 10.1038/s41550-018-0440-1
  * Drahus et al. (2020) Drahus, M., Guzik, P., Udalski, A., et al. 2020, The Astronomer’s Telegram, 13549, 1
  * Dunbar (1997) Dunbar, S. R. 1997, The College Mathematics Journal, 28, 187, doi: 10.1080/07468342.1997.11973860
  * Engelhardt et al. (2017) Engelhardt, T., Jedicke, R., Vereš, P., et al. 2017, AJ, 153, 133, doi: 10.3847/1538-3881/aa5c8a
  * Fitzsimmons et al. (2018) Fitzsimmons, A., Snodgrass, C., Rozitis, B., et al. 2018, Nature Astronomy, 2, 133, doi: 10.1038/s41550-017-0361-4
  * Fitzsimmons et al. (2019) Fitzsimmons, A., Hainaut, O., Meech, K. J., et al. 2019, ApJ, 885, L9, doi: 10.3847/2041-8213/ab49fc
  * Füglistaler & Pfenniger (2015) Füglistaler, A., & Pfenniger, D. 2015, A&A, 578, A18
  * Grude Flekkøy & Brodin (2022) Grude Flekkøy, E., & Brodin, J. 2022, arXiv e-prints, arXiv:2201.10327.  https://arxiv.org/abs/2201.10327
  * Guzik & Drahus (2021) Guzik, P., & Drahus, M. 2021, Nature, 593, 375, doi: 10.1038/s41586-021-03485-4
  * Guzik et al. (2020) Guzik, P., Drahus, M., Rusek, K., et al. 2020, Nature Astronomy, 4, 53, doi: 10.1038/s41550-019-0931-8
  * Hallatt & Wiegert (2020) Hallatt, T., & Wiegert, P. 2020, AJ, 159, 147
  * Hands & Dehnen (2020) Hands, T. O., & Dehnen, W. 2020, MNRAS, 493, L59, doi: 10.1093/mnrasl/slz186
  * Hsieh et al. (2021) Hsieh, C.-H., Laughlin, G., & Arce, H. G. 2021, ApJ, 917, 20, doi: 10.3847/1538-4357/ac0729
  * Hui et al. (2020) Hui, M.-T., Ye, Q.-Z., Föhring, D., Hung, D., & Tholen, D. J. 2020, AJ, 160, 92, doi: 10.3847/1538-3881/ab9df8
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jackson & Desch (2021) Jackson, A. P., & Desch, S. J. 2021, Journal of Geophysical Research: Planets, e2020JE006706
  * Jewitt et al. (2020a) Jewitt, D., Kim, Y., Mutchler, M., et al. 2020a, ApJ, 896, L39, doi: 10.3847/2041-8213/ab99cb
  * Jewitt & Luu (2019) Jewitt, D., & Luu, J. 2019, ApJ, 886, L29, doi: 10.3847/2041-8213/ab530b
  * Jewitt et al. (2017) Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJ, 850, L36
  * Jewitt et al. (2020b) Jewitt, D., Mutchler, M., Kim, Y., Weaver, H., & Hui, M.-T. 2020b, The Astronomer’s Telegram, 13611, 1
  * Jones & ESA Comet Interceptor Team (2019) Jones, G., & ESA Comet Interceptor Team. 2019, Comet Interceptor A Mission to a Dynamically New Solar System Object.  http://www.cometinterceptor.space/uploads/1/2/3/7/123778284/comet_interceptor_executive_summary.pdf
  * Jones et al. (2009) Jones, R. L., Chesley, S. R., Connolly, A. J., et al. 2009, Earth Moon and Planets, 105, 101, doi: 10.1007/s11038-009-9305-z
  * Jones et al. (2018) Jones, R. L., Slater, C. T., Moeyens, J., et al. 2018, Icarus, 303, 181, doi: 10.1016/j.icarus.2017.11.033
  * Kareta et al. (2020) Kareta, T., Andrews, J., Noonan, J. W., et al. 2020, ApJ, 889, L38, doi: 10.3847/2041-8213/ab6a08
  * Kim et al. (2020) Kim, Y., Jewitt, D., Mutchler, M., et al. 2020, ApJ, 895, L34, doi: 10.3847/2041-8213/ab9228
  * Knight et al. (2017) Knight, M. M., Protopapa, S., Kelley, M. S. P., et al. 2017, ApJL, 851, L31, doi: 10.3847/2041-8213/aa9d81
  * Laughlin & Batygin (2017) Laughlin, G., & Batygin, K. 2017, Research Notes of the American Astronomical Society, 1, 43, doi: 10.3847/2515-5172/aaa02b
  * Levine et al. (2021) Levine, W. G., Cabot, S. H. C., Seligman, D., & Laughlin, G. 2021, arXiv e-prints, arXiv:2108.11194.  https://arxiv.org/abs/2108.11194
  * Levine & Laughlin (2021) Levine, W. G., & Laughlin, G. 2021, ApJ, 912, 3, doi: 10.3847/1538-4357/abec85
  * Lin et al. (2020) Lin, H. W., Lee, C.-H., Gerdes, D. W., et al. 2020, ApJ, 889, L30, doi: 10.3847/2041-8213/ab6bd9
  * Linares et al. (2020) Linares, R., Landau, D., Miller, D., Weiss, B., & Lozano, P. 2020, arXiv e-prints, arXiv:2012.12935.  https://arxiv.org/abs/2012.12935
  * Luu et al. (2020) Luu, J. X., Flekkøy, E. G., & Toussaint, R. 2020, ApJ, 900, L22
  * Mamajek (2017) Mamajek, E. 2017, Research Notes of the American Astronomical Society, 1, 21, doi: 10.3847/2515-5172/aa9bdc
  * Manzini et al. (2020) Manzini, F., Oldani, V., Ochner, P., & Bedin, L. R. 2020, MNRAS, 495, L92, doi: 10.1093/mnrasl/slaa061
  * Mashchenko (2019) Mashchenko, S. 2019, MNRAS, 489, 3003, doi: 10.1093/mnras/stz2380
  * Masiero (2017) Masiero, J. 2017, arXiv e-prints, arXiv:1710.09977.  https://arxiv.org/abs/1710.09977
  * McKay et al. (2020) McKay, A. J., Cochran, A. L., Dello Russo, N., & DiSanti, M. A. 2020, ApJ, 889, L10, doi: 10.3847/2041-8213/ab64ed
  * McNeill et al. (2018) McNeill, A., Trilling, D. E., & Mommert, M. 2018, ApJL, 857, L1, doi: 10.3847/2041-8213/aab9ab
  * Meech et al. (2017) Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378, doi: 10.1038/nature25020
  * Micheli et al. (2018) Micheli, M., Farnocchia, D., Meech, K. J., et al. 2018, Nature, 559, 223, doi: 10.1038/s41586-018-0254-4
  * Moore et al. (2021) Moore, K., Courville, S., Ferguson, S., et al. 2021, Planet. Space Sci., 197, 105137, doi: 10.1016/j.pss.2020.105137
  * 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 (2018) Moro-Martín, A. 2018, ApJ, 866, 131, doi: 10.3847/1538-4357/aadf34
  * Moro-Martín (2019) —. 2019, AJ, 157, 86, doi: 10.3847/1538-3881/aafda6
  * Moro-Martín (2019) Moro-Martín, A. 2019, ApJ, 872, L32
  * Napier et al. (2021) Napier, K. J., Adams, F. C., & Batygin, K. 2021, , 2, 53, doi: 10.3847/PSJ/abe76e/53
  * Opitom et al. (2019) Opitom, C., Fitzsimmons, A., Jehin, E., et al. 2019, A&A, 631, L8, doi: 10.1051/0004-6361/201936959
  * Pau Sánchez et al. (2021) Pau Sánchez, J., Morante, D., Hermosin, P., et al. 2021, arXiv e-prints, arXiv:2107.12999.  https://arxiv.org/abs/2107.12999
  * Price et al. (2021) Price, E. M., Cleeves, L. I., Bodewits, D., & Öberg, K. I. 2021, ApJ, 913, 9, doi: 10.3847/1538-4357/abf041
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Sekanina (2019) Sekanina, Z. 2019, arXiv e-prints, arXiv:1901.08704.  https://arxiv.org/abs/1901.08704
  * Seligman & Laughlin (2018) Seligman, D., & Laughlin, G. 2018, AJ, 155, 217, doi: 10.3847/1538-3881/aabd37
  * Seligman & Laughlin (2020) —. 2020, ApJ, 896, L8
  * Seligman et al. (2019) Seligman, D., Laughlin, G., & Batygin, K. 2019, ApJ, 876, L26
  * Seligman et al. (2021) Seligman, D., Levine, W. G., Cabot, S. H. C., Laughlin, G., & Meech, K. 2021, arXiv e-prints, arXiv:2107.06834.  https://arxiv.org/abs/2107.06834
  * Siraj & Loeb (2019) Siraj, A., & Loeb, A. 2019, ApJ, 872, L10, doi: 10.3847/2041-8213/ab042a
  * Solomon (1978) Solomon, H. 1978, 6. Random Chords in the Circle and the Sphere, 127–172, doi: 10.1137/1.9781611970418.ch6
  * Solontoi et al. (2011) Solontoi, M., Ivezić, Ž., & Jones, L. 2011, in American Astronomical Society Meeting Abstracts# 217, Vol. 217, 252–11
  * Tonry et al. (2018) Tonry, J. L., Denneau, L., Heinze, A. N., et al. 2018, PASP, 130, 064505, doi: 10.1088/1538-3873/aabadf
  * Trilling et al. (2017) Trilling, D. E., Robinson, T., Roegge, A., et al. 2017, ApJ, 850, L38, doi: 10.3847/2041-8213/aa9989
  * Trilling et al. (2018) Trilling, D. E., Mommert, M., Hora, J. L., et al. 2018, AJ, 156, 261
  * Vereš & Chesley (2017a) Vereš, P., & Chesley, S. R. 2017a, AJ, 154, 13, doi: 10.3847/1538-3881/aa73d0
  * Vereš & Chesley (2017b) —. 2017b, AJ, 154, 12, doi: 10.3847/1538-3881/aa73d1
  * Xing et al. (2020) Xing, Z., Bodewits, D., Noonan, J., & Bannister, M. T. 2020, ApJ, 893, L48, doi: 10.3847/2041-8213/ab86be
  * Yang et al. (2021) Yang, B., Li, A., Cordiner, M. A., et al. 2021, Nature Astronomy, doi: 10.1038/s41550-021-01336-w
  * Ye et al. (2020) Ye, Q., Kelley, M. S. P., Bolin, B. T., et al. 2020, AJ, 159, 77, doi: 10.3847/1538-3881/ab659b
  * Ye et al. (2017) Ye, Q.-Z., Zhang, Q., Kelley, M. S. P., & Brown, P. G. 2017, ApJ, 851, L5, doi: 10.3847/2041-8213/aa9a34
  * Zhang et al. (2020) Zhang, Q., Ye, Q., & Kolokolova, L. 2020, The Astronomer’s Telegram, 13618, 1