Identifying Interstellar Objects Trapped in the Solar System through Their Orbital Parameters

The first interstellar object, 'Oumuamua, was discovered in the solar system by the Pan-STARRS telescope (Meech et al. 2017; Micheli et al. 2018). Several follow-up studies of 'Oumuamua were conducted to better understand its origin and composition (Bannister et al. 2017; Bolin et al. 2017; Gaidos et al. 2017; Jewitt et al. 2017; Mamajek 2017; Ye et al. 2017; Fitzsimmons et al. 2018). The detection of 'Oumuamua allowed for a calibration of the number density of objects of similar size, estimated to be ∼10¹⁵ pc⁻³ (Do et al. 2018). This updated number density is much higher than the previous estimate of ≲10¹² pc⁻³ (Moro-Martín et al. 2009; Engelhardt et al. 2017). Lingam & Loeb (2018) used this calibration to estimate the capture rate of 'Oumuamua-size interstellar objects by means of gravitational interactions with Jupiter and the Sun to be 1.2 × 10⁻² yr⁻¹, and the resulting number of 'Oumuamua-size interstellar objects bound to the solar system at any given time to be ∼6 × 10³. Observing or visiting such objects could facilitate the search for signs of extraterrestrial life locally, without the need to send interstellar probes (Loeb 2018). However, such a search would come with a caveat that most asteroids in the solar system reside outside the habitable zone. Here, we explore whether or not it is possible to identify trapped interstellar objects through their orbital parameters alone.

The outline of this Letter is as follows. In Section 2 we explore the orbital parameters of trapped interstellar objects by simulating their gravitational capture by the Jupiter–Sun system. In Section 3 we describe our results, and in Section 4 we explore their implications in identifying interstellar objects in our solar system through their orbital parameters. In Section 5 we identify known objects that may be of interstellar origin, and in Section 6 we investigate the future detectability of the trapped interstellar object population. Finally, Section 7 summarizes our main conclusions.

In this simulation, we consider interactions only with Jupiter as it is the most massive planet. To model the motion of interstellar objects under the gravitational influence of the Jupiter–Sun system, we developed a Python code that randomly initializes and integrates the motions of particles from their points of closest approach to Jupiter to both the past and the future, and searches for particles that are initially unbound yet end in bound orbits after their gravitational interaction with Jupiter. The Python code created for this work used the open-source N-body integrator software REBOUND¹ to trace the motions of particles under the gravitational influence of the Jupiter–Sun system (Rein & Liu 2012).

We initialize the simulation with the Sun, Jupiter, and a volume of test particles surrounding Jupiter at their distances of closest approach to the planet. The Sun and Jupiter define the ecliptic plane. We set the range of impact parameter, b, of the particles relative to Jupiter to be between 1 and $10\,{R}_{{\rm{J}}}$, where R_J is Jupiter's radius. The upper limit of $10\,{R}_{{\rm{J}}}$ was chosen after runs of the complete code with b set to vary between 0 and $\sim {10}^{3}\,{R}_{{\rm{J}}}$ resulted in no captured particles for $b\gt 10\,{R}_{{\rm{J}}}$. To choose the value of b for each particle, we draw randomly from a weighted distribution in which the likelihood of a particle having any given value of b between 1 and $10\,{R}_{{\rm{J}}}$ is proportional to b², due to the fact that each value of b represents a spherical infinitesimal shell of volume 4πb² db around Jupiter.

For each value of b that is chosen, we randomly pick an angle within the ecliptic plane between 0 and 2π, as well as a zenith angle between 0 and π. Using these two angles, we set the direction of each particle's velocity vector.

