# Modelling circumbinary protoplanetary disks

###### Key Words.:

methods: numerical – hydrodynamics – planets and satellites: formation – protoplanetary disks – binaries: close###### Abstract

Context:The Kepler mission’s discovery of a number of circumbinary planets orbiting close ( 1.1 au) to the stellar binary raises questions as to how these planets could have formed given the intense gravitational perturbations the dual stars impart on the disk. The gas component of circumbinary protoplanetary disks is perturbed in a similar manner to the solid, planetesimal dominated counterpart, although the mechanism by which disk eccentricity originates differs.

Aims:This is the first work of a series that aims to investigate the conditions for planet formation in circumbinary protoplanetary disks.

Methods:We present a number of hydrodynamical simulations that explore the response of gas disks around two observed binary systems: Kepler-16 and Kepler-34. We probe the importance of disk viscosity, aspect-ratio, inner boundary condition, initial surface density gradient, and self-gravity on the dynamical evolution of the disk, as well as its quasi-steady-state profile.

Results:We find there is a strong influence of binary type on the mean disk eccentricity, , leading to for Kepler-16 and in Kepler-34. The value of -viscosity has little influence on the disk, but we find a strong increase in mean disk eccentricity with increasing aspect-ratio due to wave propagation effects. The choice of inner boundary condition only has a small effect on the surface density and eccentricity of the disk. Our primary finding is that including disk self-gravity has little impact on the evolution or final state of the disk for disks with masses less than 12.5 times that of the minimum-mass solar nebula. This finding contrasts with the results of self-gravity relevance in circumprimary disks, where its inclusion is found to be an important factor in describing the disk evolution.

Conclusions:

## 1 Introduction

Extrasolar planets in circumbinary (p-type) orbits around short-period stellar binary systems are a recent addition to the increasingly diverse range of planetary characteristics found outside our solar system. We now know of nine short-period planets^{1}^{1}1Kepler-16(AB)b (Doyle et al., 2011), Kepler-34(AB)b 35(AB)b (Welsh et al., 2012), Kepler-38(AB)b (Orosz et al., 2012b), Kepler-47(AB)b,c (Orosz et al., 2012a), Kepler-64(AB)b (Schwamb et al., 2013), Kepler-413(AB)b (Kostov et al., 2014) and KIC 9632895 (Welsh et al., 2014). with semi-major axis 1.1 au, placing them close to the binary barycentre and under the influence of perturbative forces from the secondary star.

The stability of circumbinary planet orbits was studied well before the first confirmed exoplanet detection in 1992 (Wolszczan & Frail, 1992) by Heppenheimer (1978). More recently, Holman & Wiegert (1999) performed a parameter space study of the long-term stability of planets in P-type orbits. They calculated a range of orbital radii that allows for dynamical stability for a range of binary eccentricities and semi-major axes. Interestingly, with the exception of Kepler-47(AB)c, all known circumbinary planets lie just outside their inner most stable orbit, . One plausible explanation for this occurrence is that at some point in their evolution, these planets underwent a period of migration that was halted.

Pierens & Nelson (2008) showed that Saturn- and Jupiter-mass planets undergo different interactions with the binary that can result in either stabilisation of the planetary orbits or its removal from the system via ejection or scattering. Specifically they find that after a period of inwards migration, Saturn-mass planets have a stable evolution. A torque reversal, caused by an eccentricity increase from interaction with the binary, leads to stable outwards migration. However, Jupiter-mass planets often become trapped in 4:1 mean motion resonance with the binary, which causes significant increases in planetary eccentricity (Nelson, 2003). Eventually, close encounters with the stars at periapsis lead to outward scattering or complete ejection. The latter effect may explain why all Kepler circumbinary planets discovered so far are sub-Jupiter mass, a result difficult to explain by observational biases since larger planets are easier to detect.

Simulation | Binary System | () | Self-Gravity | Boundary | |||

A | Kepler-16 | 0.03 | No | Rigid | 1.5 | ||

B | Kepler-16 | 0.05 | No | Rigid | 1.5 | ||

C | Kepler-16 | 0.03 | No | Rigid | 1.5 | ||

D | Kepler-16 | 0.05 | No | Rigid | 1.5 | ||

E | Kepler-34 | 0.03 | No | Rigid | 1.5 | ||

F | Kepler-34 | 0.05 | No | Rigid | 1.5 | ||

G | Kepler-34 | 0.03 | No | Rigid | 1.5 | ||

H | Kepler-34 | 0.05 | No | Rigid | 1.5 | ||

I | Kepler-16 | 0.03 | Yes | Rigid | 1.5 | ||

J | Kepler-16 | 0.05 | Yes | Rigid | 1.5 | ||

K | Kepler-34 | 0.03 | Yes | Rigid | 1.5 | ||

L | Kepler-34 | 0.05 | Yes | Rigid | 1.5 | ||

M | Kepler-16 | 0.05 | No | Open | 1.5 | ||

N | Kepler-34 | 0.05 | No | Open | 1.5 | ||

O | Kepler-16 | 0.05 | Yes | Rigid | 1.5 | ||

P | Kepler-16 | 0.05 | Yes | Rigid | 1.5 | ||

