Binary Planet Formation by Gas-assisted Encounters of Planetary Embryos

To avoid confusion, let us list several definitions that we use throughout the rest of the paper.

• 1.  
  Binary is shorthand for a binary planet, not to be mistaken with binary stars etc.
• 2.  
  Transient (also transient binary or transient pair) is a binary that forms by two-body encounters of planetary embryos (e.g., Astakhov et al. 2005), in our case with the assistance of the disk gravity as we shall demonstrate later. We choose the name transient because we find the typical lifetime of these binaries to be of the order of one stellarcentric orbital period.
• 3.  
  Hardening (e.g., Hills 1975) is a process during which the orbital energy of a binary configuration is dissipated and the separation of binary components decreases.
• 4.  
  Stability of a binary planet is discussed if it can survive more than one stellarcentric orbital period, which is usual after hardening. In principle, such a binary can be observed. We characterize stability by means of the lifetime over which the binary components remain gravitationally bound.
• 5.  
  Encounter refers to a close encounter of two or more planetary embryos (single or binary), when they enter one another's Hill sphere.
• 6.  
  Merger refers to a physical collision of two embryos. In our approximation, we replace the colliding embryos by a single object, assuming perfect merger (mass and momentum conservation).
• 7.  
  We denote orbital elements in the stellarcentric frame with a subscript "s" to distinguish them from the orbital elements of one binary component with respect to another (e.g., a_s is the stellarcentric semimajor axis but a is the semimajor axis of the binary configuration).

The individual constituents of our model are as follows. First, we consider radiation transfer, which is essential to properly reproduce the disk structure (Bitsch et al. 2013) and to account for all components of the Type I torque acting on low-mass planets (e.g., Baruteau & Masset 2008; Kley & Crida 2008; Kley et al. 2009; Lega et al. 2014).

Second, we use a two-fluid approximation to include a disk of pebbles, which serves as a material reservoir for the accreting embryos (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorný 2012).

Third, we also take into account that pebbles heat the accreting embryos, which in turn heat the gas in their vicinity. The migration is then modified due to the thermal torque (Benéz-Llambay et al. 2015; Masset 2017) and its dynamical component—the hot-trail effect (Chrenko et al. 2017; Eklund & Masset 2017; Masset & Velasco Romero 2017)—which perturbs the embryos in such a way that their orbital eccentricities are excited. This is due to the epicyclic motion of the embryo, which causes variations in the azimuthally uneven distribution of the heated (and thus underdense) gas.

The numerical modeling is done with the fargo_thorin code⁴ , which was introduced and described in detail in Chrenko et al. (2017). The code is based on fargo (Masset 2000). The model is 2D (vertically averaged) but planets are evolved in 3D. A number of important vertical phenomena were implemented, although some with unavoidable approximations.

A new phenomenon implemented in this study is the vertical hot-trail effect described by Eklund & Masset (2017), which can excite orbital inclinations. Such excitation should not occur for an isolated and non-inclined orbit because it is quenched by the growth in eccentricity, which is faster (Eklund & Masset 2017); however, the vertical hot-trail effect operates when a non-negligible inclination is initially excited by some other mechanism. For example, it can become important in a system of multiple embryos where close encounters temporarily pump up the inclinations. The vertical hot trail then starts to counteract the usual inclination damping by bending waves (Tanaka & Ward 2004).

Gas is treated as a viscous Eulerian fluid described by the surface density Σ, flow velocity ${\boldsymbol{v}}$ on the polar staggered mesh, and the internal energy . The pebble disk is represented by an inviscid and pressureless fluid with its own surface density Σ_p and velocity ${\boldsymbol{u}}$. We assume two-way coupling between the two fluids by linear drag terms, with the Stokes number τ calculated for the Epstein regime.

The RHD partial differential equations read

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+{\rm{\nabla }}\cdot ({\rm{\Sigma }}{\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+{\boldsymbol{v}}\cdot {\rm{\nabla }}{\boldsymbol{v}} & = & -\displaystyle \frac{1}{{\rm{\Sigma }}}{\rm{\nabla }}P+\displaystyle \frac{1}{{\rm{\Sigma }}}{\rm{\nabla }}\cdot \overleftrightarrow{T}-\displaystyle \frac{\displaystyle \int \rho {\rm{\nabla }}\phi \,{dz}}{{\rm{\Sigma }}}\\ & & +\ \displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{p}}}}{{\rm{\Sigma }}}\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{K}}}}{\tau }({\boldsymbol{u}}-{\boldsymbol{v}}),\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \epsilon }{\partial t}+{\rm{\nabla }}\cdot (\epsilon {\boldsymbol{v}}) & = & -P{\rm{\nabla }}\cdot {\boldsymbol{v}}+{Q}_{\mathrm{visc}}+{Q}_{\mathrm{acc}}\\ & & +\ \displaystyle \frac{2{\sigma }_{{\rm{B}}}}{{\tau }_{\mathrm{eff}}}({T}_{\mathrm{irr}}^{4}-{T}^{4})-2H{\rm{\nabla }}\cdot {\boldsymbol{F}},\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{p}}}}{\partial t}+{\rm{\nabla }}\cdot ({{\rm{\Sigma }}}_{{\rm{p}}}{\boldsymbol{u}})=-{\left(\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{p}}}}{\partial t}\right)}_{\mathrm{acc}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}=-\displaystyle \frac{\int {\rho }_{{\rm{p}}}{\rm{\nabla }}\phi \,{dz}}{{{\rm{\Sigma }}}_{{\rm{p}}}}-\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{K}}}}{\tau }({\boldsymbol{u}}-{\boldsymbol{v}}),\end{eqnarray} \tag{ 5 }$$

