Modeling circumbinary planets: The case of Kepler38^{⋆}
^{1} Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
email: wilhelm.kley@unituebingen.de
^{2} Institute for Astronomy and NASA Astrobiology Institute, University of HawaiiManoa, Honolulu, HI 96825, USA
email: nader@ifa.hawaii.edu
Received: 11 December 2013
Accepted: 28 January 2014
Context. Recently, a number of planets orbiting binary stars have been discovered by the Kepler space telescope. In a few systems the planets reside close to the dynamical stability limit. Owing to the difficulty of forming planets in such close orbits, it is believed that they have formed farther out in the disk and migrated to their present locations.
Aims. Our goal is to construct more realistic models of planet migration in circumbinary disks and to determine the final position of these planets more accurately. In our work, we focus on the system Kepler38 where the planet is close to the stability limit.
Methods. The evolution of the circumbinary disk is studied using twodimensional hydrodynamical simulations. We study locally isothermal disks as well as more realistic models that include full viscous heating, radiative cooling from the disk surfaces, and radiative diffusion in the disk midplane. After the disk has been brought into a quasiequilibrium state, a 115 Earthmass planet is embedded and its evolution is followed.
Results. In all cases the planets stop inward migration near the inner edge of the disk. In isothermal disks with a typical disk scale height of H/r = 0.05, the final outcome agrees very well with the observed location of planet Kepler38b. For the radiative models, the disk thickness and location of the inner edge is determined by the mass in the system. For surface densities on the order of 3000 g/cm^{2} at 1 AU, the inner gap lies close to the binary and planets stop in the region between the 5:1 and 4:1 meanmotion resonances with the binary. A model with a disk with approximately a quarter of the mass yields a final position very close to the observed one.
Conclusions. For planets migrating in circumbinary disks, the final position is dictated by the structure of the disk. Knowing the observed orbits of circumbinary planets, radiative disk simulations with embedded planets can provide important information on the physical state of the system during the final stages of its evolution.
Key words: planetdisk interactions / binaries: general / hydrodynamics
Movies are available in electronic form at http://www.aanda.org
© ESO, 2014
1. Introduction
An interesting development in exoplanetary science has been the recent discovery of circumbinary planets (CBPs) by the Kepler space telescope. In contrast to systems where a planet revolves around one star of a binary, here the planet orbits the entire binary system. Currently known main sequence binaries with circumbinary planets are: Kepler16 (Doyle et al. 2011), Kepler34 and 35 (Welsh et al. 2012), Kepler38 (Orosz et al. 2012a), Kepler47 (Orosz et al. 2012b), and Kepler64 (Schwamb et al. 2013). In all these systems, the binary is close with an orbital period of 7 to 41 days. The orbital periods of their planets range from 50 to 300 days.
Since the discovery of the first circumbinary planet in the Kepler16 system, it has been noted that in some of these systems, the innermost planet is very close to the boundary of the dynamical stability limit (Dvorak 1986; Holman & Wiegert 1999). The observations indicate that the orbital planes of these binaries are perfectly aligned with those of their planets, and this implies that the planets are formed in flat, circumbinary disks. The proximity of the orbits of these CBPs to the stability limit raises the question as to whether these planets formed in their current orbits, or at greater distances from the binary and then migrated to their present locations.
Although not all aspects of planet formation are fully understood, it is widely accepted that planets form around single stars through a bottomup process where growth is achieved via a sequence of sticking collisions with subsequent gas accretion for the more massive objects (Armitage 2010). In a binary star system, the presence of the binary in the middle of the protoplanetary disk will alter this process and make planet formation more complicated. The perturbation of the binary and its (eccentric) disk will dynamically excite the orbits of planetesimals and hinder their growth into larger objects (Scholl et al. 2007; Meschiari 2012; Marzari et al. 2013). In a circumbinary disk the outer edge of the central cavity, created as a result of the effect of tidal forces from the binary on the disk material, coincides with the boundary of the planetary stability, and so the previous statement implies that the in situ formation of CBPs close to the stability limit may not be possible (Paardekooper et al. 2012). However, because of the small separations of CBPhosting binaries (each of the abovementioned systems fits into the orbit of Mercury), their effects on the formation of planets at great distances will be negligibly small, suggesting that these planets could have formed farther out and migrated to their current orbits (Dunhill & Alexander 2013; Marzari et al. 2013).
Planet migration is a natural consequence of planetdisk interaction (Nelson et al. 2000; Kley & Nelson 2012). To study the migration of planets in circumbinary disks, the structure of the disk has to be analyzed and compared to those around single stars. The most important dynamical difference between the two kinds of disks is the existence of a central cavity in the circumbinary disks. As shown by Artymowicz & Lubow (1994), the radial extent of this cavity is a function of the binary semimajor axis, eccentricity, and massratio as well as the disk viscosity. As indicated by these authors, for typical values of the disk viscosity, and depending on the binary eccentricity, the location of the inner edge of the disk will be about 2 to 3 times the separation of the binary.
The migration of planets in circumbinary disks was first studied by Nelson (2003) for massive, Jupitertype planets. He found that for low values of the binary eccentricity, the migrating planet will be captured in a 4:1 meanmotion resonance (MMR) with the binary whereas in more eccentric systems (e_{bin} ≳ 0.2), the planet would be captured in a stable orbit farther out (see also Pierens & Nelson 2008b). Subsequent studies by Pierens & Nelson (2007) extended these analyses to planets with masses as low as 20 M_{Earth} showing that these planets usually stop near the edge of the cavity and they suggested that planets in circumbinary disks should be predominantly found in that area. As the inner edge of the disk roughly coincides with the boundary of planetary stability, the predictions of these authors turned out to agree very closely with the orbital architecture of several Kepler CBPs. More general cases in which accretion and multiple planets were also considered were later studied by the same authors (Pierens & Nelson 2008a,b).
In a recent paper Pierens & Nelson (2013) revisited planet migration in circumbinary disks and developed models to explain the orbits of the planets around binaries Kepler16, Kepler34, and Kepler35. Similar to the majority of the simulations of this kind, these authors considered a locally isothermal approximation where the disk thermodynamics is modeled by prescribing a given temperature profile. They also considered a closed boundary condition at the edge of the cavity assuming zero mass accretion onto the central binary. Although the work by Pierens & Nelson (2013) represents significant results, their applicability may be limited because of the simplifying assumptions of an isothermal disk with a closed boundary conditions for the central cavity. Numerical simulations by Artymowicz & Lubow (1996) and Günther & Kley (2002) have shown that despite the appearance of the abovementioned cavity in the center of a circumbinary disk, material can still flow inside and onto the central binary. In other words, a more realistic boundary condition would allow for the inflow of the material through the inner edge of the disk into the disk cavity.
Here, we present an improved and extended disk model that allows for much more realistic simulations. We consider the system of Kepler38 as this system represents one of the binaries in which the circumbinary planet is in close proximity to the stability limit. We have developed an improved approach which includes a detailed balance of viscous heating and radiative cooling from the surface of the disk (D’Angelo et al. 2003), as well as additional radiative diffusion in the plane of the disk (Kley & Crida 2008). For a planet embedded in disks, this improved thermodynamics can have dramatic effects on the planet orbital dynamics such that for lowmass planets, it can even reverse the direction of migration (Kley & Crida 2008; Kley et al. 2009). In addition, unlike Pierens & Nelson (2013), we allow free inflow of material from the disk into the central cavity and construct models with net mass inflow through the disk.
This paper is organized as follows. In Sect. 2, we present the physical setup of our disk models. In Sect. 3, we study locally isothermal disks and in Sect. 4, we present the results of our full radiative models and compare them to the results from standard isothermal cases. Section 5 concludes this study by summarizing the results and discussing their implications.
2. The hydrodynamic model
As mentioned above, we consider Kepler38 to be our test binary star and carry out simulations for that system. To model the evolution of a disk around a binary star, we consider a circumbinary disk around Kepler38 and assume that the disk is vertically thin and the system is coplanar. The assumption of coplanarity is well justified as the mutual inclination of Kepler38 and its planet is smaller than 0.2 degrees. We then perform twodimensional (2D) hydrodynamical simulations in the plane of the binary.
To simulate the evolution of the disk, the viscous hydrodynamic equations are solved in a polar coordinate system (r,φ) with its origin at the center of mass of the binary. In our standard model, which we use as a basis for comparing the results of our subsequent simulations, we consider the radial extent of the disk (r) to be from 0.25 AU to 2.0 AU and φ to vary in an entire annulus of [0,360°]. This domain is covered by an equidistant grid of 256 × 512 gridcells. For testing purposes, we have also used a grid twice as fine. In all models we evolve the vertically integrated equations for the surface density Σ and the velocity components (v_{r},v_{φ}).
When simulating locally isothermal disks, we do not evolve the energy equation, and instead use an isothermal equation of state for the pressure. When considering radiative models, a vertically averaged energy equation is used, which evolves the temperature of midplane (Müller & Kley 2012). Radiative effects are included in two ways. First, a cooling term is considered to account for the radiative loss from the disk surface (D’Angelo et al. 2003). Second, we include diffusive radiative transport in the midplane of the disk using fluxlimited diffusion (Kley & Crida 2008). The use of a fluxlimiter is required because near its inner edge, the disk is very optically thin, which would lead to an unphysical high energy flux in the diffusive part of radiative transport equations. For more details on the implementation of this diffusive transport, see Müller & Kley (2013). In our radiative simulations, we consider the full, vertically averaged dissipation (Müller & Kley 2012). However, stellar irradiation is not taken into account.
To calculate the necessary height of the disk H at a position r, we first note that in a circumstellar disk (around a single star at the position r_{∗}), the vertical height of the disk is given by (1)In this equation, M_{∗} is the mass of the central star, c_{s} is the speed of sound in the disk at the location r, and G is the gravitational constant. In a binary star system, H(r) has to be calculated by taking the contributions of both stars into account (Günther & Kley 2002). That is, (2)where M_{i} = M_{1} and M_{2} are the masses of the stars of the binary. Equation (2) indicates that in a circumbinary disk, the total height H is always smaller than the individual heights H_{i}.
The equation of state of the gas in the disk is given by the ideal gas law using a mean molecular weight of μ = 2.35 (in atomic mass units) and an adiabatic exponent of γ = 1.4. For the shear viscosity we use the αparametrization with α = 0.01 for our standard model, and we set the bulk viscosity to zero. For the Rosseland opacity, we use analytic formula as provided by Lin & Papaloizou (1985) and the fluxlimiter as in Kley (1989).
To calculate the gravitational potential of the binaryplanet system at a position r in the disk, for all three bodies we use a smoothed potential function of the form (3)where M_{k} denotes the masses of the objects (stars and planet), r − r_{k} is the vector from a point in the disk to each object, and H is calculated using Eq. (2). The quantity ϵ in Eq. (3) is a smoothing parameter that accounts for the effect of the finite thickness of the disk. In our simulations, we assume ϵ = 0.6 or 0.7 (Müller et al. 2012). We note that for the planet, the smoothing length cannot be smaller than ϵ r_{H}, where r_{H} = a_{p}(m_{p}/3M_{bin})^{1/3} is the radius of the Hill sphere around the planet, a_{p} is the semimajor axis of the planet’s orbit, and M_{bin} = M_{1} + M_{2}.
The hydrodynamical equations are integrated using a mixed, explicit/implicit, scheme. The standard hydrodynamical terms are integrated explicitly using a standard Courant condition with Courant number f_{C} = 0.66. The viscous and the diffusive radiative parts are integrated implicitly to avoid instabilities. For the timeintegration of the hydrodynamical equations, we utilize the RH2Dcode with the FARGO (Masset 2000) upgrade.
The motions of the stellar binary and the planet are integrated using a fourth order RungeKutta integrator. All objects, including the two stars, feel the gravity of the disk as well as their mutual gravitation. The selfgravity of the disk is not considered. The indirect terms are taken into account through a shift of the positions of the objects such that the binary’s barycenter always coincides with the origin of the coordinate system.
To calculate the force from the disk acting on the planet, we exclude parts of the Hill sphere of the planet using a tapering function (4)In this equation d is the distance from the grid cell to the planet and p is a dimensionless parameter that is set to 0.6 for the isothermal models and 0.8 for the radiative runs. We calculate the orbital parameters of the planet using Jacobian coordinates, assuming the planet orbits a star with mass M_{bin} = M_{1} + M_{2}, at the binary barycenter.
2.1. Initial setup and boundary conditions
In all models, the disk initial surface density is chosen to have a Σ(r) = Σ_{0} r^{− 1/2} profile where r is the distance from the center of mass of the binary and Σ_{0} = 3000g/cm^{2} is the surface density at 1 AU. This radial profile for Σ corresponds to that usually adopted for the minimum solar mass nebula. As the total stellar mass in the Kepler38 system is approximately 1.2 solar masses, one can assume that, to a first order of approximation, the mass of the circumbinary disk would lie in the same range as that of the disk around the protosun. Assuming that the observed circumbinary planets reached their current positions in the late phase of their evolution, this disk surface density may be on the high side. However, together with the relatively high value of the viscosity parameter (α = 0.01), it allows for a sufficiently fast evolution of the system in our standard model such that it can be numerically simulated. We will vary these parameters for the radiative models below.
We choose the initial temperature of the disk to vary with r as T(r) ∝ r^{1} such that, assuming a central star of M_{bin}, the vertical thickness of the disk will always maintain the condition H/r = 0.05. The initial angular velocity of the disk at a distance r is chosen to be equal to the Keplerian velocity at that distance, and the radial velocity is set to zero.
Binary parameter and observed planetary parameter of the Kepler38 system used in the simulations.
The parameters of the central binary have been adopted from Orosz et al. (2012a) and can be found in Table 1, together with the planetary data. Owing to the interaction of the binary with the disk (and planet), these parameters will slowly vary during a simulation. For the models with an embedded planet, we consider the planet to have a mass of m_{p} = 3.63 × 10^{4} M_{1} which is equivalent to m_{p} = 0.34 M_{Jup} or about 115 Earth masses (see Table 1). At the beginning of each simulation, we start the binary at its periastron.
The boundary conditions of the simulations are constructed such that at the outer boundary of the disk, r_{max} = 2.0 AU, the surface density remains constant. This is achieved by using a damping boundary condition where the density is relaxed toward its initial value, Σ(r_{max}) = 2121g/cm^{2}, and the radial velocity is damped toward zero. To do this, we use the procedure specified by de ValBorro et al. (2006). The angular velocity at the outer boundary is also kept at the initial Keplerian value, and for the temperature we use a reflecting condition so that there will be no artificial radiative flux through the outer boundary. These conditions at the outer boundary lead to a disk with zero eccentricity at r_{max}. Hence, r_{max} has to be large enough such that the inner regions are not influenced.
At the inner boundary (r_{min}), we consider a boundary condition such that the inflow of disk material onto the binary is allowed. This means, for the radial boundary grid cells at r_{min}, we choose a zerogradient mass outflow condition, where the material can freely leave the grid and flow onto the binary. No mass inflow into a grid is allowed at r_{min}. The zerogradient condition is also applied to the angular velocity of the material since because of the effect of the binary, no welldefined Keplerian velocity can be found that could be used otherwise. This zerogradient condition for the angular velocity implies a physically more realistic zerotorque boundary.
With these boundary conditions, the disk can reach a quasistationary state in which there will be a constant mass flow through the disk.
3. The locally isothermal case
In this section we do not evolve the energy equation but instead leave the temperature of the disk at its initial value. This procedure has the advantage of making the simulations much faster as no heating and cooling of the disk have to be considered. This is an approximation that has often been used in planet disk simulations and gives the first order results, but it also has its limitations when planetary migration is concerned (Kley & Crida 2008). We will first show the results of our standard case (model 1), and then compare these results to those recently obtained by Pierens & Nelson (2013) who also used a locally isothermal model. Table 2 gives a very brief overview of the models and can serve as a reference.
Properties for the locally isothermal disk models without a planet.
We begin by discussing models without a planet where we bring the disk into equilibrium and analyze its structure. This equilibration process is required because the disk structure is determined by the action of the binary for which no simple analytic model can be prescribed. Additionally, the evolution of the planet in the disk occurs so slowly that the disk always remains near equilibrium.
3.1. The structure of the circumbinary disk
Fig. 1 Total mass of the disk (without the planet) in our locally isothermal standard model (model 1). After a few thousand years, the disk reaches equilibrium and a nearly constant mass is established. 

