Press Release
Open Access
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/0004-6361/202141139
Published online 17 August 2021

© K. Silsbee and R. R. Rafikov 2021

Licence Creative CommonsOpen 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 planet-hosting systems as important test beds of planet formation theory.

In particular, tight eccentric binary systems that host planets orbiting one of the stellar companions (so-called S-type systems in the classification of Dvorak 1982), such as γ Cephei (Hatzes et al. 2003), present important challenges to theories of planet formation via core accretion, which involves a steady agglomeration of planetary cores through mutual collisions of numerous small planetesimals. Indeed, gravitational perturbations from the eccentric companion in such binaries with semimajor axes ≲20–30 AU are expected to drive planetesimal eccentricities to high values. This naturally results in the destruction, rather than the growth, of planetesimals in mutual collisions. A simple calculation (Heppenheimer 1978), including only the dominant secular effects of the stellar companion on planetesimal orbits, yields planetesimal collision velocities of a few km s−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 size-dependent alignment of planetesimal orbits (Thébault et al. 2008; Rafikov & Silsbee 2015a) gives rise to a dynamical environment with a mixture of low- and high-velocity collisions, making it challenging, even qualitatively, to determine the evolution of the size distribution.

The problem is further complicated by the radial inspiral of planetesimals due to interaction with the gas. The gas disk is slightly pressure-supported and orbits ~ 30 m s−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 particle-gas eccentricity where radial migration slows down, thus accelerating coagulation. Inspiral may also help to alleviate the detrimental effect of the erosion of growing planetesimals by small bodies by rapidly removing such bodies from the system.

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

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

2 Model setup

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

A planet with mass Mpl = 1.6 MJ orbits the primary star of the γ Cephei system with the semimajor axis apl = 2 AU and eccentricity epl = 0.12, which we use to motivate our subsequent calculations. There are a number of other binaries with ab ≲ 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: (1)

