Free Access
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/0004-6361/201628349
Published online 11 May 2016

© 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 pumped-up 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 N-body 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, Kepler-16 and Kepler-34, and found that in all cases the inner disk became largely asymmetric with a build-up of gas at the disk apoapsis. Such regions of high surface density could present strong time-dependent gravitational forcing and drag on the planetesimals.

Due to the nature of performing simulations that unify both N-body 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 non-axisymmetric 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 N-body code to perform 3D, high-resolution, inter-particle gravity-enabled 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 N-body simulations are conducted using the widely adopted code PKDGRAV (Stadel 2001; Richardson et al. 2000). PKDGRAV is a highly parallelised, multi-disciplinary code capable of handling large N. An efficient tree code and multipole expansion allows for Nlog N scaling and hence the ability to run self-gravitating simulations with N105 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, non-erosive/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 semi-analytical approach to integrating the gas disk into the N-body code. The surface density data from 2D hydrodynamical simulations of circumbinary gas disks around Kepler-16 and Kepler-34 performed in Paper I are used.

thumbnail Fig. 1

Synthetic surface density maps of Kepler-16 and Kepler-34 created using the method described in Sect. 2. The orientation of the disk is chosen to match its phase at the end of the N-body simulations in Fig. 4.

Open with DEXTER

The quasi-steady-state (QSS) gas surface density data from Lines 2015 (simulations runs B and F) are used. The simulations correspond to non-self-gravitating disks since, for the disk masses explored here, the difference in gas surface density between self- and non-self-gravitating 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 self-gravity. These fluid parameters are typical of circumbinary disk simulations. Self-gravitating 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/au2. These hydrodynamical simulations are performed using FARGO-ADSG (Baruteau & Masset 2008; Baruteau 2008) over a polar mesh grid with Nr = 395 and Ns = 512 where Nr and Ns are the number of radial and azimuthal cells respectively. In Kepler-16 the gas disk has precession period of Pd = 2000 PAB and in Kepler-34 has Pd = 3000 PAB yielding Pd ≈ 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:

thumbnail Fig. 2

Map of planetesimal accelerations from an N-body test simulation calculated from the gas disk potential of Kepler-34. 50% of planetesimals in the simulation are shown and are plotted as unit vectors with the colour corresponding to their acceleration magnitude.