Open with DEXTER 
To analyse the disk’s own dynamics (i.e., without a planet) due to the effect of the central binary, we simulated the dynamics of the disk in our standard model with a zeromass planet. Because the disk inner boundary is open, and at the outer boundary its surface density is considered to be constant, the disk settles into a final quasistationary state in which the mass of the disk is constant and there is a constant mass accretion rate onto the central binary. This has been shown in Fig. 1 where the total mass in the computational domain is plotted vs. time for the locally isothermal model with no planet. After about 2000−3000 yr a quasistationary state is reached where the outflow through the inner boundary is exactly balanced by the “inflow” from outside as established by keeping the value of the density constant at r_{max}. During this process, the disk mass decreases by about 20% from 3.8 × 10^{3} M_{⊙} to 3.04 × 10^{3} M_{⊙}. The small remaining oscillations are caused by the eccentric motion of the inner binary which induces variations of the disk structure. In this simulation where we are interested in the structure of the disk without the planet, the disk does not affect the motion of the binary. Nevertheless, we numerically integrate the motion of the binary to examine the accuracy of the numerical method. We find that during the total evolution time of about 4000 years, which corresponds to over 80 000 binary orbits and 4.5 million time steps, the semimajor axis of the binary has only shrunk by 0.014%, a value low enough for our purposes.
At equilibrium, the rate of the mass flow through the disk, Ṁ_{disk}, becomes on average nearly constant. Figure 2 shows this where Ṁ_{disk}, measured at two different distances, is plotted near the end of the simulation. The resulting accretion rate equals approximately Ṁ_{disk} = 8.7 × 10^{7}M_{⊙}/yr. The oscillations are much larger at the inner parts of the disk because of the perturbation from the binary, and they occur at the binary’s orbital period. This value may be compared to the typical theoretical equilibrium disk accretion rate, Ṁ_{th} = 3πΣν. At the outer boundary, where the mass density has been kept constant, we find Ṁ_{th} = 5.4 × 10^{7}M_{⊙}/yr.
Fig. 2 Rate of mass accretion through the disk in model 1 measured at two different distances. The mean value, −8.68 × 10^{7}, is also plotted for r = 0.68 AU. As the central binary accretes material, the mass flow is directed inward and becomes negative. In the main text, however, we report the accretion rate in the form of positive absolute values. 

