Planet formation in stellar binaries: Global simulations of planetesimal growth

Planet formation around one component of a tight, eccentric binary system such as $\gamma$ Cephei (with semimajor axis around 20 AU) is theoretically challenging because of destructive high-velocity collisions between planetesimals. Despite this fragmentation barrier, planets are known to exist in such orbital configurations. Here we present a novel numerical framework for carrying out multi-annulus coagulation-fragmentation calculations of planetesimal growth, which fully accounts for the specifics of planetesimal dynamics in binaries, details of planetesimal collision outcomes, and the radial transport of solids in the disk due to the gas drag-driven inspiral. Our dynamical inputs properly incorporate the gravitational effects of both the eccentric stellar companion and the massive non-axisymmetric protoplanetary disk in which planetesimals reside, as well as gas drag. We identify a set of disk parameters that lead to successful planetesimal growth in systems such as $\gamma$ Cephei or $\alpha$ Centauri starting from $1-10$ km size objects. We identify the apsidal alignment of a protoplanetary disk with the binary orbit as one of the critical conditions for successful planetesimal growth: It naturally leads to the emergence of a dynamically quiet location in the disk (as long as the disk eccentricity is of order several percent), where favorable conditions for planetesimal growth exist. Accounting for the gravitational effect of a protoplanetary disk plays a key role in arriving at this conclusion, in agreement with our previous results. These findings lend support to the streaming instability as the mechanism of planetesimal formation. They provide important insights for theories of planet formation around both binary and single stars, as well as for the hydrodynamic simulations of protoplanetary disks in binaries (for which we identify a set of key diagnostics to verify).


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 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 axisymmet-Article number, page 1 of 29 arXiv:2107.11389v1 [astro-ph.EP] 23 Jul 2021 A&A proofs: manuscript no. finalPaper ric, 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,  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,  extended the model of  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  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 by coagulation).
Despite this work, some uncertainty remains regarding the collisional environments that would allow planetesimal growth. In particular,  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 for unimpeded 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  gives rise to a dynamical environment with a mixture of low-and highvelocity collisions, making it challenging, even qualitatively, to determine the evolution of the size distribution.
The problem is further complicated by the radial inspiral of planetesimals due to interaction with the gas. The gas disk is slightly 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  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 coagulationfragmentation code for accurately following the evolution of the planetesimal size spectrum at every radius in the disk. In addition to that, our model directly follows the exchange of mass in solid objects between the different radial locations in the disk. These components of the model use the collisional velocities and radial inspiral rates calculated in Rafikov & Silsbee (2015a,b), which fully account for the secular gravitational effects of both the stellar companion and the massive eccentric disk in which planetesimals orbit. Outcomes of planetesimal collisions are determined using the prescriptions derived in Stewart & Leinhardt (2009). With this physical model, we are able to determine whether the coagulation process as a whole is successful for a given environment (i.e., a disk model), rather than determining if it is successful based on the outcomes of individual collisions. Armed with this powerful tool, we then perform an extensive exploration of the parameter space of the disk models to determine the conditions under which planet formation is possible in compact eccentric binaries such as γ Cephei or α Centauri. This paper is organized as follows. After describing our model setup in Sect. 2, we provide an overview of the important aspects of planetesimal dynamics in 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.

