Modelling circumbinary protoplanetary disks
I. Fluid simulations of the Kepler16 and 34 systems
^{1} School of Physics, University of Bristol, H. H. Wills Physics Laboratory, Tyndall Avenue, Bristol, BS8 1TL, UK
email: stefan.lines@bristol.ac.uk; zoe.leinhardt@bristol.ac.uk; p.carter@bristol.ac.uk
^{2} CNRS, IRAP, 14 avenue Edouard Belin, 31400 Toulouse, France
^{3} Université de Toulouse, UPSOMP, IRAP, Toulouse, France
email: clement.baruteau@irap.omp.eu
^{4} Astronomy Unit, School of Physics Astronomy, Queen Mary University of London, UK
email: s.j.paardekooper@qmul.ac.uk
^{5} DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
Received: 10 April 2015
Accepted: 21 July 2015
Context. The Kepler mission’s discovery of a number of circumbinary planets orbiting close (a_{p}< 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: Kepler16 and Kepler34. We probe the importance of disk viscosity, aspectratio, inner boundary condition, initial surface density gradient, and selfgravity on the dynamical evolution of the disk, as well as its quasisteadystate profile.
Results. We find there is a strong influence of binary type on the mean disk eccentricity, e̅_{d}, leading to e̅_{d} = 0.02 − 0.08 for Kepler16 and e̅_{d} = 0.10 − 0.15 in Kepler34. The value of αviscosity has little influence on the disk, but we find a strong increase in mean disk eccentricity with increasing aspectratio 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 selfgravity has little impact on the evolution or final state of the disk for disks with masses less than 12.5 times that of the minimummass solar nebula. This finding contrasts with the results of selfgravity relevance in circumprimary disks, where its inclusion is found to be an important factor in describing the disk evolution.
Key words: hydrodynamics / methods: numerical / planets and satellites: formation / protoplanetary disks / binaries: close
© ESO, 2015
1. Introduction
Extrasolar planets in circumbinary (ptype) orbits around shortperiod 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 shortperiod planets^{1} with semimajor axis a_{p} < 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 longterm stability of planets in Ptype orbits. They calculated a range of orbital radii that allows for dynamical stability for a range of binary eccentricities and semimajor axes. Interestingly, with the exception of Kepler47(AB)c, all known circumbinary planets lie just outside their inner most stable orbit, a_{c}. 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 Jupitermass 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, Saturnmass planets have a stable evolution. A torque reversal, caused by an eccentricity increase from interaction with the binary, leads to stable outwards migration. However, Jupitermass 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 subJupiter mass, a result difficult to explain by observational biases since larger planets are easier to detect.
Fig. 1 Surface density maps, Σ_{2D}, for non selfgravitating models around Kepler16 (top row) and Kepler34 (bottom row) at 16 000 P_{AB}. Frames a)–d) and e)–h) correspond to simulations A–D and E–H in Table 1, respectively. 

Open with DEXTER 
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 kilometresized 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 Nbody 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 nonlinear mode coupling between the eccentric m = 1 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 m = 1 potential (Lubow & Artymowicz 2000). This coupling produces an m = 2 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 Fig. 1).
Parameter setup of each simulation A through U.
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 selfgravity on the evolution of the system is less well explored. Recent studies such as PN13 and Kley & Haghighipour (2014) choose not to include selfgravity in their simulations. The former justifies this approximation by confirming the Toomre parameter, (1)where c_{s} is the sound speed, Ω is the angular velocity, G is the gravitational constant, and Σ is the surface density, satisfies the Toomre stability criterion (Q> 1) at all disk radii for their chosen disk mass. However, even in selfgravitating disks of relatively low mass, lowfrequency global modes exist that do not have a counterpart in non selfgravitating disks (Papaloizou 2002). Moreover, the study by Marzari et al. (2009) shows that selfgravity in circumprimary disks can be important even in lowmass disks. Therefore, an investigation into whether or not selfgravity 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 Kepler16 and Kepler34 systems. The resulting quasisteadystate surface density profiles will be used to parameterize the gas disk potential for use in our Nbody 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 Sects. 3 and 4 we discuss the results and implications of the 21 simulations of the Kepler16 and Kepler34 systems. Finally, we summarise our findings in the conclusion in Sect. 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 q = M_{B}/M_{A} are first initialised at their respective apoapses, such that the Cartesian coordinates are initially zero in the Y direction for both stars, Y_{A}&Y_{B} = 0, and the X position is (2)and (3)where M_{A} and M_{B} are the mass of star A and B, respectively, and the distance of a star from its focus, r, is determined from the shape equation (4)where a_{b} is the binary semimajor axis, e_{b} 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 θ(t) and thus, the position vectors of each star. The stellar binary parameters are given in Table 2.
Fig. 2 Radial velocity maps at 16 000 P_{AB} for a) Kepler16 α = 10^{3}, h = 0.03; b) Kepler16 α = 10^{3}, h = 0.05; c) Kepler34 α = 10^{3}, h = 0.03 and d) Kepler34 α = 10^{3}, h = 0.05. 