Open with DEXTER 
Fig. 3 Azimuthally averaged surface density and disk eccentricity for model 1. Displayed are the profiles at different evolutionary times (in years); the red curve denotes the initial setup. 

Open with DEXTER 
During time, the surface density slowly evolves away from its initial profile until it settles in a new equilibrium state. This is shown in the upper panel of Fig. 3 where the azimuthally averaged radial surface density profile is shown at different evolutionary times. In agreement with Fig. 1, the equilibration time takes about 2000 yr. At this state, the surface density profiles will no longer change with time. Because of the tidal action of the binary on the disk, a central gap is formed with a surface density many orders of magnitude lower than inside the disk. The inner edge of the disk, which we can define, very approximately, as that radius at which the surface density is about half the maximum value, lies here at around r ≈ 0.45 AU. This is slightly larger than 3a_{bin} and in a good agreement with the findings of Artymowicz & Lubow (1994). We note that for the binary system considered here, the stability limit for planetary orbits is about 0.4 AU (Dvorak 1986; Holman & Wiegert 1999). Outside of the gap beyond r = 0.6 AU, the surface density profile becomes relatively flat. One can notice that the massflow rates shown in Fig. 2 correspond to a region deep inside the gap (r = 0.31 ≈ 2a_{bin}), and immediately outside of it (r = 0.68).
As was shown by Pierens & Nelson (2013), circumbinary disks can attain significant eccentricities. We calculate the eccentricity of the disk by treating each grid cell as an individual particle with a mass and velocity equal to the mass and velocity of the cell (Kley & Dirksen 2006). To calculate a radial dependence for the disk eccentricity, e_{disk}(r), we average over the angular direction. In our simulations, the disk eccentricity remains low, as shown in the lower panel of Fig. 3. In the regions outside of the central gap (beyond r = 0.6), the disk eccentricity is always below e_{disk} < 0.07. At radial distances r > 1.0, this eccentricity becomes lower than about 0.01. These values of the disk eccentricity are in contrast to those reported by Pierens & Nelson (2013), who found, on average, higher values. We attribute this difference to the assumption of an open inner boundary, and we will analyse this in more detail in the next section.
Fig. 4 Twodimensional density structure of an isothermal disk (model 1) around the central binary star. The graph shows a local view around the binary. The computational grid extends from r = 0.25 AU to r = 2.0 AU. The white inner region lies inside the computational domain and is not covered by the grid. The positions of the primary and secondary are indicated by the gray dots. The Roche lobe of the secondary is also shown. (Online movie) 

Open with DEXTER 
In Fig. 4, we show the 2D surface density distribution for our isothermal disk models. We note that the figure shows only the inner part of the computational domain around the central binary. The Roche lobe of the secondary star is shown as well. As shown here, an eccentric central binary strongly perturbs the disk and produces time varying patterns. As indicated by the flat surface density profile outside of the central gap, the disk shows less structure and is relatively homogeneous.
Fig. 5 Azimuthally averaged surface density and disk eccentricity profiles for two models with different inner boundary conditions. The red curve corresponds to the standard model with an open inner boundary and the blue curve represents a system with a reflecting inner boundary (model 2). In both models the disk is at equilibrium. 

Open with DEXTER 
3.2. A closed inner boundary
In this section, we study the effect of different boundary conditions at the inner edge of the computational domain on the final results. Inspired by the work of Pierens & Nelson (2013), who focused on models with closed inner boundaries, we carried out simulations using their boundary conditions. In Fig. 5, we show the surface density and disk eccentricity for two simulations; our standard model with an open inner boundary, and disk model 2 (see Table 1) with reflecting conditions at r_{min} as used by Pierens & Nelson (2013). In the latter model material is not allowed to leave the disk through the inner boundary, and so it accumulates near the outer edge of the gap, as seen by the spike in the density at r = 0.8 AU in the top panel of Fig. 5. At the same time, the eccentricity of the disk in this model becomes substantially higher than that in the model with the open boundary condition. Our results of the simulations with the closed inner boundary are very similar to those reported by Pierens & Nelson (2013), although these authors did not model the system of Kepler38, but some related ones.
Fig. 6 Azimuthally averaged surface density and eccentricity profiles for three disk models with different locations of the inner boundary, r_{rmin}. The red curve shows the standard model (model 1) with r_{min} = 0.25 AU, the blue curve is for a disk with r_{min} = 0.22 AU, and the green curve represents the disk model with r_{min} = 0.20 AU (model 3). In all three models the disk is at equilibrium. 

Open with DEXTER 
3.3. Varying the location of the inner radius
In this section, we examine the effect of the location of the inner boundary on the disk structure (model 3). Pierens & Nelson (2013) found that, at least in the case of a closed inner boundary, the amplitude of the disk eccentricity depends on the location of r_{min}. In particular, if r_{min} is so large that some MMRs between the binary and the disk fall outside the computational domain, the eccentricity of the disk will remain very low. To examine whether the disk will behave similarly for an open inner boundary, we carried out two additional simulations with slightly smaller inner radii, r_{min} = 0.22 AU and 0.20 AU. Given the binary’s semimajor, a_{bin} = 0.1469 AU, these radii correspond to 1.5 a_{bin} and 1.36 a_{bin}, respectively.
As shown in Fig. 6 (top panel), for low values of r_{min}, the density of the disk increases outside of the gap, and at the same time, the flow of the mass onto the binary is slightly reduced (to about 8 × 10^{7}M_{⊙}/yr). The profile of the edge of the gap, where the slope of the density is positive, is hardly affected by the choice of r_{min}. However, the gap is slightly wider. The eccentricity in the main part of the disk (outside r = 0.7) does not vary with the location of r_{min} (lower panel of Fig. 6). At r = 1 the disk eccentricity is about e_{disk} = 0.1 and beyond r = 1.0 it is always below 0.01. Within the gap region where the density is low, e_{disk} rises to a maximum of 0.4 at the inner boundary for the models with the smaller r_{min}. The disks with r_{min} = 0.2 AU and r_{min} = 0.22 AU show very similar behaviour. For the sake of comparison, we also ran simulations with higher spatial resolution. The results were very similar to those depicted by the two panels of Fig. 6.
Overall, the results of our simulations indicate that in a model with an open inner boundary, the disk eccentricity e_{disk} remains low in a large portion of the disk and independent of the location of the inner boundary. This is unlike the findings of Pierens & Nelson (2013) who showed that in cases with close boundaries, the disk eccentricity will rise to much higher values (Pierens & Nelson 2013) because a reflecting inner boundary creates a cavity that can sustain an eccentric global mode of the disk. Our findings are also in agreement with Kley et al. (2008) who found an eccentricity reduction for circumstellar disks in close binary systems when using outflow boundary conditions. In our opinion, open inner boundaries are more realistic because they allow for the accretion of material onto the binary stars. However, we would like to caution that in the case of more eccentric binaries, it has been found that an open inner boundary may cause the central cavity to be very large since disk material cannot reenter the computational domain once it has been lost through r_{min} (see Marzari et al. 2009). This situation occurs primarily in highly eccentric binaries such as Kepler34 which has an eccentricity of 0.52 (A. Pierens, private communication). For our system, Kepler38, we found the opposite behaviour; a reduced disk eccentricity for an open inner boundary. This follows from the fact that an open inner boundary reduces the strength of a global eccentric mode. In this paper we are more interested in the evolution of planets in circumbinary disks, and so we will defer the study of this interesting behaviour of the disk to a future time.
We note that among the three models shown in Fig. 6, the model with the smallest inner radii can have a considerably smaller numerical time step. In addition to the smaller size of the grid cells, this effect is caused by the higher deviation from a uniform Keplerian velocity as induced by the binary star, which makes the FARGO algorithm less effective. Therefore, to save computational time and resources, in the remainder of this study, we will use for the inner radius of the disk, a standard value of r_{min} = 0.25 AU, although we will present some comparison simulations using r_{min} = 0.22 AU as well.
3.4. Planets in isothermal disks
Fig. 7 Evolution of the semimajor axis of the planet in an isothermal disk (model 1). In two simulations, the planet is started in a circular orbit at a distance of a_{0} = 1.75 AU from the center of mass of the central binary. In the simulation shown in red, the planet is included in the disk from the very beginning of the simulation and it evolves as the disk evolves starting from its initial density profile. In the simulation shown in blue, the disk is first relaxed to the equilibrium and then the planet is embedded in it. In the simulation shown in green, the planet is started at a_{0} = 1.0 AU in a relaxed disk. The inset shows the results near the end of the simulation for the first model. The dashed black line in the inset corresponds to the mean value of the planet’s semimajor axis at a_{p} = 0.436 AU. 

