EDP Sciences
Free Access
Issue
A&A
Volume 582, October 2015
Article Number A5
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201526295
Published online 25 September 2015

© ESO, 2015

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 planets1 with semi-major axis ap < 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, ac. 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.

thumbnail Fig. 1

Surface density maps, Σ2D, for non self-gravitating models around Kepler-16 (top row) and Kepler-34 (bottom row) at 16 000 PAB. Frames a)d) and e)h) correspond to simulations AD and EH 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 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 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).

Table 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 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 cs 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 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 N-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 Sects. 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 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 = MB/MA are first initialised at their respective apoapses, such that the Cartesian coordinates are initially zero in the Y direction for both stars, YA&YB = 0, and the X position is (2)and (3)where MA and MB 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 ab is the binary semi-major axis, eb 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.

thumbnail Fig. 2

Radial velocity maps at 16 000 PAB for a) Kepler-16 α = 10-3, h = 0.03; b) Kepler-16 α = 10-3, h = 0.05; c) Kepler-34 α = 10-3, h = 0.03 and d) Kepler-34 α = 10-3, h = 0.05.

Open with DEXTER

Table 2

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 FARGO-ADSG (Baruteau & Masset 2008), which includes optional disk self-gravity. The hydrodynamical equations are solved on a polar mesh which is composed of nφ = 512 equally divided azimuthal sector cells, and nr = 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, fgap from Günther et al. (2004) to better model the initial surface density profile in the inner disk region: (5)where Rgap = 2.5ab 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.

thumbnail Fig. 3

Time-averaged quasi-steady state surface density profiles (a) and b)), Σ, and time evolution of disk mean eccentricity (c) and d)), , and time-averaged radial eccentricity profiles e) and f), e), for Kepler-16 (Models A–D) and Kepler-34 (Models E–H) under non self-gravitating conditions. Time-averaging is over the last 1000 binary orbits of the simulations.

Open with DEXTER

The initial surface density of the gas disk follows Σ(R) = fgapΣ0RαΣ 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 A-N and U) we assume a half minimum-mass solar nebula, thus, Σ0 = 10-4M/au2 (Hayashi 1981). In simulations OTΣ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 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, Σ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, cs = vkh, where vk 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 half-MMSN density models, even near the inner edge of the disk where the density is at a maximum, the minimum Toomre value, Qmin, 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 IL and OQ.

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, PAB, equivalent to 1700 orbits at 1 au, with the exception of M and N which run for 4000 PAB 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 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 2-dimensional 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 Rin and Rout are the disk inner and outer radii and Σc and ec are the fluid element surface density and eccentricity respectively. Note ec 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 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 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 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 Sect. 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, Pd, settles at around 2500 PAB 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 Kepler-16, 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.

thumbnail Fig. 4