Model setup
To illustrate our calculations, we consider a model S-type binary with parameters similar to that of the γ Cephei system (Hatzes et al. 2003;Chauvin et al. 2011). This binary consists of M p = 1.6M (primary) and M s = 0.4M (secondary) stars in an orbit with a semimajor axis a b = 20 AU and eccentricity e b = 0.4.
A planet with mass M pl = 1.6M J orbits the primary star of the γ Cephei system with the semimajor axis a pl = 2 AU and eccentricity e pl = 0.12, which we use to motivate our subsequent calculations. There are a number of other binaries with a b 20 AU that have similar characteristics (Chauvin et al. 2011;Marzari & Thebault 2019).
We assume the primary component of the binary to host a massive protoplanetary disk, coplanar with the binary orbit, with properties similar to the disk described in : 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: where a 0 is a reference semimajor axis. We assume the disk to be sharply truncated by the torque due to the secondary star at the outer boundary a out = 5 AU (i.e., Σ(a) = 0 for a > a out ). The apsidal angle of the gas streamlines in the disk with respect to the binary orbit is denoted d and, for simplicity, is assumed to be independent of a, that is, the disk as a whole is apsidally aligned. There is no particular reason for choosing e g (a) in the form (1), except that it scales as the forced eccentricity due to the secondary. Development of a more accurate, better motivated model for the e g (a) behavior is beyond the scope of this study. Throughout this paper, we assume a solid-to-gas ratio of 0.01, and the disk to be slightly flared with the h/r profile given by Numerical calculations self-consistently following growth of planetesimals in protoplanetary disks from small sizes, and often all the way into the planetary regime, have been routinely used to understand planet formation around single stars. A standard approach is to follow the evolution of the size (mass) spectrum of colliding objects by solving the so-called Smoluchowski (Smoluchowski 1916) or coagulation equation (Safronov 1972;Wetherill & Stewart 1989;Kenyon & Luu 1998), often with account being given to the possibility of collisional fragmentation (Wetherill & Stewart 1993;Kenyon & Luu 1999). Global coagulation simulations, incorporating radial drift of mass in the disk driven by gas drag (Spaute et al. 1991;Inaba et al. 2001;Kenyon & Bromley 2002;Kobayashi et al. 2010), represent more sophisticated versions of such calculations. Some studies have also included direct N-body integration of a modest number (∼ 10 3 ) of planetary embryos formed as a result of planetesimal coagulation (Bromley & Kenyon 2011), and accounted for the most recent results on the outcomes of planetesimal collisions (Stewart & Leinhardt 2009). However, probably the most important ingredient of such coagulation-fragmentation calculations is the proper description of planetesimal dynamics (Stewart & Wetherill 1988;Stewart & Ida 2000;Rafikov 2003aRafikov ,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.

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.  computed the potential due to an apsidally aligned disk with the surface density and eccentricity profiles given by Eq.
(1); their calculation was subsequently generalized for a disk with arbitrary profiles of Σ(a) and e g (a) by Davydenkova & Rafikov (2018).
Including disk potential, the secular disturbing function characterizing the perturbation of planetesimal motion by the external potential becomes ) where a, n, and are the semimajor axis, mean motion, and apsidal angle of planetesimal orbit (both and the disk apsidal angle, d , are measured relative to the apsidal line of the binary), B b and B d describe the torques due to the non-axisymmetric components of the potentials of the secondary and the disk, and A = A b + A d is the combined free precession rate due to the secondary (A b ) and the disk (A d ). The explicit expressions for A d , B d , A b and B b can be found in  (see their Eqs. [5]-[6] and [9]-[10], respectively). We note that while A d depends only on the Σ(a) profile, B d is determined by both Σ(a) and e g (a).
Of particular importance for planet formation in S-type binaries is the fact that the secondary companion always drives prograde apsidal precession of the planetesimal orbit (A b > 0), while the precession due to disk gravity is usually retrograde (A d < 0). Since A b and A d vary with a in different ways, one often finds that A = A b + A d → 0 at a particular semimajor axis in the disk, if the disk is sufficiently massive ( 10 −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 were previously uncovered in the context of planet formation in binaries in Rafikov (2013a), Meschiari (2014), Silsbee & Rafikov (2015b,a), and  and for massive disks with planets in Heppenheimer (1980), Ward (1981), and Sefilian et al. (2021).
In addition to the secular perturbations we have considered, the companion star also produces short-period perturbations which act on the orbital period of the planetesimals (Murray & Dermott 1999). The magnitude of the short-period perturbations is quite small for bodies which are close together. If we consider two bodies with relative eccentricity, e r , which will collide in the next orbit, then their distance is of order e r a p . Their relative acceleration due to the tidal perturbation of the secondary is GM s e r a p /R 3 s . If this acceleration acts over a time ∼ n −1 p it results in a relative velocity e r v K (a p /R s ) 3 , which is suppressed by a factor (a p /R s ) 3 relative to their typical relative velocity e r v K in the absence of such short-period perturbations. Over longer timescales such short-period perturbations average out to zero.
Because the disk we consider is small compared to the binary orbit, the resonant perturbations are also not important. Indeed, Article number, page 3 of 29 A&A proofs: manuscript no. finalPaper the only resonances present in the disk are of very high order, and therefore very narrow and unlikely to significantly affect the dynamics.
The torque due to the secondary is expected to sharply (although not as discontinuously as we assume in Eq. (1)) truncate disk density at the semimajor axis a out . We take this truncation into account when computing disk gravity. The disk gravityrelated coefficients A d and B d diverge logarithmically when a sharp disk edge is approached 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 , provided that the inner boundary is located at no more than half the semimajor axis of interest, it makes no significant contribution to the disk disturbing function R.

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  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.  found that planetesimal dynamics admits two distinct regimes, depending on the degree of coupling between the planetesimal and the gas. Strongly coupled (small) planetesimals follow orbits nearly aligned with the local eccentric gas streamline. Weakly coupled (large) planetesimals have orbits nearly aligned with the forced eccentricity determined by the gravitational perturbations due to both the companion and the disk. Separating these two regimes is the characteristic planetesimal size d c that is defined as the planetesimal size at which the damping timescale 1 τ d of the planetesimal eccentricity (relative to the eccentricity vector of the local gas streamline e g = e g (cos d , sin d )) is of order A −1 , where A is the precession rate of free eccentricity in the combined potential of the disk and perturbing star (see Eq. (2)). More specifically,  show that for a planetesimal of radius d p (see their Eq. [32]), so that |Aτ d | ≈ 1.27 for d p = d c .  also show that the relative planetesimal-gas eccentricity e r = |e − e g | is given by (see their Eq. [28]) where v c = e c v K is the characteristic planetesimal-gas velocity for planetesimals with d p ≥ d c : e c is the magnitude of the relative eccentricity between the gas streamlines and the forced eccentricity from the gravitational potential. An explicit expression 1 The damping timescale is defined such that the decay of planetesimal eccentricity vector e due to drag is described byė drag = −(e − e g )/τ d ; it should be noted that τ d ∝ |e − e g | −1 for the quadratic gas drag law. for different disk models, in which we vary disk eccentricity e 0 , surface density normalization Σ 0 , and apsidal angle d (one parameter at a time relative to the fiducial model described in Sect. 5 and shown with the thick black curve), as indicated in the legend. See the text for more details on various features of these curves.
for e c in terms of the disk and binary parameters only is provided by Eq.
[29] of . Knowing e c one can then directly compute d c through Eq.
[31] of . Both d c and v c (or e c ) are the key metrics of planetesimal dynamics in S-type systems and will frequently appear in this work. Using these two parameters, strong coupling corresponds to small bodies, d p d c , when |Aτ d | ≈ (d p /d c ) 1/2 1 and e r ≈ e c (d p /d c ) 1/2 e c . Weak coupling is realized in the opposite limit of d p d c , when |Aτ d | ≈ d p /d c 1 and e r ≈ e c . This suggests another intuitive definition of v c (or e c ) as the characteristic relative velocity (or eccentricity) between the weakly and strongly coupled (or large and small) objects.
The parameters v c and d c vary significantly throughout the disk, and with different disk models, which is illustrated in Fig.  1, where we vary (one at a time) the values of the disk eccentricity e 0 , surface density normalization Σ 0 and disk orientation d . One can see that the behavior of v c and d c directly reflects many peculiar features of planetesimal dynamics that emerge when both secondary and disk gravity are important, as described in Sect. 3.1. In particular, (i) d c and v c diverge at 2.5 AU in a model Article number, page 4 of 29 Kedron Silsbee & Roman R. Rafikov: Global simulations of planetesimal growth in stellar binaries 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 a out ); (ii) many models, in which the disk is apsidally aligned with the binary, feature a dynamically quiet location where d c , v c → 0 (with the semimajor axis varying with Σ 0 and e 0 ); (iii) dynamically quiet locations do not appear for all values of e 0 and Σ 0 or for apsidally misaligned disks (e.g., for d = 90 • ). These observations will greatly help us in interpreting the results of our detailed calculations in Sects. 5 and 6.
The concept of a dynamically quiet location in the disk will be very important for understanding the results of this work.  have demonstrated (see their Eq. [50]) that in an apsidally aligned disk, in the presence of gas drag, d c , e c → 0 at the semimajor axis a dq where 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  and used in this work assumes that planetesimal eccentricities are at their steady state values, given by Eqs.
[22]- [27] in , to which they converge on a timescale ∼ τ d due to gas drag. This is a reasonable assumption since for small objects τ d is short, 10 3 yr for d p = 1 km objects, while for bigger bodies (for which τ d ∝ d p is longer) the timescale for collisions with comparable objects (capable of significantly perturbing eccentricity) is long.

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,  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  as described in Sect. 3.2, plus a random component. We take this random component of the eccentricity vector to have standard deviation σ e which is twice σ i (Stewart & Ida 2000), but allow horizontal and vertical motions to be independent of one another. The degree of random motion can therefore be fully parameterized by σ i .
If random motions result from viscous stirring, one would typically expect them to be on the order of the escape velocity from the planetesimals. However, in the present perturbed system, because the relative motions between planetesimals are much larger than their random motions, the excitation from twobody encounters is lowered and the damping due to gas drag is enhanced. On the other hand, if present, disk turbulence could excite larger random motions. A study of the random planetesimal motions in a binary system would be subject to many uncertainties, and we do not attempt it here. We note that our typical random motions correspond to the escape velocity from roughly 10 km planetesimals.
Varying σ i has two effects. Let e i j be the relative eccentricity between two bodies in the absence of any random motions. For collisions between bodies where e i j σ i , varying σ i has little effect on the collision speed, but the frequency of such collisions varies inversely with σ i , as σ i sets the thickness of the planetesimal disk and therefore the planetesimal number density. In the limit that e i j σ i , the collision velocity is proportional to σ i and the collision rate becomes independent of σ i provided that for colliding objects with sizes d 1 , d 2 and masses m 1 , m 2 ), that is, gravitational focusing is negligible; the collision rate scales as σ −2 i in the opposite limit σ i v K v esc when gravitational focusing enhances the collision rate. The effect of random velocities on the collision velocities and rates and their implementation in our numerical framework are discussed in detail in Sects. A.3 and A.4. Figure 2 illustrates the collision velocities v coll between planetesimals as a function of their sizes. These are given by Eq. (A.16), with λ and e rand assigned their mean values of 0.81 and 2 √ πσ i , respectively. This figure is made for the fiducial system described in Sect. 5 with rather small σ i = 10 −4 , so that relative velocities are typically dominated by the secondary and disk forcing (i.e., by e i j ). The three panels differ in the assumed location a p of the colliding planetesimals in the disk (which determines their v c and d c ), as labeled on the figure. As mentioned in Sect. 3.2, v coll ∼ v c for pairs of collision partners in which one object has size d 1 d c and another has size d 2 d c . If the collision partners have similar size, or both are much larger (or both much smaller) than d c , then the collision velocities are determined predominantly by random motions and become independent of the sizes of the collision partners (since we take σ i to be size independent in this work).
As mentioned above, Fig. 2 would look very different if it were drawn assuming planetesimal dynamics typical for disks around single stars (i.e., without apsidal alignment): in that case the "valley" at d 1 ≈ d 2 would disappear and v coll would monotonically increase along this line (since smaller planetesimals have lower random velocities as a result of gas drag, Rafikov 2003aRafikov , 2004. These differences clearly demonstrate the importance of apsidal alignment naturally emerging in disks around binaries for setting the pattern of planetesimal collision velocities.

Collision outcomes
Armed with the understanding of the distribution of relative velocities of planetesimals (Sect. 3.3), we determine the outcome of 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 Article number, page 5 of 29 A&A proofs: manuscript no. finalPaper 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 , we divide the collision outcomes into three classes. We say that a collision leads to catastrophic disruption if the largest remnant contains less than half the combined mass of the two incoming planetesimals, erosion if the largest remnant is smaller than the larger of the two incoming bodies, and growth if the largest remnant is larger than either of the incoming bodies. We can then use the recipe of Stewart & Leinhardt (2009) to delineate in Fig. 2 the domains in d 1 − d 2 phase space, corresponding to each to type of the outcome: regions of catastrophic disruption lie within the solid white lines, regions of erosive collisions are outlined in dotted white, and regions leading to growth lie outside both white lines. We use the "strong rock" planetesimal composition of Stewart & Leinhardt (2009) in this calculation. One can see that, unless v c is very high, catastrophic disruption occurs only between collision partners with similar (but not equal) sizes. It is mainly relevant for objects with sizes around d c , as that is where even objects of comparable size experience high-velocity collisions (although the regions of catastrophic disruption are not exactly centered on d c because the strength of planetesimals is size-dependent, with a minimum around 0.1 km (Stewart & Leinhardt 2009). In contrast to catastrophic disruption, erosion is possible even when the sizes of the colliding bodies are very unequal.
Due to the complexity of secular dynamics in binaries outlined in Sect. 3.2, there is a significant variation in the collisional environment across the disk, as different panels of Fig. 2 demonstrate. The regions corresponding to erosion and catastrophic disruption grow larger as the characteristic velocity v c increases. The weakest planetesimals with sizes around 0.1 km may be destroyed even in the absence of any secularly excited eccentricities, in collisions arising just from random motions, which is almost what is happening in Fig. 2a. But it is still clear from panels (b) and (c) of that figure that secular pumping of planetesimal eccentricities in binaries endows these systems with much harsher collisional environments than in disks around single stars.
The heuristic picture of collision outcomes based on the relative velocity maps has been used in  and  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), Thebault (2011, and  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.

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),  found that the inspiral rateȧ p is given by (their Eq. 11): where e r and τ d are given by Eqs. (4) and (3), respectively, η ∼ (h/a p ) 2 ∼ 0.003 (h is the disk scale height) is a measure of the particle-gas velocity differential due to the pressure support in the gas disk, and α ∼ 1 is a constant that depends on the disk model as described in . From Eqs. (3) and (4), one can show ) that e r τ d ∝ d p .
As a result, one finds thatȧ p ∝ d −1 p if either (i) d p d c , when e r becomes constant, or (ii) e r η, so that the last term in brackets in Eq. (6) dominates.  Figure 3 shows the speed of radial drift for different planetesimal sizes and locations in the disk in the γ Cephei system, assuming fiducial disk parameters as in Sect. 5. We see that the curves in the upper panel match the limiting cases discussed in the previous paragraph; the regions whereȧ p ∝ d −1 p are joined by an intermediate region where the inspiral velocity depends less strongly on d p . In all cases, inspiral speed is a monotonically decreasing function of d p , meaning that small bodies are preferentially flushed out of the system, thus reducing the erosion suffered by larger objects.
In the lower panel we see that for moderately large planetesimals ( 0.1 km), there is a noticeable dip in the inspiral rate in the broad region around 1.8 AU. This is because for this disk model v c vanishes at 1.8 AU, so the only contribution to the inspiral comes from the sub-Keplerian rotation of the gas and not from any relative eccentricity between the planetesimals and gas, significantly reducingȧ p .
These considerations highlight important differences in the radial drift behavior between the binary and single stars. First, high planetesimal eccentricities driven by the secular effect of the disk and the secondary result in faster radial drift in disks in binaries compared to single systems (see Eq. (6)). Second, because higher e r increases the inspiral rate, planetesimals will be naturally concentrated in regions where e r is low, and thus where they can grow more easily. This effect does not exist for planetesimals around single stars (unless one includes additional physics).

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.

