# The Galactic Interstellar Object Population: A Framework for Prediction and
Inference

Matthew J. Hopkins Department of Physics, University of Oxford, Denys
Wilkinson Building, Keble Road, Oxford, OX1 3RH, UK  Chris Lintott Department
of Physics, University of Oxford, Denys Wilkinson Building, Keble Road,
Oxford, OX1 3RH, UK  Michele T. Bannister School of Physical and Chemical
Sciences|Te Kura Matū, University of Canterbury, Private Bag 4800,
Christchurch 8140, New Zealand  J. Ted Mackereth Banting Fellow  Just Group
plc, Enterprise House, Bancroft road, Reigate, Surrey RH2 7RP, UK  Canadian
Institute for Theoretical Astrophysics, University of Toronto, 60 St. George
Street, Toronto, ON, M5S 3H8, Canada  Dunlap Institute for Astronomy and
Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S
3H4, Canada  John C. Forbes School of Physical and Chemical Sciences|Te Kura
Matū, University of Canterbury, Private Bag 4800, Christchurch 8140, New
Zealand

###### Abstract

The Milky Way is thought to host a huge population of interstellar objects
(ISOs), numbering approximately $\qty{1e15}{\per\cubed}$ around the Sun, which
are formed and shaped by a diverse set of processes ranging from planet
formation to galactic dynamics. We define a novel framework: firstly to
predict the properties of this Galactic ISO population by combining models of
processes across planetary and galactic scales, and secondly to make
inferences about the processes modelled, by comparing the predicted population
to what is observed. We predict the spatial and compositional distribution of
the Galaxy's population of ISOs by modelling the Galactic stellar population
with data from the APOGEE survey and combining this with a protoplanetary disk
chemistry model. Selecting ISO water mass fraction as an example observable
quantity, we evaluate its distribution both at the position of the Sun and
averaged over the Galactic disk; our prediction for the Solar neighbourhood is
compatible with the inferred water mass fraction of 2I/Borisov. We show that
the well-studied Galactic stellar metallicity gradient has a corresponding ISO
compositional gradient. We also demonstrate the inference part of the
framework by using the current observed ISO composition distribution to
constrain the parent star metallicity dependence of the ISO production rate.
This constraint, and other inferences made with this framework, will improve
dramatically as the Vera C. Rubin Observatory Legacy Survey of Space and Time
(LSST) progresses and more ISOs are observed. Finally, we explore
generalizations of this framework to other Galactic populations, such as that
of exoplanets.

Interstellar objects (52), Small Solar System bodies(1469), Galaxy Evolution
(594)

††software: NumPy (Harris et al., 2020); SciPy (Virtanen et al., 2020);
Astropy (Astropy Collaboration et al., 2013, 2018, 2022); Matplotlib (Hunter,
2007)

##  1 Introduction

1I/‘Oumuamua (Meech et al., 2017) and
2I/Borisov111https://minorplanetcenter.net/mpec/K19/K19RA6.html and
https://minorplanetcenter.net/mpec/K19/K19S72.html are the first two observed
samples from a highly numerous population: interstellar objects (ISOs).
Estimated to number $\sim\qty{1e15}{\per\cubed}$ around the Sun (Engelhardt et
al., 2017; Do et al., 2018), they are implied to have a spatial distribution
spanning the entire Galaxy. This population has been predicted to exist for
decades (McGlynn & Chapman, 1989), based on models of the accretion and
migration of the giant planets, which predict that 75–85% of cometary bodies
initially in the Solar System must have been scattered into interstellar space
(Fernandez & Ip, 1984; Brasser et al., 2006). Modern exoplanet surveys
consistently find that giant planets are common across the Galaxy around stars
with a range of spectral types (Fulton et al., 2021; Sabotta et al., 2021).
This makes planetesimal scattering common across the Galaxy. A significant
number of planetesimals can also be ejected by close stellar flybys early in a
planetary system's life (e.g. Pfalzner et al., 2021). The protoplanetary disks
of other stars are therefore expected to be a source of ISOs (Stern, 1990;
Moro-Martín, 2022).

Initially it was expected that interstellar objects would display cometary
characteristics (e.g. Jewitt, 2003). The population's dominant dynamical
formation mechanisms would preferentially harvest more distant, ice-rich
planetesimals from the disks of the source systems. More of the cometary ISO
population passing through the Solar System could be detected than rocky ISOs,
as comae make cometary ISOs brighter down to smaller diameters (Engelhardt et
al., 2017). 2I/Borisov appeared relatively similar in size and composition to
Solar System comets (Jewitt & Seligman, 2022). Its distinctive features were
an exceptionally high CO and NH2 content, implying it formed on the edge of
its home system's protoplanetary disk, beyond the hypervolatile CO ice line
(Bodewits et al., 2020; Cordiner et al., 2020).

However, 1I/‘Oumuamua had a mix of observed characteristics. A 160-m scale
object, it lacked a coma in deep imaging (Jewitt et al., 2017) or detectable
outgassing in CO, C$\mathrm{O}_{2}$ (Trilling et al., 2018) or CN (Ye et al.,
2017). Despite this, it still underwent non-gravitational acceleration similar
to that experienced by Solar System comets (Micheli et al., 2018). The large
amplitude of its light curve implied a high-aspect-ratio shape (Mashchenko,
2019). At the time, this combination of factors seemed unusual, although the
surface reflectance properties were consistent with outer Solar System bodies
(e.g. Bannister et al., 2017). Hypotheses for the composition and formation of
1I/‘Oumuamua that match the limited data remain varied. It may be a
planetesimal ('Oumuamua ISSI Team et al., 2019); a fragment of a comet
devolatized by passages close to its parent star before ejection (Raymond et
al., 2018a); an icy fractal aggregate (Moro-Martín, 2019); a hydrogen iceberg
formed in a molecular cloud (Seligman & Laughlin, 2020); or a nitrogen ice
fragment from the surface of a Pluto-like dwarf planet (Jackson & Desch,
2021). Recently, observation of similarly-small near-Earth asteroids have
identified objects with the lightcurve amplitudes seen in 1I. Additionally,
Farnocchia et al. (2023) and Seligman et al. (2023) report six asteroids with
significant non-gravitational acceleration and no coma. While from the two
ISOs found so far, the compositions of interstellar objects are clearly
varied, 1I may be less extreme than it was first considered.

ISOs formed in a protoplanetary disk carry information about their home
systems in their composition: we have samples of other planetary systems
coming to us for study. The composition of a protoplanetary disk correlates
with the elemental abundances of its central star, due to their formation from
the same gas and dust in a molecular cloud core (Öberg & Bergin, 2021). We can
thus expect stars of different metallicities to produce ISOs of different
compositions. The composition of this gas and dust varies in both space and
time as the Galaxy chemically evolves due to stellar nucleosynthesis, making
the current Galactic ISO population sensitive to the entire history and
evolution of the Milky Way over cosmic time (Tinsley & Cameron, 1974; Lintott
et al., 2022). Since the occurrence of planetesimal-scattering giant planets
also has a metallicity dependence (Fischer & Valenti, 2005), the relative
occurrence of ISOs of different compositions will therefore carry information
about the Galaxy's distribution of planetary architectures. Finally, we expect
many ISOs to outlive their parent stars, as they occupy a similar environment
to Gyr-old Oort cloud comets. At minor-planet sizes, ISOs are not subject to
any known destructive forces in the ISM (Guilbert-Lepoutre et al., 2015),
other than minor erosion by dust, in their frequent passages through molecular
clouds (Pfalzner et al., 2020). They could also be disrupted or devolatilised
in rare close encounters with stars (Raymond et al., 2018a; Forbes & Loeb,
2019). ISOs thus present the possibility of studying long-lost planetary
systems.

The Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST) (Ivezić
et al., 2019) is predicted to provide a sample of tens of 1I/‘Oumuamua-like
ISOs (e.g. Levine et al., 2021), as well as any more 2I/Borisov-like cometary
ISOs that enter our Solar System. This is in addition to the continuing
contributions of the NEO surveys and other observatories that found 1I and 2I
in the first place (Meech et al., 2017). These will provide a large and varied
sample of interstellar objects, for study and comparison to predictions.

The many dependencies over planetary and Galactic scales make the observed ISO
population a fascinating tool: it has the potential to test models in both
Galaxy and planetary physics, with an entirely different set of biases to
traditional methods. Lintott et al. (2022) introduced the concept of
predicting the composition of a Galaxy's population of ISOs from the Galaxy's
stellar distribution. Using simulated Galaxies from EAGLE (Schaye et al.,
2015), they showed that the water content of ISOs was sensitive to Galactic
star formation history. In this work, we develop this method and apply it to
the stellar population of the Milky Way, estimated with data from the APOGEE
survey, to predict a broader set of properties of our own Galaxy's population
of interstellar objects. We predict the distribution of ISOs in both their
spatial position in the Galaxy and their water mass fraction. By evaluating
this distribution at the current position of the Solar System, we predict the
properties of the population of ISOs from which the observed sample of the
2020s will be drawn, and we compare this to the whole-Galaxy distribution. We
then detail a Bayesian method of comparing the predicted and observed
distributions to make inferences about the planetary and Galactic processes
modelled, and demonstrate this method by constraining the metallicity
dependence of the ISO production rate.

##  2 APOGEE and Stellar Density Modelling

To predict the distribution of ISOs in the Milky Way, we first obtain the
distribution of all stars throughout Galactic history over a large swath of
the Galactic disk, which we model by fitting simple density profiles to
debiased data from the APOGEE survey. While APOGEE's main sample is not
representative of all extant stars, we can extrapolate from it to recover the
total stellar population of the Milky Way. By design, the APOGEE main sample
mainly contains red giants: stars in a relatively short-lived phase towards
the end of their lives. Since each stellar generation forms stars with a range
of masses (Chabrier, 2003) and a corresponding range of lifespans, all but the
newest stellar populations will have some stars currently in the red giant
stage. This means multiple generations are represented in the APOGEE sample.
We then extrapolate to the entire stellar population. This reconstruction is
detailed in 2.3.

###  2.1 Observational Data: APOGEE

