11institutetext: 1 Anton Pannekoek institute for Astronomy (API), University
of Amsterdam, Science Park 904, 1098XH Amsterdam
2 Department of Astronomy, Tsinghua University, 30 Shuangqing Rd, 100084
Beĳing, PR China

# Forming equal-mass planetary binaries via pebble accretion

T.J. Konijn1 Corresponding author<EMAIL_ADDRESS>R.G. Visser1 C. Dominik1 C.W.
Ormel2

(Received 10 October 2022 / Accepted 29 November 2022)

###### Abstract

Context. Binary Solar System objects are common, ranging from satellite
systems with very large mass ratios, $M_{1}/M_{2}$, to those with mass ratios
approaching unity. One well-known example of a binary is the Pluto-Charon
system. With Charon being ’only’ eight times less massive than Pluto, the
question arises (as in the case of many other systems) as to why the mass
ratio is still close to unity. There is much evidence that (binary)
planet(esimal) formation happened early, when the protoplanetary gas disk was
still present. It is likely that (at least some of) these binaries evolved
together, as a result of pebble accretion. Pebble accretion is a new key
paradigm in planetary formation and it is believed to play a major role in
many aspects of the formation of planetary systems, from the radial transport
of material to the rapid growth of planetary embryos throughout the system.

Aims. Here, we focus on the question of how the mass arriving in the
gravitational influence zone of the binary during pebble accretion is
distributed over the binary components for a given initial mass ratio. We also
consider whether accretion over time leads to equal-mass binaries (converging
mass ratio) or to a dominant primary mass with a small moon (diverging mass
ratio).

Methods. We numerically integrated two-dimensional (2D) pebble trajectories in
the same typical fashion as for a single mass that is subject to pebble
accretion. We tracked the efficiency of accretion for the two separate binary
components, compared to a single body with the same mass. These numerical
simulations were done for a range of binary mass ratios, mutual separations,
Stokes numbers, and two orbital distances, 2.5 and 39 au.

Results. We find that in the limit where pebbles start to spiral around the
primary (this holds for relatively large pebbles), the pebble preferentially
collides with the secondary, causing the mass ratio to converge towards unity.
In our tested case, where the total binary mass is equal to that of the Pluto-
Charon system, this takes place on $\sim$Myr timescales. In this regime the
total sweep-up efficiency can lower to half that of a pebble-accreting single
body because pebbles that are thrown out of the system, after close encounters
with the system. These timescales and sweep-up efficiency are calculated under
the assumption our 2D simulations compare with the 3D reality. The results
show that systems such as Pluto-Charon and other larger equal mass binaries
may well have co-accreted by means of pebble accretion in the disk phase
without producing binaries, with highly diverging mass ratios.

###### Key Words.:

Pebble Accretion – Binary Planetesimals – Planetary formation – Streaming
instability

##  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) 111The authors classify it as ”equal size binaries”. The
internal densities are very similar for cold classicals making it equivalent
to the definition of equal mass binaries..

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;
Agarwal et al., 2017; Marchis et al., 2012)), 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, 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 &
Nesvorný, 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^{-2}$ 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_{1}$ and $M_{2}$, such that $M=M_{1}+M_{2}$ with an
increasingly diverging mass fraction, $f_{m}=M_{1}/M_{2}\in[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_{2}$) 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 $\sim$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 Section 2, we discuss the model setup
and the assumptions we use. In Section 3, we describe the results of the
numerical simulations. In Section 4, we discuss the outcomes and results,
after which we summarise the implications and present our conclusions in
Section 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):

$\displaystyle T(r_{0})$ $\displaystyle=170\ \mathrm{K}\left(\frac{r_{0}}{1\
\mathrm{au}}\right)^{-1/2},$ (1) $\displaystyle\Sigma(r_{0})$
$\displaystyle=1700\ \mathrm{g\ cm}^{-2}\left(\frac{r_{0}}{1\
\mathrm{au}}\right)^{-3/2},$ (2)

where $r_{0}$ is the radial distance to the central star. We assume
hydrostatic equilibrium, whereby the vertical density structure forms as a
Gaussian:

$\displaystyle\rho(r_{0},z)$
$\displaystyle=\frac{\Sigma(r_{0})}{\sqrt{2\pi}H(r_{0})}e^{-\frac{1}{2}\left(\frac{z}{H}\right)^{2}},$
(3)

with $H=c_{\mathrm{s}}/\Omega_{0}$ as the scale height, $\Omega_{0}$ is the
Keplerian frequency, $c_{\mathrm{s}}=\sqrt{k_{\mathrm{B}}T/\bar{m}}$ is the
local isothermal sound speed, $k_{\mathrm{B}}$ as the Boltzmann’s constant,
and $\bar{m}$ the mean molecular weight. The gas moves at a slightly lower
velocity than Keplerian, $v_{\mathrm{K}}=\sqrt{GM_{\star}/r_{0}}$ where $G$ is
the gravitational constant and $M_{\star}$ is the mass of the star.
(Weidenschilling, 1977a) gives a dimensionless constant relating
$v_{\mathrm{K}}$ and the headwind particles travelling in Keplerian orbit
feel. This dimensionless constant, $\eta,$ can be estimated by looking at the
$r$-dependencies of Eqs. 1-3:

$\displaystyle v_{\mathrm{hw}}$ $\displaystyle=\eta v_{\mathrm{K}},$ (4)
$\displaystyle\eta$
$\displaystyle=\frac{r_{0}}{2v_{\mathrm{K}}^{2}\rho}\frac{\mathrm{d}P}{\mathrm{d}r_{0}}=-\frac{13c_{\mathrm{s}}^{2}}{8v_{\mathrm{K}}^{2}},$
(5)

