Misaligned Circumbinary Disks as Efficient Progenitors of Interstellar Asteroids

Circumbinary gas disks (CBDs) with large misalignments relative to the binary orbital plane are commonly observed in nature (e.g., Chiang & Murray-Clay 2004; Köhler 2011; Andrews et al. 2014; Brinch et al. 2016; Takakuwa et al. 2017; Fang et al. 2019; Kennedy et al. 2019; Kenworthy et al. 2022; Zhu et al. 2022). The degree of CBD misalignment often increases with binary separation and eccentricity (Czekala et al. 2019). Misalignments may initially arise as a result of turbulence in the molecular gas cloud (Offner et al. 2010; Bate 2012; Tokuda et al. 2014), later accretion of material by the young binary (Bate et al. 2010; Bate 2018), warping by a tertiary companion such as a stellar flyby (Nealon et al. 2020), or if the binary forms from a cloud whose elongated axis is misaligned to its rotation axis (Bonnell & Bastien 1992). The misaligned disk may precess as a solid body if the communication timescale is shorter than the precession timescale (Papaloizou & Terquem 1995; Larwood et al. 1996). As a result of dissipation, a viscous disk evolves toward either a coplanar or polar (90°) alignment to the binary orbital plane (Martin & Lubow 2017; Lubow & Martin 2018; Zanazzi & Lai 2018; Cuello & Giuppone 2019) although, depending on the binary and disk parameters, the timescale for alignment may be longer than the disk lifetime, meaning that planet formation can take place in misaligned disks (e.g., Martin & Lubow 2018).

The late stage of terrestrial planet formation takes place after Moon-size planetesimals and Mars-size embryos have formed and the gas disk has dispersed. These solid bodies interact with one another through purely gravitational interactions to form terrestrial planets through core accretion (Artymowicz 1987; Lissauer 1993; Pollack et al. 1996). Coplanar and polar circumbinary orbits are stationary states in which the particles do not undergo significant nodal precession. As a result, terrestrial planets can efficiently form in coplanar (Quintana & Lissauer 2006; Childs & Martin 2021a) and polar aligned (Childs & Martin 2021b) circumbinary disks through core accretion. Such terrestrial circumbinary planets (CBPs) have yet to be observed, however. While this may be attributed to observational bias against such small planets in a circumbinary orbit (Windemuth et al. 2019; Martin & Fabrycky 2021; Standing et al. 2022), it may also indicate that terrestrial planets do not form through core accretion in a circumbinary disk or that current core accretion models are missing key physics.

In a disk that is misaligned from a stationary state, nodal precession can lead to large mutual misalignments and collisions with large impact velocities that may result in ejection from the system. Whether planets can form or not depends upon the misalignment and the binary eccentricity. Terrestrial CBPs that do form end up either coplanar or polar to the binary orbital plane since mergers between bodies with random nodal phase angles lead to lower inclination to the stationary states (Childs & Martin 2022). While collisions were resolved with only perfect merging in Childs & Martin (2022), in this work we consider fragmentation as a more realistic outcome from such high-energy collisions.

'Oumuamua was the first confirmed interstellar asteroid (ISA) to be observed (Chambers et al. 2016). 'Oumuamua does not exhibit comet-like features indicating that its composition is more consistent with a refractory planetoid (Jewitt et al. 2017; Meech et al. 2017; Ye et al. 2017). This ISA has an unexpectedly low velocity relative to the local standard of rest, ∼10 km s⁻¹ (Mamajek 2017), and has a highly elongated shape (Meech et al. 2017; Bolin et al. 2018). The elongated shape and tumbling motion of this body suggest that it was involved in a violent collision in its past and was sent tumbling in its parent planetary system, indicating that collisions of solid bodies in other planetary systems are not uncommon (Drahus et al. 2017; Fraser et al. 2018). The currently observed mass of 'Oumuamua is estimated to be about 10⁻¹⁷ M_⊕; however, if 'Oumuamua is composed of entirely N₂ ice, it could have lost up to 92% of its initial mass upon entering the solar system (Desch & Jackson 2021; Jackson & Desch 2021). Seligman & Laughlin (2020) proposed that if 'Oumuamua contained a significant amount of H₂ ice, it was likely pancake-shaped when it was near periapsis, and Mashchenko (2019) found that the light curve is consistent with such a shape. In these cases, the observed properties of 'Oumuamua are not representative of its origins.