Open with DEXTER 
Stellar binary parameters.
To simulate the hydrodynamical evolution of the circumbinary disks we use a modified version of the Eulerian fluid code FARGO (Masset 2000), called FARGOADSG (Baruteau & Masset 2008), which includes optional disk selfgravity. The hydrodynamical equations are solved on a polar mesh which is composed of n_{φ} = 512 equally divided azimuthal sector cells, and n_{r} = 395 logarithmically spaced radial cells. Logarithmic spacing in the radial direction is a necessary condition for the fast fourier transform algorithm in the selfgravity 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 zerogradient 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, f_{gap} from Günther et al. (2004) to better model the initial surface density profile in the inner disk region: (5)where R_{gap} = 2.5a_{b} is the estimated size of the gap (Artymowicz & Lubow 1994) and R 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.
Fig. 3 Timeaveraged quasisteady state surface density profiles (a) and b)), Σ, and time evolution of disk mean eccentricity (c) and d)), , and timeaveraged radial eccentricity profiles e) and f), e), for Kepler16 (Models A–D) and Kepler34 (Models E–H) under non selfgravitating conditions. Timeaveraging is over the last 1000 binary orbits of the simulations. 

Open with DEXTER 
The initial surface density of the gas disk follows Σ(R) = f_{gap}Σ_{0}R^{− αΣ} from PN13, where Σ_{0} is the surface density at 1 au, R is the radial position in au, and α_{Σ} defines the initial surface density gradient. For the majority of the simulations (Table 1 simulations AN and U) we assume a half minimummass solar nebula, thus, Σ_{0} = 10^{4}M_{⊙}/au^{2} (Hayashi 1981). In simulations O − TΣ_{0} 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, U, 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 Kepler34 and Kepler16. 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, Σ_{min}, is set to 10^{9} of the initial density. This means that mass is not strictly conserved in regions of very low density.
The aspect ratio, , where H 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, c_{s} = v_{k}h, where v_{k} 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^{3} and 10^{2}.
Since our disks are low in mass, we expect them all to be linearly stable according to the Toomre stability criterion, Q> 1 (Toomre 1964). For our halfMMSN density models, even near the inner edge of the disk where the density is at a maximum, the minimum Toomre value, Q_{min}, is 140. However, on consideration of the results of Marzari et al. (2009) that show selfgravity plays a significant role in circumprimary disks, selfgravity is included in simulations I–L and O–Q.
Over 21 simulations, we explore the effect that changing the stellar binary parameters, aspect ratio, αviscosity, boundary condition and inclusion of disk selfgravity has on the evolution and quasi steadystate (QSS) profile of the disk. All simulations are evolved for 16 000 binary orbits, P_{AB}, equivalent to 1700 orbits at 1 au, with the exception of M and N which run for 4000 P_{AB} to test the effect of the chosen inner boundary condition on the initial response of the disk.
3. Results
3.1. Non selfgravitating disks
3.1.1. Kepler16
We first look at the response of a non selfgravitating disk to a stellar binary with parameters chosen to match the observed properties of the system Kepler16(AB). The azimuthally and time averaged quasisteady state surface density is shown in Fig. 3a. The density maximum or peak, Σ_{peak}, indicates the location of the truncated inner edge, although the eccentric morphology of the disk requires a full 2dimensional density map, Σ_{2D} (Fig. 1), to identify the edge precisely as a function of azimuth angle. Therefore, we will consider the location of Σ_{peak} 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 (Figs. 1c and e), particularly for α = 10^{3}, 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 A 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 (Fig. 3c), , which we define according to Pierens & Nelson (2007) as (6)where R_{in} and R_{out} are the disk inner and outer radii and Σ_{c} and e_{c} are the fluid element surface density and eccentricity respectively. Note e_{c} is calculated assuming the cell is orbiting the binary barycentre alone (a two body problem). Indeed, models with h = 0.05 show a higher value of at times both early in evolution and quasisteady 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 7.5 × 10^{3} to 1.3 × 10^{2} as opposed to a consistent value of 5.0 × 10^{3} 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 B and D) have a larger average eccentricity than those with a smaller aspect ratio (models A and C) 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 Fig. 4, is plot as a function of time for Kepler16. 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 Sect. 4. The position angle of the disk centreofmass is calculated in the inertial frame and shows the circulation of the disk increasing in frequency until the disk precession period, P_{d}, settles at around 2500 P_{AB} at QSS.
The eccentric morphology of the disk is apparent through both Σ_{2D} (Fig. 1) and . Figure 2 shows the radial velocity which reveals a number of spiral features. In both the h = 0.03 and 0.05 models of Kepler16, the velocity map shows m = 2 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.
Fig. 4 Time evolution of disk mean longitude of periastron and the position angle of the disk centreofmass for Kepler16 under both non selfgravitating conditions (black solid) and selfgravitating (red dashed) conditions. 