where $v_{\mathrm{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_{\mathrm{hw}}=-32$ m s-1 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
($\tau_{\mathrm{s}}$), relating the stopping time ($t_{\mathrm{s}}$) to one
orbital timescale:

$\tau_{\mathrm{s}}=t_{\mathrm{s}}\Omega_{0}.$ (6)

Table 1: Overview of the parameter study used in the simulations. The ranges between brackets are logarithmically spaced containing 30 values for $\tau_{\mathrm{s}}$ and 10 values for $a_{\mathrm{b}}$. Since we are simulating two different mass ratios, $f_{m}$, at two different orbital distances $r_{0}$, we have 30x10x2x2 = 1200 simulations in total. The boldfaced quantities indicate the parameter set of the fiducial model. $r_{0}$ [au] | $\tau_{\mathrm{s}}$ |  $a_{\mathrm{b}}$ $[R_{\mathrm{H}}]$ | $f_{m}$ | N | $M_{\mathrm{total}}$
---|---|---|---|---|---
2.5 | [$10^{-2}$-0.4] | [0.0075-0.05] | $2,8$ | 8100 | $M_{\mathrm{Pluto}}+M_{\mathrm{Charon}}$
$\mathbf{39}$ | [$5\times 10^{-4}$-0.1], $\mathbf{10^{-2}}$ | [$10^{-3}$-$10^{-2}$], $\mathbf{0.0025}$ | $2,\mathbf{8}$ | 8100 | $M_{\mathrm{Pluto}}+M_{\mathrm{Charon}}$

###  2.2 Pebble and binary dynamics

Figure 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, $M_{1},M_{2}$, and radii, $R_{1},R_{2}$, travelling
around their barycenter in circular orbit at distance, $a_{1}$ and $a_{2}$,
respectively. Pebbles are released from a distance $y_{0}$ between the edges
of the impact range, $x_{1},x_{2}$, for a given initial phase angle of the
binary, $\varphi_{0}$, with respect to the $x$-axis. The starting velocities
of different pebbles differ depending on the $x$-coordinate as seen in
Equation 17, therefore, their velocity vectors point in different directions.
The background gas headwind is denoted with the black arrow on the y-axis.

We adopted a two-dimensional (2D) local shearing sheet box co-moving with the
binary center-of-mass (barycenter) at distance, $r_{0}$, from the central
star. In this frame the equation of motion of a pebble is given by (Ormel &
Klahr, 2010):

(7) (8) -∑j=12GMj|r-rj(t)|3(x)-xj(t) y-yj(t) +adrag,

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}\equiv M_{1}/M_{2}>1$ we say that the more massive one
with mass, $M_{1}$, is the ’primary’, and the least massive one with mass,
$M_{2}$, 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_{1}+M_{2}$, is then modelled over time, $t,$
with:

$\pmatrix{x}$ (9) y1,2 =a1,2(cos)(φ0+ωt) sin(φ0+ωt) ,

where $\omega=\sqrt{GM/a_{\mathrm{b}}^{3}}$ as the angular frequency,
$\varphi_{0}$ 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_{\mathrm{b}}=a_{1}+a_{2}$ and, if $f_{m}\neq 1$, with $a_{1}$ as the
distance of the primary with mass, $M_{1}$, and radius, $R_{1}$, and $a_{2}$
as the distance of the secondary with mass, $M_{2}$, and radius, $R_{2}$, 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:

$\mathbf{a}_{\mathrm{drag}}=-\frac{\mathbf{v}-\mathbf{v}_{\mathrm{g}}}{t_{\mathrm{s}}},$
(10)

with $\mathbf{v}_{\mathrm{g}}$ the gas velocity including shear given by:

$\mathbf{v}_{\mathrm{g}}=\left(v_{\mathrm{hw}}-\frac{3}{2}\Omega_{0}x\right)\mathbf{e}_{y}.$
(11)

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

The total mass, $M,$ of the system in our model corresponds to the mass of the
Pluto-Charon system, with an internal density
$\rho_{\bullet}=1.778\,\mathrm{g\ cm^{-3}}$ taken to be equal for both bodies.
The (shared) Hill radius is then given by:

$R_{\mathrm{H}}=r_{0}\left(\frac{M}{3M_{\mathrm{\star}}}\right)^{1/3}.$ (12)

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_{1}$ and $a_{2}$, directly follow as:

$a_{1,2}=\frac{M_{2,1}}{M}a_{\mathrm{b}}.$ (13)

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_{\bullet}}\right)^{1/3}.$ (14)

###  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
$\mathrm{tol}=10^{-6}$. The pebbles are released on a distance far away from
the binary in the y-direction:

$\displaystyle y_{0}=C\sqrt{\frac{GMt_{\mathrm{s}}}{v_{\mathrm{hw}}}},$ (15)

(as illustrated in figure 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):

$\displaystyle v_{x,\infty}$
$\displaystyle=-\frac{2\tau_{\mathrm{s}}}{\tau_{\mathrm{s}}^{2}+1}v_{\mathrm{hw}},$
(16) $\displaystyle v_{y,\infty}$
$\displaystyle=-\frac{1}{\tau_{\mathrm{s}}^{2}+1}v_{\mathrm{hw}}-\frac{3}{2}\Omega_{0}x_{0},$
(17)

