Issue 
A&A
Volume 590, June 2016



Article Number  A62  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201628349  
Published online  11 May 2016 
Modelling circumbinary protoplanetary disks
II. Gas disk feedback on planetesimal dynamical and collisional evolution in the circumbinary systems Kepler16 and 34
^{1}
School of Physics, University of Bristol, H. H. Wills Physics
Laboratory,
Tyndall Avenue,
Bristol
BS8 1TL,
UK
email:
stefan.lines@bristol.ac.uk
^{2}
CNRS, IRAP, 14
avenue Édouard Belin, 31400
Toulouse,
France
^{3}
Université de Toulouse, UPSOMP, IRAP, 31400
Toulouse,
France
^{4} Astronomy Unit, School of Physics & Astronomy, Queen
Mary University of London, UK
^{5}
DAMTP, University of Cambridge, Wilberforce Road, Cambridge
CB3 0WA,
UK
Received: 19 February 2016
Accepted: 1 April 2016
Aims. We investigate the feasibility of planetesimal growth in circumbinary protoplanetary disks around the observed systems Kepler16 and Kepler34 under the gravitational influence of a precessing eccentric gas disk.
Methods. We embed the results of our previous hydrodynamical simulations of protoplanetary disks around binaries into an Nbody code to perform 3D, highresolution, interparticle gravityenabled simulations of planetesimal growth and dynamics that include the gravitational force imparted by the gas.
Results. Including the full, precessing asymmetric gas disk generates high eccentricity orbits for planetesimals orbiting at the edge of the circumbinary cavity, where the gas surface density and eccentricity have their largest values. The gas disk is able to efficiently align planetesimal pericenters in some regions leading to phased, noninteracting orbits. Outside of these areas eccentric planetesimal orbits become misaligned and overlap leading to crossing orbits and high relative velocities during planetesimal collisions. This can lead to an increase in the number of erosive collisions that far outweighs the number of collisions that result in growth. Gravitational focusing from the static axisymmetric gas disk is weak and does not significantly alter collision outcomes from the gas free case.
Conclusions. Due to asymmetries in the gas disk, planetesimals are strongly perturbed onto highly eccentric orbits. Where planetesimals orbits are not well aligned, orbit crossings lead to an increase in the number of erosive collisions. This makes it difficult for sustained planetesimal accretion to occur at the location of Kepler16b and Kepler34b and we therefore rule out in situ growth. This adds further support to our initial suggestions that most circumbinary planets should form further out in the disk and migrate inwards.
Key words: methods: numerical / hydrodynamics / planets and satellites: formation / protoplanetary disks / binaries: close
© ESO, 2016
1. Introduction
Circumbinary planets that orbit close to their parent stars form an interesting subset of the extreme planetary systems discovered to date. The strong gravitational perturbations that the stars exert on the protoplanetary disk can significantly alter the dynamics of both planetesimals (Lines et al. 2014; Paardekooper et al. 2012; Meschiari 2012a; Thébault et al. 2006) and gas (Lines et al. 2015; Paardekooper et al. 2008; Marzari et al. 2013; Pelupessy & Portegies Zwart 2013) that orbit in close proximity to the system barycenter. Studies of the dynamics and collisional evolution of planetesimals in these hostile disks show that their eccentricities are pumpedup by these perturbations and their orbits misaligned, causing high velocity collisions during orbit crossing events (Lines et al. 2014; Paardekooper et al. 2012; Moriwaki & Nakagawa 2004). Lines et al. (2014) found that such a large number of high velocity encounters produces an overwhelming proportion of erosive collisions and showed that sustainable planetesimal growth was not possible at the orbital radius of the majority of discovered circumbinary planets.
The situation is more complex than these purely Nbody studies reveal however, since planetesimals will feel both aerodynamic drag from, and the gravitational potential of, the gaseous component of the protoplanetary disk. A number of studies on the dynamics of the circumbinary gas disks show that the fluid is also perturbed by the binary, leading to the tidal truncation of the inner disk and the generation of spiral density waves by both direct forcing from the binary on the disk and mode coupling between the binary and disk potentials (Pierens & Nelson 2013; Lubow 1991). This results in an eccentric precessing disk (Pierens & Nelson 2013). Lines et al. (2015, hereafter Paper I), explored a variety of fluid parameters for two known circumbinary planet systems, Kepler16 and Kepler34, and found that in all cases the inner disk became largely asymmetric with a buildup of gas at the disk apoapsis. Such regions of high surface density could present strong timedependent gravitational forcing and drag on the planetesimals.
Due to the nature of performing simulations that unify both Nbody and hydrodynamical effects, there is currently limited research on how the presence of a gas disk affects the dynamics of planetesimals and ultimately the collisions they undergo. Rafikov (2013) investigated how the gravitational potential from an axisymmetric disk affects embedded planetesimals, but we now know the disk is not in a static axisymmetric configuration. Additionally Marzari et al. (2013) examine both the gravitational and drag effects from a nonaxisymmetric gas disk by running hybrid simulations but, they do not consider the interaction between the planetesimals themselves. Their work shows that planetesimal eccentricities are elevated by interaction with the gas disk, in some areas of the disk by up to ten times the value seen in simulations performed without gas disk gravity enabled.
In this paper we embed the hydrodynamic results of our Paper I into an Nbody code to perform 3D, highresolution, interparticle gravityenabled simulations of planetesimal dynamics that include the gravitational forces imparted by the gas disk. In Sect. 2 we discuss our numerical method and initial conditions, in Sect. 3 we present the results of our simulations and in Sect. 4 we discuss our results in the context of previous work.
2. Numerical methods
Our Nbody simulations are conducted using the widely adopted code PKDGRAV (Stadel 2001; Richardson et al. 2000). PKDGRAV is a highly parallelised, multidisciplinary code capable of handling large N. An efficient tree code and multipole expansion allows for Nlog N scaling and hence the ability to run selfgravitating simulations with N≥10^{5} in a practical timeframe. The analytical collision model EDACM (Leinhardt & Stewart 2012) which has been previously used in simulations of planet formation (Leinhardt et al. 2015; Carter et al. 2015) is used to determine the outcome regime of a collision. EDACM takes the collision velocity, impact parameter, collider mass ratio and material parameters to assign the outcome as either perfect merging, partial accretion, nonerosive/erosive hit and run, erosive disruption and supercatastrophic erosion.
2.1. Gas potential
To avoid the computational constraints of using a fully hybridised code, we take a semianalytical approach to integrating the gas disk into the Nbody code. The surface density data from 2D hydrodynamical simulations of circumbinary gas disks around Kepler16 and Kepler34 performed in Paper I are used.
Fig. 1
Synthetic surface density maps of Kepler16 and Kepler34 created using the method described in Sect. 2. The orientation of the disk is chosen to match its phase at the end of the Nbody simulations in Fig. 4. 
The quasisteadystate (QSS) gas surface density data from Lines 2015 (simulations runs B and F) are used. The simulations correspond to nonselfgravitating disks since, for the disk masses explored here, the difference in gas surface density between self and nonselfgravitating disks is marginal. These runs correspond to a locally isothermal disk with uniform aspect ratio, h = 0.05, αviscosity, α = 10^{3}, initial surface density, (approximately half minimum mass solar nebula), a rigid inner boundary and no disk selfgravity. These fluid parameters are typical of circumbinary disk simulations. Selfgravitating gas disk simulations are ignored, since Lines et al. (2015) found that the gas disk only gravitationally interacts with itself for Σ_{disk} ≥ 2.5 × 10^{3} M_{⊙}/au^{2}. These hydrodynamical simulations are performed using FARGOADSG (Baruteau & Masset 2008; Baruteau 2008) over a polar mesh grid with N_{r} = 395 and N_{s} = 512 where N_{r} and N_{s} are the number of radial and azimuthal cells respectively. In Kepler16 the gas disk has precession period of P_{d} = 2000 P_{AB} and in Kepler34 has P_{d} = 3000 P_{AB} yielding P_{d} ≈ 250 yr for both cases. The surface density maps of these disks can be seen in Fig. 1 and are incorporated in the following way:
Fig. 2
Map of planetesimal accelerations from an Nbody test simulation calculated from the gas disk potential of Kepler34. 50% of planetesimals in the simulation are shown and are plotted as unit vectors with the colour corresponding to their acceleration magnitude. 