Various formation scenarios for 'Oumuamua and other ISAs have been proposed such as ejections from a system as a result of tidal interactions with a white dwarf (Rafikov 2018), ejections of fragments from tidally disrupted planets by a dense member of a binary system (Ćuk 2018), and ejections of a comet-like planetesimal from giant-planet interactions (Raymond et al. 2018). Binary stars have been suggested to be more likely than single stars to cause the ejection of a rocky body from planetary systems (Jackson et al. 2018). Planetesimals are ejected when they migrate inside the stability limit of the binary although this may require the presence of other planets (Fitzmaurice et al. 2022).

In this letter, we propose that, shortly after a highly misaligned circumbinary gas disk dissipates, solid bodies undergo violent collisions and become a source for ISAs such as 'Oumuamua. The highly inclined particles have no requirement to migrate close to the binary because ejections occur over a wide radial range. In Section 2 we first conduct three-body simulations to show how the initial disk misalignment and the binary eccentricity affect the particles' mutual inclinations, and thus impact velocities. In Section 3 we then conduct n-body studies of terrestrial CBP formation in highly misaligned CBDs and resolve collisions with fragmentation. We closely follow the collisions and the fate of the ejected material to better understand the nature of ISAs that are generated from misaligned CBDs. In Section 4 we discuss the implication of these results for ISAs. Lastly, we conclude with a summary of our findings in Section 5.

A particle in a circumbinary orbit around an eccentric binary can undergo two types of nodal precession depending upon its initial inclination. For low initial tilt, the orbit is circulating, meaning that the particle orbit precesses around the binary angular momentum vector. If the initial inclination is above the critical value, the orbit will be librating, meaning that it precesses about the binary eccentricity vector (Verrier & Evans 2009; Farago & Laskar 2010; Doolin & Blundell 2011; Aly et al. 2015). The critical inclination depends upon the binary eccentricity and the angular momentum of the particle (Chen et al. 2019; Martin & Lubow 2019). In the test particle limit, the minimum critical inclination that separates circulating and librating orbits is given by

$$\begin{eqnarray}&&{i}_{\mathrm{crit}}={\sin }^{-1}\sqrt{\displaystyle \frac{1-{e}_{{\rm{b}}}^{2}}{1+4{e}_{{\rm{b}}}^{2}}}\end{eqnarray} \tag{ 1 }$$

(Farago & Laskar 2010). This critical inclination occurs for a longitude of ascending node of ϕ = 90° measured in the frame of the binary (see Equation (3) in Chen et al. 2019).

We measure the particle misalignment with respect to one of the two stable configurations, coplanar or polar. The inclination of the particle orbit relative to the binary orbit is given by

$$\begin{eqnarray}&&{i}_{{\rm{b}}}={\cos }^{-1}({\hat{{\boldsymbol{l}}}}_{{\rm{b}}}\cdot {\hat{{\boldsymbol{l}}}}_{{\rm{p}}}),\end{eqnarray} \tag{ 2 }$$

and the inclination of the particle orbit relative the binary eccentricity vector is given by

$$\begin{eqnarray}&&{i}_{{\rm{e}}}={\cos }^{-1}({\hat{{\boldsymbol{e}}}}_{{\rm{b}}}\cdot {\hat{{\boldsymbol{l}}}}_{{\rm{p}}}),\end{eqnarray} \tag{ 3 }$$

where l _b and l _p are the angular momentum vectors of the binary and particle, respectively, e _b is the eccentricity vector of the binary, and $\hat{{\boldsymbol{}}}$ denotes a unit vector. For orbits with ϕ = 90° initially, if the initial particle inclination is smaller than the critical inclination (circulating orbit), we measure i_b, and if the particle inclination is larger than the critical inclination (librating orbits) we measure i_e.

