A&A 374, 1092-1099 (2001)
DOI: 10.1051/0004-6361:20010779

On disc driven inward migration of resonantly coupled planets with application to the system around GJ876

M. D. Snellgrove - J. C. B. Papaloizou - R. P. Nelson

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

Received 3 April 2001 / Accepted 28 May 2001

Abstract
We consider two protoplanets gravitationally interacting with each other and a protoplanetary disc. The two planets orbit interior to a tidally maintained disc cavity while the disc interaction induces inward migration. When the migration is slow enough, the more rapidly migrating outer protoplanet approaches and becomes locked in a 2:1 commensurability with the inner one. This is maintained in subsequent evolution. We study this evolution using a simple analytic model, full hydrodynamic 2D simulations of the disc planet system and longer time N-body integrations incorporating simple prescriptions for the effects of the disc on the planet orbits. The eccentricities of the protoplanets are found to be determined by the migration rate and circularization rate induced in the outer planet orbit by the external disc. We apply our results to the recently discovered resonant planets around GJ876. Simulation shows that a disc with parameters expected for protoplanetary discs causes trapping in the 2:1 commensurability when the planets orbit in an inner cavity and that eccentricities in the observed range may be obtained.

Key words: stars: planetary systems: formation - stars: protoplanetary discs - celestial mechanics -
stars: individual: GJ876


   
1 Introduction

The recent discovery of extrasolar giant planets orbiting around nearby solar-type stars (Marcy & Butler 1998, 2000) has stimulated renewed interest in the theory of planet formation. The objects observed so far have masses, $M_{\rm p}$, that are characteristic of giant planets ( $0.4 \; {M}_{\rm J} \; \raisebox{-.8ex}{$\buildrel{\textstyle<}\over\sim$ }\; M_{\rm p} \ \raisebox{-.8ex}{$\buildrel{\textstyle<}\over\sim$ }\; 11 \; {M}_{\rm J}$), $M_{\rm J}$ denoting a Jupiter mass. The orbital semi-major axes are in the range $0.04 \; {\rm AU} \; \raisebox{-.8ex}{$\buildrel{\textstyle<}\over\sim$ }\; a \; \raisebox{-.8ex}{$\buildrel{\textstyle<}\over\sim$ }\; 2.5 \; {\rm AU}$, and orbital eccentricities in the range $ 0.0 \; \raisebox{-.8ex}{$\buildrel{\textstyle<}\over\sim$ }\; e \; \raisebox{-.8ex}{$\buildrel{\textstyle<}\over\sim$ }\; 0.67$ (Marcy & Butler 2000).

Disc-protoplanet interactions have been invoked to explain the presence of giant planets orbiting close to their host stars through inward orbital migration induced through disc-protoplanet tidal interaction (e.g. Papaloizou et al. 1999; Lin et al. 2000). Up to now, for the most part extrasolar planets appear to be isolated. However, a few multiple systems are known. The configuration of these may contain important information about their origin and possible migration history. Of special interest is the recently discovered system around GJ876. This is found to be close to a 2:1 commensurability. Such a configuration is indeed suggestive of orbital migration. Commensurable satellite systems such as the Galilean satellites are thought to owe their origin to migration induced by tidal interaction with the central planet (e.g. Goldreich 1965).

Recent simulations of single protoplanets in the observed mass range (Kley 1999; Bryden et al. 1999; Lubow et al. 1999) interacting with a disc with parameters thought to be typical of protoplanetary discs, but constrained to be in circular orbit, indicate gap formation and upper mass limit consistent with the observations. Nelson et al. (2000) relaxed the assumption of fixed circular orbits, found inward migration and that the disc-protoplanet interaction leads to strong eccentricity damping. Due to accretion onto the central star, an inner cavity was formed in the disc interior to which the protoplanet orbited.

Simulations of two planets interacting with a disc have been performed by Kley (2000), Bryden et al. (2000), and Masset & Snellgrove (2001). So far inward migration of two planets locked into a 2:1 commensurability has not been simulated. However, Bryden et al. (2000) found a tendency for gap material between the two planets in fixed circular orbits to be cleared, ending up interior to the inner planet orbit or exterior to the outer planet orbit.

Taken together the above results suggest a natural outcome of two protoplanets interacting with a disc is that they orbit interior to an inner disc cavity while the external disc causes inward migration of the outer orbit. This catches up the inner orbit leading to the possibility of resonant interaction.

It is the purpose of this paper to investigate such resonant interaction and whether, for reasonable protoplanetary disc models, it leads to a locking such that the planets subsequently migrate maintaining the commensurability. Similar behavior occurs as a result of the tidally induced migration of the Gallilean satellites (e.g. Goldreich 1965; Lin & Papaloizou 1979).

