© The Authors 2023

1 Introduction

Roughly twenty percent of the cold classicals in the Kuiper belt are binary systems (Noll et al. 2008), with the mass ratio of the components being close to unity. A well-known example is the Pluto-Charon system (Christy & Harrington 1978) located at an orbital distance from the Sun of around 39 au. With Pluto being roughly eight times as massive as Charon, it is often categorised as a nearly equal-mass binary. Recent studies have shown, with increasing certainty, that equal-mass binaries are not rare in the Kuiper Belt and their occurrence rate is estimated to be of at least several percent (Noll et al. 2008)¹.

In addition, the asteroid belt features numerous objects with companions (Merline et al. 2002). Some of these binaries have similar masses (e.g. 90 Antiope, 2006 VW139, 2017 YE5, and 69230 Hermes Marchis et al. 2008, 2004, 2012; Agarwal et al. 2017), while others have large mass ratios, so that calling the smaller objects asteroid moons or satellites and the larger ones planet or primary would seem more correct.

The Earth-Moon system is a binary with a relatively large mass ratio as well, although here the definition of equal mass binary reveals its ambiguity. The most recent definition of ‘equal-mass’ and ‘near-equal-mass’ reveals the need for the barycenter to lie outside both companions (Stern & Levison 2002). Another way of categorising the mass equality of binaries is the ‘tug-of-war’ ratio, expressing the planetary gravity over the solar gravity for a satellite. Equal-mass planet binaries are often also called ‘double planets’, making the definition of an equal-mass binary even more vague; the Moon ticks off all the criteria for being a planet on its own and the tug-of-war ratio for the Earth-Moon system is 0.46 (Herschel 1833), but the barycenter lies within the body of the Earth making the Earth-Moon system: a double-planet system, an equal-mass binary, or a primary (Earth) and secondary (Moon) system, respectively (Russell 2017).

The formation of binary planets is still an open question and a number of channels have been proposed. In the classical model it is thought, at least in the case of Pluto-Charon (Canup 2005, Canup 2011) and Earth-Moon systems, that they formed via giant impacts on the primaries, smashing debris into stable orbiting secondaries (Cameron & Ward 1976; Canup 2012). This scenario requires specific impact parameters that ensure sufficient angular momentum in the ejected debris to allow for the formation of an object outside the Roche limit.

A second viable mechanism to create binary planet(esimals) is capture via three-body interactions in an early, dense phase of the solar system, where such interactions were much more frequent than they are today (Goldreich et al. 2002; Astakhov et al. 2005; Lee et al. 2007; Schlichting & Sari 2008). This might be the case for the dynamically hot classicals in the Kuiper Belt, which are more scattered due to perturbations by Neptune and have relatively wide mutual separations. They also show a broad and more chaotic distribution in orbital parameters, are subject to Kozai-Lidov effects (Kozai 1962; Lidov 1962), and show a much lower equal mass binary occurrence than the cold classical population because disruption is relatively easy (Morbidelli & Nesvorny 2020). The objects in this population thus have a dynamical profile that fits a random capture process which suggests a different formation process than the cold classicals. The cold classicals have strong color correlations, indicating the same composition (Benecchi et al. 2009), they have mass ratios close to unity, and show a clear preference for prograde orbits (80%) versus retrograde orbits (20%; Grundy et al. 2019), strongly disfavouring a capture event.

Finally, a recently proposed mechanism to explain the formaton of planetary binaries is the streaming instability (SI; Youdin & Goodman 2005; Johansen & Klahr 2011; Visser et al. 2021). In this scenario, dust clumps together via an interaction with the disk gas and when it becomes gravitationally bound, it collapses on a dynamical timescale to form comets, planetesimals, and protoplanets (Nesvorný et al. 2010, 2019). Depending on the amount of angular momentum present in the clump before the collapse (Nesvorný et al. 2021), and possibly added to it intrinsically during the collapse (Visser & Brouwers 2022), the forming object will become a binary in order to absorb excess angular momentum. In this scenario, the disk is expected to be loaded with pebble-sized objects, with Stokes numbers between 10⁻² to unity.

Regardless of the process of formation, if these binary systems formed early in the disk, their further evolution and growth is non-trivial. The presence of nebular gas and/or solids might cause the binary components to spiral closer and merge if they are aerodynamically small within a Myr (Lyra et al. 2021) or start wide and merge through secular effects (Rozner et al. 2020). However, if the binary components form early and start massive enough for gravity to dominate, they might continue to grow through co-accretion. The growth of such a binary system is widely known to lead to a convergence in the mass ratio – if the satellite and primary accrete from a planetesimal swarm under the right conditions. A satellite would grow significantly faster than its primary due to its higher surface-to-mass ratio and by exploiting the gravitational focusing of the primary (Harris & Kaula 1975), provided that the mutual separation is small enough (Morishima & Watanabe 2001).

While this co-accretion scenario for binaries has been studied for planetesimal accretion, it has not been applied in the case of pebble accretion. In pebble accretion, massive bodies efficiently accrete solids through the interplay of gravity and the presence of nebular gas (Ormel & Klahr 2010; Lambrechts & Johansen 2012). As the gas drag removes relative velocity and angular momentum from the solids, they are efficiently captured into the gravitational potential, leading to largely enhanced accretion cross-sections. Pebble accretion has only been invoked explicitly for a single pebble-accreting body to explain many observational difficulties in planet formation with great success, such as the size distributions of asteroids and Kuiper Belt objects (Johansen et al. 2015), formation processes in the cores of gas giants (Levison et al. 2015), TRAPPIST system architecture (Ormel et al. 2017; Schoonenberg et al. 2019), spin axis of solar system bodies (Visser et al. 2020), formation of terrestrial planets (Johansen et al. 2021), and many other applications.