The maximum impact velocity for a collision between two particles orbiting at the same semimajor axis in circular orbits with Keplerian velocity, v_K, can be estimated with

$$\begin{eqnarray}&&{v}_{max}=2{v}_{{\rm{K}}}\sin \left({i}_{max}/2\right),\end{eqnarray} \tag{ 4 }$$

where ${i}_{\max }=\max \,{i}_{{\rm{m}}}$ is the maximum value over a nodal precession period of the mutual inclination between the two particles, i_m. In a colliding system with particles of different nodal precession rates, i_max is twice the maximum inclination a particle reaches over its nodal precession period measured with respect to the stationary inclination about which it precesses (either coplanar or polar).

To probe the maximum mutual inclinations expected in a circumbinary disk as a function of binary eccentricity we conduct three-body simulations of a very low-mass particle at an orbital radius of 5 a_b from the barycenter of the three-body system. We change the eccentricity and initial inclination of the particle to the binary, i_b0, and integrate for two full nodal precession periods (see Equations (6)–(10) of Childs & Martin 2022) using the WHFAST integrator in rebound (Rein & Liu 2012). Initially the longitude of ascending node is ϕ = 90° in all cases.

Figure 1 shows the results from our three-body simulations. We plot the maximum inclination a single particle experiences relative to the axis about which it precesses (coplanar or polar), which is equivalent to i_m/2, as a function of the particle's initial inclination to the binary, i_b0, and binary eccentricity, e_b. We plot the critical inclination (Equation (1)) as the solid black line. We see that as the binary eccentricity increases, the less inclined the particle needs to be to reach the maximum inclination away from a stable configuration. This is expected since the maximum inclination a particle experiences is near the critical inclination, which decreases as binary eccentricity increases. This indicates that CBDs with relatively modest inclinations around highly eccentric binaries can harbor particles with high mutual inclinations, and thus colliding particles can experience high impact velocities. We do not expect this trend to change for binary systems with different mass ratios or different separations since these effects change only the timescale for the particle dynamics. The critical inclination and the particle dynamics depend only on the binary eccentricity and angular momentum ratio of the particle to the binary (Farago & Laskar 2010; Chen et al. 2019; Martin & Lubow 2019).

Zoom In Zoom Out Reset image size

We now explore simulations of terrestrial planet formation in the inner parts of an initially misaligned circumbinary disk including the effects of fragmentation. The fragmentation code we use is detailed in Childs & Steffen (2022). We distribute 26 Mars-sized embryos (m ≈ 0.1 M_⊕) and 260 Moon-sized planetesimals (m ≈ 0.01 M_⊕) along the fits from smoothed particle hydrodynamics simulations detailed in Childs & Martin (2021a). This bimodal mass distribution is adopted from n-body simulations that were successfully able to recover the masses of the terrestrial planets in the solar system (Chambers 2001).

Unlike WHFAST, the IAS15 integrator in rebound is a high-precision non-symplectic integrator that is able to resolve close encounters and collisions (Rein & Spiegel 2015). This feature is necessary for modeling core accretion of the terrestrial planets. To overcome the excessive CPU time associated with a non-symplectic integrator we apply an expansion factor of f = 25 to all the particles after integrating the system for 100 kyr and the phase angles of the particles have randomized. Childs & Martin (2021a) performed convergence tests with f = 25 and smaller expansion factors in n-body simulations of CBP terrestrial planet formation using the IAS15 integrator. They found that while larger expansion factors lead to some differences in system architecture, the general planet formation trends that emerge as a function of binary eccentricity and separation remain. Furthermore, Childs & Steffen (2022) studied the effects of expansion factors with fragmentation. Their findings indicate that our use of an expansion factor will lead to shorter collision timescales and more damped orbits, which are more likely to underestimate impact velocities in collisions.

We set the minimum fragment mass to half the size of the Moon (m ≈ 0.005 M_⊕). Fragments are expected to be much smaller, but we choose this value to reduce the CPU time of the simulations. While this fragment mass is orders of magnitude larger than that of 'Oumuamua, it should be viewed as an upper limit because fragment-producing collisions will produce a wide distribution of fragment sizes (Leinhardt & Stewart 2012).