where $-\frac{3}{2}\Omega_{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_{0}$, 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_{1}+M_{2}$, 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:

$\displaystyle\sqrt{(x-x_{1,2})^{2}+(y-y_{1,2})^{2}}\leq R_{1,2},$ (18)

with the subscript again referring to the corresponding body. After the cross-
section edges have been found, we release pebbles between the edges
$x_{0}\in[x_{1},x_{2}]$ accounting for an error of missing trajectories of
$\sim$ 5 percent. The resolution in pebbles released is taken to be
$N_{x_{0}}=90$ for each initial phase angle linearly spaced between
$\varphi_{0}\in[0,2\pi)$ of the binary system with resolution
$N_{\varphi_{0}}=90$. In total, we then released the amount of
$N_{x_{0}}\times N_{\varphi_{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:

$\varepsilon\equiv\frac{\dot{M_{1}}}{\dot{M_{2}}}=\frac{N_{1}}{N_{2}},$ (19)

with $N_{1}$ and $N_{2}$ 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_{0}=39$ au and separation distance $a_{\mathrm{b}}=0.0025R_{\mathrm{H}}$,
with a pebble Stokes number of $\tau_{\mathrm{s}}=10^{-2}$. 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 online222The reproduction
package can be found on doi.org/10.5281/zenodo.7324045.

Figure 2: Hitmap of pebble impacts on an equal mass binary system such that
$M_{1}=M_{2}=(M_{\mathrm{Pluto}}+M_{\mathrm{Charon}})/2,$ shown at the top.
Every pixel is a pebble falling into the system. On the $x$-axis the $x_{0}$
coordinate at moment of release is shown in Hill radii, while on the $y$-axis
the initial phase orientation of the binary $\varphi_{0}$ in radians is
depicted as illustrated in figure 1. The orange pixels show impacts on
$M_{1}$, the brown pixels on $M_{2}$ & the beige pixels indicate ejected
pebbles. The final accretion ratio $\varepsilon=1$, meaning that both bodies
accrete an equal flux of pebbles because of symmetry arguments. The dashed red
line indicates a starting angle of $\varphi_{0}=\pi$ 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 $(N_{\varphi_{0}}\times
N_{X_{0}}=500\times 1000)$. 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, $x_{0}$, 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. Figure 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 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 $x_{0}$-position where released pebbles will interact with the system
and all starting angles $\varphi_{0}\in[0,2\pi),$ as described in Section 2.3.
In Figure 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 $\varepsilon$
(Equation 19) 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_{1}=M_{2}$, 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 Figure 2. A horizontal dashed red line is added at
$\varphi_{0}=\pi$ 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 Figure 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
$\tau_{\mathrm{s}},a_{\mathrm{b}}$, etc.), we can see pebbles spiralling in
slowly and favouring hitting the outer (less massive) planetesimal, as in
Figure 3.

In the second, as shown in the middle of the hitmap (in Figure 2, specifically
between $0.28R_{\mathrm{H}}\lesssim x_{0}\lesssim 0.32R_{\mathrm{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 $\sim 0.37R_{\mathrm{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 $\tau_{\mathrm{s}},a_{\mathrm{b}}$, and
mass ratio will affect where exactly these three regimes (prograde, head on,
and retrograde) are, but the fundamental structure remains similar.

Figure 4: Growth efficiency, $\Gamma,$ as described in Equation 20. Left: Mass
ratio of $f_{m}=2$. Right: Mass ratio of $f_{m}=8$. for a range of Stokes
numbers, $\tau_{\mathrm{s}}$, mass-ratios, and binary separation,
$a_{\mathrm{b}}$, at orbital distance of $r_{0}=2.5$ au from the central star
(top panels) and $r_{0}=39$ au (bottom panels). The most opaque lines denote
an orbital separation of $a_{\mathrm{b}}=0.0075R_{\mathrm{H}}$ while the most
transparent lines shows an orbital separation of
$a_{\mathrm{b}}=0.05R_{\mathrm{H}}$, other lines indicate orbital separations
going logarithmic in between.

###  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 Figure 3. In this case,
the masses are different, specifically, $M_{1}=2M_{2}$. Integrated over all
pebbles in the impact parameter, the accretion ratio $\varepsilon$ as
described in Equation 19 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_{0}$, but
are launched with the same starting angle $\varphi_{0}$ 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.

Figure 5: Normalised growth ratios expressed in $\tilde{\varepsilon}$ (as
described in Equation 22) for a range of Stokes numbers $\tau_{\mathrm{s}}$,
mass-ratios, and binary separation, $a_{\mathrm{b}}$, at orbital distance of
$r_{0}=2.5$ au from the central star, shown at the top. The most opaque lines
denote an orbital separation of $a_{\mathrm{b}}=0.0075R_{\mathrm{H}}$, while
the most transparent lines show an orbital separation of
$a_{\mathrm{b}}=0.05R_{\mathrm{H}}$; 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=M_{1}+M_{2}$
(as described in Eqs. 26 and 27), shown at the bottom. Black and red lines
denote $M_{1}$, while the blue lines denote $M_{2}$. Figure 6: Normalised
growth ratios expressed in $\tilde{\varepsilon}$ (as described in Equation 22)
for a range of Stokes numbers $\tau_{\mathrm{s}}$, mass-ratios, and binary
separation, $a_{\mathrm{b}}$, at orbital distance of $r_{0}=39$ au from the
central star, shown at the top. The most opaque lines denote an orbital
separation of $a_{\mathrm{b}}=10^{-3}R_{\mathrm{H}}$, while the most
transparent lines show an orbital separation of
$a_{\mathrm{b}}=10^{-2}R_{\mathrm{H}}$; 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=M_{1}+M_{2}$ (as described in Eqs. 26 and 27), shown at the bottom. Black
and red lines denote $M_{1}$, while the blue lines denote $M_{2}$.

###  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 Figure 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:

$\Gamma\equiv\frac{N_{1}+N_{2}}{N_{\mathrm{all}}},$ (20)

where $N_{\mathrm{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, $\Gamma$, would be
unity. The resulting $\Gamma$ values for our parameter study are shown in
Figure 4 for 2.5 au and 39 au. The accretion efficiency can drop to as low as
$\Gamma\sim 0.4$ for the maximum Stokes numbers considered at 39 au.

Knowing the accretion efficiency $\Gamma$ 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_{\mathrm{growth}}$, of a single body with radius, $R$, internal density,
$\rho_{\bullet}$, and mass, $M$, is given by (Visser & Ormel, 2016):

$t_{\mathrm{growth}}=\frac{M}{\dot{M}}=\frac{4\rho_{\bullet}R}{3v_{\mathrm{hw}}\rho_{\mathrm{p}}f_{\mathrm{coll}}},$
(21)

where $\rho_{\mathrm{p}}$ is the spatial density of solid particles and
$f_{\mathrm{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{\varepsilon}=\frac{M_{2}}{M_{1}}\varepsilon.$ (22)

Whenever this ratio $\tilde{\varepsilon}$ dips below one, the binary will grow
towards equal mass. In the top parts of Figures 5 and 6, values of
$\tilde{\varepsilon}$ is shown as a function of the pebble Stokes number
$\tau_{\mathrm{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 Figures 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 $\Gamma$ we can derive the mass accreted by each body over a set
amount of time:

$\dot{M}_{1}+\dot{M_{2}}=\Gamma\dot{M},$ (23)

which can be rewritten as:

$\displaystyle\dot{M}_{1}$
$\displaystyle=\frac{\Gamma\dot{M}}{\left(1+\varepsilon^{-1}\right)},$ (24)
$\displaystyle\dot{M}_{2}$
$\displaystyle=\frac{\Gamma\dot{M}}{\left(1+\varepsilon\right)},$ (25)

since $\varepsilon=\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:

$\displaystyle t_{\mathrm{growth,1}}$
$\displaystyle=\frac{M_{1}}{\dot{M}_{1}}=\frac{\left(1+\varepsilon^{-1}\right)f_{m}}{\left(1+f_{m}\right)\Gamma}t_{\mathrm{growth}}\equiv\chi_{1}t_{\mathrm{growth}},$
(26) $\displaystyle t_{\mathrm{growth,2}}$
$\displaystyle=\frac{M_{2}}{\dot{M}_{2}}=\frac{\left(1+\varepsilon\right)}{\left(1+f_{m}\right)\Gamma}t_{\mathrm{growth}}\equiv\chi_{2}t_{\mathrm{growth}},$
(27)

where we have defined $\chi_{i}$ to be the scale factor that compares the
$t_{\mathrm{growth}}$ of a single body to that of both planetesimals in a
binary. The normalised accretion ratio $\tilde{\varepsilon}$ and the growth
timescales of both bodies are shown in the top and in the bottom panels of
Figures 5 and 6, respectively. At both 2.5 and 39 au (as we see in Figure 4),
the efficiency goes down with increasing, $\tau_{\mathrm{s}}$, which can also
be seen with the accretion ratio $\varepsilon$. We believe both effects have
somewhat the same cause: $\varepsilon$ favours the smaller body since it
sweeps up the inspiralling pebbles from the larger mass (as described in
Section 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.

Figure 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 $\tau_{\mathrm{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, $a_{\mathrm{b}}$, 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.

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 $\chi_{2}$ 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 $\Gamma$ 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^{3}$ to
$e^{5}$, and the mass ratio would change by a factor $e^{2}$ to e4 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.

##  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 Figure 2A, B, and F):

  1. 1.

The Stokes number, $\tau_{\mathrm{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. 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. 3.

When $\tau_{\mathrm{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_{0}$. We show the final accretion efficiency for a distance of $r_{0}=2.5$
au and have verified it shows the same trend for $r_{0}=10$ au in the same
manner as Figure 6. While the parameters shift somewhat, the results show the
same outcome.

###  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$:

$\displaystyle t_{\mathrm{em}}$
$\displaystyle=\frac{f_{m}}{\dot{f_{m}}}=\left(\frac{1}{1-\tilde{\varepsilon}^{-1}}\right)\frac{M_{1}}{\dot{M_{1}}}=\frac{\tilde{\varepsilon}}{\tilde{\varepsilon}-1}t_{\mathrm{growth,1}}$
$\displaystyle=\frac{\tilde{\varepsilon}f_{m}+1}{\left(\tilde{\varepsilon}-1\right)\left(1+f_{m}\right)}\frac{t_{\mathrm{growth}}}{\Gamma}.$
(28)

with $\tilde{\varepsilon}=\varepsilon/f_{m}<1$ and both $\dot{f_{m}}$ and
$t_{\mathrm{em}}$ being negative since we are looking at the time it takes to
grow towards equal mass. This equal mass timescale $t_{\mathrm{em}}$ is
plotted in Figure 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 $\rho_{p}$ from Equation 21 is
given by:

$\rho_{p}=0.1\rho_{g}=0.1\frac{\Sigma}{\sqrt{2\pi}H},$ (29)

where $\rho_{g}$ is the density of the gas (see Equation 3). We plot this
$t_{\mathrm{em}}$ in Figure 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_{\mathrm{em}}$ by looking
at a fixed, rather than evolving, mass ratio, from the figure we can conclude
that $t_{\mathrm{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 Figure 8 hold for
smaller binaries than Pluto-Charon. Since $t_{\mathrm{growth}}$ of a single
body at 3 AU is already within a few Myrs once its radius is larger than
$\sim$300km (Visser & Ormel, 2016), we estimate the $t_{\mathrm{em}}$ of
binaries with bodies between that and Pluto ($\sim$1200km) to be on the same
order as shown in Figure 8.

Figure 8: Timescale of $e$-folding towards equal mass $t_{\mathrm{em}}$ as
described in Equation 28 for two tested mass ratios, $f_{m}=[2,8],$ at two
different orbital distances, $r_{0}=[2.5,39]$ AU, and at different binary
separations, $a_{b}$, denoted by the opaqueness of the lines. We note only the
values, where $\tilde{\varepsilon}<1$ are plotted since we are looking at the
times it takes for the mass ratio, $f_{m}$, 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, $\tau_{\mathrm{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. 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_{\mathrm{b}}$ or larger.

  2. 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\leq\tau_{s}\sim 0.1,$ the growth timescale of the smaller body can be
three to five times shorter, changing the mass ratio quickly.

  3. 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 & M.
Rustenburg for their insightful comments and discussions. We thank everyone at
the disk meeting for feedback and improvements to the manuscript. RGV
acknowledges funding from the Dutch Research Council, ALWGO/15-01.

## References

  * Agarwal et al. (2017) Agarwal, J., Jewitt, D., Mutchler, M., Weaver, H., & Larson, S. 2017, Nature, 549, 357
  * Astakhov et al. (2005) Astakhov, S. A., Lee, E. A., & Farrelly, D. 2005, MNRAS, 360, 401
  * Bai & Stone (2010) Bai, X.-N. & Stone, J. M. 2010, ApJ, 722, 1437
  * Benecchi et al. (2009) Benecchi, S. D., Noll, K. S., Grundy, W. M., et al. 2009, Icarus, 200, 292
  * Cameron & Ward (1976) Cameron, A. G. W. & Ward, W. R. 1976, in Lunar and Planetary Science Conference, Vol. 7, Lunar and Planetary Science Conference, 120
  * Canup (2005) Canup, R. M. 2005, Science, 307, 546
  * Canup (2011) Canup, R. M. 2011, AJ, 141, 35
  * Canup (2012) Canup, R. M. 2012, Science, 338, 1052
  * Christy & Harrington (1978) Christy, J. W. & Harrington, R. S. 1978, AJ, 83, 1005
  * Es-hagh (2005) Es-hagh, M. 2005, Journal of the Earth & Space Physics, 31, 1
  * Fehlberg (1969) Fehlberg, E. 1969, NASA Technical Report, NASA-TR-R-315
  * Goldreich et al. (2002) Goldreich, P., Lithwick, Y., & Sari, R. 2002, Nature, 420, 643
  * Grundy et al. (2019) Grundy, W. M., Noll, K. S., Roe, H. G., et al. 2019, Icarus, 334, 62
  * Harris & Kaula (1975) Harris, A. W. & Kaula, W. M. 1975, Icarus, 24, 516
  * Hayashi et al. (1985) Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 1100–1153
  * Herschel (1833) Herschel, J. F. W. 1833, A Treatise on Astronomy
  * Johansen & Klahr (2011) Johansen, A. & Klahr, H. 2011, Earth Moon and Planets, 108, 39
  * Johansen et al. (2015) Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1, 1500109
  * Johansen et al. (2021) Johansen, A., Ronnet, T., Bizzarro, M., et al. 2021, Science Advances, 7, eabc0444
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Lambrechts & Johansen (2012) Lambrechts, M. & Johansen, A. 2012, A&A, 544, A32
  * Lee et al. (2007) Lee, E. A., Astakhov, S. A., & Farrelly, D. 2007, MNRAS, 379, 229
  * Levison et al. (2015) Levison, H. F., Kretke, K. A., & Duncan, M. J. 2015, Nature, 524, 322
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719
  * Lyra et al. (2021) Lyra, W., Youdin, A. N., & Johansen, A. 2021, Icarus, 356, 113831
  * Marchis et al. (2008) Marchis, F., Descamps, P., Baek, M., et al. 2008, Icarus, 196, 97
  * Marchis et al. (2004) Marchis, F., Descamps, P., Hestroffer, D., Berthier, J., & de Pater, I. 2004, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 36, AAS/Division for Planetary Sciences Meeting Abstracts #36, 46.02
  * Marchis et al. (2012) Marchis, F., Enriquez, J. E., Emery, J. P., et al. 2012, Icarus, 221, 1130
  * Merline et al. (2002) Merline, W. J., Weidenschilling, S. J., Durda, D. D., et al. 2002, in Asteroids III, 289–312
  * Morbidelli & Nesvorný (2020) Morbidelli, A. & Nesvorný, D. 2020, in The Trans-Neptunian Solar System, ed. D. Prialnik, M. A. Barucci, & L. Young, 25–59
  * Morishima & Watanabe (2001) Morishima, R. & Watanabe, S. 2001, Earth, Planets and Space, 53, 213
  * Nesvorny et al. (2020) Nesvorny, D., Li, R., Simon, J. B., et al. 2020, Binary Planetesimal Formation from Gravitationally Collapsing Pebble Clouds
  * Nesvorný et al. (2021) Nesvorný, D., Li, R., Simon, J. B., et al. 2021, Planetary Science Journal, 2, 27
  * Nesvorný et al. (2019) Nesvorný, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nature Astronomy, 3, 808
  * Nesvorný et al. (2010) Nesvorný, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 785
  * Noll et al. (2020) Noll, K., Grundy, W. M., Nesvorný, D., & Thirouin, A. 2020, in The Trans-Neptunian Solar System, ed. D. Prialnik, M. A. Barucci, & L. Young, 201–224
  * Noll et al. (2008) Noll, K. S., Grundy, W. M., Stephens, D. C., Levison, H. F., & Kern, S. D. 2008, Icarus, 194, 758
  * Ormel & Klahr (2010) Ormel, C. W. & Klahr, H. H. 2010, A&A, 520, A43
  * Ormel et al. (2017) Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1
  * Rozner et al. (2020) Rozner, M., Grishin, E., & Perets, H. B. 2020, MNRAS, 497, 5264
  * Russell (2017) Russell, D. G. 2017, International Journal of Astronomy and Astrophysics, 7, 291
  * Schlichting & Sari (2008) Schlichting, H. E. & Sari, R. 2008, ApJ, 673, 1218
  * Schoonenberg et al. (2019) Schoonenberg, D., Liu, B., Ormel, C. W., & Dorn, C. 2019, A&A, 627, A149
  * Stern & Levison (2002) Stern, S. A. & Levison, H. F. 2002, Highlights of Astronomy, 12, 205
  * Visser et al. (2020) Visser, R., Ormel, C., Dominik, C., & Ida, S. 2020, Icarus, 335, 113380
  * Visser & Brouwers (2022) Visser, R. G. & Brouwers, M. G. 2022, A&A, 663, A164
  * Visser et al. (2021) Visser, R. G., Drażkowska, J., & Dominik, C. 2021, A&A, 647, A126
  * Visser & Ormel (2016) Visser, R. G. & Ormel, C. W. 2016, Astronomy & Astrophysics, 586, A66
  * Weidenschilling (1977a) Weidenschilling, S. 1977a, Monthly Notices of the Royal Astronomical Society, 180, 57
  * Weidenschilling (1977b) Weidenschilling, S. 1977b, Astrophys. Space Sci.; (Netherlands), 51:1
  * Youdin & Goodman (2005) Youdin, A. N. & Goodman, J. 2005, The Astrophysical Journal, 620, 459–469


| (J2000) | [$\mu$Jy] | [$\mu$Jy] | [mJy] | [mJy] | [mJy] | [mJy] | [kpc] | [mJy] | [mJy]
SMM J121223.0+273351a (850.1) | 12 12 23.050 +27 33 51.89 | 20 $\pm$ 10 | 30 $\pm$ 10 |  $<$ 0.76 | 72 $\pm$ 21b [4.4] | 54.1 $\pm$ 1.5 [30.5] | 42.3 $\pm$ 1.6c | 5.6 $\pm$ 1.4 | 1.03 $\pm$ 0.08 | 0.49 $\pm$ 0.04
SMM J121220.4+273410a (850.2) | 12 12 20.484 +27 34 10.45 | 23 $\pm$ 10 | 32 $\pm$ 10 |  $<$ 0.76 | 49 $\pm$ 21d [3.5] | 21.2 $\pm$ 1.5 [11.0] | 22.5 $\pm$ 1.6e | 4.6 $\pm$ 2.4 | 0.25 $\pm$ 0.04 | 0.17 $\pm$ 0.03
SMM J121213.4+273517 (850.3) | 12 12 13.43 +27 35 17.0a | … | … | … |  $<$ 60f | 16.5 $\pm$ 2.2 [7.4] | … | … | … | …

Note. — The numbers in square brackets are the signal-to-noise of the flux
density measurements. a Position from SCUBA-2 850 $\mu$m. The positions
derived from the original M15AI29 SCUBA-2 observations of both sources were
offset from those reported here by $\sim$ 4′′ due to the use of M 87 as a
pointing source, where emission from the radio jet significantly perturbs the
position at 850 $\mu$m. This issue was identified from the SMA observations
analyzed in this work and M 87 subsequently removed from the JCMT pointing
catalog. b 850.1 un-deboosted 450-$\mu$m flux density: 93 mJy. c 850.1 peak
SMA flux density: 18.5 $\pm$ 0.8 mJy. d 850.2 un-deboosted 450-$\mu$m flux
density: 73 mJy. e 850.2 peak SMA flux density: 8.9 $\pm$ 0.8 mJy. f
2-$\sigma$ limit.

Table 4: Observed visible and near-infrared photometry

Source | F435W | F606W | F814W | F090W | F150W | F200W | F277W | F356W | F444W | Comments
---|---|---|---|---|---|---|---|---|---|---
850.1 |  $>$ 26.5 | 25.37 0.21 | 24.28 0.17 | 24.07 0.09 | 23.50 0.07 | 23.09 0.06 | 22.49 0.04 | 21.87 0.04 | 21.43 0.03 | 2$\farcs$0⌀ equivalent, including B2
|  $>$ 26.5 |  $>$ 27.2 |  $>$ 26.6 |  $>$ 27.1 |  $>$ 25.8 | 25.07 0.15 | 23.14 0.05 | 22.19 0.04 | 21.66 0.03 | Excluding B2
| 24.46 0.22 | 24.16 0.11 | 22.98 0.10 | 22.86 0.06 | 22.23 0.05 | 21.90 0.05 | 21.48 0.04 | 21.06 0.03 | 20.77 0.03 | Including B1 and B2
850.2 |  $>$ 26.5 | 23.10 0.08 | 21.82 0.05 | 21.75 0.03 | 21.50 0.03 | 21.29 0.03 | 21.03 0.02 | 20.87 0.02 | 20.83 0.02 | 2$\farcs$4⌀
B1 | 25.02 0.27 | 24.98 0.18 | 23.62 0.12 | 23.44 0.07 | 22.92 0.05 | 22.69 0.05 | 22.43 0.04 | 22.22 0.04 | 22.18 0.04 | 0$\farcs$9⌀, knot close to 850.1
B2 |  $>$ 26.5 | 25.37 0.21 | 24.28 0.17 | 24.07 0.09 | 23.50 0.07 | 23.28 0.06 | 23.35 0.07 | 23.36 0.07 | 23.23 0.06 | 0$\farcs$9⌀, knot close to 850.1
A1 |  $>$ 26.5 | 26.04 0.29 | 25.31 0.26 | 25.14 0.15 | 24.63 0.12 | 24.52 0.11 | 24.75 0.13 | 25.24 0.16 | 25.65 0.19 | 0$\farcs$9⌀, galaxy SE of 850.1
A2 |  $>$ 26.5 | 25.84 0.26 | 24.80 0.21 | 24.23 0.10 | 22.02 0.04 | 21.75 0.03 | 21.49 0.03 | 21.31 0.03 | 21.19 0.03 | 0$\farcs$9⌀, galaxy N of 850.1
A3 | 25.13 0.28 | 24.97 0.18 | 24.14 0.15 | 24.17 0.10 | 24.23 0.10 | 24.12 0.10 | 24.27 0.10 | 24.61 0.12 | 24.90 0.14 | 0$\farcs$9⌀, galaxy NW of 850.1
A4 | 23.34 0.12 | 21.52 0.04 | 20.40 0.03 | 20.17 0.02 | 19.33 0.01 | 19.08 0.01 | 19.18 0.01 | 19.67 0.01 | 19.85 0.01 | 1$\farcs$8⌀, bright spiral near 850.1
A5 | 23.87 0.16 | 22.25 0.05 | 21.36 0.04 | 21.17 0.02 | 20.51 0.02 | 20.34 0.02 | 20.52 0.02 | 21.02 0.02 | 21.25 0.03 | 2$\farcs$4⌀, galaxy S of 850.2

Note. — All limits correspond to 3 $\sigma$.

Table 3: ISM properties of sources

Source | $\mu M_{\rm d}$ | $\mu I_{\rm CO(4-3)}$ | $\mu L^{\prime}_{\rm CO(4-3)}$ | $\mu I_{\rm[CI]}$ | $\mu L^{\prime}_{\rm[CI]}$ | $\mu I_{\rm CO(5-4)}$ | $\mu L^{\prime}_{\rm CO(5-4)}$ | $\mu M_{\rm g}^{a}$ | $M_{\rm g}/(M_{\rm g}+M_{\ast})$
---|---|---|---|---|---|---|---|---|---
| [108 M⊙] | [Jy km s-1] | [1010 K km s-1 pc2] | [Jy km s-1] | [1010 K km s-1 pc2] | [Jy km s-1] | [1010 K km s-1 pc2] | [1011 M⊙] |
850.1 | 61${}^{+14}_{-11}$ | 4.57 $\pm$ 0.03 | 20.8 $\pm$ 0.2 | 1.52 $\pm$ 0.02 | 6.1 $\pm$ 0.1 | 3.85 $\pm$ 0.03 | 11.1 $\pm$ 0.1 | 4.2 $\pm$ 0.1 | 0.16 $\pm$ 0.03
850.2 | 26${}^{+3}_{-2}$ | 1.87 $\pm$ 0.04 | 8.5 $\pm$ 0.2 | 0.62 $\pm$ 0.09 | 2.5 $\pm$ 0.4 | 1.76 $\pm$ 0.04 | 5.1 $\pm$ 0.1 | 1.6 $\pm$ 0.1 | 0.62 $\pm$ 0.07

Note. — a Adopting $\alpha_{\rm CO}=$ 1.

Table 5: Observed properties of sources

Source | $z$ | $z_{\rm ph}$ |  $z_{\rm Z20}$a | $\mu M_{\ast}$ |  $\mu$ SFR | $\mu L_{\rm IR}$ | $A_{V}$ | $T_{\rm d}$ | $\mu^{0.5}R^{\rm F444W}_{\rm e}$ | Comments
---|---|---|---|---|---|---|---|---|---|---
|  |  |  | [1011 M⊙] | [103 M⊙ yr-1] | [1012 L⊙] |  | [K] | [kpc] |
850.1b | 4.26 | 4.39${}^{+0.73}_{-1.91}$ | … | 22${}^{+5}_{-2}$ | 3.7${}^{+1.6}_{-1.2}$ | 49${}^{+20}_{-13}$ | 5.0${}^{+0.3}_{-0.2}$ | 41${}^{+8}_{-6}$ | … | Sum of 850.1 & B2 fits, $\mu=$ 4.0${}^{+1.0}_{-2.2}$
| 4.26 | 4.41${}^{+0.77}_{-1.95}$ | … | 22${}^{+5}_{-2}$ | 4.0${}^{+1.6}_{-1.4}$ | 54${}^{+22}_{-15}$ | 5.1${}^{+0.2}_{-0.2}$ | 41${}^{+8}_{-5}$ | 7.5$\pm$0.8 | Excluding B2c
| 4.26 | … | … | 16${}^{+4}_{-1}$ | 4.8${}^{+1.2}_{-0.3}$ | 60${}^{+3}_{-4}$ | 3.6${}^{+1.5}_{-1.2}$d | 43${}^{+2}_{-2}$ | … |  $\Sigma_{\rm pix}$, including B2
850.2 | 4.26 | 4.28${}^{+0.08}_{-0.10}$ | 4.47${}^{+0.07}_{-0.09}$ | 1.0${}^{+0.1}_{-0.1}$ | 1.2${}^{+0.1}_{-0.1}$ | 18${}^{+2}_{-2}$ | 1.1${}^{+0.1}_{-0.1}$ | 38${}^{+2}_{-2}$ | 5.8$\pm$0.6 |  $\mu=$ 5.6${}^{+1.0}_{-3.3}$
| 4.26 | … | … | 2.2${}^{+2.0}_{-0.6}$ | 1.1${}^{+0.9}_{-0.4}$ | 11${}^{+13}_{-5}$ | 0.7${}^{+0.3}_{-0.3}$ | 44${}^{+2}_{-2}$ | … | $\Sigma_{\rm pix}$
B1 | 1.1 | 1.12${}^{+0.11}_{-0.12}$ | 0.91${}^{+0.47}_{-0.13}$ | 0.1${}^{+0.1}_{-0.0}$ | 0.01${}^{+0.01}_{-0.01}$ | 0.1${}^{+0.1}_{-0.1}$ | 1.06${}^{+0.50}_{-0.47}$ | … | … |
B2 | 4.26 | 4.17${}^{+0.16}_{-3.28}$ | 0.33${}^{+0.20}_{-0.10}$ | 0.04${}^{+0.06}_{-0.02}$ | 0.04${}^{+0.02}_{-0.03}$ | 0.53${}^{+0.78}_{-0.50}$ | 0.6${}^{+0.9}_{-0.3}$ | … | … |  $\mu\sim$ 4
A1 | 0.35 | 0.31${}^{+0.16}_{-0.04}$ | 0.53${}^{+0.72}_{-0.38}$ | … | … | … | … | … | … | Cluster member?
A2 | 1.8 | 1.77${}^{+0.28}_{-0.25}$ | 0.55${}^{+0.54}_{-0.25}$ | … | … | … | … | … | … |
A3 | 0.8 | 0.77${}^{+0.11}_{-0.12}$ | 0.54${}^{+1.57}_{-0.39}$ | … | … | … | … | … | … |
A4 | 0.35 | 0.46${}^{+0.08}_{-0.09}$ | 0.35${}^{+0.14}_{-0.05}$ | … | … | … | … | … | … | Cluster member?
A5 | 0.35 | 0.35${}^{+0.06}_{-0.03}$ | 0.29${}^{+0.03}_{-0.06}$ | … | … | … | … | … | … | Cluster member?

Note. — a BPZ-derived photometric redshifts from Zitrin et al. (2020) using
only HST photometry. b $z_{\rm ph}=$ 2.63${}^{+0.46}_{-0.23}$ when B1 and B2
are included in the photometric aperture. c A single magphys fit at $z$ $=$
4.26 to the “unresolved” photometry of 850.1 including B2 increases $\chi^{2}$
from 1.7 to 27.1 and the fit parameters become: $\mu M_{\ast}$ $=$
0.9${}^{+0.1}_{-0.1}$ $\times$ 1011 M⊙, $\mu$ SFR $=$
0.9${}_{-0.1}^{+0.1}\times$ 103 M⊙ yr-1, $\mu M_{\rm d}$ $=$
99${}_{-20}^{+20}\times$ 108 M⊙, $A_{V}$ $=$ 2.7${}_{-0.1}^{+0.1}$, $T_{\rm
d}$ $=$ 30${}^{+3}_{-1}$ K with an age for the best fit star-formation history
of 50${}^{+1}_{-1}$ Myr. d F444W-weighted mean $A_{V}$.

Table 6: Intrinsic properties of sources

Propertya | 850.1b | 850.2 | Units
---|---|---|---
$z$ | 4.2599 $\pm$ 0.0001 | 4.2593 $\pm$ 0.0003 |
$\mu$ | 4.0${}^{+1.0}_{-2.2}$ | 5.6${}^{+1.0}_{-3.3}$ |
$S_{\rm 850\mu m}$ | 13.5${}^{+0.4}_{-0.4}$ ${}^{+5.8}_{-2.7}$ | 3.8${}^{+0.3}_{-0.3}$ ${}^{+5.4}_{-0.6}$ | mJy
$M_{\ast}$ | 5.5${}^{+1.3}_{-0.5}$ ${}^{+2.4}_{-1.1}$ | 0.18${}^{+0.02}_{-0.02}$ ${}^{+0.25}_{-0.03}$ | 1011 M⊙
$M_{\rm d}$ | 15${}^{+3}_{-3}$ ${}^{+7}_{-3}$ | 4.6${}^{+0.5}_{-0.4}$ ${}^{+6.7}_{-0.7}$ | 108 M⊙
$M_{\rm g}$ | 1.05${}^{+0.03}_{-0.03}$ ${}^{+0.45}_{-0.21}$ | 0.29${}^{+0.02}_{-0.02}$ ${}^{+0.41}_{-0.05}$ | 1011 M⊙
$M_{\ast}+M_{\rm g}$ | 6.6${}^{+1.3}_{-0.5}$ ${}^{+2.4}_{-1.1}$ | 0.5${}^{+0.03}_{-0.03}$ ${}^{+0.5}_{-0.5}$ | 1011 M⊙
SFR | 900${}^{+400}_{-300}$ ${}^{+400}_{-150}$ | 210${}^{+20}_{-20}$ ${}^{+300}_{-30}$ | M⊙ yr-1
$R^{\rm F444W}_{\rm e}$ | 3.8${}^{+0.4}_{-0.4}$ ${}^{+0.7}_{-0.4}$ | 2.5${}^{+0.3}_{-0.3}$ ${}^{+1.3}_{-0.2}$ | kpc
$R^{\rm dust}_{\rm e}$ | 2.8${}^{+0.7}_{-0.7}$ ${}^{+0.5}_{-0.3}$ | 1.9${}^{+1.0}_{-1.0}$ ${}^{+1.1}_{-0.1}$ | kpc
$R^{\rm mass}_{\rm e}$ | 1.0$\pm$0.2 ${}^{+0.2}_{-0.1}$ | 2.2$\pm$0.4 ${}^{+1.1}_{-0.2}$ | kpc

Note. — a Dual error bars give first the random and then the systematic
uncertainties due to the lensing magnification. b Sum of the parameter
estimates for the fit to 850.1 and the fit to B2.