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

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

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

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

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

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

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

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

We assume that the interstellar objects within the solar neighborhood follow a Schwarzschild velocity distribution with respect to the LSR, with velocity dispersions σ_R , σ_ϕ , and σ_z . The right-handed galactic coordinate system is oriented by the vertex deviation, v_d , for the primary eigenvector of the ellipsoidal distribution function. The spatial number density of interstellar objects, n(Δ v), in a specific volume of phase space, Δ v = (Δv₁, Δv₂, Δv₃), is given by

$$\begin{eqnarray}&&n({\boldsymbol{\Delta }}v)=\displaystyle \frac{n}{8}\left(\,\prod _{i=1}^{3}\mathrm{erf}\left(\displaystyle \frac{{v}_{i}}{\sqrt{2}{\sigma }_{i}}\right){| }_{{v}_{i}-{\rm{\Delta }}{v}_{i}}^{{v}_{i}+{\rm{\Delta }}{v}_{i}}\right).\end{eqnarray} \tag{ 1 }$$

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

We adopt the values of σ_R , σ_ϕ , σ_z , and v_d for G stars (Binney & Merrifield 1998; Seligman & Laughlin 2018). It is important to note that different stellar populations correspond to different typical ages and velocity dispersions, as summarized in Table 1. The typical ages of interstellar objects are somewhat unconstrained. We verified that the results of our simulations are not sensitively dependent on the spectral type that we use. It would be informative to perform the calculations presented in Sections 4 and 5 for spectral classes representative of older and younger populations, but this is outside of the scope of this paper.

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

$$\begin{eqnarray}&&\xi \approx 1+{\left(\displaystyle \frac{42\,\mathrm{km}\,{{\rm{s}}}^{-1}}{{v}_{\infty }+\beta }\right)}^{2}.\end{eqnarray} \tag{ 2 }$$

Here v_∞ is the hyperbolic excess velocity with respect to the Sun, and β is an artificial smoothing factor that we set to β = 2 km s⁻¹ to avoid the hyperbolic case where v_∞ = 0 (Seligman & Laughlin 2018). We compute the probability that an interstellar object enters the solar system with a given velocity using Equations (1) and (2).

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

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

Zoom In Zoom Out Reset image size

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

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

$$\begin{eqnarray}&&{\boldsymbol{R}}=-k{{\boldsymbol{v}}}_{\infty },\end{eqnarray} \tag{ 3 }$$

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

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{\mathrm{sky}}=-\left(\sqrt{{R}^{2}-{b}^{2}}\right){\hat{v}}_{\infty },\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&{\boldsymbol{b}}=\left[\begin{array}{c}b\cos \phi \cos \psi \\ b\cos \phi \sin \psi \\ b\sin \phi \end{array}\right],\end{eqnarray} \tag{ 5 }$$

where the angle ϕ is given by

$$\begin{eqnarray}&&\tan \phi =-\left(\displaystyle \frac{\cos \psi \,{v}_{\infty ,x}+\sin \psi \,{v}_{\infty ,y}}{{v}_{\infty ,z}}\right).\end{eqnarray} \tag{ 6 }$$

Here the subscripts on v_∞ indicate the component of the vector. This allows us to construct b for each initial velocity vector. The initial position vector, R , is given by

$$\begin{eqnarray}&&{\boldsymbol{R}}={{\boldsymbol{R}}}_{\mathrm{sky}}+{\boldsymbol{b}}.\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

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

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

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

$$\begin{eqnarray}&&m=H+2.5{\mathrm{log}}_{10}\left(\displaystyle \frac{{d}_{\mathrm{BS}}^{2}{d}_{\mathrm{BO}}^{2}}{p(\theta ){d}_{\mathrm{OS}}^{4}}\right).\end{eqnarray} \tag{ 8 }$$

The parameters d_BO, d_OS, and d_BS represent the distances between the body and observer, Sun and observer, and body and Sun, respectively. The phase angle, θ, and the phase integral, p(θ), are defined by