The orbital elements for each body are randomly chosen in each run. We consider perfectly circular binaries with e = 0.0 and highly eccentric binaries with e = 0.8 with particle disks that are initially inclined 30° and 60° above the binary orbital plane. A particle is considered ejected from the system once its distance from the barycenter of the system exceeds 100 au. The binary consists of equal-mass stars with a total binary mass of 1 M_⊙ separated by 0.5 au. We perform 50 runs for each setup and integrate for a total of 7 Myr.

The different binary models and their corresponding binary eccentricities are listed in Table 1 and are marked by red triangles in Figure 1. The initial surface density profile is taken to be that of a steady circumbinary gas disk as described in Childs & Martin (2021b). At least initially, both particle disks around the circular orbit binary are in circulating orbits while around the eccentric orbit binary they are librating.

Our simulations reveal that the models C60 and E30 eject the most material and do so throughout the entirety of the simulation. The sustained ejection rates indicate that particles are not just being quickly ejected in the inner, unstable region of the disk close to the binary as a result of strong binary–particle interactions (e.g., Holman & Wiegert 1999; Chen et al. 2020). Particles found in larger, initially more stable, orbits also get ejected on longer timescales as a result of particle–particle interactions. On average, the C60 and E30 systems eject 80% of their disk mass and the C30 and E60 systems eject 1% of their disk mass by the end of our simulations. The large difference between these ejection percentages is the result of the different mutual inclinations the particles reach in the runs. In Figure 1 we see that the C60 and E30 particles can reach maximum mutual inclinations of i_m = 120° and i_m = 132°, respectively, for circular orbits. Such large mutual inclinations will result in high impact velocities when a collision takes place which is likely to result in mass ejection. The C30 and E60 particles reach maximum inclinations of i_m = 60° and i_m = 64°, respectively, which will result in much lower impact velocities and less mass being ejected from the system.

Because of the high ejection rates in C60 and E30, no terrestrial planets are formed. On average, the C30 and E60 systems form at least one terrestrial planet with a mass greater than 1 M_⊕ in an orbit that is nearly circular and coplanar to the circular binary or polar to the eccentric binary.

The final planetary systems with fragmentation are similar to the final planetary systems of Childs & Martin (2022), who modeled planet formation in similar misaligned CBDs but resolved collisions with only perfect merging. They did not consider a system analogous to E60, but our E60 planetary systems closely resemble those formed in the C30 runs since both are inclined by 30° to a stationary inclination. The planet eccentricities and inclinations are slightly damped, relative to the planets formed with only perfect merging, due to dynamical damping from the fragments. A notable difference between the planetary systems that formed with perfect merging and with fragmentation is the formation timescale. In agreement with Quintana et al. (2016), who compared planet formation in systems with and without fragmentation, fragmentation approximately doubles the formation and CPU time but results in similar planetary systems to those formed with only perfect merging.

The left column of Figure 2 shows the impact velocity versus the semimajor axis for all the collisions that took place in one run of each model. The pluses mark collisions that lead to ejections. We also plot the curves that represent the maximum impact velocity for a circular orbit, given by Equation (4), for each model. We use the maximum inclination the particle reached in our three-body simulations for i_m/2 to predict the maximum impact velocity in systems with the same binary setup. We see that the majority of the collisions in each model are near the impact velocity curves for a circular orbit. Particles that are on eccentric orbits can result in larger impact velocities than v_max. We see that collisions that lead to ejections can be found for a wide range of velocities but are typically found when the impact velocity is greater than ∼15 km s⁻¹.

Zoom In Zoom Out Reset image size