Of course, pebble accretion will not be relevant for all planetary binary objects. The Earth-Moon system formed later, after the disk gas had already been dispersed – so pebbles would have no longer been drifting inwards and the accretion of small rocks would not be enhanced by the interaction with the gas. Also, low-mass binary objects, as in the case of most Kuiper Belt binaries, will not have sufficient gravity to get pebble accretion started in a serious way. The process is, however, applicable for bodies that are sufficiently massive and form so early that the gas disk and pebbles are still around.

In this study, we investigate the evolution and growth of a planetary binary subject to pebble accretion. We subdivide the total mass, M, of a single body in two masses, M₁ and M₂, such that M = M₁ + M₂ with an increasingly diverging mass fraction, f_m = M₁/M₂ ∊ [1, 8]. We determine the timescale to e–fold the mass of the binary compared to the growth time of a single body with the same mass. We then quantify the e-folding time needed to converge the mass ratio to unity by looking at the individual efficiency of accretion from the flux of pebbles swept up by the system. In many cases we find that, in the ‘die-hard’ pebble accretion regime, the lower mass body (in this case M₂) grows faster than the more massive one in essence for the same reasons as in the planetesimal accretion case described above. As a result, the mass ratio can converge back to unity, well within the disk lifetime of ~1 Myr. Pebble accretion provides an explanation for the significant fraction of (near) equal mass observed planetary binaries.

The paper is organised as follows. In Sect. 2, we discuss the model setup and the assumptions we use. In Sect. 3, we describe the results of the numerical simulations. In Sect. 4, we discuss the outcomes and results, after which we summarise the implications and present our conclusions in Sect. 5.

2 Setup and assumptions

2.1 Disk model

In order to see how these binary systems evolve, we begin with a general description of the planet-forming disk in which the binary resides. For the gas temperature and surface density, we assume power-law profiles given by Weidenschilling (1977b) and Hayashi et al. (1985): $$T\left(r_0\right)=170\text{\hspace{0.17em}K}{\left(\frac{r_0}{1\text{\hspace{0.17em}au}}\right)}^{-1/2},$$(1) $$\text{Σ}\left(r_0\right)=1700\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-2}{\left(\frac{r_0}{1\text{\hspace{0.17em}au}}\right)}^{-3/2},$$(2)

where r₀ is the radial distance to the central star. We assume hydrostatic equilibrium, whereby the vertical density structure forms as a Gaussian: $$\rho \left(r_0,z\right)=\frac{\text{Σ}\left(r_0\right)}{\sqrt{2\pi }H\left(r_0\right)}{e}^{-\frac{1}{2}{\left(\frac{z}{H}\right)}^2},$$(3)

with H = c_s/Ω₀ as the scale height, Ω₀ is the Keplerian frequency, $c_{\text{s}}=\sqrt{k_{\text{B}}T/\overline{m}}$ is the local isothermal sound speed, k_B as the Boltzmann’s constant, and $\overline{m}$ the mean molecular weight. The gas moves at a slightly lower velocity than Keplerian, ${\upsilon }_{\text{K}}=\sqrt{G{M}_{\star }/r_0}$ where G is the gravitational constant and M_★ is the mass of the star. Weidenschilling (1977a) gives a dimensionless constant relating v_K and the headwind particles travelling in Keplerian orbit feel. This dimensionless constant, η, can be estimated by looking at the r-dependencies of Eqs. (1)–(3): $${\upsilon }_{\text{hw}}=\eta {\upsilon }_{\text{K}},$$(4) $$\eta =\frac{r_0}{2{\upsilon }_{\text{K}}^2\rho }\frac{\text{d}P}{\text{d}{r}_0}=-\frac{13c_{\text{s}}^2}{8{\upsilon }_{\text{K}}^2},$$(5)

where v_hw is the headwind and P the pressure governed by the ideal gas law. For the adopted disk profile, the headwind has a numerical value of v_hw = −32 m s⁻¹ downward in the co-moving local frame. In our parameter study, the Epstein regime for pebble stopping time is the relevant one and we calculated it from the Stokes number (τ_s), relating the stopping time (t_s) to one orbital timescale: $${\tau }_{\text{s}}=t_{\text{s}}{\text{Ω}}_0.$$(6)

2.2 Pebble and binary dynamics

We adopted a two-dimensional (2D) local shearing sheet box co-moving with the binary center-of-mass (barycenter) at distance, r₀, from the central star. In this frame the equation of motion of a pebble is given by (Ormel & Klahr 2010): $$\frac{\text{d}v}{\text{d}t}=\left(\begin{array}{c}2{\text{Ω}}_0{\upsilon }_y+3{\text{Ω}}_0^{2}x\\ -2{\text{Ω}}_0{\upsilon }_x\end{array}\right)-{\displaystyle \sum _{j=1}^2\frac{G{M}_j}{{\left|r-r_j\left(t\right)\right|}^3}\left(\begin{array}{c}x-x_j\left(t\right)\\ y-y_j\left(t\right)\end{array}\right)+a_{\text{drag}},}$$(7)

with the first expression in brackets consisting of the Coriolis and tidal accelerations, the second term with the two-body gravity of the separate binary components with mass, M_j, and Cartesian positions (x_j(t), y_j(t)).

From here on we use the following terminology for the binary components. We say the binary is equal mass if the mass ratio is unity and no explicit reference to one of the bodies has to be made due to symmetry. If the mass ratio changes to f_m ≡ M₁/M₂ > 1 we say that the more massive one with mass, M₁, is the ‘primary’, and the least massive one with mass, M₂, is the ‘secondary’.

We consider a binary separation that is narrow and we assume mutual circular motion of the binary system during the accretion. The circular motion of the system with total mass, M = M₁ + M₂, is then modelled over time, t, with: $$\left(\begin{array}{c}{x}_{1,2}\\ y_{1,2}\end{array}\right)=a_{1,2}\left(\begin{array}{c}\cos \left({\phi }_0+\omega t\right)\\ \sin \left({\phi }_0+\omega t\right)\end{array}\right),$$(8)

