Hidden Planets: Implications from 'Oumuamua and DSHARP

In 2017 October, 1I/2017 U1, now 'Oumuamua, was identified by the Panoramic Survey Telescope and Rapid Response System (Chambers et al. 2016, Pan-STARRS) as the first interstellar interloper observed traversing the solar system (Meech et al. 2017). 'Oumuamua's measured eccentricity e = 1.1956 ± 0.0006 places its trajectory firmly in the regime of hyperbolic orbits. It was on its way out of the solar system with heliocentric distance d_H = 1.22 au upon its discovery, having passed perihelion with d_H = 0.25 au at closest approach (Meech et al. 2017). There is no clear consensus about the exact nature of 'Oumuamua, due in part to the short window of observability after its discovery and the corresponding sparsity of data, as well as the seemingly conflicting lines of evidence that have since emerged (Sekanina 2019).

Following the detection of 'Oumuamua, several authors estimated occurrence rates for analogous interstellar objects (hereafter, ISOs) free-floating in galactic space, with all estimates falling within the range n ≈ 0.004–0.24 au⁻³ deduced by Portegies Zwart et al. (2018). Assuming a cylindrical galaxy with 10¹¹ stars, R = 3 × 10⁴ pc, and H = 10³ pc as in Laughlin & Batygin (2017), even the most conservative of these estimates, n ≈ 0.004 au⁻³, results in at minimum 0.36M_⊕ of free-floating material per star, while the highest estimate, n ≈ 0.24 au⁻³, implies over 20M_⊕ per star. Throughout this work, we adopt the fiducial estimate n = 0.2 au⁻³ from Do et al. (2018), which carefully incorporates the volume probed by Pan-STARRS to provide the most robust estimate to date.

'Oumuamua's detection appears to require a high density of material ejected into interstellar space, though both the ejection mechanism and the origins of the body remain under debate (e.g., Bannister et al. 2019). No coma, nor any cometary molecular emission bands, was detected in association with 'Oumuamua (Meech et al. 2017; Trilling et al. 2018). Micheli et al. (2018), however, found that 'Oumuamua accelerated out of the solar system at a rate that could not be explained by gravitational forces alone, suggesting a cometary nature with acceleration induced by outgassing.

Icy material is ejected from stellar systems during multiple evolutionary stages. Early on, close interactions between stars forming in open clusters can liberate icy planetesimals from their circumstellar disks (Hands et al. 2019). Later, after giant planets have formed in the system, debris is ejected through dynamical interactions with these planets (e.g., Barclay et al. 2017; Raymond et al. 2018). Much of a star's enveloping Oort cloud is ultimately shed during the post-main-sequence stage of a star's lifetime (e.g., Veras et al. 2011; Do et al. 2018; Torres et al. 2019), releasing additional volatile-rich material.

In this paper, we quantify the rate of icy planetesimal ejection induced by gravity assists from circumstellar planets, with the aim of setting expectations for future observations from the Large Synoptic Survey Telescope (LSST; Abell et al. 2009). This exercise is critical to accurately assess the arrival prospects for future ISOs, each of which provides a uniquely valuable window into the nature of other planetary systems.

To accomplish this, we study planetesimal ejection rates from the long-period giant-planet population suggested by the Disk Substructures at High Angular Resolution Project (DSHARP; Andrews et al. 2018), which surveyed 20 nearby protoplanetary disks at high resolution using the Atacama Large Millimeter Array (ALMA). Considering the case in which 'Oumuamua stems from the population of ejected volatile-rich planetesimals, we show that gravity assists by long-period giant planets are capable of reproducing the number density of interstellar asteroids implied by the detection of 'Oumuamua. We use our results to predict the detection rate by LSST for ISOs produced through this channel.

Moro-Martín (2018) previously used multi-component power-law models to explore a wide range of possible ISO size distributions consistent with 'Oumuamua's appearance. Our estimates build on this previous work in several ways:

• 1.  
  Moro-Martín (2018) adopted relatively low giant-planet occurrence rates (f_pl = 0.2 for A-K2 stars and f_pl = 0.03 for K2-M stars) based on radial velocity (RV) measurements probing planets with semimajor axes a ≲ 3 au. The DSHARP survey, however, suggests an abundance of giant planets (f_pl ∼ 0.5) at wider separation a ≳ 5 au. We focus on these long-period planets, which are substantially more effective planetesimal ejectors than their shorter-period counterparts.
• 2.  
  Changes in the boundaries of a size distribution (r_min and r_max) can dramatically alter the resulting fit. In Moro-Martín (2018), these boundaries were set as constant values r_min = 1000 km and r_min = 1 μm, motivated by solar system models. We include no assumptions for the value of r_min and leave r_max as a free parameter within our models. We thus avoid assumptions about the degree of similarity between the size distributions of solar system KBOs and free-floating ISOs, which may originate from a range of collisional histories that do not necessarily resemble that of the Kuiper Belt.
• 3.  
  Lastly, Moro-Martín (2018) assumes that 100% of solid circumstellar material around single stars and wide binaries is ejected. We instead perform N-body simulations of several representative systems from the DSHARP survey to deduce the expected mass of ejected material.

We present our simulation setup in Section 2 and subsequent results in Section 3, including our final range of possible power-law solutions. We then discuss implications of our findings for future detections in Section 4 and sources of uncertainty in Section 5 prior to concluding.

2.1. Planets as Interstellar Comet Ejectors

Ejection of planetesimals is a natural outcome of close encounters with Jupiter, as observed in N-body simulations of early solar system evolution (e.g., O'Brien et al. 2006; Levison et al. 2009). Yet, not all planets are capable of efficiently ejecting material from their circumstellar systems. To readily expel material through gravity assists, a planet must have a Safronov number Θ ≳ 1, with Θ given by

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{{v}_{\mathrm{esc}}^{2}}{{v}_{\mathrm{orb}}^{2}}=\displaystyle \frac{{M}_{p}a}{{R}_{p}{M}_{* }}.\end{eqnarray} \tag{ 1 }$$

Here, M_p, R_p, and a are the mass, radius, and semimajor axis of a planet, respectively, while M_* is the mass of the host star. Notably, hot Jupiters and super Earths cannot efficiently eject material from their systems. Long-period (a ≳ 5 au) planets of Neptune's mass or greater are the most effective ejectors; these planets, however, lie in a region of parameter space that is heavily disfavored by the detection biases of the transit and RV methods. Direct imaging results show that occurrence rates for giant planets with orbital separation 10–100 au are low, at ${9}_{-4}^{+5} \%$ between 5–13 M_Jup (Nielsen et al. 2019), indicating that the primary planetary ejectors may be less massive than the current direct imaging detection limits and/or lie within ∼10 au of the host star.

Recently, a candidate population of such planets was inferred by Zhang et al. (2018) using results from DSHARP. Disk substructures were nearly universally found in the DSHARP sample, where axisymmetric gaps and rings are most common. The DSHARP sample is biased toward relatively large disks around massive stars; however, gaps and rings at comparable radii have also been identified in samples more representative of the average protoplanetary disk population (e.g., Long et al. 2018), suggesting that they may be prevalent among the underlying population. Although these substructures could be induced by a number of mechanisms, mounting evidence for planetary companions in similar systems (e.g., PDS 70; Haffert et al. 2019) favors the hypothesis that the substructures in the DSHARP sample are caused by planets. Recent kinematic detections also point toward the presence of planets at radial locations coincident with observed gaps in protoplanetary disks (Pinte et al. 2019), further strengthening the case that some of the observed substructures are indeed induced by planets.

Consequently, the DSHARP sample suggests a ∼50% occurrence rate for giant planets with masses between ∼M_Nep and a few M_Jup orbiting their host stars at separation a ≳ 5 au (Zhang et al. 2018). This abundance of long-period giant planets is roughly in agreement with previous results from Bryan et al. (2016), which combined Keck RV measurements with NIRC2 adaptive optics imaging to obtain an occurrence rate of 52% ± 5% for planets with M = 1–20M_Jup and a = 5–20 au. An abundance of Neptune-mass planets is also supported by microlensing results, which find a peak in planet occurrence rates at a planet-to-star mass ratio q ∼ 10⁻⁴, corresponding to roughly 20M_⊕, or 1.2 M_Nep, for a typical host-star mass ∼0.6M_⊙ (Suzuki et al. 2016).

2.2. Simulation Setup

To determine how this planet population connects to the population of interstellar asteroids, we complete a detailed assessment of mass ejection rates in three DSHARP disk systems with clear radial gaps—HD 143006, AS 209, and HD 163296—using the REBOUND orbital integrator (Rein & Liu 2012). We adopt host-star masses and radial planet positions from Table 3 of Zhang et al. (2018). For each predicted planet, Zhang et al. (2018) reports several masses derived with varying assumptions for the disk dust size distribution and α viscosity parameter. Because these model parameters are not well-constrained, we randomly assign the mass of each planet from the range [M_min, M_max], where M_min and M_max are the minimum and maximum derived masses, respectively. In accordance with results from Zhang et al. (2018), the HD 143006 and AS 209 systems each include two planets, while the HD 163296 system includes three planets.

We focus on the millimeter-sized dust population in our simulations, since dust masses in this size regime are relatively well-constrained by ALMA. We estimate the total millimeter dust mass M_dust of each disk using the relation

$$\begin{eqnarray}&&\mathrm{log}({M}_{\mathrm{dust}}/{M}_{\oplus })=1.3\,\times \,\mathrm{log}({M}_{* }/{M}_{\odot })+1.1,\end{eqnarray} \tag{ 2 }$$

derived in Pascucci et al. (2016), where M_* is the mass of the host star. This derivation is based on the 887 μm flux measured from observations of the Chamaeleon I star forming region, with bandwidth spanning 867–917 μm. Thus, the dust masses in our simulations represent only the subset of grains with radii in this size range.

To reproduce the dust mass distribution in each system, we initialize 3500 equal-mass test particles radially distributed using the semimajor axis probability distributions depicted in Figure 1. These distributions are based on radially symmetric best-fitting models consisting of several superposed Gaussians for each disk, with parameters for AS 209, HD 163296, and HD 143006 given in Guzmán et al. (2018), Isella et al. (2018), and Pérez et al. (2018), respectively. Each planetesimal is initialized with orbital elements

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{\Omega }} & = & \omega =0\\ f & \in & [-\pi ,+\pi ]\\ i & \in & [0.0005,0.005]\\ e & \in & [0.04,0.06],\end{array}\end{eqnarray*}$$