When the kinetic energy (KE) of a body is equal to or less than the potential energy (PE) of the body, it will remain gravitationally bound to the binary. The right column of Figure 2 shows the KE = PE line for a range of masses and radii as a black dashed line. We expect collisions on and below this line to remain bound to the binary and collisions above this line to be ejected from the system. This line also corresponds to ${v}_{i}=2{v}_{{\rm{K}}}\sin {45}^{\circ }$. Using data from one run in each setup, we plot the PE and KE of each collision using the total mass of the colliding system, the impact velocity, and the last recorded semimajor axis of the target. Collisions that involve bodies that are eventually ejected from the system are marked with a plus. We also plot the line with slope KE/PE using the maximum impact velocity for a circular orbit for each system, over a range of radii and mass. CBDs with particles with 45° ≤ i_m/2 will result in collisions where $1\leqslant {\rm{KE}}/{\rm{PE}}$, which will lead to ejection from the system unless dynamics from the multibody system prevent this. We see that the KE versus PE line for the C60 and E30 system, which have a maximum impact velocity for a circular orbit of ${v}_{{\rm{i}}}=2{v}_{{\rm{K}}}\sin {60}^{\circ }$ and ${v}_{{\rm{i}}}=2{v}_{{\rm{K}}}\sin {66}^{\circ }$ respectively, is above the KE = PE line where ∼80% of ejections are found. The C30 and E60 systems, which have a maximum impact velocity for a circular orbit of ${v}_{{\rm{i}}}=2{v}_{{\rm{K}}}\sin {30}^{\circ }$ and ${v}_{{\rm{i}}}=2{v}_{{\rm{K}}}\sin {32}^{\circ }$ respectively, are found below the KE = PE line, explaining the low number of ejections we observe in these systems. These i_m/2 values are taken from the corresponding binary systems in Figure 1, which are marked by red triangles.

The mass of the largest remnant of a collision that results in fragmentation is M_lr. The formula for calculating M_lr is taken from Leinhardt & Stewart (2012) and is a function of impact energy and impact angle. As expected, the systems that are inclined 60° away from a stable configuration (C60 and E30) experience collisions with the highest impact velocities, which result in smaller values for M_lr. Although the resolution of our simulations is limited by the minimum mass of the fragment we define, we calculate the true M_lr for all collisions. The smallest mass of the largest remnant calculated in our simulations (although not included in the simulation due to lower limits on the fragment mass) is 2 × 10⁻⁶ M_⊕. While this is still orders of magnitude larger than the expected mass of 'Oumuamua, this is the largest remnant expected from a collision that will also produce a distribution of smaller fragments. The most massive body ejected was 0.91 M_⊕ in an E60 run, and so we expect a large distribution of fragment masses in misaligned CBDs.

Using our results, we can place an upper limit on the fraction of binaries in the Galaxy that are likely to eject the majority of their terrestrial planet-building material. Moe & Di Stefano (2017) compiled observations of early-type binaries to quantify the distributions of binary eccentricities. Close-in binaries with separations ≤0.4 au have small eccentricities ≤0.4 due to tidal circularization, while the eccentricities of more widely separated binaries are weighted to larger values. However, here we assume that binary eccentricity is uniformly distributed in the range 0.0–0.8 (Raghavan et al. 2010) and that the initial orientation of the CBD is uniformly distributed for simplicity.

A CBD ejects most of its solid material when the mutual inclination between colliding bodies becomes greater than 90°, which is ${v}_{{\rm{i}}}=2{v}_{{\rm{K}}}\sin {45}^{\circ }$, the PE = KE black line in Figure 2. The probability, p, that a binary ejects most of its solid material can be estimated with the fraction of orbits for which a particle in a CBD around the binary has a maximum inclination to either coplanar or polar that is greater than 45°.

Figure 3 shows precession paths of a particle angular momentum vector in the frame of the binary from our three-body simulations, for nine different binary eccentricities. To find p we find the fraction of the sphere's surface area that corresponds to paths with particle inclinations that are at some point greater than 45°, given by

$$\begin{eqnarray}&&p=1-\displaystyle \frac{2{A}_{{\rm{e}}1}+2{A}_{{\rm{e}}2}}{4\pi },\end{eqnarray} \tag{ 5 }$$

where A_e1 is the surface area of a bold green curved ellipse that corresponds to the phase space where the particle's inclination is always less than 45° to coplanar, and A_e2 is the surface area enclosed by the bold and dashed purple curved ellipse, which corresponds to the phase space where the particle's inclination is always less than 45° to polar. To calculate the curved area of the ellipses we consider the edge as the intersection of an elliptical cylinder with the sphere. The surface area of the curved ellipse is then given by

