© K. Silsbee and R. R. Rafikov 2021

1 Introduction

Planets have been discovered in a wide variety of stellar systems, including stellar binaries and systems of higher multiplicity (Marzari & Thebault 2019). Dynamically active environments typical of multiple stellar systems are believed to be hostile for planetary assembly, positioning such planet-hosting systems as important test beds of planet formation theory.

In particular, tight eccentric binary systems that host planets orbiting one of the stellar companions (so-called S-type systems in the classification of Dvorak 1982), such as γ Cephei (Hatzes et al. 2003), present important challenges to theories of planet formation via core accretion, which involves a steady agglomeration of planetary cores through mutual collisions of numerous small planetesimals. Indeed, gravitational perturbations from the eccentric companion in such binaries with semimajor axes ≲20–30 AU are expected to drive planetesimal eccentricities to high values. This naturally results in the destruction, rather than the growth, of planetesimals in mutual collisions. A simple calculation (Heppenheimer 1978), including only the dominant secular effects of the stellar companion on planetesimal orbits, yields planetesimal collision velocities of a few km s⁻¹ at the current location of the planet in the γ Cephei system (around 2 AU from the stellar primary). This is enough to destroy even planetesimals of hundreds of kilometers in size in catastrophic collisions, presenting an important barrier to planet formation in binaries.

A number of ideas have been advanced over the years to overcome this “fragmentation barrier.” In particular, Marzari & Scholl (2000) noticed that aerodynamic drag due to the gaseous component of the protoplanetary disk in which planetesimals move eventually apsidally aligns their orbits. While their orbits would still have eccentricities similar to those in the Heppenheimer (1978) model, their relative eccentricities (which determine the magnitude of the mutual collision velocities) would be dramatically reduced by apsidal alignment. However, it was later recognized by Thébault et al. (2006) that, even with apsidal alignment, only planetesimals of similar size would have small relative velocities as the strength of the frictional coupling to the gas disk is directly determined by the planetesimal size. As a result, planetesimals with different sizes are not well aligned with one another and therefore experience high collision velocities (Thébault et al. 2008, 2009; Thebault 2011).

On the other hand, Rafikov (2013a,b) noted that the gas disk should couple to planetesimals not only through gas drag, but also gravitationally. Neglecting the gas drag in the first instance, he showed that a massive axisymmetric disk should induce a rapid apsidal precession of planetesimal orbits, severely suppressing the eccentricity excitation due to the stellar companion. However, it is not at all clear that the protoplanetary disk orbiting one of the components of the binary should be axisymmetric, given the strong gravitational perturbation due to the external stellar companion (Paardekooper et al. 2008; Marzari et al. 2009; Regály et al. 2011). For this reason, Silsbee & Rafikov (2015b) generalized the work of Rafikov (2013a) and calculated the gravitational effect of a massive eccentric disk on the orbits of planetesimals moving inside such a disk (again, neglecting gas drag), finding that collision outcomes depend heavily on the details of the disk structure.

Subsequently, Rafikov & Silsbee (2015a) extended the model of Silsbee & Rafikov (2015b) by fully incorporating the effects of gas drag on planetesimal motion in eccentric protoplanetary disks and derived collision velocities of planetesimals as a function of their sizes. These results on planetesimal dynamics, coupled with simple prescriptions for collisional outcomes based on Stewart & Leinhardt (2009), were then used in Rafikov & Silsbee (2015b) to determine the conditions under which planetesimals could grow to sizes large enough to withstand collisional fragmentation and erosion (thus providing a pathway to planetary core assembly bycoagulation).

Despite this work, some uncertainty remains regarding the collisional environments that would allow planetesimal growth. In particular, Rafikov & Silsbee (2015b) used simple heuristics when deciding whether planetesimal growth is possible, such as assuming it to occur if there are no catastrophically disruptive collisions. However, this condition is neither necessary nor sufficient forunimpeded growth, since evolution of the planetesimal mass spectrum is affected both by coagulation at low relative velocities and by erosion and catastrophic disruption at high relative velocities. The size-dependent alignment of planetesimal orbits (Thébault et al. 2008; Rafikov & Silsbee 2015a) gives rise to a dynamical environment with a mixture of low- and high-velocity collisions, making it challenging, even qualitatively, to determine the evolution of the size distribution.

The problem is further complicated by the radial inspiral of planetesimals due to interaction with the gas. The gas disk is slightly pressure-supported and orbits ~ 30 m s⁻¹ slower than the Keplerian speed (Weidenschilling 1977). Large bodies with stopping times longer than an orbital time move at the Keplerian speed and therefore experience a headwind, which causes them to spiral into the star. The inspiral is more rapid in the case of a body having a nonzero forced eccentricity relative to the gas (Adachi et al. 1976). It was suggested in Rafikov & Silsbee (2015b) that the dependence of the inspiral rate on planetesimal eccentricity could concentrate solid bodies in the regions of the disk with low relative particle-gas eccentricity where radial migration slows down, thus accelerating coagulation. Inspiral may also help to alleviate the detrimental effect of the erosion of growing planetesimals by small bodies by rapidly removing such bodies from the system.

In this paper, we present the first realistic global model of planetesimal growth in S-type stellar binaries, which explicitly includes the aforementioned physical ingredients missing in previous studies. A specific question we address is whether the fragmentation barrier can somehow be overcome such that planetesimals can grow from some small size (that we determine as part of our analysis) to sizes of several hundred kilometers, at which point fragmentation is no longer a threat to their growth. Our modeling framework employs a multi-annulus coagulation-fragmentation code for accurately following the evolution of the planetesimal size spectrum at every radius in the disk. In addition to that, our model directly follows the exchange of mass in solid objects between the different radial locations in the disk. These components of the model use the collisional velocities and radial inspiral rates calculated in Rafikov & Silsbee (2015a,b), which fully account for the secular gravitational effects of both the stellar companion and the massive eccentric disk in which planetesimals orbit. Outcomes of planetesimal collisions are determined using the prescriptions derived in Stewart & Leinhardt (2009). With this physical model, we are able to determine whether the coagulation process as a whole is successful for a given environment (i.e., a disk model), rather than determining if it is successful based on the outcomes of individual collisions. Armed with this powerful tool, we then perform an extensive exploration of the parameter space of the disk models to determine the conditions under which planet formation is possible in compact eccentric binaries such as γ Cephei or α Centauri.

This paper is organized as follows. After describing our model setup in Sect. 2, we provide an overview of the important aspects of planetesimal dynamics in S-type binaries in Sect. 3. We briefly outline the functions of the coagulation code in Sect. 4; a detailed description of the code and its tests can be found in Appendices A and B, respectively. In Sect. 5, we describe simulations for our fiducial disk model, both in single- (Sect. 5.1) and multi-annulus (Sect. 5.2) setups. We then present the results of the model parameter space exploration in Sect. 6. In Sect. 7, we provide a discussion of our findings with implications for planet formation in binaries (Sect. 7.1), planetesimal formation (Sect. 7.2), and hydrodynamic simulations of protoplanetary disks in S-type binaries (Sect. 7.3), among other things. We summarize our main conclusions in Sect. 8.

2 Model setup

To illustrate our calculations, we consider a model S-type binary with parameters similar to that of the γ Cephei system (Hatzes et al. 2003; Chauvin et al. 2011). This binary consists of M_p = 1.6 M_⊙ (primary) and M_s = 0.4 M_⊙ (secondary) stars in an orbit with a semimajor axis a_b = 20 AU and eccentricity e_b = 0.4.

A planet with mass M_pl = 1.6 M_J orbits the primary star of the γ Cephei system with the semimajor axis a_pl = 2 AU and eccentricity e_pl = 0.12, which we use to motivate our subsequent calculations. There are a number of other binaries with a_b ≲ 20 AU that have similar characteristics (Chauvin et al. 2011; Marzari & Thebault 2019).

We assume the primary component of the binary to host a massive protoplanetary disk, coplanar with the binary orbit, with properties similar to the disk described in Silsbee & Rafikov (2015b): it is an eccentric disk with a power law radial dependence of the surface density Σ(a) and eccentricity e(a) on the semimajor axis a of the eccentric fluid trajectories. The stated dependence on a strictly applies to the surface density and eccentricity of the gas streamlines at their periastra: $$\Sigma (a)={\Sigma }_0\left(\frac{a_0}{a}\right),e_{\text{g}}(a)=e_0{\left(\frac{a_0}{a}\right)}^{-1},a<a_{\text{out}},$$(1)

where a₀ is a reference semimajor axis. We assume the disk to be sharply truncated by the torque due to the secondary star at the outer boundary a_out = 5 AU (i.e., Σ(a) = 0 for a > a_out). The apsidal angle of the gas streamlines in the disk with respect to the binary orbit is denoted ϖ_d and, for simplicity, is assumed to be independent of a, that is, the disk as a whole is apsidally aligned. There is no particular reason for choosing e_g (a) in the form (1), except that it scales as the forced eccentricity due to the secondary. Development of a more accurate, better motivated model for the e_g(a) behavior is beyond the scope of this study. Throughout this paper, we assume a solid-to-gas ratio of 0.01, and the disk to be slightly flared with the h∕r profile given by Eq. (14) from Rafikov & Silsbee (2015b) and varying between 0.03 and 0.05.

As a consequence of the generally nonzero disk eccentricity, the surface density distribution in the disk ends up being non-axisymmetric, as described in Ogilvie (2001), Statler (2001), and Silsbee & Rafikov (2015b). Since in this work we fully account for the gravitational effect of the disk on the motion of the planetesimals embedded in the disk, the non-axisymmetry of the disk generally leads to nonzero torque affecting planetesimal dynamics, in addition to the torque exerted by the secondary star. We cover the roles played by these different dynamical agents next.

3 Physical ingredients of the model

Numerical calculations self-consistently following growth of planetesimals in protoplanetary disks from small sizes, and often all theway into the planetary regime, have been routinely used to understand planet formation around single stars. A standard approach is to follow the evolution of the size (mass) spectrum of colliding objects by solving the so-called Smoluchowski (Smoluchowski 1916) or coagulation equation (Safronov 1972; Wetherill & Stewart 1989; Kenyon & Luu 1998), often with account being given to the possibility of collisional fragmentation (Wetherill & Stewart 1993; Kenyon & Luu 1999). Global coagulation simulations, incorporating radial drift of mass in the disk driven by gas drag (Spaute et al. 1991; Inaba et al. 2001; Kenyon & Bromley 2002; Kobayashi et al. 2010), represent more sophisticated versions of such calculations. Some studies have also included direct N-body integration of a modest number (~ 10³) of planetary embryos formed as a result of planetesimal coagulation (Bromley & Kenyon 2011), and accounted for the most recent results on the outcomes of planetesimal collisions (Stewart & Leinhardt 2009). However, probably the most important ingredient of such coagulation-fragmentation calculations is the proper description of planetesimal dynamics (Stewart & Wetherill 1988; Stewart & Ida 2000; Rafikov 2003a,b,c, 2004), which determines a particular mode of planetary growth (Safronov 1972; Ida & Makino 1993; Kokubo & Ida 1998; Rafikov 2003d).

While these developments significantly advanced our ideas about the origin of planets around single stars, a proper understanding of planet formation in stellar binaries requires a number of physical ingredients that are unique to these systems to be accounted for. These include: (i) excitation of planetesimal eccentricities by gravitational perturbations due to both the stellar companion and the massive, non-axisymmetric protoplanetary disk (Sect. 3.1); (ii) a unique size-dependent dynamical state (Sect. 3.2) in which planetesimals are a result of coupling between the aforementioned excitation and gas drag, strongly varying across the disk and featuring resonant locations; (iii) a nonstandard distribution of planetesimal collisional velocities (Sect. 3.3), different from that in disks around single stars; (iv) an important role played by the destructive collisional outcomes (erosion and catastrophic disruption) in determining the evolution of the planetesimal size distribution (Sect. 3.4); and (v) nonuniform radial drift of planetesimals across the disk (Sect. 3.5) strongly affected by their nontrivial dynamics, which may lead to their trapping at certain locations.

In the rest of this section, we expand on how we calculate these effects in our models.

3.1 Gravitational perturbations in S-type binaries

The motion of planetesimals around the primary in S-type binaries is affected by the gravity of both the secondary and the protoplanetary disk. Since the seminal work of Heppenheimer (1978), the direct effect of the secondary’s gravity has been accounted for using the secular approximation, which averages the (time-dependent) potential of the secondary over its orbit. However, since the gaseous protoplanetary disks in S-type binaries are themselves prone to developing eccentricity under the perturbation due to the secondary (Paardekooper et al. 2008; Regály et al. 2011), one needs to also account for the gravitational potential of a non-axisymmetric, eccentric disk. Silsbee & Rafikov (2015b) computed the potential due to an apsidally aligned disk with the surface density and eccentricity profiles given by Eq. (1); their calculation was subsequently generalized for a disk with arbitrary profiles of Σ(a) and e_g (a) by Davydenkova & Rafikov (2018).

Including disk potential, the secular disturbing function characterizing the perturbation of planetesimal motion by the external potential becomes (Silsbee & Rafikov 2015b) $$R=a^{2}n\left[\frac{1}{2}A{e}^2+B_{\text{b}}e\cos \varpi +B_{\text{d}}e\cos (\varpi -{\varpi }_{\text{d}})\right],$$(2)

where a, n, and ϖ are the semimajor axis, mean motion, and apsidal angle of planetesimal orbit (both ϖ and the disk apsidal angle, ϖ_d, are measured relative to the apsidal line of the binary), B_b and B_d describe the torques due to the non-axisymmetric components of the potentials of the secondary and the disk, and A = A_b + A_d is the combined free precession rate due to the secondary (A_b) and the disk (A_d). The explicit expressions for A_d, B_d, A_b and B_b can be found in Silsbee & Rafikov (2015b, see their Eqs. (5)–(6) and (9)–(10), respectively). We note that while A_d depends only on the Σ(a) profile, B_d is determined by both Σ(a) and e_g(a).

Of particular importance for planet formation in S-type binaries is the fact that the secondary companion always drives prograde apsidal precession of the planetesimal orbit (A_b > 0), while the precession due to disk gravity is usually retrograde (A_d < 0). Since A_b and A_d vary with a in different ways, one often finds that A = A_b + A_d → 0 at a particular semimajor axis in the disk, if the disk is sufficiently massive (≳ 10⁻³ M_⊙, Rafikov & Silsbee 2015b) but not too massive (see below). At this radial location a secular resonance emerges, at which the eccentricities of planetesimals are driven to very high values (in the absence of damping effects). Such resonances due to disk gravity werepreviously uncovered in the context of planet formation in binaries in Rafikov (2013a), Meschiari (2014), Silsbee & Rafikov (2015b,a), and Rafikov & Silsbee (2015a) and for massive disks with planets in Heppenheimer (1980), Ward (1981), and Sefilian et al. (2021).

In addition to the secular perturbations we have considered, the companion star also produces short-period perturbations which act on the orbital period of the planetesimals (Murray & Dermott 1999). The magnitude of the short-period perturbations is quite small for bodies which are close together. If we consider two bodies with relative eccentricity, e_r, which will collide in the next orbit, then their distance is of order e_ra_p. Their relative acceleration due to the tidal perturbation of the secondary is $G{M}_{\text{s}}{e}_{\text{r}}{a}_{\text{p}}/R_{\text{s}}^3$. If this acceleration acts over a time $~n_p^{-1}$ it results in a relative velocity $e_{\text{r}}{v}_{\text{K}}{(a_{\text{p}}/R_{\text{s}})}^3$, which is suppressed by a factor ${(a_{\text{p}}/R_{\text{s}})}^3$ relative to their typical relative velocity e_rv_K in the absence of such short-period perturbations. Over longer timescales such short-period perturbations average out to zero.

Because the disk we consider is small compared to the binary orbit, the resonant perturbations are also not important. Indeed, the only resonances present in the disk are of very high order, and therefore very narrow and unlikely to significantly affect the dynamics.

The torque due to the secondary is expected to sharply (although not as discontinuously as we assume in Eq. (1)) truncate disk density at the semimajor axis a_out. We take this truncation into account when computing disk gravity. The disk gravity-related coefficients A_d and B_d diverge logarithmically when a sharp disk edge is approached (Silsbee & Rafikov 2015b; Davydenkova & Rafikov 2018). For this reason, we leave a buffer of 1 AU between the outer edge of the disk and the region in which we model planetesimal growth. The inner boundary of the disk is ignored in this calculation: As shown in Fig. 10 of Silsbee & Rafikov (2015b), provided that the inner boundary is located at no more than half the semimajor axis of interest, it makes no significant contribution to the disk disturbing function $R$.

3.2 Overview of planetesimal dynamics in S-type binaries

The model for planetesimal dynamics employed in this work builds upon the understanding of the planetesimal disturbing function described in Sect. 3.1, while additionally incorporating the dissipative effect of gas drag caused by the relative motion of planetesimals with respect to the noncircular streamlines of the eccentric protoplanetary disk. It has been developed in Rafikov & Silsbee (2015a) for the realistic case of the quadratic gas drag (Adachi et al. 1976), and here we briefly summarize its main features relevant for the present work.