where all angles are in radians. We initialize each of our three disk-star-planet systems in three iterations each with a different random seed, resulting in a total of nine simulations with varying planet masses, orbital elements, and initial dust particle placements. We integrate each of these systems for one week on the Grace supercomputing cluster at Yale University.

Zoom In Zoom Out Reset image size

To encapsulate the typical behavior of all systems, we find the average millimeter mass m_ej ejected by each star as a function of time t for all nine simulations, with our results displayed in Figure 2. The resulting curve, representing the total mass of all particles with a formally positive energy relative to the host star, is well-fit by the exponential function

$$\begin{eqnarray}&&{m}_{\mathrm{ej}}=a{\mathrm{log}}_{10}t+b,\end{eqnarray} \tag{ 3 }$$

with a = 0.468 and b = −2.274. Using this profile, we extend our results to t = 10⁸ yr, at which point m_ej = 1.47M_⊕. We take this as the representative millimeter-sized dust mass ejected by each giant-planet-hosting star throughout its main-sequence lifetime, selecting a relatively late time t to compensate for mass ejected prior to the gas-clearing phase in each system. Adopting a 50% occurrence rate for this population of planets, we approximate that the average rate of millimeter-sized planetesimal ejection is 0.74 M_⊕ star⁻¹.

Zoom In Zoom Out Reset image size