Here we shall assume the prior evolution of the system leads to the orbital separation of the planets being slightly larger than that required for a strict 2:1 commensurability without considering the history in detail as it is beyond the scope of this paper. However, we comment that this might have been complicated with the planet masses varying with time through mass accretion from the disc.

Analytic methods, N-body integrations and direct two dimensional numerical simulations are used to investigate the evolution and found to give consistent results.

We consider two protoplanets gravitationally interacting with each other and a protoplanetary disc. The two planets orbit interior to a tidally maintained disc cavity while the disc interaction induces inward migration. This migration is of type II and so is regulated by the magnitude of the disc viscosity. When the migration is slow enough, the more rapidly migrating outer protoplanet approaches and becomes locked in a 2:1 mean motion commensurability with the inner one. The commensurability persists in subsequent evolution. The eccentricities of the protoplanets are increased by the resonant perturbations. In principle, tidal interaction with the disc may cause growth or decay of orbital eccentricity according to the protoplanet mass and physical conditions in the disc (Lin & Papaloizou 1993). Wide gaps, small disc viscosity and high masses favour eccentricity growth while narrow gaps, larger viscosity and lower masses lead to decay of eccentricity. For disc viscosities normally assumed in protoplanetary disc models eccentricity growth occurs only for protoplanet central star mass ratios greater than ${\sim} 0.02$ (Papaloizou et al. 2001). Thus for the conditions considered here, orbits are circularized through interaction with the disc (e.g. Goldreich & Tremaine 1981). A balance is achieved in which the eccentricities of the resonantly coupled protoplanets are determined by their migration rate and the circularization rate induced in the outer planet orbit by the external disc. We apply our results to the recently discovered resonant planets around GJ876. Simulation shows that migration induced by a disc with parameters expected for protoplanetary discs results in trapping in the 2:1 commensurability when the planets orbit in an inner cavity. Eccentricities in the observed range may be obtained. Further studies using N-body integrations indicate that the planetary system will remain stable for at least $2 \times 10^{7}$ orbits when the external disc is removed.

We note that we are restricting this analysis to mass ratios such that both protoplanets will open a gap, and maintain the cavity, corresponding to protoplanet masses above around one Jupiter mass (assuming a solar mass primary). This is in contrast to previous work by Masset & Snellgrove (2001), where the outer protoplanet had insufficient mass to fully open a gap, which resulted in a migration reversal (outward migration) when in resonance with the more massive inner protoplanet. In the case presented here, as both protoplanets are in a cavity, the migration is driven by the outer protoplanet interacting with the outer disc only, which results in inward migration.

In Sect. 2 we describe a simple analytic model of two migrating protoplanets in a 2:1 commensurability. The eccentricities of the protoplanet orbits are related to the migration rate, and circularization rate induced by the disc. In Sect. 3 we describe a simulation of two planets orbiting in an inner disc cavity. Parameters appropriate to GJ867 are adopted. This demonstrates resonant trapping and that eccentricities of the observed magnitude may be produced.

In Sect. 4 we describe N-body calculations confirming the above conclusions and indicating the long term stability of the system. Finally in Sect. 5 we discuss our results.

   
2 A simple model

We consider a system consisting of 2 planets and a primary star moving under their gravitational attraction. When there are no disc interactions and the motion is conservative the system is conveniently expressed in Hamiltonian form using Jacobi coordinates (e.g. Sinclair 1975). The coordinates, ${\vec r}_2,$ of inner planet of reduced mass m2are referred to the central star of mass M*and the coordinates of the outer planet, ${\vec r }_1,$of reduced mass m1 are referred to the center of mass of the central star and inner planet. The Hamiltonian can be written correct to second order in the planetary masses as

$\displaystyle H = {1\over 2} ( m_1 \vert \dot {\vec r}_1\vert^2 +m_2\vert \dot ... ...} + {Gm_{1}m_2 {\vec r}_1\cdot {\vec r}_2 \over \vert {\vec r}_{1}\vert^3}\cdot$     (1)

Here M*1=M*+m1, M*2= M* + m2 and $ {\vec r}_{12}= {\vec r}_{2}- {\vec r}_{1}.$The Hamiltonian can be expressed in terms of the osculating semi-major axes, eccentricities and longitudes of periastron $a_i,e_i,\varpi_i, i=1,2$ respectively as well as the longitudes $\lambda_i,$ and the time t. We recall that $\lambda_i = n_i(t-t_{0i}) + \varpi_i,$ with ni being the mean motion and t0i giving the time of periastron passage. The energy is given by Ei = -GmiM*i/(2ai),and the angular momentum $h_i = m_i\sqrt{GM_{*i}a_i(1-e_i^2)}$ which may be used to describe the motion instead of ai and ei. Only terms first order in the eccentricities and involving the resonant angles $\phi = 2\lambda_1-\lambda_2-\varpi_1, $ and $\psi = 2\lambda_1-\lambda_2-\varpi_2,$ are retained in the expansion (see Sinclair 1975). Then the perturbing part of the Hamiltonian $(\propto\! m_1 m_2)$ can be written