Open with DEXTER 
To study the evolution of the system with an embedded planet, we start our simulations with a planet at different initial distances (semimajor axes, a_{0}) from the center of mass (barycenter) of the binary and in a circular orbit. The planet’s evolution is determined by the gravitational action of the disk and the central binary. As mentioned earlier, during the evolution of the planet, its orbital elements are calculated using Jacobian coordinates with respect to the barycenter of the central binary. The planet mass is fixed to 115 Earth masses, and there is no accretion onto the planet. As a reference, we present in Table 1, the orbital parameters of the planet Kepler38, as inferred from the observations.
Figure 7 shows the evolution of the planet through the disk. We present here the results of three simulations that were carried out for different initial setups to examine the dependence of the outcome on the initial state of the disk and the planet’s orbital elements. In the first model (red line) the planet is started directly in the initial disk which has not been brought into an equilibrium (i.e., using the initial disk setup of model 1). In this simulation, the disk and the planet evolve simultaneously. In the other two simulations, the planet is started in a disk that is already in an equilibrium. The blue line corresponds to the planet starting at a_{0} = 1.75 AU and the green line is when the planet starts at a = 1.0 AU. We have shifted the last model in time to ensure an overlap with the first two.
As shown in the Fig. 7, the semimajor axis of the planet continuously shrinks (i.e., the planet migrates toward the central binary) until it reaches the inner cavity created by the binary. The strong drop in the disk surface density within the central cavity and the resulting positive density gradient act as a planet trap (Masset et al. 2006). Here the negative Lindblad torques are balanced by positive corotation torques. As a result, the planet stops its inward migration in this region. As shown in the figure, in all runs, the planet stops at approximately the same distance from the binary. The semimajor axis of the planet in the three runs has an average value of a_{p} = 0.436 AU (black line in the inset in Fig. 7). The results indicate that in all simulations, the system evolves to the same final state. At this state, the planet resides in an orbit with a period of about 100 days which is slightly outside of its observed orbit (see Table 1). This implies that the outcome of the migration process does not depend on the history of the system and is determined solely by the physical parameters of the disk. In an isothermal system and for a given binary, these parameters are determined by the values of H/r and viscosity.
Fig. 8 Evolution of the eccentricity of the planet in an isothermal disk shown in red in Fig. 7 (model 1). The planet was embedded in the disk prior to the start of the disk evolution at a semimajor of a_{0} = 1.75 AU. 

Open with DEXTER 
Figure 8 shows the evolution of the eccentricity of the planet corresponding to the red model in Fig. 7. As shown here, during the planet’s inward migration, its eccentricity increases continuously until t = 4500 when it reaches a constant value close to 0.05 and oscillates around that value for approximately 1000 years. At time (t ~ 5700), the planet is temporarily captured in a meanmotion resonance with the binary and its eccentricity undergoes a sudden jump (Murray & Dermott 1999) where it oscillates between 0.15 and 0.20. We will discuss this resonant capture in more detail in the next section.
Figure 9 shows four 2D snapshots of the evolution of the disk surface density during the inward migration of an embedded planet in one of our models. One can clearly see the spiral disturbances created by the planet during its migration. The last panel, corresponding to t = 7417, shows the final position of the planet in an orbit with a semimajor axis of a_{p} = 0.436 AU (see Fig. 7). In Fig. 10 we plot the azimuthally averaged surface density of the disk corresponding to the same times as in the snapshots in Fig. 9. The colored circles in this figure represent the planet in its actual radial distance from the center of mass of the binary. A planet of this mass does not produce a full gap, but only clears out about a third of the surface density. As the planet moves inward, the surface density of the inner part of the disk is strongly distorted. However, it regains its unperturbed profile after the planet has reached its final position (see Fig. 3).
3.4.1. Capture in resonance
The evolution of the semimajor axis of the planet and in particular, the sudden increase in its eccentricity at time t ≈ 5700 seem to point to a temporary capture in a meanmotion resonance with the central binary. To analyze this, we plotted the variations of the periodratio of the planet and binary, P_{planet}/P_{bin}, during the inward migration of the planet. Figure 11 shows the results. As shown here, as the planet migrates inward this periodratio continuously decreases until it reaches a constant value at t ≈ 5000. The inset in Fig. 11 shows the evolution of the system near the end of the simulation for one thousand years. The black dashed line in this inset denotes the arithmetic averaged value of the period ratio, P_{planet}/P_{bin} = 5.22, suggesting a possible capture into a 5:1 MMR with the binary.
Fig. 9 Twodimensional snapshots of the evolution of the disk surface density with an embedded planet. The snapshots correspond to the times 449, 729, 2807, and 7414 years. The planet, shown as the gray dot with the Roche radius, was embedded in the disk at an initial distance of r = 1.75 AU. The dashed gray line around the central binary shows the approximate boundary of stability for planetary orbits. A movie of the simulation shown here can be found online. 

Open with DEXTER 
Fig. 10 Averaged radial surface density of a disk with an embedded planet. The colors correspond to the times of the snapshots shown in Fig. 9. The big circle in each graph represents the radial position of the planet at the time corresponding to the graph. For illustrative purposes, we have moved the circles close to their corresponding curves. 

Open with DEXTER 
To determine whether the planet is truly captured in a 5:1 MMR, its corresponding resonant angles have to be analyzed. The capture in a resonance between two bodies occurs when at least one of the resonant angles of the system librates around a certain value and does not cover the full range of 0 to 2π, see Murray & Dermott (1999) for details. For a general p:q commensurability with p > q and in a planar system, the resonant angles are defined as (5)where λ_{b},λ_{p},ϖ_{b}, and ϖ_{p} denote the mean longitudes and longitude of periapse for the inner binary (b) and outer planet (p). The integer k in this equation satisfies the condition q ≤ k ≤ p and has p − q + 1 possible values. Of the p − q + 1 resonant angles, at most two are linearly independent. That means, in an actual resonant configuration, at least one of these angles will librate (Nelson & Papaloizou 2002). As our numerical simulations indicate the possibility of a capture in the 5:1 MMR, we studied the timevariations of all resonant angles of our system. We note that in this case p = 5, q = 1, and k runs from 1 to 5. Figure 12 shows the evolution of the following two resonant angles for 180 years from t = 5620 to t = 5800, (6)and (7)As shown here, before t = 5670, the evolution of these angles is irregular confirming that the planet was not in resonance. After this time, Φ_{1} shows a regular circulating behaviour, whereas the angle Φ_{2} begins to librate around zero with an amplitude of less than 40 degrees. This librating evolution of Φ_{2} indicates that the planet is indeed captured in a 5:1 MMR with the binary star. Our analysis of all other resonant angles of the system, Φ_{3}, Φ_{4}, and Φ_{5}, showed that these all have circulating behaviours as well, either prograde (as Φ_{1}) or retrograde. For all cases shown in Fig. 7, the planet remains captured in the 5:1 resonance for the duration of the simulations. We also studied the behaviour of the longitude of periastron of the binary star, ϖ_{bin}, during the slow inward migration of the planet. Our analysis indicated that prior to the capture of the planet in a 5:1 MMR, ϖ_{bin} was precessing in a prograde sense. After capture into the 5:1 resonance, this precession became retrograde.
Fig. 11 Evolution of the period ratio of the planet and binary for the model where the planet is embedded in the initial disk at a_{0} = 1.75 AU. The inset shows the period ratio for the final stage of its evolution. The dashed line in the inset corresponds to the mean value of P_{planet}/P_{bin} = 5.22. 