Basics of the code structure
Our code models the evolution of the planetesimal size distribution in N ann discrete spatial annuli placed at different radii from the primary. These should really be thought of as bins in the space of the semimajor axis, rather than radius, but we do not make that distinction in the following discussion. Within each annulus, the planetesimal size distribution is represented using N bins logarithmically spaced mass bins. At any point in time, in a given annulus, all of the information about the size distribution is encoded in a single vector n, such that n i is the number of planetesimals in mass bin i. To follow the evolution of n i in each annulus the code performs two basic functions.
First, in each of the annuli, it calculates the evolution of the size distribution due to planetesimal-planetesimal collisions. This is done in a standard fashion, effectively by solving a discretized version of the Smoluchowski equation (Smoluchowski 1916) accounting for the possibility of fragmentation (Sects. A.1, A.2, and A.6.1). Inclusion of fragmentation generally makes this an O(N 3 bins ) calculation at each time step, which is numerically expensive. To overcome this problem, our code employs the new fragmentation algorithm developed in Rafikov et al. (2020), for which the numerical cost goes only as O(N 2 bins ), 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  and accounting for both the forced eccentricity and random velocities (see Sects. 3.2 and 3.3). Our collision rates include gravitational focusing and smoothly interpolate between the shear-and dispersiondominated velocity regimes (Sect. A.4). The full distribution of collision velocities is also used to model collision outcomes following the recipes in Stewart & Leinhardt (2009)

(see Sects. 3.4 and A.2).
Our implementation of many code components has certain elements of stochasticity in it, which is quite important. It has been previously found (Windmark et al. 2012;Garaud et al. 2013) that including the distribution of collision velocities can qualitatively change the outcome of the coagulation process. In particular, in a study of dust growth, Windmark et al. (2012) found that using a Maxwellian collision velocity distribution instead of a delta function at the rms velocity of the Maxwellian can cause a few particles to experience a series of lucky lowvelocity collisions and grow to larger sizes; the larger bodies Article number, page 7 of 29 A&A proofs: manuscript no. finalPaper 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 see Sect. 3.3), as described in Sect. A.3.
The second operation performed by the code is the radial redistribution of mass between different annuli caused by the sizedependent inward migration of planetesimals due to gas drag (Sect. 3.5). The implementation of this procedure and its tests are described in Sects. A.5, A.6.2, and B.3, respectively.
Our model implicitly assumes that other than redistribution caused by radial drift, planetesimals interact only with other planetesimals within their semimajor axis bin. In practice this is not strictly true: Planetesimals in neighboring bins can also collide. However, because the growth we are modeling occurs in regions where the relative eccentricity between planetesimals of different sizes is very small (typically less than 1%), we expect the coupling of planetesimals with those in adjacent semimajor axis bins to be a minor effect.

Parameters of the numerical model
Multiple modules comprising our code employ a number of parameters, both numerical and physical. All our simulations use a standard set of these parameters described here and listed in Table C.1.
Our grid in mass space employs N bins = 100 bins; the smallest size (corresponding to the lowest bin) below which mass is supposed to be flushed out from the system due to gas drag is d p = 10 m. The largest mass bin in our simulations corresponds to an object with d p = 300 km, since beyond this point destructive collisions are no longer able to prevent growth to larger sizes. This results in a mass ratio between (logarithmically uniformly spaced) mass bins of M i+1 /M i = 1.37. To represent our global simulation domain extending from 1 AU to 4 AU we use N ann = 60 annuli, spaced in semimajor axis as described in Sect. A.5. In choosing this radial range we are motivated by the planetary semimajor axes in several S-type systems: a pl = 2 AU in γ Cephei, a pl = 2.6 AU in HD 196885, and a pl = 1.6 AU in HD 41004 (Chauvin et al. 2011).
Some additional parameters used in our simulations are as follows: the slope of the fragment mass spectrum ξ = −1 (Sect. 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 time step 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).

Results: Fiducial model
We start presentation of our results by describing in detail a particular (fiducial) simulation, which will also illustrate the metric we employ to determine the success of planet formation in a given disk model (Sect. 5.2.1). Our fiducial disk model uses the orbital parameters of the γ Cephei system and disk characteristics as described in Sect. 2. For the disk and planetesimal  parameters it adopts e 0 = 0.01, Σ 0 = 1000 g cm −2 at a 0 = 1AU, d = 0 • , σ i = 10 −4 , and a solid to gas ratio of 0.01. The total disk mass is 3.7M J . We present the results for other disk models in Sect. 6.
To assist in interpreting the outcomes of our calculations, we first describe the results of several 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).