Rafikov & Silsbee (2015a) found that planetesimal dynamics admits two distinct regimes, depending on the degree of coupling between the planetesimal and the gas. Strongly coupled (small) planetesimals follow orbits nearly aligned with the local eccentric gas streamline. Weakly coupled (large) planetesimals have orbits nearly aligned with the forced eccentricity determined by the gravitational perturbations due to both the companion and the disk. Separating these two regimes is the characteristic planetesimal size d_c that is defined as the planetesimal size at which the damping timescale¹ τ_d of the planetesimal eccentricity (relative to the eccentricity vector of the local gas streamline e_g = e_g(cosϖ_d, sinϖ_d)) is of order A⁻¹, where A is the precession rate of free eccentricity in the combined potential of the disk and perturbing star (see Eq. (2)). More specifically, Rafikov & Silsbee (2015a) show that for a planetesimal of radius d_p $$\left|A{\tau }_{\text{d}}\right|=\frac{d_{\text{p}}}{d_{\text{c}}}{\left[\frac{1}{2}+\sqrt{\frac{1}{4}+{\left(\frac{d_{\text{c}}}{d_{\text{p}}}\right)}^2}\right]}^{1/2}$$(3)

(see their Eq. (32)), so that |Aτ_d|≈ 1.27 for d_p = d_c.

Rafikov & Silsbee (2015a) also show that the relative planetesimal-gas eccentricity e_r = |e −e_g| is given by (see their Eq. (28)) $$e_{\text{r}}=\frac{v_{\text{c}}}{v_{\text{K}}}\frac{A{\tau }_{\text{d}}}{\sqrt{1+{(A{\tau }_{\text{d}})}^2}},$$(4)

where v_c = e_cv_K is the characteristic planetesimal-gas velocity for planetesimals with d_p ≥ d_c: e_c is the magnitude of the relative eccentricity between the gas streamlines and the forced eccentricity from the gravitational potential. An explicit expression for e_c in terms of the disk and binary parameters only is provided by Eq. (29) of Rafikov & Silsbee (2015a). Knowing e_c one can then directly compute d_c through Eq. (31) of Rafikov & Silsbee (2015a).

Both d_c and v_c (or e_c) are the key metrics of planetesimal dynamics in S-type systems and will frequently appear in this work. Using these two parameters, strong coupling corresponds to small bodies, d_p ≪ d_c, when $\left|A{\tau }_{\text{d}}\right|\approx {(d_{\text{p}}/d_{\text{c}})}^{1/2}\ll 1$ and $e_{\text{r}}\approx e_{\text{c}}{(d_{\text{p}}/d_{\text{c}})}^{1/2}\ll e_{\text{c}}$. Weak coupling is realized in the opposite limit of d_p ≪ d_c, when |Aτ_d|≈ d_p∕d_c ≫ 1 and e_r ≈ e_c. This suggests another intuitive definition of v_c (or e_c) as the characteristic relative velocity (or eccentricity) between the weakly and strongly coupled (or large and small) objects.

The parameters v_c and d_c vary significantly throughout the disk, and with different disk models, which is illustrated in Fig. 1, where we vary (one at a time) the values of the disk eccentricity e₀, surface density normalization Σ₀ and disk orientation ϖ_d. One can see that the behavior of v_c and d_c directly reflects many peculiar features of planetesimal dynamics that emerge when both secondary and disk gravity are important, as described in Sect. 3.1. In particular, (i) d_c and v_c diverge at 2.5 AU in a model with Σ₀ = 500 g cm⁻² because of a secular resonance, which is not present in higher-mass models as it gets pushed out of the disk (i.e., disk gravity dominates planetesimal dynamics all the way out to a_out); (ii) many models, in which the disk is apsidally aligned with the binary, feature a dynamically quiet location where d_c, v_c → 0 (with the semimajor axis varying with Σ₀ and e₀); (iii) dynamically quiet locations do not appear for all values of e₀ and Σ₀ or for apsidally misaligned disks (e.g., for ϖ_d = 90°). These observations will greatly help us in interpreting the results of our detailed calculations in Sects. 5 and 6.

The concept of a dynamically quiet location in the disk will be very important for understanding the results of this work. Rafikov & Silsbee (2015a) have demonstrated (see their Eq. (50)) that in an apsidally aligned disk, in the presence of gas drag, d_c, e_c → 0 at the semimajor axis a_dq where $$A(a_{\text{dq}})e_{\text{g}}(a_{\text{dq}})+B_{\text{d}}(a_{\text{dq}})+B_{\text{b}}(a_{\text{dq}})=0.$$(5)

At this location, the forced eccentricity is equal to the local gas eccentricity, so weakly and strongly coupled planetesimals are driven to identical orbits.

It should also be noted that the description of planetesimal dynamics provided in Rafikov & Silsbee (2015a) and used in this work assumes that planetesimal eccentricities are at their steady state values, given by Eqs. (22)–(27) in Rafikov & Silsbee (2015a), to which they converge on a timescale ~ τ_d due to gas drag. This is a reasonable assumption since for small objects τ_d is short, ≲ 10³ yr for d_p = 1 km objects, while for bigger bodies (for which τ_d ∝ d_p is longer) the timescale for collisions with comparable objects (capable of significantly perturbing eccentricity) is long.

Fig. 1 Radial profiles of characteristic velocity vc (a) and planetesimal size dc (b) for different disk models, in which we vary disk eccentricity e0, surface density normalization Σ0, and apsidal angle ϖd (one parameter at a time relative to the fiducial model described in Sect. 5 and shown with the thick black curve), as indicated in the legend. See the text for more details on various features of these curves.

3.3 Collision velocities

An important feature of planetesimal dynamics in S-type binaries is that not only the planetesimal eccentricity but also the apsidal angle takes on a unique, size-dependent value at each semimajor axis as a result of a competition between gravitational perturbations and gas drag. This is different from disks around single stars, in which a balance between gas drag and dynamical excitation due to numerous embedded objects – other planetesimals and planetary embryos – results in a size-dependent planetesimal eccentricity distribution, whereas the apsidal angles of the planetesimals are randomly distributed.

This difference becomes important when calculating the distribution of collision velocities of planetesimals, a procedure which is described in detail in Sect. A.3. Around single stars, the relative velocities of planetesimals follow a 3D-Gaussian (or Schwarzschild) distribution (Stewart & Ida 2000; Rafikov 2003c,a). But in S-type binaries, neglecting dynamical excitation by embedded objects, Rafikov & Silsbee (2015a) found a very different shape for the distribution of relative velocities, given by Eq. (A.12), with a velocity scale set by the relative eccentricity of the two colliding objects (see Eq. (62) of that paper).

In this study, we also account for random planetesimal motions due to dynamical excitation by embedded objects, in addition to those arising from the secondary and disk eccentricity. In particular, we assume that the two components of each planetesimal’s inclination vector are drawn independently from a Gaussian distribution with standard deviation σ_i, which, for simplicity, is assumed to be independent of the size of the planetesimal. We then assume each planetesimal’s eccentricity vector to consist of the equilibrium component calculated in Rafikov & Silsbee (2015a) as described in Sect. 3.2, plus a random component. We take this random component of the eccentricity vector to have standard deviation σ_e which is twice σ_i (Stewart & Ida 2000), but allow horizontal and vertical motions to be independent of one another. The degree of random motion can therefore be fully parameterized by σ_i.

If random motions result from viscous stirring, one would typically expect them to be on the order of the escape velocity from the planetesimals. However, in the present perturbed system, because the relative motions between planetesimals are much larger than their random motions, the excitation from two-body encounters is lowered and the damping due to gas drag is enhanced. On the other hand, if present, disk turbulence could excite larger random motions. A study of the random planetesimal motions in a binary system would be subject to many uncertainties, and we do not attempt it here. We note that our typical random motions correspond to the escape velocity from roughly 10 km planetesimals.

Varying σ_i has two effects. Let e_ij be the relative eccentricity between two bodies in the absence of any random motions. For collisions between bodies where e_ij ≫ σ_i, varying σ_i has little effect on the collision speed, but the frequency of such collisions varies inversely with σ_i, as σ_i sets the thickness of the planetesimal disk and therefore the planetesimal number density. In the limit that e_ij ≪ σ_i, the collision velocity is proportionalto σ_i and the collision rate becomes independent of σ_i provided that σ_iv_K ≫ v_esc (where $v_{\text{esc}}^2=2G(m_1+m_2)/(d_1+d_2)$ for colliding objects with sizes d₁, d₂ and masses m₁, m₂), that is, gravitational focusing is negligible; the collision rate scales as ${\sigma }_i^{-2}$ in the opposite limit σ_iv_K ≪ v_esc when gravitational focusing enhances the collision rate. The effect of random velocities on the collision velocities and rates and their implementation in our numerical framework are discussed in detail in Sects. A.3 and A.4.

Figure 2 illustrates the collision velocities v_coll between planetesimals as a function of their sizes. These are given by Eq. (A.16), with λ and e_rand assigned their mean values of 0.81 and $2\sqrt{\pi }{\sigma }_i$, respectively.This figure is made for the fiducial system described in Sect. 5 with rather small σ_i = 10⁻⁴, so that relative velocities are typically dominated by the secondary and disk forcing (i.e., by e_ij). The three panels differ in the assumed location a_p of the colliding planetesimals in the disk (which determines their v_c and d_c), as labeled on the figure. As mentioned in Sect. 3.2, v_coll ~ v_c for pairs of collision partners in which oneobject has size d₁ ≪ d_c and another has size d₂ ≫ d_c. If the collision partners have similar size, or both are much larger (or both much smaller) than d_c, then the collision velocities are determined predominantly by random motions and become independent of the sizes of the collision partners (since we take σ_i to be size independent in this work).

As mentioned above, Fig. 2 would look very different if it were drawn assuming planetesimal dynamics typical for disks around single stars (i.e., without apsidal alignment): in that case the “valley" at d₁ ≈ d₂ would disappear and v_coll would monotonically increase along this line (since smaller planetesimals have lower random velocities as a result of gas drag, Rafikov 2003a, 2004). These differences clearly demonstrate the importance of apsidal alignment naturally emerging in disks around binaries for setting the pattern of planetesimal collision velocities.

3.4 Collision outcomes

Armed with the understanding of the distribution of relative velocities of planetesimals (Sect. 3.3), we determine the outcomeof their physical collision using the recipes provided in the work of Stewart & Leinhardt (2009). This study assumes, quite generally, that as a result of a collision one is left with a large remnant body and a spectrum of small fragments. The size of the largest remnant is the main ingredient in their description of the collision outcome (see our Eq. (A.4)). We provide the details of our implementation of this prescription in Sect. A.2, but will briefly highlight the differences with the single-star case here.

