Astronomy Unit, Queen Mary, University of London, Mile End Rd, London, E1 4NS, UK

Abstract
Aims. We present the results of hydrodynamic simulations of the formation and subsequent orbital evolution of giant planets embedded in a circumbinary disc. The aim is to examine whether or not giant planets can be found to orbit stably in close binary systems.
Methods. We performed numerical simulations using a grid-based hydrodynamics code. We assume that a $20~M_\oplus$ core has migrated to the edge of the inner cavity formed by the binary where it remains trapped by corotation torques. This core is then allowed to accrete gas from the disc, and we study its orbital evolution as it grows in mass. For each of the two accretion time scales we considered, we performed three simulations. In two of the three simulations, we stopped the accretion onto the planet once its mass became characteristic of that of Saturn or Jupiter. In the remaining case, the planet accreted disc material freely in such a way that its mass became higher than Jupiter's.
Results. The simulations show different outcomes depending on the final mass $m_{\rm p}$ of the giant. For $m_{\rm p}=1~M_{\rm S}$ (where $M_{\rm S}$ is Saturn's mass), we find that the planet migrates inward through its interaction with the disc until its eccentricity becomes high enough to induce a torque reversal. The planet then migrates outward, and the system remains stable on long time scales. For $m_{\rm p} \ge 1~M_{\rm J}$ (where $M_{\rm J}$ is Jupiter's mass) we observed two different outcomes. In each case the planet enters the 4:1 resonance with the binary, and resonant interaction drives up the eccentricity of the planet until it undergoes a close encounter with the secondary star, leading to scattering. The result can either be ejection from the system or scattering out into the disc followed by a prolonged period of outward migration. These results suggest that circumbinary planets are more likely to be quite common in the Saturn-mass range. Jupiter-mass circumbinary planets are likely to be less common because of their less stable evolution, but if present are likely to orbit at large distances from the central binary.

Key words: accretion, accretion disks - planets and satellites: formation - stars: binaries: close - hydrodynamics - methods: numerical