Time evolution of disk mean longitude of periastron and the position angle of the disk centre-of-mass for Kepler-16 under both non self-gravitating conditions (black solid) and self-gravitating (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π/PAB. Nixon & Lubow (2015) find that eccentric binaries produce modes with lm. In particular, in their retrograde circumbinary disk simulations, highly eccentric binaries (eb ≥ 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 Kepler-16 system (eb = 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.

thumbnail Fig. 5

Fourier analysis of the surface density of all non self-gravitating Kepler-16 circumbinary disks. The transform is done on the time-averaged surface density, over a single binary orbit at 15 000 PAB, 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. 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 (Fig. 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 PAB. However, looking at the instantaneous surface density profiles for a Kepler-34 run during the last few hundred binary orbits as shown in Fig. 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 104PAB.

thumbnail Fig. 6

Instantaneous surface density profiles for Model H (Kepler-34) at three times near quasi-steady-state separated by 640 PAB.

Open with DEXTER

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 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 mass-weighted 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. Self-gravitating disks

In models IL we enable disk self-gravity (DSG) for both Kepler-16 and Kepler-34. 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 Kepler-16, as seen in Fig. 7, during the first 2000 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 Fig. 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, Σ0 is increased by 5, 10 and 25 times. Self-gravity 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 (Qmin = 28), 5 ΣMMSN (Qmin = 14) and 12.5 ΣMMSN (Qmin = 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 self-gravitating runs, at a constant 0.08 for each surface density. When self-gravity 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 self-gravitating 12.5 ΣMMSN simulation is more dynamic due to an eccentricity fall off before increasing again.

thumbnail Fig. 7

Comparison of Kepler-16 disk eccentricity evolution and time averaged surface density for 0.5 ΣMMSN density disks with (I, J, K and L) and without self-gravity (A, B, E and F). Self-gravitating runs are shown with dashed lines.

Open with DEXTER

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 PAB (0.5 ΣMMSN) to 300 PAB (12.5 ΣMMSN) 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 (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 self-gravity 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 self-gravity changes the strength of these modes. When self-gravity 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 self-gravitating disk, the profiles for both the self-gravitating m = 1 and m = 2 modes do not trace their non self-gravitating equivalents.

thumbnail Fig. 8

Evolution of the mean disk eccentricity for standard density self-gravitating Kepler-16 run J and increased density runs OT (dashed lines). The corresponding simulation with self-gravity 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 vr 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.

thumbnail Fig. 9

Fourier analysis of the surface density of Kepler-16 simulations testing the relevance of self-gravity 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 self-gravity is enabled are shown by dotted lines. Fourier analysis is done on the time-averaged surface density over the period of one binary orbit from 15 000 PAB.

Open with DEXTER

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 Σpeak. This effect, which is shown in Fig. 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 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 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.

thumbnail Fig. 10

Comparison of surface density and eccentricity evolution in Kepler-16 runs between density profiles of form Σ(r)Σ0rαΣ with αΣ = 0.5 (blue dashed) and 1.5 (black solid).

Open with DEXTER

thumbnail Fig. 11

Comparison of open (dashed) and closed (solid) boundary conditions in Kepler-16 (black) and Kepler-34 (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 PAB. 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 (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 Rgap = 2.5ab.

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 (Fig. 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 a< 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 N-body results. The planetesimal response in the N-body 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, ef, is given by (Moriwaki & Nakagawa 2004)(9)where MA and MB are the primary and secondary stellar masses, M is the total binary mass, R is the distance from the binary barycentre and ab and eb are the binary semi-major axis and eccentricity, respectively. While it is clear that ef varies with eb from Eq. (9), how ed is precisely affected by eb 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 ef and ed equal? On average we find at 1 au, ed ~ 0.2 for Kepler-16 (as per Fig. 3e), whereas ef ~ 0.02. Kepler-34 has an even larger difference between its value of ed ~ 0.45 and ef ~ 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, cs. The higher cs 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 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 Kepler-16 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 Kepler-16 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 Kepler-34 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 eb, μb and ab are the binary eccentricity, mass parameter (MB/Mt) and semi-major 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 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 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 Kepler-34 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 Kepler-34, since the stars have equal mass (μb = 0.5). A similar effect is observed in the N-body case where the mass ratio of Kepler-34 leads to the removal of secular forcing (ef = 0 and dynamical forcing eff ≠ 0) (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 Fig. 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.

thumbnail Fig. 12

Comparison of surface density and eccentricity profiles in Kepler-16 runs between self-gravitating (dashed) and non self-gravitating (solid) for 0.5 ΣMMSN (black) and 12.5 ΣMMSN (red). The 12.5 ΣMMSN surface density profile is normalised to the initial half minimum-mass solar nebula surface density.

Open with DEXTER

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 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 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 Fig. 12 where the eccentricity of the inner disk is greatly reduced with self-gravity enabled for 12.5 ΣMMSN. The eccentricity then increases over the non self-gravitating 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 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 (Fig. 9). With the exception of 12.5 ΣMMSN where including self-gravity is required to accurately model the m = 1 modes in the outer disk, it does not appear that including self-gravity 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 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 (Σ ∝ 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 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 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 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 ef.


1

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

All Tables

Table 1

Parameter setup of each simulation A through U.

Table 2

Stellar binary parameters.

All Figures

thumbnail Fig. 1

Surface density maps, Σ2D, for non self-gravitating models around Kepler-16 (top row) and Kepler-34 (bottom row) at 16 000 PAB. Frames a)d) and e)h) correspond to simulations AD and EH in Table 1, respectively.

Open with DEXTER
In the text
thumbnail Fig. 2

Radial velocity maps at 16 000 PAB for a) Kepler-16 α = 10-3, h = 0.03; b) Kepler-16 α = 10-3, h = 0.05; c) Kepler-34 α = 10-3, h = 0.03 and d) Kepler-34 α = 10-3, h = 0.05.

Open with DEXTER
In the text
thumbnail Fig. 3

Time-averaged quasi-steady state surface density profiles (a) and b)), Σ, and time evolution of disk mean eccentricity (c) and d)), , and time-averaged radial eccentricity profiles e) and f), e), for Kepler-16 (Models A–D) and Kepler-34 (Models E–H) under non self-gravitating conditions. Time-averaging is over the last 1000 binary orbits of the simulations.

Open with DEXTER
In the text
thumbnail Fig. 4

Time evolution of disk mean longitude of periastron and the position angle of the disk centre-of-mass for Kepler-16 under both non self-gravitating conditions (black solid) and self-gravitating (red dashed) conditions.

Open with DEXTER
In the text
thumbnail Fig. 5

Fourier analysis of the surface density of all non self-gravitating Kepler-16 circumbinary disks. The transform is done on the time-averaged surface density, over a single binary orbit at 15 000 PAB, 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
thumbnail Fig. 6

Instantaneous surface density profiles for Model H (Kepler-34) at three times near quasi-steady-state separated by 640 PAB.

Open with DEXTER
In the text
thumbnail Fig. 7

Comparison of Kepler-16 disk eccentricity evolution and time averaged surface density for 0.5 ΣMMSN density disks with (I, J, K and L) and without self-gravity (A, B, E and F). Self-gravitating runs are shown with dashed lines.

Open with DEXTER
In the text
thumbnail Fig. 8

Evolution of the mean disk eccentricity for standard density self-gravitating Kepler-16 run J and increased density runs OT (dashed lines). The corresponding simulation with self-gravity disabled is shown by solid lines.

Open with DEXTER
In the text
thumbnail Fig. 9

Fourier analysis of the surface density of Kepler-16 simulations testing the relevance of self-gravity 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 self-gravity is enabled are shown by dotted lines. Fourier analysis is done on the time-averaged surface density over the period of one binary orbit from 15 000 PAB.

Open with DEXTER
In the text
thumbnail Fig. 10

Comparison of surface density and eccentricity evolution in Kepler-16 runs between density profiles of form Σ(r)Σ0rαΣ with αΣ = 0.5 (blue dashed) and 1.5 (black solid).

Open with DEXTER
In the text
thumbnail Fig. 11

Comparison of open (dashed) and closed (solid) boundary conditions in Kepler-16 (black) and Kepler-34 (blue) via surface density and eccentricity profiles.

Open with DEXTER
In the text
thumbnail Fig. 12

Comparison of surface density and eccentricity profiles in Kepler-16 runs between self-gravitating (dashed) and non self-gravitating (solid) for 0.5 ΣMMSN (black) and 12.5 ΣMMSN (red). The 12.5 ΣMMSN surface density profile is normalised to the initial half minimum-mass solar nebula surface density.

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.