Open with DEXTER 
The disk surface density can be decomposed into modes described by azimuthal mode numbers m and frequency mode numbers l 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 Σ_{l,m}(r) is a complex function and Ω_{b} is the mean motion of the binary, 2π/P_{AB}. Nixon & Lubow (2015) find that eccentric binaries produce modes with l ≠ m. In particular, in their retrograde circumbinary disk simulations, highly eccentric binaries (e_{b} ≥ 0.6) produce powerful l = − 1, m = 2 modes. We remove the time dependence of the disk disturbances by since our analysis is concerned with the eccentric Kepler16 system (e_{b} = 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 l = 0. 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 Fig. 5 the strengths of the m = 1 and m = 2 modes relative to the axisymmetric component is plot as a function of orbital radius. The results confirm that for the high aspect ratio (h = 0.05) models B and D, the contribution from the m = 2 spirals is more significant in the outer disk than that seen for h = 0.03.
Fig. 5 Fourier analysis of the surface density of all non selfgravitating Kepler16 circumbinary disks. The transform is done on the timeaveraged surface density, over a single binary orbit at 15 000 P_{AB}, to consider only l = 0. The strengths of both the m = 1 and m = 2 modes are normalised against the axisymmetric (m = 0) contribution and plot as a function of radius in the disk. 

Open with DEXTER 
3.1.2. Kepler34
Disks surrounding Kepler34 are significantly more affected by the binary, with the QSS ranging from 0.1 to 0.15, typically three times more eccentric than Kepler16 (Fig. 3d). Since the massweighted 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 P_{AB}. However, looking at the instantaneous surface density profiles for a Kepler34 run during the last few hundred binary orbits as shown in Fig. 6, it appears as though a steadystate has been reached. Our Kepler34 results agree roughly with PN13, who find that typically, for the systems mutually covered, a quasi steadystate is achieved by 10^{4}P_{AB}.
Fig. 6 Instantaneous surface density profiles for Model H (Kepler34) at three times near quasisteadystate separated by 640 P_{AB}. 

Open with DEXTER 
As with Kepler16 the average eccentricity in the Kepler34 disks once they have reached QSS is larger for higher values of aspect ratio. The large difference in ( = 0.05) between models F and G is shown in Fig. 3f to originate from the disk beyond the truncation radius where the surface density is much higher, and is therefore expected since is a massweighted value. From 2.0 au, model F 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. Selfgravitating disks
In models I − L we enable disk selfgravity (DSG) for both Kepler16 and Kepler34. The αviscosity is fixed at 10^{3} but the aspect ratio is varied between 0.03 and 0.05 to adjust the Toomre value and allow for direct comparison with models A, B, E and F. For Kepler16, as seen in Fig. 7, during the first 2000 binary orbits the disk responds almost identically to that of its non selfgravitating 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 Kepler34 simulations, it is possible only to differentiate between self and non selfgravitating 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 Fig. 7 the time averaged surface density profiles also reveal that a selfgravitating disk in Kepler16 has no effect on the final disk structure, and almost no effect for disks around Kepler34.
To test at what disk mass selfgravity become important, Σ_{0} is increased by 5, 10 and 25 times. Selfgravity importance is therefore tested for Q→ 1 but Q> 1 at all times. We define standard density models (Σ_{0} = 1 × 10^{4} code units) as being 0.5 Σ_{MMSN} and those scaled up by a factor of 5, 10 and 25 as 2.5 Σ_{MMSN} (Q_{min} = 28), 5 Σ_{MMSN} (Q_{min} = 14) and 12.5 Σ_{MMSN} (Q_{min} = 5) respectively.
In Fig. 8 the evolution of the mean disk eccentricity for each of the four varying densities is shown. remains, for the non selfgravitating runs, at a constant 0.08 for each surface density. When selfgravity is enabled, with the exception of 12.5 Σ_{MMSN}, 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 selfgravitating 12.5 Σ_{MMSN} simulation is more dynamic due to an eccentricity fall off before increasing again.
Fig. 7 Comparison of Kepler16 disk eccentricity evolution and time averaged surface density for 0.5 Σ_{MMSN} density disks with (I, J, K and L) and without selfgravity (A, B, E and F). Selfgravitating runs are shown with dashed lines. 

Open with DEXTER 
At standard densities we find that the oscillation frequency increases and amplitude decreases with selfgravity enabled. Increasing the density does not change the oscillation frequency without selfgravity enabled, but the period of the oscillations decreases from 2000 P_{AB} (0.5 Σ_{MMSN}) to 300 P_{AB} (12.5 Σ_{MMSN}) when it is included.
Another measure of the effect of selfgravity 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 selfgravitating Kepler16 runs, is done. The Fourier coefficients are determined and normalised against the axisymmetric (m = 0) component. The strength of these modes normalised to the axisymmetric part, Φ_{m}/ Φ_{0}, is shown for each surface density in Fig. 9. Again, a comparison is made between simulations with and without selfgravity included.
There is a strong contribution from both the m = 1 and m = 2 modes, with the m = 1 mode always 20% or larger in the inner disk. It is less clear as to how the inclusion of selfgravity changes the strength of these modes. When selfgravity is enabled in the 0.5 Σ_{MMSN}, 2.5 Σ_{MMSN} and 5 Σ_{MMSN} surface density runs, there are only subtle differences between the strength of the modes. At 12.5 Σ_{MMSN}, due to oscillations in the strength of both modes in the selfgravitating disk, the profiles for both the selfgravitating m = 1 and m = 2 modes do not trace their non selfgravitating equivalents.
Fig. 8 Evolution of the mean disk eccentricity for standard density selfgravitating Kepler16 run J and increased density runs O–T (dashed lines). The corresponding simulation with selfgravity disabled is shown by solid lines. 

Open with DEXTER 
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 U the density gradient is made shallower by changing the surface density exponent α_{Σ} from 1.5 to 0.5. In Fig. 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 Σ_{peak} 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 v_{r} and v_{φ}. 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.
Fig. 9 Fourier analysis of the surface density of Kepler16 simulations testing the relevance of selfgravity with increasing Σ_{0}. Black lines correspond to the Fourier coefficients of the m = 1 modes, normalised to the axisymmetric m = 0 component. Red lines show the m = 2 modes. Simulations where selfgravity is enabled are shown by dotted lines. Fourier analysis is done on the timeaveraged surface density over the period of one binary orbit from 15 000 P_{AB}. 

Open with DEXTER 
We find that there are small differences in both the disk eccentricity and structure between the two. In Kepler16 the surface density profiles show that including an open boundary keeps the peak density at the same location but increases the mass interior to Σ_{peak}. This effect, which is shown in Fig. 11 is stronger in Kepler34 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 Sect. 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 Kepler16 and Kepler34 stellar binary systems. We vary αviscosity and aspect ratio, as well as probing the effects of disk selfgravity, initial surface density profile, and the inner boundary condition.
Fig. 10 Comparison of surface density and eccentricity evolution in Kepler16 runs between density profiles of form Σ(r)∝Σ_{0}r^{− αΣ} with α_{Σ} = 0.5 (blue dashed) and 1.5 (black solid). 

Open with DEXTER 
Fig. 11 Comparison of open (dashed) and closed (solid) boundary conditions in Kepler16 (black) and Kepler34 (blue) via surface density and eccentricity profiles. 

Open with DEXTER 
The majority of our simulations have reached steady state (a consistent value of ) by 16 000 P_{AB}. The steadystate 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 (Fig. 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 Σ_{peak} in Fig. 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 R_{gap} = 2.5a_{b}.
The results for the surface densities and eccentricities of our non selfgravitating disks agree with those found by PN13. Subtle differences such as the eccentricity magnitude (Fig. 3), particularly in Kepler34, 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 Kepler34, and leads to lower disk eccentricities for a< 1.5 au in Kepler16. 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 semimajor axis. Therefore, it is likely that a circumbinary gas disk is only weakly affected by the relatively low eccentricity of Kepler38 (e = 0.10) (Orosz et al. 2012b) compared to the significantly higher eccentricity of Kepler34 (e = 0.53). We would therefore expect the gas disk around Kepler34 to adopt a higher eccentricity than a similar disk around Kepler38. This theory is supported by the significantly larger disk eccentricities in simulations around the eccentric Kepler34 compared to Kepler16 where the stellar binary has a much smaller eccentricity.
Increasing disk eccentricity with increasing stellar binary eccentricity is a trend also seen in Nbody results. The planetesimal response in the Nbody scenario is dependent upon the forcing of stellar binary on the disk, and perturbation theory provides a value of eccentricity of which planetesimals tend towards. This forced eccentricity, e_{f}, is given by (Moriwaki & Nakagawa 2004)(9)where M_{A} and M_{B} are the primary and secondary stellar masses, M_{∗} is the total binary mass, R is the distance from the binary barycentre and a_{b} and e_{b} are the binary semimajor axis and eccentricity, respectively. While it is clear that e_{f} varies with e_{b} from Eq. (9), how e_{d} is precisely affected by e_{b} is not known since Kepler34 and Kepler16 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 e_{f} and e_{d} equal? On average we find at 1 au, e_{d} ~ 0.2 for Kepler16 (as per Fig. 3e), whereas e_{f} ~ 0.02. Kepler34 has an even larger difference between its value of e_{d} ~ 0.45 and e_{f} ~ 0.002.
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, c_{s}. The higher c_{s} 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 massweighted mean disk eccentricity, (see Fig. 2). This effect is seen in Fig. 5 as the enhanced propagation from a larger aspect ratio leads to a more significant contribution of m = 2 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 m = 1 modes in the outer disk, as the binary transfers energy to the fluid elements. Additionally, the radial eccentricity profiles of Kepler16 show that a high aspect ratio leads to eccentricity increases in the outer regions of the disk. This can be seen in Fig. 3e, where the eccentricity beyond 1.5 au in Kepler16 is close to 0 for the low aspect ratio runs A and C, but remains greater than 0 for both runs B and D. This is also the case for Kepler34 which has e = 0.05 at 3 au for both high aspect ratio runs F and H, but close to 0 for the low aspect ratio runs E and G.
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 Fig. 4 and the disk eccentricity in Fig. 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 m = 1 binary potential (Φ_{1,1}) on the disk can become important and the growth of the disk eccentricity is described (their Eq. (2)) as (10)where e_{b}, μ_{b} and a_{b} are the binary eccentricity, mass parameter (M_{B}/M_{t}) and semimajor axis respectively, Ω_{d} 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 Eq. (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 Kepler16. Since the binary is placed on a fixed orbit . Therefore by reading the evolution of the longitude of periastron of run B in Fig. 4, during periods where ϖ< 0, Eq. (10) gives ė_{d}> 0 and hence eccentricity growth. These periods are matched by an eccentricity increase for run B in Fig. 3c. Similarly regions where ϖ> 0 and hence ė_{d}< 0, a damping of eccentricity is seen in the eccentricity plot. The comparison of the phasing of ė_{d} with ϖ for Kepler34 is difficult since the eccentricity oscillations in disks around this binary are less pronounced. This validates the model further however, as Φ_{1,1} = 0 for Kepler34, since the stars have equal mass (μ_{b} = 0.5). A similar effect is observed in the Nbody case where the mass ratio of Kepler34 leads to the removal of secular forcing (e_{f} = 0 and dynamical forcing e_{ff} ≠ 0) (Paardekooper et al. 2012).
There is a discrepancy between the centreofmass 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 centreofmass.
We find a very minimal contribution from selfgravity to the dynamics of the disks. At standard densities, in both Kepler16 and Kepler34 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 stype 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 Fig. 4. Both these effects have only a minor impact on the disk’s density profile and its eccentricity. Therefore including disk selfgravity in simulations of disks at half minimummass solar nebula size is not necessary in modelling the disk and retrieving the correct eccentricity and structure.
Fig. 12 Comparison of surface density and eccentricity profiles in Kepler16 runs between selfgravitating (dashed) and non selfgravitating (solid) for 0.5 Σ_{MMSN} (black) and 12.5 Σ_{MMSN} (red). The 12.5 Σ_{MMSN} surface density profile is normalised to the initial half minimummass solar nebula surface density. 

Open with DEXTER 
We also find that selfgravity has no real influence on the evolution of the magnitude of the disk eccentricity or steadystate surface density at 2.5 and 5 times the minimummass solar nebula. At 12.5 times, selfgravity 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 Fig. 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 selfgravity 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 Fig. 12 where the eccentricity of the inner disk is greatly reduced with selfgravity enabled for 12.5 Σ_{MMSN}. The eccentricity then increases over the non selfgravitating run for R> 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 selfgravity therefore acts to slightly damp the overall eccentricity of the disk. Global selfgravitating 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 selfgravity included (Fig. 9). With the exception of 12.5 Σ_{MMSN} where including selfgravity is required to accurately model the m = 1 modes in the outer disk, it does not appear that including selfgravity is important when considering disks with Q> 14. Regardless of the surface density magnitude, the results reveal the importance of including up to the m = 2 contribution of the surface density in fully describing the structure of the disk since outside the cavity, within the disk, the strength of the m = 2 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 Kepler34 than Kepler16. 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 (Σ ∝ r^{− 1 / 2}) 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 Fig. 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 selfgravity on the evolution and quasi steadystate of such disks.
We found that the choice of stellar binary has a significant impact on the eccentricity of the disk with Kepler34 pumping up the mean disk eccentricity to three times the value seen in some simulations of Kepler16. The truncation radius is also strongly affected by the binary configuration, with the more eccentric Kepler34 forcing the inner edge out to an average radius of 2.0 au, compared to 1.0 au in Kepler16. 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 aspectratio. 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 Kepler34, and Kepler16, 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 selfgravity in disks near minimummass solar nebula size, and saw only a minimal decrease in mean disk eccentricity when selfgravity is included. Fourier analysis of the surface density to reveal the presence of modes also suggested that enabling selfgravity, even in disks where the surface density causes the Toomre parameter Q to approach 1, does not change the structure of the disks. We saw a subtle difference however in the m = 1 mode strength and variation in our high mass (12.5 Σ_{MMSN}) model when selfgravity is included, suggesting that work on disk masses larger than this value should include selfgravitation in order to accurately resolve the structure of the disk. This contrasts results showing the relevance of selfgravity 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 steadystate disk for both Kepler16 and Kepler34 for use in gas potential feedback on our Nbody 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 aspectratio 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,
semimajor axis and mass ratio. This may help to identify the link between and e_{f}.
Kepler16(AB)b (Doyle et al. 2011), Kepler34(AB)b 35(AB)b (Welsh et al. 2012), Kepler38(AB)b (Orosz et al. 2012b), Kepler47(AB)b,c (Orosz et al. 2012a), Kepler64(AB)b (Schwamb et al. 2013), Kepler413(AB)b (Kostov et al. 2014) and KIC 9632895 (Welsh et al. 2015).
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 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, P., & Lubow, S. H. 1994, ApJ, 421, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Günther, R., Schäfer, C., & Kley, W. 2004, A&A, 423, 559 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayasaki, K., & Okazaki, A. T. 2009, ApJ, 691, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Heppenheimer, T. A. 1978, A&A, 65, 421 [NASA ADS] [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., & Haghighipour, N. 2014, A&A, 564, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14 [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]
 Lubow, S. H. 1991, ApJ, 381, 259 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H., & Artymowicz, P. 2000, Protostars and Planets IV, 731 [Google Scholar]
 Marzari, F., Thébault, P., & Scholl, H. 2008, ApJ, 681, 1599 [NASA ADS] [CrossRef] [Google Scholar]
 Marzari, F., Scholl, H., Thébault, P., & Baruteau, C. 2009, A&A, 508, 1493 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marzari, F., Thebault, P., Scholl, H., Picogna, G., & Baruteau, C. 2013, A&A, 553, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [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]
 Nelson, R. P. 2003, MNRAS, 345, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Nixon, C., & Lubow, S. 2015, MNRAS, 448, 3472 [NASA ADS] [CrossRef] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, Science, 337, 1511 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, ApJ, 758, 87 [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]
 Papaloizou, J. C. B. 2002, A&A, 388, 615 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pelupessy, F. I., & Portegies Zwart, S. 2013, MNRAS, 429, 895 [NASA ADS] [CrossRef] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2007, A&A, 472, 993 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2008, A&A, 483, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2013, A&A, 556, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scholl, H., Marzari, F., & Thébault, P. 2007, MNRAS, 380, 1119 [NASA ADS] [CrossRef] [Google Scholar]
 Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [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]
 Welsh, W. F., Orosz, J. A., Short, D. R., et al. 2015, ApJ, 809, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Wolszczan, A., & Frail, D. A. 1992, Nature, 355, 145 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Surface density maps, Σ_{2D}, for non selfgravitating models around Kepler16 (top row) and Kepler34 (bottom row) at 16 000 P_{AB}. Frames a)–d) and e)–h) correspond to simulations A–D and E–H in Table 1, respectively. 