Coagulation in a single annulus
We ran our single-annulus calculations of planetesimal growth at three semimajor axesa p = 2 AU, 2.3 AU, and 2.5 AU -in our fiducial disk model. These are the same locations for which Fig. us in understanding the role of planetesimal dynamics on their growth. The simulations started with all planetesimals having the same size d init = 4.3 km; for the adopted bulk density ρ p = 3 g cm −3 this corresponds to planetesimal mass of m = 10 18 g. Figure 4 illustrates how the size distribution in our model system changes as a function of time. What is actually shown is mdN/d ln d p , a variable equal to the mass per unit log planetesimal size. Different curves correspond to different times since the beginning of the simulation, as indicated by the color bar. The simulation for panel (a) was stopped when the largest body went beyond the largest mass bin in the simulation (d p = 300 km). In panels (b) and (c), the simulation was stopped after 10 6 yr.
In all plots, the curves corresponding to times 10 4 years display features due to the initial conditions and finite spacing of the mass bins. Initially, all of the mass is distributed in just two mass bins surrounding the initial (seed) size d init , which is marked by the dashed vertical black line. The wiggles just to the right of d init are related to the finite spacing of the mass bins. The gap just to the left of d init exists because debris of that size is not created in collisions of objects with size d init -they produce only smaller fragments to the left of the gap. The size spectrum of objects to the left from the gap is universal at early times (∼ 10 3 −10 4 yr), as it simply reflects the size distribution of fragments formed in collisions of roughly equal mass planetesimals of size d init . The gap region is only later filled in by coagulation of the debris particles produced early on and erosion of the seed bodies by small debris.
Over time, memory of the initial conditions is erased, the gaps near d init get filled, and the size distribution becomes smooth. There are some other notable features that develop in the size distribution at late times. First, after a few times 10 4 years mdN/d ln d p starts featuring a dip around d p ∼ 0.1 km, especially pronounced in panel (b). Its location corresponds to the planetesimal size which is most easily destroyed in collisions (see . In each panel, the vertical dotted line corresponds to d c for that simulation, which varies by a factor of several across the different environments, but stays pretty close to d p ∼ 0.1 km. Particles with sizes around d c experience high-velocity collisions with almost all other bodies, and are therefore preferentially destroyed, which likely affects the depth of the dip in different panels. Second, one can see (especially in panel (b)) mdN/d ln d p developing power law behavior to the left and right of the d p ∼ 0.1 km line, with different slopes. We believe this behavior to be real but will refrain from offering an explanation for the slope of these segments. This cannot be done based on standard results regarding fragmentation (O'Brien & Greenberg 2003;Tanaka et al. 1996) since they assume a power-law dependence of the collision rate on colliding partner size. Careful examination of our simulation outputs shows that the largest objects grow primarily in collisions with objects of similar size, certainly for d p > d init . This is because most solid mass is contained in such objects, but also because growth is most favorable in collisions of such objects.
We observe, not surprisingly, that higher v c suppresses growth. In panel (a), with v c = 77 m/s, growth to hundreds of km occurs within several ×10 5 years. This can also be seen in Fig. 5 where we show the size of the largest body as a function of time for the three local simulations discussed here. Were we to simulate larger bodies, beyond 300 km, growth would continue in panel (a) until the mass reservoir is fully depleted. Because of low v c , even when a significant amount of solid mass in the system has been converted to debris (at late times), collisions of the largest objects with debris particles are not energetic enough to significantly impede their growth.
In panel (b) with v c = 265 m/s, planetesimals can only grow to a few times the initial size. After that growth of the largest objects stalls (see also Fig. 5), while the population of objects with d p ∼ d init gets gradually eroded by smaller objects produced in previous collisions (this can be seen in the decay of amplitude of dN/d ln d p around that size). If we were to run this simulation beyond 1 Myr, large objects with d p ∼ d init would eventually be eroded, followed by smaller bodies also grinding themselves down.
Finally, in panel (c) of Fig. 4, collision velocities are so high that there is barely any growth -the largest objects present in the system at early time have grown in size beyond d init only by ∼ 2. But after a population of small debris develops in the system (already by 10 3 yr), Fig. 5 shows that large bodies get eroded away by 10 5 yr, leaving only a small population of debris near the lower end of the simulation mass range (which is an artifact of our calculation, as we mentioned earlier). Figure 5 shows that in all three environments, the largest bodies initially evolve in a very similar fashion, steadily growing in size. This universality is caused by the fact that initial growth occurs by mergers of objects with d p = d init , which have low relative velocities. The differences emerge only after ∼ 10 4 yr when a sufficient amount of small debris capable of efficiently eroding big bodies in high-v c environments accumulates in the system. Article number, page 9 of 29 A&A proofs: manuscript no. finalPaper Here we again consider Fig. 2, which illustrates, in particular, collision outcomes as a function of the sizes of the collision partners, for the same environments that the simulations shown in Fig. 4 were run for. Interestingly, we see that to halt planetesimal growth it is not necessary to have catastrophic disruption of objects at d p ∼ d init . Indeed, although Fig. 4 shows that growth is halted in panels (B) and (C), Fig. 2 shows no catastrophic disruption of planetesimals at the initial size: the point (d init , d init ) is outside the solid white contours for all dynamic environments, as d init = 4.3 km considerably exceeds d c in all panels.
It may seem strange that in all three panels in Fig. 2 the point (d init , d init ) also lies outside the dotted contours, bordering the region of erosive encounters -it would seem natural that one should then find growth starting with a population of d p = d init objects in all three environments. The resolution of this apparent paradox involves two factors. First, according to our definition of erosion in Fig. 2 dotted contours correspond to collisions, in which the mass of the largest remnant is equal to the mass of the bigger object (target) involved in the collision; the mass equal to the mass of the smaller object (projectile) is reduced to debris. Similarly, even though the point (d init , d init ) sits outside the dotted contour, some amount of debris will still be produced in collisions of two objects of initial size d init ; the amount of created debris is larger the closer this point is to the dotted contour (i.e., in panels (b) and (c) relative to (a)). And once some debris particles are created, Fig. 2 shows that their collisions with seed bodies (d p = d init ) are very erosive, especially in panel (c), but only barely so in panel (a).
Second, relative velocity maps in Fig. 2 assume a particular (mean) value of encounter velocity, whereas in practice v coll has a reasonably broad distribution around this mean (see Sect. A.3). As a result, some collisions occur at higher velocities and are more destructive -a fact that is fully accounted for in our coagulation-fragmentation simulations depicted in Fig. 4. This is another clear illustration of the importance of considering the distribution of planetesimal velocities when studying their growth. The combination of the two aforementioned factors leads to erosion eventually becoming the dominant player in the collisional evolution in panels (b) and (c) of Fig. 4.

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

Outcome diagnostic
In any reasonable environment, planetesimal growth would proceed to very large sizes, eventually leading to planet formation, provided that the initial planetesimal size d init is large enough. Indeed, objects with sizes d p d c do not experience catastrophic disruption and get eroded only by much smaller planetesimals, which cannot compete with the addition of mass in collisions with larger bodies. Thus, given large enough d init , planet formation is guaranteed to be successful (e.g., Thebault 2011) even in dynamically harsh environments of the binary stars.
However, in many environments, if the initial planetesimals are too small, then they will be eroded or destroyed and large bodies will not form (like in panels (b) and (c) of Fig. 4). For this reason, the appropriate question to ask is not whether planetesimal growth can occur, as the answer would depend on the initial condition -the size of the seed planetesimals d init . Instead, one should be trying to figure out from what initial planetesimal size d init can sustained planetesimal growth occur. In a given dynamical environment, there will be a minimum size d min such that if the initial bodies are smaller than d min then planetesimal growth in the system will eventually be halted, as in Fig. 4b,c, whereas if seed objects have d init > d min , then planetesimal coagulation will eventually form large bodies, thus overcoming the fragmentation barrier.
This question, in principle, can be asked locally, at a given semimajor axis in the disk. In this case Fig. 4 makes it clear that, for a fiducial disk model, d min < 4.3 km at a p = 2 AU, whereas d min > 4.3 km at a p = 2.3 and 2.5 AU. However, the calculations presented in Sect. 5.1 do not account for radial mass transport due to gas drag and may thus be inaccurate. For this reason, we will be more interested in a general question of what d min is needed for the fragmentation barrier to be overcome and for a planet to form at some location in the whole disk, given its parameters. Once d min is determined, some other information, for example, the location where planetesimal growth is fastest, would be the outcome of our global calculations, giving them certain predictive power. Of course, d init may (and probably should) vary across the disk. Moreover, at each location, a whole spectrum of initial sizes is likely to be present. However, in this work, to reduce the number of degrees of freedom, we assume that seed planetesimals have the same size d init across the whole disk and determine a single value of d min for a given disk model based on this assumption.
To determine d min defined in this manner, we employed the following iterative algorithm. We chose a value of d min and used that as our starting size d init . The simulation was run until one of three conditions were met: (i) 1 Myr has passed, (ii) a 300 km body is formed, or (iii) there is no longer sufficient solid mass remaining in the simulation to produce a 300 km body. If a 300 km body forms, then we know d min < d init , whereas if it does not, then d min > d init . By trying several values of d init and running our global simulations for each of them, we were able to first bracket d min within some interval, and then to converge to its value with high accuracy. This procedure is demonstrated next (allowing us to also illustrate the sensitivity of global planetesimal growth to d init ).
Since the planetesimal composition is uncertain, in addition to our default simulations that use the material properties of "strong rock" bodies from Stewart & Leinhardt (2009), we also ran simulations using the strength parameters for their "rubble pile" planetesimals. We denote the minimum size assuming rubble pile planetesimals as d rubble min . Figure 6 shows the evolution of the planetesimal size distribution in the entire disk for several runs of our fiducial simulation with different initial sizes d init expressed in terms of d min = 1.6 km, as determined for this disk model. Curves of different color correspond to different times since the simulation was started. Because collisional processes proceed at different rates across the disk, reflecting the complexity of planetesimal dynamics in binaries (see Fig. 1), the evolution of the planetesimal mass spectrum looks somewhat different from that in a single annulus (see Fig. 4; e.g., the initial gaps to the left of d init are barely there, and there is no pile-up at small sizes). This is to be expected since the global size spectrum in Fig. 6 is an average of multiple singlezone size distributions (such as the ones shown in Fig. 4) with the added complication of radial inspiral of solids. Nevertheless, some key features of the size distribution evolution (e.g., a dip in mdN/d ln d p near d p ∼ 0.1 km, mass flow toward small sizes, and so on) are still preserved even in the full disk calculations.

Evolution of the system as a whole
We can see that in panels (a) and (b), essentially no growth occurs beyond the initial size. The behavior shown in these panels differs only in how rapidly the bodies are ground down. The lack of any appreciable growth in panel (b) with initial size only a factor of two less than d min shows that there is a sharp transition over less than a factor of two in d init from growth to large sizes to essentially no growth whatsoever.
In panels (c) and (d), coagulation to 300 km does occur. In both cases, though, less than half of the initial solid mass in the simulation remains at the end. The time required to produce a 300 km body is almost four times less for the simulation with d init /d min = 5 than for the one with d init /d min = 2 (8 × 10 4 years vs. 3×10 5 years), and in the former case more than twice as much solid material remains (47% of the original 7M ⊕ in panel (d) vs. 21% in panel (c)). In the simulation with d init /d min = 2 (5), 79% (88%) of the mass is lost to the creation of rubble smaller than the smallest body tracked in the simulation, with the remaining amount lost to inspiral into the central star. Figures 7 and 8 provide details on how the coagulationfragmentation process proceeds as a function of semimajor axis in the disk, for each global simulation shown in Fig. 6. Figure  7 shows dM/da -the total solid mass (in planetesimals of all sizes) per unit a, whereas Fig. 8 depicts the size of the largest body (in a given annulus), both as functions of a and time. In Fig.  7, the horizontal dashed violet line represents the initial mass distribution, which is flat because of our chosen Σ profile (see Eq. (1)). In Fig. 8 that line corresponds to the size of seed objects d init , which is constant across the disk according to our assumption (see Sect. 5.2.1).
To illustrate how the dynamics of planetesimals affect their size evolution, we also show in Fig. 9 the run of the collisional velocity v coll between d 1 = 2 and d 2 = 5 km planetesimals as a function of semimajor axis (for many of the different simulations carried out in this work, see Sect. 6). It accounts for both the secular and random velocity contributions and is calculated using whereas the black dotted line in panel (a) shows v coll computed in the absence of random motions. One can easily recognize the "valley" forming around 1.8 AU as the dynamically quiet location where v c , d c → 0 in Fig. 1. Nonzero random velocities of planetesimals "fill" this valley, not allowing v coll to drop to zero in their presence. Nevertheless, v coll is still relatively low around this region.
A large increase in v coll in the outer disk starting around 2 AU arises because for Σ 0 = 10 3 g cm −2 the free precession rate A becomes very small there, while the excitation by the companion grows large. As a result, both v c and v coll become very large, of order km s −1 , even though a true secular resonance does not appear in this disk model. We note that v coll is often several times smaller than v c , because d c is often less than d 1 = 2 km. Figure 7 shows that, in the cases with d init /d min < 1, most of the solid mass is removed from the outer and inner portions of the disk on timescales of 10 4 years, leaving the majority of the remaining solids in a narrow band around 1.8 AU. This band coincides with the dynamically quiet location in the fiducial disk model (see Fig. 9), where (i) the velocity of inward drift of planetesimals goes down (see Fig. 3) and (ii) v coll is relatively low. These factors naturally cause planetesimal accumulation (dM/da temporarily exceeds its initial value due to arrival of mass from larger radii) and promote growth at this location (see Fig. 8). But over time, the mass concentration around 1.8 AU is eroded and moves inward due to gas drag.
In the outer disk, outside 2 AU, both the amount of mass and the size of the largest object drops precipitously early on, primarily because of very high v coll ∼ 1 km s −1 in this region (see the black curve in Fig. 9) resulting in rapid planetesimal erosion down to the smallest size (10 m) tracked in our simulations 2 . Some of this debris gets swept into the dynamically quiet zone through inspiral before it is fully ground down. In the inner disk, below 1.8 AU, v coll is much lower (enabling some temporary growth for the largest objects) and mass is removed predominantly by gas drag.
In the simulations in which coagulation is successful (d init /d min > 1), the mass is again quickly removed from the outer part of the disk. But in the inner disk, particularly around 1.8 AU, the mass density (per unit semimajor axis) remains roughly constant as large objects are able to efficiently accrete smaller ones, preventing them from getting lost via inspiral or grinding down. In the run with d init /d min = 5 we find that in the outer disk, the largest bodies have size (2−3)d init for the duration of the simulation. This illustrates that larger seed objects find it much easier to survive than smaller ones even in the dynamically harsh environments.

Results: Variation of system parameters
Having discussed in detail the evolution of the global planetesimal size distribution for a particular disk model (Sect. 5), we proceed to explore the changes brought in by varying disk model parameters relative to their fiducial values. We vary several key physical parameters -surface density normalization Σ 0 , disk eccentricity normalization e 0 , disk orientation d , amplitude of random velocity σ i -which control planetesimal dynamics. We also run calculations in which we artificially turn off disk gravity (but not gas drag), and vary the inspiral rate, to explore the role played by these physical processes in determining the outcome of planetesimal evolution. In carrying out this parameter exploration we use the minimum planetesimal size d min (at which growth beyond 300 km becomes possible) as our metric to characterize the success of planet formation in a given model. Table 1 lists the parameters of the different models that have been explored, as well as their d min and d rubble min . 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).

Variation of surface density: Simulations 2-7
Somewhat unexpectedly, we find that variation in the surface density normalization Σ 0 within a broad range relative to the fiducial value of 10 3 g cm −2 has almost no effect on d min , which ranges between 1.1 and 1.7 km (see Fig. 10a and Table 1). Even though both the dynamically quiet location at which v c and d c vanish and the location (and existence) of the secular resonance 1. Disk surface density at a 0 = 1 AU 2. Disk eccentricity at a 0 = 1 AU 3. Angle between the apsidal line of the disk and that of the binary. 4. Dispersion of the distribution for each component of the planetesimal inclination vector. 5. If "yes", then disk gravity is included when calculating collision rates and outcomes, if "no", then it is ignored. 6. The inspiral rate is artificially multiplied by the inspiral rate multiplier χ ir . 7. Smallest initial planetesimals size which results in a 300 km body forming within 1 Myr. 8. Smallest initial planetesimal size which results in a 300 km body forming within 1 Myr, assuming rubble-pile planetesimals. change with Σ 0 (see Fig. 1), this does not dramatically affect d min . Figure 10a does show a weak decrease in d min as Σ 0 increases from 500 to 5000 g cm −2 . This is because, as in the fiducial model, the most favorable conditions for planetesimal growth occur at the dynamically quiet location, where v c and d c vanish, leaving only random velocity to produce relative velocity between the colliding planetesimals (i.e., v coll ∼ σ i v K ). As Σ 0 increases, the dynamically quiet location moves out in the disk. Since the Keplerian speed drops with distance, planetesimals collide at lower velocities at the more distant dynamically quiet locations (i.e., for higher Σ 0 ). As a result, within the range of Σ 0 where the dynamically quiet location exists, collisions are less destructive and planetesimal growth becomes possible starting at lower d init for larger Σ 0 .
We note that in the Σ 0 = 500 g cm 2 model a secular resonance ) appears in the disk around 2.4 AU (see Fig. 1), driving v coll to very high values 3 at moderate (2 -2.5 AU) semimajor axes. However, d min is only slightly larger than for other disk models (not featuring resonances), simply because most growth occurs at the dynamically quiet location, which is far from this resonance.
The situation changes for even more massive disks, Σ 0 > 5000 g cm −2 , in which the disk gravity is dominant over gravity from the companion throughout the whole disk. As a consequence, such disks have no dynamically quiet location 4 and we find a slight increase in d min with Σ 0 .
In Simulation 7 we also considered a model of a massive axisymmetric disk, with Σ 0 = 20,000 g cm −2 (corresponding to a total disk mass of 0.07 M interior to 5 AU) and e 0 = 0. In this setup, previously explored in Rafikov (2013b), disk gravity does not give rise to planetesimal eccentricity excitation but effectively suppresses the dynamical excitation due to the secondary's torque. As a result, in this model collision velocities are reduced even further (compared to, e.g., Simulation 6), and we find a significantly lower d min = 0.7 km (see Table 1).

