Issue 
A&A
Volume 652, August 2021



Article Number  A104  
Number of page(s)  27  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202141139  
Published online  17 August 2021 
Planet formation in stellar binaries: global simulations of planetesimal growth
^{1}
MaxPlanckInstitut für Extraterrestrische Physik,
85748
Garching,
Germany
email: ksilsbee@mpe.mpg.de
^{2}
Department of Applied Mathematics and Theoretical Physics, University of Cambridge,
Cambridge
CB3 0WA,
UK
^{3}
Institute for Advanced Study,
Einstein Drive,
Princeton,
NJ
08540,
USA
Received:
20
April
2021
Accepted:
21
June
2021
Planet formation around one component of a tight, eccentric binary system such as γ Cephei (with semimajor axis around 20 AU) is theoretically challenging because of destructive highvelocity collisions between planetesimals. Despite this fragmentation barrier, planets are known to exist in such (socalled Stype) orbital configurations. Here we present a novel numerical framework for carrying out multiannulus coagulationfragmentation calculations of planetesimal growth, which fully accounts for the specifics of planetesimal dynamics in binaries, details of planetesimal collision outcomes, and the radial transport of solids in the disk due to the gas dragdriven inspiral. Our dynamical inputs properly incorporate the gravitational effects of both the eccentric stellar companion and the massive nonaxisymmetric protoplanetary disk in which planetesimals reside, as well as gas drag. We identify a set of disk parameters that lead to successful planetesimal growth in systems such as γ Cephei or α Centauri starting from 1 to 10 km size objects. We identify the apsidal alignment of a protoplanetary disk with the binary orbit as one of the critical conditions for successful planetesimal growth: It naturally leads to the emergence of a dynamically quiet location in the disk (as long as the disk eccentricity is of order several percent), where favorable conditions for planetesimal growth exist. Accounting for the gravitational effect of a protoplanetary disk plays a key role in arriving at this conclusion, in agreement with our previous results. These findings lend support to the streaming instability as the mechanism of planetesimal formation. They provide important insights for theories of planet formation around both binary and single stars, as well as for the hydrodynamic simulations of protoplanetary disks in binaries (for which we identify a set of key diagnostics to verify).
Key words: planets and satellites: formation / protoplanetary disks / planetdisk interactions / planetstar interactions / planets and satellites: dynamical evolution and stability / binaries: general
© K. Silsbee and R. R. Rafikov 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Open Access funding provided by Max Planck Society.
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 planethosting systems as important test beds of planet formation theory.
In particular, tight eccentric binary systems that host planets orbiting one of the stellar companions (socalled Stype 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^{−1} 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 sizedependent alignment of planetesimal orbits (Thébault et al. 2008; Rafikov & Silsbee 2015a) gives rise to a dynamical environment with a mixture of low and highvelocity 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 pressuresupported and orbits ~ 30 m s^{−1} 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 particlegas 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 Stype 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 multiannulus coagulationfragmentation 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 Stype 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 multiannulus (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 Stype 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 Stype 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),\text{\hspace{1em}}{e}_{\text{g}}(a)={e}_{0}{\left(\frac{{a}_{0}}{a}\right)}^{1},\text{\hspace{1em}}a<{a}_{\text{out}},$$(1)
where a_{0} 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 solidtogas 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 nonaxisymmetric, 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 nonaxisymmetry 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 selfconsistently 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 socalled 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 Nbody integration of a modest number (~ 10^{3}) 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 coagulationfragmentation 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, nonaxisymmetric protoplanetary disk (Sect. 3.1); (ii) a unique sizedependent 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 Stype binaries
The motion of planetesimals around the primary in Stype 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 (timedependent) potential of the secondary over its orbit. However, since the gaseous protoplanetary disks in Stype 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 nonaxisymmetric, 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\mathrm{cos}\varpi +{B}_{\text{d}}e\mathrm{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 nonaxisymmetric 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 Stype 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^{−3} 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 shortperiod perturbations which act on the orbital period of the planetesimals (Murray & Dermott 1999). The magnitude of the shortperiod 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_{r}a_{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_{r}v_{K} in the absence of such shortperiod perturbations. Over longer timescales such shortperiod 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 gravityrelated 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 Stype 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^{1} τ_{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^{−1}, 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} $$\leftA{\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 planetesimalgas 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_{c}v_{K} is the characteristic planetesimalgas 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 Stype systems and will frequently appear in this work. Using these two parameters, strong coupling corresponds to small bodies, d_{p} ≪ d_{c}, when $\leftA{\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_{0}, surface density normalization Σ_{0} 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 Σ_{0} = 500 g cm^{−2} because of a secular resonance, which is not present in highermass 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 Σ_{0} and e_{0}); (iii) dynamically quiet locations do not appear for all values of e_{0} and Σ_{0} 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^{3} 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 v_{c} (a) and planetesimal size d_{c} (b) for different disk models, in which we vary disk eccentricity e_{0}, 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 Stype binaries is that not only the planetesimal eccentricity but also the apsidal angle takes on a unique, sizedependent 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 sizedependent 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 3DGaussian (or Schwarzschild) distribution (Stewart & Ida 2000; Rafikov 2003c,a). But in Stype 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 twobody 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 σ_{i}v_{K} ≫ v_{esc} (where ${v}_{\text{esc}}^{2}=2G({m}_{1}+{m}_{2})/({d}_{1}+{d}_{2})$ for colliding objects with sizes d_{1}, d_{2} and masses m_{1}, m_{2}), that is, gravitational focusing is negligible; the collision rate scales as ${\sigma}_{i}^{2}$ in the opposite limit σ_{i}v_{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^{−4}, 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_{1} ≪ d_{c} and another has size d_{2} ≫ 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_{1} ≈ d_{2} 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 singlestar 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_{1} − d_{2} 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 highvelocity collisions (although the regions of catastrophic disruption are not exactly centered on d_{c} because the strength of planetesimals is sizedependent, 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 coagulationfragmentation 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 d_{1} and d_{2} 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 v_{c} and d_{c} 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 d_{init} = 4.3 km used in the simulation described in Sect. 5.1. The point (d_{c}, d_{c}) 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 \text{\hspace{0.17em}}~\text{\hspace{0.17em}}{(h/{a}_{\text{p}})}^{2}\text{\hspace{0.17em}}~\text{\hspace{0.17em}}0.003$ (h is the disk scale height) is a measure of the particlegas 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 subKeplerian 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 e_{0} = 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 multiannulus (or multizone) coagulationfragmentation 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 planetesimalplanetesimal 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 selfsimilar (which is a standard and physically motivated assumption anyway). The coagulationfragmentation 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 dispersiondominated 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 lowvelocity 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 sizedependent 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 Stype 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^{−2} used to decide whether the largest fragment is distinct from the continuous spectrum of fragments (Sect. A.2); tolerance parameters ϵ_{1} = 0.05 and ϵ_{2} = 10^{−6} 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} = 0.01, Σ_{0} = 1000 g cm^{−2} at a_{0} = 1 AU, ϖ_{d} = 0°, σ_{i} = 10^{−4}, 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 onezone simulations of planetesimal growth at different semimajor axes in Sect. 5.1. We then move on to the global, multizone 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 dragdriven radial inspiral (naturally absent in the former case).
5.1 Coagulation in a single annulus
We ran our singleannulus 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^{−3} this corresponds to planetesimal mass of m = 10^{18} 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^{6} yr.
In all plots, the curves corresponding to times ≲10^{4} 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^{3}–10^{4} 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^{4} 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 highvelocity 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 powerlaw 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^{−1}, growth to hundreds of km occurs within several ×10^{5} 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^{−1}, 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^{3} yr), Fig. 5 shows that large bodies get eroded away by 10^{5} 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^{4} yr when a sufficient amount of small debris capable of efficiently eroding big bodies in highv_{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 coagulationfragmentation 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 singleannulus 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 v_{c} and critical size d_{c}, 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 d_{c}, and the vertical dashed lines correspond to the initial size of the planetesimals d_{init}. 
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, multiannulus simulation of planetesimal growth, fully accounting for the effects of sizedependent 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 pileup at small sizes). This is to be expected since the global size spectrum in Fig. 6 is an average of multiple singlezone 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^{4} yr vs. 3 × 10^{5} 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 coagulationfragmentation 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_{1} = 2 and d_{2} = 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 Σ_{0} = 10^{3} g cm^{−2} 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^{−1}, 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_{1} = 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^{4} 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^{−1} 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^{2}. 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 multizone model (Simulation 1). Different panels are produced for different initial planetesimal sizes d_{init} in terms of d_{min} (as labeled on the plot), which is d_{min} = 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 d_{init} 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 d_{init}. 
Fig. 9 Collision velocity between d_{1} = 2 and d_{2} = 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 d_{min} – 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 e_{0}, (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 d_{min}, 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 Σ_{0}, disk eccentricity normalization e_{0}, 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 Σ_{0} within a broad range relative to the fiducial value of 10^{3} g cm^{−2} 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 Σ_{0} (see Fig. 1), this does not dramatically affect d_{min}.
Figure 10a does show a weak decrease in d_{min} as Σ_{0} increases from 500 to 5000 g cm^{−2}. 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} ~ σ_{i}v_{K}). As Σ_{0} 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 Σ_{0}). As a result, within the range of Σ_{0} where the dynamically quiet location exists, collisions are less destructive and planetesimal growth becomes possible starting at lower d_{init} for larger Σ_{0}.
We note that in the Σ_{0} = 500 g cm^{2} 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^{3} 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, Σ_{0} > 5000 g cm^{−2}, 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^{4} and we find a slight increase in d_{min} with Σ_{0}.
In Simulation 7, we also considered a model of a massive axisymmetric disk, with Σ_{0} = 20, 000 g cm^{−2} (corresponding to a total disk mass of 0.07 M_{⊙} interior to 5 AU) and e_{0} = 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_{0}: 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_{0} changes from 0 to 0.1. This variation is nonmonotonic.
As e_{0} increases from zero to the fiducial value e_{0} = 0.01, we find d_{min} decreasing for reasons similar to the ones mentioned in Sect. 6.1: Higher e_{0} means larger B_{d} ∝ e_{0}Σ_{0} (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} ~ σ_{i}v_{K} ∝ a^{−1∕2} at this location, higher e_{0} results in lower v_{coll} (see Fig. 9c) promoting planetesimal growth and leading to lower d_{min}. In an axisymmetric disk (e_{0} = 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_{0} between 0.004 and 0.018.
For e_{0} > 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_{0} (Silsbee & Rafikov 2015b), leading to higher v_{coll} and a less growthfriendly environment as e_{0} increases (see Fig. 9c); hence the rapidly increasing values of d_{min}. Even a disk with e_{0} = 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_{0}, 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 Σ_{0} and e_{0} (Rafikov 2013b; Silsbee & Rafikov 2015b; Rafikov & Silsbee 2015a).
Parameters of our model runs.
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 antialigned 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^{−5} to almost 30 km at σ_{i} = 10^{−3}. 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 Σ_{0} and e_{0}). In this narrow zone v_{c} is very low and v_{coll} is set entirely by random velocity σ_{i}v_{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, rubblepile 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^{−4}). 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^{−4}, 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^{−3} that d_{min} is lower in the model with disk gravity than without. This reversal occurs because for σ_{i} = 10^{−3}, 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^{−5} to 6.9 km at σ_{i} = 10^{−5}. 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 rubblepile 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 onezone coagulationfragmentation 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 noinspiral 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 noinspiral simulation with fiducial parameters (χ_{ir}= 0, analogous to a simple addition of a number of multizone 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 fullinspiral 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 dragdriven 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 d_{init}∕d_{min} = 2, to be compared with Fig. 8c (see Sect. 6.7 for details). (b) Fiducial model computed for d_{init}∕d_{min} 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^{−1} regime (see Fig. 9), much higher than in disks around single stars. As a result, growth to 10^{2} 10^{3}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^{5} – 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 multidimensional 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_{0} 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 Stype 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 γ Ceplike system.
Our finding that d_{min} only weakly depends on Σ_{0} is quite important. As we vary Σ_{0} from 500 g cm^{−2} to 20,000 g cm^{−2} (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_{0} = 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 Σ_{0}, 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 Σ_{0} when a secular resonance appears in the disk; see Fig. 1).
Figure 12 shows a_{dq} as a function of e_{0} for various values of Σ_{0} in a system which has, other than the variation of e_{0} and Σ_{0}, our fiducial parameters. We have marked the fiducial model with an orange star along the Σ_{0} = 1000 g cm^{−2} 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_{0} 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 Σ_{0} = 500 g cm^{−2} 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 Σ_{0} = 500 g cm^{−2} (corresponding to 6M_{⊕} of solid material), the dynamically quiet region will not exist if e_{0} is greater than about 0.02, corresponding to a massweighted 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^{−4} featured in our fiducial model corresponds to random velocities of about 10 m s^{−1}. 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 twobody 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 dragdriven 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 nonmonotonic variation with d_{init}∕d_{min}, as well as the possibility of latetime planet migration (that must have certainly affected Stype 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^{5} yr as this ratioincreases from 1 to 5 in our fiducial model. The expectation of short disk lifetimes in Stype 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 Σ_{0} ≲ 5000 g cm^{−2} in this particular system.
Fig. 12 Location a_{dq} of the dynamically quiet zone in an apsidally aligned disk. We assume the fiducial disk plus binary parameters, except for Σ_{0} and e_{0} 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 Stype 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 multikilometer 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 metersize objects (Weidenschilling 1977). In Stype 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_{0} 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 Stype 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 (massweighted mean) disk eccentricities, ⟨e_{d} ⟩ ~ 0.1–0.3, with e_{d} profile typically exhibiting a doublepeaked 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 diskperturber 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^{6}, 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_{0}.
Similarly, accounting for the effect of disk selfgravity on its dynamics has also been shown to suppress ⟨e_{d} ⟩ (Marzari et al. 2009). Another interesting finding is that Stype 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 welldefined 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 antialignment (ϖ_{d} ≈ 180°), whereas Marzari et al. (2009) find either antialignment or precession, depending on whether selfgravity is included. Paardekooper et al. (2008) find either antialignment or precession depending on the flux limiter used in their code.
To summarize, the existing numerical studies of Stype disks in binaries have not yet arrived at a consensus regarding the values of either e_{0} or ϖ_{d}, and their dependence on various physical inputs, for example, realistic thermodynamics, disk selfgravity, viscosity, etc. Given the importance of these characteristics for understanding planet formation in Stype binaries, we encourage future numerical efforts to focus on measuring the true value of ϖ_{d} as well as disk eccentricity e_{0}, using the most realistic physical inputs, such as disk thermodynamics and selfgravity.
7.4 Comparison of detailed coagulationfragmentation calculations with simple heuristics
Full coagulationfragmentation 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^{−2}).
Now, armed with the full coagulationsimulation 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 singleannulus 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 coagulationfragmentation 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 simulationbased 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 coagulationfragmentation simulations.
Fig. 13 Minimum planetesimal size d_{min} resulting in the formation of large objects obtained with different methods for different values of v_{c} (indicated inpanels) and d_{c}. The solid blue line shows d_{min} computedusing our multiannulus simulations. The green line is the analytic estimate for d_{min} 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 d_{min} 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 Stype 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^{2} 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 coagulationfragmentation calculations fully accounting for the gas dragdriven 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 Σ_{0} 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 coagulationfragmentation framework for exploring planetesimal growth in Stype 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 massweighted mean disk eccentricity ≲0.05 (corresponding, in the context of our model, to e_{0} ≲ 0.02);

We identify disk eccentricity (e_{0}) 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 Σ_{0};

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 coagulationfragmentation 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 Ptype 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 15XRP1520139, 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 coagulationfragmentation 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 overcounting 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^{5} 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. Highvelocity collisions will result in many small fragments, and no large remnant bodies. Lowvelocity collisions will result in one body which contains nearly all of the precollision 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}}}=10.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 centerofmass 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}},\text{\hspace{1em}}{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 longterm 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 \\ \mathrm{max}({M}_{\text{tot}}b,0.5{M}_{\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}}}\text{\hspace{0.17em}}\times \text{\hspace{0.17em}}{({M}_{i}/{M}_{\text{cutoff}})}^{1+\xi}\text{\hspace{0.17em}}\times \text{\hspace{0.17em}}\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_{1}. 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 selfsimilar 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_{ij}v_{K} and v_{max} = e_{ij}v_{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}}\mathrm{exp}\left[\frac{{({k}_{\text{rand}})}^{2}}{2{\sigma}_{e}^{2}}\right],\\ \rho ({h}_{\text{rand}})=\frac{1}{{\sigma}_{e}\sqrt{2\pi}}\mathrm{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}}\mathrm{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 timestep would be too timeconsuming. 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_{1} points in vspace, 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 ${\int}_{0}^{{v}_{i}}F}(v)\text{d}v=(i+0.5)/{N}_{1$, for i ∈ [0, N_{1} − 1]. We then calculate the collision rates and outcomes for each set of collision partners at the beginning of the simulation, using the N_{1} different velocities. During the simulation, in each timestep, and for each set of collision partners, we randomly choose one of the N_{1} collisional outcomes and rates to use for that timestep.
To calculate the N_{1} values v_{i}, we initially create two vectors z_{rand} and z_{λ} of length N_{2} 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 ${\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_{3} samples of the collision velocity using Eq. (A.16) and randomly selected entries from z_{rand} and z_{λ}. We then select the N_{1} values of v_{i} so there are an equal number of the N_{3} sample velocities between consecutive values of v_{i}. We take N_{1} = 11, N_{2} = 100, and N_{3} = 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}}\mathrm{exp}\left(\frac{{q}^{2}}{2{\sigma}_{i}^{2}}\right);\text{\hspace{1em}}\rho (p)=\frac{1}{{\sigma}_{i}\sqrt{2\pi}}\mathrm{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}\mathrm{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 timestep, 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 twobody 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 sheardominated regime by replacing e_{*} in Eq. (A.23) with e_{*} + c_{1}e_{H}, and then multiplying the entire right hand side by 1 + c_{2}(e_{H}∕e_{*}), where c_{1} and c_{2} 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_{1} and c_{2} 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_{1} = 4.75, c_{2} = 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_{1} be the collision rate in the low secular velocity limit given by Eq. (A.24), and let R_{2} be the collision rate when the secular velocity dominates; R_{2} 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_{2}∕R_{1} = 2.8 in the limit that v_{esc}∕(σ_{e}v_{K}) ≪ 1, and 1.6∕λ^{2} in the limit that v_{esc}∕(σ_{e}v_{K}) ≫ 1. To make the collision rate continuous, we approximate it as $${R}_{ij}({v}_{ij})=f{R}_{1}+(1f){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_{ij}v_{K} is the relative secular velocity between bodies in massbins 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_{1} = (m_{1}− 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_{2} analogously for the second body. The relative secularly excited velocity v_{12} between particles 1 and 2 is then equalto v_{max}x_{1} − x_{2} and has a distribution (assuming x_{1} and x_{2} 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^{−4}.
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 Coagulationfragmentation time step
The time step for the coagulationfragmentation 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)}<\u03f5$$(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<{\u03f5}_{1}\text{\hspace{1em}}\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<{\u03f5}_{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 ϵ_{1}, or the decrease in mass is less than ϵ_{2} 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 ϵ_{1} and ϵ_{2} 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^{−6} have a negligible effect on d_{min}.
In the simulation used to create Fig. 4, we used a lower value of ϵ_{2} so that the distributions would remain smooth, even in regions of massspace 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}}\text{\hspace{0.17em}}\times \text{\hspace{0.17em}}\mathrm{max}\left(\frac{{v}_{\text{drift}}}{\Delta {r}_{i}}\right)={\u03f5}_{\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 = 10^{12} 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 A_{ij}, 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^{12} 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 highmass 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 wellequipped 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 timedependent behavior in fragmenting systems. However, there are some knownresults for the steadystate 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 powerlaw. 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 steadystate 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 powerlaw index of the size distribution of the steadystate 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 powerlaw 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 steadystate mass distribution with no upper or lower cutoff. 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 steadystate 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})\text{\hspace{0.17em}}~\text{\hspace{0.17em}}{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 N_{ann} = 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 N_{ann} 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 powerlaws, 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 powerlaw. Then the overabundance of small masses creates an underabundance of objects of slightly larger masses, giving rise to the wavelike pattern that we observe. This effect is reduced in panels c and d because in this example we tuned the Q_{0} 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 sizedependent 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_{0} 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_{0}, this leads to an evolution of the number density given by $$\mathfrak{n}(a,t)=\{\begin{array}{ll}{\mathfrak{n}}_{0}{s}^{2}\left(4a/({a}_{0}s)\right)\hfill & a<4{a}_{0}/(1+4{v}_{0}t/{a}_{0}),\hfill \\ 0\hfill & a>4{a}_{0}/(1+4{v}_{0}t/{a}_{0}),\hfill \end{array}$$(B.10)
where $s=1at{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 timestep, 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 inspiraldriven 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.
Variation of numerical parameters.
References
 Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progr. Theor. Phys., 56, 1756 [Google Scholar]
 Barenfeld, S. A., Carpenter, J. M., Sargent, A. I., et al. 2019, ApJ, 878, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Beaugé, C., Leiva, A. M., Haghighipour, N., & Otto, J. C. 2010, MNRAS, 408, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Belyaev, M. A., & Rafikov, R. R. 2011, Icarus, 214, 179 [CrossRef] [Google Scholar]
 Benedict, G. F., Harrison, T. E., Endl, M., & Torres, G. 2018, Res. Notes Am. Astron. Soc., 2, 7 [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [Google Scholar]
 Bromley, B. C., & Kenyon, S. J. 2011, ApJ, 731, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Chauvin, G., Beust, H., Lagrange, A. M., & Eggenberger, A. 2011, A&A, 528, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Davydenkova, I., & Rafikov, R. R. 2018, ApJ, 864, 74 [CrossRef] [Google Scholar]
 Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [NASA ADS] [CrossRef] [Google Scholar]
 Dvorak, R. 1982, Oesterreich. Akad. Wissensch. Math. Naturwissensch. Klasse Sitzungsber. Abteilung, 191, 423 [Google Scholar]
 Garaud, P., Meru, F., Galvagni, M., & Olczak, C. 2013, ApJ, 764, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Greenzweig, Y., & Lissauer, J. J. 1992, Icarus, 100, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Gyergyovits, M., Eggl, S., PilatLohinger, E., & Theis, C. 2014, A&A, 566, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harris, R. J., Andrews, S. M., Wilner, D. J., & Kraus, A. L. 2012, ApJ, 751, 115 [Google Scholar]
 Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383 [Google Scholar]
 Hayashi, C. 1981, Progr. Theor. Phys. Suppl., 70, 35 [Google Scholar]
 Heppenheimer, T. A. 1978, A&A, 65, 421 [NASA ADS] [Google Scholar]
 Heppenheimer, T. A. 1980, Icarus, 41, 76 [CrossRef] [Google Scholar]
 Ida, S., & Makino, J. 1993, Icarus, 106, 210 [Google Scholar]
 Inaba, S., Tanaka, H., Nakazawa, K., Wetherill, G. W., & Kokubo, E. 2001, Icarus, 149, 235 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 547 [Google Scholar]
 Kenyon, S. J.,& Bromley, B. C. 2002, AJ, 123, 1757 [NASA ADS] [CrossRef] [Google Scholar]
 Kenyon, S. J.,& Luu, J. X. 1998, AJ, 115, 2136 [NASA ADS] [CrossRef] [Google Scholar]
 Kenyon, S. J.,& Luu, J. X. 1999, AJ, 118, 1101 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., & Nelson, R. P. 2008, A&A, 486, 617 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kobayashi, H., Tanaka, H., Krivov, A. V., & Inaba, S. 2010, Icarus, 209, 836 [NASA ADS] [CrossRef] [Google Scholar]
 Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [Google Scholar]
 Martin, R. G., Lissauer, J. J., & Quarles, B. 2020, MNRAS, 496, 2436 [CrossRef] [Google Scholar]
 Marzari, F., & Scholl, H. 2000, ApJ, 543, 328 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., & Thebault, P. 2019, Galaxies, 7, 84 [CrossRef] [Google Scholar]
 Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marzari, F., Baruteau, C., Scholl, H., & Thebault, P. 2012, A&A, 539, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meschiari, S. 2014, ApJ, 790, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Miranda, R., & Rafikov, R. R. 2019, ApJ, 878, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, T. W. A., & Kley, W. 2012, A&A, 539, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Sol. Syst. Dyn., 253 [Google Scholar]
 Nesvold, E. R., & Kuchner, M. J. 2015, ApJ, 798, 83 [NASA ADS] [CrossRef] [Google Scholar]
 O’Brien, D. P., & Greenberg, R. 2003, icarus, 164, 334 [NASA ADS] [CrossRef] [Google Scholar]
 Ogilvie, G. I. 2001, MNRAS, 325, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Paardekooper, S.J., Thébault, P., & Mellema, G. 2008, MNRAS, 386, 973 [NASA ADS] [CrossRef] [Google Scholar]
 Picogna, G., & Marzari, F. 2013, A&A, 556, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rafikov, R. R. 2003a, AJ, 126, 2529 [CrossRef] [Google Scholar]
 Rafikov, R. R. 2003b, AJ, 125, 922 [CrossRef] [Google Scholar]
 Rafikov, R. R. 2003c, AJ, 125, 906 [CrossRef] [Google Scholar]
 Rafikov, R. R. 2003d, AJ, 125, 942 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R. 2004, AJ, 128, 1348 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R. 2013a, ApJ, 764, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R. 2013b, ApJ, 765, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R., & Silsbee, K. 2015a, ApJ, 798, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R., & Silsbee, K. 2015b, ApJ, 798, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Rafikov, R. R., Silsbee, K., & Booth, R. A. 2020, ApJS, 247, 65 [CrossRef] [Google Scholar]
 Regály, Z., Sándor, Z., Dullemond, C. P., & Kiss, L. L. 2011, A&A, 528, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Safronov, V. S. 1972, Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets. [Google Scholar]
 Schäfer, U., Yang, C.C., & Johansen, A. 2017, A&A, 597, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sefilian, A. A., Rafikov, R. R., & Wyatt, M. C. 2021, ApJ, 910, 13 [CrossRef] [Google Scholar]
 Silsbee, K., & Rafikov, R. R. 2015a, ApJ, 808, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Silsbee, K., & Rafikov, R. R. 2015b, ApJ, 798, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Smoluchowski,M. V. 1916, Zeitsch. Phys., 17, 557 [Google Scholar]
 Spaute, D., Weidenschilling, S. J., Davis, D. R., & Marzari, F. 1991, Icarus, 92, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Statler, T. S. 2001, AJ, 122, 2257 [NASA ADS] [CrossRef] [Google Scholar]
 Stewart, G. R., & Ida, S. 2000, Icarus, 143, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Stewart, S. T., & Leinhardt, Z. M. 2009, ApJ, 691, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Stewart, G. R., & Wetherill, G. W. 1988, Icarus, 74, 542 [NASA ADS] [CrossRef] [Google Scholar]
 Tanaka, H., Inaba, S., & Nakazawa, K. 1996, icarus, 123, 450 [NASA ADS] [CrossRef] [Google Scholar]
 Thebault, P. 2011, Celest. Mech. Dyn. Astron., 111, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Thebault, P., Kral, Q., & Olofsson, J. 2021, A&A, 646, A173 [CrossRef] [EDP Sciences] [Google Scholar]
 Thébault, P., Marzari, F., & Scholl, H. 2006, Icarus, 183, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Thébault, P., Marzari, F., & Scholl, H. 2008, MNRAS, 388, 1528 [NASA ADS] [CrossRef] [Google Scholar]
 Thébault, P., Marzari, F., & Scholl, H. 2009, MNRAS, 393, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Trubnikov, B. A. 1971, Sov. Phys. Dokl., 16, 124 [Google Scholar]
 Ward, W. R. 1981, Icarus, 47, 234 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
 Wetherill, G. W., & Stewart, G. R. 1989, Icarus, 77, 330 [Google Scholar]
 Wetherill, G. W., & Stewart, G. R. 1993, Icarus, 106, 190 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012, A&A, 544, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zsom, A., Sándor, Z., & Dullemond, C. P. 2011, A&A, 527, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zurlo, A., Cieza, L. A., Pérez, S., et al. 2020, MNRAS, 496, 5089 [Google Scholar]
The saturation of the size of the largest body at the level of several tens of meters in the outer disk is an artifact of our calculations (see Sect. 5.1).
Unlike v_{c} ∝ A^{−1}, neither e_{r} nor collisional velocity diverge at the resonance (where A → 0), since the product v_{c}A in Eq. (4) remains finite as A → 0.
We note that Fig. 9 does not reflect these changes well: The v_{coll} curves in its panel B look very similar to each other for Σ_{0} ≥ 2000 g cm^{−2} (as if they were fully dominated by random motions), whereas the d_{min} outcomes are not exactly the same. This is because this figure is made only for a particular pair of planetesimal sizes, not revealing the full picture of the behavior of v_{coll} across all planetesimal sizes.
Another work accounting for radiative effects by Picogna & Marzari (2013) is very different from the rest – it is a 3D Smoothed Particle Hydrodynamics (SPH) study – and found high disk eccentricity.
All Tables
All Figures
Fig. 1 Radial profiles of characteristic velocity v_{c} (a) and planetesimal size d_{c} (b) for different disk models, in which we vary disk eccentricity e_{0}, 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 d_{1} and d_{2} 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 v_{c} and d_{c} 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 d_{init} = 4.3 km used in the simulation described in Sect. 5.1. The point (d_{c}, d_{c}) 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 e_{0} = 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 singleannulus 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 v_{c} and critical size d_{c}, 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 d_{c}, and the vertical dashed lines correspond to the initial size of the planetesimals d_{init}. 

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 multizone model (Simulation 1). Different panels are produced for different initial planetesimal sizes d_{init} in terms of d_{min} (as labeled on the plot), which is d_{min} = 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 d_{init} 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 d_{init}. 

In the text 
Fig. 9 Collision velocity between d_{1} = 2 and d_{2} = 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 d_{min} – 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 e_{0}, (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 d_{min}, 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 d_{init}∕d_{min} = 2, to be compared with Fig. 8c (see Sect. 6.7 for details). (b) Fiducial model computed for d_{init}∕d_{min} just slightly exceeding unity, to be compared with Fig. 8c,d (see Sect. 7.1 for details). 

In the text 
Fig. 12 Location a_{dq} of the dynamically quiet zone in an apsidally aligned disk. We assume the fiducial disk plus binary parameters, except for Σ_{0} and e_{0} 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 d_{min} resulting in the formation of large objects obtained with different methods for different values of v_{c} (indicated inpanels) and d_{c}. The solid blue line shows d_{min} computedusing our multiannulus simulations. The green line is the analytic estimate for d_{min} 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 d_{min} 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 = 10^{12} 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 A_{ij}, 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 powerlaw 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 N_{ann} = 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 N_{ann} values, but ϵ_{drift} = 0.3. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.