where $\omega =\sqrt{GM/a_b^3}$ as the angular frequency, φ₀ as the initial phase angle with which the simulation started and the subscripts indicating to which body it applies, respectively. The binary components are assumed massive and unaffected by the gas. They are placed on a circular orbit around the barycenter with orbital separation a_b = a₁ + a₂ and, if f_m ≠ 1, with a₁ as the distance of the primary with mass, M₁, and radius, R₁, and a₂ as the distance of the secondary with mass, M₂, and radius, R₂, to the barycenter. The system rotates in a prograde orientation around its barycenter.

The Cartesian coordinates and velocities of the pebble are given by (x, y, v_x, v_y). The last term contains the acceleration due to the drag force on the pebble equal to: $$a_{\text{drag}}=-\frac{v-v_{\text{g}}}{t_{\text{s}}},$$(9)

with v_g the gas velocity including shear given by: $$v_{\text{g}}=\left({\upsilon }_{\text{hw}}-\frac{3}{2}{\text{Ω}}_{0}x\right)e_y.$$(10)

An illustrative sketch of the setup is shown in Fig. 1.

The total mass, M, of the system in our model corresponds to the mass of the Pluto-Charon system, with an internal density ρ_• = 1.778 g cm⁻³ taken to be equal for both bodies. The (shared) Hill radius is then given by: $$R_{\text{H}}=r_0{\left(\frac{M}{3M_{\star }}\right)}^{1/3}.$$(11)

In varying the mass ratio of the binary system, we then subdivided the total mass accordingly over the components from which their barycentric distance, a₁ and a₂, directly follow as: $$a_{1,2}=\frac{M_{2,1}}{M}{a}_{\text{b}}.$$(12)

From the assumed internal density the radius of each body is then calculated from: $$R_{1,2}={\left(\frac{3M_{1,2}}{4\pi {\rho }_{•}}\right)}^{1/3}.$$(13)

Fig. 1 Overview of the binary planetesimal system in the co-moving frame. The star is located on the left side of the frame. The binary components are denoted with their masses, M1, M2, and radii, R1, R2, travelling around their barycenter in circular orbit at distance, a1 and a2, respectively. Pebbles are released from a distance y0 between the edges of the impact range, x1, x2, for a given initial phase angle of the binary, φ0, with respect to the x-axis. The starting velocities of different pebbles differ depending on the x-coordinate as seen in Eq. (16), therefore, their velocity vectors point in different directions. The background gas headwind is denoted with the black arrow on the y-axis.

2.3 Numerical method and initial conditions

We are interested in the (relative) growth rate of the binary (components) subject to pebble accretion. To do so we integrate the pebble trajectories with the Runge-Kutta-Fehlberg step variable integration method (Fehlberg 1969; Es-hagh 2005), with a relative error tolerance of tol = 10⁻⁶. The pebbles are released on a distance far away from the binary in the y-direction: $$y_0=C\sqrt{\frac{GM{t}_{\text{s}}}{{\upsilon }_{\text{hw}}}},$$(14)

(as illustrated in Fig. 1) with C = 200, a safety factor. This ensures that, initially, gas drag is dominant over the gravity from the binary center-of-mass. The initial velocities of the incoming pebbles are, hence, given by the unperturbed radial and azimuthal drift velocities (Weidenschilling 1977a): $${\upsilon }_{x,\infty }=-\frac{2{\tau }_{\text{s}}}{{\tau }_{\text{s}}^2+1}{\upsilon }_{\text{hw}},$$(15) $${\upsilon }_{y,\infty }=-\frac{1}{{\tau }_{\text{s}}^2+1}{\upsilon }_{\text{hw}}-\frac{3}{2}{\text{Ω}}_0{x}_0,$$(16)

where $-\frac{3}{2}{\text{Ω}}_0{x}_0$ is the correction to the Keplerian Shear. The pebbles are released with a shear corrected sampling function, equivalent to the method used in Visser et al. (2020). In the co-moving local frame we adopt, the pebble flux released then increases for increasing x₀, to account for the increasing relative velocities of the pebbles with respect to the binary barycenter due to shear.

We then use a bisection algorithm to find the range of initial x-coordinates for which the pebbles hit one of the binary components, this range is verified to be equal to the impact parameter of a single planetesimal with mass, M = M₁ + M₂, in the parameter space we investigate. We do note that this only is valid for the close binary separations we consider. A pebble is registered as a hit when it satisfies the condition: $$\sqrt{{\left(x-x_{1,2}\right)}^2+{\left(y-y_{1,2}\right)}^2}\le R_{1.2},$$(17)

with the subscript again referring to the corresponding body. After the cross-section edges have been found, we release pebbles between the edges x₀ ∈ [x₁, x₂] accounting for an error of missing trajectories of ~5 percent. The resolution in pebbles released is taken to be $N_{x_0}=90$ for each initial phase angle linearly spaced between φ₀ ∈ [0, 2π) of the binary system with resolution $N_{{\phi }_0}=90$. In total, we then released the amount of $N_{x_0}\times N_{{\phi }_0}=8100$ pebbles per individual simulation. Finally, in each simulation, we summed the total number of accreted pebbles on each body N_1,2 over each initial phase angle and we defined the accretion ratio to be: $$\epsilon \equiv \frac{{\dot{M}}_1}{{\dot{M}}_2}=\frac{N_1}{N_2},$$(18)

with Ν₁ and Ν₂ as the total hits on the primary or secondary, respectively, but when the mass ratio is not unity.