Q | Kepler-16 | 0.05 | Yes | Rigid | 1.5 | ||

R | Kepler-16 | 0.05 | No | Rigid | 1.5 | ||

S | Kepler-16 | 0.05 | No | Rigid | 1.5 | ||

T | Kepler-16 | 0.05 | No | Rigid | 1.5 | ||

U | Kepler-16 | 0.05 | Yes | Rigid | 0.5 |

Although all the observed circumbinary planets are currently in stable orbits, it is not obvious that they formed in situ. Several studies have been carried out to try to understand how the processes behind planet formation and evolution can occur in spite of an intense gravitational influence on the protoplanetary disk from the central binary. One of the key hurdles in forming solid cores and terrestrial planets in circumbinary disks is explaining how rocky, metre- to kilometre-sized planetesimals undergo collisions that result in growth rather than erosion. It has been shown in previous work that eccentricity forcing from the binary on the protoplanetary disk can drive up relative velocities which not only inhibits accretion but can lead to energetic and highly erosive impacts. In particular, Lines et al. (2014) showed that even massive planetesimals are often disrupted despite having a large gravitational binding energy. Their work corroborates that of Paardekooper et al. (2012) and Meschiari (2012a, b) suggesting that sustained planetesimal accretion is unlikely and thus, the observed Kepler circumbinary planets did not form in situ. The most likely explanation for the presence of these planets is the migration of planetary cores that form at larger orbital radii where the gravitational forcing from the binary is of significantly lower magnitude.

A key component missing from many of the N-body studies of planetesimal evolution in circumbinary disks is the inclusion of a gas disk which is typically 100 times the mass of its solid counterpart. As is the case with the rocky planetesimals, the gas is also perturbed by the binary, raising the eccentricity of the disk. The eccentricity in the gas disk is established through parametric instabilities originating from non-linear mode coupling between the eccentric modes of both the initial disk eccentricity and the binary potential (Papaloizou, 2002; Hayasaki & Okazaki, 2009; Lubow, 1991) and well as the direct driving from the binary potential (Lubow & Artymowicz, 2000). This coupling produces an spiral wave at the 3:1 Lindblad resonance which increases disk eccentricity via the outwards transport of angular momentum. The overall result is an eccentric, asymmetric disk precessing about the binary (Kley & Haghighipour, 2014; Marzari et al., 2013; Pelupessy & Portegies Zwart, 2013; Pierens & Nelson, 2013, (hereafter referred to as PN13)). The asymmetric gas disk interacts with the solid planetesimal disk through gas drag and the time dependent gravitational potential of the gas disk. The latter effect is particularly important when large asymmetries in the gas surface density arise (see Figure 1).

The role of gas drag has been probed in several studies which found that differential orbital phasing caused relative velocities to remain small between planetesimals of comparative mass, thus, orbit crossing will occur between planetesimals of different sizes (Scholl et al., 2007; Marzari et al., 2008). The influence of disk self-gravity on the evolution of the system is less well explored. Recent studies such as PN13 and Kley & Haghighipour (2014) choose not to include self-gravity in their simulations. The former justifies this approximation by confirming the Toomre parameter,

(1) |

where is the sound speed, is the angular velocity, is the gravitational constant, and is the surface density, satisfies the Toomre stability criterion () at all disk radii for their chosen disk mass. However, even in self-gravitating disks of relatively low mass, low-frequency global modes exist that do not have a counterpart in non self-gravitating disks (Papaloizou, 2002). Moreover, the study by Marzari et al. (2009) shows that self-gravity in circumprimary disks can be important even in low-mass disks. Therefore, an investigation into whether or not self-gravity is required to fully describe a circumbinary disk, even at lower masses, is required.

In this paper, we aim to explore a number of parameters characterising a gaseous circumbinary disk in an attempt to find a suitable surface density profile for both the Kepler-16 and Kepler-34 systems. The resulting quasi-steady-state surface density profiles will be used to parameterize the gas disk potential for use in our -body simulations of planetesimal growth. This next work will be presented in a second paper, part II. Section 2 presents the numerical method and model parameters. In Section 3 and 4 we discuss the results and implications of the 21 simulations of the Kepler-16 and Kepler-34 systems. Finally, we summarise our findings in the conclusion in Section 5.

## 2 Numerical methods

The central stellar binary potential is calculated in a fixed steady state orbit centred on the binary barycentre, meaning that the stars do not respond to the gas disk (Pierens & Nelson, 2007). The stars, with mass ratio are first initialised at their respective apoapses, such that the Cartesian coordinates are initially zero in the direction for both stars, = 0, and the position is

(2) |

and

(3) |

where and are the mass of star and , respectively, and the distance of a star from its focus, r, is determined from the shape equation

(4) |

where is the binary semi-major axis, is the binary eccentricity, and the angle of the star from apoapsis is . Using Kepler’s 2nd law to determine the fractional area swept out at any given time, an iterative procedure is invoked to calculate and thus, the position vectors of each star. The stellar binary parameters are given in Table 2.