We use data from the APOGEE SDSS-IV Data Release 16 (Jönsson et al., 2020).
APOGEE is a near-infrared, high-resolution ($R\sim$22\,500$$) spectroscopic
stellar survey used to estimate high-precision chemical abundances for a
sample of over $200\,000$ Milky Way stars (Majewski et al., 2017). The
survey's simple and well-characterised selection function, based on apparent
magnitude and dereddened colour, is optimised to select red giants (Zasowski
et al., 2013, 2017). As infrared bands suffer less from dust extinction, the
red giants selected are visible across the Galactic disk. The high-precision
abundance measurements mean that we can identify monoabundance populations
with very low levels of contamination by binning the APOGEE stars in
$[\text{Fe}/\text{H}]$ and $[\alpha/\text{Fe}]$ (Bovy et al., 2016b).
Additionally, the well-characterised nature of the selection function means
that it can easily be accounted for using the method of Bovy et al. (2016a).
This makes APOGEE an ideal choice for modelling the spatial and chemical
distribution of the Milky Way's red giant population.

We obtain the chemical abundances, heliocentric distance and age of stars in
APOGEE DR16. To obtain each star's abundances, we use the calibrated ASPCAP
pipeline's (García Pérez et al., 2016) abundances of iron
$[\text{Fe}/\text{H}]$ and alpha elements $[\alpha/\text{Fe}]$, calculated
from an average of the abundances of the elements O, Mg, Si, S, and Ca, after
Bovy et al. (2016b). For each star's heliocentric distance and age, we use the
weighted_dist and age_lowess_correct estimate of AstroNN (Leung & Bovy, 2019;
Mackereth et al., 2019). weighted_dist is a weighted average of the distance
estimate from the Gaia parallax measurement of the star (Gaia Collaboration et
al., 2016) and a spectro-photometric distance estimate of the star from a
neural network trained on $265\,761$ stars surveyed in common between APOGEE
DR14 (Abolfathi et al., 2018; Holtzman et al., 2018) and Gaia DR2 (Gaia
Collaboration et al., 2018). age_lowess_correct is a measurement of the
stellar age from a neural network trained on $6676$ stars with both
spectroscopic measurement by APOGEE DR14 and asteroseismic age measurement by
the Kepler mission (Borucki et al., 2010), corrected for biases from the
neural network as described in Mackereth et al. (2017). Due to a lack of stars
with low metallicity in the training data, reliable ages aren't available for
stars with $[\text{Fe}/\text{H}]<-0.5$, so for these stars we must assume an
age distribution as described below.

The APOGEE DR16 ``statistical sample" (the sample of stars for which the
selection function can be reconstructed) contains 165,768 stars. However, this
is partly made up of dwarfs, which have higher uncertainties in their
atmospheric parameters and abundances. We thus select a subset of the
statistical sample with a calibrated ASPCAP surface gravity of $1\leq\log
g<3$. We additionally select only stars with fractional uncertainty in
heliocentric distance $D$ of less than 0.5. To restrict our sample to the
Milky Way's disk, we select stars with Galactocentric radii $R$ between 4 and
12 and height above the disk $z$ of -5 to 5. This gives us a sample of
$80\,958$ red giants.

###  2.2 Density Modelling of Red Giants across the Galaxy

To calculate the distribution of red giant stars in the Milky Way from the
APOGEE data we use the method of Bovy et al. (2016a), as this accounts for
both dust and the survey selection function simultaneously. In brief, assuming
the stars observed in a survey are distributed independently in a space of
some observables $O$ (for example position, colour and magnitude, chemical
abundances), then the positions of $N$ stars in the space of observables
$O_{1},\ldots,O_{N}$ is a realisation of an inhomogeneous Poisson point
process. This process is a random distribution of points defined by a rate
function $\lambda(O)$, such that the number of points in a given volume $V$ in
the space of the observables is a Poisson random variable with mean and
variance equal to the integral of the rate function over that volume,
$\int_{V}\lambda(O)\odif{O}$. It follows that the probability of finding a
point (i.e. an observed star) with observables in the infinitesimal volume
$\fdif{O}$ is given by $\lambda(O)\fdif{O}$, and the total number of points
(i.e. stars observed) is a Poisson random variable with mean and variance
$\Lambda=\int\lambda(O)\odif{O}$. Since the rate function is equal to the rate
of occurrence of observed stars, it can account for both the underlying true
density of stars, as well as the effect of the survey selection function and
dust which prevent all extant stars from being observed. If the rate function
is modelled with $\lambda(O\mid\theta)$, where $\theta$ parameterises the
model, then the likelihood of the model given $N$ observed stars
$O_{1},\ldots,O_{N}$ is given by

$\ln\mathcal{L}(\theta)=\sum_{i}\ln\lambda(O_{i}\mid\theta)-\int\lambda(O\mid\theta)\odif{O}\,.$
(1)

APOGEE has a selection function based on bins in dereddened colour and
apparent magnitude, so Bovy et al. (2016a) define the effective selection
function $\mathfrak{S}$, a convenient quantity equal to the fraction of a
population's stars at each heliocentric distance $D$ that will be
spectroscopically observed in each of APOGEE's fields. This is calculated for
each field by placing a tracer sample of stars in the field at that distance,
and calculating the fraction that would be observed, given a dust map and the
survey selection function in dereddened colour and apparent magnitude. We
evaluate this for each monoabundance population in each field at a range of
heliocentric distances $D$ and ages $\tau$. This allows us to treat stars as
only having three observables: the field they appear in, their distance from
the Sun $D$, and their age $\tau$. We used the effective selection function
implementation in the apogee222https://github.com/jobovy/apogee package,
described in Bovy et al. (2016a). We obtained the tracer population by
sampling PARSEC stellar model isochrones (Bressan et al., 2012; Marigo et al.,
2017) at a range of ages with a Kroupa initial mass function with a minimum
mass of $0.08M_{\odot}$ (Kroupa, 2001), cut to $1\leq\log g<3$ to match the
APOGEE red giant sample to which it was being fit. For a dust map we used a
combination of Drimmel et al. (2003), Marshall et al. (2006) and Green et al.
(2019), combined with the package mwdust333https://github.com/jobovy/mwdust,
also described in Bovy et al. (2016a).

Following Bovy et al. (2016b), we separate our red giant sample into
monoabundance populations by binning the stars in $[\text{Fe}/\text{H}]$ and
$[\alpha/\text{Fe}]$, then fit a separate number density model
$n_{\text{giants}}$ to each monoabundance population. The density model we
chose to fit to each monoabundance population was a simple axisymmetric
exponential in Galactocentric radius $R$ and height above the plane of the
disk $z$,

$n_{\text{giants}}(R,z\mid\text{logA},a_{R},a_{z})=\exp(\text{logA}-a_{R}(R-R_{0})-a_{z}\lvert
z\rvert)\,\unit{\per\kilo\cubed}\,,$ (2)

parameterised by an amplitude logA and two scale parameters $a_{R}$ and
$a_{z}$. $R_{0}=\qty{8.1}{\kilo}$ is the radial distance of the Sun from the
Galactic centre (GRAVITY Collaboration et al., 2018). Note that both $R$ and
$z$ are functions of our observables, the angular coordinates of the field
pointing and distance from the Sun $D$. We simultaneously fit the age
distribution of each monoabundance population, assumed to be a normal
distribution with mean $\tau_{0}$ and variance $1/\omega$,

$g(\tau\mid\tau_{0},\omega)=\sqrt{\frac{\omega}{2\pi}}\exp\left(-\frac{\omega}{2}(\tau-\tau_{0})^{2}\right)\,$
(3)

This form is justified by the simple narrow, monomodal age distributions for
monoabundance populations found by Lian et al. (2022), and the fact that the
effective selection function and IMF are sufficiently well-behaved that small
changes in the age distribution will not significantly change our results.

To calculate the rate function $\lambda(\text{field},D,\tau)$ for each
monoabundance population, the number density of red giants needs to be
multiplied a Jacobian factor $\lvert
J(\text{field},D)\rvert=\Omega_{\text{field}}D^{2}$ to convert it from a
density in Cartesian coordinates to a density in $D$. Multiplying by the age
distribution gives the combined underlying distribution in field, distance and
age. The observed distribution is then found by multiplying this underlying
distribution by the effective selection function
$\mathfrak{S}(\text{field},D,\tau)$. Thus the rate function for each
monoabundance population is given by

$\begin{split}\lambda(\text{field},D,\tau\mid\text{logA},a_{R},a_{z},\tau_{0},\omega)&=n_{\text{giants}}(R,z\mid\text{logA},a_{R},a_{z})\cdot
g(\tau\mid\tau_{0},\omega)\\\ &\quad\cdot\lvert
J(\text{field},D)\rvert\cdot\mathfrak{S}(\text{field},D,\tau)\,.\end{split}$
(4)

This particular form for the density profile has the advantage that the
Poisson point process likelihood takes the tractable form

$\begin{split}\ln\mathcal{L}(\text{logA},a_{R},a_{z},\tau_{0},\omega)=\text{const}&+N\left(\text{logA}-a_{R}\langle
R-R_{0}\rangle-a_{z}\langle\lvert
z\rvert\rangle+\tfrac{1}{2}\ln\omega-\frac{\omega}{2}\left(\langle\tau^{2}\rangle-2\tau_{0}\langle\tau\rangle+\tau_{0}^{2}\right)\right)\\\
&-\sum_{\text{field}}\int\int\lambda(\text{field},D\mid\text{logA},a_{R},a_{z},\tau_{0},\omega)\odif{D}\odif{\tau}\,,\end{split}$
(5)

where the sum over every star in the monoabundance population
$\sum_{i}\ln\lambda(O_{i}\mid\theta)$ in Eq. 1 is reduced to a linear
combination of the parameters with aggregates of the data in the monoabundance
population being fitted: $N$ as the number of stars observed, $\langle
R\rangle$ and $\langle\lvert z\rvert\rangle$ as the mean coordinate values,
and $\langle\tau\rangle$ and $\langle\tau^{2}\rangle$ as the mean age and age
squared. To build our model of the Milky Way disk between
$R=\qty{4}{\kilo}-\qty{12}{\kilo}$ and $\lvert z\rvert=0-\qty{5}{\kilo}$, we
found a best-fit model for each monoabundance population by maximising this
likelihood with respect to logA, $a_{R}$, $a_{z}$, $\tau_{0}$, and $\omega$,
visually checking each fit to ensure the global maximum had been found. The
results are plotted in Figure 1.