Open with DEXTER 
Fig. 12 Resonant angles Φ_{1} and Φ_{2} (see Eqs. (6) and (7)) during the capture phase in a 5:1 meanmotion resonance. 

Open with DEXTER 
Fig. 13 Evolution of the semimajor axis of a planet embedded in an isothermal disk for different values of the disk mass. The graph in red corresponds to the simulation of the standard disk model with the planet initially starting at a_{0} = 1.0 AU. The graph in blue shows a similar simulation but with half the disk mass. The graph in green corresponds to a simulation with a quarter of the initial disk mass that is continued from the model in blue at t = 4300. 

Open with DEXTER 
3.4.2. Variation of disk mass
Up to this point, the presented results and analyses were based on simulations in which, as shown in Fig. 1, the disk had a given mass. To examine whether the results will depend on the mass of the disk, we carried out simulations where the disk surface density was reduced by a factor of two and four. The results indicated that disks with reduced surface density show similar density evolution to those shown by the top panel in Fig. 3, but only scaled down by an appropriate factor. This identical density evolution can be attributed to the fact that in an isothermal disk with an αtype viscosity and no selfgravity, the surface density of the disk is scaled out of the normalized equations of state making the disk evolution independent of the actual density.
In a disk with a lower mass and surface density, the inward migration of the planet will be slower. Figure 13 shows the evolution of a planet in an isothermal disk with different masses. The red line in this figure represents the inward migration of a planet started at a_{0} = 1.0 AU in our standard disk model. The blue line corresponds to the migration of the same planet in a disk with half the mass of our standard model. As shown here and as expected, the inward migration is about two times slower than in the standard case. We started our third simulation (green line) with a quarter of the disk mass and to save computational time, continued it from the second mode (blue line) at time 4300. The planet migration rate is again reduced by a factor of two. As pointed out by Pierens & Nelson (2013) when analyzing the system of Kepler16, in a system with a slower rate of migration, a capture in higher order MMRs is more expected. These authors found that Kepler16b could have undergone a temporary capture in a 6:1 MMR with the central binary. In our simulations, the planet is always captured in a 5:1 resonance and in the same configuration, where only Φ_{2} librates. The eccentricity of the disk also reaches the same value in all three cases.
Fig. 14 Evolution of the semimajor axis of a planet embedded in an isothermal disk for different locations of the disk’s inner boundary. The graph in red corresponds to the standard model with an open inner boundary at 0.25 AU. The graph in blue represents a model with a closed inner boundary at 0.25 AU, and the graph in green is for a disk with an open inner boundary at 0.22 AU. 

Open with DEXTER 
3.4.3. Planets in disks with different inner boundary conditions at r_{min}
To test the sensitivity of our simulations to the choice of the disk boundary conditions at r_{min}, we ran three different models. In two models we considered r_{min} = 0.25 AU and once we assumed the inner boundary to be closed to mass inflow, and once we considered it to be open (see model 2 in Sect. 3.2). We also ran a third model for which we moved r_{min} to a different inner radius at 0.22 AU (see model 3 in Sect. 3.3). Figure 14 shows the evolution of the planet’s semimajor axis for these three cases. The planet was started at a_{0} = 1.0 AU in an already relaxed disk in all three cases. See Figs. 5 and 6 for the density and eccentricity distribution of models 2 and 3.
Because of the higher density in the simulation with a closed inner boundary (shown in blue), the planet migrates inward at a much faster speed. Similar to the previous simulations, in this case the planet is captured in a 5:1 MMR with the binary. At this state, the planet’s resonance angle Φ_{2} begins to librate around zero and its other resonance angles start circulating. In contrast to the previous simulations, the fast inward migration of the planet causes its orbital eccentricity to be slightly higher (e_{p} = 0.25) than in the case when the disk had an open inner boundary. This higher value of eccentricity leads to orbital instability as a result of which the planet leaves the resonance at time t ≈ 2300 and is scattered out into an orbit with a much larger semimajor axis and very high eccentricity.
In the simulations where the open inner boundary has been shifted to r_{min} = 0.22 AU, the planet migrates inward at a reduced speed and reaches a_{p} = 0.47 AU at the end of the simulation. This is in an excellent agreement with the observed orbit of Kepler38b. The eccentricity of the planet at this state oscillates between 0.01 and 0.06 with an average around 0.03. These results confirm the expectation that the final position of the planet is determined by the location of the outer edge of the central cavity.
In summary, in the model in which the outer edge of the central cavity is at a smaller distance and has the unrealistic boundary condition of being closed to mass inflow, the planet crosses the 5:1 MMR and eventually becomes unstable. When the outer edge of the central cavity is at a larger distance and has the more realistic boundary condition of being open to mass inflow, the planet stops its inward migration outside of the 5:1 MMR and maintains its orbital stability for a long time. We will return to this point when we discuss radiative models in the next section.
Fig. 15 Evolution of the semimajor axis of the binary in a system with an isothermal disk for different values of the disk mass. The graph uses the same models and colors as in Fig. 13 for the evolution of the planets. 

Open with DEXTER 
3.4.4. The orbital elements of the binary
As mentioned earlier in the description of our models and their setups, the motions of all objects (including the stars of the binary) are affected by the gravity of the disk. As a result, the orbit of the central binary will change in time as it interacts with the disk and the planet. In Fig. 15, we plot the evolution of the semimajor axis of the central binary in three systems with disks with masses as in Fig. 13. One can see from this figure that during the evolution of the system, the semimajor axis of the binary reduces with time. This can be explained by noting that while the disk interacts with the binary, angular momentum is transferred from the binary to the disk material causing the orbit of the binary to shrink. As expected, disks with lower masses lead to a slower decline in the orbit of the binary. At the point where the planet is captured into resonance, the binary loses angular momentum more rapidly.
4. The radiative case
As mentioned earlier, although the study of a locally isothermal disk can portray a general picture of the evolution of the system and its planet, it does not present a realistic model. For instance, in an isothermal disk, despite the strong heating action of the spiral waves that are induced by the binary star, the temperature does not vary and is held constant. However, because the migration of the planet and the evolution of its eccentricity depend on the disk temperature (Ward 1988), a more realistic model should allow the disk temperature to vary as well. In this section, we consider the second case. Because the evolved quantity in the simulations is the midplane temperature, the inclusion of an energy equation with viscous dissipation, radiative cooling from the disk surfaces, and diffusion in the disk’s midplane will result in a selfconsistent temperature determination in the disk and make the models much more realistic.
Similar to the case of isothermal disks, we first study the evolution of the disk without a planet, and then evolve the planet in the disk. The models are generally started from the same initial conditions as in the isothermal disks. They are then evolved to an equilibrium state that now takes about 1500 years to reach. We have compared the results to models which were started from the equilibrated isothermal states and found that they both reach identical final configurations. There was no gain in computational time when started from an equilibrium isothermal disk because the timescale for the disk adjustment was identical as it is determined by the radiative and viscous diffusion timescale.
4.1. Disk structure
Figure 16 shows a comparison between the density and temperature of a radiative disk and that of an isothermal model after they have reached their corresponding equilibrium states. The boundary conditions for the disk density are the same in both models. That means the density near the outer boundary is relaxed to the initial value which is held constant. This implies that the density for both models is identical near the outer boundary at r = 2.0 AU. In the inner part of the disk, however, the radiative model has a higher density throughout (top panel in Fig. 16). This high value of density may be due to the higher eccentricity of the disk (bottom panel in Fig. 16), similar to the evolution of an isothermal disk with the closed inner boundary. The total mass of the radiative model in the computational domain is approximately 18% higher than in isothermal cases. The oscillations in the disk mass in the quasiequilibrium state reach ~1.5% which is about ten times greater than those in the isothermal models (see Fig. 1). These wider variations are caused by the greater disk eccentricity in the radiative models, which results in large variations in the mass flow rate through the disk, with an equilibrium value of Ṁ_{disk} ≈ 1.2 × 10^{6} M_{⊙}/yr. Because of higher temperatures in the disk, in these models the inner edge of the disk is slightly closer to the binary than in the isothermal cases.
Fig. 16 Surface density (in g/cm^{2}, top panel), temperature (in K, middle panel), and eccentricity (bottom panel) of a radiative disk in the initial model (red), an isothermal case with H/r = 0.05 (model 1, blue), and a quasiequilibrium state (green). 