To simulate the hydrodynamical evolution of the circumbinary disks we use a modified version of the Eulerian fluid code FARGO (Masset, 2000), called FARGO-ADSG (Baruteau & Masset, 2008), which includes optional disk self-gravity. The hydrodynamical equations are solved on a polar mesh which is composed of = 512 equally divided azimuthal sector cells, and = 395 logarithmically spaced radial cells. Logarithmic spacing in the radial direction is a necessary condition for the fast fourier transform algorithm in the self-gravity calculation, but is also a preference for achieving the highest resolution closest to the binary (Baruteau & Masset, 2008).

The majority of our runs use a rigid, reflecting inner boundary condition which is fixed at 0.345 au with the rigid outer boundary at 4.0 au. In a couple of simulations we have used standard outflow boundaries at the grid’s inner and outer edges, where the zero-gradient condition has been extended to the gas azimuthal velocity. This is a necessary implementation since there is no well defined equilibrium between the central gravity from the binary and the opposing centrifugal force and pressure gradient. In all simulations the inner disk edge is truncated by the binary causing a steep positive density gradient which is numerically difficult to model smoothly, thus, we employ a gap function, from Günther et al. (2004) to better model the initial surface density profile in the inner disk region:

(5) |

where is the estimated size of the gap (Artymowicz & Lubow, 1994) and is the orbital radius.

The simulated gas disk is simplified by using the isothermal disk approximation, thereby minimising the number of unknowns. We leave the inclusion of an energy equation to account for the thermal evolution of the disk for a future paper.

System | (au) | ||
---|---|---|---|

Kepler-16 | 0.29 | 0.16 | |

Kepler-34 | 1.0 | 0.52 |

The initial surface density of the gas disk follows from PN13, where is the surface density at 1 au, is the radial position in au, and defines the initial surface density gradient. For the majority of the simulations (Table 1 simulations - & ) we assume a half minimum-mass solar nebula, thus, M/au (Hayashi, 1981). In simulations - is increased by a factor of 5, 10 and 25. The density gradient, is set to 1.5 in all cases except the final simulation, , where a shallower value of 0.5 is used for comparison to previous work (Marzari et al., 2013). In all simulations, the total central stellar mass is normalised to 1.0 M but the surface density of the disk is not scaled to account for the physical mass difference between Kepler-34 and Kepler-16. Therefore, although the total physical mass of the disk is fixed, the disk mass relative to combined stellar mass changes.

To avoid numerical instabilities caused by extremely low density fluid cells in the inner cavity due to the torque of the binary on the disk inner edge, a density floor is added such that the minimum value, , is set to 10 of the initial density. This means that mass is not strictly conserved in regions of very low density.

The aspect ratio, , where is the pressure scale height, is constant across the disk and varies between each model from 0.03 to 0.05. The aspect ratio does not determine the physical thickness of the disk, since the simulations are two dimensional, but defines the sound speed, , where is the Keplerian speed. Turbulence in the disk, a likely contribution of the MRI effect (Balbus & Hawley, 1991), is explored by varying an -viscosity between 10 and 10.

Since our disks are low in mass, we expect them all to be linearly stable according to the Toomre stability criterion, 1 (Toomre, 1964). For our half-MMSN density models, even near the inner edge of the disk where the density is at a maximum, the minimum Toomre value, , is 140. However, on consideration of the results of Marzari et al. (2009) that show self-gravity plays a significant role in circumprimary disks, self-gravity is included in simulations - and -.

Over 21 simulations, we explore the effect that changing the stellar binary parameters, aspect ratio, -viscosity, boundary condition and inclusion of disk self-gravity has on the evolution and quasi steady-state (QSS) profile of the disk. All simulations are evolved for 16,000 binary orbits, , equivalent to 1,700 orbits at 1 au, with the exception of which run for 4,000 to test the effect of the chosen inner boundary condition on the initial response of the disk.

## 3 Results

### 3.1 Non self-gravitating disks

#### 3.1.1 Kepler-16

We first look at the response of a non self-gravitating disk to a stellar binary with parameters chosen to match the observed properties of the system Kepler-16(AB). The azimuthally and time averaged quasi-steady state surface density is shown in Figure 3a. The density maximum or peak, , indicates the location of the truncated inner edge, although the eccentric morphology of the disk requires a full 2-dimensional density map, (Figure 1), to identify the edge precisely as a function of azimuth angle. Therefore, we will consider the location of to be at the average truncation radius, . The value of lies between 1.0 and 1.4 au which is consistent with what was found by PN13.

Increasing the aspect ratio from 0.03 to 0.05 raises the disk eccentricity (Figures 1c & 1e), particularly for , which can also be seen in the surface density profiles by an increase in the value of . The more circular form of the disk in model leads to a higher averaged surface density, since more mass orbits at a similar orbital radius.

The eccentric nature of the disk is better shown by the time evolution of the mean disk eccentricity (Figure 3c), , which we define according to Pierens & Nelson (2007) as

(6) |