1.
The disk surface density is averaged over one binary orbit to remove transient features such as short period modes. This is a valid approximation as a) the gas surface density evolves on a much longer timescale than a binary orbital period and b) short period modes are most active at the gas disk inner edge which is 0.3 au interior to the planetesimal inner edge and therefore would not have a significant impact of planetesimal dynamics.

2.
The gravitational potential for each cell Φ_{c} is calculated from the averaged cell surface density Σ_{c} using direct summation. We use the Plummer potential Φ_{c} ∝ (r^{2} + ϵ^{2})^{−1/2}, with a smoothing length ϵ = 1.2H, where H is the disk thickness, to smooth out singularities and to account for the vertical extent of the 2D disk (Müller et al. 2012).

The gravitational potential is then Fourier transformed in angle into the lowest ten azimuthal modes (enough to generate a synthetic potential map that accurately describes the real output): (1)where φ is the cell position angle and Φ_{m} is the complex form of the potential contribution, at a given radius, from each azimuthal mode m. At each radius the complex potentials are recorded in a table to be read into PKDGRAV.
Fig. 3 Time evolution of planetesimal disk locations (top row) and eccentricity (bottom row) to quasisteadystate around Kepler16, under the influence of a full asymmetric gas potential.

4.
PKDGRAV applies an acceleration on the planetesimals from the reconstructed gravitational potential of the gas disk. The gas disk potential must first be rebuilt by combining the complex modes. Mode combination is done as follows: (2)where and are real (magnitude) and imaginary (phase) components of Φ_{m} respectively, t is the simulation time and ω is the precession frequency. This method retrieves the disk gravitational potential for any value of r, φ and t.

5.
The gas potential and acceleration are evaluated at the particles location using a bilinear interpolation.