Open with DEXTER 
As shown in the middle panel of Fig. 16, the viscously heated disk is much hotter in the regions outside the inner cavity than an isothermal disk with the same initial profile. The increased temperature of the disk results in a greater vertical height of H/r ≈ 0.06 to 0.075 within the main part of the disk, compared to H/r = 0.05 for the locally isothermal disk. Inside the cavity in the radiative model, the temperature becomes very low because the reduced density leads to optically thin disks that can cool more efficiently. Furthermore, the diffusion of radiation in the midplane of the disk leads to additional loss of energy through the open inner boundary (stellar irradiation was not included in these models).
While the disk eccentricity in isothermal models was low, we find that in radiative models, the disk eccentricity becomes very high (see bottom panel in Fig. 16). This high eccentricity is caused by the higher disk temperature in the radiative case. For disks with less mass, the radiative models have lower temperatures and lower eccentricities, (see Sect. 4.3 below). The high eccentricity of the disk in radiative models is in agreement with the results of Pierens & Nelson (2013), who also consistently find higher disk eccentricities for higher H/r. Surprisingly, this behaviour is in contrast to the case where the disk surrounds one star of a binary. In those cases, disks with smaller H/r tend to have higher eccentricities (Müller & Kley 2012). We do not have a full explanation for this behaviour, but we suspect that the tidal truncation of the disk in this case, along with the eccentricity of the binary may play a role. In circumbinary disks, there is no outer radius that may handicap the growth of highly eccentric modes.
Fig. 17 Evolution of the semimajor axis and eccentricity of an embedded planet in a radiative disk. The planet was started at a_{0} = 1.0 AU (red curve). In the two other models shown in green and blue, simulations were continued at specified times with slightly different numerical parameters (see text for details). (Online movie) 

Open with DEXTER 
4.2. Planets in radiative disks
Similar to the isothermal case, after the radiative disk has reached thermal and dynamical quasistatic equilibrium, we place a m_{p} = 0.34 M_{Jup} planet in the disk and follow its evolution. In previous models the final outcome did not depend on the planet’s initial condition, and the evolution of the planet is driven by dissipative processes, and so we expect that in the radiative case, the final outcome will also be independent of the initial condition, and we place the planets in a circular orbit with a semimajor axis of a_{0} = 1.0 AU. This choice of semimajor axis will allows us to save considerable computer time.
Figure 17 shows the evolution of the semimajor axis and eccentricity of the planet. As shown here, the higher disk density in the radiative case causes the planet to migrate inward at a faster pace than in the corresponding isothermal model (see Fig. 7). Similar to the isothermal case, the inward migration slows down when the planet reaches the central cavity with the positive density slope. In the radiative case, this occurs at t ≈ 500 years. After this time, the migration continues and the planet reaches the location of 5:1 MMR with the binary at t ≈ 1400 years where it is captured in that orbit. As expected, upon reaching the resonance, the planet eccentricity increases sharply (bottom panel). Shortly after capture in resonance at about t = 1500 year, the planet’s orbit becomes unstable and the planet is scattered into a highly eccentric orbit with a much larger semimajor axis.
To examine the effects of the numerics on the results, we changed two numerical parameters that can influence the calculations of the gravitational force between planet and the disk, namely the gravitational smoothing parameter ϵ (from 0.6 as in the isothermal case to 0.7) and the force truncation parameter p (from 0.6 to 0.8). We then carried out simulations for the evolution of the disk and planet. The two runs with these new parameters were continued from the first model at t = 1150 and t = 1270 years (green and blue lines in Fig. 17), respectively, shortly before instability set in. In both models, the planet is temporarily captured in a 5:1 resonance with the binary. As shown by the blue and green curves in Fig. 17, at this state (the evolutionary times of 1700 and 1900 years) the semimajor axis of the planet remains constant and its eccentricity increases. Subsequently, the planet leaves the resonance migrating farther inward until its reaches an orbit with a semimajor axis of a ≈ 0.415 AU between the 5:1 and 4:1 resonances where it then becomes stable. At this state, the planet’s eccentricity reduces to an average value of e ≈ 0.03. Compared to the semimajor axis and eccentricity of Kepler38b (see Table 1), the planet in the radiative disk model is closer to the binary; however, its eccentricity is in a similar range. This mismatch between the observed value of the planet’s semimajor axis and that obtained from our simulations may be due to a disk mass that is too high.
Fig. 18 Surface density (in g/cm^{2}, top panel), temperature (in K, middle panel), and eccentricity (bottom panel) of two radiative disks with different masses. The green curve corresponds to the same disk model as the green curve in Fig. 16 with the standard disk mass and viscosity (α = 0.01). The orange curve shows a disk model with a mass equal to a quarter of the mass of our standard disk model and a lower disk viscosity of α = 0.004. 

Open with DEXTER 
Fig. 19 Evolution of the semimajor axis of an embedded planet in a radiative disk with different masses. The red and blue curves correspond to the standard disk mass model (as shown in Fig. 17) and the orange curve corresponds to a model with a quarter disk mass. 

Open with DEXTER 
4.3. Planets in radiative disks with reduced mass and lower viscosity
As indicated by the results of our simulations, the final location of a planet in a circumbinary disk depends on the location of the disk’s inner edge. In the radiative model with the standard disk mass, the inner edge of the disk is quite close to the binary star. As a result, the planet is able to cross the 5:1 resonance, or it becomes unstable after being temporarily captured in that resonance. A survey of the currently known circumbinary planets indicates that in contrast to these results, all the observed planets are in orbits outside the 5:1 resonance. To ensure that the final location of the planet in our simulations will be farther away, we need to develop a disk model with a large inner cavity such that the inner edge of the disk will be at a large distance from the binary star. Two effects will be able to accommodate such a disk model; a reduced temperature, and a lower viscosity.
In a radiative model, the midplane temperature of the disk is determined by the disk opacity, which can be easily adjusted by the masscontent of the disk. To simulate the evolution of the disk and its embedded planet for lower values of the disk temperature, we consider a model with a lower disk mass. Additionally, we reduce the value of the disk viscosity. Specifically, we consider a disk with a quarter of the mass (i.e. a model in which the surface density at the outer boundary was set to a quarter of the original value), and a lower viscosity parameter, α = 0.004. This value of α may be more realistic for protoplanetary disks in their final phases. In Fig. 18, we compare the structure of this disk with that of our standard radiative disk model. From top to bottom, the panels show the radial dependence of the disk surface density, its temperature, and the disk eccentricity at the equilibrium state. The green curves correspond to the quantities in Fig. 16. Because of the disk’s lower mass, the surface density of the disk is naturally lower. This reduces the disk’s optical depth and leads to a lower temperature, which subsequently results in a lower value for the disk thickness H. In agreement with previous results, a lower H results in a lower eccentricity for the disk.
Fig. 20 Evolution of the semimajor axis and eccentricity of a planet embedded in a radiative disk with a mass lower than that of our standard disk model. The orange curve corresponds to the same model as the orange curve in Fig. 19 with the planet started at a_{0} = 1.0 AU. The purple curve corresponds to a system in which the planet was embedded in the disk at a_{0} = 0.5 AU and moved only under the influence of the gravity of the binary. The light blue curve shows the same planet as in the purple curve, where at time t = 1700 the effect of the disk was also taken into account. The dashed line in the bottom panel shows the evolution of the planet’s mean eccentricity. The first equilibrium is reached at e ≈ 0.036, and under the action of the disk, this value of the eccentricity exponentially decreases to a new mean value of e ≈ 0.026 in a timespan of 600 years. 