where and are the disk inner and outer radii and and are the fluid element surface density and eccentricity respectively. Note is calculated assuming the cell is orbiting the binary barycentre alone (a two body problem). Indeed, models with = 0.05 show a higher value of at times both early in evolution and quasi-steady state (QSS), particularly for the low viscosity run. Our final disk eccentricity values lie between 0.03 and 0.08, the spread matching closely that found in PN13. However, we observe much larger values of eccentricity oscillation amplitude at QSS, which varies between models from to as opposed to a consistent value of found by PN13. This may be due to our smaller value of initial surface density. The increase in eccentricity with aspect ratio is also seen in the radial profiles of the disk’s eccentricity (Fig. 3e). The models with a larger aspect ratio (models and ) have a larger average eccentricity than those with a smaller aspect ratio (models and ) from the inner boundary until about 2.5 au.

We can also calculate the disk mean longitude of periastron, , from the mass weighted average of the cell longitudes, :

(7) |

The value oscillates about the binary periastron which itself does not change due to the analytic prescription for the stellar orbits. In Figure 4, is plot as a function of time for Kepler-16. On cross examination with the eccentricity time evolution (Fig. 3c), it is clear that there is some relation between the eccentricity and periastron longitude since the maxima in the eccentricity oscillations matches the periastron longitude maxima. We discuss this link further in Section 4. The position angle of the disk centre-of-mass is calculated in the inertial frame and shows the circulation of the disk increasing in frequency until the disk precession period, , settles at around 2,500 at QSS.

The eccentric morphology of the disk is apparent through both (Figure 1) and . Figure 2 shows the radial velocity which reveals a number of spiral features. In both the = 0.03 and 0.05 models of Kepler-16, the velocity map shows spiral waves with a larger radial extent for the higher aspect ratio. To better assess the level of asymmetry and presence of modes in the disk, a Fourier analysis of the surface density is performed.

The disk surface density can be decomposed into modes described by azimuthal mode numbers and frequency mode numbers since disk disturbances caused by the binary can have both an angular and time dependency. The surface density distribution can then be written as (Nixon & Lubow, 2015)

(8) |

where is a complex function and is the mean motion of the binary, . Nixon & Lubow (2015) find that eccentric binaries produce modes with . In particular, in their retrograde circumbinary disk simulations, highly eccentric binaries () produce powerful , modes. We remove the time dependence of the disk disturbances by since our analysis is concerned with the eccentric Kepler-16 system ( = 0.16) and because we focus only on a comparative study of the azimuthal mode strengths between simulations. We do this by choosing to discuss modes only with . This is done by performing the Fourier analysis of the surface density averaged over the period of a single binary orbit at QSS. This mode analysis allows for the identification of the mode propagation and dissipation properties with changing disk parameters.

In Figure 5 the strengths of the and modes relative to the axisymmetric component is plot as a function of orbital radius. The results confirm that for the high aspect ratio ( = 0.05) models and , the contribution from the spirals is more significant in the outer disk than that seen for = 0.03.

#### 3.1.2 Kepler-34

Disks surrounding Kepler-34 are significantly more affected by the binary, with the QSS ranging from 0.1 to 0.15, typically three times more eccentric than Kepler-16 (Figure 3d). Since the mass-weighted eccentricities decrease slightly with time by the simulation end, it is not clear from if the disks have reached a steady state by 16,000 . However, looking at the instantaneous surface density profiles for a Kepler-34 run during the last few hundred binary orbits as shown in Figure 6, it appears as though a steady-state has been reached. Our Kepler-34 results agree roughly with PN13, who find that typically, for the systems mutually covered, a quasi steady-state is achieved by .

As with Kepler-16 the average eccentricity in the Kepler-34 disks once they have reached QSS is larger for higher values of aspect ratio. The large difference in ( = 0.05) between models and is shown in Figure 3f to originate from the disk beyond the truncation radius where the surface density is much higher, and is therefore expected since is a mass-weighted value. From 2.0 au, model has a small positive eccentricity offset until 3.0 au, but there is very little difference between all models interior to this location. There is a high level of agreement between aspect ratio and -viscosity models in the surface density profiles too; our simulations see sharing a similar value of 2.0 au across all models. These results agree with PN13.

### 3.2 Self-gravitating disks

In models we enable disk self-gravity (DSG) for both Kepler-16 and Kepler-34. The -viscosity is fixed at but the aspect ratio is varied between 0.03 and 0.05 to adjust the Toomre value and allow for direct comparison with models , , and . For Kepler-16, as seen in Figure 7, during the first 2,000 binary orbits the disk responds almost identically to that of its non self-gravitating counterpart. The evolution up to QSS is subtly different however. The frequency of the eccentricity oscillations increases with gravity enabled suggesting that the disk circulation frequency increases. Due to the absence of strongly defined eccentricity oscillations in the Kepler-34 simulations, it is possible only to differentiate between self- and non self-gravitating models by the eccentricity magnitude, and not the oscillations. We find a slightly reduced disk eccentricity with DSG enabled for the larger aspect ratio simulation. In Figure 7 the time averaged surface density profiles also reveal that a self-gravitating disk in Kepler-16 has no effect on the final disk structure, and almost no effect for disks around Kepler-34.

To test at what disk mass self-gravity become important, is increased by 5, 10 and 25 times. Self-gravity importance is therefore tested for 1 but 1 at all times. We define standard density models ( = 1 10 code units) as being 0.5 and those scaled up by a factor of 5, 10 and 25 as 2.5 ( = 28), 5 ( = 14) and 12.5 ( = 5) respectively.