Open with DEXTER

  • 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 ∝ (r2 + ϵ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.

    thumbnail Fig. 3

    Time evolution of planetesimal disk locations (top row) and eccentricity (bottom row) to quasi-steady-state around Kepler-16, under the influence of a full asymmetric gas potential.

    Open with DEXTER

  • 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 inter-planetesimal and stellar-planetesimal interactions. The result of this implementation is shown in Fig. 2 where the planetesimal accelerations in response to the potential of a Kepler-34 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.

Table 1

Parameter setup of each simulation A through L.

thumbnail Fig. 4

Planetesimal disks at quasi-steady-state (end of simulation at t ≈ 40 000 PAB. A-F: Location and eccentricity of planetesimals in disk around Kepler-16. G-L: Location and eccentricity of planetesimals in disk around Kepler-34.

Open with DEXTER

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

Open with DEXTER

2.2. Initial conditions

We test the response of a circumbinary planetesimal disk around two different stellar binary systems, Kepler-16 (Doyle et al. 2011) and Kepler-34 (Welsh et al. 2012) which have a total mass of 0.89 M and 2.07 M respectively. Kepler-34 has MA = 1.05 M and MB = 1.02 M where MA and MB are the primary and secondary stellar mass, stellar separation is ab = 0.22 au and binary eccentricity is eb = 0.53. Kepler-16 has MA = 0.69 M and MB = 0.20 M, ab = 0.22 au and eb = 0.16. Each planetesimal disk has a total mass of 2.8 M spread over N = 105 planetesimals and the disk domain ranges from Rin = 0.6 to Rout = 3.0 au. The planetesimal density is 2.0 g/cm3 to match that of silicon-rich rocky bodies. Two different models are used for the planetesimal disk size distribution: unimodal and mixed. For the unimodal size distribution, mp = 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 mp = 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 N-body resolution or reducing Mdisk 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 N-body resolution for this work is necessary to evolve the disk to a quasi-steady-state in a practical timeframe. Such large planetesimals, which may be able to form at this size through fast clumping assisted by streaming-instabilities (Johansen et al. 2007; Carrera et al. 2015), provide a convenient best-case 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 N-body 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 N-body 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 e2 = 2 ⟨i2 = 0.007 (Ida & Makino 1992). We enable collisions, and gas gravity if applicable, from t = 0 and allow the disk to reach a quasi-steady-state. Only collisions that occur from after a quasi steady-state are used in the analysis and results. Each simulation runs for 3800 years in total (34 000 PAB and 50 000 PAB for Kepler-16 and Kepler-34 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 over-density 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 self-gravitating nature.

Figure 3 shows the formation and evolution of this feature along with the associated planetesimal eccentricities, for an asymmetric gas potential on a Kepler-16 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 PAB 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 PAB. 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 emax = 0.4. This is four times that found in the gas free simulation (emax = 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.

thumbnail Fig. 6

Longitude of periastron for each planetesimal (black dots) in the Kepler-16 system.

Open with DEXTER

Figure 5 assists in understanding the change in planetesimal mass distribution under the influence of the gas disk gravity. In Kepler-16 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 Kepler-34 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 Kepler-16 and Kepler-34. 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.

thumbnail Fig. 7

Change in planetesimal eccentricities by reducing planetesimal self-gravity. Snapshot at 14 000 PAB for Kepler-16 with an axisymmetric gas disk.

Open with DEXTER

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 no-collision 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 Kepler-16 in Fig. 11 and for Kepler-34 in Fig. 12. For Kepler-16 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 Kepler-16b 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 Kepler-34b 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 Kepler-16, 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 emax = 0.18 is much smaller than that found in the full asymmetric case.

An interesting observation is that in Kepler-16, 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, Kepler-16) but reducing the mass of planetesimal by a factor of 105 (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 self-gravitation of our planetesimals has a large contribution to the dynamics. Not only does reducing inter-particle gravity remove the high eccentricities at the disk outer edge, but it also removes the peak associated with the high density gas disk ring.

thumbnail Fig. 8

Planetesimal orbital velocity in Kepler-16 (top row) and Kepler-34 (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 (over-density of fluid region) is shown as a red dashed line.

Open with DEXTER

thumbnail Fig. 9

Impact velocity for collisions in disk around Kepler-16 (top row) and Kepler-34 (bottom row).

Open with DEXTER

In Kepler-34 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 emax being set by the binary forcing and not from the presence of gas. The ring is static and not eccentric, unlike that seen in Kepler-16.

thumbnail Fig. 10

Spatial distribution of collisions recorded in the simulation for the asymmetric gas potential case of both Kepler-16 (left) and Kepler-34 (right). The position of each collision is rotated by the gas disk phase to produce a relative position.

Open with DEXTER

thumbnail 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 Kepler-16. White line shows total number of collisions.

Open with DEXTER

For Kepler-16, 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 Kepler-16 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 Kepler-34 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 Kepler-16 and Kepler-34 as shown in Figs. 11 and 12.

thumbnail 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 Kepler-34. White line shows total number of collisions.

Open with DEXTER

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, multi-modal 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 mass-ratio 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 hit-and-run 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 hit-and-run encounters.

In both Kepler-16 and Kepler-34 (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 hit-and-runs. 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 hit-and-runs 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 growth-enabling collisions. Our work ties together the pure N-body 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 Kepler-16 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 105 interacting planetesimals with a radius of Rp ≈ 1000 km, as opposed to 400 Rp 525 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, Kepler-16 and Kepler-34, 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 Kepler-16 and Kepler-34; 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 quasi-steady-state 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. Kepler-16

In Kepler-16 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 au1.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 Kepler-16b 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 Kepler-16b 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. Kepler-34

Since the density peak in the axisymmetric surface density occurs further out than in Kepler-16, 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 Kepler-34 occurs further out than in Kepler-16, 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 Kepler-34b, 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. Self-gravity

The role of inter-planetesimal self-gravity 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 self-supported by inter-particle gravity. This inter-particle 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 gravity-supported structures. It should be noted however, that planetesimal self-gravity 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 self-gravity was also highlighted in Lines et al. (2014) where collision outcomes were changed by reducing self-gravity that did not model correctly the dynamics between planetesimal on close approach. Therefore we advise caution for future simulations that do not take self-gravity into account for at least large (Rp> 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 non-trivial nature of unifying these two elements. Our work has succeeded in hybridising 2D hydrodynamical simulations performed with FARGO-ADSG with 3D, high resolution, inter-particle gravity-enabled N-body simulations of circumbinary planetesimal disks. Using the state-of-the-art collision model EDACM we have identified regions of disks around Kepler-16 and Kepler-34 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 net-growth regions to r> 1.5 au for Kepler-16 and r> 2.0 au for Kepler-34 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 N-body 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 self-gravitating 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

All Tables

Table 1

Parameter setup of each simulation A through L.

All Figures

thumbnail Fig. 1

Synthetic surface density maps of Kepler-16 and Kepler-34 created using the method described in Sect. 2. The orientation of the disk is chosen to match its phase at the end of the N-body simulations in Fig. 4.

Open with DEXTER
In the text
thumbnail Fig. 2

Map of planetesimal accelerations from an N-body test simulation calculated from the gas disk potential of Kepler-34. 50% of planetesimals in the simulation are shown and are plotted as unit vectors with the colour corresponding to their acceleration magnitude.

Open with DEXTER
In the text
thumbnail Fig. 3

Time evolution of planetesimal disk locations (top row) and eccentricity (bottom row) to quasi-steady-state around Kepler-16, under the influence of a full asymmetric gas potential.

Open with DEXTER
In the text
thumbnail Fig. 4

Planetesimal disks at quasi-steady-state (end of simulation at t ≈ 40 000 PAB. A-F: Location and eccentricity of planetesimals in disk around Kepler-16. G-L: Location and eccentricity of planetesimals in disk around Kepler-34.

Open with DEXTER
In the text
thumbnail 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.

Open with DEXTER
In the text
thumbnail Fig. 6

Longitude of periastron for each planetesimal (black dots) in the Kepler-16 system.

Open with DEXTER
In the text
thumbnail Fig. 7

Change in planetesimal eccentricities by reducing planetesimal self-gravity. Snapshot at 14 000 PAB for Kepler-16 with an axisymmetric gas disk.

Open with DEXTER
In the text
thumbnail Fig. 8

Planetesimal orbital velocity in Kepler-16 (top row) and Kepler-34 (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 (over-density of fluid region) is shown as a red dashed line.

Open with DEXTER
In the text
thumbnail Fig. 9

Impact velocity for collisions in disk around Kepler-16 (top row) and Kepler-34 (bottom row).

Open with DEXTER
In the text
thumbnail Fig. 10

Spatial distribution of collisions recorded in the simulation for the asymmetric gas potential case of both Kepler-16 (left) and Kepler-34 (right). The position of each collision is rotated by the gas disk phase to produce a relative position.

Open with DEXTER
In the text
thumbnail 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 Kepler-16. White line shows total number of collisions.

Open with DEXTER
In the text
thumbnail 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 Kepler-34. White line shows total number of collisions.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.