Following Rafikov & Silsbee (2015b), we divide the collision outcomes into three classes. We say that a collision leads to catastrophicdisruption if the largest remnant contains less than half the combined mass of the two incoming planetesimals, erosion if the largest remnant is smaller than the larger of the two incoming bodies, and growth if the largest remnant is larger than either of the incoming bodies. We can then use the recipe of Stewart & Leinhardt (2009) to delineate in Fig. 2 the domains in d₁ − d₂ phase space,corresponding to each to type of the outcome: regions of catastrophic disruption lie within the solid white lines, regions of erosive collisions are outlined in dotted white, and regions leading to growth lie outside both white lines. We use the “strong rock” planetesimal composition of Stewart & Leinhardt (2009) in this calculation. One can see that, unless v_c is very high, catastrophicdisruption occurs only between collision partners with similar (but not equal) sizes. It is mainly relevant for objects with sizes around d_c, as that is where even objects of comparable size experience high-velocity collisions (although the regions of catastrophic disruption are not exactly centered on d_c because the strength of planetesimals is size-dependent, with a minimum around 0.1 km (Stewart & Leinhardt 2009). In contrast to catastrophic disruption, erosion is possible even when the sizes of the colliding bodies are very unequal.

Due to the complexity of secular dynamics in binaries outlined in Sect. 3.2, there is a significant variation in the collisional environment across the disk, as different panels of Fig. 2 demonstrate. The regions corresponding to erosion and catastrophic disruption grow larger as the characteristic velocity v_c increases. The weakest planetesimals with sizes around 0.1 km may be destroyed even in the absence of any secularly excited eccentricities, in collisions arising just from random motions, which is almost what is happening in Fig. 2a. But it is still clear from panels (b) and (c) of that figure that secular pumping of planetesimal eccentricities in binaries endows these systems with much harsher collisional environments than in disks around single stars.

The heuristic picture of collision outcomes based on the relative velocity maps has been used in Rafikov & Silsbee (2015a) and Silsbee & Rafikov (2015a) to understand the possibility of planetary growth, which was assumed to take place when neither catastrophic disruption nor substantial erosion were taking place at certain locations in the disk. Obviously, such a simplistic approach cannot provide a perfect characterization of the conditions necessary for sustained planetesimal growth. Indeed, persistent planetesimal accretion may be possible even if some collisions are destructive; however, one does not know a priori how harsh of a collisional environment can be tolerated. Thébault et al. (2008, 2009), Thebault (2011), and Rafikov & Silsbee (2015b) all had to make certain assumptions when drawing their conclusions on the likelihood of planet formation in binaries (see Sect. 7.4 for an assessment of their realism). For this reason, in this work we follow the evolution of the planetesimal size distribution in full, using detailed coagulation-fragmentation calculations with realistic physical inputs as described in Sect. 4 and Appendix A.

Fig. 2 Collision velocities between two planetesimals as a function of their sizes d1 and d2 in a disk with the fiducial set of parameters given in Sect. 5 and σi = 10−4. Each panel corresponds to a different location in the disk, which leads to the varying values of vc and dc labeled on the panels. Regions of erosive collisions are outlined in dotted white lines, and regions of catastrophic disruption are enclosed by solid white lines. The horizontal and vertical yellow lines correspond to the initial planetesimal size dinit = 4.3 km used in the simulation described in Sect. 5.1. The point (dc, dc) is shown as a red dot.

3.5 Radial inspiral

In this work, we also explicitly include the radial inspiral of planetesimals, which we expect to impact planetesimal growth in two different ways. First, inspiral drains the disk of solid material. If this happens more rapidly than coagulation, then large bodies will not form due to lack of solid material needed for their growth. Second, since small bodies inspiral more rapidly than large ones, one might expect that, once small bodies are flushed out, the effect of erosive collisions would be reduced, thus favoring the growth of large planetesimals.

Using the results of Adachi et al. (1976), Rafikov & Silsbee (2015b) found that the inspiral rate ȧ_p is given by (their Eq. (11)): $${\dot{a}}_{\text{p}}=-\pi \frac{a_{\text{p}}}{e_{\text{r}}{\tau }_{\text{d}}}{\left(\frac{5}{8}{e}_{\text{r}}^2+{\eta }^2\right)}^{1/2}\left[\left(\frac{\alpha }{4}+\frac{5}{16}\right)e_{\text{r}}^2+\eta \right],$$(6)

where e_r and τ_d are given by Eqs. (4) and (3), respectively, $\eta ~{(h/a_{\text{p}})}^2~0.003$ (h is the disk scale height) is a measure of the particle-gas velocity differential due to the pressure support in the gas disk, and α ~ 1 is a constant that depends on the disk model as described in Rafikov & Silsbee (2015b). From Eqs. (3) and (4), one can show (Rafikov & Silsbee 2015b) that e_r τ_d ∝ d_p. As a result, one finds that ${\dot{a}}_{\text{p}}\propto d_{\text{p}}^{-1}$ if either (i) d_p ≫ d_c, when e_r becomes constant, or (ii) e_r ≪ η, so that the last term in brackets in Eq. (6) dominates.

Figure 3 shows the speed of radial drift for different planetesimal sizes and locations in the disk in the γ Cephei system, assuming fiducial disk parameters as in Sect. 5. We see that the curves in the upper panel match the limiting cases discussed in the previous paragraph; the regions where ${\dot{a}}_{\text{p}}\propto d_{\text{p}}^{-1}$ are joined by an intermediate region where the inspiral velocity depends less strongly on d_p. In all cases, inspiral speed is a monotonically decreasing function of d_p, meaning that small bodies are preferentially flushed out of the system, thus reducing the erosion suffered by larger objects.

In the lower panel we see that for moderately large planetesimals (≳ 0.1 km), there is a noticeable dip in the inspiral rate in the broad region around 1.8 AU. This is because for this disk model v_c vanishes at 1.8 AU, so the only contribution to the inspiral comes from the sub-Keplerian rotation of the gas and not from any relative eccentricity between the planetesimals and gas, significantly reducing ȧ_p.

These considerations highlight important differences in the radial drift behavior between the binary and single stars. First, high planetesimal eccentricities driven by the secular effect of the disk and the secondary result in faster radialdrift in disks in binaries compared to single systems (see Eq. (6)). Second, because higher e_r increases the inspiral rate, planetesimals will be naturally concentrated in regions where e_r is low, and thus where they can grow more easily. This effect does not exist for planetesimals around single stars (unless one includes additional physics).

Fig. 3 Inspiral speed as a function of planetesimal size at different locations in the disk (a) (see legend), and as a function of location in the disk for different planetesimal sizes (b) (see legend). The figure was made assuming the disk model described in Sect. 5, with e0 = 0.01, Σ0 = 1000 g cm−2, and ϖd= 0.

4 Numerical ingredients of the model

Our calculations of planetesimal growth in binaries employ a multi-annulus (or multi-zone) coagulation-fragmentation code specifically designed for this task. Here we briefly summarize its main features and provide references to relevant sections in Appendices A and B.

4.1 Basics of the code structure

Our code models the evolution of the planetesimal size distribution in N_ann discrete spatial annuli placed at different radii from the primary. These should really be thought of as bins in the space of the semimajor axis, rather than radius, but we do not make that distinction in the following discussion. Within each annulus, the planetesimal size distribution is represented using N_bins logarithmically spaced massbins. At any point in time, in a given annulus, all of the information about the size distribution is encoded in a single vector n, such that n_i is the number of planetesimals in mass bin i. To follow the evolution of n_i in each annulus the code performs two basic functions.

First, in each of the annuli, it calculates the evolution of the size distribution due to planetesimal-planetesimal collisions. This is done in a standard fashion, effectively by solving a discretized version of the Smoluchowski equation (Smoluchowski 1916) accounting for the possibility of fragmentation (Sects. A.1, A.2, and A.6.1). Inclusion of fragmentation generally makes this an $O(N_{\text{bins}}^3)$ calculation at each time step, which is numerically expensive. To overcome this problem, our code employs the new fragmentation algorithm developed in Rafikov et al. (2020), for which the numerical cost goes only as $O(N_{\text{bins}}^2)$, as long as the size distribution of fragments formed in a collision of two objects is self-similar (which is a standard and physically motivated assumption anyway). The coagulation-fragmentation component of the code has been extensively tested against the known analytical solutions as described in Sects. B.1 and B.2.

To compute the number of collisions between different mass bins we use collision rates calculated using the relative velocities from Rafikov & Silsbee (2015a) and accounting for both the forced eccentricity and random velocities (see Sects. 3.2 and 3.3). Our collision rates include gravitational focusing and smoothly interpolate between the shear- and dispersion-dominated velocity regimes (Sect. A.4). The full distribution of collision velocities is also used to model collision outcomes following the recipes in Stewart & Leinhardt (2009; see Sects. 3.4 and A.2).

Our implementation of many code components has certain elements of stochasticity in it, which is quite important. It has been previously found (Windmark et al. 2012; Garaud et al. 2013) that including the distribution of collision velocities can qualitatively change the outcome of the coagulation process. In particular, in a study of dust growth, Windmark et al. (2012) found that using a Maxwellian collision velocity distribution instead of a delta function at the rms velocity of the Maxwellian can cause a few particles to experience a series of lucky low-velocity collisions and grow to larger sizes; the larger bodies formed via lucky collisions are more resistant to destruction in subsequent collisions and would continue to grow to arbitrarily large sizes. In our case, the strong variation in collision velocity with collision partner size provides the dominant source of such randomness. In addition, even for bodies with given sizes, there is a range in their relative collision velocity (see Sect. A.3). To account for such a possibility, we also draw collision velocity between two bodies from a physically motivated distribution (Rafikov & Silsbee 2015a; see Sect. 3.3), as described in Sect. A.3.

The second operation performed by the code is the radial redistribution of mass between different annuli caused by the size-dependent inward migration of planetesimals due to gas drag (Sect. 3.5). The implementation of this procedure and its tests are described in Sects. A.5, A.6.2, and B.3, respectively.

Our model implicitly assumes that other than redistribution caused by radial drift, planetesimals interact only with other planetesimals within their semimajor axis bin. In practice this is not strictly true: Planetesimals in neighboring bins can also collide. However, because the growth we are modeling occurs in regions where the relative eccentricity between planetesimals of different sizes is very small (typically less than 1%), we expect the coupling of planetesimals with those in adjacent semimajor axis bins to be a minor effect.

4.2 Parameters of the numerical model

Multiple modules comprising our code employ a number of parameters, both numerical and physical. All our simulations use a standard set of these parameters described here and listed in Table C.1.

Our grid in mass space employs N_bins = 100 bins; the smallest size (corresponding to the lowest bin) below which mass is supposed to be flushed out from the system due to gas drag is d_p = 10 m. The largest mass bin in our simulations corresponds to an object with d_p = 300 km, since beyond this point destructive collisions are no longer able to prevent growth to larger sizes. This results in a mass ratio between (logarithmically uniformly spaced) mass bins of M_i+1∕M_i = 1.37. To represent our global simulation domain extending from 1 AU to 4 AU we use N_ann = 60 annuli, spaced in semimajor axis as described in Sect. A.5. In choosing this radial range we are motivated by the planetary semimajor axes in several S-type systems: a_pl = 2 AU in γ Cephei, a_pl = 2.6 AU in HD 196885, and a_pl = 1.6 AU in HD 41004 (Chauvin et al. 2011).

Some additional parameters used in our simulations are as follows: the slope of the fragment mass spectrum ξ = −1 (Sect. sectionlinking A.2); parameter b = 10⁻² used to decide whether the largest fragment is distinct from the continuous spectrum of fragments (Sect. A.2); tolerance parameters ϵ₁ = 0.05 and ϵ₂ = 10⁻⁶ used in timestep determination (Sect. A.6).

We studied the sensitivity of the code to changes of these parameters and found only slight variations in the outcomes (see Appendix C).

5 Results: Fiducial model

We start presentation of our results by describing in detail a particular (fiducial) simulation, which will also illustrate the metric we employ to determine the success of planet formation in a given disk model (Sect. 5.2.1). Our fiducial disk model uses the orbital parameters of the γ Cephei system and disk characteristics as described in Sect. 2. For the disk and planetesimal parameters it adopts e₀ = 0.01, Σ₀ = 1000 g cm⁻² at a₀ = 1 AU, ϖ_d = 0°, σ_i = 10⁻⁴, and a solid to gas ratio of 0.01. The total disk mass is 3.7M_J. We present the results for other disk models in Sect. 6.

To assist in interpreting the outcomes of our calculations, we first describe the results of several one-zone simulations of planetesimal growth at different semimajor axes in Sect. 5.1. We then move on to the global, multi-zone calculation of planet formation in the fiducial model in Sect. 5.2. This approach not only allows us to separate local and global factors affecting planet formation but also highlights the role played by the gas drag-driven radial inspiral (naturally absent in the former case).

5.1 Coagulation in a single annulus

We ran our single-annulus calculations of planetesimal growth at three semimajor axes – a_p = 2 AU, 2.3 AU, and 2.5 AU – in our fiducial disk model. These are the same locations for which Fig. 2 illustrates the distribution of collision velocities, which helps us in understanding the role of planetesimal dynamics on their growth. The simulations started with all planetesimals having the same size d_init = 4.3 km; for the adopted bulk density ρ_p = 3 g cm⁻³ this corresponds to planetesimal mass of m = 10¹⁸ g.

Figure 4 illustrates how the size distribution in our model system changes as a function of time. What is actually shown is mdN∕d ln d_p, a variable equal to the mass per unit log planetesimal size. Different curves correspond to different times since the beginning of the simulation, as indicated bythe color bar. The simulation for panel (a) was stopped when the largest body went beyond the largest mass bin in the simulation (d_p = 300 km). In panels (b) and (c), the simulation was stopped after 10⁶ yr.

In all plots, the curves corresponding to times ≲10⁴ yr display features due to the initial conditions and finite spacing of the mass bins. Initially, all of the mass is distributed in just two mass bins surrounding the initial (seed) size d_init, which is marked by the dashed vertical black line. The wiggles just to the right of d_init are related to the finite spacing of the mass bins. The gap just to the left of d_init exists because debris of that size is not created in collisions of objects with size d_init – they produce only smaller fragments to the left of the gap. The size spectrum of objects to the left from the gap is universal at early times (~ 10³–10⁴ yr), as it simply reflects the size distribution of fragments formed in collisions of roughly equal mass planetesimals of size d_init. The gap region is only later filled in by coagulation of the debris particles produced early on and erosion of the seed bodies by small debris.

Over time, memory of the initial conditions is erased, the gaps near d_init get filled, and the size distribution becomes smooth. There are some other notable features that develop in the size distribution at late times. First, after a few times 10⁴ yr mdN∕d ln d_p starts featuring a dip around d_p ~ 0.1 km, especially pronounced in panel (b). Its location corresponds to the planetesimal size which is most easily destroyed incollisions (see Rafikov & Silsbee 2015b). In each panel, the vertical dotted line corresponds to d_c for that simulation, which varies by a factor of several across the different environments, but stays pretty close to d_p ~ 0.1 km. Particles with sizes aroundd_c experience high-velocity collisions with almost all other bodies, and are therefore preferentially destroyed, which likely affects the depth of thedip in different panels.

Second, one can see (especially in panel (b)) mdN∕d ln d_p developing power law behavior to the left and right of the d_p ~ 0.1 km line, with different slopes. We believe this behavior to be real but will refrain from offering an explanation for the slope of these segments. This cannot be done based on standard results regarding fragmentation (O’Brien & Greenberg 2003; Tanaka et al. 1996) since they assume a power-law dependence of the collision rate on colliding partner size.

Third, panel c exhibits an accumulation of particles of 10–20 m in size; a similar accumulation is also present in panel b as a small bump in the same size range. This feature is artificial and arises because in our calculations we do not follow the fate of debris smaller than 10 m in size, simply removing it from the system. As a result, 10–20 m objects do not have smaller bodies to erode them, and their own relative velocities are too low to result in catastrophic disruption: because of their comparable size their collisional velocities are small, thanks to the peculiarity of planetesimal dynamics in binaries (see Fig. 2). Nevertheless, the emergence of this feature near the bottom of our mass grid does not affect the outcome of planetesimal evolution for d_p ≳ d_init.

Careful examination of our simulation outputs shows that the largest objects grow primarily in collisions with objects of similar size, certainly for d_p > d_init. This is because most solid mass is contained in such objects, but also because growth is most favorable in collisions of such objects.

We observe, not surprisingly, that higher v_c suppresses growth. In panel (a), with v_c = 77 m s⁻¹, growth to hundreds of km occurs within several ×10⁵ yr. This can also be seen in Fig. 5 where we show the size of the largest body as a function of time for the three local simulations discussed here. Were we to simulate larger bodies, beyond 300 km, growth would continue in panel (a) until the mass reservoir is fully depleted. Because of low v_c, even when asignificant amount of solid mass in the system has been converted to debris (at late times), collisions of the largest objects withdebris particles are not energetic enough to significantly impede their growth.

In panel b with v_c = 265 m s⁻¹, planetesimals can only grow to a few times the initial size. After that growth of the largest objects stalls (see also Fig. 5), while the population of objects with d_p ~ d_init gets gradually eroded by smaller objects produced in previous collisions (this can be seen in the decay of amplitude of dN∕d ln d_p around that size). If we were to run this simulation beyond 1 Myr, large objects with d_p ~ d_init would eventually be eroded, followed by smaller bodies also grinding themselves down.

Finally, in panel c of Fig. 4, collision velocities are so high that there is barely any growth – the largest objects present in the system at early time have grown in size beyond d_init only by ~2. But after a population of small debris develops in the system (already by 10³ yr), Fig. 5 shows that large bodies get eroded away by 10⁵ yr, leaving only a small population of debris near the lower end of the simulation mass range (which is an artifact of our calculation, as we mentioned earlier).

Figure 5 shows that in all three environments, the largest bodies initially evolve in a very similar fashion, steadily growing in size. This universality is caused by the fact that initial growth occurs by mergers of objects with d_p = d_init, which have low relative velocities. The differences emerge only after ~10⁴ yr when a sufficient amount of small debris capable of efficiently eroding big bodies in high-v_c environments accumulates in the system.

Here we again consider Fig. 2, which illustrates, in particular, collision outcomes as a function of the sizes of the collision partners, for the same environments that the simulations shown in Fig. 4 were run for. Interestingly, we see that to halt planetesimal growth it is not necessary to have catastrophic disruption of objects at d_p ~ d_init. Indeed, although Fig. 4 shows that growth is halted in panels (B) and (C), Fig. 2 shows no catastrophic disruption of planetesimals at the initial size: the point (d_init, d_init) is outside the solid white contours for all dynamic environments, as d_init = 4.3 km considerably exceeds d_c in all panels.

It may seem strange that in all three panels in Fig. 2 the point (d_init, d_init) also lies outside the dotted contours, bordering the region of erosive encounters – it would seem natural that one should then find growth starting with a population of d_p = d_init objects in all three environments. The resolution of this apparent paradox involves two factors. First, according to our definition of erosion in Fig. 2 dotted contours correspond to collisions, in which the mass of the largest remnant is equal to the mass of the bigger object (target) involved in the collision; the mass equal to the mass of the smaller object (projectile) is reduced to debris. Similarly, even though the point (d_init, d_init) sits outside the dotted contour, some amount of debris will still be produced in collisions of two objects of initial size d_init; the amount of created debris is larger the closer this point is to the dotted contour (i.e., in panels b and c relative to (a)). And once some debris particles are created, Fig. 2 shows that their collisions with seed bodies (d_p = d_init) are very erosive, especially in panel (c), but only barely so in panel a.

Second, relative velocity maps in Fig. 2 assume a particular (mean) value of encounter velocity, whereas in practice v_coll has a reasonably broad distribution around this mean (see Sect. A.3). As a result, some collisions occur at higher velocities and are more destructive – a fact that is fully accounted for in our coagulation-fragmentation simulations depicted in Fig. 4. This is another clear illustration of the importance of considering the distribution of planetesimal velocities when studying their growth. The combination of the two aforementioned factors leads to erosion eventually becoming the dominant player in the collisional evolution in panels b and c of Fig. 4.

Fig. 4 Evolution of the size distribution as a function of time in a set of single-annulus simulations, described in Sect. 5.1, which should be consulted for details. The three panels correspond to environments with the different values of the critical velocity vc and critical size dc, previously shown in Fig. 2, at different locations (as labeled) in the disk around the primary star in the γ Cephei system. The vertical dotted lines correspond to the critical size dc, and the vertical dashed lines correspond to the initial size of the planetesimals dinit.

Fig. 5 Time evolution of the size of the biggest object in the system, shown for the three simulations depicted in Fig. 4 and labeled according to their semimajor axis.

5.2 Full global calculation with inspiral

We now proceed to describe a fully global, multi-annulus simulation of planetesimal growth, fully accounting for the effects of size-dependent radial inspiral (see Sect. 3.5) and using the fiducial disk parameters. This simulation is also used to motivate our choice of a particular metric determining whether planetesimal growth in the system successfully overcomes the fragmentation barrier (see next section).

5.2.1 Outcome diagnostic

In any reasonable environment, planetesimal growth would proceed to very large sizes, eventually leading to planet formation, provided that the initial planetesimal size d_init is large enough. Indeed, objects withsizes d_p ≫ d_c do not experience catastrophic disruption and get eroded only by much smaller planetesimals, which cannot compete with the addition of mass in collisions with larger bodies. Thus, given large enough d_init, planet formation is guaranteed to be successful (e.g., Thebault 2011) even in dynamically harsh environments of the binary stars.

However, in many environments, if the initial planetesimals are too small, then they will be eroded or destroyed and large bodies will not form (like in panels b and c of Fig. 4). For this reason, the appropriate question to ask is not whether planetesimal growth can occur, as the answer would depend on the initial condition – the size of the seed planetesimals d_init. Instead, one should be trying to figure out from what initial planetesimal size d_init can sustained planetesimal growth occur. In a given dynamical environment, there will be a minimum size d_min such that if the initial bodies are smaller than d_min then planetesimal growth in the system will eventually be halted, as in Fig. 4b,c, whereas if seed objects have d_init > d_min, then planetesimal coagulation will eventually form large bodies, thus overcoming the fragmentation barrier.

This question, in principle, can be asked locally, at a given semimajor axis in the disk. In this case Fig. 4 makes it clear that, for a fiducial disk model, d_min < 4.3 km at a_p = 2 AU, whereas d_min > 4.3 km at a_p = 2.3 and 2.5 AU. However, the calculations presented in Sect. 5.1 do not account for radial mass transport due to gas drag and may thus be inaccurate. For this reason, we will be more interested in a general question of what d_min is needed for the fragmentation barrier to be overcome and for a planet to form at some location in the whole disk, given its parameters. Once d_min is determined, some other information, for example, the location where planetesimal growth is fastest, would be the outcome of our global calculations, giving them certain predictive power. Of course, d_init may (and probably should) vary across the disk. Moreover, at each location, a whole spectrum of initial sizes is likely to be present. However, in this work, to reduce the number of degrees of freedom, we assume that seed planetesimals have the same size d_init across the whole disk and determine a single value of d_min for a given disk model based on this assumption.

To determine d_min defined in this manner, we employed the following iterative algorithm. We chose a value of d_min and used that as our starting size d_init. The simulation was run until one of three conditions were met: (i) 1 Myr has passed, (ii) a 300 km body is formed, or (iii) there is no longer sufficient solid mass remaining in the simulation to produce a 300 km body. If a 300 km body forms, then we know d_min < d_init, whereas if it does not, then d_min > d_init. By trying several values of d_init and running our global simulations for each of them, we were able to first bracket d_min within someinterval, and then to converge to its value with high accuracy. This procedure is demonstrated next (allowing us to also illustrate the sensitivity of global planetesimal growth to d_init).

Since the planetesimal composition is uncertain, in addition to our default simulations that use the material properties of “strong rock" bodies from Stewart & Leinhardt (2009), we also ran simulations using the strength parameters for their “rubble pile” planetesimals. We denote the minimum size assuming rubble pile planetesimals as $d_{\text{min}}^{\text{rubble}}$.

5.2.2 Evolution of the system as a whole

Figure 6 shows the evolution of the planetesimal size distribution in the entire disk for several runs of our fiducial simulation with different initial sizes d_init expressed in terms of d_min = 1.6 km, as determined for this disk model. Curves of different color correspond to different times since the simulation was started.Because collisional processes proceed at different rates across the disk, reflecting the complexity of planetesimal dynamics in binaries (see Fig. 1), the evolution of the planetesimal mass spectrum looks somewhat different from that in a single annulus (see Fig. 4; e.g., the initial gaps to the left of d_init are barely there, and there is no pile-up at small sizes). This is to be expected since the global size spectrum in Fig. 6 is an average of multiple single-zone size distributions (such as the ones shown in Fig. 4) with the added complication of radial inspiral of solids. Nevertheless, some key features of the size distribution evolution (e.g., a dip in mdN∕d ln d_p near d_p ~ 0.1 km, mass flow toward small sizes, and so on) are still preserved even in the full disk calculations.

We can see that in panels a and b, essentially no growth occurs beyond the initial size. The behavior shown in these panels differs only in how rapidly the bodies are ground down. The lack of any appreciable growth in panel b with initial size only a factor of two less than d_min shows that there is a sharp transition over less than a factor of two in d_init from growthto large sizes to essentially no growth whatsoever.

In panels c and d, coagulation to 300 km does occur. In both cases, though, less than half of the initial solid mass in the simulation remains at the end. The time required to produce a 300 km body is almost four times less for the simulation with d_init∕d_min = 5 than for the one with d_init∕d_min = 2 (8 × 10⁴ yr vs. 3 × 10⁵ yr), and in the former case more than twice as much solid material remains (47% of the original 7 M_⊕ in panel d vs. 21% in panel c). In the simulation with d_init∕d_min = 2 (5), 79% (88%) of the mass is lost to the creation of rubble smaller than the smallest body tracked in the simulation, with the remaining amount lost to inspiral into the central star.

Figures 7 and 8 provide details on how the coagulation-fragmentation process proceeds as a function of semimajor axis in the disk, for each global simulation shown in Fig. 6. Figure 7 shows dM∕da – the total solid mass (in planetesimals of all sizes) per unit a, whereas Fig. 8 depicts the size of the largest body (in a given annulus), both as functions of a and time. In Fig. 7, the horizontal dashed violet line represents the initial mass distribution, which is flat because of our chosen Σ profile (see Eq. (1)). In Fig. 8 that line corresponds to the size of seed objects d_init, which is constant across the disk according to our assumption (see Sect. 5.2.1).

To illustrate how the dynamics of planetesimals affect their size evolution, we also show in Fig. 9 the run of the collisional velocity v_coll between d₁ = 2 and d₂ = 5 km planetesimals as a function of semimajor axis (for many of the different simulations carried out in this work, see Sect. 6). It accounts for both the secular and random velocity contributions and is calculated using Eq. (A.16) with parameters λ and e_rand taking their mean values. The black solid line corresponds to the fiducial disk model, whereas the black dotted line in panel (a) shows v_coll computed in the absence of random motions. One can easily recognize the “valley” forming around 1.8 AU as the dynamically quiet location where v_c, d_c → 0 in Fig. 1. Nonzero random velocities of planetesimals “fill" this valley, not allowing v_coll to drop to zero in their presence. Nevertheless, v_coll is still relatively low around this region.

A large increase in v_coll in the outer disk starting around 2 AU arises because for Σ₀ = 10³ g cm⁻² the free precession rate A becomes very small there, while the excitation by the companion grows large. As a result, both v_c and v_coll become very large, of order km s⁻¹, even though a true secular resonance does not appear in this disk model. We note that v_coll is often several times smaller than v_c, because d_c is often less than d₁ = 2 km.

Figure 7 shows that, in the cases with d_init∕d_min < 1, most of the solid mass is removed from the outer and inner portions of the disk on timescales of 10⁴ yr, leaving the majority of the remaining solids in a narrow band around 1.8 AU. This band coincides with the dynamically quiet location in the fiducial disk model (see Fig. 9), where (i) the velocity of inward drift of planetesimals goes down (see Fig. 3) and (ii) v_coll is relativelylow. These factors naturally cause planetesimal accumulation (dM∕da temporarily exceeds its initial value due to arrival of mass from larger radii) and promote growth at this location (see Fig. 8). But over time, the mass concentration around 1.8 AU is eroded and moves inward due to gas drag.

In the outer disk, outside 2 AU, both the amount of mass and the size of the largest object drops precipitously early on, primarily because of very high v_coll ~ 1 km s⁻¹ in this region (see the black curve in Fig. 9) resulting in rapid planetesimal erosion down to the smallest size (10 m) tracked in our simulations². Some of thisdebris gets swept into the dynamically quiet zone through inspiral before it is fully ground down. In the inner disk, below 1.8 AU, v_coll is much lower (enabling some temporary growth for the largest objects) and mass is removed predominantly by gas drag.

In the simulations in which coagulation is successful (d_init∕d_min > 1), the mass is again quickly removed from the outer part of the disk. But in the inner disk, particularly around 1.8 AU, the mass density (per unit semimajor axis) remains roughly constant as large objects are able to efficiently accrete smaller ones, preventing them from getting lost via inspiral or grinding down. In the run with d_init∕d_min = 5 we find that in the outer disk, the largest bodies have size (2−3)d_init for the duration of the simulation. This illustrates that larger seed objects find it much easier to survive than smaller ones even in the dynamically harsh environments.

Fig. 6 Evolution of the mass distribution in the entire disk in our fiducial multi-zone model (Simulation 1). Different panels are produced for different initial planetesimal sizes dinit in terms of dmin (as labeled on the plot), which is dmin = 1.6 km for this simulation.

Fig. 7 Mass in planetesimals of all sizes (tracked in our mass grid) in the disk per unit semimajor axis dM∕da = 2πaΣ(a) as a functionof a in our fiducial model at different moments of time (see color bar). The different panels correspond to the same initial planetesimalsizes dinit as in Fig. 6. The dashed violet horizontal line corresponds to the mass distribution at t = 0.

Fig. 8 Similar to Fig. 7 but showing the size of the largest body as a function of the semimajor axis. The dashed violet horizontal line now shows the size of seed planetesimals dinit.

Fig. 9 Collision velocity between d1 = 2 and d2 = 5 km planetesimals as a function of the semimajor axis a for the different values of model parameters varied one at a time (with respect to the fiducial set of parameters) and shown in the subplot legends. The fiducial value of the varied parameter in each subplot is shown with the solid black curve. The dotted line in panel (a) shows the collision velocities in the limit that σi = 0.

Fig. 10 Dependence of our key metric dmin – a minimum size of seed planetesimals that results in their eventual growth beyond 300 km – on some parameters of our simulations: (a) surface density normalization Σ0, (b) disk eccentricity normalization e0, (c) apsidal angle of the disk ϖd, (d) amplitude of the random velocity of planetesimals σi, and (e) parameter χir (artificially) regulating the magnitude of planetesimal inspiral. Solid lines and diamond points correspond to dmin, whereas dashed lines and circles are for $d_{\text{min}}^{\text{rubble}}$ – the minimum size for collisionally weak planetesimals. The causes of the trends shown here are discussed in Sects. 6.1–6.7. In each panel, the fiducial simulation is highlighted in red.

6 Results: Variation of system parameters

Having discussed in detail the evolution of the global planetesimal size distribution for a particular disk model (Sect. 5), we proceed to explore the changes brought in by varying disk model parameters relative to their fiducial values. We vary several key physical parameters – surface density normalization Σ₀, disk eccentricity normalization e₀, disk orientation ϖ_d, amplitude of random velocity σ_i – which control planetesimal dynamics. We also run calculations in which we artificially turn off disk gravity (but not gas drag), and vary the inspiral rate, to explore the role played by these physical processes in determining the outcome of planetesimal evolution. In carrying out this parameter exploration we use the minimum planetesimal size d_min (at which growth beyond 300 km becomes possible) as our metric to characterize the success of planet formation in a given model.

Table 1 lists the parameters of the different models that have been explored, as well as their d_min and $d_{\text{min}}^{\text{rubble}}$. These two simulation metrics are also shown in Fig. 10 and are discussed in the rest of this section. Simulation 1 is our fiducial model, and all other simulations vary one or two of the parameters while leaving the rest fixed. In each row, the parameters that differ from Simulation 1 are highlighted in bold. We also vary a number of numerical parameters, and show that these choices make little difference to the outcome (see Appendix C and Table C.1 for details).

6.1 Variation of surface density: Simulations 2–7

Somewhat unexpectedly, we find that variation in the surface density normalization Σ₀ within a broad range relative to the fiducial value of 10³ g cm⁻² has almost no effect on d_min, which ranges between 1.1 and 1.7 km (see Fig. 10a and Table 1). Even though both the dynamically quiet location at which v_c and d_c vanish and the location (and existence) of the secular resonance change with Σ₀ (see Fig. 1), this does not dramatically affect d_min.

Figure 10a does show a weak decrease in d_min as Σ₀ increases from 500 to 5000 g cm⁻². This is because, as in the fiducial model, the most favorable conditions for planetesimal growth occur at the dynamically quiet location, where v_c and d_c vanish, leaving only random velocity to produce relative velocity between the colliding planetesimals (i.e., v_coll ~ σ_iv_K). As Σ₀ increases, the dynamically quiet location moves out in the disk. Since the Keplerian speed drops with distance, planetesimals collide at lower velocities at the more distant dynamically quiet locations (i.e., for higher Σ₀). As a result, within the range of Σ₀ where the dynamically quiet location exists, collisions are less destructive and planetesimal growth becomes possible starting at lower d_init for larger Σ₀.

We note that in the Σ₀ = 500 g cm² model a secular resonance (Silsbee & Rafikov 2015b) appears in the disk around 2.4 AU (see Fig. 1), driving v_coll to very high values³ at moderate (2–2.5 AU) semimajor axes. However, d_min is only slightly larger than for other disk models (not featuring resonances), simply because most growth occurs at the dynamicallyquiet location, which is far from this resonance.

The situation changes for even more massive disks, Σ₀ > 5000 g cm⁻², in which thedisk gravity is dominant over gravity from the companion throughout the whole disk. As a consequence, such disks have no dynamically quiet location⁴ and we find a slight increase in d_min with Σ₀.

In Simulation 7, we also considered a model of a massive axisymmetric disk, with Σ₀ = 20, 000 g cm⁻² (corresponding to a total disk mass of 0.07 M_⊙ interior to 5 AU) and e₀ = 0. In this setup, previously explored in Rafikov (2013b), disk gravity does not give rise to planetesimal eccentricity excitation but effectively suppresses the dynamical excitation due to the secondary’s torque. As a result, in this model collision velocities are reduced even further (compared to, e.g., Simulation 6), and we find a significantly lower d_min = 0.7 km (see Table 1).

6.2 Variation of e₀: Simulations 8–12

Figure 10b shows that disk eccentricity has a dramatic effect on planetesimal growth: d_min varies by 1.4 dex as the disk eccentricity normalization e₀ changes from 0 to 0.1. This variation is non-monotonic.

As e₀ increases from zero to the fiducial value e₀ = 0.01, we find d_min decreasing for reasons similar to the ones mentioned in Sect. 6.1: Higher e₀ means larger |B_d |∝ e₀Σ₀ (Silsbee & Rafikov 2015b; Rafikov & Silsbee 2015a), shifting the dynamically quiet location further from the primary (here we again assume ϖ_d = 0°; see Fig. 1). As v_coll ~ σ_iv_K ∝ a^−1∕2 at this location, higher e₀ results in lower v_coll (see Fig. 9c) promoting planetesimal growth and leading to lower d_min. In an axisymmetric disk (e₀ = 0), B_d = 0 and cancellation of the secondary’s torque does not occur, leading to a rather high d_min = 4.2 km. The dynamically quiet location exists in our simulation domain for values of e₀ between 0.004 and 0.018.

For e₀ > 0.018 the torque due to disk gravity starts to dominate planetesimal eccentricity excitation, and v_c is nonzero throughout the whole simulated portion of the disk; the dynamically quiet zone disappears. In this regime v_c ∝|B_d|∝ e₀ (Silsbee & Rafikov 2015b), leading to higher v_coll and a less growth-friendly environment as e₀ increases (see Fig. 9c); hence the rapidly increasing values of d_min. Even a disk with e₀ = 0.05 requires seed planetesimals with d_min ≈ 22 km to ensure their growth into planetary regime.

We note that while the location of the dynamically quiet region is sensitive to e₀, variation of disk eccentricity does not affect the presence (or location) of the secular resonance in the disk. The latter is determined by the disk mass only, while the former depends on both Σ₀ and e₀ (Rafikov 2013b; Silsbee & Rafikov 2015b; Rafikov & Silsbee 2015a).

6.3 Variation of ϖ_d: Simulations 13–15

We also find substantial sensitivity of d_min to ϖ_d. Figure 10c shows that increasing the deviation of the disk orientation from apsidal alignment with the binary orbit leads to larger d_min: Keeping everything else fixed, d_min grows from 1.6 km for an aligned disk (ϖ_d = 0°) to 9.7 km for an anti-aligned disk (ϖ_d = 180°). Even a ϖ_d = 10° misalignment is sufficient to increase d_min to 2.5 km.

This trend, again, owes its origin to the presence (or absence) of a dynamically quiet location in the disk. Secularly forced v_c can vanish only in apsidally aligned disks (see Fig. 1a). A small misalignment makes v_c nonzero everywhere (Rafikov & Silsbee 2015a,b) but only slightly modifies the behavior of v_coll (see, e.g., the ϖ_d = 10° curve in Fig. 9d), because (i) v_c is still dramatically reduced at the location where Eq. (5) is satisfied and (ii) random velocity σ_i is nonzero. However, for a significant misalignment (ϖ_d ~ 1), collisional velocities end up being much higher than in the aligned disk (see ϖ_d = 90° and ϖ_d = 180° curves in Fig. 9d), requiring larger seed planetesimals to overcome the fragmentation barrier. Therefore, it is reasonable that d_min increases with the degree of misalignment, as our results show.

6.4 Variation of σ_i: Simulations 16–19

Figure 10d demonstrates a strong, monotonic dependence of d_min on the magnitude of random motions σ_i: for collisionally strong planetesimals d_min rises from 0.3 km at σ_i = 10⁻⁵ to almost 30 km at σ_i = 10⁻³. This is to be expected for disk models in Simulations 16–19, which all feature a dynamically quiet zone around 1.8 AU (set by the values of Σ₀ and e₀). In this narrow zone v_c is very low and v_coll is set entirely by random velocity σ_iv_K (see the black solid and dotted curves in Fig. 9a). As a result, even at this favorable location in the disk, d_min is limited by a combination of radial inspiral and random motions that can destroy small planetesimals if σ_i gets large, requiring larger seed bodies to ensure steady growth. We anticipate that σ_i would be less important in models which do not have a dynamically quiet region (see Sect. 6.6).

6.5 Strong versus weak planetesimals

Dashed curves in Fig. 10 show our results for collisionally weak, rubble-pile planetesimals as defined in Stewart & Leinhardt (2009). Rather unsurprisingly, we find that for such planetesimals $d_{\text{min}}^{\text{rubble}}$ is larger than d_min for collisionally stronger objects, typically by a factor of 1.5–3. Thus, more massive seed bodies would be necessary to enable sustained growth if they are collisionally weak.

6.6 The case without disk gravity: Simulations 20–24

In order to better understand the impact of properly accounting for disk gravity and to isolate its effect on planetesimal growth, we also ran simulations where disk gravity was artificially turned off in our secular solutions for planetesimal eccentricity, that is, we set B_d = A_d = 0. In this case, the secular eccentricity of planetesimals is set only by the gravity of the companion and gas drag (Thébault et al. 2008). We carried out this exercise while varying σ_i, with the results presented by orange curves and points in Fig. 10d.

Turning off disk gravity eliminates the region of vanishing v_c (see Fig. 1a) and collision velocities around 1.8 AU go up substantially (see Fig. 9a). Therefore, as one might expect, planetesimal growth requires significantly higher d_min = 6.9 km, compared to 1.6 km when disk gravity is fully accounted for (assuming the fiducial σ_i = 10⁻⁴). We note that without diskgravity v_coll in the outer disk is lower than with it, but this still does not help to keep d_min low.

In the absence of disk gravity, reducing σ_i has very little effect on d_min because for σ_i ≲ 10⁻⁴, destructive collisions arise mostly due to secularly excited eccentricities and are therefore unaffected by further reducing σ_i. But at sufficiently high σ_i, the collision velocities are dominated by random motions, and d_min is no longeraffected by the disk model. We do not consider such high values of σ_i, but we do find at σ_i = 10⁻³ that d_min is lower in the model with disk gravity than without. This reversal occurs because for σ_i = 10⁻³, the outer part of the disk becomes the most favorable environment for planet formation in the absence of disk gravity because the random velocities are lower there. With disk gravity the collision velocities in the outer disk are still dominated by the secular eccentricities, so planet formation must occur in the inner disk where the random velocities are higher.

Interestingly, d_min increases very slightly as we go to very low values of σ_i, going from 6.6 km at σ_i = 3 × 10⁻⁵ to 6.9 km at σ_i = 10⁻⁵. Possibly that is because decreasing σ_i decreases the timescale for collisional evolution (see Eq. (A.25)) while leaving the timescale for radial inspiral unaffected. As a result, erosion in collisions is more of an issue at low σ_i because small bodies have less time to be removed from the system before they collide with larger ones. But quantitatively, this effect is rather small, and in fact not present for the rubble-pile planetesimals.

6.7 The effect of artificially altering the inspiral rate: Simulations 25–36

The key factor not allowing us to treat the global size evolution of planetesimals as a simple superposition of one-zone coagulation-fragmentation simulations such as the one presented in Sect. 5.1, is the mass exchange between the different annuli caused by radial inspiral due to gas drag. Thus, to determine the degree to which the growth process is affected by the mass exchange, we also carried out some simulations with the inspiral rateartificially increased or decreased compared to the value given by Eq. (6) by a factor χ_ir as indicated in Table 1. Simulations 25–31 explore what happens in the fiducial model when the inspiral rate is varied with χ_ir changing from 0 to 100.

One possibility that we wanted to test through this exercise is that the removal of small planetesimals by inspiral should reduce the detrimental effect of erosion on planetesimal growth, with the expectation that high χ_ir > 1 (i.e., faster removal of small fragments) would facilitate growth, while χ_ir < 1 (weakened inspiral) would suppress it. However, Fig. 10e reveals a pattern of d_min behavior that runs counter to this expectation. In fact, in the no-inspiral case (χ_ir= 0) we find d_min ≈ 1 km to be lower than in the case with full inspiral (χ_ir = 1) for which d_min ≈ 1.6 km. And χ_ir > 1 results in even larger d_min.

We believe that there are two reasons for this behavior. First, in the outer disk, inspiral does not help much with cleaning out small fragments, simply because they get primarily ground down below the smallest size captured in our simulations. Second, the radial inspiral negatively impacts the growth of large planetesimals by shifting their semimajor axis away from the dynamically quiet zone where their growth is most favorable. Indeed, Fig.11a shows the results of a no-inspiral simulation with fiducial parameters (χ_ir= 0, analogous to a simple addition of a number of multi-zone simulations not interacting with each other) for d_init∕d_min = 2, which can be directly compared to Fig. 8c (for which χ_ir = 1). One can see that without inspiral, the largest objects form closer to a_dq = 1.8 AU where v_c = 0, which certainly facilitates their growth compared to the full-inspiral calculation, in which the 300 km object emerges interior to 1.4 AU. But the general conclusion that we can draw from this exercise is that the gas drag-driven inspiral has a rather modest effect on the outcome of planetesimal growth.

Fig. 11 Evolution of the size of the largest object at different locations in the disk for two models used to illustrate certain aspects of our calculations. (a) Fiducial model with gas drag turned off, χir = 0, computed for dinit∕dmin = 2, to be compared with Fig. 8c (see Sect. 6.7 for details). (b) Fiducial model computed for dinit∕dmin just slightly exceeding unity, to be compared with Fig. 8c,d (see Sect. 7.1 for details).

7 Discussion

The results of the previous sections highlight the important differences between planet formation around single stars and in binaries. In the binary case, erosion plays a much more important role in planetesimal size evolution, which is caused by the obvious differences in underlying dynamics: whereas in disks around single stars planetesimal motions are excited only by gravitational perturbations between numerous planetesimals (and, possibly, turbulence in the disk), planetesimals in binaries get additionally (and very strongly) excited by the gravity of the companion star and the disk. This can easily raise the collision speeds of planetesimals of comparable size into the km s⁻¹ regime (see Fig. 9), much higher than in disks around single stars. As a result, growth to 10² -10³-km size objects in binaries is possible only if planetesimals start already rather large, with d_init ~ 1−10 km.

This places planets in binaries such as γ Cep in a unique category of dynamical extremophiles or hyperdynamophiles⁵ – objects capable of surviving and growing despite the harsh dynamical environment in which they reside. The very existence of such planets puts planet formation theories to a very stringent test.

Our parameter space exploration in Sect. 6 is not entirely exhaustive: because of the multi-dimensional space of physical parameters we could vary at most one or two of them at once with respect to our fiducial model, meaning that some domains of the parameter space are left unexplored. Also, our models rely on certain simplifying assumptions that were necessary to reduce the number of model inputs, for example, the initial planetesimal size d_init is the same everywhere in the disk, the amplitude of random motions σ_i is the same regardless of location and object’s size, etc. Nevertheless, our results do allow us to robustly single out the most important physical ingredients determining the outcome of planetesimal growth in binaries.

The gravitational effect of the protoplanetary disk is the main such ingredient and should not be overlooked. The inclusion of disk gravity significantly reduces d_min in our fiducial case, by about a factor of four compared to the case in which disk gravity is ignored (see Fig. 10d). Eccentricity e₀ and orientation ϖ_d of the protoplanetary disk are two other key parameters determining the outcome of planetesimal growth. Acting in conjunction with disk gravity, they can lead to the emergence of dynamically quiet locations in the disk, providing favorable conditions for coagulation.

Inclusion of radial inspiral of solids is what makes our models fully global. We find that inspiral somewhat complicates planet formation by causing the migration of growing planetesimals out of the dynamically quiet zone in the disk and raising d_min compared to the case in which inspiral is ignored (see Fig. 10e). But the overall effect is not as large as one might have expected (see Sect. 6.7). Next we discuss some other consequences of our results.

7.1 Implications for planet formation in binaries

About a dozen binaries are currently known to harbor planets in S-type configurations (Chauvin et al. 2011; Marzari & Thebault 2019), some of them with binary semimajor axis even smaller than that of γ Cep (i.e., a_b < 20 AU). Our results provide a pathway to understanding the origin of such planets. Our previous attempt to explain their existence (Rafikov & Silsbee 2015b) focused on explaining the formation of planets in some of these systems (Chauvin et al. 2011) at their current locations. But in this work we are more interested in the overall possibility of planet formation at one to a few AU separation in a γ Cep-like system.

Our finding that d_min only weakly depends on Σ₀ is quite important. As we vary Σ₀ from 500 g cm⁻² to 20,000 g cm⁻² (i.e., between 0.3 and 12 times the minimum mass solar nebula density at 1 AU Hayashi 1981), implying the variation of the protoplanetary disk mass from 1.8 M_J to 70 M_J (for the assumed Σ profile (1) with a₀ = 1 AU and a_out = 5 AU), d_min changes by less than a factor of two. What this means physically is that, even with the variation of Σ₀, the dynamically quiet region can still exist somewhere in the disk, ensuring that planetesimal growth can be achieved in the disk globally (things may get more complicated at very low Σ₀ when a secular resonance appears in the disk; see Fig. 1).

Figure 12 shows a_dq as a function of e₀ for various values of Σ₀ in a system which has, other than the variation of e₀ and Σ₀, our fiducial parameters. We have marked the fiducial model with an orange star along the Σ₀ = 1000 g cm⁻² curve. This shows that over a wide range of disk masses, there is a dynamically quiet region in the disk. For each disk mass however, there is only a relatively narrow (approximately factor of two) range in e₀ over which the dynamically quiet region exists in the part of the disk outside 1 AU.

At low disk masses, a secular resonance appears in the disk (see the curve for Σ₀ = 500 g cm⁻² in Fig. 1). This resonance dramatically increases collision velocities in the outer disk as illustrated in Fig. 9b. Collision velocities between d_p > 10 km planetesimals are even higher than in that figure, since their orbits are less affected by gas drag.

In a static disk model such as ours, this resonance is localized, and planet formation can still proceed in other regions of the disk. In reality though, the disk evaporates over time, which changes the resonance location and causes the resonance to “sweep” through the disk. This effect is similar to the sweeping secular resonances involving planets (Ward 1981).

On the other hand, the dynamically quiet location also sweeps through the disk, meaning that over the disk lifetime, planetesimal growth may be promoted in many parts of the disk. Exactly what is the effect of this combined sweeping dynamically quiet location and sweeping resonance will depend on the details and timing of the disk dispersal.

It is worth noting that Fig. 12 shows that even for a quite low surface density disk with Σ₀ = 500 g cm⁻² (corresponding to 6M_⊕ of solid material), the dynamically quiet region will not exist if e₀ is greater than about 0.02, corresponding to a mass-weighted mean disk eccentricity of 0.05. This is within the range found in numerical simulations, but on the low end of that range (see Sect. 7.3 for details).

The low levels of planetesimal random motions required to drive d_min into the kilometer size range may seem problematic: σ_i = 10⁻⁴ featured in our fiducial model corresponds to random velocities of about 10 m s⁻¹. However, even this is above the level of the escape speed from the surface of our initial planetesimals in most cases, which is the most the two-body relaxation would pump up the random motions to. Moreover, given the high relative velocities (up to ~ v_c) between the planetesimal populations of different sizes, both viscous stirring and dynamical friction are likely to be very inefficient (Stewart & Ida 2000; Rafikov 2003a), further reducing the amplitude of planetesimal random motions.

The precise location and time at which planetesimals reach our assumed threshold size of 300 km are important outcomes of our models. They turn out to be rather strong functions of d_init∕d_min ≥ 1. If the size of seed planetesimals d_init is close to d_min, planetesimals grow rather slowly, causing their substantial gas drag-driven inspiral, so that d_p reaches 300 km somewhat closer to the star than a_dq (see Fig. 11b). As d_init∕d_min increases, planetesimals tend to grow to large sizes closer to a_dq (see Fig. 8c for d_init∕d_min = 2). However, further increase in d_init∕d_min causes the location of preferential planetesimal growth to shift inward again (see Fig. 8d), presumably because of faster coagulation there. Such non-monotonic variation with d_init∕d_min, as well as the possibility of late-time planet migration (that must have certainly affected S-type systems harboring hot Jupiters), complicate the effort to constrain the initial planetesimal size d_init based on the semimajor axes of observed planets within binaries.

At the same time, the time to reach 300 km is a monotonic function of d_init∕d_min, steadily falling from ~1 Myr to < 10⁵ yr as this ratioincreases from 1 to 5 in our fiducial model. The expectation of short disk lifetimes in S-type binaries (Harris et al. 2012; Barenfeld et al. 2019; Zurlo et al. 2020) may then favor larger d_init, for which the growth time is shorter.

As previously stated, our fiducial disk model contains less than four Jupiter masses of material. The planet in the γ Cephei system was determined to have around ten Jupiter masses (Benedict et al. 2018). Unless the disk is being replenished by material infalling from the circumbinary disk, this would rule out all disk models with less than Σ₀ ≲ 5000 g cm⁻² in this particular system.

Fig. 12 Location adq of the dynamically quiet zone in an apsidally aligned disk. We assume the fiducial disk plus binary parameters, except for Σ0 and e0 which are varied as shown on the plot. The fiducial model is marked by an orange star along the curve corresponding to Σ0 = 1000 g cm−2.

7.2 Implications for planetesimal origin

Our finding that some minimum size d_min ~ 1–10 km of seed planetesimals is necessary for ensuring the success of planetary genesis in S-type binary configurations has important consequences for theories of planetesimal formation. Two main ideas for the origin of planetesimals – steady growth by coagulation and rapid formation via streaming instability (Johansen et al. 2014) – predict very different formation pathways.

Growth by particle sticking must necessarily cover all sizes in the range from dust grains to multi-kilometer objects (subject to substantial gravitational focusing). This formation mode is known to be problematic because of the bouncing barrier around millimeter to centimeter sizes (Blum & Wurm 2008) and the rapid radial inspiral threatening meter-size objects (Weidenschilling 1977). In S-type binaries this planetesimal formation mechanism encounters another major bottleneck – efficient growth is impossible until seed objects are 1–10 km in size, substantially larger than any of the aforementioned barriers. This additional obstacle significantly reduces chances of slow planetesimal growth by coagulation in binaries.

On the contrary, the competing mechanism relying on the operation of streaming instability appears to naturallyproduce large planetesimals with sizes in the range from tens to hundreds of kilometers (Simon et al. 2016; Schäfer et al. 2017). This mechanism is also rather fast (with a typical timescale of tens of local orbital periods), resulting in a rapid production of the population of seed planetesimals with d_init exceeding d_min, which would enable further growth into the planetary regime as we have demonstrated here.

Thus, results of our work provide indirect support for the idea that streaming instability is a pathway for forming planetesimals in binaries (similar reasoning can be found in Thebault 2011). Of course, whether and how this instability would operate in an eccentric disk in presence of strong perturbations from the secondary, is an open question. But if it works in binaries, then there is no reason for it not to be the mechanism for forming planetesimals also in disks around single stars.

7.3 Implications for simulations of gaseous disks in binaries

Given the strong sensitivity of planetesimal growth to disk eccentricity e₀ and orientation ϖ_d that we identify in this work, it is imperative to better constrain these variables. Unfortunately, observations are not yet at a stage where one can reliably determine eccentricity and orientation of protoplanetary disks in binaries with small separations a_b ≲ 30 AU. Instead, onehas to rely on numerical simulations for guidance in that matter.

A number of studies have explored the amplitude of disk eccentricity excited by the gravity of an eccentric companion in a coplanar S-type configuration using 2D hydrodynamic simulations. The conclusions that can be drawn from these works very strongly depend on the assumed physical inputs: equation of state, treatment of disk gravity, inner boundary conditions, and so on.

For example, studies using a locally isothermal equation of state (Kley & Nelson 2008; Paardekooper et al. 2008; Marzari et al. 2009; Zsom et al. 2011; Regály et al. 2011; Martin et al. 2020) typically find rather high (mass-weighted mean) disk eccentricities, ⟨e_d ⟩ ~ 0.1–0.3, with e_d profile typically exhibiting a double-peaked structure as a function of semimajor axis and with high values of e_d close to the disk center. These results may need to be viewed with caution in light of the recently identified issues with using the locally isothermal equation of state for numerical modeling of the disk-perturber tidal interaction (Miranda & Rafikov 2019). On the contrary, studies that attempt to represent disk thermodynamics in a more realistic way (Marzari et al. 2012; Gyergyovits et al. 2014) find a dramatic reduction of disk eccentricity⁶, down to ⟨e_d⟩ ~ 0.02–0.05. As our results demonstrate, disk eccentricity at this level would allow planetesimal growth to proceed successfully starting from objects of a few to 10 km in size (see Fig. 1b). We note that in our model, ⟨e_d ⟩ = 2.5e₀.

Similarly, accounting for the effect of disk self-gravity on its dynamics has also been shown to suppress ⟨e_d ⟩ (Marzari et al. 2009). Another interesting finding is that S-type disks tend to be less eccentric in more eccentric binaries (Marzari et al. 2009, 2012; Regály et al. 2011), presumably because of their more compact size.

Also important for the possibility of efficient planetesimal growth is the orientation of the disk, since the dynamically quiet zone, greatly facilitating planetesimal growth, arises only in disks which are close to apsidal alignment (ϖ_d ≈ 0) with the binary orbit (see Fig. 1c). Unfortunately, at the moment there is little clarity regarding this issue. Marzari et al. (2009) find the disk to be precessing, that is, there is not a well-defined value for ϖ_d. In some other studies, the disk is aligned with the binary orbit (Müller & Kley 2012; Martin et al. 2020). Gyergyovits et al. (2014) find the disk to be close to alignment, with |ϖ_d |≈ 0.6 rad. Marzari et al. (2012) find anti-alignment (ϖ_d ≈ 180°), whereas Marzari et al. (2009) find either anti-alignment or precession, depending on whether self-gravity is included. Paardekooper et al. (2008) find either anti-alignment or precession depending on the flux limiter used in their code.

To summarize, the existing numerical studies of S-type disks in binaries have not yet arrived at a consensus regarding the values of either e₀ or ϖ_d, and their dependence on various physical inputs, for example, realistic thermodynamics, disk self-gravity, viscosity, etc. Given the importance of these characteristics for understanding planet formation in S-type binaries, we encourage future numerical efforts to focus on measuring the true value of ϖ_d as well as disk eccentricity e₀, using the most realistic physical inputs, such as disk thermodynamics and self-gravity.

7.4 Comparison of detailed coagulation-fragmentation calculations with simple heuristics

Full coagulation-fragmentation simulations like the ones featured in this work are the most accurate predictor of the outcome of planetesimal growth. However, they are numerically expensive. For that reason, in our previous work (e.g., Rafikov & Silsbee 2015b), we employed simple local (i.e., neglecting radial inspiral) heuristic criteria to determine conditions under which growth would occur.

The simplest of these is the assumption that growth occurs if and only if there are no catastrophically disruptive collisions. However, as shown in Sect. 5.1, it is also possible that erosion can prevent growth even in the absence of catastrophicdisruption. For that reason Rafikov & Silsbee (2015b) considered an alternative criterion, in which growth happens provided that no erosive collisions occur between planetesimals of “similar" sizes (i.e., bodies with a mass ratio greater than 10⁻²).

Now, armed with the full coagulation-simulation code, we can test the performance of each of these growth criteria in predicting the success of planetesimal evolution. To that effect we used our simulations in a single-annulus mode to calculate d_min for a variety of models featuring different values of d_c and v_c. As discussed in Rafikov & Silsbee (2015a), the collision rates and velocities are entirely determined by d_c and v_c (for the assumed level of random motions), regardless of what combination of disk parameters lead to the specific d_c and v_c (we set v_K as appropriate at 2 AU in the disk around γ Cephei). Therefore, by exploring the range of plausible v_c and d_c we check the validity of these heuristics for all values of the disk parameters in our model.

Figure 13 displays the results of this calculation (each panel corresponds to a different value of v_c), with the solid blue line showing d_min as a function of d_c, computed by our code using the procedure described in Sect. 5.2.1. The red line shows the minimum value of d_init such that no catastrophic disruption occurs for planetesimals of size d_init or larger. The green line shows the minimum value of d_init such that no body larger than d_init suffers an erosive collision with a partner with mass ratio greater than 1%. The same criteria were used to determine the gray and black regions, respectively, in Fig. 4 of Rafikov & Silsbee (2015b).

This figure shows that the criterion of no catastrophic disruption is too optimistic, typically predicting d_min between two and five times smaller (depending on d_c and v_c) than the one determined by the coagulation-fragmentation simulation. On the contrary, the criterion of no erosion between bodies with a mass ratio greater than 1% is too restrictive, with the green curve exceeding the blue curve by a factor two to five.

Interestingly, the blue curve is often roughly equally separated from both red and green ones (in log space), which motivated us to also show the dashed blue line computed as the geometric mean of the predictions from the “no catastrophic disruption” and “no erosion” models. As one can see, it matches the simulation-based d_min (blue curve) surprisingly well, within a factor of two or even better (although we do not claim this to have any deep meaning). Thus, the geometric mean of the predictions for d_min based on the two criteria used in Rafikov & Silsbee (2015b) can be employed as a rough approximation for predicting the success of planetesimal growth in the absence of sophisticated coagulation-fragmentation simulations.

Fig. 13 Minimum planetesimal size dmin resulting in the formation of large objects obtained with different methods for different values of vc (indicated inpanels) and dc. The solid blue line shows dmin computedusing our multi-annulus simulations. The green line is the analytic estimate for dmin assuming that growth can occur starting from the smallest size that does not suffer erosion by any particle greater than 1% of its mass, according to the Stewart & Leinhardt (2009) prescription. The red line is the analytic estimate for dmin assuming that growth can occur as long as catastrophic disruptions are absent. The dashed blue line shows the geometric mean of the red and green lines.

7.5 Comparison with previous work

A number of past studies of planet formation in S-type binaries have arrived at rather pessimistic conclusions regarding planetesimal growth starting from a population of 1–10 km objects (Thébault et al. 2006, 2008, 2009; Beaugé et al. 2010). In particular, Thebault (2011) found that even d_init ≳ 10² km planetesimals do not ensure planetesimal growth in HD 196885 system at the location of the known planet. All these studies included the effect of gas drag on planetesimal dynamics, and Paardekooper et al. (2008) and Beaugé et al. (2010) even considered the possibility of the disk being eccentric.

However, as noted already in Rafikov (2013b) and Rafikov & Silsbee (2015a,b), all these calculations ignored the effect of disk gravity, which dramatically changes planetesimal dynamics. For example, Beaugé et al. (2010) find low planetesimal collisional velocities in the outer disk, whereas we find a very hostile dynamical environment there (see Fig. 9). The existence of the dynamically quiet location in apsidally aligned disks substantially alleviates the fragmentation bottleneck, allowing in many cases planetesimal growth to proceed even starting at d_init ~ 1 km (see Fig. 10).

Another factor that distinguishes our study from others in terms of realism is our reliance on the global coagulation-fragmentation calculations fully accounting for the gas drag-driven inspiral of planetesimals. Whereas Thébault et al. (2006, 2008, 2009) and Rafikov & Silsbee (2015b) drew their conclusions based on individual collisional outcomes (see Sect. 7.4), we are able to model the evolution of the planetesimal size distribution in its entirety, fully accounting for the fine details of planetesimal dynamics atevery step.

Compared to Rafikov & Silsbee (2015b), we find less need for a massive (≳ 30M_J) protoplanetary disk: the lowest value of Σ₀ allowing growth to proceed starting at d_init ~ 1 km in Fig. 10a corresponds to a disk mass of just 2 M_J. This difference arises primarily from our current full treatment of the planetesimal size evolution, but also from focusing on planet formation not at a given location but globally, in the whole disk.

8 Conclusions

We present the first fully global coagulation-fragmentation framework for exploring planetesimal growth in S-type binary systems, such as α Cen and γ Cep. Our framework fully accounts for the specifics of planetesimal dynamics in binaries, including perturbations not only due to the eccentric stellar companion, but also due to the gravity of the eccentric protoplanetary disk in which planetesimals orbit. We include the effects of gas drag on the eccentricity dynamics of planetesimals, as well as on their radial inspiral. Thus, this framework brings studies of planet formation in binaries to an entirely new level of realism.

By running a suite of simulations with varied model inputs, we are able to delineate the conditions under which planetesimals can grow to large sizes (hundreds of kilometers) despite the adverse effects of fragmentation. The difficulty of planet formation for a given disk model was measured via the parameter d_min – the smallest size of seed planetesimals that allows their growth into the planetary regime – with the smaller d_min implying more favorable conditions for planet formation. Based on the results of these calculations, we are able to draw the following general conclusions.

• For most disk parameters considered in this paper, planet formation in binaries such as γ Cephei can successfully occur provided that the initial planetesimal size is ≳10 km; however, for favorable disk parameters, this minimum initial size can go down to ≲1 km;

• The gravitational effect of the protoplanetary disk plays the key role in lowering the minimum initial planetesimal size necessary for sustained growth by a factor of four. This reduction can be achieved in protoplanetary disks apsidally aligned with the binary, in which a dynamically quiet zone appears within the disk provided that the mass-weighted mean disk eccentricity ≲0.05 (corresponding, in the context of our model, to e₀ ≲ 0.02);

• We identify disk eccentricity (e₀) and the apsidal misalignment angle (ϖ_d) as the parameters to which planetesimal growth is most sensitive. Whenever the dynamically quiet location exists in the disk, growth is also strongly dependent on the level of random motion of planetesimals;

• For the models explored in this work, the outcome of planetesimal evolution is only weakly dependent on the disk surface density Σ₀;

• Radial inspiral of solid material appears to make little difference to the critical planetesimal size from which collisional agglomeration can occur;

• The planetesimal strength, as parameterized by the two models of Stewart & Leinhardt (2009), makes only about a factor of two difference to the critical planetesimal size;

• Planetesimal growth to large sizes can occur even in the face of some erosion, but in the environments with frequent collisions leading to catastrophic disruption, growth will not occur. We used our coagulation-fragmentation simulations to provide a new heuristic criterion for successful planetesimal evolution in Sect. 7.4;

• Our results provide indirect support for models in which planetesimals are born large (≳10 km), such as those based on streaming instability.

The general framework for the realistic treatment of planetesimal growth in binaries developed in this work can be applied to study a range of other problems: circumbinary planet formation in P-type systems (Silsbee & Rafikov 2015a), collisional evolution of debris disks perturbed by planetary (Nesvold & Kuchner 2015; Sefilian et al. 2021) or stellar (Thebault et al. 2021) companions,and so on.

Acknowledgements

The authors acknowledge helpful discussion with Lucas Jordan. K.S. acknowledges the support of the Max Planck Society. R.R.R. thanks NASA grant 15-XRP15-2-0139, STFC grant ST/T00049X/1, and John N. Bahcall Fellowship for financial support.

Appendix A Detailed description of the numerical methods

Here we describe a number of technical details of the implementation of our code, including both numerical details and physical inputs.

A.1 Mass redistribution due to collisions

In a given annulus, our code solves the standard coagulation-fragmentation problem, which effectively is a version of the coagulation equation of Smoluchowski (1916) modified to account for fragmentation (see Eq. (1) of Rafikov et al. 2020). Our framework employs mass space discretization in the form of logarithmically spaced mass bins with constant M_i+1∕M_i. The grid extends from the minimum mass corresponding to a 10 m object to a maximum mass corresponding to a 300 km object. Fragments with sizes below 10 m that form in planetesimal collisions are simply removed from the system (i.e., total mass is not conserved for this reason alone). This leads to certain numerical artifacts discussed in Sect. 5.1.

At each time step, the code calculates the mean collision rate R_ij between planetesimals in bin i and bin j for all i and j in [1, N_bins]. This calculation is described in Appendix A.4. The number of collisions Z_ij in a given time step of length Δt between bodies in mass bins i and j is then drawn randomly from a Poisson distribution with the mean R_ijΔt: $$Z_{ij}=\text{Poisson}(R_{ij}\Delta t).$$(A.1)

We define an operator $\widehat{F}$ such that F_ijk is the number of fragments inmass bin i produced by a collision of a particle in bin j with a particle in bin k (see Sect. A.2).

The mass distribution n is updated as $$n_i=n_i+{\Sigma }_{j=1}^{N_{\text{bins}}}{\Sigma }_{k=1}^{N_{\text{bins}}}{C}_{ijk},$$(A.2)

where $$C_{ijk}=Z_{jk}\frac{1+{\delta }_{jk}}{2}\left[F_{ijk}-{\delta }_{ij}-{\delta }_{ik}\right].$$(A.3)

This takes into account both addition of particles to bin i from collisional fragments (the F_ijk operator), as well as collisional destruction of particles in bin i (the − δ_ij and − δ_ik terms). The factor of (1 + δ_jk)∕2 is to prevent over-counting when summing over both j and k.

At times, due to a fluctuation, a number of particles in a particular mass bin may become negative. In this case we reset the number of particles in such a bin to zero; this effectively leads to addition of mass to the system. However, we checked several runs of our fiducial simulation, and found that relative to the total mass in the simulation, less than 1 part in 10⁵ was added in this way.

A.2 Collisional outcomes

This section describes the physics which goes into the determination of F_ijk. The outcome of a collision between a body in mass bin i and one in mass bin j is dependent on the collision velocity. High-velocity collisions will result in many small fragments, and no large remnant bodies. Low-velocity collisions will result in one body which contains nearly all of the pre-collision mass and only very small fragments. Calculation of collision velocity v_coll is described in Sect. A.3.

Given a value of v_coll, we use the recipe provided in Stewart & Leinhardt (2009) to determine the mass of the largest remnant M_lr. We take the mass of the largest remnant to be $$\frac{M_{\text{lr}}}{M_{\text{tot}}}=1-0.5\frac{Q_R}{Q_{\text{RD}}^{*}},$$(A.4)

where M_tot = m_i + m_j is the total mass of the two incoming bodies, Q_R is the specific center-of-mass energy of the collision given by $$Q_R=0.5\frac{m_i{m}_j{v}_{\text{coll}}^2}{M_{\text{tot}}},$$(A.5)

and $Q_{\text{RD}}^{*}$ is the critical value of Q_R which leads to catastrophic disruption (which corresponds to M_lr < 0.5M_tot); $Q_{\text{RD}}^{*}$ is a function of the sizes and composition of the planetesimals, and v_coll, which is given by Eq. (2) from Stewart & Leinhardt (2009). As our default, we used the parameters appropriate for their “strong rock” planetesimals; however we also ran simulations using their “rubble pile” planetesimals.

When M_lr is calculated, it will in general be between M_i and M_i+1 for some i. We then add a fraction P_i of the mass to bin i and P_i+1 to bin i + 1 given by $$P_i=\frac{M_{i+1}-M_{\text{lr}}}{M_{i+1}-M_i},P_{i+1}=\frac{M_{\text{lr}}-M_i}{M_{i+1}-M_i}.$$(A.6)

We assume that the rest of the fragments follow a mass distribution with a power law index ξ. The precise value of ξ does not affect the long-term outcome of fragmentation cascade (O’Brien & Greenberg 2003). The upper cutoff mass for this distribution is given by $$M_{\text{cutoff}}=\{\begin{array}{ll}0.5(M_{\text{tot}}-M_{\text{lr}})\hfill & M_{\text{lr}}/M_{\text{tot}}\ge 0.5,\hfill \\ \max (M_{\text{tot}}b,0.5M_{\text{lr}})\hfill & M_{\text{lr}}/M_{\text{tot}}<\mathrm{0.5.}\hfill \end{array}$$(A.7)

Here b is an adjustable parameter which determines the cutoff mass for catastrophic collisions: whenever Eq. (A.4) yields M_lr< 2bM_tot, we put all of the mass into fragments (as described below) without leaving a single largest remnant. This situation often arises in very energetic collisions when M_lr is small.

The fragment mass is assigned to mass bins with 0 < i ≤ i_cutoff, where i_cutoff is the index of the largest bin with mass less than M_cutoff. This is donesuch that $$F_{ijk}={{\rm Y}}_{jk}{M}_i^{1+\xi },$$(A.8)

where ϒ_jk is a normalization constant. The exponent of M_i is 1 + ξ rather than ξ because of the logarithmic spacing of the mass bins.

In order to solve for ϒ_jk, we note that $${\sum }_{i=-\infty }^{i_{\text{cutoff}}}{{\rm Y}}_{jk}{M}_i^{2+\xi }=M_{\text{tot}}-M_{\text{lr}},$$(A.9)

where nonpositive values of i correspond to (hypothetical) mass bins with M_i less than the minimum size considered in our simulations. Letting ℜ ≡ M_i+1∕M_i, we can then solve for ϒ_jk and write $$F_{ijk}=\frac{M_{\text{tot}}-M_{\text{lr}}}{M_{\text{cutoff}}}\times {(M_i/M_{\text{cutoff}})}^{1+\xi }\times \left(1-{\Re }^{-\xi -2}\right).$$(A.10)

Some mass is lost from the simulation because it is contained in bodies with mass less than M₁. This lost mass is equal to $$M_{jk}^{\text{lost}}={{\rm Y}}_{jk}{\sum }_{i=-\infty }^0{M}_i^{2+\xi }=(M_{\text{tot}}-M_{\text{lr}}){\left[\frac{M_0}{M_{\text{cut}}}\right]}^{2+\xi }.$$(A.11)

Because of the self-similar form of the fragment mass spectrum, we were able to implement this fragmentation algorithm with numerical cost scaling as $O(N_{\text{bins}}^2)$ using explicit time stepping (see Rafikov et al. 2020 for details). This considerably reduces the computational burden of our calculations.

A.3 Collision velocity

We calculate the collision velocity v_coll by considering two separate dynamical regimes. In the first, the relative motion of planetesimals is dominated by their secularly excited eccentricities. Equation (64) of RS15a then gives the relative eccentricity e_ij between planetesimals as a function of their masses. Given e_ij, there is a distribution of relative encounter velocities v_ij between v_min = 0.5e_ijv_K and v_max = e_ijv_K, where v_K is the local Keplerian speed. This distribution is given by Eq. (62) of RS15a as $$\frac{\text{d}{f}_{ij}}{\text{d}{v}_{ij}}=\frac{v_{\text{max}}^{-1}}{E(\sqrt{3}/2)}\frac{v_{ij}^2}{\left[(v_{\text{max}}^2-v_{ij}^2)(v_{ij}^2-v_{\text{min}}^2)\right]}.$$(A.12)

At each time step, for each set of collision partners, a random number λ between 0.5 and 1 is drawn based on the distribution in Eq. (A.12), and the collision velocity in this regime is taken to be $$v_{\text{coll}}=\lambda e_{ij}{v}_{\text{K}}.$$(A.13)

In the second dynamical regime, we assume that random motions dominate over the secularly excited eccentricities. We assume that planetesimals have eccentricity vectors (k, h) = (k_s + k_rand, h_s + h_rand) where (k_s, h_s) and (k_rand, h_rand) are the secular and random components, respectively; k_rand and h_rand are assumed to be drawn from independent Gaussian distributions $$\begin{array}{l}\rho (k_{\text{rand}})=\frac{1}{{\sigma }_e\sqrt{2\pi }}\exp \left[-\frac{{(k_{\text{rand}})}^2}{2{\sigma }_e^2}\right],\\ \rho (h_{\text{rand}})=\frac{1}{{\sigma }_e\sqrt{2\pi }}\exp \left[-\frac{{(h_{\text{rand}})}^2}{2{\sigma }_e^2}\right],\end{array}$$(A.14)

where we take σ_e = 2σ_i. This leads to a distribution of the relative random eccentricity $e_{\text{rand}}=\sqrt{{(k_{\text{rand},\text{1}}-k_{\text{rand},\text{2}})}^2+{(h_{\text{rand},\text{1}}-h_{\text{rand},\text{2}})}^2}$ given by theRayleigh distribution $$\rho (e_{\text{rand}})=\frac{e_{\text{rand}}}{2{\sigma }_e^2}\exp \left[-\frac{{(e_{\text{rand}})}^2}{4{\sigma }_e^2}\right].$$(A.15)

In principle,for each set of collision partners, we should then take the collision velocity to be given by $$v_{\text{coll}}=v_{\text{K}}\sqrt{{\left(\lambda e_{\text{ij}}\right)}^2+e_{\text{rand}}^2},$$(A.16)

where e_rand is drawn from the distribution in Eq. (A.15), and λ is drawn from the distribution in Eq. (A.12).

In practice, doing this at each time-step would be too time-consuming. Instead, we use the following algorithm. Let F(v) be the true collision velocity distribution for a given set of collision partners in a given annulus. We approximately calculate N₁ points in v-space, such that there is a constant fraction of the distribution F(v) between consecutive points: this means that the i^th point v_i satisfies the relation ${\displaystyle {\int }_0^{v_i}F}(v)\text{d}v=(i+0.5)/N_1$, for i ∈ [0, N₁ − 1]. We then calculate the collision rates and outcomes for each set of collision partners at the beginning of the simulation, using the N₁ different velocities. During the simulation, in each timestep, and for each set of collision partners, we randomly choose one of the N₁ collisional outcomes and rates to use for that timestep.

To calculate the N₁ values v_i, we initially create two vectors z_rand and z_λ of length N₂ with a constant fraction of the distribution for e_rand between each entry of z_rand and a constant fraction of the distribution of λ between each entry of z_λ. This means that ${\displaystyle {\int }_{z_{\text{rand}}^i}^{z_{\text{rand}}^{i+1}}\rho }(e_{\text{rand}})=1/N_2$, and a similar relation exists between z_λ and the distribution for λ implied by Eq. (A.12).

We then draw N₃ samples of the collision velocity using Eq. (A.16) and randomly selected entries from z_rand and z_λ. We then select the N₁ values of v_i so there are an equal number of the N₃ sample velocities between consecutive values of v_i. We take N₁ = 11, N₂ = 100, and N₃ = 990.

We ignore shear motion in calculating the collision velocity because if the shear motion is the dominant effect, then the planetesimals will completely merge in collisions, so there is no reason to consider shear motion except in how it affects the collision rate.

We note that the differential radial drift described in Sect. 3.5 is not included in our calculation of v_coll as it does not typically contribute substantially to the collision velocities between bodies of the sizes that we consider in this paper.

A.4 Collision rate

Next we describe our calculation of the collision rate. In the first dynamical regime (where the relative motion of planetesimals is dominated by their secularly excited eccentricities), we use Eq. (63) from RS15a (corrected so that e_max is replaced with v_max∕a_p). This gives a flux F_ij of planetesimals from bin j past a planetesimalin bin i of $$F_{ij}=\frac{E(\sqrt{3}/2){\Sigma }_j{v}_{\text{max}}}{{\pi }^2{a}_{\text{p}}{m}_j{I}_{ij}},$$(A.17)

where I_ij is the relative inclination between the two groups of particles. We assume that all of the planetesimal inclination is due to random motions. We assume all planetesimals have inclination vectors (q, p) with q and p drawn from independent Gaussian distributions: $$\rho (q)=\frac{1}{{\sigma }_i\sqrt{2\pi }}\exp \left(-\frac{q^2}{2{\sigma }_i^2}\right);\rho (p)=\frac{1}{{\sigma }_i\sqrt{2\pi }}\exp \left(-\frac{p^2}{2{\sigma }_i^2}\right).$$(A.18)

The distribution of $I_{\text{ij}}=\sqrt{{(q_i-q_j)}^2+{(p_i-p_j)}^2}$ is therefore given [analogous to Eq. (A.15)] by $$\rho (I_{ij})=\frac{1}{2{\sigma }_i^2}{I}_{ij}\exp \left(-\frac{I_{ij}^2}{4{\sigma }_i^2}\right).$$(A.19)

Then, averaged over I_ij, the average value of F_ij is given by $$\langle F_{ij}\rangle ={\displaystyle {\int }_{I_{ij}=0}^{\infty }\rho }(I_{ij})F_{ij}d{I}_{ij}=\frac{E(\sqrt{3}/2){\Sigma }_j{v}_{\text{max}}}{2{\pi }^{3/2}{\sigma }_i{m}_j{a}_{\text{p}}}.$$(A.20)

The collision rate for a given particle in bin i with the population in bin j is given by R_ij = ⟨F_ij⟩σ_ij, where σ_ij is the collision cross section including gravitational focusing: $${\sigma }_{ij}=A_{\text{g}}\left[1+{\left(\frac{v_{\text{esc}}}{\lambda e_{ij}{v}_{\text{K}}}\right)}^2\right],$$(A.21)

where $A_{\text{g}}=\pi {(d_i+d_j)}^2$ is the geometric cross section, and $v_{\text{esc}}=\sqrt{2G(m_i+m_j)/(d_i+d_j)}$. As discussed in Sect. A.3, at each time-step, we assume a particular value of the collision velocity drawn from the distribution F(v). For the purposes of calculating the collision rate, we assume λ to be in the same quantile of the distribution for λ as v is for F(v).

In the second dynamical regime of dominant random motions ($\sqrt{\pi }{\sigma }_e>e_{ij}$) we use the two-body accretion formula from Eq. (17) of Greenzweig & Lissauer (1992): $$R_{ij}=\frac{{\Sigma }_j{A}_{\text{g}}}{2{\pi }^2{m}_j}\left[\mathcal{F}(I_{*})+\frac{v_{\text{esc}}^2}{v_{\text{K}}^2}\frac{\mathcal{G}(I_{*})}{e_{*}^2}\right].$$(A.22)

Here, $\mathcal{F}$ and $\mathcal{G}$ are functions of I_* = σ_i∕σ_e. We assume that I_* = 1∕2, in which case $\mathcal{F}(1/2)=17.3$, and $\mathcal{G}(1/2)=38.2$. e_* is the rms random eccentricity. Greenzweig & Lissauer (1992) considered a planetesimal accreting on a circular orbit, but in our case, we should use the rms relative random eccentricity which can be shown from Eq. (A.15) to be e_* = 2σ_e. Additionally, we must divide by the orbital period as Eq. (A.22) is in units where the orbital period is unity: $$R_{ij}=\frac{{\Sigma }_j{A}_{\text{g}}{v}_{\text{K}}}{4{\pi }^3{a}_{\text{p}}{m}_j}\left[\mathcal{F}(I_{*})+\frac{v_{\text{esc}}^2}{v_{\text{K}}^2}\frac{\mathcal{G}(I_{*})}{e_{*}^2}\right].$$(A.23)

We approximate the transition to the shear-dominated regime by replacing e_* in Eq. (A.23) with e_* + c₁e_H, and then multiplying the entire right hand side by 1 + c₂(e_H∕e_*), where c₁ and c₂ are interpolation constants. This gives us a formula $$R_{ij}=\frac{{\Sigma }_j{A}_{\text{g}}{v}_{\text{K}}}{4{\pi }^3{a}_{\text{p}}{m}_j}\left[\mathcal{F}(I_{*})+\frac{v_{\text{esc}}^2}{v_{\text{K}}^2}\frac{\mathcal{G}(I_{*})}{{(e_{*}+c_1{e}_{\text{H}})}^2}\right]\left(1+c_2\frac{e_{\text{H}}}{e_{*}}\right).$$(A.24)

We estimatec₁ and c₂ by comparing the rates given by Eq. (A.24) with those shown in Fig. 5 of Greenzweig & Lissauer (1992). From these data we estimate c₁ = 4.75, c₂ = 22.6, which we use in our work.

We expect a transition between the two dynamical regimes to occur when $e_{\text{ij}}=\langle e_{\text{r}}\rangle =\sqrt{\pi }{{\sigma }^{\prime }}_e$, where ${{\sigma }^{\prime }}_e={\sigma }_e+c_1{e}_{\text{H}}/2$ (the factor of 2 arising from the fact that e_* = 2σ_e). The collision rates in the two regimes are not exactly equal at the transition point. Let R₁ be the collision rate in the low secular velocity limit given by Eq. (A.24), and let R₂ be the collision rate when the secular velocity dominates; R₂ is given by $$R_2=\frac{E(\sqrt{3}/2){\Sigma }_j{A}_{\text{g}}{v}_{\text{K}}{e}_{ij}}{2{\pi }^{3/2}{\sigma }_i{a}_{\text{p}}{m}_{\text{p}}}\left[1+{\left(\frac{v_{\text{esc}}}{v_{ij}}\right)}^2\right].$$(A.25)

Assuming σ_i∕σ_e = 1∕2, then at $e_{ij}=\sqrt{\pi }{\sigma }_e$, ignoring shear, R₂∕R₁ = 2.8 in the limit that v_esc∕(σ_ev_K) ≪ 1, and 1.6∕λ² in the limit that v_esc∕(σ_ev_K) ≫ 1. To make the collision rate continuous, we approximate it as $$R_{ij}(v_{ij})=f{R}_1+(1-f)R_2,$$(A.26)

with $f=1/(1+{[e_{ij}/(\sqrt{\pi }{{\sigma }^{\prime }}_e)]}^2)=1/(1+[v_{ij}/{(\sqrt{\pi }\lambda {{\sigma }^{\prime }}_e{v}_{\text{K}})}^2])$. v_ij = λe_ijv_K is the relative secular velocity between bodies in mass-bins i and j for that time step.

The above prescription needs to be modified in the case of collisions between two bodies within the same mass bin, where the focusing factor can vary significantly depending on the location within the mass bin of the two colliding partners. For example, two bodies with masses right at the center of the mass bin will have a higher focusing factor than one body at the bottom end of the mass bin and one body at the top end; in the latter case, the approach velocity will be higher.

Let v_max,i be the relative (secularly excited) velocity between the heaviest body in the ith mass bin and the lightest. Let x₁ = (m₁− m_bottom)∕(m_top − m_bottom) be the position of body 1 in the mass bin, where m_top and m_bottom are the masses of the heaviest and lightest bodies in the bin. We define x₂ analogously for the second body. The relative secularly excited velocity v₁₂ between particles 1 and 2 is then equalto v_max|x₁ − x₂| and has a distribution (assuming x₁ and x₂ to be uniformly distributed in [0, 1]) $$\rho (v_{12})=\frac{2}{v_{\text{max}}}\left(1-\frac{v_{12}}{v_{\text{max}}}\right).$$(A.27)

Therefore, ⟨R_ii⟩, the rate at which a particle in bin i experiences collisions with other particles in bin i, is given by $$\langle R_{ii}\rangle ={\displaystyle {\int }_0^{v_{\text{max}}}\rho }(v_{12})R_{ii}(v_{12})\text{d}{v}_{12}.$$(A.28)

A.5 Radial inspiral

Using the results of Rafikov & Silsbee (2015a,b) we obtain the radial drift velocity of particles as a function of their size and semimajor axis (see Sect. 3.5 for details). Let Δr_i be the radial width of annulus i, and let the drift rate of a particle in annulus i and mass bin j be v_i,j. We then let the change in the number of particles in annulus i and mass bin j in a time interval Δt_insp be given by $$\Delta n_{i,j}=-\frac{n_{i,j}{v}_{i,j}\Delta t_{\text{insp}}}{\Delta r_i}+\frac{n_{i+1,j}{v}_{i+1,j}\Delta t_{\text{insp}}}{\Delta r_{i+1}},$$(A.29)

We use N_ann = 60 annuli. These were spaced so as to minimize changes in e_c or a within one bin. The exact bin boundaries depend on the run of e_c(a) which varies between our different disk models. For the fiducial model, a changed by less than 16% between adjacent bins and e_c by less than 41% in the cases where e_c > 10⁻⁴.

A.6 Time step determination

In this section, we describe how the time steps for the coagulation process and the radial inspiral are determined.

A.6.1 Coagulation-fragmentation time step

The time step for the coagulation-fragmentation process must be short enough that the distribution does not change appreciably during the length of a time step. In principle, this means that we would like $$\frac{|n_i(t)-n_i(t+\Delta t_{\text{coag}})|}{n_i(t)}<ϵ$$(A.30)

for all i, where ϵ is some small number. In practice, this is complicated by the fact that mass is added to bins in discrete quantities, and that some bins have a small number of particles. In a bin with a small number of particles, even adding one additional particle will make a large fractional change to the amount of mass in that bin. Additionally, collisions are a random process, so a fluctuation can change the mass in a bin by more than expected. To address these concerns, we instead pick the largest time step Δt_coag such that for all i $$\begin{array}{l}\left|\langle \frac{n_i(t)-n_i(t+\Delta t_{\text{coag})}}{n_i(t)}\rangle \right|<{ϵ}_1\text{or}\\ \left|\langle \frac{m_i(n_i(t)-n_i(t+\Delta t_{\text{coag}}))}{{\Sigma }_{i=1}^{N_{\text{bins}}}{m}_i{n}_i}\rangle \right|<{ϵ}_2,\end{array}$$(A.31)

where the angle brackets refer to averaging over the distribution of encounter velocity in Eq. (A.12). To save computational time, we approximate the average over the encounter velocity distribution by evaluating the collisional outcomes using the median encounter velocity. Here, as before, a time step is acceptable if for each bin, either the fractional decrease in mass is less than ϵ₁, or the decrease in mass is less than ϵ₂ times the total mass in the system. Time step Δt_coag is recalculated every time we update the mass distribution, separately for each annulus, and the minimum value among all the annuli is used for each time step. It is not obvious from first principles which values of ϵ₁ and ϵ₂ are appropriate. But rows 10 and 11 of Table C.1 (which explores sensitivity of our calculations to the numerical parameters used; see Appendix C) show that variation of these values by a factor of two (or more) from our fiducial values of 0.05 and 10⁻⁶ have a negligible effect on d_min.

In the simulation used to create Fig. 4, we used a lower value of ϵ₂ so that the distributions would remain smooth, even in regions of mass-space containing a negligible amount of mass. The success of the coagulation process, and the shape of the distribution at sizes greater than d_init in those simulations is not affected by this change.

A.6.2 Radial inspiral time step

The code uses a different time step for the radial drift. Mass transport due to radial drift is recalculated using Eq. (A.29) every interval of time Δt_insp such that $$\Delta t_{\text{insp}}\times \max \left(\frac{v_{\text{drift}}}{\Delta r_i}\right)={ϵ}_{\text{drift}},$$(A.32)

which guarantees that the inspiraling particles do not cross the full radial bin width in time Δt_insp (v_drift is the inspiral speed). While we see from Eq. (A.29) that ϵ_drift must not be larger than 1, it is apparent from the results of Sect. B.3 below that there is little advantage to picking ϵ_drift less than 1. Therefore, we use ϵ_drift = 1.

Appendix B Tests of the code

In this section, we describe the tests we did of the code against known analytic solutions in certain limiting cases. We independently verify three code functionalities – coagulation, fragmentation, and radial inspiral – and find that the code does a reasonable job of reproducing known solutions to the accuracy warranted given the imperfect knowledge of the physical processes involved.

Fig. B.1 Test of our coagulation module: comparison of simulation and analytic results for different coagulation kernels starting with N = 1012 particles of unit mass. Here η = Nt is a time coordinate, rescaled by the number of particles initially in the system. Each panel is made assuming a different coagulation kernel Aij, labeled in the panels and discussed in Sect. B.1.

B.1 Tests against known solutions of the coagulation equation

In this section we test the coagulation module of our code. We describe the results of simulations of the coagulation alone, without fragmentation, and with simplified collision rates and initial conditions that have analytic solutions we can compare with. There are three known analytic solutions to the Smoluchowski coagulation equation $$\frac{\text{d}{n}_k}{\text{d}t}=\frac{1}{2}{\sum }_{i+j=k}{A}_{ij}{n}_i{n}_j-n_k{\sum }_{i=1}^{\infty }{A}_{ik}{n}_i,$$(B.1)

obtained for kernels A_ij = 1 (Smoluchowski 1916), A_ij = i + j and A_ij = ij (Trubnikov 1971), with the initial condition of N particles of a single mass bin. Our code successfully reproduces these solutions with small deviations as discussed below.

The top panel of Fig. B.1 shows the number of particles per unit mass as a function of mass for three different times for the kernel A_ij = 1. In this figure we have used a mass bin spacing M_i+1∕M_i = 1.15. Initially there were N = 10¹² particles all of mass 1. Because the distribution evolves more rapidly if there are more particles, we express time in terms of η = Nt. In the limit of large N, the shape of the distribution depends only on η, rather than on N or t independently. The large scatter at low masses is due to the logarithmic spacing of the mass bins and the initial condition that all the mass is in bodies with unit mass. Since we are assuming perfect sticking in collisions for this exercise, only mass bins containing integer masses will contain a nonzero number of particles, leading to an uneven distribution at low masses. Other than this discrepancy at low mass, there is nearly perfect agreement between our code and the analytic solution.

The middle panel of Fig. B.1 shows the number of particles per unit mass as a function of mass for three different values of η for the kernel A_ij = i + j. The distribution evolves more rapidly than in the case of the constant kernel because of the higher collision rates. In this case, in addition to the scatter at low masses, we see that the numerical solution is slightly behind the analytic solution at the high-mass end of the distribution.

The bottom panel of Fig. B.1 shows the number of particles per unit mass as a function of mass for three different values of η for the kernel A_ij = ij. Here, the system evolves in a more complicated manner than the previous two cases (Trubnikov 1971). For η < 1, orderly growth from smaller to larger particles progresses as in the previous cases. However, for η > 1, one body becomes much larger than the others, and eventually consumes all the mass in the system. Our simulation is not well-equipped to handle this runaway growth phase, so we restrict ourselves to reproducing the analytic solutions for η < 1.

The discrepancy between the numerical and analytic solutions at η = 0.99 looks somewhat alarming. This is because the distribution function is evolving quite rapidly at this time. We also plotted the analytic solution for η = 0.997, and we see that it matches the numerical solution well, implying that this difference only corresponds to a time lag of just 0.3% in the numerical solution, which is not critical for our calculations.

B.2 Tests of the fragmentation module

Next we checked the ability of our code to properly handle fragmentation. Unlike the coagulation problem, there are no known analytical solutions describing time-dependent behavior in fragmenting systems. However, there are some knownresults for the steady-state fragmentation cascades to which collisionally evolving systems tend to converge. Inthis section we use these results to verify that the fragmentation module of our code is operating correctly.

O’Brien & Greenberg (2003) present a model of collisional fragmentation that leads to a steady state size distribution given by a power-law. Their model has one free parameter s, which parametrizes the strength of particles of size d_p in catastrophic disruption events as $$Q_{\text{RD}}^{*}=Q_0{d}_{\text{p}}^s.$$(B.2)

They show that the distribution of debris sizes D_f in the steady-state collisional cascade is given by $$\frac{\text{d}{n}_{\text{f}}(D_{\text{f}})}{\text{d}{D}_{\text{f}}}=B{D}_{\text{f}}^{-y},$$(B.3)

where n_f is the number of objects of size D_f (per unit D_f), and B is an overall scaling constant. The power-law index of the size distribution of the steady-state collisional cascade is given by $$y=\frac{21+s}{6+s}.$$(B.4)

O’Brien & Greenberg (2003) have assumed that the collisional cross section between two bodies is proportional to the square of the size of the target body (i.e., the geometric cross section). If we relax this assumption, as was done in Tanaka et al. (1996), and instead let the collision cross section be proportional to $d_{\text{p}}^{\alpha }$, one can show that y gets modified to $$y=\frac{15+3\alpha +s}{6+s}.$$(B.5)

Fig. B.2 Test of our fragmentation module: comparison of our simulation results (blue) with the theoretical predictions of O’Brien & Greenberg (2003, green; they fall on top of the simulation outputs in the top panels). The upper panel shows the number of particles per mass bin as a function of mass. The bottom panel shows the power-law slope of this mass distribution. The rightmost 40% of the mass bins are held fixed, so as not to introduce errors due to the finite extent of our simulation in mass space and time. The values of α, s, and z are labeled on the panels. See Sect. B.2 for details.

Equation (B.5) was derived assuming a steady-state mass distribution with no upper or lower cut-off. Since we are only able to simulate a finite range in masses, we must have a way of supplying mass to the system to achieve a steady state; otherwise, all mass will eventually be contained in particles smaller than the lower cutoff size for the simulation. We do this by fixing the counts in the top 40% of the mass bins, and letting the collisional cascade provide the mass to the mass bins corresponding to smaller particles. By doing this we continuously resupply mass to the large mass bins.

Figure B.2 shows a comparison of the steady-state result of the simulation against the expected power law distribution from Eq. (B.5) for three different values of s and α. The simulations were run with bin spacing M_i+1∕M_i = 1.15. In each case, the top panel shows the number of bodies in each bin as a function of bin mass. The blue line shows the particle numbers from the simulation, while the green line shows the expected number counts based on Eq. (B.5). The bottom panels show the logarithmic slope d ln n_f∕d ln m as a function of mass.

Panels a and b show the case where s = 0 and α = 2 (Dohnanyi 1969). In this case Eq. (B.5) predicts y = 7∕2. Assuming constant bulk density, so that $m\propto d_{\text{p}}^3$, we can show that $$\frac{\text{d}{n}_{\text{f}}}{\text{d}m}=B^{\prime }{m}^{-\frac{y+2}{3}},$$(B.6)

where B′ is the appropriate normalization. The number of particles contained in the bin with bin mass M_i scales as $\Delta n(M_i)~M_i^{-z}$, with z = (y − 1)∕3 = (3 + α)∕(6 + s), due to the logarithmic spacing of the mass bins. For s = 0 and α = 2, z = 5∕6.

Fig. B.3 Test of the radial inspiral module: a comparison of the evolution of the surface density profile due to radial inspiral in the simulation (dotted lines) vs. the analytic result (solid lines) at different moments of time for the inspiral rate given by Eq. (B.7). Panels a, c, and e correspond to the cases with Nann = 30, 60, and 240 radial annuli between 1 and 4 AU, and ϵdrift = 1.0. The bottom panels, b, d, and f, have the same Nann values, but ϵdrift = 0.3.

Panels c and d are the same as a and b except that s = −1.5 and α = 0, leading to a slope d ln n_f∕d ln m = −z = −2∕3. The slope is reproduced faithfully except right near the lower end of the mass distribution. This time, there is a regular small oscillation of the slope that does not go away with decreasing time step – a behavior that often emerges in fragmentation simulations (Belyaev & Rafikov 2011).

Finally, in panels e and f we let s = 0 and α = 1. This leads to the expectation that z = 2∕3.

In these three plots, we see that our simulation has nearly perfect agreement with the expected power-laws, except for understandable deviations at specific points. Objects at the small end of the mass distribution, have nothing smaller than themselves in the simulation, so they are destroyed less frequently than the model would predict, hence the steepening of the power-law. Then the over-abundance of small masses creates an under-abundance of objects of slightly larger masses, giving rise to the wave-like pattern that we observe. This effect is reduced in panels c and d because in this example we tuned the Q₀ in Eq. (B.2) so that objects at the lower end of the simulated mass range would not be destroyed by smaller particles, even if they were in the simulation. The wiggle at the lowest masses in this plot is due to the finite spacing of the mass bins.

B.3 Tests of the radial inspiral module

In this section, we describe our test of the code module handling the size-dependent radial inspiral. We made this test for a simple inspiral rate given by $$v_{\text{drift}}=v_0{\left(\frac{a}{a_0}\right)}^2.$$(B.7)

Let 𝔫(a, t) be the number of particles per unit semimajor axis at time t. Under the influence of inspiral, 𝔫 solves the equation $$\frac{\partial \mathfrak{n}(a,t)}{\partial t}-\frac{\partial }{\partial a}\left[v_{\text{drift}}(a)\mathfrak{n}(a,t)\right]=0.$$(B.8)

We note that the minus sign in Eq. (B.8) arises from the fact that we adopt a sign convention where positive v₀ corresponds to inward motion, toward smaller a.

Adopting the initial planetesimal distribution in the disk in the form $$\mathfrak{n}(a,t=0)={\mathfrak{n}}_0\left(4-\frac{a}{a_0}\right),$$(B.9)

for 0 < a < 4a₀, this leads to an evolution of the number density given by $$\mathfrak{n}(a,t)=\{\begin{array}{ll}{\mathfrak{n}}_0{s}^{-2}\left(4-a/(a_{0}s)\right)\hfill & a<4a_0/(1+4v_{0}t/a_0),\hfill \\ 0\hfill & a>4a_0/(1+4v_{0}t/a_0),\hfill \end{array}$$(B.10)

where $s=1-at{v}_0/a_0^2$.

Figure B.3 displays the evolution of the planetesimal radial distribution obtained by our code, in comparison with the analytical solution (B.10). Panels a, c, and e correspond to the cases with N_ann = 30, 60, and 240 radial annuli between 1 and 4 AU. The annuli are located at the same values of a as in Simulations 5, 1, and 7, respectively, and we take ϵ_drift = 1.0. The corresponding bottom panels b, d, and f are for ϵ_drift = 0.3.

We see that generally our numerical solutions stay within several percent of the analytic solutions, except the locations where the analytic solution tends to zero. The simulation does not do a good job of resolving the sharp outer edge of the distribution because particles are effectively redistributed through each annulus every time-step, leading to substantial numerical diffusion. The time step (i.e., the choice of ϵ_drift) makes little difference to the results, but the agreement is a strong function of the number of radial annuli N_ann used in our calculations (see Fig. B.3 and Table C.1).

Appendix C Sensitivity to numerical parameters

Our code employs a number of numerical parameters to set various time steps, properly resolve mass and inspiral-driven evolution, and so on. We varied these numerical code inputs to see the effect on the calculation outcomes. Table C.1 shows the results of this exercise, in which we use d_min and $d_{\text{min}}^{\text{rubble}}$ as our sensitivity metrics.

Simulation 1 uses the fiducial set of numerical parameters that we employ in all our production runs. As in Table 1, in all other runs the parameter that is different from Simulation 1 is highlighted in bold in each row. We see that d_min does not depend strongly on the numerical parameters, which justifies our particular choice of their values.

References

All Tables

All Figures

Fig. 1 Radial profiles of characteristic velocity vc (a) and planetesimal size dc (b) for different disk models, in which we vary disk eccentricity e0, surface density normalization Σ0, and apsidal angle ϖd (one parameter at a time relative to the fiducial model described in Sect. 5 and shown with the thick black curve), as indicated in the legend. See the text for more details on various features of these curves.

In the text

Fig. 2 Collision velocities between two planetesimals as a function of their sizes d1 and d2 in a disk with the fiducial set of parameters given in Sect. 5 and σi = 10−4. Each panel corresponds to a different location in the disk, which leads to the varying values of vc and dc labeled on the panels. Regions of erosive collisions are outlined in dotted white lines, and regions of catastrophic disruption are enclosed by solid white lines. The horizontal and vertical yellow lines correspond to the initial planetesimal size dinit = 4.3 km used in the simulation described in Sect. 5.1. The point (dc, dc) is shown as a red dot.

In the text

	Fig. 3 Inspiral speed as a function of planetesimal size at different locations in the disk (a) (see legend), and as a function of location in the disk for different planetesimal sizes (b) (see legend). The figure was made assuming the disk model described in Sect. 5, with e0 = 0.01, Σ0 = 1000 g cm−2, and ϖd= 0.
In the text

Fig. 4 Evolution of the size distribution as a function of time in a set of single-annulus simulations, described in Sect. 5.1, which should be consulted for details. The three panels correspond to environments with the different values of the critical velocity vc and critical size dc, previously shown in Fig. 2, at different locations (as labeled) in the disk around the primary star in the γ Cephei system. The vertical dotted lines correspond to the critical size dc, and the vertical dashed lines correspond to the initial size of the planetesimals dinit.

In the text

	Fig. 5 Time evolution of the size of the biggest object in the system, shown for the three simulations depicted in Fig. 4 and labeled according to their semimajor axis.
In the text

	Fig. 6 Evolution of the mass distribution in the entire disk in our fiducial multi-zone model (Simulation 1). Different panels are produced for different initial planetesimal sizes dinit in terms of dmin (as labeled on the plot), which is dmin = 1.6 km for this simulation.
In the text

	Fig. 7 Mass in planetesimals of all sizes (tracked in our mass grid) in the disk per unit semimajor axis dM∕da = 2πaΣ(a) as a functionof a in our fiducial model at different moments of time (see color bar). The different panels correspond to the same initial planetesimalsizes dinit as in Fig. 6. The dashed violet horizontal line corresponds to the mass distribution at t = 0.
In the text

	Fig. 8 Similar to Fig. 7 but showing the size of the largest body as a function of the semimajor axis. The dashed violet horizontal line now shows the size of seed planetesimals dinit.
In the text

Fig. 9 Collision velocity between d1 = 2 and d2 = 5 km planetesimals as a function of the semimajor axis a for the different values of model parameters varied one at a time (with respect to the fiducial set of parameters) and shown in the subplot legends. The fiducial value of the varied parameter in each subplot is shown with the solid black curve. The dotted line in panel (a) shows the collision velocities in the limit that σi = 0.

In the text

Fig. 10 Dependence of our key metric dmin – a minimum size of seed planetesimals that results in their eventual growth beyond 300 km – on some parameters of our simulations: (a) surface density normalization Σ0, (b) disk eccentricity normalization e0, (c) apsidal angle of the disk ϖd, (d) amplitude of the random velocity of planetesimals σi, and (e) parameter χir (artificially) regulating the magnitude of planetesimal inspiral. Solid lines and diamond points correspond to dmin, whereas dashed lines and circles are for $d_{\text{min}}^{\text{rubble}}$ – the minimum size for collisionally weak planetesimals. The causes of the trends shown here are discussed in Sects. 6.1–6.7. In each panel, the fiducial simulation is highlighted in red.

In the text

Fig. 11 Evolution of the size of the largest object at different locations in the disk for two models used to illustrate certain aspects of our calculations. (a) Fiducial model with gas drag turned off, χir = 0, computed for dinit∕dmin = 2, to be compared with Fig. 8c (see Sect. 6.7 for details). (b) Fiducial model computed for dinit∕dmin just slightly exceeding unity, to be compared with Fig. 8c,d (see Sect. 7.1 for details).

In the text

	Fig. 12 Location adq of the dynamically quiet zone in an apsidally aligned disk. We assume the fiducial disk plus binary parameters, except for Σ0 and e0 which are varied as shown on the plot. The fiducial model is marked by an orange star along the curve corresponding to Σ0 = 1000 g cm−2.
In the text

Fig. 13 Minimum planetesimal size dmin resulting in the formation of large objects obtained with different methods for different values of vc (indicated inpanels) and dc. The solid blue line shows dmin computedusing our multi-annulus simulations. The green line is the analytic estimate for dmin assuming that growth can occur starting from the smallest size that does not suffer erosion by any particle greater than 1% of its mass, according to the Stewart & Leinhardt (2009) prescription. The red line is the analytic estimate for dmin assuming that growth can occur as long as catastrophic disruptions are absent. The dashed blue line shows the geometric mean of the red and green lines.

In the text

	Fig. B.1 Test of our coagulation module: comparison of simulation and analytic results for different coagulation kernels starting with N = 1012 particles of unit mass. Here η = Nt is a time coordinate, rescaled by the number of particles initially in the system. Each panel is made assuming a different coagulation kernel Aij, labeled in the panels and discussed in Sect. B.1.
In the text

Fig. B.2 Test of our fragmentation module: comparison of our simulation results (blue) with the theoretical predictions of O’Brien & Greenberg (2003, green; they fall on top of the simulation outputs in the top panels). The upper panel shows the number of particles per mass bin as a function of mass. The bottom panel shows the power-law slope of this mass distribution. The rightmost 40% of the mass bins are held fixed, so as not to introduce errors due to the finite extent of our simulation in mass space and time. The values of α, s, and z are labeled on the panels. See Sect. B.2 for details.

In the text

Fig. B.3 Test of the radial inspiral module: a comparison of the evolution of the surface density profile due to radial inspiral in the simulation (dotted lines) vs. the analytic result (solid lines) at different moments of time for the inspiral rate given by Eq. (B.7). Panels a, c, and e correspond to the cases with Nann = 30, 60, and 240 radial annuli between 1 and 4 AU, and ϵdrift = 1.0. The bottom panels, b, d, and f, have the same Nann values, but ϵdrift = 0.3.

In the text