In Figure 8 the evolution of the mean disk eccentricity for each of the four varying densities is shown. remains, for the non self-gravitating runs, at a constant 0.08 for each surface density. When self-gravity is enabled, with the exception of 12.5 , falls off slowly with time and the rate of this fall off increases with increasing surface density. This leads to disks with higher surface densities having a slightly (5-25%) smaller eccentricity at the end of the simulation. The evolution of the eccentricity in the self-gravitating 12.5 simulation is more dynamic due to an eccentricity fall off before increasing again.

At standard densities we find that the oscillation frequency increases and amplitude decreases with self-gravity enabled. Increasing the density does not change the oscillation frequency without self-gravity enabled, but the period of the oscillations decreases from 2000 (0.5 ) to 300 (12.5 ) when it is included.

Another measure of the effect of self-gravity is to identify the strength of the modes present in the disk. A fourier analysis of the surface density, identical to that done for the non self-gravitating Kepler-16 runs, is done. The Fourier coefficients are determined and normalised against the axisymmetric () component. The strength of these modes normalised to the axisymmetric part, , is shown for each surface density in Figure 9. Again, a comparison is made between simulations with and without self-gravity included.

There is a strong contribution from both the and modes, with the mode always 20 or larger in the inner disk. It is less clear as to how the inclusion of self-gravity changes the strength of these modes. When self-gravity is enabled in the 0.5 , 2.5 and 5 surface density runs, there are only subtle differences between the strength of the modes. At 12.5 , due to oscillations in the strength of both modes in the self-gravitating disk, the profiles for both the self-gravitating and modes do not trace their non self-gravitating equivalents.

### 3.3 Surface density profile gradient

The gradient of the initial density profile is a value which often varies between studies, justified by the lack of observational data. In simulation the density gradient is made shallower by changing the surface density exponent from 1.5 to 0.5. In Figure 10, the radially varying surface density is shown alongside the time evolution of the mean disk eccentricity. Changing the exponent has little effect on the truncation of the inner disk edge since adopts the same value. However, there is a significant difference in the mean disk eccentricity at all times, with = 0.5 obtaining a much lower eccentricity throughout.

### 3.4 Boundary conditions

Previous work on circumbinary disks sees the implementation of two boundary conditions for the inner edge; namely rigid/reflecting and open. Open boundary conditions are chosen to be the more realistic scenario since matter can flow through the boundary and accrete onto the star(s) as would occur in a real system. The difficulty in these cases however is setting a well defined value for and . Rigid boundary conditions are chosen to remove the enhanced mass loss experienced through the inner disk from an eccentric disk. The rigid approximation is, however, unrealistic and is often modified to include routines that remove reflected waves at the boundaries. We explore the influence of the choice of boundary condition on the disk by comparing the surface density and eccentricity between identical simulations, but with different inner boundary conditions.

We find that there are small differences in both the disk eccentricity and structure between the two. In Kepler-16 the surface density profiles show that including an open boundary keeps the peak density at the same location but increases the mass interior to . This effect, which is shown in Figure 11 is stronger in Kepler-34 where the density plateaus at around 1.1 au to 1.5 au before further increasing to a peak at 2 au. We present our explanation for this in the Section 4. The disk eccentricity reveals less of a difference between the two boundary conditions. In the main part of the disk the eccentricity is the same regardless of boundary condition although the mass weighted eccentricity is higher in the cavity.

## 4 Discussion

In this paper, we have shown the results of a series of hydrodynamical simulations that explore the response of a circumbinary gas disk to both the Kepler-16 and Kepler-34 stellar binary systems. We vary -viscosity and aspect ratio, as well as probing the effects of disk self-gravity, initial surface density profile, and the inner boundary condition.

The majority of our simulations have reached steady state (a consistent value of ) by 16,000 . The steady-state nature of the disk is also reflected in the stationary form of the instantaneous, azimuthally averaged surface density during the last few thousand binary orbits (Figure 6), showing the balance between the truncation effects from tidal torques on the inner edge of the disk and the viscous replenishment of mass to the inner disk. The formation of the central cavity, which is shown by the movement of in Figure 1, causes the inner edge of the disk to move outwards. By taking the cavity edge to be the location at which the positive density gradient begins, we show good agreement with the Artymowicz & Lubow (1994) gap estimate .

The results for the surface densities and eccentricities of our non self-gravitating disks agree with those found by PN13. Subtle differences such as the eccentricity magnitude (Figure 3), particularly in Kepler-34, may be explained by a number of factors that include our reduced outer boundary (from 5 au to 4 au) and the difference in initial surface density profiles. We certainly retrieve higher disk eccentricities than found by Kley & Haghighipour (2014) who attribute their lower values in comparison with PN13 to the choice of boundary condition. Kley & Haghighipour (2014) claim that the rigid boundary condition used by PN13 leads to a higher disk eccentricity. We find that the choice of a rigid boundary condition leads to only the slightest increase in disk eccentricity over the open boundary scenario in Kepler-34, and leads to lower disk eccentricities for ¡ 1.5 au in Kepler-16. Therefore, we suggest that it might be more appropriate to link the discrepancy in disk eccentricity to the choice of stellar binary for which we have shown significantly affects the disk structure and dynamics.