To structure the results from our parameter study, we defined a fiducial model based on the Pluto-Charon mass(-ratio), radii, orbital distance from the star r₀ = 39 au and separation distance a_b = 0.0025 R_H, with a pebble Stokes number of τ_s = 10⁻². An overview of the fiducial model and the rest of our parameter study is provided in Table 1. A reproduction package of the code can be found online².

3 Results

3.1 General interpretation

We describe above the way in which a single pebble interacts with a binary planetesimal system, however, that does not explain how that system of planetesimals is expected to grow over time. For that purpose, we need every starting x0-position where released pebbles will interact with the systemand all starting angles φ₀ ∈ [0, 2π), as described in Sect. 2.3. In Fig. 2, we plot a hitmap of 500.000 pebbles using our fiducial system from Table 1, with the only exception the mass ratio is set to f_m = 1. Because of symmetry considerations, the accretion ratio ɛ (Eq. (18)) in this run is therefore unity.

Every pixel in the hitmap is a single pebble falling into the binary system, set by the x-coordinate when released, and the phase of the binary at that moment. When a pebble falls toward the system, there are three possible outcomes; the pebble can fall onto either of the planetesimals or it can fly by without hitting any of the two bodies and continue on its path around the central star. In the first two cases, the pixel is coloured either brown or orange depending on onto which body it accretes. If the trajectory does not lead to accretion, it is coloured in beige. This specific case describes an equal-mass binary, M₁ = M₂, the planetesimals are moving as a mirror image of each other around the barycenter. Therefore, we expect that the top and the bottom half of the hitmap are identical, with inverted colours. This is indeed the case in Fig. 2. A horizontal dashed red line is added at φ₀ = π to emphasise this symmetry. The exact same image can be seen above and below this line, only with the brown and orange colours swapped. In a system of a non-equal mass binary, this symmetry is expected to disappear.

By looking at the different trajectories found in the hitmap shown above, we can identify (at the simplest) three regions for all pebbles falling into a binary system, ranging from the most outer left to the most outer right. These three regions are as follows.

In the first, most pebbles falling in from the left side of the map enter the system circling in a prograde direction around the binary, in the same direction as the planetesimals are moving. This can result in a pebble spiralling into the system slowly (as in trajectory E), but as seen in the example of Fig. 2, it mostly results in chaotic paths as seen in trajectories A and B. In this case, we believe these chaotic paths happen because when a pebble gets close to one planetesimal it stays in its vicinity for longer and has a high probability of getting slung around, which would result in a chaotic path. In other scenarios (different τ_s, a_b, etc.), we can see pebbles spiralling in slowly and favouring hitting the outer (less massive) planetesimal, as in Fig. 3.

In the second, as shown in the middle of the hitmap (in Fig. 2, specifically between 0.28R_H ≲ x₀ ≲ 0.32R_H), the pebbles fall in from straight above the system. Falling in head on, with great velocity, either hitting one of the two bodies or going right past them and missing both. This region is exemplified by trajectories C and D.

In the third, shown on the right side, the pebbles come into the system circling around in retrograde fashion, opposing the movement of the binary. In the hitmap, this results in a large region of pebbles spiralling in slowly as seen in E. This happens all the way to ~0.37R_H, from which the last trajectories are mostly chaotic paths. We believe this to result from the fact that the pebbles swirl around the binary and come back into the system almost from below. In that case, it would not have the same momentum as when it comes in from above since it now is moving opposite to the shearing wind. This lack of momentum makes it easy to get swung around when actually getting near the two bodies and do more chaotic paths, as exemplified in F.

Modifying the fundamental parameters τ_s, a_b, and mass ratio will affect where exactly these three regimes (prograde, head on, and retrograde) are, but the fundamental structure remains similar.

Fig. 2 Hitmap of pebble impacts on an equal mass binary system such that M1 = M2 = (MPluto + MCharon)/2, shown at the top. Every pixel is a pebble falling into the system. On the x-axis the x0 coordinate at moment of release is shown in Hill radii, while on the y-axis the initial phase orientation of the binary φ0 in radians is depicted as illustrated in Fig. 1. The orange pixels show impacts on M1, the brown pixels on M2 and the beige pixels indicate ejected pebbles. The final accretion ratio ε = 1, meaning that both bodies accrete an equal flux of pebbles because of symmetry arguments. The dashed red line indicates a starting angle of φ0 = π rad. From 6 results in this graph the trajectories have been plotted below, those are denoted as A-F in this hitmap. The resolution has been taken $\left(N_{x_0}\times N_{{\phi }_0}=500\times 1000\right)$. Bottom panel: six trajectories from the hitmap above. Shown is a top-down view of the binary system with a pebble “falling” in. In each of the trajectories, an arrow denotes the pebble coming into the frame. Pebble trajectories A through F are sorted with their starting distance, x0, going from left to right, corresponding to closer radial distance from the star to further away, respectively. We note that even though the trajectory for F appears to start on the bottom left, it originated on the right outside of the frame only entering the image below.

3.2 Spiralling pebbles

The situation where a pebble is caught by the binary into a spiralling motion is particularly interesting. For an equal-mass binary, both bodies are then equally likely to accrete the pebble. However, if the binary has one body more massive than the other and if the distance between the two bodies is smaller than the initial size of the spiral, the lower-mass body becomes the more likely target. This important effect is illustrated in Fig. 3. In this case, the masses are different, specifically, M₁ = 2M₂. Integrated over all pebbles in the impact parameter, the accretion ratio ε as described in Eq. (18) of this system is less than one, meaning that the smaller body is accreting more mass per unit time than its more massive counterpart. This may seem counter-intuitive as the less massive body has a smaller gravitational pull. However, it can be understood by the in-spiralling motion of the pebbles that are bound to the binary system long before they actually get accreted. Figure 3 shows a number of pebbles from a horizontal row in a hitmap. These pebbles have different starting positions x₀, but are launched with the same starting angle φ₀ of the binary. Many of the orbits turn into a spiral that slowly closes in on the binary. The small body moves like a vacuum cleaner through its orbit around the massive body (in fact, the barycenter, which will be very close to or even inside the larger body). As long as the settling velocity on the spiral path is low, the smaller body is likely to sweep up the incoming pebble before it can reach the massive body.