We randomly draw the initial speed of each particle, v_∞, from a Maxwellian distribution with velocity dispersion 40 km s⁻¹, because this is the approximate characteristic distribution of the velocity dispersion for stars in the Solar neighborhood (Binney & Tremaine 2008; Li & Adams 2016). From the constructed velocity vector, we subtract the motion of the Sun relative to the local standard of rest (LSR), ${(U,V,W)}_{\mathrm{LSR}}=(11.1,12.2,7.3)$ $\mathrm{km}\ {{\rm{s}}}^{-1}$ (Schonrich et al. 2010). Finally, we pick a random position in Jupiter's orbit and account for its velocity vector. We subsequently calculate the speed of each particle using conservation of energy at closest approach, using the equation

$$\begin{eqnarray}&&{{\boldsymbol{E}}}_{k}={{\boldsymbol{E}}}_{k,\infty }-{{\boldsymbol{E}}}_{g,\odot }-{{\boldsymbol{E}}}_{g,\mathrm{Jupiter}}.\end{eqnarray} \tag{ 1 }$$

While the gravitational influence of the Sun is not very large at the distance of Jupiter, we have included it in our calculation of each interstellar object's velocity vector. The resulting velocities approximate the characteristic velocity of interstellar objects; however, any other object with the same velocity, such as an Oort cloud object ejected at a particular speed, could be misconstrued with interstellar objects in this analysis.

To ensure that each particle is at its distance of closest approach, we require the position vector to lie in the plane perpendicular to the velocity vector relative to Jupiter. For each particle, we pick a random angle between 0 and 2π to determine in which direction the position vector points within this plane. Using the impact parameter, b, and the angle within the plane perpendicular to the velocity vector, we construct each particle's position vector. At this point, we have fully initialized each particle with both a position and a velocity vector.

In the first stage of the simulation, we integrate all of the particles backward in time and determine which ones are initially unbound. We use the WHFast² integrator in REBOUND to trace each particle from t = 0 to an earlier time −t_i (Rein & Tamayo 2015), where t_i is an amount of time to sample either side of the closest approach to Jupiter. The only constraint on t_i is that it is a time interval at and above which the results of the simulation do not change, on the order of Jupiter's orbital period. We then compute the escape velocity relative to the Jupiter–Sun system for each particle using the expression below and label all particles with speed $v(-{t}_{i})\gt {v}_{\mathrm{esc}}$ as initially unbound, where

$$\begin{eqnarray}&&{v}_{\mathrm{esc}}=\sqrt{\displaystyle \frac{2{{GM}}_{\odot }}{{d}_{\odot }}+\displaystyle \frac{2{{GM}}_{J}}{{d}_{J}}},\end{eqnarray} \tag{ 2 }$$

with d_⊙ and d_J being the distances from the Sun and Jupiter, respectively. In the second stage of the simulation, we integrate the particles with unbound initial conditions forward in time and find which ones will end in bound orbits around the Sun. We use WHFast to integrate each particle from t = 0 to t_i. We again compute each particle's escape velocity and determine that all particles with speed $v({t}_{i})\lt {v}_{\mathrm{esc}}$ become gravitationally bound and therefore satisfy the condition of capture.

Finally, our Python code calculates initial conditions and final orbital parameters for the captured particles. At $-{t}_{i}$, we calculate each particle's incoming zenith angle, θ, where 90° corresponds to the plane of the ecliptic. We then use Orbital,³ an open-source orbital mechanics package, to compute each particle's semimajor axis a, eccentricity e, and inclination i, at t_i. We subsequently compute orbital period, T, perihelion distance, q, and aphelion distance, Q.

We ran our Python code for 10⁵ particles, which resulted in 3202 instances of capture. We verified that the percentage of interstellar objects captured by Jupiter remained at ∼3% for 2× and 5× particles. The distribution of semimajor axis, a, is positively skewed with a median of 13 au, 25th and 75th percentiles of 4.3 au and 16 au, respectively, and a range of 0.93 au to ∼5 × 10³ au. The distribution of eccentricity, e, is negatively skewed with a median of 0.74, 25th and 75th percentiles of 0.53 and 0.89, respectively, and a range of up to ∼1. The distribution of inclination, i, is positively skewed with a median of 82°, 25th and 75th percentiles of 41° and 124°, respectively, and a range of up to 180°. The probability distribution of inclination, i, is displayed in Figure 1, and is normalized to unit area. Owing to the nearly flat distribution of inclinations in Figure 1 that results from the broad initial velocity distribution, the population of interstellar objects dominates over the Centaur inclination distribution at an inclination angle of 48°.