Open with DEXTER  
In the text 
Fig. 2 Radial velocity maps at 16 000 P_{AB} for a) Kepler16 α = 10^{3}, h = 0.03; b) Kepler16 α = 10^{3}, h = 0.05; c) Kepler34 α = 10^{3}, h = 0.03 and d) Kepler34 α = 10^{3}, h = 0.05. 

Open with DEXTER  
In the text 
Fig. 3 Timeaveraged quasisteady state surface density profiles (a) and b)), Σ, and time evolution of disk mean eccentricity (c) and d)), , and timeaveraged radial eccentricity profiles e) and f), e), for Kepler16 (Models A–D) and Kepler34 (Models E–H) under non selfgravitating conditions. Timeaveraging is over the last 1000 binary orbits of the simulations. 

Open with DEXTER  
In the text 
Fig. 4 Time evolution of disk mean longitude of periastron and the position angle of the disk centreofmass for Kepler16 under both non selfgravitating conditions (black solid) and selfgravitating (red dashed) conditions. 

Open with DEXTER  
In the text 
Fig. 5 Fourier analysis of the surface density of all non selfgravitating Kepler16 circumbinary disks. The transform is done on the timeaveraged surface density, over a single binary orbit at 15 000 P_{AB}, to consider only l = 0. The strengths of both the m = 1 and m = 2 modes are normalised against the axisymmetric (m = 0) contribution and plot as a function of radius in the disk. 