Figure 1: The best-fit values of the density modelling parameters for each
monoabundance population with 20 or more observed stars, and
$\rho_{\text{sm}}/n_{\text{giants}}$, the ratio between the sine morte stellar
mass density and the red giant number density, explained in section 2.3.

Reliable ages are not available for stars $[\text{Fe}/\text{H}]<-0.5$ in
AstroNN, so for these populations we fit only for logA, $a_{R}$ and $a_{z}$.
We use an effective selection function dependent only on field and $D$,
calculated assuming a uniform age distribution. While a uniform age
distribution is an inaccurate description of a monoabundance population, it is
a non-informative assumption that ensures that the effective selection
function is never zero where a star may actually be observed. On testing, we
found that the effective selection function did not vary strongly with the age
distribution assumed. Even then, as discussed in the section below, we do not
expect stars with metallicity $[\text{Fe}/\text{H}]<-0.5$ to contribute a
significant number of ISOs, making our conclusions independent of the age
distribution assumed for these stars.

The two main distinct chemodynamical populations of the Milky Way are clearly
shown in our model (Fig. 1). The young, high $[\text{Fe}/\text{H}]$, low
$[\alpha/\text{Fe}]$ monoabundance populations have the low vertical scale
lengths and high radial scale lengths which form the thin disk, whereas the
old low $[\text{Fe}/\text{H}]$, high $[\alpha/\text{Fe}]$ monoabundance
populations have the opposite trend in scale lengths, forming the thick disk
(Mashonkina et al., 2019).

This approach gives us simple but accurate models of the trends of the Milky
Way disk's stellar population. These results agree with the results of Bovy et
al. (2016b) which also models monoabundance populations in APOGEE, fitting
more complex models to only stars in the red clump, using their highly
consistent absolute magnitudes as an accurate distance measurement.

###  2.3 The Sine Morte Stellar Population

Having obtained a model for the distribution of red giants in the Milky Way,
we then infer the distribution of all stars throughout Galactic history. As
described in § 3.2, ISOs form mostly at the start of their parent star's life
and then outlive their parent star. Under this constraint, the population of
currently living stars is too limited to use to predict the ISO population.
Instead, we must consider what the stellar population would be at the present
time if stars did not die, but instead continued to orbit around the Milky Way
indefinitely — while being subjected to the same dynamical effects that
affected both their longer-lived companion stars and the ISOs they had
released on similar orbits. We introduce this as the sine morte444Latin for
“without death”, [sine mrte] IPA pronunciation, or ‘seen-ay mort-ay’. stellar
population.

First, we calculate the sine morte number density of stars in each
monoabundance population, from the red giant number density. For this, we use
the same PARSEC stellar model isochrones and Kroupa initial mass function with
a minimum mass of $0.08M_{\odot}$ as in § 2. Using the isochrones with
metallicities corresponding to each monoabundance population, weighted by the
fitted age distribution $g(\tau\mid\tau_{0},\omega)$ for that monoabundance
population, we calculate $N_{\text{giants}}/N_{\text{sm}}$: the fraction of
all stars ever created that are currently in the red giant phase, which is
again defined by $1\leq\log g<3$. Dividing the number density of giants by
this fraction gives us the sine morte stellar number density. This assumes
that the age distribution of each monoabundance population does not vary
significantly across the Milky Way disk, which is reasonable for the range of
$R$ and $z$ we are considering (Lian et al., 2022).

Next, we calculate the sine morte stellar mass density by multiplying the sine
morte stellar number density by the average star's initial mass, $\langle
M_{\text{int}}\rangle$. By definition, the average mass of stars in the sine
morte population is the average initial mass, which is only dependent on the
initial mass function. Thus we calculate $\langle M_{\text{int}}\rangle$ from
the same Kroupa IMF. All combined, the sine morte stellar mass density is
given by

$\rho_{\text{sm}}(\mathbf{x})=\frac{\langle
M_{\text{int}}\rangle}{N_{\text{giants}}/N_{\text{{sm}}}}\cdot
n_{\text{giants}}(\mathbf{x})\,.$ (6)

The difference between the two populations at the position of the Sun is
illustrated in Fig. 2. The ratio between the sine morte stellar mass density
and the red giant number density is plotted in the top right panel in Fig. 1.

Figure 2: Comparison of the $[\text{Fe}/\text{H}]$ distribution of red giants
and the $[\text{Fe}/\text{H}]$ distribution of the sine morte stellar
population around the Sun, both normalised. Though similar, there are subtle
differences largely caused by the fact that the distribution of ages changes
as a function of $[\text{Fe}/\text{H}]$, as shown in Fig. 1.

For the monoabundance populations with $[\text{Fe}/\text{H}]<-0.5$, without
reliable age measurements, we calculate an upper limit for the sine morte mass
density. We assume an old age distribution, with mean
$\tau_{0}=\qty{12}{\giga}$ and standard deviation
$\omega^{-\frac{1}{2}}=\qty{1}{\giga}$, which minimises
$N_{\text{giants}}/N_{\text{sm}}$. Even with this upper limit, below we find
that stars with $[\text{Fe}/\text{H}]<-0.5$ contribute a very small number of
ISOs, making our conclusions independent of this approximation.

The chemical model we use in § 3.1 to link the composition of ISOs to the
composition of stars depends only on stellar metallicity
$[\text{Fe}/\text{H}]$: so when we evaluate the model we sum the sine morte
distribution over the bins in $[\alpha/\text{Fe}]$ to get a distribution in
only $[\text{Fe}/\text{H}]$ bins. To ensure accurately fitted models, we
include only monoabundance populations with 20 or more observed stars. We then
smooth this binned distribution by taking the derivative of a spline fit to
the cumulative distribution, knotted at the edges of the bins. Plotted in Fig.
3 is the sine morte metallicity distribution
$\rho_{\text{sm}}([\text{Fe}/\text{H}])$ evaluated at the position of the Sun
($R=\qty{8.1}{\kilo}$, $z=\qty{0.021}{\kilo}$, GRAVITY Collaboration et al.
(2018); Bennett & Bovy (2019)), and integrated over the broader
$R=\qty{4}{\kilo}-\qty{12}{\kilo}$ and $\lvert z\rvert=0-\qty{5}{\kilo}$ range
of the Galactic disk we are modelling.

Figure 3: The normalised mass-weighted sine morte stellar metallicity
distribution, evaluated at the position of the Sun and integrated over the
range of the Milky Way disk we are modelling. Vertical lines show the range in
$[\text{Fe}/\text{H}]$ for which we model variations in the water mass
fraction $f_{\mathrm{H}_{2}\mathrm{O}}$ of ISOs

##  3 Predicting the Interstellar Object Distribution

In the previous section, we calculated the sine morte stellar population of
the Milky Way from the APOGEE survey. In this section we describe how to
predict an example physical property of the Galactic ISO population — the ISO
water mass fraction distribution — from this stellar population.

###  3.1 Protoplanetary Disk Model