Zoom In Zoom Out Reset image size

Is it possible to identify trapped interstellar objects through their orbital parameters alone? To explore this question, we compare our orbital element distributions of trapped ∼100 m size interstellar objects to those of ∼100 m size background objects that originated in the solar system. The population of solar system objects with orbital elements similar to our theoretical population are the Centaurs, which are defined by Grav et al. (2011) as objects that have $5.5\leqslant a\leqslant 30\,\mathrm{au}$ and $q\gt 5.5\,\mathrm{au}$. We will call our theoretical population "Centaur-like" as the objects have $\sim 3\leqslant a\leqslant 30\,\mathrm{au}$ and $q\gt 2\,\mathrm{au}$.

The size distribution of Centaurs is not well constrained, and ∼100 m diameter Centaurs have yet to be discovered. We will estimate the number density of ∼100-m-sized Centaurs by using a model derived by Bauer et al. (2013): $\ n(\gt D)\propto {D}^{-1.7}$, where D is the diameter of the object, because size distribution tends to follow a power-law distribution for a given population (Trilling et al. 2017). Extrapolating this model to $D\sim 100\,{\rm{m}}$ gives us an estimate of ∼10⁷ Centaurs.

Jedicke & Herron (1997) described the distributions of a, e, and i for Centaurs using analytic expressions. We normalize these expressions, as corrected by Grav et al. (2011), to the respective ranges for each parameter, forming the following probability distributions:

$$\begin{eqnarray}&&{P}_{\mathrm{Centaur}}(a)=\displaystyle \frac{1}{6.9\sqrt{2\pi }}\ \exp \ \left[-\displaystyle \frac{{\left(a-32\mathrm{au}\right)}^{2}}{2\cdot {\left(6.9\mathrm{au}\right)}^{2}}\right]\,{\mathrm{au}}^{-1}\,,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{Centaur}}(e)=2.26\ \exp \ \left[-\displaystyle \frac{{\left(e-0.21\right)}^{2}}{2\cdot {\left(0.21\right)}^{2}}\right]\,,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}{P}_{\mathrm{Centaur}}(i) & = & 3.96\times {10}^{-3}\ (i+2.7)\ \\ & & \times \exp \ \left[-\displaystyle \frac{{\left(i+2\buildrel{\circ}\over{.} 7\right)}^{2}}{2\cdot {\left(15^\circ \right)}^{2}}\right]\,{\deg }^{-1}.\end{array}\end{eqnarray} \tag{ 5 }$$

To mathematically describe the orbital element distributions of trapped interstellar objects, we fit Gaussian distributions to the segments of distributions that deviate most significantly from their Centaur counterparts. We describe $0\leqslant a\leqslant 10\,\mathrm{au}$, 0 ≤ e ≤ 1, and 0 ≤ i ≤ 180° with Gaussians, and express the resulting functions as normalized probability distributions:

$$\begin{eqnarray}&&{P}_{\mathrm{interstellar}}(a)=1.5\times {10}^{-1}\ \exp \ \left[-\displaystyle \frac{{\left(a-4\mathrm{au}\right)}^{2}}{2\cdot {\left(3\mathrm{au}\right)}^{2}}\right]\,{\mathrm{au}}^{-1},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{interstellar}}(e)=2.2\ \exp \ \left[-\displaystyle \frac{{\left(e-1\right)}^{2}}{2\cdot {\left(0.36\right)}^{2}}\right],\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{interstellar}}(i)=6.4\times {10}^{-3}\ \exp \ \left[-\displaystyle \frac{{\left(i-68^\circ \right)}^{2}}{2\cdot {\left(100^\circ \right)}^{2}}\right]\,{\deg }^{-1}.\end{eqnarray} \tag{ 8 }$$