Fig. 3 Selection of pebble trajectories for a given initial phase angle for the fiducial model with a mass ratio of 2, to illustrate the accretion advantage for the secondary. Left panel: pebbles enter the Hill radius (dashed circle) from the top and approach the binary center of mass in the middle of the panel. Gray trajectories miss the entire system, blue ones accrete onto the secondary and red ones onto the primary. Right panel: zoom in on the trajectories shown in the left panel. The hits on the secondary are indicated with the blue crosses at point of impact. It clearly illustrates the accretion advantage for the smaller body due to close encounters and spirals. Additionally, we show the more massive body in the center (dark red solid circle) with a direct hit of a pebble (dark red trajectory). The trajectories are from different times and are plotted in one panel only for illustration.

3.3 Accretion efficiency

The impact parameter range of pebbles reaching the binary turns out to be the same as in a computation where the total mass of the binary is present as a single, larger body. This is the case when the system has a small binary separation compared to the range of impact parameters in which pebble accretion is relevant. Nevertheless, not all pebbles within the impact parameter range are actually accreted by the binary. The pebbles that are lost escape through scattering by close encounters with the primary and the secondary mass, as illustrated in Fig. 2 (A, D, and F). This is different from the single mass case, in which, for the considered parameter space, all pebbles in the release range would be accreted. Therefore, we define the fraction of pebbles accreted on the binary, compared to the single-body case: $$\text{Γ}\equiv \frac{N_1+N_2}{N_{\text{all}}},$$(19)

where N_all is the total number of released pebbles and N_1,2 are the number of pebbles accreted on bodies 1 and 2, respectively. For a single body with the total binary mass this fraction, Γ, would be unity. The resulting Γ values for our parameter study are shown in Fig. 4 for 2.5 and 39 au. The accretion efficiency can drop to as low as Γ ~ 0.4 for the maximum Stokes numbers considered at 39 au.

Knowing the accretion efficiency Γ of the binary system, we can relate the growth time of the binary system to the growth time of a single body with the same mass. That is, in the pebble accretion regime, the growth timescale, t_growth, of a single body with radius, R, internal density, ρ_•, and mass, M, is given by (Visser & Ormel 2016): $$t_{\text{growth}}=\frac{M}{\dot{M}}=\frac{4{\rho }_{•}R}{3{\upsilon }_{\text{hw}}{\rho }_{\text{p}}{f}_{\text{coll}}},$$(20)

where ρ_p is the spatial density of solid particles and f_coll the collision factor. In the next section, we determine the relative growth timescales of the binary components based on this analysis.

3.4 Relative growth rates

To see if the binary mass ratio will shift towards unity during growth by pebble accretion, we need to look at the normalised growth ratio $$\tilde{\epsilon }=\frac{M_2}{M_1}\epsilon .$$(21)

Whenever this ratio $\tilde{\epsilon }$ dips below one, the binary will grow towards equal mass. In the top parts of Figs. 5 and 6, values of $\tilde{\epsilon }$ is shown as a function of the pebble Stokes number τ_s for various combinations of our key parameters. There is a clear regime, namely, an interval of Stokes numbers, where the mass ratio is indeed moving towards equal mass. This is both true in the asteroid belt at 2.5 au, and in the Kuiper belt at 39 au. This effect is seen not only for a moderate mass ratio of f_m = 2, but also in the more extreme case of f_m = 8. It becomes more prevalent when the two objects in the binary are closer together. In Figs. 5 and 6 it can be seen that while this is not always the case, it is true for a substantial part of the parameter study, in particular for large Stokes numbers. In this parameter regime, pebble accretion works toward reducing the mass ratio of the binary. The question arises as to how relevant this effect is. We consider whether bodies spend enough time in this regime to have their mass ratio significantly altered. In fact, it comes down to the question of what the growth timescales of these bodies are, where the growth timescale is the time to e-fold the mass of a component.

With the growth timescale derived for a single body and the accretion efficiency Γ we can derive the mass accreted by each body over a set amount of time: $${\dot{M}}_1+{\dot{M}}_2=\text{Γ}\dot{M},$$(22)

which can be rewritten as: $${\dot{M}}_1=\frac{\text{Γ}\dot{M}}{\left(1+{\epsilon }^{-1}\right)},$$(23) $${\dot{M}}_2=\frac{\text{Γ}\dot{M}}{\left(1+\epsilon \right)},$$(24)

since $\epsilon ={\dot{M}}_1/{\dot{M}}_2$. We now can use the mass ratio, f_m, such that $M_1=\frac{f_m}{1+f_m}M$ and $M_2=\frac{1}{1+f_m}M$. In this way, we find the growth timescales of both planetesimals: $$t_{\text{growth,}1}=\frac{M_1}{{\dot{M}}_1}=\frac{\left(1+{\epsilon }^{-1}\right)f_m}{\left(1+f_m\right)\text{Γ}}{t}_{\text{growth}}\equiv {\chi }_1{t}_{\text{growth}},$$(25) $$t_{\text{growth,2}}=\frac{M_2}{{\dot{M}}_2}=\frac{\left(1+\epsilon \right)}{\left(1+f_m\right)\text{Γ}}{t}_{\text{growth}}\equiv {\chi }_2{t}_{\text{growth}},$$(26)