We make the foundational assertion that all ISOs we consider form as
planetesimals in a protoplanetary disk ('Oumuamua ISSI Team et al., 2019). A
protoplanetary disk has to first order the same composition as the star it
forms around, since they both form from the same molecular cloud core. Under
this assumption, Bitsch & Battistini (2020) predict the composition of
planetesimals formed around stars of different metallicities. They do this for
stars with metallicities in the range $-0.4\leq[\text{Fe}/\text{H}]\leq 0.4$,
using the average elemental composition of stars at each value of
$[\text{Fe}/\text{H}]$ in the GALAH DR2 (Buder et al., 2018), and for
planetesimals that form both interior and exterior to the water ice line: the
inner edge of the region of the protoplanetary disk where it is cool enough to
form water ice. We assume that each star in the sine morte population produces
only ISOs with their composition set by the Bitsch & Battistini (2020) formula
for that metallicity, exterior to the water ice line. While in reality, stars
will each produce a distribution of ISOs that formed at different positions in
their protoplanetary disk and thus have a range of compositions, this
simplification of only modelling planetesimals which form exterior to the
water ice line is justified by the proportionally greater reservoir of
snowline-exterior planetesimals, and the higher efficiencies of formation
mechanisms dynamically stripping them into the interstellar population
(Fitzsimmons et al., 2023). In our Solar System, the vast majority of Oort
cloud objects are ice-rich (Meech et al., 2016); therefore both these and the
majority of ISOs produced by the Solar System must have formed outside of the
water ice line (Filacchione et al., 2022). Additionally, planetesimals beyond
the water ice line are more loosely bound to their parent stars, so will be
more easily ejected (e.g. Moro-Martín (2018)).

We focus on the mass fraction of water, $f_{\mathrm{H}_{2}\mathrm{O}}$, as
this varies significantly and decreases monotonically with
$[\text{Fe}/\text{H}]$ in the models of Bitsch & Battistini (2020). To obtain
a smooth map from $[\text{Fe}/\text{H}]$ to $f_{\mathrm{H}_{2}\mathrm{O}}$, we
fit a third-order polynomial to the water ice mass fraction data points in
figure 10 of Bitsch & Battistini (2020). For metallicities outside the range
of $-0.4\leq[\text{Fe}/\text{H}]\leq 0.4$ we assume that the relationship
between $f_{\mathrm{H}_{2}\mathrm{O}}$ and $[\text{Fe}/\text{H}]$ continues to
be monotonic, with $f_{\mathrm{H}_{2}\mathrm{O}}$ remaining high beyond the
low $[\text{Fe}/\text{H}]$ limit and remaining low beyond the high
$[\text{Fe}/\text{H}]$ limit: allowing us to track the fraction of ISOs with
high, low, and intermediate water mass fraction. This range corresponds to
$0.07\leq f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.51$ in water mass fraction.

During the writing of this paper, Cabral et al. (2023) was published,
containing an updated version of the Bitsch & Battistini (2020) chemical model
which included data from GalaH DR3 (Buder et al., 2021) and APOGEE DR17
(Abdurro'uf et al., 2022), and allowed for variation in $[\alpha/\text{Fe}]$
as well as $[\text{Fe}/\text{H}]$. They found that the water mass fraction of
planetesimals, consistent with the Bitsch & Battistini (2020) assumptions we
use, is most dependent on $[\text{Fe}/\text{H}]$. They also confirmed the
trend between $[\text{Fe}/\text{H}]$ and $f_{\mathrm{H}_{2}\mathrm{O}}$ found
in Bitsch & Battistini (2020). However, they also found that for the APOGEE
survey there was a smaller variation in $f_{\mathrm{H}_{2}\mathrm{O}}$ over
the same range in $[\text{Fe}/\text{H}]$ compared to the GalaH data used in
their work. Cabral et al. (2023) notes that this means that the trends seen
are robust. The findings of Cabral et al. (2023) do mean that the predictions
of the ISO distributions in our work may overestimate the width of the
distribution in $f_{\mathrm{H}_{2}\mathrm{O}}$.

We assume that every star produces ISOs, and the number of ISOs produced by
each star depends on its mass and metallicity. Lu et al. (2020) argue that the
mass of planet-forming material in a protoplanetary disk is proportional to
both the mass of the host star mass $M_{\ast}$ and its metal mass fraction $Z$
— well approximated by $Z_{\odot}10^{[\text{Fe}/\text{H}]}$ for small values
of Z ($Z_{\odot}=0.0153$, Caffau et al. (2011)). In the absence of confirmed
and comprehensive knowledge of ISO formation mechanisms, we use this as a
reasonable proxy for the number of ISOs produced by each star. However, the
number of ISOs produced by a star may not be simply proportional to the mass
of planet-forming material, because ISO production also requires the ejection
of planetesimals — which is dependent on system dynamical architecture. One
major ISO ejection pathway is scattering by giant planets, the occurrence of
which has its own metallicity dependence (Osborn & Bayliss, 2020).
Additionally, scattering by giant planets may form ISOs with extra
fragmentation, reweighting the number distribution towards lower masses and
increasing the number of ISOs produced from the same mass of planet-forming
material (Raymond et al., 2018b).

We therefore assume the number of ISOs produced by each star is proportional
to its mass, while incorporating into the model a power law dependence on
metallicity mass fraction: thus the number of ISOs produced is proportional to
$M_{*}\cdot 10^{\beta[\text{Fe}/\text{H}]}$. Here $\beta=1$ corresponds to the
simple assumption that the number of ISOs produced is proportional to the mass
of planet-forming material. $\beta=0$ corresponds to no metallicity dependence
at all. However, we expect $\beta>0$: planetesimals require some fraction of
dust and ice in order to exist, i.e. at the elemental level, C/N/O, Al/Si, &c.
must be present, so a minimum metallicity constraint must exist (Johnson & Li,
2012). We do not model the constant of proportionality here, as this depends
on the size distribution of ISOs, which remains observationally unconstrained
with only two ISOs ('Oumuamua ISSI Team et al., 2019; Jewitt & Seligman,
2022). We predict only the normalised distribution of ISO water mass
fractions.

###  3.2 Predicting the ISO Population From the Stellar Population

Since the majority of ISOs are expected to be ejected by dynamical processes
within several hundred Myr of a star's formation (Pfalzner & Bannister, 2019;
Fitzsimmons et al., 2023), we assume that the birth of a star and release of
its ISOs are contemporaneous, and omit accounting for any delay. As Oort cloud
comets, which also occupy an interstellar environment (Kaib & Volk, 2022), do
not exhibit substantive erosion after 4.5 Gyr (Stern & Shull, 1988), we
continue with the expectation that ISO erosional processes are broadly similar
to those of Solar System comets (e.g. Guilbert-Lepoutre et al., 2015). This
implies ISOs from stars prior to the Sun will outlive their parent stars. We
thus use the sine morte stellar distribution established in § 2.3. For
assessing the water mass fraction of the resulting population, we assume
processing in the interstellar medium has a negligible effect555While this is
valid for water, it may not hold for hypervolatiles; e.g. Seligman et al.
(2022) suggest that processing in the interstellar medium will remove CO and
C$\mathrm{O}_{2}$ relative to $\mathrm{H}_{2}$O. However, this is countered by
the CO-rich nature of 2I., and that ISOs broadly represent the water mass
fraction of their source planetesimal populations.

To link the spatial distribution of stars throughout the Galaxy to the ISOs
they produced, we need to model their combined Galactic dynamics. ISOs are not
expected to remain near to their parent stars on Gyr timescales, and stars do
die. However, it is helpful to consider what the behaviour of the ISO
population is in the case that neither of these are true. If ISOs did remain
near their parent stars and both ISOs and stars existed for an infinite length
of time, the number density distribution of ISOs with a given water mass
fraction $f_{\mathrm{H}_{2}\mathrm{O}}$ and position in the Galaxy
$\mathbf{x}$ in the present day, denoted
$n_{\text{ISO}}(\mathbf{x},f_{\mathrm{H}_{2}\mathrm{O}})$, would be equal to
the number of ISOs produced by stars currently at the same position and of the
corresponding metallicity $[\text{Fe}/\text{H}]$. Under the model of § 3.1,
the number of ISOs produced by a star of mass $M_{*}$ and metallicity
$[\text{Fe}/\text{H}]$ is proportional to $M_{*}\cdot
10^{\beta[\text{Fe}/\text{H}]}$, thus

$n_{\text{ISO}}(\mathbf{x},f_{\mathrm{H}_{2}\mathrm{O}}\mid\beta)\propto
10^{\beta[\text{Fe}/\text{H}]}\cdot\odv{[\text{Fe}/\text{H}]}{f_{\mathrm{H}_{2}\mathrm{O}}}\cdot\rho(\mathbf{x},[\text{Fe}/\text{H}])$
(7)

where $\rho(\mathbf{x},[\text{Fe}/\text{H}])$ is the mass density distribution
of stars at position $\mathbf{x}$ with the metallicity $[\text{Fe}/\text{H}]$
corresponding to the ISO water mass fraction $f_{\mathrm{H}_{2}\mathrm{O}}$,
and $\odif{[\text{Fe}/\text{H}]}/\odif{f_{\mathrm{H}_{2}\mathrm{O}}}$ is the
gradient of the relationship between $[\text{Fe}/\text{H}]$ and
$f_{\mathrm{H}_{2}\mathrm{O}}$ described in § 3.1.

The fact that stars do in fact die is then corrected for by replacing
$\rho(\mathbf{x},[\text{Fe}/\text{H}])$ in this equation with
$\rho_{\text{sm}}(\mathbf{x},[\text{Fe}/\text{H}])$, the sine morte stellar
mass density introduced in § 2.3. Correcting for the fact that ISOs disperse
from near their parent stars is more complex, so we proceed with some
simplifying assumptions. Once ejected, unless its resultant velocity exceeds
the Galactic escape velocity, an ISO will orbit the Galactic centre, as stars
do. We assume ISOs are ejected relatively slowly from their parent planetary
systems compared to the stellar velocity dispersion. This is justified
assuming ejection velocities $<\qty{10}{\kilo\per}$, the maximum ejection
velocity from a planetary system under an expected suite of scattering
mechanisms (Pfalzner & Bannister, 2019; Fitzsimmons et al., 2023), and the
velocity dispersions $\gtrsim\qty{20}{\kilo\per}$ measured in the Solar
Neighbourhood (Anguiano et al., 2020). Stars form on near-circular orbits
around the Galactic centre (Frankel et al., 2020), so a cloud of recently
ejected ISOs will all have similar orbits to their parent star: nearly
circular with similar ranges of oscillation in Galactocentric $R$ and $z$.
However, the slight differences in their orbits will give the ISOs different
orbital periods around the Galactic centre, meaning they will disperse along
their similar near-circular orbital paths. Therefore, though ISOs do not stay
near their parent star, we assume here that they only disperse in the
azimuthal direction, and remain in the same $R$ and $z$ range as their parent
star. Under this assumption equation 7 still holds if the stellar density
model is axisymmetric, depending only on $R$ and $z$, as does ours:

$n_{\text{ISO}}(R,z,f_{\mathrm{H}_{2}\mathrm{O}}\mid\beta)\propto
10^{\beta[\text{Fe}/\text{H}]}\cdot\odv{[\text{Fe}/\text{H}]}{f_{\mathrm{H}_{2}\mathrm{O}}}\cdot\rho_{\text{sm}}(R,z,[\text{Fe}/\text{H}])$
(8)

Orbits around the Galactic centre can evolve with time, due to the influence
of perturbing potentials such as spiral arms, the bar and molecular clouds.
These effects can cause both dynamical `heating', an increase in the size of
radial and vertical excursions from a circular orbit, and `migration', changes
to the radius of an orbit while it remains nearly circular (Sellwood & Binney,
2002). We introduced Eq. 8 with the assumption that an ISO will stay in the
same range of $R$ and $z$ as its parent star. Due to their azimuthal
separation, the star and ISOs will experience slightly different perturbing
potentials, causing their orbits to evolve adjacently but independently.
However, if the Galactic stellar and ISO distributions are sufficiently
axisymmetric, perturbations will consistently change together both the orbits
of stars and the orbits of a corresponding number of ISOs. Thus, Eq. 8 holds
for our model.

In this work we predict the distribution of ISOs in
$f_{\mathrm{H}_{2}\mathrm{O}}$ both at particular values of $R$ and $z$ and
integrated over the whole Milky Way. Since we do not model the total number of
ISOs, we remove the need for the constant of proportionality in Eq. 8 by
normalising each $n_{\text{ISO}}$ distribution we calculate with

$p(f_{\mathrm{H}_{2}\mathrm{O}}\mid\beta)=\frac{n_{\mathrm{ISO}}(f_{\mathrm{H}_{2}\mathrm{O}}\mid\beta)}{\int
n_{\mathrm{ISO}}(f_{\mathrm{H}_{2}\mathrm{O}}\mid\beta)\odif{f_{\mathrm{H}_{2}\mathrm{O}}}}\,.$
(9)

This gives us the distribution of ISOs within the bounds of the protoplanetary
disk chemical model: $0.07\leq f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.51$. Outside
of this range, we can calculate the fraction of ISOs with
$f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.07$ and $f_{\mathrm{H}_{2}\mathrm{O}}\geq
0.51$, by assuming that the relation between $[\text{Fe}/\text{H}]$ and
$f_{\mathrm{H}_{2}\mathrm{O}}$ remains monotonic. Thus, all stars with
$[\text{Fe}/\text{H}]\leq-0.4$ contribute ISOs with
$f_{\mathrm{H}_{2}\mathrm{O}}\geq 0.51$, and all stars with
$[\text{Fe}/\text{H}]\geq 0.4$ contribute ISOs with
$f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.07$.

##  4 Results

In this section we demonstrate the prediction framework of section 3 by making
two different predictions. We demonstrate two different example values for a
stellar metallicity dependence for ISO production, $\beta$, with $\beta=1$ for
our principal prediction and $\beta=0$ as an alternate prediction.

###  4.1 Principal Prediction: $\beta=1$

First, we predict the distribution of ISOs assuming that the number produced
by each star is proportional to the star's metal mass fraction $Z$, by setting
$\beta=1$ in equation 8. As described in section 3.1, in the absence of
concrete knowledge of ISO formation mechanisms this is a reasonable value to
assume, and thus we consider this our principal prediction.

ISO $f_{\mathrm{H}_{2}\mathrm{O}}$ range | Fraction of ISOs around Sun | Fraction of ISOs in Milky Way Disk
---|---|---
$f_{\mathrm{H}_{2}\mathrm{O}}<0.07$ | 0.017 | 0.024
$0.07\leq f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.51$ | 0.955 | 0.915
$0.51<f_{\mathrm{H}_{2}\mathrm{O}}$ | 0.027 | 0.061
Table 1: Primary prediction for the fraction of ISOs in each
$f_{\mathrm{H}_{2}\mathrm{O}}$ range evaluated at the position of the Sun and
integrated over whole Milky Way, with $\beta=1$.

Table 1 lists the the fraction of ISOs within and outside either end of the
Bitsch & Battistini (2020) protoplanetary disk chemical model
$f_{\mathrm{H}_{2}\mathrm{O}}$ range, $0.07\leq
f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.51$. We assess both the distribution of
ISOs at the position of the Sun, at $R=\qty{8.1}{\kilo}$,
$z=\qty{0.021}{\kilo}$ (GRAVITY Collaboration et al., 2018; Bennett & Bovy,
2019), and the distribution of ISOs integrated over the region of the Milky
Way disk we are modelling. In both cases, the vast majority of ISOs lie within
the range of the model. The significant mass of stars beyond the lower
$[\text{Fe}/\text{H}]$ limit of the $\rho_{\text{sm}}([\text{Fe}/\text{H}])$
distribution in Fig. 3, both at the position of the Sun and over the whole
Milky Way disk, do not contribute significantly to the fractions of ISOs in
the high $f_{\mathrm{H}_{2}\mathrm{O}}$ range. This is because of the
exponential dependence of the number of ISOs produced by each star on
$[\text{Fe}/\text{H}]$ in Eq. 8.

Figure 4: Primary prediction for the distribution of ISO water mass fractions,
evaluated at the position of the Sun and integrated over the Milky Way disk,
for $\beta=1$.

Within the range of the chemical model, we can plot the distribution function
of the ISO water mass fractions. Figure 4 shows both the population of ISOs
around the Sun and over the whole Milky Way. The different shapes of the two
metallicity distribution functions in Figure 3 is also apparent here, as the
wider Milky Way metallicity distribution function results in a wider ISO water
mass fraction distribution.

The distributions of ISOs around the Sun and averaged over the Milky Way disk
are remarkably similar. We explore this in Figure 5, through the sine morte
stellar mass $[\text{Fe}/\text{H}]$ distribution and the ISO
$f_{\mathrm{H}_{2}\mathrm{O}}$ distribution within the range of the chemical
model at a range of values of $R$, integrated over $z$. Also plotted are the
median values of $[\text{Fe}/\text{H}]$ and $f_{\mathrm{H}_{2}\mathrm{O}}$ at
each value of $R$. Clear in the left-hand panel of Fig. 5 is the well-studied
Galactic metallicity gradient (Cheng et al., 2012) — but additionally in the
right-hand panel is a corresponding gradient in ISO water mass fraction. Since
the composition of ISOs depends on the chemical makeup of the stars that they
form around, we expect trends in the chemical abundances of stars to be
accompanied by equivalent trends in the compositions of ISOs. Figure 5 also
shows why the Solar neighbourhood distributions are similar to the whole-
Galaxy integrated distributions in Figures 3 and 4: The Solar neighbourhood
happens to be at an intermediate value of $R$ (8.1), where both the stellar
$[\text{Fe}/\text{H}]$ distribution and therefore the ISO
$f_{\mathrm{H}_{2}\mathrm{O}}$ distributions are approximately midway between
the high and low extremes.

Figure 5: Sine morte stellar mass distributions (left) and ISO water mass
fraction distributions (right), integrated over $z$, at each distance to the
Galactic Centre $R$. Both distributions are normalised such that at each value
of $R$ the integral over $[\text{Fe}/\text{H}]$ or
$f_{\mathrm{H}_{2}\mathrm{O}}$ is unity. Dashed lines show the median value of
each distribution at each $R$.

In these plots, we have normalised the distributions in $[\text{Fe}/\text{H}]$
and $f_{\mathrm{H}_{2}\mathrm{O}}$ such that at each value of $R$ the integral
over $[\text{Fe}/\text{H}]$ or $f_{\mathrm{H}_{2}\mathrm{O}}$ is unity.
However, it is worth noting that our model implies the densities of both the
stellar and ISO populations will decrease with distance from the Galactic
centre. The exponential stellar density profile is well established (Jurić et
al., 2008). Here we predict that the Galactic ISO population density profile
will decrease faster than that of the stars: due to the metallicity dependence
of the number of ISOs produced by each star, the higher-metallicity stars in
the inner disk will produce more ISOs per unit of stellar mass than the lower-
metallicity stars of the outer disk.

###  4.2 Alternate Prediction: $\beta=0$

Setting $\beta=1$ is just a choice in a basic and observationally
unconstrained model. Therefore, we explore how changing its value affects the
predicted ISO distribution. Lintott et al. (2022) predicted the ISO
populations of Galaxies from the EAGLE hydrodynamical simulation using the
same protoplanetary disk chemical model to map stellar metallicities
$[\text{Fe}/\text{H}]$ to ISO water mass fractions
$f_{\mathrm{H}_{2}\mathrm{O}}$. However, Lintott et al. (2022) assumed that
the number of ISOs produced by each star was independent of the star's
metallicity. This is equivalent to setting $\beta=0$ in Eq. 8 in this work,
and we use this to make our alternate prediction. The resulting water mass
fraction distributions for ISOs at the position of the Sun and integrated over
the whole Milky Way are plotted in Figure 6 and tabulated in Table 2.

Figure 6: Alternate prediction for the distribution of ISO water mass fractions, evaluated at the position of the Sun and integrated over the Milky Way disk with $\beta=0$: assuming the number of ISOs produced by each star is independent of the star's metallicity. ISO $f_{\mathrm{H}_{2}\mathrm{O}}$ range | Fraction of ISOs around Sun | Fraction of ISOs in Milky Way Disk
---|---|---
$f_{\mathrm{H}_{2}\mathrm{O}}<0.07$ | 0.007 | 0.009
$0.07\leq f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.51$ | 0.903 | 0.798
$0.51<f_{\mathrm{H}_{2}\mathrm{O}}$ | 0.090 | 0.194
Table 2: Alternate prediction for fraction of ISOs in each
$f_{\mathrm{H}_{2}\mathrm{O}}$ range, evaluated at the position of the Sun and
integrated over the whole Milky Way with $\beta=0$: assuming the number of
ISOs produced by each star is independent of the star's metallicity.

Assuming $\beta=0$ causes the ISO distribution to be more weighted towards a
higher water mass fraction. These high water mass fraction ISOs come from low
metallicity stars, which in this scenario no longer have their contributions
to the ISO population suppressed by the exponential dependence on
$[\text{Fe}/\text{H}]$. This means that the low metallicity tails of the
$\rho_{\text{sm}}([\text{Fe}/\text{H}])$ distributions in Figure 3 now
contribute a significant number of ISOs with
$0.51<f_{\mathrm{H}_{2}\mathrm{O}}$: almost 20% when averaged over the Milky
Way disk. A comparison of these results to those of Lintott et al. (2022) is
discussed in section 6.1.

##  5 Inference Framework

###  5.1 Bayesian Inference

In this work so far we have detailed a method of predicting the distribution
of ISOs around the Sun, and demonstrated that this method can be used with
varying models to get different predictions. These predictions can be compared
to an observed sample of ISOs in order to make inferences about the models
which went into producing them. In this section, we detail how to do this in a
Bayesian manner.

The predictions we make in this work are normalised ISO water mass fraction
distributions, $p(f_{\mathrm{H}_{2}\mathrm{O}}\mid\beta)$, parameterised by
the ISO production power law slope $\beta$. This however could be generalised
to a distribution in any set of observable properties $\mathbf{O}$, and
parameterised by any set of parameters $\theta$. These distributions are also
probability density functions for the value of $\mathbf{O}$ of a randomly
selected ISO, so $p(\mathbf{O}\mid\theta)\fdif{\mathbf{O}}$ is the probability
of an ISO having properties in the infinitesimal volume $\fdif{\mathbf{O}}$.
Extending this to multiple ISOs, a sample of ISOs with properties
${\mathbf{O}}_{1},\ldots,{\mathbf{O}}_{N}$ will have likelihood

$p({\mathbf{O}}_{1},\ldots,{\mathbf{O}}_{N}\mid\theta)=\prod_{i}p({\mathbf{O}}_{i}\mid\theta)$

assuming ISOs arrive independently with independent compositions. This is a
valid assumption if the ISOs we observe are samples from a large Galactic
population with contributions from stars across the Galaxy and cosmic time.
This is because the probability that two randomly selected ISOs have the same
parent star is approximately equal to one over the total number of stars which
could contribute ISOs to the population around the Sun. The ISO sample
likelihood can then be combined with priors to form a posterior distribution
on the model parameters:

$p(\theta\mid{\mathbf{O}}_{1},\ldots,{\mathbf{O}}_{N})\propto
p(\theta)\cdot\prod_{i}p({\mathbf{O}}_{i}\mid\theta)$ (10)

This posterior can be used to calculate estimates and confidence intervals for
the values of the parameters $\theta$. An example of this calculation is given
in section 6.3.

###  5.2 Additional Properties of ISOs

As with all minor planets, there are many observable quantities for ISOs that
could be used in this framework in the set of observable properties
$\mathbf{O}$. If the distribution of ISOs in these properties can be
predicted, the inference method of § 5.1 can be used to compare the predicted
distribution to the observed distribution in these properties, to make
inferences about the models used. If multiple properties are included in
$\mathbf{O}$, then the joint distribution of ISOs in these properties can be
predicted and used in inference. Including different ISO properties in the
framework will allow inferences to be made about different processes affecting
the ISO population. We briefly consider several such properties here.

The Milky Way stellar population is broadly divided between different chemo-
dynamical populations: the thin disk (high $[\text{Fe}/\text{H}]$, low
velocity dispersion), thick disk (low $[\text{Fe}/\text{H}]$, high velocity
dispersion), and the halo (very low $[\text{Fe}/\text{H}]$, radial orbits)
(Recio-Blanco et al., 2014; Horta et al., 2023). Since the composition of an
ISO depends on the metallicity of the star it formed around, each chemo-
dynamical stellar population will contribute ISOs to the Galactic population
with a distinct joint distribution in both composition and velocity.
Therefore, including the velocities of ISOs in this framework could be used to
tie ISOs to the chemo-dynamical stellar populations their parent star belongs
to (Eubanks et al., 2021).

As an alternative measurement of composition to water mass fraction, the
carbon-to-oxygen ratio of a cometary ISO can be estimated from its coma. This
contains information about an ISO's formation location in the protoplanetary
disk relative to the H2O, CO2 and CO ice lines; modelling suggests it would be
a useful measure for future ISOs (Seligman et al., 2022).

Further additions to our framework could include predictions of the size
distribution and aspect ratio distribution of ISOs, which would be especially
interesting due to 1I/‘Oumuamua's extreme shape. The prediction of these
distributions would depend on the stellar population via models of
planetesimal formation which could link the size and shape distribution of
planetesimals to the properties of their natal star. Finally, including the
binarity rate of ISOs would test models of the applicability of formation
mechanism of planetesimals to other protoplanetary disks, as the compact
binarity of trans-Neptunian objects does for the Solar System protoplanetary
disk (e.g. Nimmo et al., 2018). This could also constrain the ejection
mechanisms of ISOs: loosely-bound wide planetesimal binaries would not survive
scattering by a giant planet, but would survive a more gentle gravitational
interaction with a stellar flyby.

###  5.3 Generalisation to Other Galactic Populations

Our method could easily be generalised to any Galaxy-wide population with a
dependence on the properties of stars, simply by replacing the ``ISO recipe''
of section 3.

For example, the distribution of planets through the Galaxy could be predicted
by substituting a model of the occurrence rate of planets as a function of
their host star's metallicity. This is based on the planet-metallicity
correlation (Fischer & Valenti, 2005; Osborn & Bayliss, 2020): as noted
earlier that dust is necessary for planets, stars with higher metallicity are
more likely to host planets. Here, the Milky Way metallicity gradient means
that planets are more common towards the Galactic centre. This is a testable
prediction with microlensing surveys such as OGLE (Udalski et al., 2015),
which continue to find exoplanets between the Solar System and the Galactic
centre (Tsapras, 2018) with distances estimated from followup observations
(Vandorou et al., 2023). If ISOs seed planet formation as hypothesised by
Pfalzner & Bannister (2019), the ISO gradient we predict in § 4.1 could also
produce a signature in the planetary population.

##  6 Discussion

###  6.1 Comparison to Previous Work

As described in section 3.1, Cabral et al. (2023) updates the protoplanetary
disk chemical model of Bitsch & Battistini (2020) that we use, and note that
the trend of planetesimal water mass fraction decreasing with stellar
metallicity is robust. They do however find that for the APOGEE survey, there
was a smaller variation in $f_{\mathrm{H}_{2}\mathrm{O}}$ over the same range
in $[\text{Fe}/\text{H}]$ compared to the GalaH data used in Bitsch &
Battistini (2020). If the true variation in $f_{\mathrm{H}_{2}\mathrm{O}}$ is
smaller than in the model of this work, our predictions may overestimate the
width of the ISO water mass fraction distribution.

Lintott et al. (2022) made a prediction of the ISO population of a simulated
Milky Way Galaxy from the EAGLE hydrodynamical cosmological simulation (Schaye
et al., 2015), using a model equivalent to that of this work with $\beta=0$.
Whereas we predict a single-peaked $f_{\mathrm{H}_{2}\mathrm{O}}$
distribution, they predicted an ISO distribution with a significant number of
ISOs with water mass fraction both below and above the range of the
protoplanetary disk model, which they interpret as a bimodal distribution in
ISO composition.

There are expected reasons for the difference between the prediction here,
based on the observed Milky Way stellar population, and the prediction by
Lintott et al. (2022) based on the simulated EAGLE Galaxy. The EAGLE Galaxy
has a much wider $[\text{Fe}/\text{H}]$ distribution than the Milky Way, with
many more stars outside of the $[\text{Fe}/\text{H}]$ range of the
protoplanetary disk model. As a smoothed particle hydrodynamics simulation,
EAGLE is susceptible to producing Galaxies with $[\text{Fe}/\text{H}]$
distributions wider than those observed in nature, due to underestimating
metal mixing between particles (Wiersma et al., 2009). Therefore the results
of this work, based on the observed stellar population of the Milky Way,
should give a much more accurate prediction for the Milky Way's population of
ISOs.

###  6.2 Distinguishing Local and Galactic Populations of ISOs

The results of § 4.1 show that the well-studied metallicity gradient of the
Milky Way has a corresponding ISO composition gradient, with ISOs generally
having a higher water mass fraction at larger Galactocentric radii. Although
we can only observe the compositions of ISOs which pass through the inner
Solar System, it is still instructive to model how the distribution of ISOs
varies across a wider portion of the Galactic disk. This is because we have
made assumptions about the Galactic dynamics of ISOs, under which the
distribution of ISOs at a point in the Galaxy corresponds to the distribution
of stars at that same point. If these assumptions break down, then the
particular way in which they break down will affect the population of ISOs
detectable in the Solar System in a related, calculable way.

For example, radial migration, caused by the non-axisymmetric potential of
spiral arms, flattens the the Milky Way metallicity gradient by blurring the
metallicity distribution in the radial direction (Vickers et al., 2021). This
widens the stellar metallicity distribution around the Sun as stars migrate in
from the metal-poor outer disk and out from metal-rich inner disk. However,
ISOs may undergo less radial migration than stars, due to the random motion
given to them by their ejection from their home planetary system (Daniel &
Wyse, 2015). The stars currently in the Solar neighbourhood may thus have a
wider range of Galactocentric radii of origin than the ISOs. This would make
the distribution of observable ISO compositions narrower than would be
predicted from the distribution of stars.

Additionally, the low velocity of 1I/‘Oumuamua relative to the local standard
of rest implies that it — and therefore a large fraction of observable ISOs —
could come from local star-forming regions (e.g. Hallatt & Wiegert, 2020).
This is also testable: their compositions will match those predicted from the
metallicities of the nearest star forming regions. In addition, if the
velocity distribution of ISOs is included in future work as described in 5.2,
it may be possible to trace individual ISOs back to the star forming regions
they came from by matching them up with both composition and velocity. This
would be hugely advantageous for studying planetesimal formation, as it would
allow us to pair the properties of detected ISOs directly with observations of
their parent planetary systems.

###  6.3 An Estimation of the ISO Production Metallicity Dependence

The two different predictions made in § 4 demonstrate that the Galactic
population of ISOs is sensitive to small changes to the processes that affect
their formation and evolution. This means that if models of these processes
can be combined to make accurate predictions of the ISO population, then the
framework described in § 5 can be used to make inferences about those
processes.

In particular, the predictions of sections 4.1 and 4.2 show that the ISO
distribution around the Sun is sensitive to the metallicity dependence of the
number of ISOs produced by each star, $\beta$. As outlined in section 5.1,
these two predictions can then be compared to a sample of ISOs. In our example
physical property, this is done by calculating the likelihood of each of these
predictions producing the observed distribution of water mass fractions, such
as for the ISOs expected to be found by the Legacy Survey of Space and Time
(LSST) of the Vera C. Rubin Observatory. For a more general result, this
likelihood can be combined with a prior on $\beta$ and Bayes' theorem to
calculate a posterior distribution for $\beta$.

Though the work of this paper has been carried out in expectation of a larger
sample of ISOs being known in the future, we do already have a preliminary
estimate of the ISO water mass fraction distribution. Due to the unknown
composition of 1I/‘Oumuamua, we can only estimate the water mass fraction of
2I/Borisov. We adopt a value of $f_{\mathrm{H}_{2}\mathrm{O}}=0.3$, after that
of Seligman et al. (2022) inferred from the production rates of 2I's coma.
However, this value is more useful for demonstrative purposes than as a finely
constrained estimate; estimating a comet's bulk composition from the
composition of its coma is challenging. For instance, Seligman et al. (2022)
note that the compositions of interstellar comets calculated from production
rates are affected by preferential desorbtion of CO and $\textrm{CO}_{2}$
relative to $\mathrm{H}_{2}$O.

With this observed distribution of $f_{\mathrm{H}_{2}\mathrm{O}}$ we can then
use the Bayesian framework of § 5.1 to calculate a posterior distribution for
$\beta$. We use our model for the distribution of ISOs around the Sun, and
taking a uniform prior on $\beta$ means that the posterior is simply
proportional to the likelihood, equal to the value of the
$f_{\mathrm{H}_{2}\mathrm{O}}$ distribution. For
$f_{\mathrm{H}_{2}\mathrm{O}}=0.3$ this posterior is maximised by a value of
$\beta=1.33$. The distribution of ISOs at this value is listed and plotted in
Table 3 and Figure 7. The symmetric 90% confidence limit for this estimation
of $\beta$ is $\beta\in(-1.3,7.2)$. This is of course a wide interval
containing physically implausible values. The physically implausible values
could be removed from the confidence interval by using a physically motivated
prior on $\beta$ — but with one known value of $f_{\mathrm{H}_{2}\mathrm{O}}$
this would make the posterior dominated by the prior. Using the sample of ISOs
that the LSST is expected to find will much better constrain the posterior for
$\beta$ and other parameters used in models with this framework.

Figure 7: Distribution of ISO water mass fractions, evaluated at the position of the Sun and integrated over the Milky Way disk, with maximum a posteriori estimate of $\beta=1.33$. At $f_{\mathrm{H}_{2}\mathrm{O}}=0.3$ a vertical line marks the measured water mass fraction of 2I/Borisov.  ISO $f_{\mathrm{H}_{2}\mathrm{O}}$ range | Fraction of ISOs around Sun | Fraction of ISOs in Milky Way Disk
---|---|---
$f_{\mathrm{H}_{2}\mathrm{O}}<0.07$ | 0.022 | 0.031
$0.07\leq f_{\mathrm{H}_{2}\mathrm{O}}\leq 0.51$ | 0.960 | 0.930
$0.51<f_{\mathrm{H}_{2}\mathrm{O}}$ | 0.018 | 0.040
Table 3: Predicted fraction of ISOs in each $f_{\mathrm{H}_{2}\mathrm{O}}$
range, evaluated at the position of the Sun and integrated over the Milky Way
dist, with maximum a posteriori estimate of $\beta=1.33$.

###  6.4 Selection Effects

It should be noted that the predictions in this work ignore some specific
effects which will influence the population of ISOs we expect to observe from
the Earth. The predictions in § 4 are the distribution of ISOs in a smooth
Galaxy-wide distribution, evaluated at the location of the Solar System.
Gravitational focussing by the Sun will increase the density of ISOs in the
inner solar system in a velocity-dependent manner (e.g. Engelhardt et al.,
2017; Forbes & Loeb, 2019; Dehnen & Hands, 2022), and therefore also in a
composition-dependent manner. This is due to the fact that we expect
compositionally and dynamically distinct populations of ISOs to come from the
chemo-dynamically distinct stellar populations in the Milky Way thin and thick
disks (c.f. Eubanks et al., 2021). This could handily be modelled when
incorporating ISO velocities into the predictions of this framework.

Additionally, the Pan-STARRS near-Earth object survey (Chambers et al., 2016)
that detected 1I/‘Oumuamua, the observations of amateur astronomers such as
Gennadiy Borisov who discovered 2I/Borisov, and the LSST which will discover
tens more ISOs: all have highly non-trivial selection functions. Since in
order to be discovered an ISO needs to be detected in multiple observations
which can be linked as the same object (e.g. Meech et al., 2017; Schwamb et
al., 2023), these selection functions are dependent on ISO size, approach
velocity, perihelion, and composition. Composition has a direct link to the
detectability of ISOs, as 2I-sized and cometary ISOs will be more likely to be
detected, as they will form a coma as they approach the Sun. These selection
effects will need to be accurately accounted for, in order for the Bayesian
framework to produce accurate inferences about the processes affecting the
Galactic ISO population.

##  7 Conclusion

In advance of the Vera C. Rubin Observatory Legacy Survey of Space and Time
(LSST), we lay out a framework to predict the Galactic distribution of ISOs,
using the stellar population of the Milky Way. Using the method of Bovy et al.
(2016a), we fit simple density models to a sample of red giants in APOGEE
binned in $[\text{Fe}/\text{H}]$ and $[\alpha/\text{Fe}]$, and use these to
evaluate the ``sine morte" metallicity distribution of stars throughout the
Galaxy's integrated history, across the Galactic disk. Under the assumption
that the spatial distribution of a population of ISOs will be the same as the
sine morte distribution of stars which formed them, we use the protoplanetary
disk model of Bitsch & Battistini (2020) to map the metallicity distribution
of stars to the distribution of ISO water mass fractions. Localising our model
to the Solar neighbourhood, we predict that 95% of ISOs around the Sun have
water mass fraction $f_{\mathrm{H}_{2}\mathrm{O}}$ between 0.07 and 0.51, with
a peak around 0.35.

By considering the distribution of ISOs over the Galactic disk, we show that
the well-studied Milky Way metallicity gradient has an equivalent gradient in
ISO composition, with the median ISO water mass fraction increasing with
distance from the Galactic Centre as the median stellar metallicity decreases.
This causes the ISO water mass fraction distribution averaged over the Milky
Way disk to be wider that the distribution around the Sun. Since we also
predict higher-metallicity stars produce more ISOs than lower-metallicity
stars, the Milky Way metallicity gradient implies that the radial ISO density
profile is steeper than the exponential radial stellar density profile, making
ISOs much more common in the inner Galactic disk than in the outer disk.

We also set out a Bayesian inference framework, which can compare predictions
of the ISO distribution to the sample observed by the LSST, in order to make
inferences about the many astrophysical processes which influence the ISO
population. To demonstrate its use, we use the composition measurement of
2I/Borisov to calculate a maximum a posteriori estimate for the power law
slope of the ISO production metallicity dependence, $\beta$, to be 1.33, with
a symmetric 90% confidence interval of $(-1.3,7.2)$.

Encoded in the population of ISOs we observe is a wealth of information about
processes through Galactic history, from the evolution of the Milky Way to
planet formation. The framework set out in this work is a novel approach that
will allow us to appreciate how these treasures can further our understanding
on both planetary and Galactic scales.

MJH acknowledges support from the Science and Technology Facilities Council through grant ST/W507726/1. MJH also thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation; his participation in the program has benefitted from this work. MTB appreciates support by the Rutherford Discovery Fellowships from New Zealand Government funding, administered by the Royal Society Te Apārangi. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

## References

  * Abdurro'uf et al. (2022) Abdurro'uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35, doi: 10.3847/1538-4365/ac4414
  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42, doi: 10.3847/1538-4365/aa9e8a
  * Anguiano et al. (2020) Anguiano, B., Majewski, S. R., Hayes, C. R., et al. 2020, AJ, 160, 43, doi: 10.3847/1538-3881/ab9813
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Bannister et al. (2017) Bannister, M. T., Schwamb, M. E., Fraser, W. C., et al. 2017, ApJ, 851, L38, doi: 10.3847/2041-8213/aaa07c
  * Bennett & Bovy (2019) Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417, doi: 10.1093/mnras/sty2813
  * Bitsch & Battistini (2020) Bitsch, B., & Battistini, C. 2020, A&A, 633, A10, doi: 10.1051/0004-6361/201936463
  * 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
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Bovy et al. (2016a) Bovy, J., Rix, H.-W., Green, G. M., Schlafly, E. F., & Finkbeiner, D. P. 2016a, ApJ, 818, 130, doi: 10.3847/0004-637X/818/2/130
  * Bovy et al. (2016b) Bovy, J., Rix, H.-W., Schlafly, E. F., et al. 2016b, ApJ, 823, 30, doi: 10.3847/0004-637X/823/1/30
  * Brasser et al. (2006) Brasser, R., Duncan, M. J., & Levison, H. F. 2006, Icarus, 184, 59, doi: 10.1016/j.icarus.2006.04.010
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Buder et al. (2018) Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513, doi: 10.1093/mnras/sty1281
  * Buder et al. (2021) Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150, doi: 10.1093/mnras/stab1242
  * Cabral et al. (2023) Cabral, N., Guilbert-Lepoutre, A., Bitsch, B., Lagarde, N., & Diakite, S. 2023, A&A, 673, A117, doi: 10.1051/0004-6361/202243882
  * Caffau et al. (2011) Caffau, E., Ludwig, H. G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255, doi: 10.1007/s11207-010-9541-4
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560, doi: 10.48550/arXiv.1612.05560
  * Cheng et al. (2012) Cheng, J. Y., Rockosi, C. M., Morrison, H. L., et al. 2012, ApJ, 746, 149, doi: 10.1088/0004-637X/746/2/149
  * 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
  * Daniel & Wyse (2015) Daniel, K. J., & Wyse, R. F. G. 2015, MNRAS, 447, 3576, doi: 10.1093/mnras/stu2683
  * Dehnen & Hands (2022) Dehnen, W., & Hands, T. O. 2022, MNRAS, 512, 4062, doi: 10.1093/mnras/stab3670
  * Do et al. (2018) Do, A., Tucker, M. A., & Tonry, J. 2018, ApJ, 855, L10, doi: 10.3847/2041-8213/aaae67
  * Drimmel et al. (2003) Drimmel, R., Cabrera-Lavers, A., & López-Corredoira, M. 2003, A&A, 409, 205, doi: 10.1051/0004-6361:20031070
  * Engelhardt et al. (2017) Engelhardt, T., Jedicke, R., Vereš, P., et al. 2017, AJ, 153, 133, doi: 10.3847/1538-3881/aa5c8a
  * Eubanks et al. (2021) Eubanks, T. M., Hein, A. M., Lingam, M., et al. 2021, arXiv e-prints, arXiv:2103.03289, doi: 10.48550/arXiv.2103.03289
  * Farnocchia et al. (2023) Farnocchia, D., Seligman, D. Z., Granvik, M., et al. 2023, , 4, 29, doi: 10.3847/PSJ/acb25b
  * Fernandez & Ip (1984) Fernandez, J. A., & Ip, W. H. 1984, Icarus, 58, 109, doi: 10.1016/0019-1035(84)90101-5
  * Filacchione et al. (2022) Filacchione, G., Ciarniello, M., Fornasier, S., & Raponi, A. 2022, arXiv e-prints, arXiv:2210.02741, doi: 10.48550/arXiv.2210.02741
  * Fischer & Valenti (2005) Fischer, D. A., & Valenti, J. 2005, ApJ, 622, 1102, doi: 10.1086/428383
  * Fitzsimmons et al. (2023) Fitzsimmons, A., Meech, K., Matrà, L., & Pfalzner, S. 2023, arXiv e-prints, arXiv:2303.17980, doi: 10.48550/arXiv.2303.17980
  * Forbes & Loeb (2019) Forbes, J. C., & Loeb, A. 2019, ApJ, 875, L23, doi: 10.3847/2041-8213/ab158f
  * Frankel et al. (2020) Frankel, N., Sanders, J., Ting, Y.-S., & Rix, H.-W. 2020, ApJ, 896, 15, doi: 10.3847/1538-4357/ab910c
  * Fulton et al. (2021) Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, ApJS, 255, 14, doi: 10.3847/1538-4365/abfcc1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * García Pérez et al. (2016) García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2016, AJ, 151, 144, doi: 10.3847/0004-6256/151/6/144
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15, doi: 10.1051/0004-6361/201833718
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93, doi: 10.3847/1538-4357/ab5362
  * Guilbert-Lepoutre et al. (2015) Guilbert-Lepoutre, A., Besse, S., Mousis, O., et al. 2015, Space Sci. Rev., 197, 271, doi: 10.1007/s11214-015-0148-9
  * Hallatt & Wiegert (2020) Hallatt, T., & Wiegert, P. 2020, AJ, 159, 147, doi: 10.3847/1538-3881/ab7336
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Holtzman et al. (2018) Holtzman, J. A., Hasselquist, S., Shetrone, M., et al. 2018, AJ, 156, 125, doi: 10.3847/1538-3881/aad4f9
  * Horta et al. (2023) Horta, D., Schiavon, R. P., Mackereth, J. T., et al. 2023, MNRAS, 520, 5671, doi: 10.1093/mnras/stac3179
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * 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), 126, e06706, doi: 10.1029/2020JE006706
  * Jewitt (2003) Jewitt, D. 2003, Earth Moon and Planets, 92, 465, doi: 10.1023/B:MOON.0000031961.88202.60
  * Jewitt et al. (2017) Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJ, 850, L36, doi: 10.3847/2041-8213/aa9b2f
  * Jewitt & Seligman (2022) Jewitt, D., & Seligman, D. Z. 2022, arXiv e-prints, arXiv:2209.08182, doi: 10.48550/arXiv.2209.08182
  * Johnson & Li (2012) Johnson, J. L., & Li, H. 2012, ApJ, 751, 81, doi: 10.1088/0004-637X/751/2/81
  * Jönsson et al. (2020) Jönsson, H., Holtzman, J. A., Allende Prieto, C., et al. 2020, AJ, 160, 120, doi: 10.3847/1538-3881/aba592
  * Jurić et al. (2008) Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864, doi: 10.1086/523619
  * Kaib & Volk (2022) Kaib, N. A., & Volk, K. 2022, arXiv e-prints, arXiv:2206.00010, doi: 10.48550/arXiv.2206.00010
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Leung & Bovy (2019) Leung, H. W., & Bovy, J. 2019, MNRAS, 489, 2079, doi: 10.1093/mnras/stz2245
  * Levine et al. (2021) Levine, W. G., Cabot, S. H. C., Seligman, D., & Laughlin, G. 2021, ApJ, 922, 39, doi: 10.3847/1538-4357/ac1fe6
  * Lian et al. (2022) Lian, J., Zasowski, G., Mackereth, T., et al. 2022, MNRAS, 513, 4130, doi: 10.1093/mnras/stac1151
  * Lintott et al. (2022) Lintott, C., Bannister, M. T., & Mackereth, J. T. 2022, ApJ, 924, L1, doi: 10.3847/2041-8213/ac41d5
  * Lu et al. (2020) Lu, C. X., Schlaufman, K. C., & Cheng, S. 2020, AJ, 160, 253, doi: 10.3847/1538-3881/abb773
  * Mackereth et al. (2017) Mackereth, J. T., Bovy, J., Schiavon, R. P., et al. 2017, MNRAS, 471, 3057, doi: 10.1093/mnras/stx1774
  * Mackereth et al. (2019) Mackereth, J. T., Bovy, J., Leung, H. W., et al. 2019, MNRAS, 489, 176, doi: 10.1093/mnras/stz1521
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94, doi: 10.3847/1538-3881/aa784d
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77, doi: 10.3847/1538-4357/835/1/77
  * Marshall et al. (2006) Marshall, D. J., Robin, A. C., Reylé, C., Schultheis, M., & Picaud, S. 2006, A&A, 453, 635, doi: 10.1051/0004-6361:20053842
  * Mashchenko (2019) Mashchenko, S. 2019, MNRAS, 489, 3003, doi: 10.1093/mnras/stz2380
  * Mashonkina et al. (2019) Mashonkina, L. I., Neretina, M. D., Sitnova, T. M., & Pakhomov, Y. V. 2019, Astronomy Reports, 63, 726, doi: 10.1134/S1063772919090063
  * McGlynn & Chapman (1989) McGlynn, T. A., & Chapman, R. D. 1989, ApJ, 346, L105, doi: 10.1086/185590
  * Meech et al. (2016) Meech, K. J., Yang, B., Kleyna, J., et al. 2016, Science Advances, 2, e1600038, doi: 10.1126/sciadv.1600038
  * 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
  * Moro-Martín (2018) Moro-Martín, A. 2018, ApJ, 866, 131, doi: 10.3847/1538-4357/aadf34
  * Moro-Martín (2019) —. 2019, ApJ, 872, L32, doi: 10.3847/2041-8213/ab05df
  * Moro-Martín (2022) —. 2022, arXiv e-prints, arXiv:2205.04277, doi: 10.48550/arXiv.2205.04277
  * Nimmo et al. (2018) Nimmo, F., Kretke, K., Ida, S., Matsumura, S., & Kleine, T. 2018, Space Sci. Rev., 214, 101, doi: 10.1007/s11214-018-0533-2
  * Öberg & Bergin (2021) Öberg, K. I., & Bergin, E. A. 2021, Phys. Rep., 893, 1, doi: 10.1016/j.physrep.2020.09.004
  * Osborn & Bayliss (2020) Osborn, A., & Bayliss, D. 2020, MNRAS, 491, 4481, doi: 10.1093/mnras/stz3207
  * 'Oumuamua ISSI Team et al. (2019) 'Oumuamua ISSI Team, Bannister, M. T., Bhandare, A., et al. 2019, Nature Astronomy, 3, 594, doi: 10.1038/s41550-019-0816-x
  * Pfalzner et al. (2021) Pfalzner, S., Aizpuru Vargas, L. L., Bhandare, A., & Veras, D. 2021, A&A, 651, A38, doi: 10.1051/0004-6361/202140587
  * Pfalzner & Bannister (2019) Pfalzner, S., & Bannister, M. T. 2019, ApJ, 874, L34, doi: 10.3847/2041-8213/ab0fa0
  * Pfalzner et al. (2020) Pfalzner, S., Davies, M. B., Kokaia, G., & Bannister, M. T. 2020, ApJ, 903, 114, doi: 10.3847/1538-4357/abb9ae
  * Raymond et al. (2018a) Raymond, S. N., Armitage, P. J., & Veras, D. 2018a, ApJ, 856, L7, doi: 10.3847/2041-8213/aab4f6
  * Raymond et al. (2018b) Raymond, S. N., Armitage, P. J., Veras, D., Quintana, E. V., & Barclay, T. 2018b, MNRAS, 476, 3031, doi: 10.1093/mnras/sty468
  * Recio-Blanco et al. (2014) Recio-Blanco, A., de Laverny, P., Kordopatis, G., et al. 2014, A&A, 567, A5, doi: 10.1051/0004-6361/201322944
  * Sabotta et al. (2021) Sabotta, S., Schlecker, M., Chaturvedi, P., et al. 2021, A&A, 653, A114, doi: 10.1051/0004-6361/202140968
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Schwamb et al. (2023) Schwamb, M. E., Jones, R. L., Yoachim, P., et al. 2023, ApJS, 266, 22, doi: 10.3847/1538-4365/acc173
  * Seligman & Laughlin (2020) Seligman, D., & Laughlin, G. 2020, ApJ, 896, L8, doi: 10.3847/2041-8213/ab963f
  * Seligman et al. (2022) Seligman, D. Z., Rogers, L. A., Cabot, S. H. C., et al. 2022, , 3, 150, doi: 10.3847/PSJ/ac75b5
  * Seligman et al. (2023) Seligman, D. Z., Farnocchia, D., Micheli, M., et al. 2023, , 4, 35, doi: 10.3847/PSJ/acb697
  * Sellwood & Binney (2002) Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785, doi: 10.1046/j.1365-8711.2002.05806.x
  * Stern (1990) Stern, S. A. 1990, PASP, 102, 793, doi: 10.1086/132704
  * Stern & Shull (1988) Stern, S. A., & Shull, J. M. 1988, Nature, 332, 407, doi: 10.1038/332407a0
  * Tinsley & Cameron (1974) Tinsley, B. M., & Cameron, A. G. W. 1974, Ap&SS, 31, 31, doi: 10.1007/BF00642598
  * Trilling et al. (2018) Trilling, D. E., Mommert, M., Hora, J. L., et al. 2018, AJ, 156, 261, doi: 10.3847/1538-3881/aae88f
  * Tsapras (2018) Tsapras, Y. 2018, Geosciences, 8, 365, doi: 10.3390/geosciences8100365
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1, doi: 10.48550/arXiv.1504.05966
  * Vandorou et al. (2023) Vandorou, A., Dang, L., Bennett, D. P., et al. 2023, arXiv e-prints, arXiv:2302.01168, doi: 10.48550/arXiv.2302.01168
  * Vickers et al. (2021) Vickers, J. J., Shen, J., & Li, Z.-Y. 2021, ApJ, 922, 189, doi: 10.3847/1538-4357/ac27a9
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wiersma et al. (2009) Wiersma, R. P. C., Schaye, J., Theuns, T., Dalla Vecchia, C., & Tornatore, L. 2009, MNRAS, 399, 574, doi: 10.1111/j.1365-2966.2009.15331.x
  * 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
  * Zasowski et al. (2013) Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, AJ, 146, 81, doi: 10.1088/0004-6256/146/4/81
  * Zasowski et al. (2017) Zasowski, G., Cohen, R. E., Chojnowski, S. D., et al. 2017, AJ, 154, 198, doi: 10.3847/1538-3881/aa8df9