The total number of ∼100-m-sized interstellar objects is estimated to be ∼6 × 10³ (Lingam & Loeb 2018), whereas the total number of ∼100-m-sized Centaurs is estimated to be ∼10⁷. The best way to distinguish between the populations is through their orbital inclination. We derive N_Centaur(i) and N_interstellar(i) as follows:

$$\begin{eqnarray}&&{N}_{\mathrm{Centaur}}(i)={10}^{7}\cdot {P}_{\mathrm{Centaur}}(i),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{N}_{\mathrm{interstellar}}(i)=6\times {10}^{3}\cdot {P}_{\mathrm{interstellar}}(i).\end{eqnarray} \tag{ 10 }$$

For the number of interstellar objects to dominate over Centaurs, we require that $\tfrac{{N}_{\mathrm{interstellar}}(i)}{{N}_{\mathrm{Centaur}}(i)}\gt 10$. We find that this condition holds for i ≳ 48°, and we integrate ${N}_{\mathrm{interstellar}}(i)\,{di}$ from 48° to 180°, resulting in an expected identifiable trapped population of ${N}_{\mathrm{interstellar}}\sim 4000$.

While the size distribution for interstellar objects is yet unconstrained, we will use the size distribution of Centaurs as an estimate, namely $n(\gt D)\propto {D}^{-1.7}$, where D is the diameter of the object, as Centaurs are the most dynamically similar population. Approximating $n(\gt D)\sim 6\,\times \,{10}^{3}$ for $D=100\,{\rm{m}}$ (Lingam & Loeb 2018), we find that $n(\gt D)\gtrsim 10$ for $D\gtrsim 4\,\mathrm{km}$.

We assume a very low albedo (.05) to find a conservative cutoff value for absolute magnitude of 4-km-sized objects, H > 16. We then filter all objects in the MPC database to those with values of semimajor axis, a, and eccentricity, e, within our established 25th and 75th percentile bounds, as well as inclination, i ≳ 48°, and absolute magnitude, H > 16. We set 5.5 au as the lower bound for semimajor axis, a, as our analysis only includes Centaurs. We find eight potential trapped interstellar objects, enumerated in Table 1.

Future surveys, such as the Large Synoptic Survey Telescope (LSST⁴ ), will be able to search for our predicted population of trapped interstellar objects. While the Hyper Suprime-Cam Survey (HSCS⁵ ) has comparable sensitivity, with the capability of reaching a magnitude of 24.5 with a 10-to-1 signal-to-noise ratio, LSST has a far larger field of view, covering 38 times more area on the sky, which will be more advantageous for the search. While HSC could reach a magnitude of 26 at a less demanding signal-to-noise constraint, the fact that the flat distribution of inclinations of trapped interstellar objects would not help with a targeted search, combined with HSC's relatively small field of view, makes LSST the preferred instrument for such a search in the near future.

An 'Oumuamua-sized object would have an apparent magnitude of 24.5, LSST's limit, at a distance of ∼2 au from the Sun. We will again use the size distribution of Centaurs as an estimate for that of trapped interstellar objects, namely $n(\gt D)\propto {D}^{-1.7}$, where D is the diameter of the object. The maximum distance at which an object can be observed scales as, $d\propto {D}^{\tfrac{1}{2}}$. Using the approximations, $d\sim 2\,\mathrm{au}$ and $n(\gt D)\sim 6\,\times \,{10}^{3}$, for $D=100\,{\rm{m}}$, we express n( > D) as the following function of d,