6.
This additional acceleration is added onto the existing acceleration of the planetesimal calculated in PKDGRAV’s other functions which evaluate the force on the planetesimals from interplanetesimal and stellarplanetesimal interactions. The result of this implementation is shown in Fig. 2 where the planetesimal accelerations in response to the potential of a Kepler34 gas disk are shown.
For simulations where we test the static, symmetric circular disk only (m = 0), the same procedure applies but values of the azimuthal mode number greater than zero are ignored. The suite of simulations defined in Table 1, explores the effects of three different gas profile types and varying the size distribution of the planetesimal disk.
Parameter setup of each simulation A through L.
Fig. 4
Planetesimal disks at quasisteadystate (end of simulation at t ≈ 40 000 P_{AB}. AF: Location and eccentricity of planetesimals in disk around Kepler16. GL: Location and eccentricity of planetesimals in disk around Kepler34. 
Fig. 5
Change in the distribution of planetesimal mass between initial conditions and simulation end for the asymmetric potential only cases. Planetesimals are placed into bins of width 0.1 au and the increase or decrease of mass in each bin is plotted as a percentage change. 
2.2. Initial conditions
We test the response of a circumbinary planetesimal disk around two different stellar binary systems, Kepler16 (Doyle et al. 2011) and Kepler34 (Welsh et al. 2012) which have a total mass of 0.89 M_{⊙} and 2.07 M_{⊙} respectively. Kepler34 has M_{A} = 1.05 M_{⊙} and M_{B} = 1.02 M_{⊙} where M_{A} and M_{B} are the primary and secondary stellar mass, stellar separation is a_{b} = 0.22 au and binary eccentricity is e_{b} = 0.53. Kepler16 has M_{A} = 0.69 M_{⊙} and M_{B} = 0.20 M_{⊙}, a_{b} = 0.22 au and e_{b} = 0.16. Each planetesimal disk has a total mass of 2.8 M_{⊕} spread over N = 10^{5} planetesimals and the disk domain ranges from R_{in} = 0.6 to R_{out} = 3.0 au. The planetesimal density is 2.0 g/cm^{3} to match that of siliconrich rocky bodies. Two different models are used for the planetesimal disk size distribution: unimodal and mixed. For the unimodal size distribution, m_{p} = 8.3 × 10^{11} M_{⊙}. For the mixed size distribution planetesimals can occupy one of ten mass bins that are equal in their total mass and range from m_{p} = 3.8 × 10^{11}−9.7 × 10^{9} M_{⊙}. The large size of these planetesimals is a necessary condition; simulating much smaller bodies would require either the impractical computational demand of increasing the Nbody resolution or reducing M_{disk} for constant N and thus reducing the collision rate. We want to ensure realistic collision timescales so that any net accretion rate can be compared to the overall lifetime of the protoplanetary disk. Additionally our planetesimal masses are an order of magnitude larger than in Lines et al. (2014), as a decrease in the Nbody resolution for this work is necessary to evolve the disk to a quasisteadystate in a practical timeframe. Such large planetesimals, which may be able to form at this size through fast clumping assisted by streaminginstabilities (Johansen et al. 2007; Carrera et al. 2015), provide a convenient bestcase scenario as they have a large gravitational binding energy that makes them difficult to disrupt. Should it be found that even they are not able to undergo accretion, smaller bodies would only disrupt more easily. The planetesimal surface density follows Σ(r) ∝ r^{1.5} which is consistent with our previous Nbody study of circumbinary planetesimal disks (Lines et al. 2014).
A small value for the timestep is required (0.0025 yr) to accurately resolve the binary, which is modelled as two interacting Nbody particles, and maintain stability over thousands of binary orbits. As per Lines et al. (2014) we begin each simulation with an unperturbed planetesimal disk with the inclinations and eccentricities assigned from a Rayleigh distribution with dispersions of ⟨e⟩^{2} = 2 ⟨i⟩^{2} = 0.007 (Ida & Makino 1992). We enable collisions, and gas gravity if applicable, from t = 0 and allow the disk to reach a quasisteadystate. Only collisions that occur from after a quasi steadystate are used in the analysis and results. Each simulation runs for 3800 years in total (34 000 P_{AB} and 50 000 P_{AB} for Kepler16 and Kepler34 respectively), a small snapshot of the million year timescale over which planet formation occurs.
3. Results
3.1. Gas potential
3.1.1. Full asymmetric disk potential
We start by looking at how the presence of the gas disk potential modifies the dynamics of the planetesimals. The strongest gravitational influence comes from areas of highest surface density, most significant around the gas disk inner edge apoapsis. At apoapsis of an eccentric orbit the gas has its lowest velocity which leads to an inevitable overdensity of fluid. The small region over which this over density exists means that it acts almost like a massive body which stirs up planetesimal eccentricities as they traverse though it. Planetesimals that orbit exterior to the density peak are attracted towards it and move inwards while those that orbit interior to it move outwards. This can be seen in Fig. 2 by the strong acceleration directed towards the dense fluid peak (the peak itself contains planetesimals with very low accelerations due to being at the point of highest potential). This ultimately leads to the formation of a narrow annulus of planetesimals that can cause a separate interaction with surrounding planetesimals due to their selfgravitating nature.
Figure 3 shows the formation and evolution of this feature along with the associated planetesimal eccentricities, for an asymmetric gas potential on a Kepler16 disk. In panel A the narrow annulus or arm takes on a spiral form due to the planetesimals experiencing a shift from positive to negative torque as they transit across the density peak. Since ω ≪ ω_{p} where ω_{p} is the planetesimal orbital frequency, planetesimals make two orders of magnitude more orbits in the time it takes for the eccentric pattern in the gas disk to make one full precession. The planetesimals are thus concentrated quickly with the transient spiral disappearing after only a few thousand binary orbits and before the gas disk has rotated twice.
In panel D at 9000 P_{AB} the planetesimal disk has almost reached a quasi steady state with an eccentric planetesimal ring that passes through the density peak and more noticeably an asymmetric cavity interior to the ring almost devoid of planetesimals due to their migration outwards. Planetesimals interior to around 0.9 au do not strongly feel the gas gravity and lead relatively unperturbed orbits, with their eccentricities primarily set by the binary forcing (for a full description of planetesimal dynamics around binaries, see Lines et al. 2014). The presence of the gas potential leads to eccentricity waves that are launched outside of the density peak region and move inwards, bunching up at the peak at 1.2 au. These eccentricity waves eventually compress together to form a smooth distribution of planetesimal eccentricities which can be seen at the end of the simulation in panel F of Fig. 4 at 40 000 P_{AB}. The planetesimal disk under the influence of the full asymmetric gas potential has a much higher mean disk eccentricity than the gas free case. Around the density peak, eccentricities reach a maximum of e_{max} = 0.4. This is four times that found in the gas free simulation (e_{max} = 0.1), which occurs at the innermost planetesimal disk edge and is set by the binary forcing. At the end of the simulation as can be seen in panel F of Fig. 4 planetesimal eccentricities are enhanced from 1.0 to 1.7 au.
Fig. 6
Longitude of periastron for each planetesimal (black dots) in the Kepler16 system. 
Figure 5 assists in understanding the change in planetesimal mass distribution under the influence of the gas disk gravity. In Kepler16 the cavity is shown by the negative mass change between 1.1 and 1.5 au. The high density planetesimal ring appears as a small increase in planetesimal mass just exterior to this deleted region. In Kepler34 there is a more defined planetesimal ring shown by the large mass increase at 2.6 au, corresponding to the fluid density peak. The cavity stretches from 1.6 to 2.5 au but is not fully depleted in the middle since planetesimals in this region has a low acceleration from the gas potential due to their distance from any significant fluid mass.
Alignment of planetesimal pericenters (ϖ; see Fig. 6) is another effect contributed by the presence of gas disk gravity. The gas free case shows no alignment at all resulting in overlapping eccentric orbits throughout the disk, but the eccentricity of the planetesimals themselves is low. The orientation of the orbits becomes much more important when the eccentricities are raised to those values seen in the presence of the full asymmetric disk gravity. The asymmetric cavity is represented by the lack of planetesimals on one side of the disk, with low occupancy between ϖ ≈ 2 and ϖ ≈ 5 in the region of a = 1.0−1.2 au. The density peak that focuses planetesimals onto highly eccentric orbits also neatly aligns them with Δϖ ≈ 0.6 between 1.4 au and 1.6 au. The gas fails to align planetesimals between 1.2 au and 1.4 au however and leads to a wide range of orbit orientations, Δϖ ≈ 3.
Planetesimals in the region around and interior to the density peak experience a departure from their Keplerian velocity. It can be seen in Fig. 8 that the velocity distribution becomes winged for both Kepler16 and Kepler34. Near the peak some planetesimals obtain a much lower orbital velocity − these are the planetesimals captured onto an eccentric ring focused by the density peak, with low orbital velocities calculated at apoapsis. When these planetesimals transit through periapsis their velocities are increased, shown by the positive velocity wing.
Fig. 7
Change in planetesimal eccentricities by reducing planetesimal selfgravity. Snapshot at 14 000 P_{AB} for Kepler16 with an axisymmetric gas disk. 
A combination of eccentricity waves generated by the asymmetric potential and inflated orbit velocities leads to significantly higher impact velocities between planetesimals during collisions. In Fig. 9 the impact velocities can be seen to exceed ten times the typical value seen for the gas free case. The impact velocity distributions, shown in Fig. 9 are enhanced around the position of the orbital velocity wings since this is the location where planetesimals can obtain either the Keplerian value or the inflated or deflated value caused by the potential of the disk, and hence the relative velocities during collisions can be large.
Spatial collision maps (Fig. 10) can be used to show that the asymmetric cavities chiseled out by the strong gravitational pull of the high density ring exterior to it, are as devoid of collisions as they are of planetesimals. Collision rates drop off with distance from the binary even without gas, due to the drop in planetesimal number density, so the addition of this nocollision zone increases the difficulty for planetesimals to grow.
The most damning data in terms of the ability for planetesimals to growth is the collision type occurrence. A collision is put into one of 24 bins with a width of 0.1 au, and at the end of the simulation the percentage of each collision type can be displayed as a function of orbital radius. This can be seen for Kepler16 in Fig. 11 and for Kepler34 in Fig. 12. For Kepler16 the presence of the full asymmetric gas gravity causes a significant enhancement in the number of erosive collisions from the disk inner edge at 0.6 au to 1.8 au. Considering the mixed size distribution at the location of the planet Kepler16b at 0.7 au, 40% of all collisions are disruptive to one of the colliders, with this increasing to 70% between 1.0 au and 1.3 au. For Kepler34b at 1.1 au, 30% of collisions are disruptive with this increasing to 75% between 1.4 au and 1.6 au. Only exterior to 2.1 au does the number of growth enabling collisions exceed that of destructive ones.
3.1.2. Axisymmetric gas potential
For the axisymmetric case, the disk is static in time with the gas surface density smeared out azimuthally. This means that the density peak in the full asymmetric description is averaged out over an annulus leading to a lower density ring. However the focusing effect of the ring on the planetesimals is comparatively low with respect to the asymmetric case since the high density has been averaged out. For Kepler16, despite the dense gas ring having e = 0.0, the gravitational forcing from the binary causes the planetesimals to form a static double ring with e> 0.0 which is clearly seen in panel B of Fig. 4. The eccentricity distribution in panel E shows that e_{max} = 0.18 is much smaller than that found in the full asymmetric case.
An interesting observation is that in Kepler16, planetesimal eccentricities obtain much higher values than the gas free case, and particularly for the region beyond where the gas potential is strong (>2 au). This is contradictory to the results of Rafikov (2013) who find that the gravitational contribution from an axisymmetric disk helps to reduce excited planetesimals. To try and understand this further we run an additional simulation, identical to that of Run B (axisymmetric, Kepler16) but reducing the mass of planetesimal by a factor of 10^{5} (decreasing a planetesimals radius by a factor of ≈50). This reduces the gravitational interaction between them, essentially leaving us with the motion of test particles about the binary and subject to the potential of the gas. The results, shown in Fig. 7, clearly show that the selfgravitation of our planetesimals has a large contribution to the dynamics. Not only does reducing interparticle gravity remove the high eccentricities at the disk outer edge, but it also removes the peak associated with the high density gas disk ring.
Fig. 8
Planetesimal orbital velocity in Kepler16 (top row) and Kepler34 (bottom row) for both axisymmetric and asymmetric gas disk gravity inclusion. The Keplerian velocity is shown as a dotted curve and the location of the gas surface density peak at apoapsis (overdensity of fluid region) is shown as a red dashed line. 
Fig. 9
Impact velocity for collisions in disk around Kepler16 (top row) and Kepler34 (bottom row). 
In Kepler34 seen in panel H of Fig. 4 the effects are far less noticeable since the gas ring is further out and hence less dense. This means the gravitational effect is weaker, but also acts on a smaller number of planetesimals since the number density of planetesimals falls off with distance from the binary. The eccentricity distribution barely diverges from the gas free case with e_{max} being set by the binary forcing and not from the presence of gas. The ring is static and not eccentric, unlike that seen in Kepler16.
Fig. 10
Spatial distribution of collisions recorded in the simulation for the asymmetric gas potential case of both Kepler16 (left) and Kepler34 (right). The position of each collision is rotated by the gas disk phase to produce a relative position. 
Fig. 11
Collision type occurrence for gas free, static axisymmetric gas potential and asymmetric gas potential simulations for a unimodal (upper) and mixed (lower) size distribution disk around Kepler16. White line shows total number of collisions. 
For Kepler16, the gas surface density maximum, shown in Fig. 4, at 1.5 au is efficient at aligning planetesimal pericenters even in the axisymmetric case. As can be seen in Fig. 6, alignment is at a maximum at 1.5 au but begins from 1.2 and continues to the outer edge of the planetesimal disk. The axisymmetric case sees a dephasing of orbits at around 2.0 au, similar to that seen when considering the full asymmetric gas potential.
The formation of a static eccentric planetesimal ring in Kepler16 leads to a similar but suppressed wing effect as compared to the full asymmetric case. Planetesimals orbiting near the ring apoapsis have a reduced orbital velocity and those at periapsis have an increased velocity. Since Kepler34 has a circular ring at the density peak there is no change to the orbital velocities of the planetesimals and they follow the Keplerian velocity. The lack of departure from the Keplerian orbital velocity and the static nature of the planetesimal rings lead to no significant difference to the gas free case when looking at the impact velocities. It it unsurprising then that the collision type occurrence show no significant change from the gas free case for both Kepler16 and Kepler34 as shown in Figs. 11 and 12.
Fig. 12
Collision type occurrence for gas free, static axisymmetric gas potential and asymmetric gas potential simulations for a unimodal (upper) and mixed (lower) size distribution disk around Kepler34. White line shows total number of collisions. 
3.2. Role of size distribution
A common simplification to many models of planetesimal growth is the assumption of a unimodal size distribution. One reason for removing the complication of a realistic, multimodal size distribution is the difficulty in understanding what this should be. This is particularly true for circumbinary disks where strong perturbations on the inner disk possibly invalidate commonly used power laws that apply for our unperturbed solar system. Our simplification to this problem, where planetesimals occupy one of ten discrete size bins, is described in Sect. 2.
Removing the assumption of same size planetesimals is important in our simulations for three main reasons:

1.
Colliders of equal mass will disrupt more easily due to enhanced momentum coupling. This will lead to a bias in the number of erosive collisions. EDACM uses massratio as an input to collision type determination and thus is susceptible to this bias.

2.
Colliders with different sizes will lead to larger range of impact velocities since encounter velocities between planetesimals are set by gravitational focusing as well as binary/gas forces.

3.
EDACM uses impact parameter as a quantifier for collision determination. Equal size colliders are more likely to undergo hitandrun collisions since the critical angle is more likely to be exceeded when the target and projectile have comparable radii. Therefore introducing a more realistic size distribution will help to remove the bias in hitandrun encounters.
In both Kepler16 and Kepler34 (see Figs. 11 and 12) we find that introducing this mixed size distribution reduces the number of supercatastrophic collisions and increases the number of erosive events and erosive hitandruns. This is consistent for all gas implementation cases. This is expected since planetesimals previously classified as involved in highly erosive (supercatastrophic) events are now less easy to disrupt but still have large values of impact parameter from the orbit crossing events. They therefore become relegated to standard erosion and bounce like hitandruns with an erosive element since they erode the smaller projectile in the process.
4. Discussion
In this paper we have performed hybrid simulations of circumbinary protoplanetary disks in an attempt to determine the importance of gas disk gravity on the ability for planetesimals to undergo growthenabling collisions. Our work ties together the pure Nbody planetesimal dynamics and growth of Lines et al. (2014) and the hydrodynamical simulations of circumbinary gas disks of Lines et al. (2015). The work is comparable to that of Marzari et al. (2013) who investigate the effects of gas disk gravity on the ability for planetesimal to accumulate in the Kepler16 system. Our work differs to Marzari et al. (2013) in a number of ways:

1.
We do not consider the effects of thermal evolution on the gas disks structure and evolution.

2.
Our work uses the full implementation of the collision model EDACM Leinhardt & Stewart (2012) which allows for the accurate determination of collision outcomes based on multiple factors, and not just impact velocity.

3.
We consider 10^{5} interacting planetesimals with a radius of R_{p} ≈ 1000 km, as opposed to 400 R_{p} ≈ 5−25 km test particles. We thus probe an entirely different size regime that considers planetesimals to have somehow grown to this size prior to the start of the simulation.
Our results show that it is imperative to include the full precessing asymmetric gas disk potential when considering the evolution of planetesimal dynamics in circumbinary disks. For both systems we have explored, Kepler16 and Kepler34, gravitational effects from the asymmetric gas disk drive up eccentricities with values often exceeding four times that of the background eccentricity set by the dynamical and secular forcing from the binary. The gas surface density, shown in Fig. 1, takes a similar shape for both Kepler16 and Kepler34; an eccentric ring of material with a high density build up around the inner disk apoapsis and an asymmetric cavity interior to this. Once the planetesimals have reached quasisteadystate they begin to adopt this shape as those orbiting within the cavity are drawn outwards by the gravitational attraction of the high density gas disk inner edge. Gas gravity not only focuses planetesimals onto eccentric orbits, but particularly around the density peak also efficiently aligns planetesimal pericenters. This means that orbit crossings are far less common in this region since planetesimals are well phased.
4.1. Kepler16
In Kepler16 there is an order of magnitude drop in impact velocity at around 1.5 au. This corresponds to a significant drop in the proportion of erosive collision types from this orbital radius, as can be seen in Fig. 11. Particularly problematic however is the region interior to this (1.0 au−1.4 au) where planetesimals have largely inflated eccentricities but do not undergo strong orbital alignment. Here, eccentric orbits can overlap contributing to large collision velocities during physical interactions. This corresponds to an overwhelming proportion of erosive collisions that reaches 60% at 1.2 au for the unimodal size case and 70% for the mixed size distribution.
Erosive collisions either match, or more often, surpass the number of growth enabling (perfect merging and partial accretion) collisions from the inner edge of the disk at 0.6 au to the region of efficient planetesimal pericenter alignment at 1.5 au. Beyond this radius collisions quickly become predominantly accretion based, consistently above 50% out to the disk edge at 3.0 au. This would suggest that Kepler16b at 0.7 au could not have formed in situ, supporting previous hypotheses (Lines et al. 2014; Paardekooper et al. 2012; Meschiari 2012a) that circumbinary planets must have formed further away from the binary barycenter and then migrated to their present location. In this case we find that Kepler16b couldn’t have formed interior to 1.5 au.
There is a noticeable change in planetesimal dynamics and collision outcomes when using only a static axisymmetric gas potential. Planetesimals are pulled onto an eccentric ring at the location of the gas density peak, which despite raising eccentricities across the disk does not lead to a change in collision outcomes over the gas free case. This can likely be attributed to the way planetesimal percenters are aligned over most of the disk, which phases planetesimal orbits such that relative velocities and hence collision velocities are kept low.
4.2. Kepler34
Since the density peak in the axisymmetric surface density occurs further out than in Kepler16, it overlays in the planetesimal disk in a region of low number density. This means that the gravitational focusing effect of the gas gravity is weak. For this reason we again find very little modification to the planetesimal velocities and the collision occurrences remain similar to that of the gas free case. For the asymmetric case, similarly because the density peak in Kepler34 occurs further out than in Kepler16, it leads to a more eccentric ring. This results in planetesimal being perturbed by the gas over a greater radial extent of the disk and can be seen to increase eccentricities to e> 0.2 between 1.4 au and 2.2 au (Fig. 4).
Consequently, collisions are dominantly erosive from 1.0 au to 2.0 au with erosion typically accounting for over 50% of all collisions within this range (Fig. 12). Interior to 1.0 au the low gas surface density has little impact on the collisions and instead the encounter velocities and hence collision occurrences are set by the binary forcing. This is clear from Fig. 12 as the high proportion of disruptive collisions from 0.6 au to 0.7 au is also seen in the gas free scenario. In a small window from 0.8 au to 1.0 au erosive collisions are equally matched by growth enabling ones. However the majority of these erosive collisions are supercatastrophic meaning that planetesimals are being ground back down into much smaller bodies and eventually dust. Therefore planetesimals would not only struggle to accrete as much as they lose in this region, but they would likely continue to undergo net erosion. From our results, we suspect that Kepler34b, which lies at 1.1 au could not have formed in situ and would have migrated from outside of 2.0 au.
It is important to note that the collision occurrence values in the gas free case are much more favourable to planetesimal growth than those found in Lines et al. (2014) since the lower resolution of these simulations, and hence increased size of planetesimals, means that planetesimals can grow more easily. This is a simple consequence of larger planetesimals having a larger gravitational binding energy, supporting them against disruption. In this vein it is also important to highlight that our results indicate an extremely hostile disk in the presence of the asymmetric gas gravity despite these large initial planetesimals.
4.3. Selfgravity
The role of interplanetesimal selfgravity should not be underestimated. We have shown in Fig. 7 that including full gravitational effects between disk bodies appears to change, significantly, the eccentricity distribution. In the axisymmetric gas disk case, the high density ring focuses planetesimals onto a narrow annulus which becomes selfsupported by interparticle gravity. This interparticle gravity supported ring then has its eccentricity increased by interaction with the binary. It is not surprising that particle gravity plays such an important role in our simulations as the initial size and mass of our planetesimals is large. Full numerical simulations are necessary to correctly account for gravitational focusing effects and gravitysupported structures. It should be noted however, that planetesimal selfgravity may not be the sole factor in the eccentricity inflation. When the mass of the planetesimal disk is reduced, it is likely the stellar binary ceases to precess, and hence it could be the precession of the binary that drives the differences observed in Fig. 7.
The role of selfgravity was also highlighted in Lines et al. (2014) where collision outcomes were changed by reducing selfgravity that did not model correctly the dynamics between planetesimal on close approach. Therefore we advise caution for future simulations that do not take selfgravity into account for at least large (R_{p}> 100 km) planetesimals.
4.4. Gas drag
Despite the large size of the planetesimals, they still feel gas drag due to the large eccentricities gained through interaction with gas disk gravity. Using a simple analytical prescription for the implementation of aerodynamic drag from an axisymmetric disk, we perform preliminary simulations with drag enabled in an attempt to quantify the relevance of the drag force in future work. We find that the inclusion of axisymmetric drag has only a marginal impact on the dynamics of collisions between planetesimals, with gravitational effects dominating over drag forces by a few orders of magnitude.
5. Summary and further work
Gas feedback on planetesimals is often a missing element from simulations of planetesimal growth in protoplanetary disks, due to the nontrivial nature of unifying these two elements. Our work has succeeded in hybridising 2D hydrodynamical simulations performed with FARGOADSG with 3D, high resolution, interparticle gravityenabled Nbody simulations of circumbinary planetesimal disks. Using the stateoftheart collision model EDACM we have identified regions of disks around Kepler16 and Kepler34 where planetesimal growth can occur. We highlight that neither of these systems could have formed in situ in the presence of a precessing asymmetric disk, even for a disk half the mass of the minimum mass solar nebula. We confine netgrowth regions to r> 1.5 au for Kepler16 and r> 2.0 au for Kepler34 placing both planets at half this critical radius. Our work corroborates that of Marzari et al. (2013) who find that the eccentric gas disk gravity inhibits accretion by increasing planetesimal eccentricities and removing their perihelia alignment, and also Meschiari (2012b) who find that a turbulence gas disk prevents planetesimal accumulation by a similar mechanism.
Several improvements can still be made and will be addressed in future work. These include a) performing further
hydrodynamical simulations to solidify our understanding of the structure and evolution of the gas disk by including the heating and cooling of the gas b) increasing the resolution of the Nbody disk to probe alternate planetesimal size regimes c) improving the handling of collision dust produced during erosive collisions d) performing further investigations into the importance of selfgravitating planetesimals and e) including drag from a precessing asymmetric disk using a method similar to that used for reconstructing the potential in this paper.
Acknowledgments
S.L and Z.M.L are supported by the STFC. P.J.C. is grateful to NERC Grant NE/K004778/1. S.J.P. is supported by a Royal Society University Research Fellowship. The authors acknowledge the University of Bristol Advanced Computing Research Centre facilities (https://www.acrc.bris.ac.uk/), and the use of Bluecrystal Phase III which was used to carry out this work. We thank an anonymous reviewer for their useful suggestions.
References
 Baruteau, C. 2008, Ph.D. Thesis, CEA Saclay, GifsurYvette [Google Scholar]
 Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carter, P. J., Leinhardt, Z. M., Elliott, T., Walter, M. J., & Stewart, S. T. 2015, ApJ, 813, 72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Ida, S., & Makino, J. 1992, Icarus, 96, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Leinhardt, Z. M., Dobinson, J., Carter, P. J., & Lines, S. 2015, ApJ, 806, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Lines, S., Leinhardt, Z. M., Paardekooper, S., Baruteau, C., & Thebault, P. 2014, ApJ, 782, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Lines, S., Leinhardt, Z. M., Baruteau, C., Paardekooper, S.J., & Carter, P. J. 2015, A&A, 582, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lubow, S. H. 1991, ApJ, 381, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., Thebault, P., Scholl, H., Picogna, G., & Baruteau, C. 2013, A&A, 553, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meschiari, S. 2012a, ApJ, 752, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Meschiari, S. 2012b, ApJ, 761, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Moriwaki, K., & Nakagawa, Y. 2004, ApJ, 609, 1065 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paardekooper, S.J., Thébault, P., & Mellema, G. 2008, MNRAS, 386, 973 [NASA ADS] [CrossRef] [Google Scholar]
 Paardekooper, S.J., Leinhardt, Z. M., Thébault, P., & Baruteau, C. 2012, ApJ, 754, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Pelupessy, F. I., & Portegies Zwart, S. 2013, MNRAS, 429, 895 [NASA ADS] [CrossRef] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2013, A&A, 556, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rafikov, R. R. 2013, ApJ, 764, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Richardson, D. C., Quinn, T., Stadel, J., & Lake, G. 2000, Icarus, 143, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Stadel, J. 2001, Ph.D. Thesis, University of Washington [Google Scholar]
 Thébault, P., Marzari, F., & Scholl, H. 2006, Icarus, 183, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
All Tables
All Figures
Fig. 1
Synthetic surface density maps of Kepler16 and Kepler34 created using the method described in Sect. 2. The orientation of the disk is chosen to match its phase at the end of the Nbody simulations in Fig. 4. 

In the text 
Fig. 2
Map of planetesimal accelerations from an Nbody test simulation calculated from the gas disk potential of Kepler34. 50% of planetesimals in the simulation are shown and are plotted as unit vectors with the colour corresponding to their acceleration magnitude. 

In the text 
Fig. 3
Time evolution of planetesimal disk locations (top row) and eccentricity (bottom row) to quasisteadystate around Kepler16, under the influence of a full asymmetric gas potential. 

In the text 
Fig. 4
Planetesimal disks at quasisteadystate (end of simulation at t ≈ 40 000 P_{AB}. AF: Location and eccentricity of planetesimals in disk around Kepler16. GL: Location and eccentricity of planetesimals in disk around Kepler34. 

In the text 
Fig. 5
Change in the distribution of planetesimal mass between initial conditions and simulation end for the asymmetric potential only cases. Planetesimals are placed into bins of width 0.1 au and the increase or decrease of mass in each bin is plotted as a percentage change. 

In the text 
Fig. 6
Longitude of periastron for each planetesimal (black dots) in the Kepler16 system. 

In the text 
Fig. 7
Change in planetesimal eccentricities by reducing planetesimal selfgravity. Snapshot at 14 000 P_{AB} for Kepler16 with an axisymmetric gas disk. 

In the text 
Fig. 8
Planetesimal orbital velocity in Kepler16 (top row) and Kepler34 (bottom row) for both axisymmetric and asymmetric gas disk gravity inclusion. The Keplerian velocity is shown as a dotted curve and the location of the gas surface density peak at apoapsis (overdensity of fluid region) is shown as a red dashed line. 

In the text 
Fig. 9
Impact velocity for collisions in disk around Kepler16 (top row) and Kepler34 (bottom row). 

In the text 
Fig. 10
Spatial distribution of collisions recorded in the simulation for the asymmetric gas potential case of both Kepler16 (left) and Kepler34 (right). The position of each collision is rotated by the gas disk phase to produce a relative position. 

In the text 
Fig. 11
Collision type occurrence for gas free, static axisymmetric gas potential and asymmetric gas potential simulations for a unimodal (upper) and mixed (lower) size distribution disk around Kepler16. White line shows total number of collisions. 

In the text 
Fig. 12
Collision type occurrence for gas free, static axisymmetric gas potential and asymmetric gas potential simulations for a unimodal (upper) and mixed (lower) size distribution disk around Kepler34. White line shows total number of collisions. 

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