where we have defined χ_i to be the scale factor that compares the i_growth of a single body to that of both planetesimals in a binary. The normalised accretion ratio $\tilde{\epsilon }$ and the growth timescales of both bodies are shown in the top and in the bottom panels of Figs. 5 and 6, respectively. At both 2.5 and 39 au (as we see in Fig. 4), the efficiency goes down with increasing, τ_s, which can also be seen with the accretion ratio ɛ. We believe both effects have somewhat the same cause: ɛ favours the smaller body since it sweeps up the inspiralling pebbles from the larger mass (as described in Sect. 3.2), while the efficiency will also go down as the smaller mass throws a significant fraction of these inspiralling pebbles out of the system in close encounters.

What can also be clearly seen is that the relative growth timescale of the larger body is greater than one in almost all of the parameter space, indicating that the growth of a massive body in the pebble accretion regime is slowed down when it has a secondary. In a large part of the parameter space, the χ₂ value of the smaller body drops below 1, indicating that its growth will be accelerated by the presence of the larger partner. Even in the region where the efficiency Γ decreases substantially, the growth of the smaller body is faster. That, combined with the fact the larger body grows more slowly, shows that we have a stabilising effect on the mass ratio. Similar-mass binaries are made more common at both 2.5 and 39 au by this process.

For Stokes numbers close to 0.1, which happens to be also the preferred Stokes number regime for triggering the streaming instability (Bai & Stone 2010), the growth timescales of the smaller body can easily be a factor of 3-5 less than that of the larger body. So in the time the mass of the large body grows by a factor e, the mass of the smaller body may grow by a factor e³ to e⁵, and the mass ratio would change by a factor e² to e⁴ or about 7 to 40. This clearly demonstrates the strong potential of pebble accretion to create equal-mass binaries, even if the initial binary had demonstrated a high mass ratio.

Fig. 4 Growth efficiency, Γ, as described in Eq. (19). Left: mass ratio of fm = 2. Right: mass ratio of fm = 8. for a range of Stokes numbers, τs, mass-ratios, and binary separation, ab, at orbital distance of r0 = 2.5 au from the central star (top panels) and r0 = 39 au (bottom panels). The most opaque lines denote an orbital separation of ab = 0.0075RH while the most transparent lines shows an orbital separation of ab = 0.05RH, other lines indicate orbital separations going logarithmic in between.

Fig. 5 Normalised growth ratios expressed in $\tilde{\epsilon }$ (as described in Eq. (21)) for a range of Stokes numbers τs, mass-ratios, and binary separation, ab, at orbital distance of r0 = 2.5 au from the central star, shown at the top. The most opaque lines denote an orbital separation of ab = 0.0075RΗ, while the most transparent lines show an orbital separation of ab = 0.05RH; other lines indicate orbital separations going logarithmically in between. Relative growth timescales of both bodies as a fraction of the growth timescale of a single body with mass, M = M1 + M2 (as described in Eqs. (25) and (26)), shown at the bottom. Black and red lines denote M1, while the blue lines denote M2.

Fig. 6 Normalised growth ratios expressed in $\tilde{\epsilon }$ (as described in Eq. (21)) for a range of Stokes numbers τs, mass-ratios, and binary separation, ab, at orbital distance of r0 = 39 au from the central star, shown at the top. The most opaque lines denote an orbital separation of ab = 10−3RH, while the most transparent lines show an orbital separation of ab = 10−2RH; other lines indicate orbital separations going logarithmically in between. Relative growth timescales of both bodies as a fraction of the growth timescale of a single body with mass, Μ = M1 + M2 (as described in Eqs. (25) and (26)), shown at the bottom. Black and red lines denote M1, while the blue lines denote M2.

4 Discussion

4.1 Pebble trajectory regimes

We have identified three different main regimes that determine the enhancement in the accretion rate of the secondary (excluding chaotic paths that pebbles may take, as shown in Figs. 2 A, B, and F):

1. The Stokes number, τ_s, is too low and pebbles are tightly coupled to the gas. The pebbles do not spiral but make a direct impact or typically one loop before impact (Fig. 7a). The accretion ratio is in favour of the primary.

2. When the binary separation is too wide, leading to the massive body accreting pebbles in isolation from the smaller one (Fig. 7b), the accretion ratio is in favour of the primary.

3. When τ_s is in the regime of pebble accretion and the separation falls within or on the typical radial extent of the spirals. The pebble undergoes tightly wound long-term spiralling inward and inevitably encounters the secondary, before reaching the more massive one (Fig. 7c). Here, the accretion ratio is in favour of the secondary.

Unsurprisingly, we find the same trend persists for varied orbital distance r₀. We show the final accretion efficiency for a distance of r₀ = 2.5 au and have verified it shows the same trend for r₀ = 10 au in the same manner as Fig. 6. While the parameters shift somewhat, the results show the same outcome.

Fig. 7 Sketch of three distinctive cases that characterise the accretion efficiency of the binary components. The red dot indicates the binary barycenter, the small and large black dots the secondary and primary, respectively, and circular velocities shown by the black arrows (illustrative, and not shown to scale). Case a: pebble Stokes number τs is low and strongly coupled to the gas. It impacts the primary after typically a single encounter and renders the probability low that it will impact the secondary. Case b: orbital separation, ab, is too wide for the secondary to benefit from accretion on the massive one. Case c: narrow binary separation combined with the pebble accretion regime in which pebble spirals are ample and dominant. The probability of impact on the secondary is high due to the repeating encounters before the pebble sinks toward the massive body.

4.2 Timescale for developing toward equal mass