Variation of e 0 : Simulations 8-12
Figure 10b shows that disk eccentricity has a dramatic effect on planetesimal growth: d min varies by 1.4 dex as the disk eccentricity normalization e 0 changes from 0 to 0.1. This variation is non-monotonic.
As e 0 increases from zero to the fiducial value e 0 = 0.01, we find d min decreasing for reasons similar to the ones mentioned in Sect. 6.1: Higher e 0 means larger |B d | ∝ e 0 Σ 0 (Silsbee & Rafikov 2015b; Rafikov & Silsbee 2015a), shifting the dynamically quiet location further from the primary (here we again assume d = 0 • ; see Fig. 1). As v coll ∼ σ i v K ∝ a −1/2 at this location, higher e 0 results in lower v coll (see Fig. 9c) promoting planetesimal growth and leading to lower d min . In an axisymmetric disk (e 0 = 0), B d = 0 and cancellation of the secondary's torque does not occur, leading to a rather high d min = 4.2 km. The dynamically quiet location exists in our simulation domain for values of e 0 between 0.004 and 0.018.
For e 0 > 0.018 the torque due to disk gravity starts to dominate planetesimal eccentricity excitation, and v c is nonzero throughout the whole simulated portion of the disk; the dynamically quiet zone disappears. In this regime v c ∝ |B d | ∝ e 0 (Silsbee & Rafikov 2015b), leading to higher v coll and a less growthfriendly environment as e 0 increases (see Fig. 9c); hence the rapidly increasing values of d min . Even a disk with e 0 = 0.05 requires seed planetesimals with d min ≈ 22 km to ensure their growth into planetary regime.
We note that while the location of the dynamically quiet region is sensitive to e 0 , variation of disk eccentricity does not affect the presence (or location) of the secular resonance in the disk. The latter is determined by the disk mass only, while the former depends on both Σ 0 and e 0 (Rafikov 2013b;  -the minimum size for collisionally weak planetesimals. The causes of the trends shown here are discussed in Sects. 6.1-6.7. In each panel, the fiducial simulation is highlighted in red.

Variation of d : Simulations 13-15
We also find substantial sensitivity of d min to d . Figure 10c shows that increasing the deviation of the disk orientation from apsidal alignment with the binary orbit leads to larger d min : Keeping everything else fixed, d min grows from 1.6 km for an aligned disk ( d = 0 • ) to 9.7 km for an anti-aligned disk ( d = 180 • ). Even a d = 10 • misalignment is sufficient to increase d min to 2.5 km.
This trend, again, owes its origin to the presence (or absence) of a dynamically quiet location in the disk. Secularly forced v c can vanish only in apsidally aligned disks (see Fig. 1a). A small misalignment makes v c nonzero everywhere (Rafikov & Silsbee 2015a,b) but only slightly modifies the behavior of v coll (see, e.g., the d = 10 • curve in Fig. 9d), because (i) v c is still dramatically reduced at the location where Eq. (5) is satisfied and (ii) random velocity σ i is nonzero. However, for a significant misalignment ( d ∼ 1), collisional velocities end up being much higher than in the aligned disk (see d = 90 • and d = 180 • curves in Fig.  9d), requiring larger seed planetesimals to overcome the fragmentation barrier. Therefore, it is reasonable that d min increases with the degree of misalignment, as our results show.

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

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

The case without disk gravity: Simulations 20-24
In order to better understand the impact of properly accounting for disk gravity and to isolate its effect on planetesimal growth, we also ran simulations where disk gravity was artificially turned off in our secular solutions for planetesimal eccentricity, that is, we set B d = A d = 0. In this case, the secular eccentricity of planetesimals is set only by the gravity of the companion and gas drag . We carried out this exercise while varying σ i , with the results presented by orange curves and points in Fig. 10d.
Turning off disk gravity eliminates the region of vanishing v c (see Fig. 1a) and collision velocities around 1.8 AU go up substantially (see Fig. 9a). Therefore, as one might expect, planetesimal growth requires significantly higher d min = 6.9 km, compared to 1.6 km when disk gravity is fully accounted for (assuming the fiducial σ i = 10 −4 ). We note that without disk gravity v coll in the outer disk is lower than with it, but this still does not help to keep d min low.
In the absence of disk gravity, reducing σ i has very little effect on d min because for σ i 10 −4 , destructive collisions arise mostly due to secularly excited eccentricities and are therefore unaffected by further reducing σ i . But at sufficiently high σ i , the collision velocities are dominated by random motions, and d min is no longer affected by the disk model. We do not consider such high values of σ i , but we do find at σ i = 10 −3 that d min is lower in the model with disk gravity than without. This reversal occurs because for σ i = 10 −3 , the outer part of the disk becomes the most favorable environment for planet formation in the absence of disk gravity because the random velocities are lower there. With disk gravity the collision velocities in the outer disk are still dominated by the secular eccentricities, so planet forma- tion must occur in the inner disk where the random velocities are higher.
Interestingly, d min increases very slightly as we go to very low values of σ i , going from 6.6 km at σ i = 3×10 −5 to 6.9 km at σ i = 10 −5 . Possibly that is because decreasing σ i decreases the timescale for collisional evolution (see Eq. (A.25)) while leaving the timescale for radial inspiral unaffected. As a result, erosion in collisions is more of an issue at low σ i because small bodies have less time to be removed from the system before they collide with larger ones. But quantitatively, this effect is rather small, and in fact not present for the 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 rate artificially increased or decreased compared to the value given by Eq. (6) by a factor χ ir as indicated in Table 1. Simulations 25-31 explore what happens in the fiducial model when the inspiral rate is varied with χ ir changing from 0 to 100. One possibility that we wanted to test through this exercise is that the removal of small planetesimals by inspiral should reduce the detrimental effect of erosion on planetesimal growth, with the expectation that high χ ir > 1 (i.e., faster removal of small fragments) would facilitate growth, while χ ir < 1 (weakened inspiral) would suppress it. However, Fig. 10e reveals a pattern of d min behavior that runs counter to this expectation. In fact, in the no-inspiral case (χ ir = 0) we find d min ≈ 1 km to be lower than in the case with full inspiral (χ ir = 1) for which d min ≈ 1.6 km. And χ ir > 1 results in even larger d min .
We believe that there are two reasons for this behavior. First, in the outer disk, inspiral does not help much with cleaning out small fragments, simply because they get primarily ground down below the smallest size captured in our simulations. Second, the radial inspiral negatively impacts the growth of large planetesimals by shifting their semimajor axis away from the dynamically quiet zone where their growth is most favorable. Indeed, Fig.  11a shows the results of a no-inspiral simulation with fiducial parameters (χ ir = 0, analogous to a simple addition of a number of multi-zone simulations not interacting with each other) for d init /d min = 2, which can be directly compared to Fig. 8c (for which χ ir = 1). One can see that without inspiral, the largest objects form closer to a dq = 1.8 AU where v c = 0, which certainly facilitates their growth compared to the full-inspiral calculation, in which the 300 km object emerges interior to 1.4 AU. But the general conclusion that we can draw from this exercise is that the gas drag-driven inspiral has a rather modest effect on the outcome of planetesimal growth.

Discussion
The results of the previous sections highlight the important differences between planet formation around single stars and in binaries. In the binary case, erosion plays a much more important role in planetesimal size evolution, which is caused by the obvious differences in underlying dynamics: whereas in disks around single stars planetesimal motions are excited only by gravitational perturbations between numerous planetesimals (and, possibly, turbulence in the disk), planetesimals in binaries get additionally (and very strongly) excited by the gravity of the companion star and the disk. This can easily raise the collision speeds of planetesimals of comparable size into the km s −1 regime (see Fig. 9), much higher than in disks around single stars. As a result, growth to 10 2 − 10 3 -km size objects in binaries is possible only if planetesimals start already rather large, with d init ∼ 1 − 10 km.
This places planets in binaries such as γ Cep in a unique category of dynamical extremophiles or hyperdynamophiles 5objects capable of surviving and growing despite the harsh dynamical environment in which they reside. The very existence of such planets puts planet formation theories to a very stringent test.
Our parameter space exploration in Sect. 6 is not entirely exhaustive: because of the multi-dimensional space of physical parameters we could vary at most one or two of them at once with respect to our fiducial model, meaning that some domains of the parameter space are left unexplored. Also, our models rely on certain simplifying assumptions that were necessary to reduce the number of model inputs, for example, the initial planetesimal size d init is the same everywhere in the disk, the amplitude of random motions σ i is the same regardless of location and object's size, etc. Nevertheless, our results do allow us to robustly single out the most important physical ingredients determining the outcome of planetesimal growth in binaries.
The gravitational effect of the protoplanetary disk is the main such ingredient and should not be overlooked. The inclusion of disk gravity significantly reduces d min in our fiducial case, by about a factor of four compared to the case in which disk gravity is ignored (see Fig. 10d). Eccentricity e 0 and orientation d of the protoplanetary disk are two other key parameters determining the outcome of planetesimal growth. Acting in conjunction with disk gravity, they can lead to the emergence of dynamically quiet locations in the disk, providing favorable conditions for coagulation.
Inclusion of radial inspiral of solids is what makes our models fully global. We find that inspiral somewhat complicates planet formation by causing the migration of growing planetesimals out of the dynamically quiet zone in the disk and raising d min compared to the case in which inspiral is ignored (see Fig.  10e). But the overall effect is not as large as one might have expected (see Sect. 6.7). Next we discuss some other consequences of our results.

Implications for planet formation in binaries
About a dozen binaries are currently known to harbor planets in S-type configurations (Chauvin et al. 2011;Marzari & Thebault 2019), some of them with binary semimajor axis even smaller than that of γ Cep (i.e., a b < 20 AU). Our results provide a pathway to understanding the origin of such planets. Our previous attempt to explain their existence ) focused on explaining the formation of planets in some of these systems (Chauvin et al. 2011) at their current locations. But in this work we are more interested in the overall possibility of planet formation at one to a few AU separation in a γ Cep-like system.
Our finding that d min only weakly depends on Σ 0 is quite important. As we vary Σ 0 from 500 g cm −2 to 20,000 g cm −2 (i.e., between 0.3 to 12 times the minimum mass solar nebula density at 1 AU (Hayashi 1981)), implying the variation of the protoplanetary disk mass from 1.8M J to 70M J (for the assumed Σ profile (1) with a 0 = 1 AU and a out = 5 AU), d min changes by less than a factor of two. What this means physically is that, even with the variation of Σ 0 , the dynamically quiet region can still exist somewhere in the disk, ensuring that planetesimal growth can be achieved in the disk globally (things may get more complicated at very low Σ 0 when a secular resonance appears in the disk; see Fig. 1). Figure 12 shows a dq as a function of e 0 for various values of Σ 0 in a system which has, other than the variation of e 0 and Σ 0 , our fiducial parameters. We have marked the fiducial model with an orange star along the Σ 0 = 1000 g cm −2 curve. This shows that over a wide range of disk masses, there is a dynamically quiet region in the disk. For each disk mass however, there is only a relatively narrow (approximately factor of two) range in e 0 over which the dynamically quiet region exists in the part of the disk outside 1 AU.
At low disk masses, a secular resonance appears in the disk (see the curve for Σ 0 = 500 g cm −2 in Fig. 1). This resonance dramatically increases collision velocities in the outer disk as il- Location a dq of the dynamically quiet zone in an apsidally aligned disk. We assume the fiducial disk plus binary parameters, except for Σ 0 and e 0 which are varied as shown on the plot. The fiducial model is marked by an orange star along the curve corresponding to Σ 0 = 1000 g cm −2 . lustrated in Fig. 9b. Collision velocities between d p > 10 km planetesimals are even higher than in that figure, since their orbits are less affected by gas drag.
In a static disk model such as ours, this resonance is localized, and planet formation can still proceed in other regions of the disk. In reality though, the disk evaporates over time, which changes the resonance location and causes the resonance to "sweep" through the disk. This effect is similar to the sweeping secular resonances involving planets (Ward 1981).
On the other hand, the dynamically quiet location also sweeps through the disk, meaning that over the disk lifetime, planetesimal growth may be promoted in many parts of the disk. Exactly what is the effect of this combined sweeping dynamically quiet location and sweeping resonance will depend on the details and timing of the disk dispersal.
It is worth noting that Fig. 12 shows that even for a quite low surface density disk with Σ 0 = 500 g cm −2 (corresponding to 6M ⊕ of solid material), the dynamically quiet region will not exist if e 0 is greater than about 0.02, corresponding to a massweighted mean disk eccentricity of 0.05. This is within the range found in numerical simulations, but on the low end of that range (see Sect. 7.3 for details).
The low levels of planetesimal random motions required to drive d min into the kilometer size range may seem problematic: σ i = 10 −4 featured in our fiducial model corresponds to random velocities of about 10 m s −1 . However, even this is above the level of the escape speed from the surface of our initial planetesimals in most cases, which is the most the two-body relaxation would pump up the random motions to. Moreover, given the high relative velocities (up to ∼ v c ) between the planetesimal populations of different sizes, both viscous stirring and dynamical friction are likely to be very inefficient (Stewart & Ida 2000;Rafikov 2003a), further reducing the amplitude of planetesimal random motions.
The precise location and time at which planetesimals reach our assumed threshold size of 300 km are important outcomes of our models. They turn out to be rather strong functions of d init /d min ≥ 1. If the size of seed planetesimals d init is close to Article number, page 17 of 29 A&A proofs: manuscript no. finalPaper d min , planetesimals grow rather slowly, causing their substantial gas drag-driven inspiral, so that d p reaches 300 km somewhat closer to the star than a dq (see Fig. 11b). As d init /d min increases, planetesimals tend to grow to large sizes closer to a dq (see Fig. 8c for d init /d min = 2). However, further increase in d init /d min causes the location of preferential planetesimal growth to shift inward again (see Fig. 8d), presumably because of faster coagulation there. Such non-monotonic variation with d init /d min , as well as the possibility of late-time planet migration (that must have certainly affected S-type systems harboring hot Jupiters), complicate the effort to constrain the initial planetesimal size d init based on the semimajor axes of observed planets within binaries.
At the same time, the time to reach 300 km is a monotonic function of d init /d min , steadily falling from ∼ 1 Myr to < 10 5 yr as this ratio increases from 1 to 5 in our fiducial model. The expectation of short disk lifetimes in S-type binaries (Harris et al. 2012;Barenfeld et al. 2019;Zurlo et al. 2020) may then favor larger d init , for which the growth time is shorter.
As previously stated, our fiducial disk model contains less than four Jupiter masses of material. The planet in the γ Cephei system was determined to have around ten Jupiter masses (Benedict et al. 2018). Unless the disk is being replenished by material infalling from the circumbinary disk, this would rule out all disk models with less than Σ 0 5000 g cm −2 in this particular system.