In order to connect the mass ejected per star to an ISO spatial number density n, we must first determine the background stellar number density n_*. We adopt the local midplane stellar mass density ρ_* = 0.043 ± 0.004 M_⊙ pc⁻³ found in McKee et al. (2015), then estimate the average stellar mass in the solar neighborhood using the sample of 75 known stars (including brown dwarfs) within 5 pc of the Sun. For stars without measured masses, we apply the piecewise mass–luminosity relation from Eker et al. (2015). We obtain an average stellar mass M ∼ 0.35M_⊙, which translates to stellar number density n_* ≈ 0.12 stars pc⁻³.

Combining n_* with our average planetesimal mass ejection rate per star, we find projected ISO mass density ρ_ISO ∼ 0.09 M_⊕ pc⁻³ free-floating in space. This means that, for number density distribution dn/dr with boundaries r_min,obs = 867 μm and r_max,obs = 917 μm,

$$\begin{eqnarray}&&{\int }_{{r}_{\min ,\mathrm{obs}}}^{{r}_{\max ,\mathrm{obs}}}\displaystyle \frac{{dn}}{{dr}}m(r){dr}=0.09\,{M}_{\oplus }\,{\mathrm{pc}}^{-3}.\end{eqnarray} \tag{ 4 }$$

We use $m(r)=\tfrac{4}{3}\pi {r}^{3}\rho$ with ρ = 1 g cm⁻³ and fit dn/dr with a power-law radius distribution of the form

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dr}}={{Cr}}^{-q},\end{eqnarray} \tag{ 5 }$$

where q is the power-law index and C is constant. We interpret previous estimates of n obtained from the detection rate implied by 'Oumuamua as n(r ≳ r_O), where r_O is 'Oumuamua's effective radius, as larger ISOs would also be detectable when passing through the solar system at 'Oumuamua's distance. Equivalently,

$$\begin{eqnarray}&&{\int }_{{r}_{{\rm{O}}}}^{{r}_{\max }}\displaystyle \frac{{dn}}{{dr}}{dr}=n.\end{eqnarray} \tag{ 6 }$$

Combining and integrating Equations (4)–(6), we parameterize our solution space with the function

$$\begin{eqnarray*}&&f=\displaystyle \frac{4}{3}\pi \rho n\left(\displaystyle \frac{1}{0.09{M}_{\oplus }\,{\mathrm{pc}}^{-3}}\right)\left(\displaystyle \frac{1-q}{4-q}\right)\left[\displaystyle \frac{{r}_{\max ,\mathrm{obs}}^{4-q}-{r}_{\min ,\mathrm{obs}}^{4-q}}{{r}_{\max }^{1-q}-{r}_{{\rm{O}}}^{1-q}}\,\right],\end{eqnarray*}$$

which is equal to unity when Equations (4)–(6) are simultaneously satisfied. Using f, we solve for power-law indices q over a range of r_max, with r_O = 55 m from Jewitt et al. (2017).

The resulting solutions are displayed in Figure 3, indicating that, for a given n, our power-law solutions are robust to changes in r_max for an ISO distribution in which ${r}_{\max }\,\rlap{/}{\approx }\,{r}_{O}$. We take a vertical slice through Figure 3 at r_max = 10⁶ cm to display q as a function of n in Figure 4. Ultimately, we find that high power-law indices q > 4 are required to fit both the n values implied by 'Oumuamua and our simulated mass ejection rates. Adopting q = 4.09, this distribution corresponds to a total mass of ejected ISO material ∼24M_⊕ per DSHARP-type system for ISOs in the size range 10⁻³ cm ≤ r ≤ 10⁶ cm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another potential constraint for this distribution is the number density ${n}_{\mathrm{CNEOS}}={10}^{{6}_{-1.5}^{+0.75}}$ au⁻³ of r_CNEOS ≳ 0.45 m ISOs implied by the candidate interstellar meteor identified from the Center for Near-Earth Object Studies (CNEOS) bolide catalog in Siraj & Loeb (2019). We note that Devillepoix et al. (2019) reported that the United States government sensors used to discover this object are generally unreliable for orbit calculations. We include the candidate CNEOS meteor for completeness and primarily focus on determining whether the object is consistent with our analysis up to this point. If the CNEOS meteor originates from the same underlying distribution as 'Oumuamua, it would require