It is informative to see non-equal mass binaries growing towards a mass ratio of unity, but it does not give us a definitive timescale of how long it will take to get there. We consider now the timescale to change the mass ratio, f_m, by a factor e: $$\begin{array}{ll}{t}_{\text{em}}\hfill & =\frac{f_m}{{\dot{f}}_m}=\left(\frac{1}{1-{\tilde{\epsilon }}^{-1}}\right)\frac{M_1}{{\dot{M}}_1}=\frac{\tilde{\epsilon }}{\tilde{\epsilon }-1}{t}_{\text{growth,1}}\hfill \\ \hfill & =\frac{\tilde{\epsilon }{f}_m+1}{\left(\tilde{\epsilon }-1\right)\left(1+f_m\right)}\frac{t_{\text{growth}}}{\text{Γ}}.\hfill \end{array}$$(27)

with $\tilde{\epsilon }=\epsilon /f_m<1$ and both ${\dot{f}}_m$ and t_em being negative since we are looking at the time it takes to grow towards equal mass. This equal mass timescale t_em is plotted in Fig. 8 for our parameter study. We are in an environment where the streaming instability is active or has been active recently, resulting in significant dust settling in the disk. We therefore take for the solid-to-gas ratio a conservative value of 0.1 such that ρ_p from Eq. (20) is given by: $${\rho }_p=0.1{\rho }_g=0.1\frac{\text{Σ}}{\sqrt{2\pi }H},$$(28)

where ρ_g is the density of the gas (see Eq. (3)). We plot this t_em in Fig. 8. For a binary with a total mass equal to that of the Pluto-Charon system at the location of the asteroid belt at 2.5 AU, the timescales to e-fold towards unity are well within a Myr. Even at a greater distance of 39 AU, in a substantial part of our parameter space, we find timescales shorter than 1 Myr, again implying that pebble accretion has true potential as a formation method of (nearly) equal-mass binaries. Even though we have taken the static approach of calculating t_em by looking at a fixed, rather than evolving, mass ratio, from the figure we can conclude that t_em becomes lower with increasing mass ratio, f_m. Hence, the mass-equality timescale slows down while this process is happening, but initially the masses are rapidly converging.

A question arises regarding whether the timescales in Fig. 8 hold for smaller binaries than Pluto-Charon. Since t_growth of a single body at 3 AU is already within a few Myrs once its radius is larger than ~300 km (Visser & Ormel 2016), we estimate the t_em of binaries with bodies between that and Pluto (~1200 km) to be on the same order as shown in Fig. 8.

Fig. 8 Timescale of e-folding towards equal mass tem as described in Eq. (27) for two tested mass ratios, fm = [2, 8], at two different orbital distances, r0 = [2.5, 39] AU, and at different binary separations, ab, denoted by the opaqueness of the lines. We note only the values, where $\tilde{\epsilon }<1$ are plotted since we are looking at the times it takes for the mass ratio, fm, to grow towards unity, not away from it.

4.3 Limitations of the model

One important simplification we have made is the assumption of 2D accretion, only looking at particles located in the midplane and ignoring the effects of a vertical distribution among incoming pebbles. We can speculate that dropping this assumption may reduce the efficiency of the ‘vacuum cleaning’ by the quickly orbiting secondary. An exploration of this effect would be very useful and should be done. For now, we hypothesise that even though the specific values of growth ratios may differ from the 3D results, the general trends will remain. The parameter space of Stokes numbers, τ_s, where the growth ratio flips in favour of the smaller body, will also remain largely the same. Therefore, we think the three main regimes of determining accretion rates we found in binaries will also remain the same as in our 2D case.

We are ignoring the effect of gas drag on the binary itself. While this assumption is certainly correct during the computation of the trajectory of a single pebble, over long evolutionary times it may affect the orbit of the binary by removing angular momentum from it.

In this work, we only considered purely circular orbits in our binary systems. Since truly circular orbits are rarely ever seen, our results might differ if we also simulate elliptical orbits. Most observed binaries that are in close orbits have been tidally circularised (Noll et al. 2020) and we find that the assumption for circular orbits is justified. We estimate, however, that the trends will again stay largely the same for eccentric binaries due to the spiralling nature of pebbles in the ‘die-hard’ pebble accretion regime.

We have only considered prograde orbital orientations for the binary mutual motion. The reason for this is that binaries forming in the context of the streaming instability show a clear preference for prograde orbits and is also the predominant preference in observations of the cold classical Kuiper belt population (Grundy et al. 2019; Nesvorny et al. 2020; Visser & Brouwers 2022).

5 Conclusions

In this paper, we study the growth and evolution of a binary planet (or planetesimal) system subject to pebble accretion. Specifically, we have investigated the accretion efficiency of the individual binary companions and the binary system as a whole for a variety of parameters such as Stokes number, binary separation, and binary mass ratio. If the binary components are pebble accreting with a mass ratio of unity, they accrete pebbles at an equal rate, as expected from symmetry arguments. If the mass ratio is not unity, we arrive at the following conclusions:

1. A single planetary body that is accreting pebbles will generally accrete more efficiently in a binary system with a more massive companion, if pebbles are slowly spiralling toward the binary center-of-mass and these spirals extend to the length of the binary separation distance of a_b or larger.

2. In that case, the mass ratio starts off above unity and the growth rate of the smaller body can be accelerated whilst the more massive one grows more slowly, pushing the binary system toward a mass ratio of unity. For Stokes numbers 0.01 ≤ τ_s ~ 0.1, the growth timescale of the smaller body can be three to five times shorter, changing the mass ratio quickly.

3. Generally, three regimes can be distinguished: (i) Stokes numbers are too small and pebbles are tightly coupled to the gas. These pebbles impact the more massive body, typically after a single encounter, with no accretion advantage for the secondary. (ii) The same holds when Stokes numbers are too large and pebbles meet the binary on ballistic orbits. (iii) Pebbles are captured by the primary in orbits wider than the binary separation. During the in-spiralling, the pebble has a high probability to be swept up by the secondary, resulting in a converging mass ratio.

Acknowledgements