The shape and strength of the binaries potential invariably changes with its orbital properties: mass ratio, eccentricity and semi-major axis. Therefore, it is likely that a circumbinary gas disk is only weakly affected by the relatively low eccentricity of Kepler-38 (e = 0.10) (Orosz et al., 2012b) compared to the significantly higher eccentricity of Kepler-34 (e = 0.53). We would therefore expect the gas disk around Kepler-34 to adopt a higher eccentricity than a similar disk around Kepler-38. This theory is supported by the significantly larger disk eccentricities in simulations around the eccentric Kepler-34 compared to Kepler-16 where the stellar binary has a much smaller eccentricity.

Increasing disk eccentricity with increasing stellar binary eccentricity is a trend also seen in -body results. The planetesimal response in the N-body scenario is dependent upon the of stellar binary on the disk, and perturbation theory provides a value of eccentricity of which planetesimals tend towards. This forced eccentricity, , is given by (Moriwaki & Nakagawa, 2004)

(9) |

where and are the primary and secondary stellar masses, is the total binary mass, is the distance from the binary barycentre and and are the binary semi-major axis and eccentricity, respectively. While it is clear that varies with from Equation 9, how is precisely affected by is not known since Kepler-34 and Kepler-16 differ in mass ratio as well as stellar binary eccentricity. The remaining question is do the gas and planetesimal disks adopt the same eccentricity; are and equal? On average we find at 1 au, for Kepler-16 (as per Figure 3e), whereas . Kepler-34 has an even larger difference between its value of and .

There is a noticeable trend in with aspect ratio but no clear change in the disk eccentricity or structure with varying -viscosity. We find considerably higher eccentricities in disks with a higher aspect ratio. This is explained by the increased level of communication in the disk fluid due to the larger value of the sound speed, . The higher allows waves in the disk to propagate over longer radial distances before being damped. This results in an increased pitch angle of the spiral arms and increases the radial extent to which the waves can travel. Waves which unfold further out can deposit their energy over a larger percentage of the disk which increases the mass-weighted mean disk eccentricity, (see Fig. 2). This effect is seen in Figure 5 as the enhanced propagation from a larger aspect ratio leads to a more significant contribution of modes in the outer regions of the disk. The deposition of energy as a result of these unfurling spirals in the large aspect ratio runs, is also indirectly observed through the enhanced eccentric modes in the outer disk, as the binary transfers energy to the fluid elements. Additionally, the radial eccentricity profiles of Kepler-16 show that a high aspect ratio leads to eccentricity increases in the outer regions of the disk. This can be seen in Figure 3e, where the eccentricity beyond 1.5 au in Kepler-16 is close to 0 for the low aspect ratio runs and , but remains greater than 0 for both runs and . This is also the case for Kepler-34 which has e = 0.05 at 3 au for both high aspect ratio runs and , but close to 0 for the low aspect ratio runs and .

There is a correlation between oscillations in the disk eccentricity and the rotational dynamics of the disk from cross examining the longitude of periastron in Figure 4 and the disk eccentricity in Figure 3. Pierens & Nelson (2013) suggest that the amplitude of the oscillations is linked to pressure effects. This is reflected in our results by the increasing amplitude of eccentricity oscillations with larger aspect ratio. For a disk around an eccentric binary Lubow & Artymowicz (2000) find the direct forcing from the binary potential () on the disk can become important and the growth of the disk eccentricity is described (their Equation 2) as

(10) |

where , and are the binary eccentricity, mass parameter () and semi-major axis respectively, is the Keplerian velocity of the disk and is the longitude of periastron of the disk relative to the binary. We do not retrieve the same magnitude of eccentricity growth as predicted by Equation 10 but this is likely due to the difference between the model’s dependency on ballistic particles in contrast to our fluid. The eccentricity growth and damping is well phased though with , mapping the eccentricity oscillations in Kepler-16. Since the binary is placed on a fixed orbit . Therefore by reading the evolution of the longitude of periastron of run B in Figure 4, during periods where , Equation 10 gives and hence eccentricity growth. These periods are matched by an eccentricity increase for run B in Figure 3c. Similarly regions where and hence , a damping of eccentricity is seen in the eccentricity plot. The comparison of the phasing of with for Kepler-34 is difficult since the eccentricity oscillations in disks around this binary are less pronounced. This validates the model further however, as for Kepler-34, since the stars have equal mass (). A similar effect is observed in the -body case where the mass ratio of Kepler-34 leads to the removal of secular forcing ( and dynamical forcing ) (Paardekooper et al., 2012).

There is a discrepancy between the centre-of-mass position angle showing perfect circulation of the disk and the longitude of periastron showing a combination of imperfect circulation and libration. This may be explained by the mass weighted value of the longitude of periastron becoming contaminated by outer regions of the disk where higher order disturbances can cause the value of to flip such that the orientation of the eccentric disk can vary between the inner and outer regions. If remained on the same hemisphere of the disk for all radii then would extend fully from - to + to show full circulation as per the centre-of-mass.