Open with DEXTER  
In the text 
Fig. 6 Instantaneous surface density profiles for Model H (Kepler34) at three times near quasisteadystate separated by 640 P_{AB}. 

Open with DEXTER  
In the text 
Fig. 7 Comparison of Kepler16 disk eccentricity evolution and time averaged surface density for 0.5 Σ_{MMSN} density disks with (I, J, K and L) and without selfgravity (A, B, E and F). Selfgravitating runs are shown with dashed lines. 

Open with DEXTER  
In the text 
Fig. 8 Evolution of the mean disk eccentricity for standard density selfgravitating Kepler16 run J and increased density runs O–T (dashed lines). The corresponding simulation with selfgravity disabled is shown by solid lines. 

Open with DEXTER  
In the text 
Fig. 9 Fourier analysis of the surface density of Kepler16 simulations testing the relevance of selfgravity with increasing Σ_{0}. Black lines correspond to the Fourier coefficients of the m = 1 modes, normalised to the axisymmetric m = 0 component. Red lines show the m = 2 modes. Simulations where selfgravity is enabled are shown by dotted lines. Fourier analysis is done on the timeaveraged surface density over the period of one binary orbit from 15 000 P_{AB}. 

Open with DEXTER  
In the text 
Fig. 10 Comparison of surface density and eccentricity evolution in Kepler16 runs between density profiles of form Σ(r)∝Σ_{0}r^{− αΣ} with α_{Σ} = 0.5 (blue dashed) and 1.5 (black solid). 

Open with DEXTER  
In the text 
Fig. 11 Comparison of open (dashed) and closed (solid) boundary conditions in Kepler16 (black) and Kepler34 (blue) via surface density and eccentricity profiles. 

Open with DEXTER  
In the text 
Fig. 12 Comparison of surface density and eccentricity profiles in Kepler16 runs between selfgravitating (dashed) and non selfgravitating (solid) for 0.5 Σ_{MMSN} (black) and 12.5 Σ_{MMSN} (red). The 12.5 Σ_{MMSN} surface density profile is normalised to the initial half minimummass solar nebula surface density. 

Open with DEXTER  
In the text 