where a0 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 aout = 5 AU (i.e., Σ(a) = 0 for a > aout). 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 eg (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 eg(a) behavior is beyond the scope of this study. Throughout this paper, we assume a solid-to-gas ratio of 0.01, and the disk to be slightly flared with the hr profile given by Eq. (14) from Rafikov & Silsbee (2015b) and varying between 0.03 and 0.05.

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

3 Physical ingredients of the model

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

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

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

3.1 Gravitational perturbations in S-type binaries

The motion of planetesimals around the primary in S-type binaries is affected by the gravity of both the secondary and the protoplanetary disk. Since the seminal work of Heppenheimer (1978), the direct effect of the secondary’s gravity has been accounted for using the secular approximation, which averages the (time-dependent) potential of the secondary over its orbit. However, since the gaseous protoplanetary disks in S-type binaries are themselves prone to developing eccentricity under the perturbation due to the secondary (Paardekooper et al. 2008; Regály et al. 2011), one needs to also account for the gravitational potential of a non-axisymmetric, eccentric disk. Silsbee & Rafikov (2015b) computed the potential due to an apsidally aligned disk with the surface density and eccentricity profiles given by Eq. (1); their calculation was subsequently generalized for a disk with arbitrary profiles of Σ(a) and eg (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) (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), Bb and Bd describe the torques due to the non-axisymmetric components of the potentials of the secondary and the disk, and A = Ab + Ad is the combined free precession rate due to the secondary (Ab) and the disk (Ad). The explicit expressions for Ad, Bd, Ab and Bb can be found in Silsbee & Rafikov (2015b, see their Eqs. (5)–(6) and (9)–(10), respectively). We note that while Ad depends only on the Σ(a) profile, Bd is determined by both Σ(a) and eg(a).

Of particular importance for planet formation in S-type binaries is the fact that the secondary companion always drives prograde apsidal precession of the planetesimal orbit (Ab > 0), while the precession due to disk gravity is usually retrograde (Ad < 0). Since Ab and Ad vary with a in different ways, one often finds that A = Ab + Ad → 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 short-period perturbations which act on the orbital period of the planetesimals (Murray & Dermott 1999). The magnitude of the short-period perturbations is quite small for bodies which are close together. If we consider two bodies with relative eccentricity, er, which will collide in the next orbit, then their distance is of order erap. Their relative acceleration due to the tidal perturbation of the secondary is . If this acceleration acts over a time it results in a relative velocity , which is suppressed by a factor relative to their typical relative velocity ervK in the absence of such short-period perturbations. Over longer timescales such short-period perturbations average out to zero.

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

The torque due to the secondary is expected to sharply (although not as discontinuously as we assume in Eq. (1)) truncate disk density at the semimajor axis aout. We take this truncation into account when computing disk gravity. The disk gravity-related coefficients Ad and Bd 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 .

3.2 Overview of planetesimal dynamics in S-type binaries

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

Rafikov & Silsbee (2015a) found that planetesimal dynamics admits two distinct regimes, depending on the degree of coupling between the planetesimal and the gas. Strongly coupled (small) planetesimals follow orbits nearly aligned with the local eccentric gas streamline. Weakly coupled (large) planetesimals have orbits nearly aligned with the forced eccentricity determined by the gravitational perturbations due to both the companion and the disk. Separating these two regimes is the characteristic planetesimal size dc that is defined as the planetesimal size at which the damping timescale1 τd of the planetesimal eccentricity (relative to the eccentricity vector of the local gas streamline eg = eg(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 dp (3)

(see their Eq. (32)), so that |d|≈ 1.27 for dp = dc.

Rafikov & Silsbee (2015a) also show that the relative planetesimal-gas eccentricity er = |eeg| is given by (see their Eq. (28)) (4)

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

Both dc and vc (or ec) are the key metrics of planetesimal dynamics in S-type systems and will frequently appear in this work. Using these two parameters, strong coupling corresponds to small bodies, dpdc, when and . Weak coupling is realized in the opposite limit of dpdc, when |d|≈ dpdc ≫ 1 and erec. This suggests another intuitive definition of vc (or ec) as the characteristic relative velocity (or eccentricity) between the weakly and strongly coupled (or large and small) objects.

The parameters vc and dc 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 e0, surface density normalization Σ0 and disk orientation ϖd. One can see that the behavior of vc and dc 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) dc and vc diverge at 2.5 AU in a model with Σ0 = 500 g cm−2 because of a secular resonance, which is not present in higher-mass models as it gets pushed out of the disk (i.e., disk gravity dominates planetesimal dynamics all the way out to aout); (ii) many models, in which the disk is apsidally aligned with the binary, feature a dynamically quiet location where dc, vc → 0 (with the semimajor axis varying with Σ0 and e0); (iii) dynamically quiet locations do not appear for all values of e0 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, dc, ec → 0 at the semimajor axis adq where (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, ≲ 103 yr for dp = 1 km objects, while for bigger bodies (for which τddp is longer) the timescale for collisions with comparable objects (capable of significantly perturbing eccentricity) is long.

thumbnail Fig. 1

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

3.3 Collision velocities

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

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

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

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

Varying σi has two effects. Let eij be the relative eccentricity between two bodies in the absence of any random motions. For collisions between bodies where eijσ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 eijσi, the collision velocity is proportionalto σi and the collision rate becomes independent of σi provided that σivKvesc (where for colliding objects with sizes d1, d2 and masses m1, m2), that is, gravitational focusing is negligible; the collision rate scales as in the opposite limit σivKvesc 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 vcoll between planetesimals as a function of their sizes. These are given by Eq. (A.16), with λ and erand assigned their mean values of 0.81 and , 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 eij). The three panels differ in the assumed location ap of the colliding planetesimals in the disk (which determines their vc and dc), as labeled on the figure. As mentioned in Sect. 3.2, vcoll ~ vc for pairs of collision partners in which oneobject has size d1dc and another has size d2dc. If the collision partners have similar size, or both are much larger (or both much smaller) than dc, 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 d1d2 would disappear and vcoll would monotonically increase along this line (since smaller planetesimals have lower random velocities as a result of gas drag, Rafikov 2003a, 2004). These differences clearly demonstrate the importance of apsidal alignment naturally emerging in disks around binaries for setting the pattern of planetesimal collision velocities.

3.4 Collision outcomes

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

Following Rafikov & Silsbee (2015b), we divide the collision outcomes into three classes. We say that a collision leads to catastrophicdisruption if the largest remnant contains less than half the combined mass of the two incoming planetesimals, erosion if the largest remnant is smaller than the larger of the two incoming bodies, and growth if the largest remnant is larger than either of the incoming bodies. We can then use the recipe of Stewart & Leinhardt (2009) to delineate in Fig. 2 the domains in d1d2 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 vc is very high, catastrophicdisruption occurs only between collision partners with similar (but not equal) sizes. It is mainly relevant for objects with sizes around dc, as that is where even objects of comparable size experience high-velocity collisions (although the regions of catastrophic disruption are not exactly centered on dc because the strength of planetesimals is size-dependent, with a minimum around 0.1 km (Stewart & Leinhardt 2009). In contrast to catastrophic disruption, erosion is possible even when the sizes of the colliding bodies are very unequal.

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

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

thumbnail Fig. 2

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

3.5 Radial inspiral

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

Using the results of Adachi et al. (1976), Rafikov & Silsbee (2015b) found that the inspiral rate ȧp is given by (their Eq. (11)): (6)

where er and τd are given by Eqs. (4) and (3), respectively, (h is the disk scale height) is a measure of the particle-gas velocity differential due to the pressure support in the gas disk, and α ~ 1 is a constant that depends on the disk model as described in Rafikov & Silsbee (2015b). From Eqs. (3) and (4), one can show (Rafikov & Silsbee 2015b) that er τddp. As a result, one finds that if either (i) dpdc, when er becomes constant, or (ii) erη, 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 are joined by an intermediate region where the inspiral velocity depends less strongly on dp. In all cases, inspiral speed is a monotonically decreasing function of dp, 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 vc vanishes at 1.8 AU, so the only contribution to the inspiral comes from the sub-Keplerian rotation of the gas and not from any relative eccentricity between the planetesimals and gas, significantly reducing ȧp.

These considerations highlight important differences in the radial drift behavior between the binary and single stars. First, high planetesimal eccentricities driven by the secular effect of the disk and the secondary result in faster radialdrift in disks in binaries compared to single systems (see Eq. (6)). Second, because higher er increases the inspiral rate, planetesimals will be naturally concentrated in regions where er 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).

thumbnail Fig. 3

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

4 Numerical ingredients of the model

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

4.1 Basics of the code structure

Our code models the evolution of the planetesimal size distribution in Nann 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 Nbins 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 ni is the number of planetesimals in mass bin i. To follow the evolution of ni in each annulus the code performs two basic functions.

First, in each of the annuli, it calculates the evolution of the size distribution due to planetesimal-planetesimal collisions. This is done in a standard fashion, effectively by solving a discretized version of the Smoluchowski equation (Smoluchowski 1916) accounting for the possibility of fragmentation (Sects. A.1, A.2, and A.6.1). Inclusion of fragmentation generally makes this an 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 , as long as the size distribution of fragments formed in a collision of two objects is self-similar (which is a standard and physically motivated assumption anyway). The coagulation-fragmentation component of the code has been extensively tested against the known analytical solutions as described in Sects. B.1 and B.2.

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

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

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

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

4.2 Parameters of the numerical model

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

Our grid in mass space employs Nbins = 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 dp = 10 m. The largest mass bin in our simulations corresponds to an object with dp = 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 Mi+1Mi = 1.37. To represent our global simulation domain extending from 1 AU to 4 AU we use Nann = 60 annuli, spaced in semimajor axis as described in Sect. A.5. In choosing this radial range we are motivated by the planetary semimajor axes in several S-type systems: apl = 2 AU in γ Cephei, apl = 2.6 AU in HD 196885, and apl = 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 e0 = 0.01, Σ0 = 1000 g cm−2 at a0 = 1 AU, ϖd = 0°, σi = 10−4, and a solid to gas ratio of 0.01. The total disk mass is 3.7MJ. We present the results for other disk models in Sect. 6.

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

5.1 Coagulation in a single annulus

We ran our single-annulus calculations of planetesimal growth at three semimajor axes – ap = 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 dinit = 4.3 km; for the adopted bulk density ρp = 3 g cm−3 this corresponds to planetesimal mass of m = 1018 g.

Figure 4 illustrates how the size distribution in our model system changes as a function of time. What is actually shown is mdNd ln dp, 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 (dp = 300 km). In panels (b) and (c), the simulation was stopped after 106 yr.

In all plots, the curves corresponding to times ≲104 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 dinit, which is marked by the dashed vertical black line. The wiggles just to the right of dinit are related to the finite spacing of the mass bins. The gap just to the left of dinit exists because debris of that size is not created in collisions of objects with size dinit – 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 (~ 103–104 yr), as it simply reflects the size distribution of fragments formed in collisions of roughly equal mass planetesimals of size dinit. 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 dinit 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 104 yr mdNd ln dp starts featuring a dip around dp ~ 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 dc for that simulation, which varies by a factor of several across the different environments, but stays pretty close to dp ~ 0.1 km. Particles with sizes arounddc experience high-velocity collisions with almost all other bodies, and are therefore preferentially destroyed, which likely affects the depth of thedip in different panels.

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

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

Careful examination of our simulation outputs shows that the largest objects grow primarily in collisions with objects of similar size, certainly for dp > dinit. 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 vc suppresses growth. In panel (a), with vc = 77 m s−1, growth to hundreds of km occurs within several ×105 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 vc, 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 vc = 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 dp ~ dinit gets gradually eroded by smaller objects produced in previous collisions (this can be seen in the decay of amplitude of dNd ln dp around that size). If we were to run this simulation beyond 1 Myr, large objects with dp ~ dinit 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 dinit only by ~2. But after a population of small debris develops in the system (already by 103 yr), Fig. 5 shows that large bodies get eroded away by 105 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 dp = dinit, which have low relative velocities. The differences emerge only after ~104 yr when a sufficient amount of small debris capable of efficiently eroding big bodies in high-vc 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 dp ~ dinit. 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 (dinit, dinit) is outside the solid white contours for all dynamic environments, as dinit = 4.3 km considerably exceeds dc in all panels.

It may seem strange that in all three panels in Fig. 2 the point (dinit, dinit) 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 dp = dinit 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 (dinit, dinit) sits outside the dotted contour, some amount of debris will still be produced in collisions of two objects of initial size dinit; 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 (dp = dinit) 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 vcoll has a reasonably broad distribution around this mean (see Sect. A.3). As a result, some collisions occur at higher velocities and are more destructive – a fact that is fully accounted for in our coagulation-fragmentation simulations depicted in Fig. 4. This is another clear illustration of the importance of considering the distribution of planetesimal velocities when studying their growth. The combination of the two aforementioned factors leads to erosion eventually becoming the dominant player in the collisional evolution in panels b and c of Fig. 4.

thumbnail Fig. 4

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

thumbnail Fig. 5

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

5.2 Full global calculation with inspiral

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

5.2.1 Outcome diagnostic

In any reasonable environment, planetesimal growth would proceed to very large sizes, eventually leading to planet formation, provided that the initial planetesimal size dinit is large enough. Indeed, objects withsizes dpdc 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 dinit, 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 dinit. Instead, one should be trying to figure out from what initial planetesimal size dinit can sustained planetesimal growth occur. In a given dynamical environment, there will be a minimum size dmin such that if the initial bodies are smaller than dmin then planetesimal growth in the system will eventually be halted, as in Fig. 4b,c, whereas if seed objects have dinit > dmin, 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, dmin < 4.3 km at ap = 2 AU, whereas dmin > 4.3 km at ap = 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 dmin 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 dmin 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, dinit 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 dinit across the whole disk and determine a single value of dmin for a given disk model based on this assumption.

To determine dmin defined in this manner, we employed the following iterative algorithm. We chose a value of dmin and used that as our starting size dinit. 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 dmin < dinit, whereas if it does not, then dmin > dinit. By trying several values of dinit and running our global simulations for each of them, we were able to first bracket dmin 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 dinit).

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 .

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 dinit expressed in terms of dmin = 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 dinit are barely there, and there is no pile-up at small sizes). This is to be expected since the global size spectrum in Fig. 6 is an average of multiple single-zone size distributions (such as the ones shown in Fig. 4) with the added complication of radial inspiral of solids. Nevertheless, some key features of the size distribution evolution (e.g., a dip in mdNd ln dp near dp ~ 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 dmin shows that there is a sharp transition over less than a factor of two in dinit 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 dinitdmin = 5 than for the one with dinitdmin = 2 (8 × 104 yr vs. 3 × 105 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 dinitdmin = 2 (5), 79% (88%) of the mass is lost to the creation of rubble smaller than the smallest body tracked in the simulation, with the remaining amount lost to inspiral into the central star.

Figures 7 and 8 provide details on how the coagulation-fragmentation process proceeds as a function of semimajor axis in the disk, for each global simulation shown in Fig. 6. Figure 7 shows dMda – 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 dinit, 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 vcoll between d1 = 2 and d2 = 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 erand taking their mean values. The black solid line corresponds to the fiducial disk model, whereas the black dotted line in panel (a) shows vcoll computed in the absence of random motions. One can easily recognize the “valley” forming around 1.8 AU as the dynamically quiet location where vc, dc → 0 in Fig. 1. Nonzero random velocities of planetesimals “fill" this valley, not allowing vcoll to drop to zero in their presence. Nevertheless, vcoll is still relatively low around this region.

A large increase in vcoll in the outer disk starting around 2 AU arises because for Σ0 = 103 g cm−2 the free precession rate A becomes very small there, while the excitation by the companion grows large. As a result, both vc and vcoll become very large, of order km s−1, even though a true secular resonance does not appear in this disk model. We note that vcoll is often several times smaller than vc, because dc is often less than d1 = 2 km.

Figure 7 shows that, in the cases with dinitdmin < 1, most of the solid mass is removed from the outer and inner portions of the disk on timescales of 104 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) vcoll is relativelylow. These factors naturally cause planetesimal accumulation (dMda 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 vcoll ~ 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 simulations2. 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, vcoll 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 (dinitdmin > 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 dinitdmin = 5 we find that in the outer disk, the largest bodies have size (2−3)dinit 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.

thumbnail Fig. 6

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

thumbnail Fig. 7

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

thumbnail Fig. 8

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

thumbnail Fig. 9

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

thumbnail Fig. 10

Dependence of our key metric dmin – a minimum size of seed planetesimals that results in their eventual growth beyond 300 km – on some parameters of our simulations: (a) surface density normalization Σ0, (b) disk eccentricity normalization e0, (c) apsidal angle of the disk ϖd, (d) amplitude of the random velocity of planetesimals σi, and (e) parameter χir (artificially) regulating the magnitude of planetesimal inspiral. Solid lines and diamond points correspond to dmin, whereas dashed lines and circles are for – the minimum size for collisionally weak planetesimals. The causes of the trends shown here are discussed in Sects. 6.16.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 e0, 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 dmin (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 dmin and . 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 103 g cm−2 has almost no effect on dmin, which ranges between 1.1 and 1.7 km (see Fig. 10a and Table 1). Even though both the dynamically quiet location at which vc and dc vanish and the location (and existence) of the secular resonance change with Σ0 (see Fig. 1), this does not dramatically affect dmin.

Figure 10a does show a weak decrease in dmin 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 vc and dc vanish, leaving only random velocity to produce relative velocity between the colliding planetesimals (i.e., vcoll ~ σivK). 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 dinit for larger Σ0.

We note that in the Σ0 = 500 g cm2 model a secular resonance (Silsbee & Rafikov 2015b) appears in the disk around 2.4 AU (see Fig. 1), driving vcoll to very high values3 at moderate (2–2.5 AU) semimajor axes. However, dmin 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 location4 and we find a slight increase in dmin 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 e0 = 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 dmin = 0.7 km (see Table 1).

6.2 Variation of e0: Simulations 8–12

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

As e0 increases from zero to the fiducial value e0 = 0.01, we find dmin decreasing for reasons similar to the ones mentioned in Sect. 6.1: Higher e0 means larger |Bd |∝ e0Σ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 vcoll ~ σivKa−1∕2 at this location, higher e0 results in lower vcoll (see Fig. 9c) promoting planetesimal growth and leading to lower dmin. In an axisymmetric disk (e0 = 0), Bd = 0 and cancellation of the secondary’s torque does not occur, leading to a rather high dmin = 4.2 km. The dynamically quiet location exists in our simulation domain for values of e0 between 0.004 and 0.018.

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

We note that while the location of the dynamically quiet region is sensitive to e0, 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 e0 (Rafikov 2013b; Silsbee & Rafikov 2015b; Rafikov & Silsbee 2015a).

Table 1

Parameters of our model runs.

6.3 Variation of ϖd: Simulations 13–15

We also find substantial sensitivity of dmin to ϖd. Figure 10c shows that increasing the deviation of the disk orientation from apsidal alignment with the binary orbit leads to larger dmin: Keeping everything else fixed, dmin grows from 1.6 km for an aligned disk (ϖd = 0°) to 9.7 km for an anti-aligned disk (ϖd = 180°). Even a ϖd = 10° misalignment is sufficient to increase dmin 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 vc can vanish only in apsidally aligned disks (see Fig. 1a). A small misalignment makes vc nonzero everywhere (Rafikov & Silsbee 2015a,b) but only slightly modifies the behavior of vcoll (see, e.g., the ϖd = 10° curve in Fig. 9d), because (i) vc 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 dmin 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 dmin on the magnitude of random motions σi: for collisionally strong planetesimals dmin 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 e0). In this narrow zone vc is very low and vcoll is set entirely by random velocity σivK (see the black solid and dotted curves in Fig. 9a). As a result, even at this favorable location in the disk, dmin is limited by a combination of radial inspiral and random motions that can destroy small planetesimals if σi gets large, requiring larger seed bodies to ensure steady growth. We anticipate that σi would be less important in models which do not have a dynamically quiet region (see Sect. 6.6).

6.5 Strong versus weak planetesimals

Dashed curves in Fig. 10 show our results for collisionally weak, rubble-pile planetesimals as defined in Stewart & Leinhardt (2009). Rather unsurprisingly, we find that for such planetesimals is larger than dmin 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 Bd = Ad = 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 vc (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 dmin = 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 vcoll in the outer disk is lower than with it, but this still does not help to keep dmin low.

In the absence of disk gravity, reducing σi has very little effect on dmin 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 dmin 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 dmin 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, dmin 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 rubble-pile planetesimals.

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

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

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

We believe that there are two reasons for this behavior. First, in the outer disk, inspiral does not help much with cleaning out small fragments, simply because they get primarily ground down below the smallest size captured in our simulations. Second, the radial inspiral negatively impacts the growth of large planetesimals by shifting their semimajor axis away from the dynamically quiet zone where their growth is most favorable. Indeed, Fig.11a shows the results of a no-inspiral simulation with fiducial parameters (χir= 0, analogous to a simple addition of a number of multi-zone simulations not interacting with each other) for dinitdmin = 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 adq = 1.8 AU where vc = 0, which certainly facilitates their growth compared to the full-inspiral calculation, in which the 300 km object emerges interior to 1.4 AU. But the general conclusion that we can draw from this exercise is that the gas drag-driven inspiral has a rather modest effect on the outcome of planetesimal growth.

thumbnail 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 dinitdmin = 2, to be compared with Fig. 8c (see Sect. 6.7 for details). (b) Fiducial model computed for dinitdmin 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 102 -103-km size objects in binaries is possible only if planetesimals start already rather large, with dinit ~ 1−10 km.

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

Our parameter space exploration in Sect. 6 is not entirely exhaustive: because of the multi-dimensional space of physical parameters we could vary at most one or two of them at once with respect to our fiducial model, meaning that some domains of the parameter space are left unexplored. Also, our models rely on certain simplifying assumptions that were necessary to reduce the number of model inputs, for example, the initial planetesimal size dinit 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 dmin in our fiducial case, by about a factor of four compared to the case in which disk gravity is ignored (see Fig. 10d). Eccentricity e0 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 dmin compared to the case in which inspiral is ignored (see Fig. 10e). But the overall effect is not as large as one might have expected (see Sect. 6.7). Next we discuss some other consequences of our results.

7.1 Implications for planet formation in binaries

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

Our finding that dmin 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 MJ to 70 MJ (for the assumed Σ profile (1) with a0 = 1 AU and aout = 5 AU), dmin 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 adq as a function of e0 for various values of Σ0 in a system which has, other than the variation of e0 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 e0 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 dp > 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 e0 is greater than about 0.02, corresponding to a mass-weighted mean disk eccentricity of 0.05. This is within the range found in numerical simulations, but on the low end of that range (see Sect. 7.3 for details).

The low levels of planetesimal random motions required to drive dmin 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 two-body relaxation would pump up the random motions to. Moreover, given the high relative velocities (up to ~ vc) 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 dinitdmin ≥ 1. If the size of seed planetesimals dinit is close to dmin, planetesimals grow rather slowly, causing their substantial gas drag-driven inspiral, so that dp reaches 300 km somewhat closer to the star than adq (see Fig. 11b). As dinitdmin increases, planetesimals tend to grow to large sizes closer to adq (see Fig. 8c for dinitdmin = 2). However, further increase in dinitdmin causes the location of preferential planetesimal growth to shift inward again (see Fig. 8d), presumably because of faster coagulation there. Such non-monotonic variation with dinitdmin, as well as the possibility of late-time planet migration (that must have certainly affected S-type systems harboring hot Jupiters), complicate the effort to constrain the initial planetesimal size dinit 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 dinitdmin, steadily falling from ~1 Myr to < 105 yr as this ratioincreases from 1 to 5 in our fiducial model. The expectation of short disk lifetimes in S-type binaries (Harris et al. 2012; Barenfeld et al. 2019; Zurlo et al. 2020) may then favor larger dinit, 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.

thumbnail Fig. 12

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

7.2 Implications for planetesimal origin

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

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

On the contrary, the competing mechanism relying on the operation of streaming instability appears to naturallyproduce large planetesimals with sizes in the range from tens to hundreds of kilometers (Simon et al. 2016; Schäfer et al. 2017). This mechanism is also rather fast (with a typical timescale of tens of local orbital periods), resulting in a rapid production of the population of seed planetesimals with dinit exceeding dmin, 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 e0 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 ab ≲ 30 AU. Instead, onehas to rely on numerical simulations for guidance in that matter.

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

For example, studies using a locally isothermal equation of state (Kley & Nelson 2008; Paardekooper et al. 2008; Marzari et al. 2009; Zsom et al. 2011; Regály et al. 2011; Martin et al. 2020) typically find rather high (mass-weighted mean) disk eccentricities, ⟨ed ⟩ ~ 0.1–0.3, with ed profile typically exhibiting a double-peaked structure as a function of semimajor axis and with high values of ed close to the disk center. These results may need to be viewed with caution in light of the recently identified issues with using the locally isothermal equation of state for numerical modeling of the disk-perturber tidal interaction (Miranda & Rafikov 2019). On the contrary, studies that attempt to represent disk thermodynamics in a more realistic way (Marzari et al. 2012; Gyergyovits et al. 2014) find a dramatic reduction of disk eccentricity6, down to ⟨ed⟩ ~ 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, ⟨ed ⟩ = 2.5e0.

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

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

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

7.4 Comparison of detailed coagulation-fragmentation calculations with simple heuristics

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

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

Now, armed with the full coagulation-simulation code, we can test the performance of each of these growth criteria in predicting the success of planetesimal evolution. To that effect we used our simulations in a single-annulus mode to calculate dmin for a variety of models featuring different values of dc and vc. As discussed in Rafikov & Silsbee (2015a), the collision rates and velocities are entirely determined by dc and vc (for the assumed level of random motions), regardless of what combination of disk parameters lead to the specific dc and vc (we set vK as appropriate at 2 AU in the disk around γ Cephei). Therefore, by exploring the range of plausible vc and dc 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 vc), with the solid blue line showing dmin as a function of dc, computed by our code using the procedure described in Sect. 5.2.1. The red line shows the minimum value of dinit such that no catastrophic disruption occurs for planetesimals of size dinit or larger. The green line shows the minimum value of dinit such that no body larger than dinit 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 dmin between two and five times smaller (depending on dc and vc) than the one determined by the coagulation-fragmentation simulation. On the contrary, the criterion of no erosion between bodies with a mass ratio greater than 1% is too restrictive, with the green curve exceeding the blue curve by a factor two to five.

Interestingly, the blue curve is often roughly equally separated from both red and green ones (in log space), which motivated us to also show the dashed blue line computed as the geometric mean of the predictions from the “no catastrophic disruption” and “no erosion” models. As one can see, it matches the simulation-based dmin (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 dmin based on the two criteria used in Rafikov & Silsbee (2015b) can be employed as a rough approximation for predicting the success of planetesimal growth in the absence of sophisticated coagulation-fragmentation simulations.

thumbnail Fig. 13

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

7.5 Comparison with previous work

A number of past studies of planet formation in S-type binaries have arrived at rather pessimistic conclusions regarding planetesimal growth starting from a population of 1–10 km objects (Thébault et al. 2006, 2008, 2009; Beaugé et al. 2010). In particular, Thebault (2011) found that even dinit ≳ 102 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 dinit ~ 1 km (see Fig. 10).

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

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

8 Conclusions

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

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

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

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

  • We identify disk eccentricity (e0) 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 coagulation-fragmentation simulations to provide a new heuristic criterion for successful planetesimal evolution in Sect. 7.4;

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

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

Acknowledgements

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

Appendix A Detailed description of the numerical methods

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

A.1 Mass redistribution due to collisions

In a given annulus, our code solves the standard coagulation-fragmentation problem, which effectively is a version of the coagulation equation of Smoluchowski (1916) modified to account for fragmentation (see Eq. (1) of Rafikov et al. 2020). Our framework employs mass space discretization in the form of logarithmically spaced mass bins with constant Mi+1Mi. 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 Rij between planetesimals in bin i and bin j for all i and j in [1, Nbins]. This calculation is described in Appendix A.4. The number of collisions Zij 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 RijΔt: (A.1)

We define an operator such that Fijk 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 (A.2)

where (A.3)

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

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

A.2 Collisional outcomes

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

Given a value of vcoll, we use the recipe provided in Stewart & Leinhardt (2009) to determine the mass of the largest remnant Mlr. We take the mass of the largest remnant to be (A.4)

where Mtot = mi + mj is the total mass of the two incoming bodies, QR is the specific center-of-mass energy of the collision given by (A.5)

and is the critical value of QR which leads to catastrophic disruption (which corresponds to Mlr < 0.5Mtot); is a function of the sizes and composition of the planetesimals, and vcoll, 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 Mlr is calculated, it will in general be between Mi and Mi+1 for some i. We then add a fraction Pi of the mass to bin i and Pi+1 to bin i + 1 given by (A.6)

We assume that the rest of the fragments follow a mass distribution with a power law index ξ. The precise value of ξ does not affect the long-term outcome of fragmentation cascade (O’Brien & Greenberg 2003). The upper cutoff mass for this distribution is given by (A.7)

Here b is an adjustable parameter which determines the cutoff mass for catastrophic collisions: whenever Eq. (A.4) yields Mlr< 2bMtot, 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 Mlr is small.

The fragment mass is assigned to mass bins with 0 < iicutoff, where icutoff is the index of the largest bin with mass less than Mcutoff. This is donesuch that (A.8)

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

In order to solve for ϒjk, we note that (A.9)

where nonpositive values of i correspond to (hypothetical) mass bins with Mi less than the minimum size considered in our simulations. Letting ℜ ≡ Mi+1Mi, we can then solve for ϒjk and write (A.10)

Some mass is lost from the simulation because it is contained in bodies with mass less than M1. This lost mass is equal to (A.11)

Because of the self-similar form of the fragment mass spectrum, we were able to implement this fragmentation algorithm with numerical cost scaling as 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 vcoll 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 eij between planetesimals as a function of their masses. Given eij, there is a distribution of relative encounter velocities vij between vmin = 0.5eijvK and vmax = eijvK, where vK is the local Keplerian speed. This distribution is given by Eq. (62) of RS15a as (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 (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) = (ks + krand, hs + hrand) where (ks, hs) and (krand, hrand) are the secular and random components, respectively; krand and hrand are assumed to be drawn from independent Gaussian distributions (A.14)

where we take σe = 2σi. This leads to a distribution of the relative random eccentricity given by theRayleigh distribution (A.15)

In principle,for each set of collision partners, we should then take the collision velocity to be given by (A.16)

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

In practice, doing this at each time-step would be too time-consuming. Instead, we use the following algorithm. Let F(v) be the true collision velocity distribution for a given set of collision partners in a given annulus. We approximately calculate N1 points in v-space, such that there is a constant fraction of the distribution F(v) between consecutive points: this means that the ith point vi satisfies the relation , for i ∈ [0, N1 − 1]. We then calculate the collision rates and outcomes for each set of collision partners at the beginning of the simulation, using the N1 different velocities. During the simulation, in each timestep, and for each set of collision partners, we randomly choose one of the N1 collisional outcomes and rates to use for that timestep.

To calculate the N1 values vi, we initially create two vectors zrand and zλ of length N2 with a constant fraction of the distribution for erand between each entry of zrand and a constant fraction of the distribution of λ between each entry of zλ. This means that , and a similar relation exists between zλ and the distribution for λ implied by Eq. (A.12).

We then draw N3 samples of the collision velocity using Eq. (A.16) and randomly selected entries from zrand and zλ. We then select the N1 values of vi so there are an equal number of the N3 sample velocities between consecutive values of vi. We take N1 = 11, N2 = 100, and N3 = 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 vcoll 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 emax is replaced with vmaxap). This gives a flux Fij of planetesimals from bin j past a planetesimalin bin i of (A.17)

where Iij 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: (A.18)

The distribution of is therefore given [analogous to Eq. (A.15)] by (A.19)

Then, averaged over Iij, the average value of Fij is given by (A.20)

The collision rate for a given particle in bin i with the population in bin j is given by Rij = ⟨Fijσij, where σij is the collision cross section including gravitational focusing: (A.21)

where is the geometric cross section, and . As discussed in Sect. A.3, at each time-step, we assume a particular value of the collision velocity drawn from the distribution F(v). For the purposes of calculating the collision rate, we assume λ to be in the same quantile of the distribution for λ as v is for F(v).

In the second dynamical regime of dominant random motions () we use the two-body accretion formula from Eq. (17) of Greenzweig & Lissauer (1992): (A.22)

Here, and are functions of I* = σiσe. We assume that I* = 1∕2, in which case , and . 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: (A.23)

We approximate the transition to the shear-dominated regime by replacing e* in Eq. (A.23) with e* + c1eH, and then multiplying the entire right hand side by 1 + c2(eHe*), where c1 and c2 are interpolation constants. This gives us a formula (A.24)

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

We expect a transition between the two dynamical regimes to occur when , where (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 R1 be the collision rate in the low secular velocity limit given by Eq. (A.24), and let R2 be the collision rate when the secular velocity dominates; R2 is given by (A.25)

Assuming σiσe = 1∕2, then at , ignoring shear, R2R1 = 2.8 in the limit that vesc∕(σevK) ≪ 1, and 1.6∕λ2 in the limit that vesc∕(σevK) ≫ 1. To make the collision rate continuous, we approximate it as (A.26)

with . vij = λeijvK is the relative secular velocity between bodies in mass-bins i and j for that time step.

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

Let vmax,i be the relative (secularly excited) velocity between the heaviest body in the ith mass bin and the lightest. Let x1 = (m1mbottom)∕(mtopmbottom) be the position of body 1 in the mass bin, where mtop and mbottom are the masses of the heaviest and lightest bodies in the bin. We define x2 analogously for the second body. The relative secularly excited velocity v12 between particles 1 and 2 is then equalto vmax|x1x2| and has a distribution (assuming x1 and x2 to be uniformly distributed in [0, 1]) (A.27)

Therefore, ⟨Rii⟩, the rate at which a particle in bin i experiences collisions with other particles in bin i, is given by (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 Δri be the radial width of annulus i, and let the drift rate of a particle in annulus i and mass bin j be vi,j. We then let the change in the number of particles in annulus i and mass bin j in a time interval Δtinsp be given by (A.29)

We use Nann = 60 annuli. These were spaced so as to minimize changes in ec or a within one bin. The exact bin boundaries depend on the run of ec(a) which varies between our different disk models. For the fiducial model, a changed by less than 16% between adjacent bins and ec by less than 41% in the cases where ec > 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 Coagulation-fragmentation time step

The time step for the coagulation-fragmentation process must be short enough that the distribution does not change appreciably during the length of a time step. In principle, this means that we would like (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 Δtcoag such that for all i (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 Δtcoag 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 dmin.

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 mass-space containing a negligible amount of mass. The success of the coagulation process, and the shape of the distribution at sizes greater than dinit 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 Δtinsp such that (A.32)

which guarantees that the inspiraling particles do not cross the full radial bin width in time Δtinsp (vdrift 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.

thumbnail Fig. B.1

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

B.1 Tests against known solutions of the coagulation equation

In this section we test the coagulation module of our code. We describe the results of simulations of the coagulation alone, without fragmentation, and with simplified collision rates and initial conditions that have analytic solutions we can compare with. There are three known analytic solutions to the Smoluchowski coagulation equation (B.1)

obtained for kernels Aij = 1 (Smoluchowski 1916), Aij = i + j and Aij = 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 Aij = 1. In this figure we have used a mass bin spacing Mi+1Mi = 1.15. Initially there were N = 1012 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 Aij = i + j. The distribution evolves more rapidly than in the case of the constant kernel because of the higher collision rates. In this case, in addition to the scatter at low masses, we see that the numerical solution is slightly behind the analytic solution at the high-mass end of the distribution.

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

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

B.2 Tests of the fragmentation module

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

O’Brien & Greenberg (2003) present a model of collisional fragmentation that leads to a steady state size distribution given by a power-law. Their model has one free parameter s, which parametrizes the strength of particles of size dp in catastrophic disruption events as (B.2)

They show that the distribution of debris sizes Df in the steady-state collisional cascade is given by (B.3)

where nf is the number of objects of size Df (per unit Df), and B is an overall scaling constant. The power-law index of the size distribution of the steady-state collisional cascade is given by (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 , one can show that y gets modified to (B.5)

thumbnail Fig. B.2

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

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

Figure B.2 shows a comparison of the steady-state result of the simulation against the expected power law distribution from Eq. (B.5) for three different values of s and α. The simulations were run with bin spacing Mi+1Mi = 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 nfd 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 , we can show that (B.6)

where B′ is the appropriate normalization. The number of particles contained in the bin with bin mass Mi scales as , 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.

thumbnail Fig. B.3

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

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

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

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

B.3 Tests of the radial inspiral module

In this section, we describe our test of the code module handling the size-dependent radial inspiral. We made this test for a simple inspiral rate given by (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 (B.8)

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

Adopting the initial planetesimal distribution in the disk in the form (B.9)

for 0 < a < 4a0, this leads to an evolution of the number density given by (B.10)

where .

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 Nann = 30, 60, and 240 radial annuli between 1 and 4 AU. The annuli are located at the same values of a as in Simulations 5, 1, and 7, respectively, and we take ϵdrift = 1.0. The corresponding bottom panels b, d, and f are for ϵdrift = 0.3.

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

Appendix C Sensitivity to numerical parameters

Our code employs a number of numerical parameters to set various time steps, properly resolve mass and inspiral-driven evolution, and so on. We varied these numerical code inputs to see the effect on the calculation outcomes. Table C.1 shows the results of this exercise, in which we use dmin and 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 dmin does not depend strongly on the numerical parameters, which justifies our particular choice of their values.

Table C.1

Variation of numerical parameters.

References

  1. Adachi, I., Hayashi, C., & Nakazawa, K. 1976, Progr. Theor. Phys., 56, 1756 [Google Scholar]
  2. Barenfeld, S. A., Carpenter, J. M., Sargent, A. I., et al. 2019, ApJ, 878, 45 [NASA ADS] [CrossRef] [Google Scholar]
  3. Beaugé, C., Leiva, A. M., Haghighipour, N., & Otto, J. C. 2010, MNRAS, 408, 503 [NASA ADS] [CrossRef] [Google Scholar]
  4. Belyaev, M. A., & Rafikov, R. R. 2011, Icarus, 214, 179 [CrossRef] [Google Scholar]
  5. Benedict, G. F., Harrison, T. E., Endl, M., & Torres, G. 2018, Res. Notes Am. Astron. Soc., 2, 7 [Google Scholar]
  6. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [Google Scholar]
  7. Bromley, B. C., & Kenyon, S. J. 2011, ApJ, 731, 101 [NASA ADS] [CrossRef] [Google Scholar]
  8. Chauvin, G., Beust, H., Lagrange, A. M., & Eggenberger, A. 2011, A&A, 528, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Davydenkova, I., & Rafikov, R. R. 2018, ApJ, 864, 74 [CrossRef] [Google Scholar]
  10. Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [Google Scholar]
  11. Dvorak, R. 1982, Oesterreich. Akad. Wissensch. Math. Naturwissensch. Klasse Sitzungsber. Abteilung, 191, 423 [Google Scholar]
  12. Garaud, P., Meru, F., Galvagni, M., & Olczak, C. 2013, ApJ, 764, 146 [NASA ADS] [CrossRef] [Google Scholar]
  13. Greenzweig, Y., & Lissauer, J. J. 1992, Icarus, 100, 440 [NASA ADS] [CrossRef] [Google Scholar]
  14. Gyergyovits, M., Eggl, S., Pilat-Lohinger, E., & Theis, C. 2014, A&A, 566, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Harris, R. J., Andrews, S. M., Wilner, D. J., & Kraus, A. L. 2012, ApJ, 751, 115 [Google Scholar]
  16. Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Hayashi, C. 1981, Progr. Theor. Phys. Suppl., 70, 35 [Google Scholar]
  18. Heppenheimer, T. A. 1978, A&A, 65, 421 [NASA ADS] [Google Scholar]
  19. Heppenheimer, T. A. 1980, Icarus, 41, 76 [CrossRef] [Google Scholar]
  20. Ida, S., & Makino, J. 1993, Icarus, 106, 210 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Inaba, S., Tanaka, H., Nakazawa, K., Wetherill, G. W., & Kokubo, E. 2001, Icarus, 149, 235 [NASA ADS] [CrossRef] [Google Scholar]
  22. 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]
  23. Kenyon, S. J.,& Bromley, B. C. 2002, AJ, 123, 1757 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kenyon, S. J.,& Luu, J. X. 1998, AJ, 115, 2136 [NASA ADS] [CrossRef] [Google Scholar]
  25. Kenyon, S. J.,& Luu, J. X. 1999, AJ, 118, 1101 [NASA ADS] [CrossRef] [Google Scholar]
  26. Kley, W., & Nelson, R. P. 2008, A&A, 486, 617 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Kobayashi, H., Tanaka, H., Krivov, A. V., & Inaba, S. 2010, Icarus, 209, 836 [NASA ADS] [CrossRef] [Google Scholar]
  28. Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [Google Scholar]
  29. Martin, R. G., Lissauer, J. J., & Quarles, B. 2020, MNRAS, 496, 2436 [CrossRef] [Google Scholar]
  30. Marzari, F., & Scholl, H. 2000, ApJ, 543, 328 [NASA ADS] [CrossRef] [Google Scholar]
  31. Marzari, F., & Thebault, P. 2019, Galaxies, 7, 84 [CrossRef] [Google Scholar]
  32. Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Marzari, F., Baruteau, C., Scholl, H., & Thebault, P. 2012, A&A, 539, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Meschiari, S. 2014, ApJ, 790, 41 [NASA ADS] [CrossRef] [Google Scholar]
  35. Miranda, R., & Rafikov, R. R. 2019, ApJ, 878, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Müller, T. W. A., & Kley, W. 2012, A&A, 539, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Murray, C. D., & Dermott, S. F. 1999, Sol. Syst. Dyn., 253 [Google Scholar]
  38. Nesvold, E. R., & Kuchner, M. J. 2015, ApJ, 798, 83 [NASA ADS] [CrossRef] [Google Scholar]
  39. O’Brien, D. P., & Greenberg, R. 2003, icarus, 164, 334 [NASA ADS] [CrossRef] [Google Scholar]
  40. Ogilvie, G. I. 2001, MNRAS, 325, 231 [NASA ADS] [CrossRef] [Google Scholar]
  41. Paardekooper, S.-J., Thébault, P., & Mellema, G. 2008, MNRAS, 386, 973 [NASA ADS] [CrossRef] [Google Scholar]
  42. Picogna, G., & Marzari, F. 2013, A&A, 556, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Rafikov, R. R. 2003a, AJ, 126, 2529 [CrossRef] [Google Scholar]
  44. Rafikov, R. R. 2003b, AJ, 125, 922 [CrossRef] [Google Scholar]
  45. Rafikov, R. R. 2003c, AJ, 125, 906 [CrossRef] [Google Scholar]
  46. Rafikov, R. R. 2003d, AJ, 125, 942 [NASA ADS] [CrossRef] [Google Scholar]
  47. Rafikov, R. R. 2004, AJ, 128, 1348 [NASA ADS] [CrossRef] [Google Scholar]
  48. Rafikov, R. R. 2013a, ApJ, 764, L16 [NASA ADS] [CrossRef] [Google Scholar]
  49. Rafikov, R. R. 2013b, ApJ, 765, L8 [NASA ADS] [CrossRef] [Google Scholar]
  50. Rafikov, R. R., & Silsbee, K. 2015a, ApJ, 798, 69 [NASA ADS] [CrossRef] [Google Scholar]
  51. Rafikov, R. R., & Silsbee, K. 2015b, ApJ, 798, 70 [NASA ADS] [CrossRef] [Google Scholar]
  52. Rafikov, R. R., Silsbee, K., & Booth, R. A. 2020, ApJS, 247, 65 [CrossRef] [Google Scholar]
  53. Regály, Z., Sándor, Z., Dullemond, C. P., & Kiss, L. L. 2011, A&A, 528, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Safronov, V. S. 1972, Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets. [Google Scholar]
  55. Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Sefilian, A. A., Rafikov, R. R., & Wyatt, M. C. 2021, ApJ, 910, 13 [CrossRef] [Google Scholar]
  57. Silsbee, K., & Rafikov, R. R. 2015a, ApJ, 808, 58 [NASA ADS] [CrossRef] [Google Scholar]
  58. Silsbee, K., & Rafikov, R. R. 2015b, ApJ, 798, 71 [NASA ADS] [CrossRef] [Google Scholar]
  59. Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55 [Google Scholar]
  60. Smoluchowski,M. V. 1916, Zeitsch. Phys., 17, 557 [Google Scholar]
  61. Spaute, D., Weidenschilling, S. J., Davis, D. R., & Marzari, F. 1991, Icarus, 92, 147 [NASA ADS] [CrossRef] [Google Scholar]
  62. Statler, T. S. 2001, AJ, 122, 2257 [NASA ADS] [CrossRef] [Google Scholar]
  63. Stewart, G. R., & Ida, S. 2000, Icarus, 143, 28 [NASA ADS] [CrossRef] [Google Scholar]
  64. Stewart, S. T., & Leinhardt, Z. M. 2009, ApJ, 691, L133 [NASA ADS] [CrossRef] [Google Scholar]
  65. Stewart, G. R., & Wetherill, G. W. 1988, Icarus, 74, 542 [NASA ADS] [CrossRef] [Google Scholar]
  66. Tanaka, H., Inaba, S., & Nakazawa, K. 1996, icarus, 123, 450 [NASA ADS] [CrossRef] [Google Scholar]
  67. Thebault, P. 2011, Celest. Mech. Dyn. Astron., 111, 29 [NASA ADS] [CrossRef] [Google Scholar]
  68. Thebault, P., Kral, Q., & Olofsson, J. 2021, A&A, 646, A173 [CrossRef] [EDP Sciences] [Google Scholar]
  69. Thébault, P., Marzari, F., & Scholl, H. 2006, Icarus, 183, 193 [NASA ADS] [CrossRef] [Google Scholar]
  70. Thébault, P., Marzari, F., & Scholl, H. 2008, MNRAS, 388, 1528 [NASA ADS] [CrossRef] [Google Scholar]
  71. Thébault, P., Marzari, F., & Scholl, H. 2009, MNRAS, 393, L21 [NASA ADS] [CrossRef] [Google Scholar]
  72. Trubnikov, B. A. 1971, Sov. Phys. Dokl., 16, 124 [Google Scholar]
  73. Ward, W. R. 1981, Icarus, 47, 234 [Google Scholar]
  74. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
  75. Wetherill, G. W., & Stewart, G. R. 1989, Icarus, 77, 330 [Google Scholar]
  76. Wetherill, G. W., & Stewart, G. R. 1993, Icarus, 106, 190 [Google Scholar]
  77. Windmark, F., Birnstiel, T., Ormel, C. W., & Dullemond, C. P. 2012, A&A, 544, A16 [CrossRef] [EDP Sciences] [Google Scholar]
  78. Zsom, A., Sándor, Z., & Dullemond, C. P. 2011, A&A, 527, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Zurlo, A., Cieza, L. A., Pérez, S., et al. 2020, MNRAS, 496, 5089 [Google Scholar]

1

The damping timescale is defined such that the decay of planetesimal eccentricity vector e due to drag is described by ; it should be noted that τd ∝|eeg|−1 for the quadratic gas drag law.

2

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).

3

Unlike vcA−1, neither er nor collisional velocity diverge at the resonance (where A → 0), since the product vcA in Eq. (4) remains finite as A → 0.

4

We note that Fig. 9 does not reflect these changes well: The vcoll 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 dmin 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 vcoll across all planetesimal sizes.

5

By analogy with hyperthermophiles – extremophilic organisms thriving in extremely hot environments.

6

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

Table 1

Parameters of our model runs.

Table C.1

Variation of numerical parameters.

All Figures

thumbnail Fig. 1

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

In the text
thumbnail Fig. 2

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

In the text
thumbnail Fig. 3

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

In the text
thumbnail Fig. 4

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

In the text
thumbnail 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
thumbnail Fig. 6

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

In the text
thumbnail Fig. 7

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

In the text
thumbnail Fig. 8

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

In the text
thumbnail Fig. 9

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

In the text
thumbnail Fig. 10

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

In the text
thumbnail 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 dinitdmin = 2, to be compared with Fig. 8c (see Sect. 6.7 for details). (b) Fiducial model computed for dinitdmin just slightly exceeding unity, to be compared with Fig. 8c,d (see Sect. 7.1 for details).

In the text
thumbnail Fig. 12

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

In the text
thumbnail Fig. 13

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

In the text
thumbnail Fig. B.1

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

In the text
thumbnail Fig. B.2

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

In the text
thumbnail Fig. B.3

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

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.