$$\begin{eqnarray}&&{\int }_{{r}_{\mathrm{CNEOS}}}^{{r}_{\max }}\displaystyle \frac{{dn}}{{dr}}{dr}={n}_{\mathrm{CNEOS}}.\end{eqnarray} \tag{ 7 }$$

Adopting q = 4.09, corresponding to n = 0.2 au⁻³ from Do et al. (2018), we obtain the number density n = 5.63 × 10⁵ au⁻³ for objects with r_CNEOS ≤ r ≤ r_max. This value is firmly within the error bars for n_CNEOS from Siraj & Loeb (2019), indicating that our results are consistent with those implied by the candidate CNEOS meteor.

Although the first ISO was not confirmed until 2017, anticipatory estimates for the detectability of such objects have been made on several earlier occasions. To account for the prior non-detection of ISOs, early assessments predicted an interstellar comet number density n ∼ 10¹³ pc⁻³ for objects with r ≥ 1 km (McGlynn & Chapman 1989). Later, Jewitt (2003) projected that if all stars eject 10¹³ comets with r ≥ 1 km, approximately 0.3 such interstellar interlopers should pass within 5 au of the Sun—the approximate detection limit of Pan-STARRS—each year.

Moro-Martín et al. (2009) followed up this study by calculating the expected detection rate by LSST of inactive extrasolar comets passing through the solar system at distances greater than 5 au. They found that the probability of LSST detecting an inactive interloper during its 10 yr lifetime is ∼0.01%–1%, and they ultimately concluded that LSST will likely not observe even one such object. Cook et al. (2016) then updated prior estimates by numerically simulating the detectability of interstellar asteroids within 5 au of the Sun, concluding that an optimistic estimate would result in 1 interstellar object detected during LSST's lifetime. Since 'Oumuamua's discovery with Pan-STARRS, Trilling et al. (2017) used a scaling argument to project that the detection rate of analogous interlopers will increase from 0.2 yr⁻¹ with Pan-STARRS to 1 yr⁻¹ with the advent of LSST.

Leveraging the power-law radius distributions obtained in Section 3, we independently estimate the expected detection rate by LSST for protostellar disk-ejected ISOs. Adopting a single-frame magnitude limit m_LSST ∼ 24 (LSST Science Collaboration et al. 2009) and cometary geometric albedo p_V = 0.04 (Lamy et al. 2004), we calculate the smallest ISO radius observable with LSST at a given distance d, given in au, using

$$\begin{eqnarray}&&{r}_{\min ,\mathrm{LSST}}=\displaystyle \frac{1\,\mathrm{au}}{\sqrt{{p}_{v}}}\,{10}^{\tfrac{1}{5}({m}_{\odot }-H)},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&H={m}_{\mathrm{LSST}}-5{\mathrm{log}}_{10}d,\end{eqnarray} \tag{ 9 }$$

where m_⊙ = −26.74 is the magnitude of the Sun and H is the absolute magnitude of the ISO. Our calculation implicitly assumes that all ISOs are observed at opposition; however, in practice the solar angle and frequency of sky coverage will also play an important role in the final ISO detection rate, meaning that our rates should be treated as upper limits.

At each of 10⁵ evenly spaced circular annuli around the Earth, from 0–40 au, we sum over the number density profile from r_min,LSST to r_max to obtain the total observable ISO number density n_ISO within that annulus. Combining n_ISO with the surface area of each annulus σ, we calculate the total detection rate of ISOs as

$$\begin{eqnarray}&&R=\sum {n}_{\mathrm{ISO}}\sigma v,\end{eqnarray} \tag{ 10 }$$

summing over all annuli. We take 'Oumuamua's velocity at infinity, 26 km s⁻¹ (Mamajek 2017), as our representative v. Completing this calculation with reference minimum radius thresholds r > 1 m, r > 10 m, and r > 55 m, we obtain the results outlined in Figure 5.

Zoom In Zoom Out Reset image size

