The formation of planets in circumbinary disks
We examine the formation of planets around binary stars in light of the recently discovered systems Kepler 16, 34 and 35. We conduct hydrodynamical simulations of self gravitating disks around binary systems. The selected binary and disk parameters are chosen consistent with observed systems. The disks are evolved until they settle in a quasi-equilibrium and the resulting systems are compared with the parameters of Kepler 16, 34 and 35. We find a close correspondence of the peak density at the inner disk gap and the orbit of the observed planets. We conclude, based on our simulations, that the orbits of the observed Kepler planets are determined by the size of the inner disk gap which for these systems results from the binary driving. This mediates planet formation either through the density enhancement or through planetary trapping at the density gradient inversion in the inner disk. For all three systems the current eccentricity of the planetary orbit is less than the disk eccentricity in the simulations. This, together with the long term stability of the orbits argues against in situ formation (e.g. a direct collapse scenario of the material in the ring). Conducting additional simulations of systems with a wider range of parameters (taken from a survey of eclipsing binaries), we find that the planet semi-major axis and binary eccentricity in such a scenario should be tightly correlated providing an observational test of this formation mechanism.
keywords:planetary systems: protoplanetary discs methods: numerical
The recent discovery of a population of planets in orbit around binary stars (Doyle et al., 2011; Welsh et al., 2012; Orosz et al., 2012a, b) adds an extra layer of complexity to planet formation theories, which are already struggling to explain the wide variety of systems uncovered by the current explosion of planet discoveries (Baraffe et al., 2010). Around a single star the time required to grow a giant planet is expected to be at least comparable, and probably exceeds, the lifetime of the protoplanetary disk (Yin et al., 2002) and is subject to a number of bottlenecks. After the formation of the gaseous protoplanetary disk (Wyatt, 2008) giant planets form through a number of subsequent phases each of which is expected to last at least a million years. The first planetesimal seeds are thought to form from the dense dustlayer in the midplane of the disk, either through gravitational instabilities directly(Michikoshi et al., 2010), or by the fast collapse (Johansen et al., 2007) and/or slow sedimentation of aggregations that form in a turbulent gas velocity field(Cuzzi et al., 2008). The subsequent growth of these planetesimals to rocky proto planet is mediated by oligarchic accretion (Ida & Makino, 1993; Kokubo & Ida, 1998), followed by the accretion of a giant envelope by capturing gas from the remaining disk (Fortier et al., 2009). The growth of the rock and the accretion of a gaseous envelope are not necessarily the most time consuming phases, and each of these can be realized in a few Myr (Fortier et al., 2009).
The most difficult part in the planet formation process appears to be the earlier phase in which the disk develops a highly contrasted density structure (Michikoshi et al., 2010), that eventually collapses under its own gravity to form the first planetesimals. In a disk around a single star the gravitational instability grows slower than the time scale for the dust to settle in a thin disk (Yamoto & Sekiya, 2006). Forming planets in such an environment is a non-trivial endeavour, and it may be mediated by turbulence in the dust disk (Johansen et al., 2007; Cuzzi et al., 2008). When the inner object is a binary rather than a single star, the orbital motion excites perturbations in the disk. These perturbations may hinder the growth of planetesimals by increasing their orbital eccentricities and hence their relative impact velocities (e.g. Marzari et al., 2008; Paardekooper et al., 2012).
In addition, the core accretion theory suffers from the problem that proto-planetary cores are expected to migrate inwards before they have a chance to grow in size. This planet migration model would naturally explain the presence of planets in close orbits around their parent stars (Ward, 1997; Tanaka et al., 2002), but have difficulty to reproduce the observed wide orbits. On the other hand, the direct formation scenario for gas giants via the collapse of disc instabilities, can only explain a limited set of planets in remote orbits (Boley, 2009; Janson et al., 2012).
Recently the concept of ’planet traps’ was introduced to overcome the problems caused by (type I) planet migration: it turns out that in regions of positive (increasing density with radius) radial density slope the net effect of the corotation and Lindblad resonance torques flips sign (Masset et al., 2006). This leads to a (possible) pile up of proto-planetary cores at the radii of positive slope (Note that additionally, these regions can have an inward pressure force leading to super Keplerian gas motions. As a consequence the gas drag on planetesimals also flips sign here and these regions are expected to show enhanced density for this reason; Masset et al., 2006).
The details of this process remain uncertain. The principal problem is that the actual radius where planets are trapped is ill constrained for disks around single stars (Morbidelli et al., 2008). If the trapping happens at the transition between the inner turbulent region and the ’dead zone’ it is uncertain at what distance and at what time the process is operating. The distance for this process could be as close to the parent star as AU (Ilgner & Nelson, 2006), up to a distance of AU (Turner et al., 2007), and the effective distance may even move with time. This is not necessarily a problem for this mechanism, but it makes it difficult to identify whether or not the planetary trapping process has taken place.
Circumbinary protoplanetary discs on the other hand have a definite scale and will have a more precisely defined region where the density gradient is positive. Circumbinary discs therefore do not suffer from the aforementioned finetuning problem. The inner regions are cleared out by the binary torque to radii (Artymowicz & Lubow, 1994, 1996). The precise location of the inner boundary is a weak function of the system parameters (semi-major axis , eccentricity and mass ratio ).
Hence circumbinary systems provide a test bed for planet formation
theories. The recent detection of circumbinary planets around Kepler
16, 34 and 35
The formation of planets around binaries may be quite common, and their apparent lack can easily be the result of biases in the techniques used to search for planets. The majority of solar type stars are in binary systems Duquennoy & Mayor (1991); Ghez et al. (1993), and each may have formed with a circum binary disk. However, the lifetime of disks around close binaries may present a constraint for planet formation: (Kraus et al., 2012) found for the Taurus Auriga star forming region that of close binaries (semi-major axis AU ) have dispersed their disk by 1 Myr, whereas 80%–90% of wide binaries and single stars retain their disk for more than 2-3 Myr. Even so, 1/3 of close binaries in the Kraus et al. study retain a disk for up to 10 Myr, so circumbinary planet formation does not seem to be precluded. If survival of these disks allows for the formation of circumbinary planets in a seizable fraction of the binary population, which seems to be confirmed by the detection of Kepler 16, 34 and 35 (Welsh et al., 2012), the prospects for detection by transit events are best in eclipsing binary systems (e.g. Devor et al., 2008) (assuming that circumbinary disks preferentially align with the binary, which is the case for the Kepler systems.) .
In this manuscript we study the formation of perturbations in hydrodynamic circumbinary disks and study the evolution of the density structure. We aim our study at the recently discovered systems Kepler 16, 34 and 35, with the aim of making a detailed comparisons with the observed planetary parameters and assess whether these provide constraints on the formation channel of these systems.
We simulate the structural evolution of a gaseous circumbinary disk
around a binary star system. The numerical framework consists of a
variety of ingredients brought together in the Astrophysics Multipurpose
Software Environment (AMUSE) (Portegies Zwart et al., 2009, 2012).
The AMUSE package combines well tested simulation codes into a software suite which can be used to perform individual tasks, or reassemble the parts into a new application that combines a wide variety of solvers. The interfaces of codes within a common domain are designed to be as homogeneous as possible. The AMUSE application consists of a user script, written in python, that controls the community modules. The user script specifies the initial conditions, manages the calling sequence and data flow between the community modules, controls the runtime error handling, checks for energy conservation and other runtime diagnostics and performs a primary analysis on the raw simulation data. In the AMUSE philosophy we use Python as a glue language to bind the community modules together. The relatively low speed of this high-level language is not an issue, because most of the work is done in the community codes.
For the simulations presented in this manuscript we co-evolved binary systems and a gaseous disc. The binary solver employed was a Kepler solver in universal variables (e.g. Bate et al., 1971). The gaseous disk is simulated using a self-gravitating smooth particle hydrodynamics solver (Pelupessy, 2005). The mutual gravitational interaction between the gas particles and the stars are implemented using the BRIDGE solver (Fujii et al., 2007). BRIDGE provides a symplectic mapping for gravitational evolution in cases where the dynamics of a system can be split into two (or more) distinct regimes. In the application presented here the internal dynamics of the binary evolves on a relatively short timescale compared to the dynamics of the circumbinary disk. We adopt the bridge scheme to couple the gravitational dynamics of the inner binary with the hydrodynamics and self gravity of the circumbinary disk. The actual mutual forces are calculated by direct summation, the self gravity of the disc is calculated using a Barnes-Hut tree (Barnes & Hut, 1986). The gas is evolved with an isothermal equation of state. We use the standard SPH viscosity formulation (Monaghan & Gingold, 1983) with and resulting in a Reynolds number of . Additionally, the stars can act as sink particles to the gas. If a gas particles passes closer than 0.05 AU by a star it is accreted onto that star, transferring its mass and momentum to the star.
3.1 Initial conditions
The initial conditions of the models consist of a binary star embedded in a circumbinary disk. For the parameters of the binary star we adopt the observed stellar masses, semi-major axis and eccentricity of Kepler 16, 34 and 35 systems. In addition to this we also run models for 9 known eclipsing binary stars taken from the Devor et al. (2008) catalogue to examine a wider range of parameters. We will discuss the latter runs in Section 3.3.
The initial protoplanetary disk is set-up as an axi-symmetric disc with mass constructed to be in equilibrium with a central mass . The disk mass adopted is . The disk is assumed to be aligned with the binary orbit, which is plausible given the alignment of the binary and planetary orbits (within ). The discs have an inner boundary at and an outer boundary (0.5 AU and 8 AU for the Kepler 16 model). The inner boundaries were chosen such that a limited number of gas particles are accreted onto the stars during the initial stages of the simulation. In this way the results do not depend on spurious accretion, as any smaller inner boundary would see the additional mass quickly accreted onto the stars. The density profile of the disk is , with and temperature profile , normalized such that the Toomre Q parameter is equal to at the disk edge, resulting in K at 1 AU (for the Kepler 16 model). The parameters of the binaries and discs for the Kepler 16, 34 and 35 systems are listed in Table 1.
The Kepler 16, 34 and 35 systems turn out to sample different combinations of mass ratio q and eccentricity , where Kepler 34 and 35 have and Kepler 16 has (see Table 1). Kepler 16 and 35 have similar low while Kepler 34 has a high eccentricity .
We run the simulations for 1000 binary orbits. A snapshot of the resulting gas distribution after approximately 300 orbits for the Kepler 16, 34 and 35 runs are plotted in figure 1. At this point the distribution is in quasi-stationary equilibrium where the distribution of gas changes little qualitatively afterwards. Compared with the initial distribution of gas, the inner disc gap has expanded, and an eccentric overdense ring has formed (in agreement with Artymowicz & Lubow, 1994, 1996). The gas disk has been pumped into an eccentric orbit by high order Lindblad resonances (Pierens & Nelson, 2007). The interaction is highly non-linear though and the inside of the gas flow exhibits strong periodic tidal streams. Close inspection of the time lapses of the simulation shows that the onset of these tidal streams coincides with the start of the eccentricity pumping. The action of the resonances excites strong waves in the gas distribution. The azimuthally averaged density plots in Figure 1 (right panels) are for this reason time averaged over 10 snapshots (spanning the 100 last orbits). For all three systems there is a close correspondence of the location of the dense ring and the planet.
In order to test the robustness of our main conclusions we varied the parameters of the Kepler 16 run to test the effects of (small) disk inclinations (), variations in Toomre Q parameter (), () and disk mass fraction (), as well as numerical resolution, running at and . The results presented here are insensitive to these parameter variations.
Figure 2 shows the distribution of the (instantaneous) semi major axis and eccentricity of the disk material for the three runs. As can be seen the gas follows a relatively narrow relation in the plane where the material in the inner disk follows increasingly eccentric orbits, while the outer disc material follows circular orbits. In these figures the points with bars indicate the position of the observed planets in the diagram. The bars do not indicate the uncertainty in the observations (these are very small) but indicate the secular variation of the orbital elements in time, derived from long-term three body integrations (Doyle et al., 2011; Welsh et al., 2012). For Kepler 16, most of the variation consists of slow (on a timescale of yr) oscillations around the forced eccentricity due to secular pertubations from the binary, while for Kepler 34 and 35 most of the variation in the eccentricity is due to fast oscillations on the orbital timescale (as expected for equal mass binaries; Moriwaki & Nakagawa, 2004). For all models the eccentricities in the gas disk is greater than the observed planetary orbit even accounting for the variations in the orbital elements, although for Kepler 34 the disk material falls only just outside the range of the long-term integrations. Although we only have three systems, and the eccentricities are poorly matched, the system with the larger disc eccentricity coincides with the larger planetary ecccentricities.
An interesting question is whether any planets forming in the disk would be stable (if they inherited the orbital elements from the disk). In order to investigate this we have run a grid of three body models for each of the Kepler binaries with the central binary and semi-major axis and choosing a grid of 40x40 models. For the simulations we employed the Huayno solver (Pelupessy et al., 2012) in the AMUSE package, using its order 10 shared variable timestep (SHARED10) integrator (Sofroniou & Spaletta, 2005). We ran models at each starting from and going to lower until a stable model was encountered (a model was deemed unstable if a deviation of in was encountered - models were run for 1 Myr). Considering coplanar orbits the most important additional parameter is the angle between binary and planetary pericenter. For a given there is a small range in where a model can be stable or unstable depending on the pericenter angle. The results are summarized in figure 2, where models which are always unstable are marked with stars, while models which can be stable depending on pericenter angle are marked with plusses or crosses. In the case of Kepler 16, varying the pericenter angle shows the biggest change in the unstable , because of the larger mass ratio between the binary components. Note that all three observed planets lie within the stable region. Also, the disc material at the radius of the planet is in the unstable or partially stable region for all the models. This makes a scenario whereby the disc material is converted in situ into a planet (inhereting both semi-major axis and eccentricity, assuming that the dissolution of the gas disk does not alter the orbits of the remaining planetesimal material) followed by slow evolution of the orbital elements unlikely, unless inward radial migration has occurred.
3.3 Additional models
In order to examine a wider variety of parameters we run in the same way a set of models varying the mass and eccentricity of the primary stars, which we chose from a catalog of 773 eclipsing binaries (Devor et al., 2008). The systems were selected from the Devor et al. (2008) survey on the basis of their orbital period ( days) and eccentricity (). This ensures that we select systems that are unlikely to have undergone tidal interactions, in which case they are unsuitable to our analysis. This results in a set consisting of 9 systems (note this does not exclude the possibility of planets around the remaining 764 systems). An overview of the initial conditions of this “survey run” is given in Table 2.
The resulting gas disks from these runs were analyzed in the same way as in Sect. 3.2, i.e. we calculate the semi-major axis and eccentricity of the disk material. In Figure 3 we plot the resulting density distributions as a function of semi-major axis. If planets commonly form around close binaries we could expect them to form preferentially close to the peak of the density distribution with an eccentricity smaller than that of the disk material, as they have been seen for the hitherto detected circumbinary Kepler systems. If this is indeed the case we can derive the expected location of any planets in these systems from the simulations. We present in Table 2 fiducial planetary elements and . The eccentricity could lie anywhere between 0 and , whereas the semi-major axis would be expected to lie within of the given.
In Figure 4 we show a number of scatter plots of binary
parameters (mass ratio and eccentricity ) versus
‘planetary’ orbit semi-major axis and eccentricity
. For comparison we include the observed planets
Kepler 16, 34, 35, 38 and 47
which is expected in case circumbinary planets form preferentially at the peak density in the disk, should be testable as a more statistically significant number of planets around binary systems are found, especially as so far only one system with binary eccentricity is found.
The simulations presented here favor the hypothesis that the structure of a circumbinary disc and the planetary orbit of the observed systems Kepler 16,34 and 35 are closely related. The agreement in the correlation between binary eccentricity and planet semi-major axis shown in Figure 4 of the observed and simulated systems favour this simple scenario.
In principle there are three possibilities to account for such a correspondence: 1) The planet has formed in situ in the overdense region, 2) the planet has formed further out in the disk and migrated inward while the circumbinary disc was still present until it encountered the region of positive density slope at the gap ( in this case the the planet occupies an orbit close to the peak density just because the surface density rises quickly at the inner edge), and 3) the planet has formed further out in the disc and migrated by secular evolution (e.g. scattering by planetesimals) and stalled at the radius previously occupied by the inner disk edge for reasons that are not directly related to the previously existing gas disk.
Upon naive inspection of figure 1 alternative (1) seems preferable. However in this case the planet eccentricities are smaller than expected from figure 2 and furthermore planetesimals originating from the gas disk at the current semi-major axis of the planet would be in marginally stable orbits. This problem could be remedied by a mechanism whereby the eccentricity of a planet would decrease (e.g. by forming multiple proto-planets and having them collide). Given the fact that at the radii of the planets orbits of modestly higher eccentricities are already unstable, it seems that a modest amount of radial migration (i.e. moving to the left in the diagram of fig. 2) is preferred.
Scenario 3) - migration after the disc has cleared - suffers from the problem of fine-tuning the resulting orbit to match the location of peak density at the inner disc edge. In this respect 2) provides a more natural migration scenario. Additionally, gravitational scattering by larger planets (which may be supported by the detection of a binary system with multiple planets, Kepler 47) is unlikely given the fact that the eccentricities of most of the Kepler circumbinary planets are quite low. These planets are too far out for tidal circularisation (such as in the case of hot jupiters).
Of course the above conclusions are based on a simplified model where we have ignored the subsequent phases of planet formation (the build up of planetesimals and their growth). We argue, based on the close comparison, that these may not be essential in determining the planetary orbit characteristics for the special case of circumbinary planet formation. If indeed this is the case a more careful consideration of the latter phases in the disk evolution where the remnant material is cleared out while the planet emerges from the disc may still be needed to show that the planets end up in their long-term stable configuration.
We investigated the formation mechanism of the recently discovered population of planets orbiting both components of a binary by studying hydrodynamic simulations of circumbinary disks. We chose binary parameters according to the observed systems Kepler 16, 34 and 35 and an initial disk with canonical initial conditions. The disks were evolved until they settled in quasi equilibrium and the resulting systems were compared with the observed systems. We also run three-body simulations to investigate the long term stability of planetary orbits around the binaries. We performed 9 additional simulations for which the initial conditions are taken from a survey of eclipsing binaries. In Tab. 2 we present, for these observed binary systems, the predicted semi-major axis and eccentricity of the circum binary planet. We further conclude that:
Planets observed are found inside the theoretically expected overdensity at the inner edge of the circumbinary disk,
The material in the circumbinary disks itself settles in a narrow range of eccentricities,
The eccentricity ranges of the observed planets are smaller than that of the disk material, with the possible exception of the Kepler 34 model.
A relatively tight positive relation between planet semi-major axis and binary eccentricity is expected if planets form preferentially at the density peak just outside the inner edge of the circum binary disk.
These results suggest that planet formation in these systems and therefore the orbital parameters of these planets, are determined for a large part by the binary driving of the proto-planetary gas disk. In addition it seems necessary that the systems we have modelled in detail, Kepler 16,34,35, have experienced at least some planetary migration. Formation in the outer parts of the disk followed by migration inward is possible, but the naive expectation that the planet would be found inside the disk gap is not entirely fulfilled. This formation mechanism would need some fine-tuning but not as severely as for planetary formation models around single stars (in so far as these invoke planetary trapping). A direct relation between disk material, e.g. through gravitational collapse, is probably not the formation route, as this would imply that the planets would form with an eccentricity close to the material of the disk. If this were the case planets with this eccentricity could in principle also form, since most of the high eccentricity disk material is well within the region of stability. Planetesimal formation from the ring material followed by aggregation is not excluded. In this scenario the eccentricity of the planet is lowered through repeated collisions. In this case the close correspondence between the gas ring in our simulations and the planets is not accidental.
Acknowledgements This work was supported by the Netherlands Research Council NWO (grants #643.200.503, #639.073.803 and #614.061.608) and by the Netherlands Research School for Astronomy (NOVA).
- during revision of the manuscript additional circumbinary planet systems, Kepler 38 and 47, where published.
- The framework and the source codes of the scripts which were used to run the simulations presented here can be downloaded from www.amusecode.org.
- Note Kepler 47 is actually a two planet system, we consider the parameters of the inner planet.
- Artymowicz P., Lubow S. H., 1994, \apj, 421, 651
- Artymowicz P., Lubow S. H., 1996, \apjl, 467, L77
- Baraffe I., Chabrier G., Barman T., 2010, Reports on Progress in Physics, 73, 016901
- Barnes J., Hut P., 1986, \nat, 324, 446
- Bate R. R., Mueller D. D., White J. E., 1971, Fundamentals of astrodynamics.
- Boley A. C., 2009, \apjl, 695, L53
- Cuzzi J. N., Hogan R. C., Shariff K., 2008, \apj, 687, 1432
- Devor J., Charbonneau D., O’Donovan F. T., Mandushev G., Torres G., 2008, \aj, 135, 850
- Doyle L. R., Carter J. A., Fabrycky D. C. e. a., 2011, Science, 333, 1602
- Duquennoy A., Mayor M., 1991, \aap, 248, 485
- Fortier A., Benvenuto O. G., Brunini A., 2009, \aap, 500, 1249
- Fujii M., Iwasawa M., Funato Y., Makino J., 2007, \pasj, 59, 1095
- Ghez A. M., Neugebauer G., Matthews K., 1993, \aj, 106, 2005
- Ida S., Makino J., 1993, \icarus, 106, 210
- Ilgner M., Nelson R. P., 2006, \aap, 445, 205
- Janson M., Bonavita M., Klahr H., Lafrenière D., 2012, \apj, 745, 4
- Johansen A., Oishi J. S., Mac Low M.-M., Klahr H., Henning T., Youdin A., 2007, \nat, 448, 1022
- Kokubo E., Ida S., 1998, \icarus, 131, 171
- Kraus A. L., Ireland M. J., Hillenbrand L. A., Martinache F., 2012, \apj, 745, 19
- Marzari F., Thébault P., Scholl H., 2008, \apj, 681, 1599
- Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006, \apj, 642, 478
- Michikoshi S., Kokubo E., Inutsuka S.-i., 2010, \apj, 719, 1021
- Monaghan J. J., Gingold R. A., 1983, Journal of Computational Physics, 52, 374
- Morbidelli A., Crida A., Masset F., Nelson R. P., 2008, \aap, 478, 929
- Moriwaki K., Nakagawa Y., 2004, \apj, 609, 1065
- Orosz J. A., Welsh W. F., Carter J. A. e. a., 2012, ArXiv e-prints
- Orosz J. A., Welsh W. F., Carter J. A. e. a., 2012, ArXiv e-prints
- Paardekooper S.-J., Leinhardt Z. M., Thébault P., Baruteau C., 2012, \apjl, 754, L16
- Pelupessy F. I., 2005, PhD thesis, Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands
- Pelupessy F. I., Jänes J., Portegies Zwart S., 2012, \na, 17, 711
- Pichardo B., Sparke L. S., Aguilar L. A., 2005, \mnras, 359, 521
- Pierens A., Nelson R. P., 2007, \aap, 472, 993
- Portegies Zwart S., McMillan S., Harfst S. e. a., 2009, \na, 14, 369
- Portegies Zwart S., McMillan S., Pelupessy I., van Elteren A., 2012, in Capuzzo-Dolcetta R., Limongi M., Tornambè A., eds, Advances in Computational Astrophysics: Methods, Tools, and Outcome Vol. 453 of Astronomical Society of the Pacific Conference Series, Multi-physics Simulations using a Hierarchical Interchangeable Software Interface. p. 317
- Sofroniou M., Spaletta G., 2005, Optimization Methods and Software, 20, 597
- Tanaka H., Takeuchi T., Ward W. R., 2002, \apj, 565, 1257
- Turner N. J., Sano T., Dziourkevitch N., 2007, \apj, 659, 729
- Ward W. R., 1997, \icarus, 126, 261
- Welsh W. F., Orosz J. A., Carter e. a., 2012, \nat, 481, 475
- Winn J. N., Albrecht S., Johnson J. A. e. a., 2011, \apjl, 741, L1
- Wyatt M. C., 2008, \araa, 46, 339
- Yamoto F., Sekiya M., 2006, \apjl, 646, L155
- Yin Q., Jacobsen S. B., Yamashita K., Blichert-Toft J., Télouk P., Albarède F., 2002, \nat, 418, 949