Implications for planetesimal origin
Our finding that some minimum size d min ∼ 1 − 10 km of seed planetesimals is necessary for ensuring the success of planetary genesis in S-type binary configurations has important consequences for theories of planetesimal formation. Two main ideas for the origin of planetesimals -steady growth by coagulation and rapid formation via streaming instability (Johansen et al. 2014)-predict very different formation pathways.
Growth by particle sticking must necessarily cover all sizes in the range from dust grains to multi-kilometer objects (subject to substantial gravitational focusing). This formation mode is known to be problematic because of the bouncing barrier around millimeter to centimeter sizes (Blum & Wurm 2008) and the rapid radial inspiral threatening meter-size objects (Weidenschilling 1977). In S-type binaries this planetesimal formation mechanism encounters another major bottleneck -efficient growth is impossible until seed objects are 1-10 km in size, substantially larger than any of the aforementioned barriers. This additional obstacle significantly reduces chances of slow planetesimal growth by coagulation in binaries.
On the contrary, the competing mechanism relying on the operation of streaming instability appears to naturally produce large planetesimals with sizes in the range from tens to hundreds of kilometers (Simon et al. 2016;Schäfer et al. 2017). This mechanism is also rather fast (with a typical timescale of tens of local orbital periods), resulting in a rapid production of the population of seed planetesimals with d init exceeding d min , which would enable further growth into the planetary regime as we have demonstrated here.
Thus, results of our work provide indirect support for the idea that streaming instability is a pathway for forming planetesimals in binaries (similar reasoning can be found in Thebault 2011). Of course, whether and how this instability would operate in an eccentric disk in presence of strong perturbations from the secondary, is an open question. But if it works in binaries, then there is no reason for it not to be the mechanism for forming planetesimals also in disks around single stars.