Among the approximately 260 extrasolar planets known at the time of writing (e.g. http://exoplanet.eu/), about 30 reside in binary or multiple-star systems, with most of them orbiting one stellar component in so-called S-type orbits (Eggenberger et al. 2004; Mugrauer et al. 2007). The majority of binary stars hosting planets have orbital separation $a_{\rm b} \ge$ 100 AU. However, there are a few cases of short period binary systems with $a_{\rm b} \sim$ 20 AU like Gliese 86, $\gamma$ Cephei and HD 41004 in which planets have been discovered to orbit at 1-2 AU from the primary (Eggenberger et al. 2004; Mugrauer & Neuhauser 2005). As suggested by studies of the long-term stability of planets in binary systems (e.g. Holman & Wiegert 1999), close binaries with $a_{\rm b} \sim$ 1 AU can, in principle, harbour planets evolving on a P-type orbit which encircles the two components of the binary system. One such circumbinary planet with mass of $m_{\rm p}=2.5~M_{\rm J}$ has been detected orbiting at 23 AU from the radio pulsar binary PSR 1620-26. Another with mass of $m_{\rm p}=2.44~M_{\rm J}$ was found evolving around a system which consists of the star HD 202206 and its 17.4 $M_{\rm J}$ brown dwarf companion (Udry et al. 2002). Because short period binaries are often rejected from observational surveys, circumbinary planets have not yet been observed in binary systems composed of two main sequence stars.

Several circumbinary discs, however, have been detected around spectroscopic binaries like DQ Tau, AK Sco and GW Ori. In GG Tau, the circumbinary disc has been directly imaged and has revealed the presence of an inner disc cavity due to the tidal torques exerted by the central binary (Dutrey et al. 1994). The existence of these circumbinary discs combined with the fact that $\sim$ 50% of solar-type stars are members of binaries (Duquennoy & Mayor 1991) suggests that circumbinary planets may be common provided that planet formation can occur inside such discs.

To date, several theoretical studies focused on planet formation in close binary systems indicate that planetesimal accretion should be possible within circumbinary discs. Moriwaki & Nakagawa (2004) found that planetesimals can grow in gas-free circumbinary discs only in regions farther out that $\sim$ 13 AU from a binary with orbital separation $a_{\rm b}=1$ AU, eccentricity $e_{\rm b}=0.1$ and mass ratio $q_{\rm b}=0.2$ . The influence of gas drag was studied recently by Scholl et al. (2007) who showed that for a binary with the same parameters, the eccentricity damping provided by the disc can enable planetesimal accretion to occur in regions located below $\sim$ 4 AU from the central binary. The later stages of planet formation in which Earth-mass planets form by accumulation of embryos was investigated by Quintana & Lissauer (2006). These authors found that planetary systems similar to those around single stars can be formed around binaries, provided that the ratio of the binary apocentre distance to planetary orbit is $\le$ 0.2. In general binaries with larger maximum separations lead to planetary systems with fewer planets.

Recently, the evolution of Earth-mass bodies embedded in a circumbinary disc and undergoing type I migration because of disc torques (e.g. Ward 1997) was examined by Pierens & Nelson (2007, hereafter referred to as Paper I). In this work, it was found that the inward drift of a protoplanet can be stopped near the edge of the cavity formed by the binary. Such an effect arises because in this region, the gradient of the disc surface density is such that the planet experiences strong positive corotation torques which can eventually counterbalance the negative differential Lindblad torque, thereby leading to the halting of migration (Masset et al. 2006). In a subsequent paper, Pierens & Nelson (2008a) extended this work by investigating the issue of how multiple protoplanets interact with each other if they form at large distance from the binary and successively migrate toward the cavity edge. The simulations performed by Pierens & Nelson (2008a) of pairs of planets interacting with each other indicated different outcomes such as resonant trapping or orbital exchange, depending on the ratio between the masses of the planets. Interestingly, in simulations involving more than two planets, planetary growth resulting from scattering and collisions between protoplanets was found to be feasible. This implies that giant cores might be formed in circumbinary discs, resulting eventually in a gas giant planet orbiting near the cavity edge.

The orbital evolution of Jovian mass planets embedded in circumbinary discs was studied by Nelson (2003). This work showed that giant planets undergoing type II migration (e.g. Lin & Papaloizou 1993; Nelson 2000) are likely to enter the 4:1 resonance with the binary. In that case, the subsequent evolution of the giant can be twofold, depending on whether or not the resonance is stable. For systems in which the 4:1 resonance is stable, the planet remains near or at the resonance. However, it appears that there is a finite probability for the system to be unstable, in which case the planet can be ejected from the system due to close encounters with the binary.

In this paper, we extend the work of Nelson (2003) by studying the whole evolution of a circumbinary planet during its growth from a core into a gas giant. To address this issue, we consider a scenario in which a 20 $M_\oplus$ core initially trapped at the edge of the cavity can slowly accrete gas from the disc. As the planet grows, the onset of non-linear effects as well as gap formation can exclude gas material from the coorbital region, thereby cancelling the effects of corotation torques. This can subsequently lead to the planet migrating inward again, until it becomes trapped eventually in a mean motion resonance with the binary. Here, we present the results of hydrodynamic calculations aimed at simulating such a scenario. In particular, we want to examine how the final outcome of the system depends on the accretion rate onto the planet as well as on the final mass of the giant. Interestingly, the results of the simulations indicate that only Saturn-mass giant planets can evolve stably in circumbinary discs. Most of the calculations of embedded giants with masses of $m_{\rm p} \ge 1~M_{\rm J}$ resulted in close encounters between the planet and the binary, leading eventually to the planet being completely ejected from the system.

This paper is organized as follows. In Sect. 2, we describe the hydrodynamical model. The results of the simulations are discussed in Sect. 3. We finally summarise and present our conclusions in Sect. 4.

2 Hydrodynamical model

2.1 Numerical method

We consider a 2D disc model in which all physical quantities are vertically averaged and we work in polar coordinates $(r,\phi)$ with the origin located at the centre of mass of the binary. The equations governing the disc evolution as well as the equations of motion for the binary plus planet system can be found in Paper I. The equations for the disc are solved using the hydrocode Genesis, which has been tested extensively against other codes (De Val-Borro et al. 2006), and employs a second order numerical scheme based on the monotonic transport algorithm (Van Leer 1977). Included in this code is a fifth-order Runge-Kutta integrator (Press et al. 1992) used to compute the evolution of the planet and binary orbits.

As in Paper I, we use $N_{\rm r}=256$ radial grid cells uniformly distributed between $r_{\rm in}=0.5$ and $r_{\rm out}=6$ and $N_{\phi}=380$ azimuthal grid cells. We adopt also the same computational units in which the mass of the binary is $M_\star=1$ , the gravitational constant is G=1 and the radius r=2 in the computational domain corresponds to 5 AU. The unit of time is $\Omega^{-1}=\sqrt{GM_*/a_{\rm b}^3}$ , where $a_{\rm b}=0.4$ is the initial binary separation. In the following, we report our results in units of the initial orbital period of the binary $P=2\pi\Omega^{-1}$ .

In this paper, the planet is allowed to accrete gas from the disc. Accretion by the protoplanet can be modelled by removing at each time-step a fraction of the gas located inside the Roche lobe of the planet and then adding the corresponding amount of matter to the mass of the planet (e.g. Kley 1999; Nelson et al. 2000). Here, the removal rate is chosen such that the accretion time scale onto the planet is $t_{\rm acc}=t_{\rm dyn}$ , where $t_{\rm dyn}$ is the orbital period of the planet. This corresponds to the maximum rate at which the planet can accrete gas material (Kley 1999). In order to examine how the final outcome of the system depends on the accretion rate, we have also performed simulations with $t_{\rm acc}$ = $10~t_{\rm dyn}$ .

When calculating the gravitational force between the disc and planet we adopt a gravitational softening parameter b=0.6H, where H is the disc height at the planet location. Material contained within the planet Hill sphere does not contribute to the gravitational acceleration experienced by the planet.

2.2 Initial conditions

In Paper I, we presented the results of simulations of protoplanets with masses of $m_{\rm p}=$ 5, 10 and 20 $M_\oplus$ embedded in circumbinary discs. We found that in each case, the inward migration of the protoplanet is halted due to the action of corotation torques which operate strongly at the cavity edge. Here, we extend the model presented in this earlier paper and examine how such a trapped protoplanet evolves as it accretes gas from the disc and grows to become a giant planet. In order to investigate this issue, we restarted one of the simulations presented in Paper I for which $m_{\rm p}=20~M_\oplus$ at a point in time when the planet is trapped at the cavity edge. Therefore, in this new series of simulations, the protoplanet evolves initially on an orbit with $a_{\rm p}\sim 1.2$ and $e_{\rm p}\sim 0.02$ , which are the values for the semi-major axis and eccentricity that a 20 $M_\oplus$ body finally attains once its migration has been stopped (see Paper I).

The initial semi-major axis and eccentricity of the binary are $a_{\rm b}\sim 0.39$ and $e_{\rm b}\sim 0.08$ respectively. These values correspond to the ones that an initially circular binary with $a_{\rm b}=0.4$ and mass ratio $q_{\rm b}=0.1$ ultimately attains as it interacts with a circumbinary disc. In Paper I, we indeed showed that the evolution outcome of such a system is an equilibrium configuration for which the disc structure as well as the binary eccentricity remain unchanged. Interestingly, from the time this quasi-steady state is reached, we found that the apsidal lines of the disc and binary are almost aligned.

The disc model used in this work is the same as that of Paper I. Accordingly, the aspect ratio is constant and equal to H/r=0.05. The initial disc surface density profile is presented in Fig. 1. For numerical reasons (see Sect. 2.3), we use a low-density region from r=4 to r=6. From $r\sim 2$ to r=4, the surface density is $\Sigma(r)=\Sigma_0\;r^{-1/2}$ where $\Sigma_0$ is chosen such that the unperturbed disc would contain 0.01 $M_\odot$ inside 10 AU (we assume that the initial binary separation corresponds to 1 AU). We further note that such a density profile corresponds to the equilibrium density for an accretion disc with a constant kinematic viscosity (e.g. Gunther & Kley 2002). In the inner parts, the disc presents an inner cavity which is created self-consistently from simulations of binary-disc interactions (see Paper I). It is worthwhile to notice that here, the onset of non-linear effects due to the presence of a 20 $M_\oplus$ initially located at $r\sim 1.2$ can modify slightly the gap profile.

We model the disc turbulent viscosity using the ``alpha'' prescription for the effective kinematic viscosity $\nu=\alpha c_{\rm s} H$ (Shakura & Sunyaev 1973). In Paper I, we set $\alpha=10^{-4}$ because using larger values caused the binary separation to decrease too rapidly to allow an equilibrium configuration to be obtained (see Paper I for details). In this work however, we examine the evolution of giant planets undergoing type II migration. Because in that case the migration rate of the planet is controlled by the disc viscous evolution, we decided here to use a more realistic $\alpha$ value, which probably lies in the range 10^-3-10^-2 in real circumstellar dics. In order to obtain a disc model in which $\alpha=10^{-3}$ , the calculations were started with $\alpha$ increasing slightly from 10^-4 to the desired value in $\sim$ 1500 binary orbits. Although not sufficient for a new equilibrium configuration to be established, this gives a sufficient time to the binary plus disc system to readjust to the new $\alpha$ value.

$\begin{figure} \par\includegraphics[width=7.3cm,clip]{9453f1.eps} \end{figure}$	Figure 1: This figure shows the initial disc surface density profile. The initial position of the planet with mass of $m_{\rm p}=20~M_\oplus$ is represented by a black circle.
Open with DEXTER

2.3 Boundary conditions

In order to avoid any wave reflection at the outer edge of the computational domain, we impose a low-density region between r=4 and r=6 using a taper function.

At the inner boundary, we model the accretion onto the central star by setting the radial velocity in the inner ghost zones to $v_{\rm r}=\beta v_{\rm r}(r_{\rm in})$ , where $v_{\rm r}(r_{\rm in})=-3\nu/2r_{\rm in}$ is the gas drift velocity at the disc inner edge due to viscous evolution and where $\beta$ is a free parameter. Following Pierens & Nelson (2008b), we set $\beta=5$ in the simulations presented here.

3 Results

$\begin{figure} \par\includegraphics[width=7.3cm,clip]{9453f2.eps}\vspace*{3mm} \includegraphics[width=7.3cm,clip]{9453f3.eps} \end{figure}$

Figure 2: This figure shows the mass of the planet as a function of time for simulations in which $t_{\rm acc}=t_{\rm dyn}$ ( upper panel) and $t_{\rm acc}=10~t_{\rm dyn}$ ( lower panel). The dotted (resp. dashed) line corresponds to simulations in which the final mass of the planet is 1 $M_{\rm s}$ (resp. 1 $M_{\rm J}$ ). In simulations represented by the solid line, the final mass of the planet is not imposed.

Open with DEXTER

$\begin{figure} \par\includegraphics[width=7.6cm,clip]{9453f4.eps}\hspace*{4mm} \includegraphics[width=7.6cm,clip]{9453f5.eps} \end{figure}$

Figure 3: This figure shows the evolution of the system for models in which $m_{\rm p}=1~M_{\rm S}$ . Here, the left and right panels correspond to $t_{\rm acc}=t_{\rm dyn}$ and $t_{\rm acc}=10~t_{\rm dyn}$ respectively. Top: evolution of the semi-major axis of the planet. Middle: evolution of the eccentricities for the planet (black), binary (red) and disc (green). Bottom: evolution of the longitudes of pericentre for the planet (black), binary (red) and disc (green).

Open with DEXTER

For each value of the accretion rate onto the planet we consider, we have performed three simulations which differ in the final mass attained by the giant. In the first and second calculations, we prevent further growth of the planet once its mass becomes $m_{\rm p}=1~M_{\rm s}$ (where $M_{\rm s}$ is Saturn's mass) and $m_{\rm p}=1~M_{\rm J}$ respectively. In the last run however, the final mass of the planet is not restricted and we allow the latter to accrete gas freely from the disc. The mass of the planet as a function of time is presented in Fig. 2 for each run.

3.1 Evolution of Saturn-mass planets

For both values of $t_{\rm acc}$ , the semi-major axis evolution of a 20 $M_\oplus$ body which grows to become a Saturn-mass planet is displayed in the upper panel of Fig. 3. In the model with $t_{\rm acc}=t_{\rm dyn}$ , the planet mass reaches $m_{\rm p}=1~M_{\rm s}$ at $t\sim 300~P$ , and this mass is attained at $t\sim 2800~P$ in the model with $t_{\rm acc}=10~t_{\rm dyn}$ (see Fig. 2). The orbital evolution of the planet is, however, very similar in both cases and typically proceeds as follows. At the beginning of the simulation, gap formation due to the growth of the protoplanet can exclude material from the coorbital region, which results in the (negative) differential Lindblad torques being no longer couterbalanced by the (positive) corotation torques. The planet thus migrates inward again. As it migrates, the interaction with the central binary makes $e_{\rm p}$ increase slowly, as observed in the middle panel of Fig. 3 which shows the evolution of the eccentricities of the binary, disc and planet for both values of $t_{\rm acc}$ . We define the disc eccentricity $e_{\rm d}$ by:

$\begin{displaymath}% e_{\rm d}=\frac{\int_0^{2\pi}\int_{r_{\rm in}}^{r_{\rm max}... ...int_0^{2\pi}\int_{r_{\rm in}}^{r_{\rm max}}\Sigma\; {\rm d}S}; \end{displaymath}$

(1)

Once the planet eccentricity reaches $e_{\rm p}\sim 0.05$ , which occurs at $t\sim 500~P$ , the migration of the planet stops and then reverses. Previous work of Papaloizou & Larwood (2000) indicated that an eccentric low-mass protoplanet can experience net positive disc torques when $e_{\rm p}\sim 1.1 H/R$ , owing to the latter crossing resonances in the disc that do not overlap the orbit at low-eccentricities. As pointed out by these authors, this occurs because a planet evolving on a high eccentric orbit rotates more slowly than the disc at apocentre, which results in the outer disc exerting a positive torque on the planet. Although here the planet mass is in the Saturnian mass range, it appears that the observed migration reversal is caused by a similar phenomenon. In order to demonstrate that such an effect is at work here, we have plotted in Fig. 4 the evolution of both the disc torques and planet orbital position $r_{\rm p}$ over a few orbital periods of the planet. In agreement with Papaloizou & Larwood (2000), the disc torques exerted on the planet are positive (negative) when the latter is at apocentre (pericentre). Clearly, there is a slight imbalance between the positive and negative torques, due to the planet orbiting at a cavity edge, which combined with the fact that the planet spends more time at apocentre, favours a net positive torque on the planet and outward migration. We note that in Fig. 4, there is a slight phase shift between the two curves which corresponds to the time needed for the planet to create an inner (outer) wake at apocentre (pericentre).

$\begin{figure} \par\includegraphics[width=7.3cm,clip]{9453f6.eps} \end{figure}$	Figure 4: This figure shows, over a few planetary orbital periods, the evolution of the torques exerted by the disc on the planet (solid line) as well as the time evolution of the orbital the position of the planet $r_{\rm p}$ (dashed line).
Open with DEXTER

Whereas this reversed migration proceeds very slowly until $t\sim 4.5$ $\times$ 10⁴ P, subsequent evolution involves exponential growth of the planet migration rate, which is characteristic of an episode of runaway migration (Masset & Papaloizou 2003). Saturn-mass planets are known to be good candidates for such a migration regime in massive protoplanetary discs, due to their ability to create a coorbital mass deficit larger than their own mass. As noticed by Masset & Papaloizou (2003), outward runaway migration can eventually occur in discs with shallow surface density profiles, provided that the planet initially migrates with a significantly positive drift rate. Here, two effects contribute to make such a phenomenon possible. First, prior to the period of runaway migration, the planet migrates outward over $\sim$ 2 $\times$ 10⁴ binary orbits, which is clearly larger than the planet libration timescale. Second, the planet evolves in a region of strong positive surface density gradient such that the coorbital mass deficit increases as the planet migrates outward. This can be seen for example in Fig. 5 which shows, for the model in which $t_{\rm acc}=t_{\rm dyn}$ , a series of surface density plots at different times.

$\begin{figure} \par\includegraphics[width=5.2cm,clip]{9453f7.eps}\hspace{3mm} \... ...f9.eps}\hspace{3mm} \includegraphics[width=5.2cm,clip]{9453f10.eps}\end{figure}$	Figure 5: This figure shows, for the model in which $m_{\rm p}=1~M_{\rm S}$ and $t_{\rm acc}=t_{\rm dyn}$ , snapshots of the disc surface density at times shown above the plots. Here, the planet is represented by a white circle.
Open with DEXTER

Interestingly, this period of outward runaway migration breaks at $t\sim 8$ $\times$ 10⁴ P, which corresponds to the time when the planet passes through the 5:1 resonance with the binary. This can excite the planet eccentricity to such values that during the course of an orbit, the radial excursion of the planet can exceed its coorbital width, leading to the loss of the coorbital mass deficit. The fourth panel in Fig. 5 displays the disc surface density just after the 5:1 resonance crossing. Comparing this panel with the second one which corresponds to a time when the planet is about to undergo runaway migration, it is clear that the coorbital region is no longer depleted, which can prevent runaway migration being sustained. This result is in broad agreement with the calculations performed by Masset & Papaloizou (2003) which indicate a similar tendency for the runaway migration to not be maintained.

After the 5:1 resonance crossing, continuation of the runs indicates that the planets undergo slow outward migration for $\sim$ 5 $\times$ 10⁴ binary orbits and then migrate inward again. Due to the very long run times required by these simulations, the final fate of the planets remains uncertain. Nevertheless, it seems likely that continued inward migration will proceed until the eccentricities of the planets are high enough for the disc torques to become positive, resulting eventually in the planets migrating outward again. Thus we expect the long term evolution to consist of periods of inward followed by outward migration until disc dispersal, resulting in a stable Saturn-mass circumbinary planet.

The evolution of the longitudes of pericentre of the binary, disc and planet is depicted in the lower panel of Fig. 3. As in Paper I, we compute the disc longitude of pericentre according to the following definition:

$\begin{displaymath}% w_{\rm d}=\frac{\int_0^{2\pi}\int_{r_{\rm in}}^{r_{\rm max}... ...int_0^{2\pi}\int_{r_{\rm in}}^{r_{\rm max}}\Sigma\; {\rm d}S}; \end{displaymath}$

(2)

3.2 Evolution of Jupiter-mass planets

$\begin{figure} \par\includegraphics[width=7.6cm,clip]{9453f11.eps}\hspace*{4mm} \includegraphics[width=7.6cm,clip]{9453f12.eps} \end{figure}$

Figure 6: This figure shows the evolution of the system for models in which $m_{\rm p}=1~M_{\rm J}$ . Here, the left and right panels correspond to $t_{\rm acc}=t_{\rm dyn}$ and $t_{\rm acc}=10~t_{\rm dyn}$ respectively. Top: evolution of the semi-major axis of the planet. Middle: evolution of the eccentricities for the planet (black), binary (red) and disc (green). Bottom: evolution of the longitudes of pericentre for the planet (black), binary (red) and disc (green).

Open with DEXTER

For the two accretion time scales we consider, the orbital evolution of a planet with final mass $m_{\rm p}=1~M_{\rm J}$ is illustrated in the upper panel of Fig. 6. In the model with $t_{\rm acc}=t_{\rm dyn}$ the planet reaches one Jupiter mass in $\sim$ 5.5 $\times$ 10³ binary orbits while it reaches this mass in $\sim$ 1.4 $\times$ 10⁴ orbits in the model with $t_{\rm acc}=10~t_{\rm dyn}$ (see Fig. 2). The early evolution of the planet is similar to that described for Saturn-mass planets, involving inward migration of the planet and continued growth of its eccentricity. However, contrary to models in which $m_{\rm p}=1~M_{\rm S}$ , simulations with $m_{\rm p}=1~M_{\rm J}$ resulted in the temporary formation of the 4:1 resonance between the planet and binary (we note that migration reversal occurs for the planets with $m_{\rm p}=1~M_{\rm S}$ before they reach the 4:1 resonance).

For the calculation with $t_{\rm acc}=t_{\rm dyn}$ , the time evolution of the resonant angles $\psi _1$ , $\psi _2$ , $\psi _3$ and $\psi _4$ associated with the 4:1 resonance is displayed in Fig. 7. These are given by:

$\begin{displaymath}% \begin{array}{l} \psi_1=4\lambda_{\rm b}-\lambda_{\rm p}-3\... ...-\lambda_{\rm p}-2\omega_{\rm p}-\omega_{\rm b},\\ \end{array}\end{displaymath}$

(3)

$\begin{figure} \par\includegraphics[width=7.1cm,clip]{9453f13.eps} \end{figure}$	Figure 7: This figure shows, for the simulation with $m_{\rm p}=1~M_{\rm J}$ and $t_{\rm acc}=t_{\rm dyn}$ , the evolution of the resonant angles $\psi _1$ , $\psi _2$ , $\psi _3$ and $\psi _4$ associated with the 4:1 resonance.
Open with DEXTER

A similar outcome is observed in the simulation with $t_{\rm acc}=10~t_{\rm dyn}$ . Here the planet is scattered at $t\sim 4$ $\times$ 10⁴ P and evolves subsequently on an orbit with significantly larger semimajor axis ( $a_{\rm p}\sim 1.6$ ) and eccentricity ( $e_{\rm p}\sim 0.3)$ . Since a close encounter is generally a chaotic dynamical process, these values are quite different from the ones observed in the simulation with $t_{\rm acc}=t_{\rm dyn}$ . In fact, these strongly depend on the encounter geometry and therefore may significantly differ from one simulation to the other.

$\begin{figure} \par\includegraphics[width=5.2cm,clip]{9453f14.eps}\hspace{3mm} ... ....eps}\hspace{3mm} \includegraphics[width=5.2cm,clip]{9453f17.eps} \end{figure}$	Figure 8: This figure shows, for the model in which $m_{\rm p}=1~M_{\rm J}$ and $t_{\rm acc}=t_{\rm dyn}$ , snapshots of the disc surface density at times shown above the plots. Here, the planet is represented by a white circle.
Open with DEXTER

Prior to this close encounter however, the evolution of the planet differed noticeably from that corresponding to the model with $t_{\rm acc}=t_{\rm dyn}$ . It appears that in the calculation with $t_{\rm acc}$ = $10~t_{\rm dyn}$ , the 4:1 resonance becomes rapidly undefined after its formation, which arises at $t\sim 4$ $\times$ 10⁴. This is because, compared with the run in which $t_{\rm acc}=t_{\rm dyn}$ , the planet mass is significantly lower when the 4:1 resonance is established, leading to a weaker resonant locking. Once the resonance is broken, the planet is located just outside of the 4:1 resonance and migrates slightly outward until $t \sim 1.7$ $\times$ 10⁴ P. Outward migration is induced by positive disc torques due to the planet eccentricity having reached $e_{\rm p}\sim 0.15$ after resonance breaking. Then, the planet migrates inward and approaches the 4:1 resonance again, but its eccentricity remains sufficiently high for there to be a close encounter between the planet and central binary system leading to the planet being scattered by the binary. The planet and binary undergo multiple encounters during this phase, with significant changes in their orbital elements occuring (see middle-right panel of Fig. 6). The planet ends up orbiting interior to the 4:1 resonance where it migrates inward toward the 3:1 resonance.

In both simulations, the evolution of the planet subsequent to this initial scattering is quite similar. In the following, we use the results of the simulation with $t_{\rm acc}=t_{\rm dyn}$ to discuss in more details how the evolution of the system proceeds after this event. At $t\sim 10^4~P$ , interaction with the disc and central binary finally results in the planet settling into an orbit further out in the disc with eccentricity $e_{\rm p} \simeq$ 0.1-0.15. The planet now forms a gap in the disc and begins to migrate inward slowly under type II migration. The evolution of the disc and planet plus binary system for this run is presented in Fig. 8, which shows snapshots of the disc surface density at different times. The first panel corresponds to a time shortly before the formation of the 4:1 resonance while the three other ones display the state of the system just after the initial scattering. These show clearly that the inner disc is progressively lost through the inner boundary as a result of viscous evolution, thereby leading to a reduction of the (positive) inner disc torques and consequently to the inward migration of the planet.

$\begin{figure} \par\includegraphics[width=7.6cm,clip]{9453f18.eps}\hspace*{4mm} \includegraphics[width=7.6cm,clip]{9453f19.eps} \end{figure}$

Figure 9: This figure shows the evolution of the system for models in which $m_{\rm p} > 1~M_{\rm J}$ . Here, the left and right panels correspond to $t_{\rm acc}=t_{\rm dyn}$ and $t_{\rm acc}=10~t_{\rm dyn}$ respectively. Top: evolution of the semi-major axis of the planet. Middle: evolution of the eccentricities for the planet (black), binary (red) and disc (green). Bottom: evolution of the longitudes of pericentre for the planet (black), binary (red) and disc (green).

Open with DEXTER

At the end of the run with $t_{\rm acc}=t_{\rm dyn}$ , the final fate of the planet is still uncertain despite the very long time scale covered by the simulation. However, it is likely that the continued inward migration of the planet, combined with the growth of the eccentricities of both the binary and planet (see middle panel of Fig. 6), will result in further close encounters and eventually in the planet being completely ejected from the system. Such an outcome is observed in the calculation with $t_{\rm acc}=10~t_{\rm dyn}$ . In that case, this arises because inward migration of the planet causes the temporary formation of the 3:1 resonance with the binary at $t\sim 1.35$ $\times$ 10⁵ P, which makes the planet eccentricity increase up to $e_{\rm p}\sim 0.15$ . This leads to close encounters between the planet and the secondary star, resulting in the scattering and ejection of the planet from the system.

The lower panel of Fig. 6 shows the evolution of the disc, binary and planet longitudes of pericentre for both models. Interestingly, we see that here the planet is aligned with neither the disc nor the binary. Therefore, relative to simulations with $m_{\rm p}=1~M_{\rm S}$ , the probability of close encounters between the planet and binary is increased, thereby leading to a system which is more likely to become unstable.

Lastly, we notice that the results of these runs are broadly consistent with previous hydrodynamical simulations of jupiter-mass planets embedded in circumbinary discs (Nelson 2003) which indicated that indeed, trapping into 4:1 resonance followed by a scattering through a close encounter with the binary is a possible outcome of such systems.

3.3 Evolution of giants with m_p > 1 M_J

The evolution of accreting bodies with final masses $m_{\rm p} > 1~M_{\rm J}$ is depicted, for both values of the accretion time scale, in the upper panel of Fig. 9. With respect to the simulations presented in Sect. 3.2, here the calculations differ in that the planet can continue to accrete gas once its mass has attained $m_{\rm p}=1~M_{\rm J}$ . Although subsequent evolution may differ, this implies that while $m_{\rm p} < 1~M_{\rm J}$ the evolution of the planet is the same as the planets considered in Sect. 3.2. Therefore, when discussing the results below, we consider only the evolution of the system from the time when the planet mass has reached $m_{\rm p}=1~M_{\rm J}$ .

For the model with $t_{\rm acc}=t_{\rm dyn}$ , the planet is in 4:1 resonance with the binary when its mass reaches and exceeds that of Jupiter, which arises at $t\sim 5.5$ $\times$ 10³ P (see Fig. 2). Figure 10 displays the evolution of the resonant angles associated with the 4:1 resonance. Comparing this figure with Fig. 7, we see that the 4:1 resonance breaks at $t\sim 8.5$ $\times$ 10³ P with $m_{\rm p}=1~M_{\rm J}$ whereas it breaks at $t\sim 1.15$ $\times$ 10⁴ P in the one with $m_{\rm p} > 1~M_{\rm J}$ . This supports the idea that the resonant interaction is stronger for higher planet masses. Once again, the resonance drives the planet eccentricity up to $e_{\rm p}\sim 0.2$ , until the planet undergoes a close encounter with the binary and is subsequently ejected from the system at $t\sim 1.15$ $\times$ 10⁴ P. Here, it is worth noting that the planet mass is $m_{\rm p}\sim 1.7~M_{\rm J}$ when this occurs (see Fig. 2).

For the model with $t_{\rm acc}=10~t_{\rm dyn}$ , the planet orbits just beyond the location of the 4:1 resonance with the binary when its mass attains $m_{\rm p}=1~M_{\rm J}$ , which corresponds to $t\sim 1.4$ $\times$ 10⁴ P. From this point in time until $t\sim 2.1$ $\times$ 10⁴ P, the orbital evolution of the planet is similar to that found in the model with $m_{\rm p}=1~M_{\rm J}$ and $t_{\rm acc}=10~t_{\rm dyn}$ (see Sect. 3.2). Then, the large value of its eccentricity ( $e_{\rm p}\sim 0.13$ ) causes the planet to undergo a close encounter with the binary and to be scattered out. Interestingly, this scattering leads to the temporary formation of the 6:1 resonance between the planet and the binary at $t\sim 2.15$ $\times$ 10⁴ P. As can be seen in the middle panel of Fig. 9, which shows the eccentricities of the planet, disc and binary, trapping into 6:1 resonance makes the planet eccentricity increase to $e_{\rm p} \ge 0.2$ . At $t\sim 2.3$ $\times$ 10⁴ P, a close encounter between the binary and planet occurs, and the planet is consequently scattered out on a high eccentricity orbit with $a_{\rm p}\sim 1.8$ .

Just after this scattering, the planet mass has reached $m_{\rm p}\sim 3~M_{\rm J}$ and the eccentricities of the disc and planet have attained $e_{\rm d}\sim 0.08$ and $e_{\rm p}\sim 0.1$ respectively. In agreement with Nelson (2003), we find that the interaction between the eccentric disc and the eccentric planet induces outward migration of the latter. The total torque exerted by the disc on the planet, as well as the torques due to the disc interior and exterior to $a_{\rm p}$ , are presented in Fig. 12. We see that the time-averaged total disc torque is clearly positive from $t\sim 2.3$ $\times$ 10⁴ P onward, and that the torques oscillate with a period of $\sim$ 3 $\times$ 10³ binary orbits. Interestingly, these time variations correspond closely to the precession of the planet relative to that of the disc. The reason for this is simply that the eccentric disc exerts an orbit-averaged positive torque on the eccentric planet when the orbits are misaligned, and exerts a negative torque when the orbits are aligned. An antialigned configuration causes the planet to interact with matter whose angular velocity is greater than that of the planet at apocentre, and this leads to the planet experiencing a strong positive torque.

At later times, the simulation indicates that the outward migration of the planet can be sustained over long time scales. There are a number of reasons for this. First, the planet maintains an eccentric orbit and experiences a strong positive torque at apocentre when the disc and planet orbits are antialigned; this positive torque experienced by the planet implies a negative torque experienced by the outer disc material, causing gas to flow through the planet orbit to form an inner disc. Figure 11 shows the state of the disc, planet and binary at different times. The first panel shows the disc just before the planet enters the 4:1 resonance, and the second panel shows the system just after the planet has been scattered outward. The third and fourth panels show the increase in inner disc size and the outward migration of the planet. During its transition from outer to inner disc, the disc material continues to exert a positive torque on the planet, manifested through a corotation torque. Moreover, the existence of an inner disc provides another source of positive Lindblad torques and assists the planet in maintaining outward migration.

Examination of the torques exerted by the disc on the planet (see Fig. 12) shows that they weaken from $t\sim 4$ $\times$ 10⁴ P. This is due mainly to the increase in gap size generated by the planet as it grows in mass, combined with the fact that the disc contains a finite reservoir of gas in our model. The increase in gap size clearly affects the torques due to both outer and inner discs.

At the end of the run, the planet has migrated to the outer edge of the disc and its mass has attained $m_{\rm p}\sim 5~M_{\rm J}$ . This suggests that another avenue for evolution of giant planets whose mass exceeds 1 $M_{\rm J}$ is long-term outward migration, in addition to the possibility of being scattered out of the system as displayed by the runs described previously. Establishing the ratio of these two outcomes will require a large suite of simulations with slightly varying initial conditions, task which is beyond the scope of the paper. It appears from our results, however, that the probability of scattering and ejection is greater than that of prolonged outward migration.

$\begin{figure} \par\includegraphics[width=7.2cm,clip]{9453f20.eps} \end{figure}$	Figure 10: This figure shows, for the simulation with $m_{\rm p} > 1~M_{\rm J}$ and $t_{\rm acc}=t_{\rm dyn}$ , the evolution of the resonant angles $\psi _1$ , $\psi _2$ , $\psi _3$ and $\psi _4$ associated with the 4:1 resonance.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=5.2cm,clip]{9453f21.eps}\hspace{3mm} ... ....eps}\hspace{3mm} \includegraphics[width=5.2cm,clip]{9453f24.eps} \end{figure}$	Figure 11: This figure shows, for the model in which $m_{\rm p} > 1~M_{\rm J}$ and $t_{\rm acc}=10~t_{\rm dyn}$ , snapshots of the disc surface density at times shown above the plots. Here, the planet is represented by a white circle.
Open with DEXTER