We find that, if 'Oumuamua is an ejected cometary planetesimal drawn from an isotropic power-law distribution, LSST will find several ISOs each year of 'Oumuamua's size or larger, as well as up to hundreds of ISOs per year with r > 1 m.

Zhu et al. (2019) recently found that protoplanetary disk masses may be heavily underestimated due to the assumption that the disks are optically thin. As a result, in the case that these disks are optically thick, our millimeter mass ejection rates may be underestimated, leading to an underestimate in the power-law slope.

However, it may instead be the case that our disk masses are overestimated. Our three sample disks are located around relatively massive stars (0.83M_⊙, 1.78M_⊙, and 2.04M_⊙ for AS 209, HD 143006, and HD 163296, respectively), meaning that treating these systems as "typical" may instead lead to an overestimate of the total mass of ejected material. This would accordingly correspond to an overestimate of the LSST detection rate, particularly on the lower end of the ISO size distribution.

Furthermore, the number density of ISOs, n, may be locally enhanced at the present time due to the Sun's current proximity to the Galactic midplane: z_⊙ ∼ 17 ± 5 pc (Karim & Mamajek 2017), as compared to its maximum height above the midplane 49–93 pc (Bahcall & Bahcall 1985). If the Sun is currently passing through a low-dispersion, "cold" population of ISOs, the apparent number density as observed from the solar system may be higher than the bulk average density throughout the Galaxy. Because we do not yet have strong constraints on the scale height of ISOs, we cannot rule out this possibility, which suggests that the true n may be lower than most current estimates. This could result in a steeper power-law size distribution and a correspondingly lower LSST detection rate for large ISOs.

It is also quite possible that the distribution of interstellar comets does not follow a simple power law. Moro-Martín (2018), after sampling a wide range of possible models, found that 'Oumuamua is likely not representative of an isotropic background population. Though our present model is already overconstrained and therefore does not warrant the inclusion of additional free parameters, we acknowledge that a power-law fit may not fully capture the true characteristics of the ejected ISO population. We may be overestimating the number density of large ISOs if, for example, the collisional evolution of solid material in extrasolar systems typically grinds the vast majority of larger bodies into dust prior to ejection.

Lastly, the detection rates reported in Section 4 may be modestly increased by gravitational focusing and comet brightening from passage near the Sun, each of which enhances the observability of the ISO population.

The discovery of the first interstellar interloper, 'Oumuamua, has already provided exceptional insight to inform our understanding of planetary systems while simultaneously presenting new puzzles. In our work, we have reconciled simulations of the observed DSHARP planet-disk systems with the detection of 'Oumuamua to constrain the range of possible size distributions for ISOs ejected through interactions with circumstellar giant planets. We conclude that the population of long-period giant planets suggested by the DSHARP sample is capable of ejecting the population of free-floating planetesimals implied by 'Oumuamua's detection. Furthermore, 'Oumuamua is consistent with an origin within a population of ISOs following a single power-law radius distribution, resulting in an anticipated LSST detection rate ranging from just a few detections per year for 'Oumuamua-sized ISOs to over 100 yr⁻¹ for ISOs with r > 1 m. Future observations of ISOs could potentially answer long-standing questions about not only the range of processes taking place in extrasolar systems, but also the population statistics of long-period giant exoplanets, the diversity of small-body populations throughout the Galaxy, and the evolutionary path of the solar system itself.

We thank the anonymous referee for helpful comments that substantially improved the manuscript. We also thank Matt Holman for insightful inputs as the manuscript was being finalized. M.R. is supported by the National Science Foundation Graduate Research Fellowship Program under grant No. DGE-1752134. This material is also based upon work supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under Cooperative Agreement Notice NNH13ZDA017C issued through the Science Mission Directorate. We acknowledge support from the NASA Astrobiology Institute through a cooperative agreement between the NASA Ames Research Center and Yale University. Simulations in this paper made use of the REBOUND code, which can be downloaded freely at http://github.com/hannorein/rebound. This research also made use of the numpy (Oliphant 2006), matplotlib (Hunter 2007), and seaborn (Waskom et al. 2014) Python packages.