where the individual quantities are the pressure P, viscous stress tensor $\overleftrightarrow{T}$ (e.g., Masset 2002), gas volume density ρ, volume density of pebbles ρ_p, coordinate z perpendicular to the midplane, gravitational potential ϕ of the primary and the planets, Keplerian angular frequency Ω_K, viscous heating term Q_visc (Mihalas & Weibel Mihalas 1984), accretion heating term Q_acc related to the accretion sink term (∂Σ_p/∂t)_acc, Stefan–Boltzmann constant σ_B, effective vertical optical depth τ_eff (Hubeny 1990), irradiation temperature T_irr (Chiang & Goldreich 1997; Menou & Goodman 2004; Baillié & Charnoz 2014), midplane gas temperature T (the term ∝T⁴ describes vertical cooling), vertical pressure scale height H, and radiative flux ${\boldsymbol{F}}$.

The ideal-gas state equation is used as the thermodynamic closing relation:

$$\begin{eqnarray}&&P=(\gamma -1)\epsilon ={\rm{\Sigma }}\displaystyle \frac{R}{\mu }T,\end{eqnarray} \tag{ 6 }$$

where γ is the adiabatic index, R is the universal gas constant, and μ is the mean molecular weight. Equation (6) has been widely used in numerical models to relax inferior isothermal approximations and to account for the disk thermodynamics (through the energy equation), which is important for accurate migration rates (e.g., Kley & Crida 2008; Kley et al. 2009; Lega et al. 2014). We also point out that the given state equation neglects the radiation pressure and phase transitions, which is a valid assumption in low-temperature disks (D'Angelo et al. 2003).

Finally, flux-limited diffusion and the one-temperature approximation are utilized to describe the in-plane radiation transport, leading to

$$\begin{eqnarray}&&{\boldsymbol{F}}=-D{\rm{\nabla }}T=-{\lambda }_{\mathrm{lim}}\displaystyle \frac{16{\sigma }_{{\rm{B}}}{T}^{3}}{{\rho }_{0}\kappa }{\rm{\nabla }}T,\end{eqnarray} \tag{ 7 }$$

where D denotes the diffusion coefficient, λ_lim is the flux limiter according to Kley (1989), ρ₀ is the midplane volume density, and κ(ρ, T) is the material opacity. We use the opacity by Bell & Lin (1994) for both the Rosseland and Planck opacities.

For completeness, we provide the accretion heating formula

$$\begin{eqnarray}&&{Q}_{\mathrm{acc}}=\displaystyle \sum _{i}\displaystyle \frac{{{GM}}_{\mathrm{em},i}{\dot{M}}_{\mathrm{em},i}}{{R}_{i}{S}_{\mathrm{cell}}},\end{eqnarray} \tag{ 8 }$$

where the sum goes over all embryos with indices i, masses M_em,i, self-consistently calculated pebble accretion rates ${\dot{M}}_{\mathrm{em},i}$, and physical radii R_i. G is the gravitational constant and S_cell is the surface area of the cell that contains the respective embryo and in which the heat is liberated (Q_acc is zero in other cells).

Planets are evolved on 3D orbits using the ias15 integrator (Rein & Liu 2012; Rein & Spiegel 2015). Planetary collisions are treated as perfect mergers. The planet–disk interactions are calculated by means of the vertical averaging procedure of Müller et al. (2012). The planetary potential is adopted from Klahr & Kley (2006) and has the smoothing length r_sm = 0.5R_H, where R_H is the planet's Hill sphere.

When computing the torque acting on an embryo, we do not exclude any part of the Hill sphere because we focus on low-mass embryos (Lega et al. 2014). Such an exclusion is required only when embryos exceed masses ≃10 M_⊕ and form a circumplanetary disk. This disk should not contribute to the gas-driven torque because it comoves with the embryo (Crida et al. 2008). Our model ignores the torques from pebbles (Benéz-Llambay & Pessah 2018) because we assume relatively low pebble-to-gas mass ratios (less than 0.001). But we point out that during accretion of pebbles, we account for the transfer of their mass and linear momentum onto the embryo.

An important ingredient when investigating 3D planetary orbits is the inclination damping (e.g., Cresswell et al. 2007). We include the damping by using the formula from Tanaka & Ward (2004). In our case, the damping acceleration perpendicular to the disk plane reads

$$\begin{eqnarray}{a}_{z}=\left\{\begin{array}{ll}\beta \displaystyle \frac{{\rm{\Sigma }}{{\rm{\Omega }}}_{{\rm{K}}}}{{c}_{{\rm{s}}}^{4}}(2{A}_{z}^{c}{v}_{z}+{A}_{z}^{s}z{{\rm{\Omega }}}_{{\rm{K}}}), & I\gt {I}_{0}\\ 0, & I\leqslant {I}_{0}\end{array}\right.\end{eqnarray} \tag{ 9 }$$

where β = 0.3 (e.g., Pierens et al. 2013), ${A}_{z}^{c}=-1.088$, and ${A}_{z}^{s}=-0.871$ are fixed coefficients, c_s is the sound speed, v_z is the embryo's vertical velocity, and z is its vertical separation from the midplane. Σ and Ω_K are evaluated along the embryo's orbit.

In writing Equation (9), we introduce a simple modification of Tanaka & Ward's formula. We assume that the inclination damping does not operate when the orbital inclination I is below a certain critical value I₀. The motivation for this modification stems from the findings of Eklund & Masset (2017), who investigated the orbital evolution of a hot (accreting) planet in a 3D radiative disk. Not only did they find the excitation of the eccentricity due to the hot-trail effect, but they also showed that the effect has a vertical component that can excite the inclinations.

In our 2D model with accretion heating, the excitation of eccentricity is reproduced naturally. Regarding the inclinations, we simply assume that the hot-trail effect operates vertically as well and balances the inclination damping up to the value I₀. For larger inclinations, the damping takes over and the standard formula for a_z applies. We consider I₀ to be a free parameter of the model and choose I₀ = 10⁻³ rad ≃ 0057, which is at the lower end of the results of Eklund & Masset (2017).

The protoplanetary disk model used in all our RHD simulations is exactly the same as in Chrenko et al. (2017), including the initial and boundary conditions (de Val-Borro et al. 2006). The parameters⁵ characterizing the initial gas disk are the surface density Σ = 750(r/(1 au))^−0.5 g cm⁻², kinematic viscosity ν = 5 × 10¹⁴ cm² s⁻¹, adiabatic index γ = 1.4, mean molecular weight μ = 2.4 g mol⁻¹, vertical opacity drop c_κ = 0.6, and disk albedo A = 0.5. The central star has an effective temperature T_⋆ = 4370 K, stellar radius R_⋆ = 1.5 R_⊙, and mass M_⋆ = 1 M_⊙. The domain stretches from r_min = 2.8 au to r_max = 14 au in radius and spans the entire azimuth, having a grid resolution of 1024 × 1536 (rings × sectors).

The gas disk is numerically evolved to its thermal equilibrium and only after that is the coupled pebble disk introduced. Pebbles are parameterized by the radial pebble mass flux ${\dot{M}}_{{\rm{F}}}=2\times {10}^{-4}\,{M}_{\oplus }\,{\mathrm{yr}}^{-1}$ (Lambrechts & Johansen 2014), the Schmidt number Sc = 1, coagulation efficiency _p = 0.5 (Lambrechts & Johansen 2014), bulk density ρ_b = 1 g cm⁻³, and turbulent stirring efficiency α_p = 1 × 10⁻⁴ (e.g., Youdin & Lithwick 2007). To infer the pebble sizes, we assume the drift-limited growth regime (Birnstiel et al. 2012; Lambrechts & Johansen 2014), leading to pebble sizes of several centimeters.

For reference, Figure 1 shows the radial profiles of the aspect ratio h(r) = H/r and temperature T(r) of the equilibrium gas disk. The h(r) profile has a maximum near 4 au, where the opacity peaks just before the sublimation of water ice (Bell & Lin 1994). Therefore the heat produced by viscous heating is not easily radiated away from this opaque region and the disk puffs up.

Zoom In Zoom Out Reset image size

Near 7 au, the disk starts to flare because the vertical optical depth is small enough for the incoming radiation to penetrate deeper and heat the disk. The transition to the flared outer parts produces a zone of convergent Type I migration for planetary embryos (Bitsch et al. 2013; Pierens 2015).

Let us now discuss and analyze our nominal simulation in which binary planets were found to form. Initially, we set four embryos on circular orbits with stellarcentric semimajor axes a_s = 5, 6.7, 8.4, and 10.1 au, inclinations I_s = I₀ = 10⁻³ rad ≃ 0057, and randomized longitudes. The initial mass of each embryo is M_em = 3 M_⊕ and their orbital separations are equal to 16 mutual Hill radii:

$$\begin{eqnarray}&&{R}_{\mathrm{mH},{ij}}=\displaystyle \frac{{a}_{{\rm{s}},i}+{a}_{{\rm{s}},j}}{2}{\left(\displaystyle \frac{{M}_{\mathrm{em},i}+{M}_{\mathrm{em},j}}{3{M}_{\star }}\right)}^{1/3}.\end{eqnarray} \tag{ 10 }$$

Embryos are numbered 1, 2, 3, and 4 from the innermost outward.

Figure 2 shows the evolution of embryos over 140 kyr of the full RHD simulation with pebble accretion and respective heating. At first, the embryos undergo convergent migration toward their zero-torque radius. Without the heating torques, the embryos would concentrate near 7 au thanks to the contribution of the entropy-related corotation torque (Paardekooper & Mellema 2008), which is positive from ≃4 to 7 au in this particular disk model.

Zoom In Zoom Out Reset image size

With the heating torques, however, the zero-torque radius is shifted further out because these torques are always positive (Benéz-Llambay et al. 2015). Moreover, the hot-trail effect quickly excites orbital eccentricities. Within ≃2 kyr, the eccentricities reach e_s ≃ 0.02 for the innermost embryo and e_s ≃ 0.04 for the outermost.

The inclinations first remain constant near the prescribed I₀ value, with only small temporal excitations not exceeding I_s ≃ 02. Even these initially small inclinations are enough to modify the encounter geometries in such a way that the system becomes gradually stirred in the vertical direction. Once the system becomes closely packed, at about ≃35 kyr into the simulation, the mutual close encounters pump the inclinations significantly, typically to I_s ≃ 1° and even up to I_s ≃ 2°. Planets can pass above or below each other because their vertical excursions are comparable to (or larger than) their Hill spheres. For example, the maximal vertical excursion of a 5 M_⊕ embryo at 10 au is z_max ≃ 2R_H when I_s = 2°.

Due to excited eccentricities, the embryos never form a stable resonant chain, in accordance with Chrenko et al. (2017). Moreover, the excited inclinations help the embryos to avoid collisions and mergers for a long period of time. Consequently, close encounters of embryos are frequent in the system. Between ≃98 kyr and ≃114 kyr, strong unphysical oscillations of the stellarcentric orbital elements appear in Figure 2 for some of the embryos, indicating formation of gravitationally bound binary planets.⁶

To identify the events related to the binary formation in the simulation described above, we computed orbits of the relative motion among all possible pairs of embryos and selected those with a < R_mH, e < 1. The results are shown in Figure 3. Throughout the simulation, we found 65 time intervals during which at least two embryos are captured on a mutual elliptical orbit, changing their relative orbital energy from initially positive to negative (and back to positive when the capture terminates). Subsequently, we also scanned the sample of bound pairs and looked for cases when a third embryo has its distance from a pair d < R_mH to identify three-body (pair–embryo) encounters. These are highlighted in Figure 3 by arrows.

Zoom In Zoom Out Reset image size

Analyzing the lifetime of the bound pairs, we found that most of them dissolve before finishing one stellarcentric orbit. However, we found a single case when several binary configurations existed consecutively for a prolonged period of time between ≃98 kyr and ≃114 kyr. This time interval is bordered by three-body encounters, which usually cause the binary separation to drop.

In summary, bound pairs can form in two-body encounters but they quickly dissolve unless they undergo a three-body encounter with one of the remaining embryos. The latter process is known as binary hardening (e.g., Hills 1975, 1990; Goldreich et al. 2002; Astakhov et al. 2005) and it occurs when an external perturber removes energy from a binary system, which then becomes more tightly bound. To distinguish between two types of events contributing to the formation of bound pairs, we call those formed in two-body encounters transient binaries because of their typically short dynamical lifetime.

Transient binaries have been the subject of many different studies, for example, in the three-body Hill problem (e.g., Simó & Stuchi 2000; Astakhov et al. 2005). Their formation is possible because the orbital energy of two bodies in no longer conserved when additional perturbers (e.g., the central star) are present (e.g., Cordeiro et al. 1999; Araujo et al. 2008). Our system is of course more complicated because additional gravitational perturbations arise from the gas disk. We will demonstrate in Section 3.4 that the gas indeed facilitates formation of transients.

One last question we address here is whether or not the occurrence of bound pairs in Figure 3 is related to the vertically stirred orbits of embryos. To find an answer, we looked for bound pairs (with a < R_mH, e < 1) in one of our previous simulations reported in Chrenko et al. (2017) (dubbed Case III), where the inclinations were damped in the standard way (Tanaka & Ward 2004). We found only six bound pairs (all transients) compared to 65 cases in our nominal simulation presented here. Therefore, the excited inclinations importantly change the outcomes of embryo encounters in the gas disk and help to form transients.

Here we investigate the formation of a transient pair of embryos 2 and 4, which precedes the binary hardening events in our nominal simulation. The pair forms in a two-body interaction at t ≃ 97.76 kyr. To see whether the embryo–disk interactions assist in the process, we show in Figure 4 the evolution of the perturbed gas surface density Σ during the encounter.

Before the encounter (panel (a)), the usual structures can be seen in the disk. The hot (underdense) trail of the outer embryo 2 can be seen as a dark oval spot in the bottom right quadrant (it moves down and gets larger in time). The hot trail of the inner embryo 4 is less prominent and looks like an underdense gap attached to the embryo from inside.

As the embryos approach in panel (b), the outer spiral arm of the inner embryo 4 and the inner spiral arm of the outer embryo 2 join and overlap. The overlap forms a strong density wave positioned between the embryos. The overdensity increases as part of the wave becomes trapped between the Hill spheres of the two embryos in panel (c). From panel (c) to (d), the embryos cross this shared density wave.

In panel (d), the embryos are so close to each other that they effectively act on the disk as a single mass, and the previously shared spiral arm splits into an inner and an outer component with a small pitch angle. There are two more spirals with a larger pitch angle that are leftovers of the initial wakes launched by the embryos. In panel (e), all spirals blend into a single pair of arms. The embryos enter one another's Hill sphere between panels (d) and (e) through the vicinity of the Lagrange points L₁ for the outer and L₂ for the inner embryo (Astakhov et al. 2005). The embryos are captured on a prograde binary orbit (in panel (e), embryo 4 orbits the central embryo 2 counterclockwise).

The spiral arm crossing that appears during the encounter is known to produce strong damping effects on the embryos (Papaloizou & Larwood 2000). It is thus likely that the gas supports the gravitational capture by dissipating the orbital energy. To quantify this effect, we measured the total gravitational force ${{\boldsymbol{F}}}_{{\rm{g}},i}$ exerted by the disk on each embryo and we calculated the mechanical power

$$\begin{eqnarray}&&{P}_{i}={\int }_{\mathrm{disk}}{{\boldsymbol{v}}}_{i}\cdot {{\boldsymbol{F}}}_{{\rm{g}},i}{dS},\end{eqnarray} \tag{ 11 }$$

where ${{\boldsymbol{v}}}_{i}$ is the velocity vector of the embryo and the integral goes over the entire disk. P_i directly determines the rate of change of the orbital energy ${\dot{E}}_{i}$ of each embryo (e.g., Cresswell et al. 2007).

Panel (f) of Figure 4 shows the total balance of the energy subtraction (or addition) for embryos 2 and 4, P = P₂ + P₄. The energy is transferred to the disk when P < 0 and subtracted from the disk when P > 0. It is obvious that throughout the closest approach (panels (b)–(e)), the orbital energy of the embryos dissipates and thus the influence of the gas disk on the formation of transients is confirmed.

The transient pair of embryos 2 and 4 does not dissolve; instead it is further stabilized in three-body encounters with the remaining embryos. Figure 5 shows these encounters in detail. First, the transient pair encounters embryo 1 at t ≃ 97.85 kyr. During this encounter, one component of the binary (embryo 4) becomes unbound and is deflected away, but the incoming embryo 1 takes its place in an exchange reaction so the binary continues to exist.

Similar situations occur at about ≃98.69 kyr, when the configuration 1 + 2 changes to 1 + 3, and at ≃103.4 kyr, when the configuration 1 + 3 changes to 3 + 4. Figure 5 also reveals that the binary becomes hardened with each encounter since the overlap between the r(t) curves for the binary components becomes tighter.

An even clearer indication of binary hardening is provided by Figure 6, where we plot the temporal evolution of the orbital elements of the binary planet. The color of the curves changes each time there is a change in the composition of the binary. One can further notice that the exchange interactions produce sudden decreases of a and also of e. Between the exchange interactions, a decreases smoothly whereas e generally increases in an oscillatory manner. We will describe these variations later.

Zoom In Zoom Out Reset image size

Considering the binary planet inclination, the transient binary 2 + 4 forms with a prograde orbit and a relatively low inclination. During the first three-body encounter, the binary is reconfigured to a retrograde orbit with the inclination oscillating between 100° and 170°. The inclination then slowly evolves toward 180° regardless of the swap encounters, which only diminish the amplitude of the oscillation.

The binary planet does not survive to the end of our simulation. At ≃111.26 kyr, it undergoes a three-body exchange interaction with embryo 1 during which their Hill spheres overlap for a prolonged period of time (see the spike in Figure 3 at ≃111.26 kyr). Consequently, the binary inclination is flipped from the retrograde configuration to I = 80°. In this configuration, the binary undergoes a fast decrease of a, accompanied by an equally fast increase of e. Consequently, the binary planet ends its life in a merger into a single body.

The lifetime of the hardened binary in our nominal simulation is long enough (∼10⁴ yr) to be interesting. There are two basic questions that we shall now address. First, what is the evolution of such a binary if the surrounding embryos and their perturbations are ignored? And second, what causes the changes in the binary orbital elements between the three-body encounters in Figure 6?

To answer these questions, we discard the non-binary embryos and restart the simulation from the configuration of the binary planet,⁷ gas, and pebbles corresponding to t = 100 kyr. Three models are numerically evolved for 45 kyr. The first one has the same setup as the initial simulation (apart from the ignored non-binary embryos). In the second one, the accretion heating is disabled but the mass of the binary components can still grow by pebble accretion. In the third one, we again switch off accretion heating and discard the pebble disk; the binary mass therefore remains constant.

The orbital evolution of the binary in these three cases is shown in Figure 7. The evolution of the inclination is more or less the same, regardless of the model, and converges toward a fully retrograde configuration. The semimajor axis decreases as a consequence of pebble accretion, which transports the linear momentum and mass onto the binary components, thus changing their orbital angular momentum. It is worth noting that if the pebble accretion and accretion heating are ignored, the isolated binary planet evolving in the radiative disk exhibits only minor orbital changes (once it adjusts to the removal of the surrounding embryos at the beginning of the restart).

Zoom In Zoom Out Reset image size

The eccentricity substantially changes only in the model with accretion heating, otherwise it oscillates around its initial value or exhibits a slow secular variation. In other words, not only is the hot-trail effect important for exciting the eccentricities of individual embryos before the encounter phase, but it is also responsible for pumping the eccentricity of the binary up to an asymptotic value e ≃ 0.75.

These findings justify our incorporation of pebble accretion and accretion heating into the model because both phenomena affect the rate of change of the binary orbital elements. Pebble accretion diminishes the semimajor axis and accretion heating excites the eccentricity.

When protoplanetary disks undergo dispersal due to photoevaporation, the emerging planetary systems may become unstable (e.g., Lega et al. 2013). Here we test whether the hardened binary planet could survive the gas removal phase and the subsequent orbital instabilities. Since the binary undergoes three-body encounters during the disk phase, they can also be expected after removal of the disk.

To investigate the evolution after the photoevaporation, we remove the fluid part of the model (i.e., gas and pebbles) instantly and continue with a pure N-body simulation. The orbits are integrated for an additional 10 Myr. To account for the chaotic nature of an N-body system with close encounters, we extract 48 orbital configurations of the embryos from between ≃99.7 kyr and ≃102.4 kyr of the nominal simulation and use them as the initial conditions for 48 independent integrations.

Our aim is to quantify the survival rate of binary planets. Figure 8 shows the evolution of the fraction N_surv/N_tot, where N_surv is the number of N-body systems still containing a binary planet at simulation time t and N_tot is the total number of systems (48). The dependence exhibits a steep decrease—the binary dissolves before 0.1 Myr of evolution in ≃45% of cases and before 1 Myr in ≃76% of cases. However, the trend for t ≥ 1 Myr becomes rather flat. The binary planet survives the whole integration timespan in 15% of our runs.

Zoom In Zoom Out Reset image size

We estimate the fraction of planetary systems f_surv in which the binary can survive the orbital instabilities. An estimate for young systems can be readily made by taking the final fraction of our integrations conserving the binary, yielding f_{surv,10 Myr} ≃ 15%. To make an estimate for older systems, we performed a power-law extrapolation of the flat tail of the distribution in Figure 8. The resulting extrapolation 0.22(t/(1 Myr))^−0.21 yields f_{surv,4.5 Gyr} ≃ 4% for t = 4.5 Gyr (i.e., comparable to the age of the solar system).

Let us briefly discuss the orbital architecture of the individual systems at the end of our integrations (Figure 9). The systems that retained the binary can be divided into two classes. The more common first class comprises five systems (simulation numbers 5, 12, 28, 34, 37) in which one of the binary components undergoes an early collision (at about ≃0.1 Myr) with one of the remaining embryos while maintaining the binary configuration. The collision reduces the multiplicity of the system and changes the mass ratio of the binary components from 1/1 to 2/1. The stability of such systems is obvious from Figure 9 because eccentricities are only marginally excited and the planets do not undergo orbital crossings.

Zoom In Zoom Out Reset image size

The less frequent second class of orbital architectures (simulation numbers 17 and 38) includes systems in which no collision occurred yet the binary managed to survive close encounters. Orbits of the single embryos in these simulations are moderately eccentric and have orbital crossings with the binary. It is very likely that these systems would reconfigure and the binary would dissolve when integrating for t > 10 Myr. However, we believe that these two cases are accounted for by the power-law extrapolation when estimating f_surv for old systems.

So far the analysis has been based on our nominal simulation. Although a broader study of the parametric space is difficult using RHD simulations, it is important to quantify how common it is for transients to become hardened and stabilized in three-body encounters. Also, it is desirable to test whether the result is sensitive to the choice of initial separations, embryo masses, and embryo multiplicity.

We perform four additional full RHD simulations in which we vary the initial conditions for embryos (the disk remains the same as in Section 3.1). Embryos start with different random longitudes, and inclinations are I_s = I₀. Two simulations (denoted I and II) are run with four embryos, each having M_em = 3 M_⊕ again, but the innermost embryo is initially placed at a_s = 6 au and the remaining ones are spaced by 10 mutual Hill radii. Two simulations (denoted III and IV) include eight embryos with M_em = 1.7 M_⊕, the innermost one being placed at a_s = 5 au and the others having initial separations of eight mutual Hill radii.

Simulations I and II cover 120 kyr of evolution. A common feature of these runs is a merger occurring relatively early (at ≃40 kyr and ≃51 kyr in simulations I and II, respectively), followed by a second merger (at ≃81 kyr and ≃94 kyr, respectively). Simulation III covers 180 kyr of evolution. There is a late violent sequence at ≃150 kyr during which two mergers occur and three embryos are scattered out of the simulated part of the disk.

In simulations I–III, only transient binaries are formed. However, we detect one case of a hardened binary in simulation IV, which covers 120 kyr of evolution. The simulation also contains 190 transients compared to 65 cases found in our nominal simulation, which implies that the increased multiplicity of the system (eight embryos instead of four) logically increases the frequency of embryo encounters.

Figure 10 shows the evolution of orbital elements of the binary configurations participating in binary hardening. First, a transient consisting of embryos 2 and 7 forms at about 122.5 kyr. After ≃1 kyr, it undergoes a three-body encounter during which the configuration changes to embryos 2 and 6 and the semimajor axis of the binary decreases. Then the binary evolves for about 5 kyr due to pebble accretion and additional three-body encounters. The secular rate of change of a is clearly related to the value of I, suggesting that the deposition of pebbles onto a prograde binary causes the separation to decrease faster than in a retrograde case. As in the nominal simulation, the binary separation shrinks until the binary merges.

Zoom In Zoom Out Reset image size

Although the statistics of our five (one nominal and four additional) RHD simulations is still poor, we can nevertheless conclude that a binary planet with a considerable lifetime can form in at least some cases. We can also make a crude estimate of the formation efficiency f_form, which we define as the fraction of simulations in which a binary planet was formed and then hardened, obtaining f_form ≃ 0.4.

For simulations that formed such binaries, we also define the total time τ_bp for which the binary planet existed in the system (regardless of which embryos were bound). The binary in our nominal simulation existed for ≃16 kyr and the binary in simulation IV existed for ≃6 kyr in total. Taking the arithmetic mean of these two values, we obtain τ_bp ≃ 11 kyr.

The choice of the RHD model was not a priori motivated by studying the formation of binary planets. The initial motivation was to check how the evolution described in Chrenko et al. (2017) changes if the orbits of embryos become inclined due to the vertical hot-trail effect (only the horizontal hot trail was modeled in Chrenko et al. 2017). We found binary planet formation as an unexpected yet natural outcome of the model.

It is possible that a less complex model could be applied to scan the parametric space, e.g., by N-body integrations. But although there are many state-of-the-art N-body models of migration in multiplanet systems with disks (e.g., Cossou et al. 2013, 2014; Izidoro et al. 2017), none of them (to our knowledge) identified the formation of binary planets. This suggests that there might be issues preventing binary planet formation in such models.

We demonstrated in Section 3.4 that hydrodynamic effects are important during the formation of transient binaries by two-body encounters. The pair of approaching embryos creates perturbations in the disk that differ from perturbations arising from isolated embryos. The prescriptions for the disk torques that are currently used in N-body models (Cresswell & Nelson 2008; Paardekooper et al. 2010, 2011; Fendyke & Nelson 2014) cannot account for such effects because they were derived from models containing a single embryo, not an interacting pair. Moreover, such prescriptions do not account for the thermal torque and hot-trail effect.

We point out that our RHD model has some limitations as well. First, although the grid resolution leads to numerical convergence of the migration rate for low-mass planets (see, e.g., Lega et al. 2014; Brož et al. 2018, for discussions of resolution), it is not tuned to resolve binary configurations with the smallest orbital separations well enough. However, looking at the turquoise curve in Figure 7 (i.e., the case without additional perturbers), the marginal change of orbital elements indicates a very low level of numerical dissipation, safely negligible over a typical binary lifetime during the disk phase. Second, the vertical hot trail cannot be implemented in our 2D model in a self-consistent way. 3D simulations would be required to assess the importance of the vertical dimension, which we neglect here. Finally, the magnitude of the thermal torques depends on the thermal diffusivity and therefore on the opacity (Masset 2017). However, we used a single opacity prescription and it remains unclear whether the described effects work the same way, e.g., in a low-density disk undergoing photoevaporation.

We found that two-body encounters of planetary embryos in the gas disk can establish transient binaries. A binary planet with a considerable lifetime can form from a transient by a three-body encounter that provides the necessary energy dissipation to make the pair more tightly bound. Here we discuss whether there are other mechanisms suitable for formation of binary planets.

An additional possibility exists during the reaccumulation phase after a large impact of two embryos approaching on initially unbound trajectories, as discussed by Ryan et al. (2014). But this situation is highly unlikely. Head-on collisions usually disrupt the protoplanets in such a way that the reaccumulation forms a large primary and a low-mass disk, from which a satellite can be assembled but not a binary companion. Only a special grazing geometry with a large impact parameter can be successful (Ryan et al. 2014) and it can only produce binaries with separations of a few planetary radii due to the angular momentum deficit of such an encounter.

Finally, planets can be captured in a binary configuration by means of tidal dissipation. Ochiai et al. (2014) studied the evolution of three hot Jupiters around a host star and discovered that binary gas giants can form in ∼10% of systems that undergo orbital crossings.

In our hydrodynamic simulations, the components of binaries have comparable masses of several M_⊕ before they merge. But as we found in the follow-up gas-free N-body simulations, the system often stabilizes by a collision of one of the embryos with the binary. If the binary survives the collision, the mass ratio of the components increases to 2/1. Therefore it seems that, if born from a population of equal-mass embryos (as obtained in the oligarchic growth scenarios), binary planets would preferentially exist with the component mass ratio 1/1 or 2/1. This aspect of our model is related to the choice of the initial embryo masses and is of course too simplified to capture the outcome of models where accretion creates a range of embryo masses.

Although simulations with gas accretion onto the planets are beyond the scope of our paper, we believe that runaway accretion of gas onto the binary would disrupt it. We thus expect that the binary planets formed in three-body encounters cannot exceed the masses of giant planet cores. This could be ensured by the mechanism of pebble isolation (Lambrechts et al. 2014; Bitsch et al. 2018) or simply because binaries could form late, just before dispersal of the gas disk. However, it is possible that binary giant planets form later by the mechanism of tidal capture (Ochiai et al. 2014).

Orbital evolution due to tidal dissipation is without doubt an important factor for the stability of binary planets. However, it is difficult to assess the tidal effects at this stage because there are many unknown parameters. A large uncertainty lies in the parameter k₂/Q, where k₂ is the degree 2 Love number and Q is the tidal quality factor (e.g., Harris & Ward 1982). These parameters reflect the interior structure and thus they depend on the planetary composition (water-rich versus silicate-rich) and state (cold versus magma worlds), the latter of which changes on an uncertain timescale. Moreover, our model treats the planets as point-mass objects, and therefore we have no information about their rotation, which is important to determine the level of spin–orbital synchronicity. Additionally, similar masses of the binary components and their (possibly) retrograde and highly eccentric orbit make the analysis of tides even more complicated. For all these reasons, the model of tidal evolution of a binary planet should be sufficiently complex and should account for the internal structure and rheology (e.g., Boué et al. 2016; Walterová & Běhounková 2017).

We define the occurrence rate of binary planets f_bp as the fraction of planetary systems that are expected to contain at least one binary planet hardened by three-body encounters. An order-of-magnitude estimate can be obtained by dividing the timescale τ_bp ≃ 10⁴ yr for which these binary planets are typically present in our simulations (Section 3.8) by the lifetime of protoplanetary disks τ_disk ≃ 10⁷ yr (e.g., Fedele et al. 2010), correcting for the formation efficiency f_form ≃ 0.4 (Section 3.8) and for the fraction of the emerging planetary systems f_stab in which the binary planet can survive after dispersal of the gas disk (Section 3.7). This leads to

$$\begin{eqnarray}&&{f}_{\mathrm{bp}}={f}_{\mathrm{form}}\displaystyle \frac{{\tau }_{\mathrm{bp}}}{{\tau }_{\mathrm{disk}}}{f}_{\mathrm{stab}}.\end{eqnarray} \tag{ 12 }$$

We quantify two characteristic values for young planetary systems (using f_{stab,10 Myr} ≃ 0.15) and for old planetary systems (using f_{stab,4.5 Gyr} ≃ 0.04), leading to f_{bp,10 Myr} ≃ 6 × 10⁻⁵ and f_{bp,4.5 Gyr} ≃ 2 × 10⁻⁵. In other words, one out of ≃(2–5) × 10⁴ planetary systems should contain a binary planet formed by two- and three-body encounters.

We emphasize that the estimate is highly uncertain because it is inferred from a small number of tests. Moreover, our simulations cover only a small region of the protoplanetary disk, the evolution timescale is still short compared to the disk's lifetime, and the number of embryos is relatively low. Last but not least, our estimate assumes that binaries only form by the processes identified in this work.

As pointed out in Section 1, the estimated occurrence rate should motivate a systematic search for binary planets in the observational data. For example, binary planets should be detectable by the Hunt for Exomoons with Kepler (HEK) project (Kipping et al. 2012). The sensitivity of this survey is ≃40% for binaries with Pluto–Charon mass ratios (Kipping et al. 2015).

The possible existence of binary planets opens new avenues in planetary sciences. A study of long-term orbital dynamics and stability of binaries in various systems is needed (e.g., Donnison 2010), including an assessment of their tidal evolution. Binary planets are challenging for hydrodynamic modeling as well. Local high-resolution simulations of interactions with the disk, pebble accretion, pebble isolation (e.g., Bitsch et al. 2018), and gas accretion (e.g., Lambrechts & Lega 2017) should be performed for binary planets (preferably in 3D) to understand the impact of these processes in detail.

In this paper, we reported various fates of binary planets (considering the set of gas-free N-body simulations). Frequently, one of the binary components underwent a collision with an equally large impactor, or the binary merged. The merger of the binary planet is a process that should be investigated (e.g., by the smoothed-particle hydrodynamics (SPH) method, Jutzi 2015).

A merger can occur in a situation when the binary orbit is inclined with respect to the global orbital plane, and the resulting body would then retain the initial angular momentum of the binary, forming a planet (or a giant planet core) with an angular tilt of the rotational axis. It might be worth investigating the relation of such an event to the origin of Uranus.

Moreover, the collision of the binary components would statistically occur at high impact angles. It is interesting that such impact angles and similar masses of the target and the impactor (although on unbound trajectories) were also used for a successful explanation of the impact origin of the Earth–Moon system (Canup 2012).