We would like to thank J. Groot, G. Koning, R. M. van der Linden and M. Rustenburg for their insightful comments and discussions. We thank everyone at the disk meeting for feedback and improvements to the manuscript. R.G.V. acknowledges funding from the Dutch Research Council, ALWGO/15-01.

References

All Tables

All Figures

Fig. 1 Overview of the binary planetesimal system in the co-moving frame. The star is located on the left side of the frame. The binary components are denoted with their masses, M1, M2, and radii, R1, R2, travelling around their barycenter in circular orbit at distance, a1 and a2, respectively. Pebbles are released from a distance y0 between the edges of the impact range, x1, x2, for a given initial phase angle of the binary, φ0, with respect to the x-axis. The starting velocities of different pebbles differ depending on the x-coordinate as seen in Eq. (16), therefore, their velocity vectors point in different directions. The background gas headwind is denoted with the black arrow on the y-axis.

In the text

Fig. 2 Hitmap of pebble impacts on an equal mass binary system such that M1 = M2 = (MPluto + MCharon)/2, shown at the top. Every pixel is a pebble falling into the system. On the x-axis the x0 coordinate at moment of release is shown in Hill radii, while on the y-axis the initial phase orientation of the binary φ0 in radians is depicted as illustrated in Fig. 1. The orange pixels show impacts on M1, the brown pixels on M2 and the beige pixels indicate ejected pebbles. The final accretion ratio ε = 1, meaning that both bodies accrete an equal flux of pebbles because of symmetry arguments. The dashed red line indicates a starting angle of φ0 = π rad. From 6 results in this graph the trajectories have been plotted below, those are denoted as A-F in this hitmap. The resolution has been taken $\left(N_{x_0}\times N_{{\phi }_0}=500\times 1000\right)$. Bottom panel: six trajectories from the hitmap above. Shown is a top-down view of the binary system with a pebble “falling” in. In each of the trajectories, an arrow denotes the pebble coming into the frame. Pebble trajectories A through F are sorted with their starting distance, x0, going from left to right, corresponding to closer radial distance from the star to further away, respectively. We note that even though the trajectory for F appears to start on the bottom left, it originated on the right outside of the frame only entering the image below.

In the text

Fig. 3 Selection of pebble trajectories for a given initial phase angle for the fiducial model with a mass ratio of 2, to illustrate the accretion advantage for the secondary. Left panel: pebbles enter the Hill radius (dashed circle) from the top and approach the binary center of mass in the middle of the panel. Gray trajectories miss the entire system, blue ones accrete onto the secondary and red ones onto the primary. Right panel: zoom in on the trajectories shown in the left panel. The hits on the secondary are indicated with the blue crosses at point of impact. It clearly illustrates the accretion advantage for the smaller body due to close encounters and spirals. Additionally, we show the more massive body in the center (dark red solid circle) with a direct hit of a pebble (dark red trajectory). The trajectories are from different times and are plotted in one panel only for illustration.

In the text

Fig. 4 Growth efficiency, Γ, as described in Eq. (19). Left: mass ratio of fm = 2. Right: mass ratio of fm = 8. for a range of Stokes numbers, τs, mass-ratios, and binary separation, ab, at orbital distance of r0 = 2.5 au from the central star (top panels) and r0 = 39 au (bottom panels). The most opaque lines denote an orbital separation of ab = 0.0075RH while the most transparent lines shows an orbital separation of ab = 0.05RH, other lines indicate orbital separations going logarithmic in between.

In the text

Fig. 5 Normalised growth ratios expressed in $\tilde{\epsilon }$ (as described in Eq. (21)) for a range of Stokes numbers τs, mass-ratios, and binary separation, ab, at orbital distance of r0 = 2.5 au from the central star, shown at the top. The most opaque lines denote an orbital separation of ab = 0.0075RΗ, while the most transparent lines show an orbital separation of ab = 0.05RH; other lines indicate orbital separations going logarithmically in between. Relative growth timescales of both bodies as a fraction of the growth timescale of a single body with mass, M = M1 + M2 (as described in Eqs. (25) and (26)), shown at the bottom. Black and red lines denote M1, while the blue lines denote M2.

In the text

Fig. 6 Normalised growth ratios expressed in $\tilde{\epsilon }$ (as described in Eq. (21)) for a range of Stokes numbers τs, mass-ratios, and binary separation, ab, at orbital distance of r0 = 39 au from the central star, shown at the top. The most opaque lines denote an orbital separation of ab = 10−3RH, while the most transparent lines show an orbital separation of ab = 10−2RH; other lines indicate orbital separations going logarithmically in between. Relative growth timescales of both bodies as a fraction of the growth timescale of a single body with mass, Μ = M1 + M2 (as described in Eqs. (25) and (26)), shown at the bottom. Black and red lines denote M1, while the blue lines denote M2.

In the text

Fig. 7 Sketch of three distinctive cases that characterise the accretion efficiency of the binary components. The red dot indicates the binary barycenter, the small and large black dots the secondary and primary, respectively, and circular velocities shown by the black arrows (illustrative, and not shown to scale). Case a: pebble Stokes number τs is low and strongly coupled to the gas. It impacts the primary after typically a single encounter and renders the probability low that it will impact the secondary. Case b: orbital separation, ab, is too wide for the secondary to benefit from accretion on the massive one. Case c: narrow binary separation combined with the pebble accretion regime in which pebble spirals are ample and dominant. The probability of impact on the secondary is high due to the repeating encounters before the pebble sinks toward the massive body.

In the text

Fig. 8 Timescale of e-folding towards equal mass tem as described in Eq. (27) for two tested mass ratios, fm = [2, 8], at two different orbital distances, r0 = [2.5, 39] AU, and at different binary separations, ab, denoted by the opaqueness of the lines. We note only the values, where $\tilde{\epsilon }<1$ are plotted since we are looking at the times it takes for the mass ratio, fm, to grow towards unity, not away from it.

In the text