\begin{displaymath}H _{12}= -{Gm_1m_2\over a_1}\left({Be_1\cos\phi} -{Ce_2\cos\psi}\right), \end{displaymath} (2)

where $C=2b_{1/2}^{(2)}({\alpha}) +(1/2)\alpha {\rm d} b_{1/2}^{(2)}({\alpha})/{\rm d}\alpha$ and $B=1.5b_{1/2}^{(1)}({\alpha}) +(1/2)\alpha {\rm d} b_{1/2}^{(1)}({\alpha})/{\rm d}\alpha -2\alpha.$Here b1/2(n) denotes the usual Laplace coefficient and $\alpha= a_2/a_1.$ From now on we replace M*i by M*.

2.1 Orbital precession

We may also take into account additional gravitational forces that may produce precession of the planetary orbits. These could result from the mass distribution provided by the disc exterior to the outer planet or the time averaged mass distribution of the planets themselves. These effects are not included in the resonant Hamiltonian given above. Accordingly we add to it another Hamiltonian

 \begin{displaymath}% H_{\rm int} = - {1\over 2}m_1 n_1 a_1^2 \omega_{\rm pr1}\,e_1^2 - {1\over 2}m_2 n_2 a_2^2 \omega_{\rm pr2}\,e_2^2 . \end{displaymath} (3)

Here, as can be verified from the equations of motion, we have adopted a parameterization such that the precession frequency induced in the orbit of mi is $\omega_{\rm pri}(a_i).$ This prescription does not allow the precession frequency to depend on other quantities. However, appropriate matching can be carried out for a particular case under consideration.

2.2 Basic equations

The equations of motion are derived from:

${{\rm d}E_i}/{{\rm d}t} =-n_i\partial H' / \partial \lambda_i -(n_1 T +D)\delta_{i1}, $
${{\rm d}h_i}/{{\rm d}t} =-{\partial H' / \partial \lambda_i} -{\partial H' / \partial \varpi_i} -T\delta_{i1},$
$ {{\rm d}\lambda_i / {\rm d}t} = n_i +n_i {\partial H' / \partial E_i} + {\partial H' / \partial h_i}, $
${{\rm d}\varpi_i / {\rm d}t} ={\partial H' / \partial h_i}, $

with $H' = H_{12} + H_{\rm int}.$ These can be obtained from Hamilton's equations (e.g. Brouwer & Clemence 1961) to which we have added, for the outer planet m1,an additional external torque - Twith an associated orbital energy loss rate n1 T together with additional orbital energy dissipation rate D. The torque and dissipation rate could be produced by tidal interaction with the disc leading to inward migration and orbital circularization. We thus obtain to lowest order in the planetary eccentricities and perturbing masses.
 
$\displaystyle % {{\rm d}n_1\over {\rm d}t} = {6n_1^2 m_2 \over M_* }\left(Be_1\sin\phi -Ce_2\sin\psi\right) + {3 n_1 a_1\over G M_*m_1}\left(n_1 T +D\right)$     (4)


\begin{displaymath}{{\rm d}n_2\over {\rm d}t} =-{3n_2^2 m_1 a_2\over M_* a_1 }\left(Be_1\sin\phi -Ce_2\sin\psi\right)\end{displaymath} (5)


\begin{displaymath}{{\rm d}e_1\over {\rm d}t} =-{ m_2 n_1\over M_* }B\sin\phi -{D a_1\over GM_* m_1 e_1}\end{displaymath} (6)


\begin{displaymath}{{\rm d}e_2\over {\rm d}t} ={ m_1 n_2 a_2\over M_* a_1 }C\sin\psi \end{displaymath} (7)


 \begin{displaymath}{{\rm d}\phi \over {\rm d}t} = 2n_1- n_2 - {m_2\over e_1 M_* }n_1 B\cos\phi -\omega_{\rm pr1} \end{displaymath} (8)


 \begin{displaymath}{{\rm d}\psi \over {\rm d}t} = 2n_1- n_2 + { m_1 a_2\over e_2 a_1 M_* }n_2 C\cos\psi -\omega_{\rm pr2}. \end{displaymath} (9)

2.3 Stationary solutions

When no migration or circularization occurs (T=D=0) equilibrium solutions may exist such that $\psi $ and $\phi $ are either zero or $\pi.$ Each of n1,n2,e1,e2 are then constant. A relation between the eccentricities then follows from Eqs. (8) and (9) in the form

 
$\displaystyle e_2a_1( e_1 M_* \omega_{pr1}+ m_2 n_1 B \cos \phi) = e_1 (-m_1 a_2 n_2 C\cos\psi+ e_2 a_1 M_* \omega_{pr2}).$     (10)

This condition matches the precession rates of the orbits of the two planets. Also 2n1=n2.Noting that the eccentricities are positive, when they are of very small magnitude, the precessional terms become negligible and there is a solution with $\psi=0, \phi =\pi$or $\psi=0, \phi =\pi.$ In either case we have

 
m2 a1 e2 n1 B = m1 a2 e1n2 C. (11)

For larger eccentricities the precessional terms may become important in Eq. (10) and then solutions with $\psi=0, \phi =0,$ may occur. Then Eq. (10) gives

 \begin{displaymath}m_2 a_1 e_2 n_1 B + m_1 a_2 e_1n_2 C= e_1 e_2 a_1 M_* ( \omega_{\rm pr2} -\omega_{\rm pr1}). \end{displaymath} (12)

For stable solutions, when perturbed, the angles may undergo librations about their equilibrium points (e.g. Sinclair 1975). There are two frequencies of oscillation $\nu_1,\nu_2$ being given for any $\psi, \phi $ in the limit of small eccentricities by $\nu_1^2 = (m_2 n_1 B)^2/(M_*e_1)^2,$and $\nu_2^2 =(m_1n_2a_2C)^2/(a_1M_* e_2)^2.$We look for a solutions with migration which are close to stable solutions of this type.

2.4 Resonant migration

We look for solutions of Eqs. (4-9) corresponding to the situation where the two planets migrate inwards locked in resonance with n1/n2 maintained nearly equal to 1/2 while the eccentricities remain nearly constant. The tendency for the resonant coupling to excite the eccentricities is counterbalanced by circularization through the action of $D \equiv (GM_* m_1 e_1^2) /(a_1 t_c)$ which defines the circularization time for e1. Similarly Eq. (4) defines an inward migration timescale $t_{\rm mig}= GM_* m_1/(3 T a_1n_1).$ We begin by supposing that the angle $\psi $ executes a libration about zero such that the mean rate of change of e2 is zero. Similarly the mean rates of change of n1 and n2 induced by $\psi $ are zero. Such a libration is seen in simulations. We also suppose the angle $\phi $ either librates or circulates but in such a way that the correspondingly induced mean rates of change are not zero. The simplest example is when the angle executes a very small or even zero amplitude libration about a value slightly offset from zero or $\pi$ (e.g. Lin & Papaloizou 1979). We suppose the circulation/libration periods to be short compared to the timescale of migration so that averaging is possible. We denote the average of $n_1 e_1\sin\phi $ by $\delta.$Then averaging Eqs. (4-9) gives

 \begin{displaymath}% {\overline{{1\over n_1}{{\rm d}n_1\over {\rm d}t}}}={6 m_2 B\delta\over M_* } +{3 a_1 \over G M_*m_1}\left(n_1 T +D\right) \end{displaymath} (13)


\begin{displaymath}% {\overline{{1\over n_2}{{\rm d}n_2\over {\rm d}t}}}=-{6 m_1 B a_2 \delta\over M_* a_1}\end{displaymath} (14)


 \begin{displaymath}% {\overline {{{\rm d}e_1^2\over {\rm d}t}}}=-{2 m_2 B \delta\over M_* } -{2D a_1 \over GM_* m_1 } \equiv 0 \end{displaymath} (15)


 \begin{displaymath}% {\overline{{{\rm d}\phi \over {\rm d}t}}} = -{m_2\over e_1 ... ... \over e_2 M_* a_1 }n_2 C -\omega_{\rm pr1} +\omega_{\rm pr2}. \end{displaymath} (16)

Here on the left hand sides the overline denotes the time average and for simplicity of notation this has been dropped from the right hand sides. The resonance condition also implies ${\overline{{1\over n_2}{{\rm d}n_2\over {\rm d}t} }}= {\overline{{1\over n_1}{{\rm d}n_1\over {\rm d}t} }}.$Using this and eliminating $\delta$ from the above, we obtain for the rate of increase of n1 through migration

\begin{displaymath}% {\overline{{1\over n_1}{{\rm d}n_1\over {\rm d}t} }} ={3 a_... ...\left(n_1 T +D\right) \left({m_1a_2\over m_1a_2+m_2a_1}\right).\end{displaymath} (17)

Also (15) gives for the eccentricity balance

 \begin{displaymath}% e_1^2 = {t_{\rm c} m_2a_1\over 3 t_{\rm mig}(2 m_1 a_2 + m_2a_1)}\cdot \end{displaymath} (18)

The above determines the eccentricity of the outer planet e1 as a function of $t_{\rm c}$ and $t_{\rm mig}.$For a system with m1/m2 =3, we get $e_1 \sim \sqrt{0.07t_{\rm c}/t_{\rm mig}}.$

For e1 in the 0.01 range we need $t_{\rm c} \sim10 $ orbits if $t_{\rm mig} \sim 10^4$ orbits. The eccentricity of the inner planet is determined by Eq. (16). For small amplitude librations this is given by Eq. (10). That would still apply when $\phi $is circulating provided the cosines are time averaged and the mean circulation rate is small.

  
3 A simulation of two migrating resonantly coupled planets

The protoplanetary disc is numerically simulated using an Eulerian 2D hydrodynamic code. The code used is a modified version of NIRVANA, which has been described, tested and used successfully elsewhere on a similar problem involving interacting planets (Masset & Snellgrove 2001). Incorporated with the hydrodynamic code is a 4th-order Runge-Kutta integrator which is used to evolve the orbits of the two planets. The gravitational forces calculated from the disc model are used in the equations of motion of the planets, and the disc itself responds to the planetary potential. Hence the system evolves in a self-consistent fashion. In order to obtain the long integration times needed for simulations of this type the FARGO algorithm (Masset 2000) is applied. However, tests have shown that the results are not affected by this.

We use a 2D cylindrical $(r,\varphi)$ grid with 200 radial zones distributed uniformly between r=0.4 and r=3.47 in dimensionless units and 300 azimuthal zones. We apply outflow conditions at the inner boundary to simulate the accretion of disc material onto the central star.

3.1 Physical setup

We attempt to simulate the resonance locking of the system GJ876 via tidally induced migration of the planets, using plausible values of the disc parameters. The disc is assumed to be thin and isothermal, with constant aspect ratio h/r = 0.07, and a constant Shakura & Sunyaev (1973) $\alpha$-viscosity prescription with ${\alpha} = 2 \times 10^{-3}$ is adopted. The two planets are initially in circular orbits coplanar with the disc, at radial locations r1 =1.0 and r2 =0.6. Hence the outer planet is located outside the exact 2:1 commensurability ( r2 =0.63). The planet masses are chosen to correspond to the minimum mass ratios obtained from observations (Marcy et al. 2001). With masses normalised so that stellar mass is M* = 1, this corresponds to $m_{1} = 6 \times 10^{-3}$ and $m_{2} = 1.8 \times 10^{-3}$. The planet masses are fixed as the planets are assumed to be no longer accreting material from the disc.

The disc is prescribed an initial surface density $\Sigma_{0}$ corresponding to what would give a disc mass of $2 \times 10^{-3}$ within the orbit of the outer planet. However we assume that both planets are located inside a tidally truncated cavity located at r < 1.3, with low surface density $\Sigma_{\mathrm{cavity}} = 0.01\Sigma_{0}$. This cavity is supposed to have already been cleared by the tidal action of the two planets. Between 1.3 < r < 1.5 the surface density is prescribed such that $\ln \Sigma$ linearly joins to $\Sigma_{0}$.

3.2 Results


  \begin{figure} \par\includegraphics[width=8.8cm,height=8.8cm,clip]{MS1339f1.ps}\par\end{figure} Figure 1: The lower plots show the semi major axes a (lower left) and eccentricities e (lower right) for both planets. The upper plots show the resonance angles $\phi $ and $\psi $. The time unit is measured in orbit periods of the disc material at r=1.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7.9cm,height=8.7cm,clip]{MS1339f2.ps}\par\end{figure} Figure 2: Surface density plot for inner regions of the disc at a time corresponding to 300 orbits. Darker areas on the plot correspond to regions of lower surface density. The disc cavity is clearly visible, as are the density waves outside the cavity. Inside the cavity are the wakes of the two planets.
Open with DEXTER

The tidal interaction of the planets with the disc material causes the planets to migrate inwards (see Fig. 1). The inner planet is deep within the cavity and only interacts with low surface density material and thus migrates slowly. The outer planet has its outer 2:1 Lindblad resonance located outside the cavity and in the body of the disc. Hence there is more material exerting a negative torque on the planet, and therefore it migrates faster, despite its larger mass. The ratio of semi-major axes a1 / a2 of the planets decreases until the planets "lock'' into a 2:1 commensurability with $n_{2} \approx 2n_{1}$ at a time $t \approx 400$ orbits. Both planets then subsequently migrate inwards a further 10% maintaining this ratio, showing the resonance to be robust.

Figure 1 also shows the calculated values of the resonance angles $\phi $ and $\psi $. Once the commensurability lock has occurred, these are both librating about zero. The resonant interaction causes an eccentricity growth of both planets, the growth halting at around average values of e1 = 0.06 and e2 = 0.34 although both eccentricities exhibit variations around these average values. The periastron angles of the two planets oscillate around the alignment position, being the natural stable state, and this is in agreement with model parameter fits to observations (Laughlin & Chambers 2001).

The disc cavity remains at low density, and the cavity edge diffuses slowly in on the viscous diffusion timescale, following the outer planet. Figure 2 shows a surface density plot of the inner regions of the disc.

3.3 Relation to analytic model and observations

The simulation confirms that tidally induced migration of two planets with the given mass ratios leads to trapping into the 2:1 commensurability. The resonance angles both undergo stable periodic librations as per Sect. 2 and the resonance is robust for the length of the simulation. We can conclude that stable resonance trapping would be the probable outcome of the evolution of such a system. The finding of the system GJ876 near to or in such a resonance is probably simply a consequence of past tidally induced migration of the two planets into such a state. The magnitudes of the eccentricities are in broad agreement with fitted model parameters to the observations of GJ876 with $e_{1} = 0.12 \pm 0.10$ and $e_{2} = 0.31 \pm 0.08$ (Marcy et al. 2001; Laughlin & Chambers 2001). The eccentricity ratio $e_{2} / e_{1} \approx 5.7.$ We comment that Eq. (11) which applies in the small eccentricity limit when one of $\phi, \psi$ is zero and the other $\pi$gives $e_1/e_2 \sim 1/11.$However, the fact that both angles librate about zero indicates, within the context of the simple model, that non resonant orbital precession needs to be incorporated and we should use Eq. (10) to make a comparison with the simulation. This is

$\displaystyle e_1 e_2 a_1 M_* \omega_{pr1}+ m_2 a_1 e_2 n_1 B \cos \phi = - m_1 a_2 e_1n_2 C\cos\psi +e_1 e_2 a_1 M_* \omega_{\rm pr2}.$     (19)

Since the librations of the angles are about a value close to zero, we replace them by zero to make simple estimates. The largest non central mass in the system is m1 so we include its effect in causing non resonant orbital precession for m2by including a non zero $\omega_{\rm pr2}.$Also we include the corresponding effect on m1 due to m2which produces a non zero value of $\omega_{\rm pr1}.$Then we have

\begin{displaymath}% \omega_{\rm pr1}+ { m_2 n_1 B\over e_1 M_*} = \omega_{\rm pr2}- {n_2 m_1a_2 C\over e_2 a_1 M_*}\cdot \end{displaymath} (20)

This condition equates the precession rate of the orbit of m2on the right with that of m1 on the left. In the former case the first contribution comes from the time averaged orbit and the second is the resonant contribution. In fact with $e_2 \sim 0.34,$ the orbits approach each other quite closely to within 22 percent of the outer semi major axis. But the resonant configuration avoids close encounters keeping the planets apart. Thus we expect the resonant and non resonant contributions to the precession rate to show significant cancellation. The large value of e2 probably makes strict comparison inaccurate. Nonetheless the simulation gives a precession period of about 70orbits in a retrograde sense. We can estimate that this to be of the same magnitude as would be predicted from the second term on the left hand side which however gives prograde precession. Thus the non resonant precession rate due to m2 needs to be be comparable to the resonant effect but of opposite sign.

Considering the precession rate of m2, we estimate that $\omega_{\rm pr2}$ corresponds to 70 orbits in the prograde sense while the resonant term corresponds to 70/2 orbits in the retrograde sense. Thus the combined effect produces a precession period of 70 orbits in the retrograde sense as seen in the simulation.

In summary our simulation gives plausible eccentricity values for the two planets, that can be understood in outline by use of a simplified analytic theory and are consistent with the current observations.

  
4 Orbit integrations

In addition to the analytic model presented in Sect. 2 and the hydrodynamic simulations presented in Sect. 3, we have also performed three-body orbit integrations using a fifth-order Runge-Kutta scheme (e.g. Press et al. 1993).

The basic assumptions of the model are that the two planets exist within the inner cavity of a tidally truncated disc that lies exterior to the outer planet. Tidal interaction with this disc causes inwards migration of the outer planet, and also leads to eccentricity damping of the outer planet. It is further assumed that as the planets migrate inwards and approach their final semi-major axes, the disc material disperses exponentially on some prescribed decay time $t_{\rm disp}$, the process being initiated such that the required final semi-major axis of the outer planet was asymptotically attained. In our numerical calculations, a torque was applied to the outermost planet such that it migrated inwards on a time scale of $t_{\rm mig}$ local orbital periods as defined in Sect. 2, and a damping force was applied in the radial direction to damp the eccentricity on a time scale of $t_{\rm c}$ local orbital periods also as defined in Sect. 2.

 

 
Table 1: This table shows the parameters used for the orbital integration runs. The first column gives the run label, the second column the migration time scale, the third column the circularization time scale, and the fourth column gives the disc dispersal time. The fifth, sixth and seventh columns give the final eccentricities for the outer and inner planets, and the ratio of these eccentricities.
Run $t_{\rm mig}$ $t_{\rm c}$ $t_{\rm disp}$ e1 e2 e2/e1
R1 $3.33\times 10^3 \;\;\;$ 360 $8 \times 10^3 \;\;$ 0.095 0.41 4.3
R2 $3.33\times 10^3 \;\;\;$ 100 $8 \times 10^3 \;\;$ 0.05 0.3 6
R3 $ 3.33 \times 10^2 \;\;\;$ 360 $8 \times 10^2 \;\;$ 0.34 0.72 2.1
R4 $ 10^3\;\;\;$ 108 $ 2.4 \times 10^3$ 0.095 0.41 4.3
R5 $3.33\times 10^3 \;\;\;$ 360 $0.0 \;\;\;\;\;$ 0.095 0.41 4.3


These integrations used initial conditions corresponding to the more massive, outermost planet being located initially at 5 AU, with the lighter inner-most planet located initially at 2.5 AU. The planet masses adopted for the orbit integrations are the same as the minimum masses reported for the planets in the system around GJ876 by Marcy et al. (2001) (i.e. 1.87 $M_{\rm J}$and 0.56 $M_{\rm J}$). The stellar mass is taken to be 0.32 $M_{\odot}$. Whilst these calculations provide only a crude approximation to the detailed physics of disc-companion interactions, their simplicity allows us to perform many calculations, covering a wide area of parameter space, and also to run for much longer time scales than is possible for simulations of the type described in Sect. 3.

A number of calculations have been performed to examine the relationship between the final values of e1, e2, and their ratio e1/e2, to the various input parameters $t_{\rm mig}$, $t_{\rm c}$, and $t_{\rm disp}$. The results of some of these calculations are presented in Table 1, and are discussed below. The unit of time used in the abscissa of the Figs. 3 to 5 is the orbital period for an object at 1 AU in orbit around a star with mass 0.32 $M_{\odot}$, and is denoted as P(1 AU).

  \begin{figure} \par\includegraphics[width=8.2cm,clip]{MS1339f3.ps}\end{figure} Figure 3: This figure shows the evolution of the planet semi-major axes and eccentricities for the run R1 shown in Table 1.
Open with DEXTER

4.1 Dependence on migration and circularization times

Equation (18) shows that the eccentricity of the outer planet, e1, depends on the ratio of $t_{\rm c}/t_{\rm mig}$. Here we present results of simulations that explore how the eccentricity ratio e2/e1 depends on $t_{\rm c}$ and $t_{\rm mig}.$

Figure 3 shows the evolution of the semi-major axes and eccentricities for the run R1, whose model parameters are described in Table 1. This figure shows the inward migration of the outer planet that subsequently locks to the inner planet as it reaches the 2:1 commensurability. The subsequent evolution is such that the two planets, now resonantly locked, migrate inwards. The eccentricities result from the balance between eccentricity driving through the resonant interaction, and eccentricity damping due to the disc interaction. As the planets approach their final semi-major axes, the effects of migration and eccentricity damping are removed exponentially with decay time $t_{\rm disp}=8 \times 10^3$ local orbits, causing them to cease migration at semi-major axes $a_1 \simeq 0.2$ and $a_2 \simeq 0.126$ which are values appropriate to GJ876. The subsequent evolution beyond a time of $t \simeq 0.45$ P(1 AU) in Fig. 3 occurs in the absence of disc effects, and suggests a long-term stability of the system given that it remains stable for $2 \times 10^{7}$ orbits of the outer planet at its final semi-major axis. Figure 4 shows the evolution of the resonant angles $\phi $ and $\psi $, which librate about $\phi =\psi =0$ in agreement with the results presented in Sect. 3 for the full hydrodynamic simulations.

Figure 5 shows the results from run R3 described in Table 1, and illustrates the effect of reducing $t_{\rm mig}$ while keeping $t_{\rm c}$ constant. As expected from Eq. (18), the eccentricity of the outer planet increases from $e_1 \simeq 0.1$ to $e_1 \simeq 0.35$ when $t_{\rm mig}$changes from $ 3.33\times10^3$ to $3.33\times10^2$ local orbital periods. However, the ratio e2/e1 has also changed from $e_2/e_1 \simeq 4$ to $e_2/e_1 \simeq 2$.

Calculation R2 illustrates the effect of keeping $t_{\rm mig}$ constant while changing $t_{\rm c}$. As expected from Eq. (18), reducing $t_{\rm c}$leads to a reduction in e1, since the disc model damps the eccentricity more effectively. We also find that the ratio e2/e1 again changes, it now being $e_2/e_1 \simeq 6$, as compared to $e_2/e_1 \simeq 4$ for run R1.

  \begin{figure} \par\includegraphics[width=8cm,clip]{MS1339f4.ps}\end{figure} Figure 4: This figure shows the evolution the resonant angles $\phi $ and $\psi $for the run R1 described in Table 1. These show libration around $\phi =\psi =0$, in agreement with the results from the hydrodynamic simulation shown in Fig. 1.
Open with DEXTER

Equation (18) predicts that keeping the ratio $t_{\rm mig}/t_{\rm c}$ constant, but changing both $t_{\rm mig}$ and $t_{\rm c}$ independently, should leave e1unchanged. Calculation R4 indicates that this is what happens, and also shows that the ratio e2/e1 remains unchanged.

Overall the results are entirely consistent with the analytic predictions presented in Sect. 2 and with the hydrodynamical simulations presented in Sect. 3. Furthermore, they indicate that the ratio e2/e1 scales rather weakly with $t_{\rm mig}/t_{\rm c}$. We comment that from Sect. 2 Eq. (10) we expect the eccentricity ratio to reach ${\sim} 11$ as $e_1 \rightarrow 0$. Long-term stability of two-planet systems that become locked due to disc-induced orbital migration is also indicated by our calculations. In particular run R1 covers a time scale corresponding to $2 \times 10^{7}$ orbits of the outer planet in its final configuration. We find that it is possible to arrange e1to match the observed value of the outer planet in the GJ876 system by fixing $t_{\rm mig}$ and choosing $t_{\rm c}$ appropriately, but it is difficult to then obtain a value for e2 that matches the reported value of e2=0.27.

4.2 Dependence on disc dispersal time scale

We have performed simulations to examine the effect of removing the disc on different time scales on the stability of the system. We find that the stability is largely unaffected by the rate at which the disc is removed. Figure 3 shows the evolution of the semi-major axes and eccentricities from a run in which the disc was removed on a time scale of $t_{\rm disp}=8 \times 10^3$ local orbital periods. In this case the disc dispersal was switched on once the outer planet semi-major axis was a1 < 0.5 AU. The final semi-major axis of the planet is determined by the radius at which the disc dispersal is initiated, and the time scale over which the disc is removed. Calculation R5 was similar to R1 except that the disc was removed instantaneously. The results presented in Table 1 show that this has little effect on the final outcome. A slight increase in the scatter of the temporal distribution of eccentricities was observed.

  \begin{figure} \par\includegraphics[width=8.2cm,clip]{MS1339f5.ps}\end{figure} Figure 5: This figure shows the evolution of the planet semi-major axes and eccentricities for the run R3 shown in Table 1.
Open with DEXTER

   
5 Summary and discussion

In this paper we have considered two protoplanets gravitationally interacting with each other and a protoplanetary disc. The two planets orbit interior to a tidally maintained disc cavity while the disc interaction induces inward migration.

We have supposed that the previous evolution of the system results in the the planets getting into a configuration with an orbital separation just exceeding that required for a 2:1 commensurability. This evolution is likely to have involved both accretion and migration. Given the set up we consider a natural evolution is towards both planets orbiting in a cavity outside of which orbits the protoplanetary disc. Tidal interaction results in the inner disc material either being expelled into the outer disc or accreting onto the central star. Subsequently the outer planet migrates towards the inner one as a result of interacting tidally with the exterior disc.

When the migration is slow enough, we found that the outer protoplanet approached and became locked into a 2:1 commensurability with the inner one. This was maintained in subsequent evolution. We studied the nature of these interactions using a simple analytic model, hydrodynamic 2D simulations and longer time N-body integrations. These all gave consistent results.

The magnitude of the stabilized eccentricities was found to be determined by the ratio of the migration rate to the circularization rate induced in the outer planet orbit by the external disc. The eccentricity ratio e1/e2 was found to vary with the magnitude of the ratio being more extreme for smaller eccentricities. We have applied our results to the recently discovered resonant planets around GJ876. Simulation shows that a disc with parameters expected for protoplanetary discs causes trapping in the 2:1 commensurability when the planets orbit in an inner cavity and that eccentricities in the observed range may be obtained. In this case the orbits were found to be aligned with both resonant angles librating about zero. The whole system then precessed in a retrograde sense. Finally when the disc is removed on a range of timescales the orbital configuration has been found to be stable for up to $2 \times 10^{7}$ orbits subsequently.

Acknowledgements
This work was supported by PPARC grant number PPARC GR/L 39094 and the computational resources of the GRAND HPC consortium. We thank Udo Ziegler for making a FORTRAN Version of his code NIRVANA publicly available.

References

 

Copyright ESO 2001