$$\begin{eqnarray}&&{A}_{{\rm{e}}}=2\pi {a}^{2}-4{a}^{2}\,{\sin }^{-1}\left(\displaystyle \frac{\sqrt{{a}^{2}+{b}^{2}}}{a}\right),\end{eqnarray} \tag{ 6 }$$

where a and b are the semimajor and semiminor axes, respectively, of the elliptical cylinder.

Zoom In Zoom Out Reset image size

In Figure 3 we list the p value for each binary eccentricity. Since we assume a uniform distribution of binary eccentricity, we take the average of all the p values to find the fraction of binaries in the Galaxy that are likely to eject most of their disk material. We find $\bar{p}=0.59$, meaning that with these assumptions, more than half of the binaries in the Galaxy are more likely to eject their terrestrial planet disk mass, and produce ISAs, than to form terrestrial planets. Additionally, we find that binaries with eccentricities greater than 0.4 have the same probability, p = 0.54, of ejecting most of their disk mass.

This is a strict upper limit to the fraction because we have assumed an initially isotropic orientation of the particle disk, which is equal to the orientation of the gas disk at the time of disk dispersal. However, the gas disk may evolve toward either coplanar or polar alignment for nonzero disk viscosity. Depending on the binary and disk parameters, the alignment timescale can be longer than the disk lifetime, and so the planetesimal disk may form in a misaligned disk. We note that the two currently observed polar circumbinary gas disks have external companions (Kennedy et al. 2019; Kenworthy et al. 2022) that truncate the outer part of the circumbinary disk, leading to a radially narrow disk and a short alignment timescale (Martin et al. 2022). The timescale for gas disk alignment also depends upon a number of other binary and disk parameters including the binary semimajor axis, binary eccentricity, disk viscosity, and disk temperature (e.g., Lubow & Martin 2018). Gas disk alignment does not proceed in a strictly linear fashion, and even a small initial misalignment can lead to very large misalignments during the disk evolution (Smallwood et al. 2019). Because of the complexity of this problem, we leave a more detailed investigation of this to future work.

We conducted a suite of n-body simulations around circular and highly eccentric, equal-mass binaries separated by 0.5 au. The circumbinary particle disks were inclined by either 30° or 60° to the binary orbital plane. We resolved all collisions with fragmentation, which allowed us to analyze the collisions and better understand the post-collision mass and velocity distributions.

We found that around highly eccentric binaries, CBDs with mild initial misalignment can result in large mutual inclinations between the particles. The impact velocity between two bodies in circular Keplerian orbits with mutual inclination i_m is given by ${v}_{{\rm{i}}}=2{v}_{{\rm{K}}}\sin ({i}_{{\rm{m}}}/2)$. CBDs that harbor particles with mutual inclinations greater than 90° have particle collisions with kinetic energies greater than the potential energy between the colliding system and the central binary, which results in ejection from the system, unless dynamics from the multiparticle disk inhibit this. This mechanism is an efficient source of ISAs with a wide range of sizes and velocities. These ISAs will have characteristics consistent with a violent past, such as the small size, elongated shape, and tumbling motion of 'Oumuamua. ISAs formed in this way will mostly be rocky in composition since the terrestrial planets are expected to form inside the snowline radius.

Assuming a uniform distribution of binary eccentricity and an isotropic distribution of the disk orientation relative to the binary, we find that 59% of binaries in the Galaxy are more likely to eject their CBD terrestrial material through high-velocity particle–particle collisions than to retain their material and build terrestrial planets. This is an upper limit to the fraction since a nonzero disk viscosity during the gas disk phases causes the alignment to tend toward either coplanar or polar depending upon the initial misalignment. These findings can help place constraints on occurrence rates for both ISAs and terrestrial CBPs.

We thank the anonymous referee for useful comments that improved the manuscript. We thank Stephen Lepp, Charlie Abod, and Ian Rabago for help with Figure 3. We acknowledge support from NASA through grant 80NSSC21K0395. Computer support was provided by UNLVs National Supercomputing Center.