We find a very minimal contribution from self-gravity to the dynamics of the disks. At standard densities, in both Kepler-16 and Kepler-34 the eccentricity is marginally lower with gravity enabled, an effect observed by Marzari et al. (2009) in their study of circumstellar disks perturbed by an s-type binary (circumprimary) configuration. Additionally, the oscillation frequency in increases with gravity enabled due to the increase in precession rate of the disk as seen in Figure 4. Both these effects have only a minor impact on the disk’s density profile and its eccentricity. Therefore including disk self-gravity in simulations of disks at half minimum-mass solar nebula size is not necessary in modelling the disk and retrieving the correct eccentricity and structure.

We also find that self-gravity has no real influence on the evolution of the magnitude of the disk eccentricity or steady-state surface density at 2.5 and 5 times the minimum-mass solar nebula. At 12.5 times, self-gravity starts to become a necessary addition to the simulations to resolve the correct morphology of the disk. Enabling it at this high disk mass changes the structure of the disk, as can be seen in the surface density plot of Figure 12. A series of small density humps in the outer disk is observed that occur at regions of increased disk eccentricity. A noticeable difference in the disk evolution when self-gravity is included is the change in frequency of the eccentricity oscillations. The trend in final mass weighted mean disk eccentricity is to decrease with increasing surface density, with the exception of the highest density run. The disk eccentricity is probed further in Figure 12 where the eccentricity of the inner disk is greatly reduced with self-gravity enabled for 12.5 . The eccentricity then increases over the non self-gravitating run for 2.0 au. This may explain why the mass weighted value is larger at the end of the simulation, since it takes time to raise the eccentricity of the outer disk. Including self-gravity therefore acts to slightly damp the overall eccentricity of the disk. Global self-gravitating modes can exist in the disk if precession due to the companion star can be ignored (Papaloizou, 2002). However, it is likely that disk precession due to a massive binary companion cannot be ignored, in which case the global modes play no role of importance.

The Fourier analysis performed to identify differences in the strength of the low frequency modes did not reveal any significant trends with increasing surface density, with or without self-gravity included (Figure 9). With the exception of 12.5 where including self-gravity is required to accurately model the modes in the outer disk, it does not appear that including self-gravity is important when considering disks with 14. Regardless of the surface density magnitude, the results reveal the importance of including up to the contribution of the surface density in fully describing the structure of the disk since outside the cavity, within the disk, the strength of the modes is often 20% of the background axisymmetric surface density.

Our results indicate that the choice of boundary condition does not reflect strongly in the dynamical evolution of disks around both binary types. In both binary configurations, allowing material to flow through an open boundary increases the mass interior to the density peak. This effect is stronger in Kepler-34 than Kepler-16. Explanations for this include advanced wave damping routines operating in the rigid boundary case that could change the structure of the inner disk, or that the open boundary allows for a more realistic truncation timescale. There is little difference in the disk eccentricities beyond the truncation radius. These results are in contrast with Kley & Haghighipour (2014) who find much larger disk eccentricities when using a rigid boundary. This is probably to be expected since their shallower initial surface density profile () and omission of a gap function to model the expected inner cavity means that a large amount of disk mass is located just exterior to the inner boundary. This means a large disk mass is perturbed close to the binary, increasing disk eccentricity. In this case, an open boundary is quick to remove any highly eccentric material close to the boundary itself.

The final result from this work is the effect the initial surface density gradient, , has on the eccentricity. By changing from 1.5 to a shallower gradient of 0.5, we retrieve a much smaller eccentricity at all times. Using = 0.5 reduces the final disk eccentricity by a half. This result is not surprising since the calculation of the mean disk eccentricity is mass weighted and therefore biased towards areas of high mass. For = 0.5, more mass is in the outer disk, as seen in the surface density profiles of Figure 10. Since the outer disk is less perturbed by the binary and has a lower eccentricity, disks with a shallower density profile will have a lower mass weighted eccentricity. This explains why, alongside the varying binary parameters, Kley & Haghighipour (2014) find considerably smaller values for their disk eccentricities than found by us and PN13.

## 5 Summary and further work

In this paper we have performed a number of hydrodynamical simulations of circumbinary gas disks that explore the binary configuration, disk aspect ratio, -viscosity, inner boundary condition, initial surface density profile and inclusion of self-gravity on the evolution and quasi steady-state of such disks.

We found that the choice of stellar binary has a significant impact on the eccentricity of the disk with Kepler-34 pumping up the mean disk eccentricity to three times the value seen in some simulations of Kepler-16. The truncation radius is also strongly affected by the binary configuration, with the more eccentric Kepler-34 forcing the inner edge out to an average radius of 2.0 au, compared to 1.0 au in Kepler-16. There is no significant correlation between the -viscosity and the final surface density but there is a noticeable trend for the disk eccentricity to increase with increasing aspect-ratio. The eccentricity of the disk is strongly linked to the initial surface density gradient, which must therefore be taken into account when comparing work of a similar nature. At least for our initial surface density profile gradient value of = 1.5, the choice of inner boundary condition did not have a substantial effect on the evolution of the disk for both Kepler-34, and Kepler-16, but may play a role in resolving the disk interior to the density peak around the truncation radius.