Implications for simulations of gaseous disks in binaries
Given the strong sensitivity of planetesimal growth to disk eccentricity e 0 and orientation d that we identify in this work, it is imperative to better constrain these variables. Unfortunately, observations are not yet at a stage where one can reliably determine eccentricity and orientation of protoplanetary disks in binaries with small separations a b 30 AU. Instead, one has to rely on numerical simulations for guidance in that matter.
A number of studies have explored the amplitude of disk eccentricity excited by the gravity of an eccentric companion in a coplanar S-type configuration using 2D hydrodynamic simulations. The conclusions that can be drawn from these works very strongly depend on the assumed physical inputs: equation of state, treatment of disk gravity, inner boundary conditions, and so on.
For example, studies using a locally isothermal equation of state (Kley & Nelson 2008;Paardekooper et al. 2008;Marzari et al. 2009;Zsom et al. 2011;Regály et al. 2011;Martin et al. 2020) typically find rather high (mass-weighted mean) disk eccentricities, e d ∼ 0.1 − 0.3, with e d profile typically exhibiting a double-peaked structure as a function of semimajor axis and with high values of e d close to the disk center. These results may need to be viewed with caution in light of the recently identified issues with using the locally isothermal equation of state for numerical modeling of the disk-perturber tidal interaction (Miranda & Rafikov 2019). On the contrary, studies that attempt to represent disk thermodynamics in a more realistic way (Marzari et al. 2012;Gyergyovits et al. 2014) find a dramatic reduction of disk eccentricity 6 , down to e d ∼ 0.02 − 0.05. As our results demonstrate, disk eccentricity at this level would allow planetesimal growth to proceed successfully starting from objects of a few to 10 km in size (see Fig. 1b). We note that in our model, e d = 2.5e 0 .
Similarly, accounting for the effect of disk self-gravity on its dynamics has also been shown to suppress e d ). Another interesting finding is that S-type disks tend to be less eccentric in more eccentric binaries (Marzari et al. , 2012Regály et al. 2011), presumably because of their more compact size.
Also important for the possibility of efficient planetesimal growth is the orientation of the disk, since the dynamically quiet zone, greatly facilitating planetesimal growth, arises only in disks which are close to apsidal alignment ( d ≈ 0) with the binary orbit (see Fig. 1c). Unfortunately, at the moment there is little clarity regarding this issue. Marzari et al. (2009) find the disk to be precessing, that is, there is not a well-defined value for d . In some other studies, the disk is aligned with the binary orbit (Müller & Kley 2012;Martin et al. 2020). Gyergyovits et al. (2014) find the disk to be close to alignment, with | d | ≈ 0.6 rad. Marzari et al. (2012) find anti-alignment ( d ≈ 180 • ), whereas Marzari et al. (2009) find either anti-alignment or precession, depending on whether self-gravity is included. Paardekooper et al. (2008) find either anti-alignment or precession depending on the flux limiter used in their code.
To summarize, the existing numerical studies of S-type disks in binaries have not yet arrived at a consensus regarding the values of either e 0 or d , and their dependence on various physical inputs, for example, realistic thermodynamics, disk self-gravity, viscosity, etc. Given the importance of these characteristics for understanding planet formation in S-type binaries, we encourage future numerical efforts to focus on measuring the true value of d as well as disk eccentricity e 0 , using the most realistic physical inputs, such as disk thermodynamics and self-gravity.

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 catastrophic disruption. For that reason  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 d min for a variety of models featuring different values of d c and v c . As discussed in , the collision rates and velocities are entirely determined by d c and v c (for the as-sumed level of random motions), regardless of what combination of disk parameters lead to the specific d c and v c (we set v K as appropriate at 2 AU in the disk around γ Cephei). Therefore, by exploring the range of plausible v c and d c we check the validity of these heuristics for all values of the disk parameters in our model. Figure 13 displays the results of this calculation (each panel corresponds to a different value of v c ), with the solid blue line showing d min as a function of d c , computed by our code using the procedure described in Sect. 5.2.1. The red line shows the minimum value of d init such that no catastrophic disruption occurs for planetesimals of size d init or larger. The green line shows the minimum value of d init such that no body larger than d init suffers an erosive collision with a partner with mass ratio greater than 1%. The same criteria were used to determine the gray and black regions, respectively, in Fig. 4 of .
This figure shows that the criterion of no catastrophic disruption is too optimistic, typically predicting d min between two and five times smaller (depending on d c and v c ) than the one determined by the coagulation-fragmentation simulation. On the contrary, the criterion of no erosion between bodies with a mass ratio greater than 1% is too restrictive, with the green curve exceeding the blue curve by a factor two to five.
Interestingly, the blue curve is often roughly equally separated from both red and green ones (in log space), which motivated us to also show the dashed blue line computed as the geometric mean of the predictions from the "no catastrophic disruption" and "no erosion" models. As one can see, it matches the simulation-based d min (blue curve) surprisingly well, within a factor of two or even better (although we do not claim this to have any deep meaning). Thus, the geometric mean of the predictions for d min based on the two criteria used in  can be employed as a rough approximation for predicting the success of planetesimal growth in the absence of sophisticated coagulation-fragmentation simulations.

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. 2006Beaugé et al. 2010). In particular, Thebault (2011) found that even d init 10 2 km planetesimals do not ensure planetesimal growth in HD 196885 system at the location of the known planet. All these studies included the effect of gas drag on planetesimal dynamics, and Paardekooper et al. (2008) and Beaugé et al. (2010) even considered the possibility of the disk being eccentric.
However, as noted already in Rafikov (2013b) and Rafikov & Silsbee (2015a,b), all these calculations ignored the effect of disk gravity, which dramatically changes planetesimal dynamics. For example, Beaugé et al. (2010) find low planetesimal collisional velocities in the outer disk, whereas we find a very hostile dynamical environment there (see Fig. 9). The existence of the dynamically quiet location in apsidally aligned disks substantially alleviates the fragmentation bottleneck, allowing in many cases planetesimal growth to proceed even starting at d init ∼ 1 km (see Fig. 10).
Another factor that distinguishes our study from others in terms of realism is our reliance on the global coagulationfragmentation calculations fully accounting for the gas dragdriven inspiral of planetesimals. Whereas Thébault et al. (2006Thébault et al. ( , 2008Thébault et al. ( , 2009) and   we are able to model the evolution of the planetesimal size distribution in its entirety, fully accounting for the fine details of planetesimal dynamics at every step.
Compared to , we find less need for a massive ( 30M J ) protoplanetary disk: the lowest value of Σ 0 allowing growth to proceed starting at d init ∼ 1 km in Fig.  10a corresponds to a disk mass of just 2M J . This difference arises primarily from our current full treatment of the planetesimal size evolution, but also from focusing on planet formation not at a given location but globally, in the whole disk.

Conclusions
We present the first fully global coagulation-fragmentation framework for exploring planetesimal growth in S-type binary systems, such as α Cen and γ Cep. Our framework fully accounts for the specifics of planetesimal dynamics in binaries, including perturbations not only due to the eccentric stellar companion, but also due to the gravity of the eccentric protoplanetary disk in which planetesimals orbit. We include the effects of gas drag on the eccentricity dynamics of planetesimals, as well as on their radial inspiral. Thus, this framework brings studies of planet formation in binaries to an entirely new level of realism.
By running a suite of simulations with varied model inputs, we are able to delineate the conditions under which planetesimals can grow to large sizes (hundreds of kilometers) despite the adverse effects of fragmentation. The difficulty of planet formation for a given disk model was measured via the parameter d min -the smallest size of seed planetesimals that allows their growth into the planetary regime -with the smaller d min implying more favorable conditions for planet formation. Based on the results of these calculations, we are able to draw the following general conclusions.
-For most disk parameters considered in this paper, planet formation in binaries such as γ Cephei can successfully occur provided that the initial planetesimal size is 10 km; however, for favorable disk parameters, this minimum initial size can go down to 1 km. -The gravitational effect of the protoplanetary disk plays the key role in lowering the minimum initial planetesimal size necessary for sustained growth by a factor of four. This reduction can be achieved in protoplanetary disks apsidally aligned with the binary, in which a dynamically quiet zone appears within the disk provided that the mass-weighted mean disk eccentricity 0.05 (corresponding, in the context of our model, to e 0 0.02). -We identify disk eccentricity (e 0 ) and the apsidal misalignment angle ( d ) as the parameters to which planetesimal growth is most sensitive. Whenever the dynamically quiet location exists in the disk, growth is also strongly dependent on the level of random motion of planetesimals. -For the models explored in this work, the outcome of planetesimal evolution is only weakly dependent on the disk surface density Σ 0 . -Radial inspiral of solid material appears to make little difference to the critical planetesimal size from which collisional agglomeration can occur. -The planetesimal strength, as parameterized by the two models of Stewart & Leinhardt (2009), makes only about a factor of two difference to the critical planetesimal size. -Planetesimal growth to large sizes can occur even in the face of some erosion, but in the environments with frequent collisions leading to catastrophic disruption, growth will not occur. We used our 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 , 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.
In a given annulus, our code solves the standard coagulationfragmentation problem, which effectively is a version of the coagulation equation of Smoluchowski (1916) modified to account for fragmentation (see Eq. [1] of Rafikov et al. 2020). Our framework employs mass space discretization in the form of logarithmically spaced mass bins with constant M i+1 /M i . The grid extends from the minimum mass corresponding to a 10 m object to a maximum mass corresponding to a 300 km object. Fragments with sizes below 10 m that form in planetesimal collisions are simply removed from the system (i.e., total mass is not conserved for this reason alone). This leads to certain numerical artifacts discussed in Sect. 5.1.
At each time step, the code calculates the mean collision rate R i j between planetesimals in bin i and bin j for all i and j in [1, N bins ]. This calculation is described in Appendix A.4. The number of collisions Z i j in a given time step of length ∆t between bodies in mass bins i and j is then drawn randomly from a Poisson distribution with the mean R i j ∆t: We define an operatorF such that F i jk is the number of fragments in mass 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 This takes into account both addition of particles to bin i from collisional fragments (the F i jk operator), as well as collisional destruction of particles in bin i (the −δ i j and −δ ik terms). The factor of (1 + δ jk )/2 is to prevent over-counting when summing over both j and k. At times, due to a fluctuation, a number of particles in a particular mass bin may become negative. In this case we reset the number of particles in such a bin to zero; this effectively leads to addition of mass to the system. However, we checked several runs of our fiducial simulation, and found that relative to the total mass in the simulation, less than 1 part in 10 5 was added in this way.

Appendix A.2: Collisional outcomes
This section describes the physics which goes into the determination of F i jk . The outcome of a collision between a body in mass bin i and one in mass bin j is dependent on the collision velocity. High-velocity collisions will result in many small fragments, and no large remnant bodies. Low-velocity collisions will result in one body which contains nearly all of the pre-collision mass and only very small fragments. Calculation of collision velocity v coll is described in Sect. A.3. Given a value of v coll , we use the recipe provided in Stewart & Leinhardt (2009) to determine the mass of the largest remnant M lr . We take the mass of the largest remnant to be where M tot = m i + m j is the total mass of the two incoming bodies, Q R is the specific center-of-mass energy of the collision given by and Q * RD is the critical value of Q R which leads to catastrophic disruption (which corresponds to M lr < 0.5M tot ); Q * RD is a function of the sizes and composition of the planetesimals, and v coll , which is given by Eq. 2 from Stewart & Leinhardt (2009). As our default, we used the parameters appropriate for their "strong rock" planetesimals; however we also ran simulations using their "rubble pile" planetesimals.
When M lr is calculated, it will in general be between M i and M i+1 for some i. We then add a fraction P i of the mass to bin i and P i+1 to bin i + 1 given by 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 Here b is an adjustable parameter which determines the cutoff mass for catastrophic collisions: whenever Eq. (A.4) yields M lr < 2bM tot , we put all of the mass into fragments (as described below) without leaving a single largest remnant. This situation often arises in very energetic collisions when M lr is small. The fragment mass is assigned to mass bins with 0 < i ≤ i cutoff , where i cutoff is the index of the largest bin with mass less than M cutoff . This is done such that where Υ jk is a normalization constant. The exponent of M i is 1 + ξ rather than ξ because of the logarithmic spacing of the mass bins. In order to solve for Υ jk , we note that where nonpositive values of i correspond to (hypothetical) mass bins with M i less than the minimum size considered in our simulations. Letting R ≡ M i+1 /M i , we can then solve for Υ jk and write Some mass is lost from the simulation because it is contained in bodies with mass less than M 1 . This lost mass is equal to Because of the self-similar form of the fragment mass spectrum, we were able to implement this fragmentation algorithm with numerical cost scaling as O(N 2 bins ) using explicit time stepping (see Rafikov et al. 2020 for details). This considerably reduces the computational burden of our calculations.

Appendix A.3: Collision velocity
We calculate the collision velocity v coll by considering two separate dynamical regimes. In the first, the relative motion of planetesimals is dominated by their secularly excited eccentricities. Equation (64) of RS15a then gives the relative eccentricity e i j between planetesimals as a function of their masses. Given e i j , there is a distribution of relative encounter velocities v i j between v min = 0.5e i j v K and v max = e i j v K , where v K 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 In the second dynamical regime, we assume that random motions dominate over the secularly excited eccentricities. We assume that planetesimals have eccentricity vectors (k, h) = (k s + k rand , h s + h rand ) where (k s , h s ) and (k rand , h rand ) are the secular and random components, respectively; k rand and h rand are assumed to be drawn from independent Gaussian distributions ρ(k rand ) = 1 where we take σ e = 2σ i . This leads to a distribution of the relative random eccentricity e rand = (k rand,1 − k rand,2 ) 2 + (h rand,1 − h rand,2 ) 2 given by the Rayleigh distribution In principle, for each set of collision partners, we should then take the collision velocity to be given by where e rand is drawn from the distribution in Eq. (A.15), and λ is drawn from the distribution in Eq. (A.12).
In practice, doing this at each time-step would be too timeconsuming. Instead, we use the following algorithm. Let F(v) be the true collision velocity distribution for a given set of collision partners in a given annulus. We approximately calculate N 1 points in v-space, such that there is a constant fraction of the distribution F(v) between consecutive points: this means that the i th point v i satisfies the relation v i 0 F(v)dv = (i + 0.5)/N 1 , for i ∈ [0, N 1 − 1]. We then calculate the collision rates and outcomes for each set of collision partners at the beginning of the simulation, using the N 1 different velocities. During the simulation, in each timestep, and for each set of collision partners, we randomly choose one of the N 1 collisional outcomes and rates to use for that timestep.
To calculate the N 1 values v i , we initially create two vectors z rand and z λ of length N 2 with a constant fraction of the distribution for e rand between each entry of z rand and a constant fraction of the distribution of λ between each entry of z λ . This means that z i+1 rand z i rand ρ(e rand ) = 1/N 2 , and a similar relation exists between z λ and the distribution for λ implied by Eq. (A.12).
We then draw N 3 samples of the collision velocity using Eq. (A.16) and randomly selected entries from z rand and z λ . We then select the N 1 values of v i so there are an equal number of the N 3 sample velocities between consecutive values of v i . We take N 1 = 11, N 2 = 100, and N 3 = 990.
We ignore shear motion in calculating the collision velocity because if the shear motion is the dominant effect, then the planetesimals will completely merge in collisions, so there is no reason to consider shear motion except in how it affects the collision rate.
We note that the differential radial drift described in Sect. 3.5 is not included in our calculation of v coll as it does not typically contribute substantially to the collision velocities between bodies of the sizes that we consider in this paper.

Appendix A.4: Collision rate
Next we describe our calculation of the collision rate. In the first dynamical regime (where the relative motion of planetesimals is dominated by their secularly excited eccentricities), we use Eq. (63) from RS15a (corrected so that e max is replaced with v max /a p ). This gives a flux F i j of planetesimals from bin j past a planetesimal in bin i of F i j = E( √ 3/2)Σ j v max π 2 a p m j I i j , (A.17) where I i j 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: The distribution of I ij = (q i − q j ) 2 + (p i − p j ) 2 is therefore given [analogous to Eq. (A.15)] by Then, averaged over I i j , the average value of F i j is given by ρ(I i j )F i j dI i j = E( √ 3/2)Σ j v max 2π 3/2 σ i m j a p . (A.20) The collision rate for a given particle in bin i with the population in bin j is given by R i j = F i j σ i j , where σ i j is the collision cross section including gravitational focusing: To address these concerns, we instead pick the largest time step ∆t coag such that for all i n i (t) − n i (t + ∆t coag) n i (t) < 1 or m i (n i (t) − n i (t + ∆t coag )) where the angle brackets refer to averaging over the distribution of encounter velocity in Eq. (A.12). To save computational time, we approximate the average over the encounter velocity distribution by evaluating the collisional outcomes using the median encounter velocity. Here, as before, a time step is acceptable if for each bin, either the fractional decrease in mass is less than 1 , or the decrease in mass is less than 2 times the total mass in the system. Time step ∆t coag is recalculated every time we update the mass distribution, separately for each annulus, and the minimum value among all the annuli is used for each time step. It is not obvious from first principles which values of 1 and 2 are appropriate. But rows 10 and 11 of Table 2 (which explores sensitivity of our calculations to the numerical parameters used; see Appendix C) show that variation of these values by a factor of two (or more) from our fiducial values of 0.05 and 10 −6 have a negligible effect on d min .
In the simulation used to create Fig. 4, we used a lower value of 2 so that the distributions would remain smooth, even in regions of mass-space containing a negligible amount of mass. The success of the coagulation process, and the shape of the distribution at sizes greater than d init in those simulations is not affected by this change.
Appendix A.6.2: Radial inspiral time step The code uses a different time step for the radial drift. Mass transport due to radial drift is recalculated using Eq. (A.29) every interval of time ∆t insp such that (A.32) which guarantees that the inspiraling particles do not cross the full radial bin width in time ∆t insp (v drift is the inspiral speed). While we see from Eq. (A.29) that drift must not be larger than 1, it is apparent from the results of Sect. B.3 below that there is little advantage to picking drift less than 1. Therefore, we use drift = 1.

Appendix B: Tests of the code
In this section we describe the tests we did of the code against known analytic solutions in certain limiting cases. We independently verify three code functionalities -coagulation, fragmentation, and radial inspiral -and find that the code does a reasonable job of reproducing known solutions to the accuracy warranted given the imperfect knowledge of the physical processes involved.
Appendix 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 obtained for kernels A i j = 1 (Smoluchowski 1916), A i j = i + j and A i j = i j (Trubnikov 1971), with the initial condition of N particles of a single mass bin. Our code successfully reproduces these solutions with small deviations as discussed below. The top panel of Fig. B.1 shows the number of particles per unit mass as a function of mass for three different times for the kernel A i j = 1. In this figure we have used a mass bin spacing at low masses. Other than this discrepancy at low mass, there is nearly perfect agreement between our code and the analytic solution.
The middle panel of Fig. B.1 shows the number of particles per unit mass as a function of mass for three different values of η for the kernel A i j = i + j. The distribution evolves more rapidly than in the case of the constant kernel because of the higher collision rates. In this case, in addition to the scatter at low masses, we see that the numerical solution is slightly behind the analytic solution at the high-mass end of the distribution.
The bottom panel of Fig. B.1 shows the number of particles per unit mass as a function of mass for three different values of η for the kernel A i j = i j. 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. effect on the calculation outcomes. Table C.1 shows the results of this exercise, in which we use d min and d rubble min as our sensitivity metrics.
Simulation 1 uses the fiducial set of numerical parameters that we employ in all our production runs. As in Table 1, in all other runs the parameter that is different from Simulation 1 is highlighted in bold in each row. We see that d min does not depend strongly on the numerical parameters, which justifies our particular choice of their values.
Article number, page 28 of 29 Kedron Silsbee & Roman R. Rafikov: Global simulations of planetesimal growth in stellar binaries