$$\begin{eqnarray}&&n(\gt D)\sim 6.3\times {10}^{4}\cdot {\left(\displaystyle \frac{d}{\mathrm{au}}\right)}^{-3.4}\,.\end{eqnarray} \tag{ 11 }$$

We then calculate the time-averaged distance of each particle in the population, $d=a(1+\tfrac{{e}^{2}}{2})$. Finally, we multiply each value of n(> D) by the proportion of the population that is at an average distance d, to obtain the number of objects expected to be detectable by LSST as a function of d, shown in Figure 2. We expect several hundreds of trapped interstellar objects to be detectable by LSST, with diameters ranging from ∼100 m to $\sim 10\,\mathrm{km}$. The stability lifetime of the trapped orbits is taken into account in this prediction, as it is based on the Lingam & Loeb (2018) treatment of the steady-state population of trapped objects. The lifetime over which an object will receive a "kick" from Jupiter is on the same order as the stability lifetime.

Zoom In Zoom Out Reset image size

We simulated the capture of interstellar objects, with initial speeds corresponding to the velocity distribution of stars in the Solar neighborhood, by means of three-body interactions with the Jupiter–Sun system. This was done by first populating the volume surrounding Jupiter with individual particles at their distances of closest approach to the planet. We then integrated the particle trajectories backward in time to determine which ones would be unbound from the Jupiter–Sun system, and subsequently integrated those particles forward in time to find which ones end in bound orbits. We then compared the resulting distributions of orbital elements to those of Centaurs to identify parameters for which interstellar objects should dominate, and therefore be easily identifiable. Using these results, we identified known objects that may be of interstellar origin. Finally, we estimated the number of trapped interstellar objects detectable by LSST.

Our calculations indicate that out of the trapped population of ∼6 × 10³ 'Oumuamua-sized interstellar objects, there should be ∼4000 trapped objects in our solar system identifiable by Centaur-like orbits at i ≳ 48°. We find that eight known objects have orbital parameters indicating their possible interstellar origin: 2018 WB1, 2013 YG48, 2018 TL6, 2017 SV13, 2013 JD4, 2018 AS18, 2008 WA95, and 2016 WS1, as listed in Table 1. Out of the entire population of trapped interstellar objects, we estimate that there are several hundreds of interstellar objects, ranging in diameter from ∼100 m to ∼10 km, detectable by LSST.

The orbit of asteroid (514107) 2015 BZ509, which was proposed to have an interstellar origin, is reproducible within the orbital parameter distributions of our trapped population (Namouni & Morais 2018). However, its orbital parameters are relatively unlikely, falling in the bottom quartile of both semimajor axis and eccentricity distributions. This could potentially be attributed to detection bias, as objects that are closer in are more likely to be detected by telescopes.

There is strong scientific motivation for investigating interstellar objects, including the potential to gain a deeper understanding of planetary system formation (Seligman & Laughlin 2018). Follow-up observations can help facilitate missions to probe objects of interstellar origin in this orbital parameter space. High-resolution spectroscopy could also be used to measure oxygen isotope ratios for objects with cometary tails in this orbital parameter space, as such ratios are expected to be markedly different for objects of interstellar origin compared to those that originated within the solar system (Lingam & Loeb 2018). The exploration of trapped interstellar objects could potentially help reveal the prospects of life in other star systems, as well as extraterrestrial artifacts (Freitas 1983; Haqq-Misra & Kopparapu 2012; Bialy & Loeb 2018; Lingam & Loeb 2019; Wright 2018). However, such a search would come with a caveat that most asteroids in the solar system reside outside the habitable zone.

We thank Martin Elvis, John Forbes, Manasvi Lingam, Roman Rafikov, and an anonymous referee for helpful comments on the manuscript. We also thank Peter Veres for helping us estimate the population of ∼100-m Centaurs in the solar system. This research has made use of data and/or services provided by the International Astronomical Union's Minor Planet Center. This work was supported in part by a grant from the Breakthrough Prize Foundation.