Open with DEXTER 
We followed our standard procedure, i.e. after allowing the new disk (with lower mass and viscosity) to reach equilibrium, we placed a planet in a circular orbit at a_{0} = 1 AU. Figure 19 shows the evolution of the semimajor axis of the planet for the new quartermass model (orange line) in comparison to the previous simulations. As shown here, in both models, the planet migrates inward and the rate of migration is initially similar in both cases. However, although the reduced density would slow down the migration and the lower temperature would speed it up, in the long run, the planet drifted inward at a much lower rate than in the previous simulation. In this reducedmass disk model with the lower temperature and viscosity, the planet carves a deeper gap where the central density in the vicinity of the planet is only about a quarter to a third of the ambient density compared to the twothirds in the previous model (see Fig. 10). This deeper gap brings the planet closer to the Type II migration regime, and eventually leads to the very slow migration rate observed in Fig. 19. In this simulation, the migration is so slow that we could not follow the model to the final equilibrium state.
To overcome the problem with the slow migration, we ran an additional model where we placed the planet initially in an orbit with a smaller semimajor axis of a_{0} = 0.5 AU. In contrast to the previous models, here we first evolve the system without the gravitational backreaction of the disk, i.e. the planet only feels the gravitational force of the two stars. In this way, we calculate the approximate equilibrium state of the combined disk+planet system, to avoid initial transients. Once this equilibrium is reached, we include the gravitational backreaction of the disk to allow for possible orbital evolution of the planet; Fig. 20 shows the results. The evolution of the planet’s semimajor axis is shown in the top panel and the bottom panel shows its eccentricity. The orange curve refers to the previous model where the planet was started at a_{0} = 1.0 AU under the full action of the disk. As shown in the figure, while the planet’s semimajor axis is continuously shrinking, its eccentricity is slowly increasing. For the model started at a_{0} = 0.5, the planet remains first in its orbit while it does not feel the gravity of the disk (purple curve). At this state, the eccentricity of the planet oscillates between 0.0 and 0.07 with an average of about 0.036. We note that these variations are solely due to the effect of the binary star and are always present for an orbiting planet. Once the backreaction of the disk is included, the orbital elements of the planet begin to slowly vary (light blue curve). The planet begins to move slightly outward while its eccentricity becomes lower because of the damping effect of the disk. For illustration purposes, we have superimposed the black dashed line which in the first phase shows the equilibrium eccentricity, and after the release of the planet, its exponential decline with a timescale of 600 yr. After its initial outward migration, the planet moves slowly inward and finally reaches a stable orbit exterior to the 5:1 resonance. From the results shown in Fig. 20, we infer that the final position of the planet is very close to a_{p} = 0.5 AU, just slightly outside its observed orbit.
5. Summary and discussion
Using 2D hydrodynamical simulations, we studied the evolution of a planet embedded in a circumbinary disk. We carried out our simulations for the planetary system of Kepler38 as our test binary. In our analysis, we adopted a twostep approach in which we first studied the structure of a circumbinary disk without a planet, and then included a planet in the disk and followed its evolution. To connect to previous studies by other authors, we began by considering a locally isothermal equation of state in our disk models. However, in order to extended the analysis to more realistic systems, we considered radiative models that include viscous dissipation, vertical cooling, and radiative diffusion. In the following, we briefly present our most important results.

a)
Inner boundary. To enable mass accretion into the inner cavityof the disk, we used an open boundary condition for the inner edgeof the disk. This boundary condition allowed for free inflow ofmaterial onto the binary star. At the outer edge of the disk, we considered a fixeddensity boundary condition which, combinedwith the open boundary at the inner edge, resulted in disk modelswith constant massaccretion rates. Compared to the disks withclosed inner boundaries or those with slow viscous inflows, asused for example by Pierens & Nelson (2013), ouropen boundary condition lead to the disk models with reducedand broader maximum densities and lower disk eccentricities. Initially, we considered the location of the inner edge of thecomputational grid, (i.e., r_{min}), to be at a distance approximately 1.7 times the binary semimajor axis. Following Pierens & Nelson (2013), who showed that the location of the disk’s inner edge can influence the structure of the inner part of the disk, we carried out simulations for different values of r_{min} and determined that a value of r_{min} ≈ 1.5a_{bin} (equal to 0.22 AU for the Kepler38 system) would be sufficient to capture all necessary physics. A disk model with r_{min} = 0.20 AU also resulted in the same gap profile. We note that from a numerical standpoint, the highest possible value of r_{min} should be considered because lower values of this quantity will lead to significantly reduced time steps in which case the FARGOalgorithm is no longer applicable. Therefore, in systems where the (expected) inner cavity is larger, a higher value of r_{min} is certainly preferable.

b)
Planets in isothermal disks. We studied locally isothermal disks with a viscosity parameter α = 0.01 and vertical scale height of H/r = 0.05. Before embedding the planet in the disk, we allowed the disk to reach a quasiequilibrium state, which for our standard disk model took approximately 2000 yr (~39 000 binary orbits). We then started the planet at various distances from the binary and showed that as expected for a dissipative viscous disk, the planet stops its radial migration and settles in the same orbit near the inner edge of the disk regardless of its starting point. For the case with r_{min} = 0.25 AU, we found that the planet was captured in a 5:1 resonance with the binary and remained in that orbit until the end of the simulation. The stability of the resonant configuration may be caused by the damping action of the disk, because circumbinary resonances such as 5:1 are known to be unstable for planetary orbits (Dvorak 1986). Hence, we expect unstable behaviour in this case once the disk has been dissipated, possibly leading to encounters with the binary and departure from a closed orbit. For a disk with r_{min} = 0.22 AU (a more realistic case), the inner cavity was larger and as a result, the planet ceased its migration slightly outside the 5:1 resonance at a ≈ 0.47 AU, which is in a very good agreement with the observed location of circumbinary planet Kepler38b. The mean eccentricity of this planet, e ≈ 0.03, also agrees with that of Kepler38b. Our results indicated that in the case of isothermal disks, the disk mass does not influence the final position of the planet because the zerotorque position only depends on the density profile and not the absolute value, as the density scales out of the equations.

c)
Planets in radiative disks. In a final set of simulations, we developed more realistic disk models by including viscous dissipation, radiative cooling, and diffusion. We found that when we used the same absolute density as in the standard isothermal disk (≈2100 g/cm^{2} at 2 AU), the viscous heating made the disk hotter than in the isothermal cases. This additional heating caused the disk to become more eccentric and the inner cavity to become smaller. As a result, the embedded planet was able to cross the 5:1 resonance and reach its final position at a ≈ 0.415 AU (slightly outside the instability region) with a low eccentricity of e ≈ 0.03. This orbit will again be unstable upon dissipation of the disk. A reduction of the disk density combined with a lower viscosity leads to a wider inner cavity with a lower disk temperature. In this case, the disk mass becomes very low and as a result, the migration time becomes too long to follow the planet’s evolution to the end, as found by Pierens & Nelson (2013) as well. To save computational effort, we started a model with a planet initially much farther in, and obtained an equilibrium state where the planet stopped its radial migration just inside of a ≈ 0.5 AU, again very close to the observed location.
In agreement with previous simulations (Pierens & Nelson 2007, 2013), our study indicates that planets migrate inward in circumbinary disks. The final position of a planet depends on the size of the inner cavity. The results of isothermal runs suggest that the disk height has to be around H/r = 0.05, however, no other information on the disk mass can be inferred from these simulations. In fully radiative simulations, the mass of the disk
will determine the temperature and hence H/r, as in these systems the cooling depends on the vertical column density. From our simulations we infer that the required surface density is probably around a quarter of the minimum solar mass nebula.
To improve the models, the stellar irradiation from the two stars of the binary can be included. The combined irradiation will heat up the inner parts of the disk causing the disk’s thickness to increase which itself results in an increase in the disk’s eccentricity. What effects this may have on the final outcome of the simulations remains to be seen. A further improvement can be obtained with full 3D simulations including radiative transport. However, considering the necessarily long evolution time where tenthousands of binary orbits have to be followed, 3D simulations are certainly still challenging.
Acknowledgments
We thank Arnaud Pierens for stimulating discussions. N.H. acknowledges support from the NASA ADAP grant NNX13AF20G, and the Alexander von Humboldt Foundation. N.H. would also like to thank the Computational Physics Group at the University of Tübingen for their kind hospitality during the course of this project.
References
 Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge: Cambridge University Press) [Google Scholar]
 Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651 [NASA ADS] [CrossRef] [Google Scholar]
 Artymowicz, P., & Lubow, S. H. 1996, ApJ, 467, L77 [NASA ADS] [CrossRef] [Google Scholar]
 D’Angelo, G., Henning, T., & Kley, W. 2003, ApJ, 599, 548 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [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]
 Dunhill, A. C., & Alexander, R. D. 2013, MNRAS, 435, 2328 [NASA ADS] [CrossRef] [Google Scholar]
 Dvorak, R. 1986, A&A, 167, 379 [NASA ADS] [Google Scholar]
 Günther, R., & Kley, W. 2002, A&A, 387, 550 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
 Kley, W., & Crida, A. 2008, A&A, 487, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., & Dirksen, G. 2006, A&A, 447, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., & Nelson, R. P. 2012, ARA&A, 50, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., Papaloizou, J. C. B., & Ogilvie, G. I. 2008, A&A, 487, 671 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1985, in Protostars and Planets II, eds. D. C. Black, & M. S. Matthews, 981 [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]
 Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Meschiari, S. 2012, ApJ, 761, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, T. W. A., & Kley, W. 2012, A&A, 539, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, T. W. A., & Kley, W. 2013, A&A, 560, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge: Cambridge University Press) [Google Scholar]
 Nelson, R. P. 2003, MNRAS, 345, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, R. P., & Papaloizou, J. C. B. 2002, MNRAS, 333, L26 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, ApJ, 758, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, Science, 337, 1511 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Paardekooper, S.J., Leinhardt, Z. M., Thébault, P., & Baruteau, C. 2012, ApJ, 754, L16 [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. 2008a, A&A, 478, 939 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pierens, A., & Nelson, R. P. 2008b, 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]
 Ward, W. R. 1988, Icarus, 73, 330 [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]