4 Summary and conclusion

In this paper we have presented the results of hydrodynamic simulations aimed at studying the formation and evolution of giant planets embedded in circumbinary discs.

We focused on a model in which a $20~M_\oplus$ core initially trapped at the edge of the cavity formed by the binary can slowly accrete material from the disc. We examined how the final outcome of the system depends on the accretion rate onto the planet and on the final planet mass attained. For each value of the accretion rate considered, we performed three calculations. In two of the three simulations, we assumed that accretion stops when the mass of the planet has reached either $m_{\rm p}=1~M_{\rm S}$ or $m_{\rm p}=1~M_{\rm J}$ . In the remaining case, we allowed the planet to accrete gas freely from the disc in such a way that its final mass was $m_{\rm p} > 1~M_{\rm J}$ . The simulations show different outcomes, depending on the final mass of the planet:

$\begin{figure} \par\includegraphics[width=7.5cm,clip]{9453f25.eps} \end{figure}$	Figure 12: This figure shows, for the model with $m_{\rm p} > 1~M_{\rm J}$ and $t_{\rm acc}=10~t_{\rm dyn}$ , the evolution of the torques exerted by the disc on the planet. The blue (red) line corresponds to the inner (outer) disc torques and the black line corresponds to the net torques exerted on the planet.
Open with DEXTER

A number of other issues remain to be resolved when considering the early stages of planet formation in circumbinary discs. For example, the question of whether or not planetary cores can grow due to planetesimal accretion needs to be examined in more detail. An eccentric binary can lead to the formation of an eccentric disc, in such a way that planetesimal accretion is prevented except in the outer regions of the disc. The ability of planetesimals to accrete, with care being taken to simulate the structure of the circumbinary disc, will be the subject of a future paper.

Another issue relates to the fact that we have only simulated a two dimensional system in this paper. The close encounters experienced by planets with the central binary are likely to produce significantly inclined orbits, and the resulting disc-planet interaction is likely to be modified (weakened) by this, leading to potentially different outcomes to those observed in our 2D runs. Given the very long evolution times involved, however, performing full 3D simulations is beyond current computational capability. There remains scope, however, for a more approximate 3D treatment of this problem.

On the formation and migration of giant planets in circumbinary discs

1 Introduction

2 Hydrodynamical model

2.1 Numerical method

2.2 Initial conditions

2.3 Boundary conditions

3 Results

3.1 Evolution of Saturn-mass planets

3.2 Evolution of Jupiter-mass planets

3.3 Evolution of giants with m_p > 1 M_J

4 Summary and conclusion

References