$$\begin{eqnarray}&&\cos \theta =\frac{{d}_{BO}^{2}+{d}_{\mathrm{BS}}^{2}-{d}_{\mathrm{OS}}^{2}}{2{d}_{\mathrm{BO}}{d}_{\mathrm{BS}}}\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&p(\theta )=\displaystyle \frac{2}{3}\left(\left(1-\displaystyle \frac{\theta }{\pi }\right)\cos \theta +\displaystyle \frac{1}{\pi }\sin \theta \right).\end{eqnarray} \tag{ 10 }$$

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

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

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

• 1.  
  Its apparent magnitude is less than 24 and higher than 16. We consider m = 24 to be the apparent magnitude of the dimmest objects that the LSST is capable of observing. Because of the LSST's 30 s exposures, sources brighter than 16th magnitude would saturate, meaning they would not be detected by the LSST.
• 2.  
  Its altitude is higher than 30° relative to the latitude/longitude of Rubin Observatory, and its decl. is less than +12°. This ensures that the object is within the nominal field of view of the LSST, as the Wide-Fast-Deep is likely to extend to +12° decl. for most right ascensions, according to Figure 1 of Bianco et al. (2022).
• 3.  
  The Sun's altitude is lower than −18° relative to the latitude/longitude of the LSST. This ensures that the object is observed before/after astronomical twilight, allowing for detailed observations of faint point sources.

These criteria ensure that when each object reaches m ≤ 24, the angular distance between the object and the Sun is large. We find that out of the ∼10⁶ simulated interstellar objects, 6.956% of them reach a minimum apparent magnitude brighter than 24, and 2.828% of them are detectable by the LSST. We define a "reachable" object as one that can be intercepted with a rocket sent from Earth with an impulsive change of velocity, Δv < 20 km s⁻¹. We further explore the reachable population in Section 5.

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

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

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

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Table 2, we show the fractional percentages of objects that will be detected by the LSST and attainable targets for a range of limiting Δv. To convert these fractions to detection rates, we consider a typical interstellar object traveling within the 5 au sphere centered at the Sun. We assume a spatial number density of n_o ∼ 1 × 10⁻¹ au⁻³, lower than that calculated by Do et al. (2018) given the nondetection of an additional 'Oumuamua-like object since 2017 (Levine et al. 2021). With this assumed spatial density, there are ∼50 interstellar objects within the 5 au sphere at any given time.

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

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

$$\begin{eqnarray}&&| | d| | =\displaystyle \frac{1}{{\left(2\pi \right)}^{2}}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{2\pi }{\left(2-2\cos \left({\theta }_{1}-{\theta }_{2}\right)\right)}^{1/2}r\,d{\theta }_{1}d{\theta }_{2},\end{eqnarray} \tag{ 11 }$$

which yields ∣∣d∣∣ = 4/π r. For a sphere,

$$\begin{eqnarray}\begin{array}{rcl}| | d| | & = & \displaystyle \frac{1}{{\left(4\pi \right)}^{2}}{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{\pi }{\displaystyle \int }_{0}^{2\pi }{\displaystyle \int }_{0}^{\pi }\left(2-2\cos {\phi }_{1}\cos {\phi }_{2}\right.\\ & & -2\cos {\theta }_{1}\cos {\theta }_{2}\sin {\phi }_{1}\sin {\phi }_{2}\\ & & {\left.-2\sin {\theta }_{1}\sin {\theta }_{2}\sin {\phi }_{1}\sin {\phi }_{2}\right)}^{1/2}\\ & & \times \,r\,\sin {\phi }_{1}\sin {\phi }_{2}d{\phi }_{1}d{\theta }_{1}d{\phi }_{2}d{\theta }_{2},\end{array}\end{eqnarray} \tag{ 12 }$$

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

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

These values correspond to crossing times of ∼0.8 and ∼1.6 yr. To convert the fractions to detection rates, we use the following equation:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{rate} & = & 33.11\,{\mathrm{yr}}^{-1}\left(\displaystyle \frac{{n}_{o}}{1\times {10}^{-1}\,{\mathrm{au}}^{-3}}\right)\\ & & \times \left(\displaystyle \frac{V}{4/3\pi {(5\,\mathrm{au})}^{3}}\right)\left(\displaystyle \frac{v}{20\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{(4/3)5\,\mathrm{au}}{| | d| | }\right),\end{array}\end{eqnarray} \tag{ 13 }$$

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

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

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

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

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

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

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

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

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

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

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