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

The Milky Way is thought to host a huge population of interstellar objects (ISOs), numbering approximately 1015 pc−3 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, first to predict the properties of this Galactic ISO population by combining models of processes across planetary and Galactic scales, and second to make inferences about the processes being modeled, by comparing the predicted population to what is observed. We predict the spatial and compositional distribution of the Galaxy’s population of ISOs by modeling the Galactic stellar population with data from the APOGEE survey and combining this with a protoplanetary disk chemistry model. Selecting the 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 neighborhood 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 progresses and more ISOs are observed. Finally, we explore generalizations of this framework to other Galactic populations, such as that of exoplanets.

1. INTRODUCTION 1I/'Oumuamua (Meech et al. 2017) and 2I/Borisov 1 are the first two observed samples from a highly numerous population: interstellar objects (ISOs).Estimated to number ∼ 10 15 pc −3 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 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.

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.

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 ∼ 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(Zasowski et al. , 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 [Fe/H] and [α/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 [Fe/H] and alpha elements [α/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 [Fe/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 ≤ 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 kpc and 12 kpc and height above the disk z of −5 kpc to 5 kpc.This gives us a sample of 80 958 red giants.

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 , . . ., O N is a realisation of an inhomogeneous Poisson point process.This process is a random distribution of points defined by a rate function λ(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, V λ(O) dO.It follows that the probability of finding a point (i.e. an observed star) with observables in the infinitesimal volume δO is given by λ(O) δO, and the total number of points (i.e.stars observed) is a Poisson random variable with mean and variance Λ = λ(O) dO.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 λ(O | θ), where θ parameterises the model, then the likelihood of the model given (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 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 τ .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 τ .We used the effective selection function implementation in the apogee2 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 ⊙ (Kroupa 2001), cut to 1 ≤ 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 mwdust3 , 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 [Fe/H] and [α/Fe], then fit a separate number density model n 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, parameterised by an amplitude logA and two scale parameters a R and a z .R 0 = 8.1 kpc 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 τ 0 and variance 1/ω, 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 λ(field, D, τ ) for each monoabundance population, the number density of red giants needs to be multiplied a Jacobian factor |J(field, D)| = Ω 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 S(field, D, τ ).Thus the rate function for each monoabundance population is given by This particular form for the density profile has the advantage that the Poisson point process likelihood takes the tractable form where the sum over every star in the monoabundance population i ln λ(O i | θ) 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, ⟨R⟩ and ⟨|z|⟩ as the mean coordinate values, and ⟨τ ⟩ and ⟨τ 2 ⟩ as the mean age and age squared.To build our model of the Milky Way disk between R = 4 kpc − 12 kpc and |z| = 0 − 5 kpc, we found a best-fit model for each monoabundance population by maximising this likelihood with respect to logA, a R , a z , τ 0 , and ω, visually checking each fit to ensure the global maximum had been found.The results are plotted in Figure 1.
Reliable ages are not available for stars [Fe/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 [Fe/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 [Fe/H], low [α/Fe] monoabundance populations have the low vertical scale lengths and high radial scale lengths which form the thin disk, whereas the old low [Fe/H], high [α/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.

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 morte4 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 ⊙ as in § 2. Using the isochrones with metallicities corresponding to each monoabundance population, weighted by the fitted age distribution g(τ | τ 0 , ω) for that monoabundance population, we calculate N giants /N sm : the fraction of all stars ever created that are currently in the red giant phase, which is again defined by 1 ≤ 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, ⟨M int ⟩.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 ⟨M int ⟩ from the same Kroupa IMF.All combined, the sine morte stellar mass density is given by 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.
For the monoabundance populations with [Fe/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 τ 0 = 12 Gyr and standard deviation ω − 1 2 = 1 Gyr, which minimises N giants /N sm .Even with this upper limit, below we find that stars with [Fe/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 [Fe/H]: so when we evaluate the model we sum the sine morte distribution over the bins in [α/Fe] to get a distribution in only [Fe/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

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.

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  (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 H2O , as this varies significantly and decreases monotonically with [Fe/H] in the models of Bitsch & Battistini (2020).To obtain a smooth map from [Fe/H] to f H2O , 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 ≤ [Fe/H] ≤ 0.4 we assume that the relationship between f H2O and [Fe/H] continues to be monotonic, with f H2O remaining high beyond the low [Fe/H] limit and remaining low beyond the high [Fe/H] limit: allowing us to track the fraction of ISOs with high, low, and intermediate water mass fraction.This range corresponds to 0.07 ≤ f H2O ≤ 0.51 in water mass fraction.They also confirmed the trend between [Fe/H] and f H2O found in Bitsch & Battistini (2020).However, they also found that for the APOGEE survey there was a smaller variation in f H2O over the same range in [Fe/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 H2O .
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 * and its metal mass fraction Z -well approximated by Z ⊙ 10 [Fe/H] for small values of Z (Z ⊙ = 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 * • 10 β[Fe/H] .Here β = 1 corresponds to the simple assumption that the number of ISOs produced is proportional to the mass of planet-forming material.β = 0 corresponds to no metallicity dependence at all.However, we expect β > 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 2023).We predict only the normalised distribution of ISO water mass fractions.

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.A fraction of a planetary system's bound planetesimals will become unbound later in a system's life, whether from its Oort cloud's continuous interaction with the Galactic environment, or in post-main-sequence escape (Veras et al. 2011(Veras et al. , 2014;;Levine et al. 2023); for simplicity we omit this smaller population here.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 effect5 , 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 H2O and position in the Galaxy x in the present day, denoted n ISO (x, f H2O ), would be equal to the number of ISOs produced by stars currently at the same position and of the corresponding metallicity [Fe/H].Under the model of § 3.1, the number of ISOs produced by a star of mass M * and metallicity [Fe/H] is proportional to M * • 10 β[Fe/H] , thus where ρ(x, [Fe/H]) is the mass density distribution of stars at position x with the metallicity [Fe/H] corresponding to the ISO water mass fraction f H2O , and d[Fe/H]/ df H2O is the gradient of the relationship between [Fe/H] and f H2O described in § 3.1.
The fact that stars do in fact die is then corrected for by replacing ρ(x, [Fe/H]) in this equation with ρ sm (x, [Fe/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 < 10 km s −1 , 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 ≳ 20 km s −1 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: 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 H2O 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 ISO distribution we calculate with This gives us the distribution of ISOs within the bounds of the protoplanetary disk chemical model: 0.07 ≤ f H2O ≤ 0.51.Outside of this range, we can calculate the fraction of ISOs with f H2O ≤ 0.07 and f H2O ≥ 0.51, by assuming that the relation between [Fe/H] and f H2O remains monotonic.Thus, all stars with [Fe/H] ≤ −0.4 contribute ISOs with f H2O ≥ 0.51, and all stars with [Fe/H] ≥ 0.4 contribute ISOs with f H2O ≤ 0.07.

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, β, with β = 1 for our principal prediction and β = 0 as an alternate prediction.

Principal Prediction: β = 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 β = 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.
Table 1 lists the the fraction of ISOs within and outside either end of the Bitsch & Battistini (2020) protoplanetary disk chemical model f H2O range, 0.07 ≤ f H2O ≤ 0.51.We assess both the distribution of ISOs at the position of the Sun, at R = 8.1 kpc, z = 0.021 kpc (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 [Fe/H] limit of the ρ sm ([Fe/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 H2O range.This is because of the exponential dependence of the number of ISOs produced by each star on [Fe/H] in Eq. 8.  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.

Local Milky Way Disk
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 [Fe/H] distribution and the ISO f H2O 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 [Fe/H] and f H2O 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 kpc), where both the stellar [Fe/H] distribution and therefore the ISO f H2O distributions are approximately midway between the high and low extremes.In these plots, we have normalised the distributions in [Fe/H] and f H2O such that at each value of R the integral over [Fe/H] or f H2O 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.

Alternate Prediction: β = 0
Setting β = 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 [Fe/H] to ISO water mass fractions f H2O .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 β = 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.
Assuming β = 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 [Fe/H].This means that the low metallicity tails of the ρ sm ([Fe/H]) distributions in Figure 3 now contribute a significant number of ISOs with 0.51 < f H2O : almost 20% when averaged over the Milky Way disk.A comparison of these results to those of Lintott et al. (2022)  Table 2. Alternate prediction for fraction of ISOs in each fH 2 O range, evaluated at the position of the Sun and integrated over the whole Milky Way with β = 0: assuming the number of ISOs produced by each star is independent of the star's metallicity.

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 H2O | β), parameterised by the ISO production power law slope β.This however could be generalised to a distribution in any set of observable properties O, and parameterised by any set of parameters θ.These distributions are also probability density functions for the value of O of a randomly selected ISO, so p(O | θ) δO is the probability of an ISO having properties in the infinitesimal volume δO.Extending this to multiple ISOs, a sample of ISOs with properties O 1 , . . ., O N will have likelihood 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: This posterior can be used to calculate estimates and confidence intervals for the values of the parameters θ.An example of this calculation is given in section 6.3.

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 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 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 [Fe/H], low velocity dispersion), thick disk (low [Fe/H], high velocity dispersion), and the halo (very low [Fe/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).A prediction and discussion of the chemodynamical distribution of ISOs using stellar velocities from the Gaia satellite is currently in preparation.
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 H 2 O, CO 2 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.

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 would mean 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.
Other variables affecting planet formation however could flatten or even reverse this gradient.Based on 28 planetary microlensing events, Koshimoto et al. (2021) find tentative evidence of a decrease in the planet occurrence rate towards the Galactic centre, possibly due to the increased rate of close stellar encounters in the Galactic bulge.

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 H2O over the same range in [Fe/H] compared to the GALAH data used in Bitsch & Battistini (2020).If the true variation in f H2O 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-like galaxy from the EAGLE hydrodynamical cosmological simulation (Schaye et al. 2015), using a model equivalent to that of this work with β = 0. Whereas we predict a single-peaked f H2O 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 [Fe/H] distribution than the Milky Way, with many more stars outside of the [Fe/H] range of the protoplanetary disk model.As a smoothed particle hydrodynamics simulation, EAGLE is susceptible to producing galaxies with [Fe/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.

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.

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, β.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 β and Bayes' theorem to calculate a posterior distribution for β.
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 H2O = 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 CO 2 relative to H 2 O.
With this observed distribution of f H2O we can then use the Bayesian framework of § 5.1 to calculate a posterior distribution for β.We use our model for the distribution of ISOs around the Sun, and taking a uniform prior on β means that the posterior is simply proportional to the likelihood, equal to the value of the f H2O distribution.For f H2O = 0.3 this posterior is maximised by a value of β = 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 β is β ∈ (−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 β -but with one known value of f H2O 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 β and other parameters used in models with this framework.Table 3. Predicted fraction of ISOs in each fH 2 O range, evaluated at the position of the Sun and integrated over the Milky Way dist, with maximum a posteriori estimate of β = 1.33.

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 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.

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 [Fe/H] and [α/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 H2O 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.

Figure 1 .
Figure 1.The best-fit values of the density modelling parameters for each monoabundance population with 20 or more observed stars, and ρsm/ngiants, the ratio between the sine morte stellar mass density and the red giant number density, explained in section 2.3.

Figure 2 .
Figure 2. Comparison of the [Fe/H] distribution of red giants and the [Fe/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 [Fe/H], as shown in Fig. 1.

Figure 3 .
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 [Fe/H] for which we model variations in the water mass fraction fH 2 O of ISOs 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 [α/Fe] as well as [Fe/H].They found that the water mass fraction of planetesimals, consistent with the Bitsch & Battistini (2020) assumptions we use, is most dependent on [Fe/H].

Figure 4 .
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 β = 1.

Figure 5 .
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 [Fe/H] or fH 2 O is unity.Dashed lines show the median value of each distribution at each R.

Figure 6 .
Figure6.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 β = 0: assuming the number of ISOs produced by each star is independent of the star's metallicity.

ISO fH 2
O range Fraction of ISOs around Sun Fraction of ISOs in Milky Way Disk fH 2 O < 0

Figure 7 .
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 β = 1.33.At fH 2 O = 0.3 a vertical line marks the measured water mass fraction of 2I/Borisov.

Table 1 .
Primary prediction for the fraction of ISOs in each fH 2 O range evaluated at the position of the Sun and integrated over whole Milky Way, with β = 1.