Online material
Movie 1: The central part of the disk. Locally isothermal Model disk model without a planet, the standard model in the paper, see Fig. 4. (Access here)
Movie 2: Moving planet around binary. Isothermal model: locally isothermal Model disk model with an embedded planet, corresponding to the final phase where the planet is in the 5:1 resonance with the binary, see Fig. 9. (Access here)
Movie 3: Moving planet around binary. Radiative model: radiative disk model with an embedded planet, corresponding to the final phase, of the blue/green curve in Fig. 17. (Access here)
All Tables
Binary parameter and observed planetary parameter of the Kepler38 system used in the simulations.
All Figures
Fig. 1 Total mass of the disk (without the planet) in our locally isothermal standard model (model 1). After a few thousand years, the disk reaches equilibrium and a nearly constant mass is established. 

Open with DEXTER  
In the text 
Fig. 2 Rate of mass accretion through the disk in model 1 measured at two different distances. The mean value, −8.68 × 10^{7}, is also plotted for r = 0.68 AU. As the central binary accretes material, the mass flow is directed inward and becomes negative. In the main text, however, we report the accretion rate in the form of positive absolute values. 

Open with DEXTER  
In the text 
Fig. 3 Azimuthally averaged surface density and disk eccentricity for model 1. Displayed are the profiles at different evolutionary times (in years); the red curve denotes the initial setup. 

Open with DEXTER  
In the text 
Fig. 4 Twodimensional density structure of an isothermal disk (model 1) around the central binary star. The graph shows a local view around the binary. The computational grid extends from r = 0.25 AU to r = 2.0 AU. The white inner region lies inside the computational domain and is not covered by the grid. The positions of the primary and secondary are indicated by the gray dots. The Roche lobe of the secondary is also shown. (Online movie) 

Open with DEXTER  
In the text 
Fig. 5 Azimuthally averaged surface density and disk eccentricity profiles for two models with different inner boundary conditions. The red curve corresponds to the standard model with an open inner boundary and the blue curve represents a system with a reflecting inner boundary (model 2). In both models the disk is at equilibrium. 

Open with DEXTER  
In the text 
Fig. 6 Azimuthally averaged surface density and eccentricity profiles for three disk models with different locations of the inner boundary, r_{rmin}. The red curve shows the standard model (model 1) with r_{min} = 0.25 AU, the blue curve is for a disk with r_{min} = 0.22 AU, and the green curve represents the disk model with r_{min} = 0.20 AU (model 3). In all three models the disk is at equilibrium. 

Open with DEXTER  
In the text 
Fig. 7 Evolution of the semimajor axis of the planet in an isothermal disk (model 1). In two simulations, the planet is started in a circular orbit at a distance of a_{0} = 1.75 AU from the center of mass of the central binary. In the simulation shown in red, the planet is included in the disk from the very beginning of the simulation and it evolves as the disk evolves starting from its initial density profile. In the simulation shown in blue, the disk is first relaxed to the equilibrium and then the planet is embedded in it. In the simulation shown in green, the planet is started at a_{0} = 1.0 AU in a relaxed disk. The inset shows the results near the end of the simulation for the first model. The dashed black line in the inset corresponds to the mean value of the planet’s semimajor axis at a_{p} = 0.436 AU. 

Open with DEXTER  
In the text 
Fig. 8 Evolution of the eccentricity of the planet in an isothermal disk shown in red in Fig. 7 (model 1). The planet was embedded in the disk prior to the start of the disk evolution at a semimajor of a_{0} = 1.75 AU. 

Open with DEXTER  
In the text 
Fig. 9 Twodimensional snapshots of the evolution of the disk surface density with an embedded planet. The snapshots correspond to the times 449, 729, 2807, and 7414 years. The planet, shown as the gray dot with the Roche radius, was embedded in the disk at an initial distance of r = 1.75 AU. The dashed gray line around the central binary shows the approximate boundary of stability for planetary orbits. A movie of the simulation shown here can be found online. 

Open with DEXTER  
In the text 
Fig. 10 Averaged radial surface density of a disk with an embedded planet. The colors correspond to the times of the snapshots shown in Fig. 9. The big circle in each graph represents the radial position of the planet at the time corresponding to the graph. For illustrative purposes, we have moved the circles close to their corresponding curves. 

Open with DEXTER  
In the text 
Fig. 11 Evolution of the period ratio of the planet and binary for the model where the planet is embedded in the initial disk at a_{0} = 1.75 AU. The inset shows the period ratio for the final stage of its evolution. The dashed line in the inset corresponds to the mean value of P_{planet}/P_{bin} = 5.22. 

Open with DEXTER  
In the text 
Fig. 12 Resonant angles Φ_{1} and Φ_{2} (see Eqs. (6) and (7)) during the capture phase in a 5:1 meanmotion resonance. 

Open with DEXTER  
In the text 
Fig. 13 Evolution of the semimajor axis of a planet embedded in an isothermal disk for different values of the disk mass. The graph in red corresponds to the simulation of the standard disk model with the planet initially starting at a_{0} = 1.0 AU. The graph in blue shows a similar simulation but with half the disk mass. The graph in green corresponds to a simulation with a quarter of the initial disk mass that is continued from the model in blue at t = 4300. 

Open with DEXTER  
In the text 
Fig. 14 Evolution of the semimajor axis of a planet embedded in an isothermal disk for different locations of the disk’s inner boundary. The graph in red corresponds to the standard model with an open inner boundary at 0.25 AU. The graph in blue represents a model with a closed inner boundary at 0.25 AU, and the graph in green is for a disk with an open inner boundary at 0.22 AU. 

Open with DEXTER  
In the text 
Fig. 15 Evolution of the semimajor axis of the binary in a system with an isothermal disk for different values of the disk mass. The graph uses the same models and colors as in Fig. 13 for the evolution of the planets. 

Open with DEXTER  
In the text 
Fig. 16 Surface density (in g/cm^{2}, top panel), temperature (in K, middle panel), and eccentricity (bottom panel) of a radiative disk in the initial model (red), an isothermal case with H/r = 0.05 (model 1, blue), and a quasiequilibrium state (green). 

Open with DEXTER  
In the text 
Fig. 17 Evolution of the semimajor axis and eccentricity of an embedded planet in a radiative disk. The planet was started at a_{0} = 1.0 AU (red curve). In the two other models shown in green and blue, simulations were continued at specified times with slightly different numerical parameters (see text for details). (Online movie) 

Open with DEXTER  
In the text 
Fig. 18 Surface density (in g/cm^{2}, top panel), temperature (in K, middle panel), and eccentricity (bottom panel) of two radiative disks with different masses. The green curve corresponds to the same disk model as the green curve in Fig. 16 with the standard disk mass and viscosity (α = 0.01). The orange curve shows a disk model with a mass equal to a quarter of the mass of our standard disk model and a lower disk viscosity of α = 0.004. 

Open with DEXTER  
In the text 
Fig. 19 Evolution of the semimajor axis of an embedded planet in a radiative disk with different masses. The red and blue curves correspond to the standard disk mass model (as shown in Fig. 17) and the orange curve corresponds to a model with a quarter disk mass. 

Open with DEXTER  
In the text 
Fig. 20 Evolution of the semimajor axis and eccentricity of a planet embedded in a radiative disk with a mass lower than that of our standard disk model. The orange curve corresponds to the same model as the orange curve in Fig. 19 with the planet started at a_{0} = 1.0 AU. The purple curve corresponds to a system in which the planet was embedded in the disk at a_{0} = 0.5 AU and moved only under the influence of the gravity of the binary. The light blue curve shows the same planet as in the purple curve, where at time t = 1700 the effect of the disk was also taken into account. The dashed line in the bottom panel shows the evolution of the planet’s mean eccentricity. The first equilibrium is reached at e ≈ 0.036, and under the action of the disk, this value of the eccentricity exponentially decreases to a new mean value of e ≈ 0.026 in a timespan of 600 years. 

Open with DEXTER  
In the text 