We did not see a strong influence of self-gravity in disks near minimum-mass solar nebula size, and saw only a minimal decrease in mean disk eccentricity when self-gravity is included. Fourier analysis of the surface density to reveal the presence of modes also suggested that enabling self-gravity, even in disks where the surface density causes the Toomre parameter to approach 1, does not change the structure of the disks. We saw a subtle difference however in the mode strength and variation in our high mass (12.5 ) model when self-gravity is included, suggesting that work on disk masses larger than this value should include self-gravitation in order to accurately resolve the structure of the disk. This contrasts results showing the relevance of self-gravity in the circumprimary case since we did not find that it is a necessary inclusion in circumbinary disks that lie above the Toomre stability limit.

The next goal is to take surface density and velocity profiles of a steady-state disk for both Kepler-16 and Kepler-34 for use in gas potential feedback on our N-body simulations of planetesimal evolution and growth. From this work we found that many of our parameter choices have little influence on the final product, with the exception of aspect-ratio and initial surface density gradient which both have a significant impact on disk eccentricity. This will be an important aspect to consider in our future work since the eccentricity of the gas disk is a measure of the level of asymmetry in the gas disk potential.

Future work on the hydrodynamical evolution of circumbinary disks will include the thermal response of the gas by including an energy equation. Additionally, since the evolution of the disks is sensitive to binary configuration, further investigations should look at changing the binary parameter space: eccentricity, semi-major axis and mass ratio. This may help to identify the link between and .

###### Acknowledgements.

S.L and Z.M.L are supported by the STFC. P.J.C is grateful to NERC Grant NE/K004778/1. SJP 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 DICE which were both used to carry out this work. The authors would also like to thank the referee S. Lubow for his insightful report.## References

- Artymowicz & Lubow (1994) Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651
- Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
- Baruteau & Masset (2008) Baruteau, C., & Masset, F. 2008, ApJ, 678, 483
- Doyle et al. (2011) Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602
- Günther et al. (2004) Günther, R., Schäfer, C., & Kley, W. 2004, A&A, 423, 559
- Hayasaki & Okazaki (2009) Hayasaki, K., & Okazaki, A. T. 2009, ApJ, 691, L5
- Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35
- Heppenheimer (1978) Heppenheimer, T. A. 1978, AA, 65, 421
- Holman & Wiegert (1999) Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621
- Kley & Haghighipour (2014) Kley, W., & Haghighipour, N. 2014, A&A, 564, A72
- Kostov et al. (2014) Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14
- Lines et al. (2014) Lines, S., Leinhardt, Z. M., Paardekooper, S., Baruteau, C., & Thebault, P. 2014, ApJ, 782, L11
- Lubow (1991) Lubow, S. H. 1991, ApJ, 381, 259
- Lubow & Artymowicz (2000) Lubow, S. H., & Artymowicz, P. 2000, Protostars and Planets IV, 731
- Marzari et al. (2009) Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493
- Marzari et al. (2008) Marzari, F., Thébault, P., & Scholl, H. 2008, ApJ, 681, 1599
- Marzari et al. (2013) Marzari, F., Thebault, P., Scholl, H., Picogna, G., & Baruteau, C. 2013, A&A, 553, A71
- Masset (2000) Masset, F. 2000, A&AS, 141, 165
- Meschiari (2012a) Meschiari, S. 2012a, ApJ, 752, 71
- Meschiari (2012b) —. 2012b, ApJL, 761, L7
- Moriwaki & Nakagawa (2004) Moriwaki, K., & Nakagawa, Y. 2004, ApJ, 609, 1065
- Nelson (2003) Nelson, R. P. 2003, MNRAS, 345, 233
- Nixon & Lubow (2015) Nixon, C., & Lubow, S. 2015, MNRAS, 448, 3472
- Orosz et al. (2012a) Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, Science, 337, 1511
- Orosz et al. (2012b) —. 2012b, ApJ, 758, 87
- Paardekooper et al. (2012) Paardekooper, S.-J., Leinhardt, Z. M., Thébault, P., & Baruteau, C. 2012, ApJL, 754, L16
- Papaloizou (2002) Papaloizou, J. C. B., 2002, A&A, 388, 615
- Pelupessy & Portegies Zwart (2013) Pelupessy, F. I., & Portegies Zwart, S. 2013, MNRAS, 429, 895
- Pierens & Nelson (2007) Pierens, A., & Nelson, R. P. 2007, AA, 472, 993
- Pierens & Nelson (2008) —. 2008, AA, 483, 633
- Pierens & Nelson (2013) —. 2013, A&A, 556, A134
- Scholl et al. (2007) Scholl, H., Marzari, F., & Thébault, P. 2007, MNRAS, 380, 1119
- Schwamb et al. (2013) Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127
- Toomre (1964) Toomre, A. 1964, ApJ, 139, 1217
- Welsh et al. (2012) Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475
- Welsh et al. (2014) Welsh, W. F., Orosz, J. A., Short, D. R., et al. 2014, ArXiv e-prints, arXiv:1409.1605
- Wolszczan & Frail (1992) Wolszczan, A., & Frail, D. A. 1992, Nature, 355, 145