A&A 483, 325-337 (2008)
DOI: 10.1051/0004-6361:20079291

Influence of an inner disc on the orbital evolution of massive planets migrating in resonance

A. Crida¹ - Z. Sándor^2,3 - W. Kley¹

1 - Institut für Astronomie & Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
2 - Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
3 - Department of Astronomy, Eötvös Loránd University, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

Received 20 December 2007 / Accepted 13 February 2008

Abstract
Context. The formation of resonant planet pairs in exoplanetary systems involves planetary migration inside the protoplanetary disc. After a resonant capture, subsequent migration leads to a large increase in planetary eccentricities, if no damping mechanism is applied. This has led to the conclusion that the migration of resonant planetary systems cannot take place across large radial distances, but must be terminated rapidly by disc dissipation.
Aims. We investigate if the presence of an inner disc could supply eccentricity damping to the inner planet, and if this effect could explain observed eccentricities in some planetary systems.
Methods. We compute hydrodynamic simulations of giant planets, in orbits of a given eccentricity about an inner gas disc, and measure the effect of the gas disc on the planetary orbital parameters. We perform detailed long-term calculations of the GJ 876 system. We complete N-body simulations, which include artificial forces on the planets that recreate the effect of the inner and outer discs.
Results. We find that we cannot neglect the influence of the inner disc, and that the disc could explain the observed eccentricities. In particular, we reproduce the orbital parameters of a few systems engaged in 2:1 mean motion resonances: GJ 876, HD 73 526, HD 82 943, and HD 128 311. Analytically, we derive the effect that the inner disc should have on the inner planet to reach a specific orbital configuration, for any given damping effect of the outer disc on the outer planet.
Conclusions. We conclude that an inner disc, even though difficult to model properly in hydrodynamical simulations, should be taken into account because of its damping effect on the eccentricity of the inner planet. By including this effect, we can explain quite naturally the observed orbital elements of the pairs of known resonant exoplanets.

Key words: accretion, accretion discs - planets and satellites: formation - celestial mechanics

  
1 Introduction

The orbital evolution of a system consisting of young protoplanets is governed by disc-planet and mutual gravitational interactions. In the case of differential migration, the semi-major axis ratio of two planets varies with time and - in the situation of convergent migration - capture in a resonant configuration can occur. A large fraction of the observed multi-planet systems, contain a pair of planets engaged in a resonance. Here, we are interested in mean motion resonances (MMR), where the ratio of the (mean) orbital periods of the outer to the inner planet, equals that of two small integers. Among the 6 systems known to contain an MMR, 4 have a ratio of 2:1: GJ 876, HD 82 943, HD 128 311, and HD 73 526. The system first discovered to be in a 2:1 configuration (GJ 876) is interesting in several aspects. The planets are both massive (0.56 and 1.94 Jupiter mass), in contrast to the small mass of the central star, which is only 0.33 $M_\odot$. The short orbital periods of the planets ($\approx$30 and $\approx$60 days) allow the accurate determination of their orbital elements, which are provided in Table 1. In GJ 876, the two outer planets are``deep'' in the 2:1 MMR: the apsidal lines of the two osculating, orbital ellipses, are always aligned and librate with small amplitudes only (so-called apsidal resonance or corotation). As a consequence of the apsidal resonance, the planetary eccentricities show only small variation with time. A resonant configuration such as that within GJ 876, can be established only by the action of dissipative effects, such as disc-planet interaction. The mere existence of systems engaged in MMRs, is one of the strongest indication that planetary migration has indeed occurred, during the early evolution of planetary systems.

The first detailed modelling of GJ 876 was conducted by Lee & Peale (2002), who performed customised 3-body simulations of a central host star and two planets, with additional (dissipative) forces that reproduced the effects of disc-planet interaction. In such simulations, it is assumed that a pair of planets is embedded in the protoplanetary disc; this disc consists of an outer disc only, while the inner disc, which is inside the orbit of the inner planet, has already been lost due to accretion onto the star and planets, or final evaporation. In such a configuration, only the outer planet is in contact with an even more distant protoplanetary disc, and experiences typically negative Lindblad torques, and migrates inward. In contrast, the inner planet has no ambient material and does not feel any disc torque. In terms of the 3-body simulations by Lee & Peale (2002), this implies that additional forces reduce the semi-major axis and eccentricity of the outer planet, while the inner planet feels only the direct gravitational forces of the star and outer planet. Capture into a resonant configuration can occur when, during the inward migration, the outer planet crosses the location at which the mean orbital periods have a ratio of two small integers.

The eccentricity damping action of the outer disc, onto the outer planet, is typically parameterised by the migration rate, i.e.

$$\frac{\dot{e}}{e} = - K \left \vert\frac{\dot{a}}{a} \right \vert$$ (1)

where a and e are the semi-major axis, and eccentricity of the planet, respectively, while the constant factor K relates the damping rates. For low-mass planets in the linear regime, it is known that the eccentricity damping timescale, $\tau_e = -e/\dot{e}$, is about h² times shorter than the migration timescale $\tau_a = -a/\dot{a}$ (Tanaka & Ward 2004), where h is the aspect ratio of the disc (generally a few percent). This short timescale was confirmed using fully-nonlinear, hydrodynamical simulations in two and three dimensions (Cresswell et al. 2007). Hence, for low-mass planets one may safely assume that they migrate inward, on nearly circular orbits, unless disturbed by additional bodies in the system. For higher planetary masses, the value (and also the sign) of K, is not precisely known; due to the gap opening, it is however expected that the damping on the eccentricity will be reduced, with respect to the linear case. The simulations by Lee & Peale (2002) indicate that for values of K equal to unity, the resonant action leads to a strong growth of eccentricity of both the inner and outer planets, much larger than the observed values ( $e_1 \approx 0.03$ and $e_2 \approx 0.22$). Lee & Peale (2002) pointed out that the final state of the system is determined by K. For a vanishing K, i.e. no eccentricity damping, the planetary eccentricities continue to grow. To match the observations, Lee & Peale (2002) had to assume a value of K=100, a value that appears to be large, considering the high masses of the planets.

The convergent migration of two massive planets was demonstrated in a variety of hydrodynamical simulations (Kley 2000; Snellgrove et al. 2001; Papaloizou 2003; Bryden et al. 2000). Using multi-dimensional hydrodynamical simulations of resonant planetary systems, it was shown that for masses in the Jupiter regime, the value of K is approximately 1-10 at most (Kley et al. 2004). In a detailed study of the system GJ 876, Kley et al. (2005) have shown that, in hydrodynamical simulations where the inner disc was depleted, the final eccentricities of the planets were always much larger than those observed, unless one assumes that the outer disc dissipates rapidly, on the viscous time scale. Hence, this scenario does not allow for the migration of the resonant planets over a large radial distance.

In hydrodynamical simulations of discs with embedded planets, it is often found that the inner disc is depleted rapidly; this may not however always be the case, and may instead be an artefact of inappropriate inner boundary conditions. As shown by Crida & Morbidelli (2007), an inner disc will survive for much longer than believed previously, even in the presence of a planet. In this situation, the inner disc should have a dynamical influence on the inner planet, and induce possibly some additional damping of the eccentricity. The effect of such eccentricity damping on the inner planet, was described first by Lee & Peale (2002), for the particular case of GJ 876. They found that a value of K=10, for both planets, yields final eccentricities in the observed range. The first full hydrodynamical study in this direction - including an inner disc - was completed for the resonant system HD 73 526, by Sándor et al. (2007), in their study, they showed that the inclusion of an inner disc leads to an eccentricity damping of the inner planet, and allows more extended radial migration with reasonable final eccentricities.

In this paper, we analyse the effect in more detail, and investigate the dynamical influence of an inner disc on a planet. In Sect. 2, we perform a sequence of hydrodynamical simulations, measure the torque and power exerted by the disc on the planet, and evaluate the change in eccentricity and semi-major axis as a function of the planetary eccentricity. In Sect. 3, we complete a full time-evolution of a pair of planets, embedded in a protoplanetary disc, for realistic parameters, with specific application to GJ 876. We show that the torque and power generated by the inner disc, yield an effective damping of e that results in moderate final eccentricities, even for extended radial migration. In Sect. 4, we apply an eccentricity damping to the inner planet, in N-body simulations, with artificial non-conservative forces that reproduce the effect of the disc. We apply this analysis to a few exoplanetary systems, and try to recover the observed orbital elements with a realistic migration scenario. Section 5 summarises our results.

   
2 Effect of a gaseous inner disc on the orbital elements of a planet on an eccentric orbit

We perform a suite of hydrodynamical simulations to measure the influence of an inner disc on the orbital elements of a planet, which has an eccentric orbit. The disc is treated as a two-dimensional gas that occupies the orbital plane of the planet. The disc material is present only inside the planetary orbit, so that the effect of the inner disc can be isolated. The planet has a mass $M_{\rm p}=2.14\times 10^{-3}\ M_*$, and a fixed orbit with semi-major axis a=1, and an eccentricity e that is varied between one simulation and another.

The disc is treated as a non-self-gravitating gas that can nevertheless interact gravitationally with the planet. Using the expressions for the planetary energy $E=-G(M_*+M_{\rm p})M_{\rm p}/2a$ and angular momentum $H=M_{\rm p}\sqrt{G(M_*+M_{\rm p})a(1-e^2)}$, one can derive:

  

$$\dot{a}/a = -\dot{E}/E$$ (2)

$$\dot{e}/e = \frac{e^2-1}{2e^2}\left(\frac{\dot{E}}{E} + 2\frac{\dot{H}}{H}\right)\cdot$$ (3)

In the present simulations we fix the orbit of the planet, but monitor the torque ($\dot{H}$) and power ($\dot{E}$) acting on the planet, averaged over one orbit. We check that if the planet is released from its fixed orbit, its eccentricity and semi-major axis follow the expected evolution.

We use the FARGO code developed by Masset (2000b,a), which is a two-dimensional hydro-code in cylindrical coordinates ($r,\varphi$), with an isothermal equation of state. Thus, the sound speed is given by $c_{\rm s}=hr\Omega$, where $\Omega$ is the local angular velocity, r is the distance to the central star, and h is the aspect ratio, which is here 0.05. The gas viscosity $\nu$ is given by an $\alpha$-prescription ( $\nu=\alpha c_{\rm s} h r$, Shakura & Sunyaev 1973), with $\alpha=10^{-2}$. The grid covers the region from 0.4 to 1.62 in radius, divided into 112 elementary rings (logarithmically spaced), themselves divided into 500 sectors. The inner boundary condition is non-reflecting; at every time step, the density in the zeroth ring is set to be equal to the density of the first ring, which is rotated to simulate wave propagation and avoid wave reflection; the density of the zeroth ring is then shifted to ensure that its azimuthally-averaged density is similar to its initial value. The outer boundary is open, which means that outflow of gas out of the grid is permitted.

	Figure 1: Gas density profiles after 400 orbits for a few planet eccentricities, and common initial profile (plain line).
Open with DEXTER

The coordinate system is centred on the star. It can be non-rotating, corotating with the planet, or rotating at a constant angular velocity with a period equal to that of the planet

$$\Omega_0 = \left(\frac{G(M_*+M_{\rm p})}{a^3}\right)^{1/2}\cdot$$ (4)

The choice of a different coordinate rotation-rate, can alter the gas dynamics in the vicinity of the planet. This may occur even for the FARGO algorithm, which treats each ring in a corotating frame. We perform comparison simulations in all three frames and show our results. In the corotating frame, gas in the Hill sphere rotates about the planet, which moves radially from apoastron to periastron and vice versa. In the non-rotating frame, the motion of the planet about the star in the grid appears to spread the Hill sphere. For large eccentricities, the co-rotating frame is disadvantageous because it does not rotate at a constant angular velocity; this implies that additional terms exist in the equations of motion. We therefore apply the second case, which is a frame rotating with a constant speed $\Omega _0$, within which the planet describes an epicycle.

The initial density profile corresponds to an approximate gap opened by the planet, without the outer disc, such that an equilibrium profile and a stationary regime are quickly reached. The initial profile is plotted as a solid line in Fig. 1. The initial total mass of gas in the disc is $1.53\times 10^{-3} M_*=71\%~M_{\rm p}$. For various planet eccentricities, we display the final density profiles, after 400 orbits, in Fig. 1. The density measurements are those when the planet is at apoastron, and the density spike at r=1+e corresponds to the Hill sphere of the planet. As the eccentricity increases, the inner disc becomes more depleted. The outer disc is clearly empty. In each of the five cases considered, the simulation was completed inside a non-rotating frame.

2.1 Average effect of the inner disc

The power and torque from the disc on the planet become almost constant after about 200 orbits, and we measure both quantities after 400 orbits. To measure the force exerted by the gas on the planet, we exclude material inside the Hill sphere, using a tapering function given by:

$$f({ d})=\left[\exp\left(-\frac{{ d}/r_{\rm H}-p}{p/10}\right)+1\right]^{-1}$$ (5)

where d is the distance to the planet, $r_{\rm H}$ is its Hill radius $r_{\rm H}=r_{\rm p}(M_{\rm p}/3M_*)^{1/3}$ ($r_{\rm p}$ being the distance between the planet and the star), and p is a dimensionless parameter set to be equal to 0.8. This function is plotted in Fig. 2 for different values of p. We stress that the use of such a tapering is important: the gas close to the planets exerts a strong force, which should not be taken into account because this gas is gravitationally bound to the planet, and should be considered as a part of the planet mass. With self-gravity, this gas should feel the same force as and naturally follow the planet; tapering frees the planet from the need to carry artificially this material. The influence of the shape of the tapering function will be discussed below.

	Figure 2: Tapering function f( d) defined by Eq. (5) for four values of the parameter p.
Open with DEXTER

For the three above mentioned options for the rotation rate of the coordinate system, the results are shown in Fig. 3. The quantities are normalised by $\mu=M_{\rm disc}/M_{\rm p}$, because the force felt by the planet is proportional to the gas density. The two top panels display the influence of the inner disc on the orbital elements e and a. For $e\geqslant 0.1$, the eccentricity is effectively damped on a timescale of about two thousands of orbits, as can be seen in the bottom left panel displaying $\tau_e=-e/(\dot{e}/\mu)$; the dispersion in the values of $\tau_e$ for $0.1\leqslant e \leqslant 0.15$ is due to the fact that $\dot{e}$ is close to zero. This confirms the idea presented by Lee & Peale (2002) and Sándor et al. (2007) that an inner disc has a damping effect on the planetary eccentricity. For low eccentricities ( $e\leqslant 0.1$), the influence of the disc is small ( $\vert\dot{e}\vert<4\times 10^{-5}$ orbit^-1), and could even lead to a small excitation of the eccentricity.

	Figure 3: Influence of the inner disc on a planet on a fixed orbit, as a function of the eccentricity. NR (+ symbols): in the non-rotating frame; CF ($\times$ symbols): in the corotating frame; $\Omega$F (open squares): in the frame rotating monotonically with the same period as the planet.
Open with DEXTER

It is well-known that a planet on a circular orbit exerts a negative torque on the inner disc (Goldreich & Tremaine 1980; Lin & Papaloizou 1979). It repels the gas, leading to the opening of a gap. When the density gradient at the gap edge is sufficiently steep, an equilibrium is achieved (see for instance Crida et al. 2006). The planet then feels a positive torque from the inner disc. For e=0, we find that the planet feels a positive torque (which is equal to the power), as expected; consequently, $\dot{a}>0$. But as e increases, the power decreases, and for $e\gtrsim 0.3$, the semi-major axis variation expected if one releases the planet is negative (top right panel). The inner disc appears to attract the planet more and more as the eccentricity grows.

The variations of $\dot{a}$ and $\dot{e}$ with e, show that the simple model represented by Eq. (1), for which K is a constant, is an over simplification. The bottom right panel of Fig. 3 displays $(\dot{e}/e)/(\dot{a}/a)=K$ as a function of e, which varies significantly. The dispersion about e=0.3-0.35 is due to the change of sign of $\dot{a}$. For 0.1<e<0.25, a K value of minus a few appears reasonable, but a precise value cannot be stated. However, it is clear that the effect of the inner disc on the planet cannot be neglected. The left panels in particular show that a significant damping of the eccentricity has to be taken into account.

2.2 More detailed analysis

To investigate the physical mechanism more thoroughly, we focus on the case where e=0.15 and where the frame is rotating at a constant angular velocity $\Omega _0$. In Fig. 4, we plot the specific torque and power felt by the planet during one orbit, starting at apoastron. The units are normalised such that a=G=M_*=1, and thus $\Omega_0 = \sqrt{1.00214}$. Four curves are presented, corresponding to measurements during the same orbit for four different values of the parameter p in Eq. (5), in the range [0.6;0.9]. The larger the p factor becomes, the smaller is the amplitude of the oscillation. This oscillation is due to gas in the Hill sphere. For each case, the average is plotted as a straight horizontal line. Both the mean torque and power vary by 10 percent for different values of p. The corresponding $\dot{a}$ also varies by 10 percent, while $\dot{e}$ varies by about $5\%$.

A fifth curve, labelled p=0, is drawn on each panel of Fig. 4, that corresponds to the case where f=1 (no tapering). The measured torque and power differ and are quite noisy with respect to the cases with tapering. This is because the density structure is non-stationary, even within the Hill sphere, due to the eccentric orbit. The deviations are so large that the mean values for the p=0 case, averaged over one orbit, give $\dot{a}/\mu= -10^{-3}$, instead of $2.8\times 10^{-4}$ for p=0.8, and $\dot{e}/\mu= -7.14\times 10^{-4}$, instead of - $7.95\times 10^{-5}$. This clear difference of the p=0 case, and the good agreement between the other four cases $p\in[~0.6;~0.9~]$, highlights the importance and necessity of tapering.

	Figure 4: Specific torque ( left panel) and power ( right panel) acting on the planet during one orbit for a fixed semi-major axis and eccentricity, e=0.15. The different curves correspond to measures using different values of the p parameter in Eq. (5). The horizontal lines correspond to the average in each of the four cases $p\neq 0$.
Open with DEXTER

Figure 5: Evolution during one orbit of the vectors ${\vec r}_{\rm p}$: position of the planet in the frame rotating at velocity $\Omega _0$ ( right); ${\vec v}_{\rm p}\vert _{\rm NR}$: velocity of the planet in the non-rotating frame ( top); ${\vec v}_{\rm p}\vert _{\Omega\rm F}$: velocity of the planet in the frame rotating at constant velocity $\Omega _0$ ( middle); $\vec{F}$: force from the disc on the planet ( left, arbitrary unit).

Open with DEXTER

	Figure 6: Density maps, in linear grey scale, in the frame rotating at constant velocity $\Omega _0$ at four different moments of the orbit corresponding to Fig. 4.
Open with DEXTER

For the frame rotating at constant angular velocity $\Omega _0$, we draw, in Fig. 5, the traces during one orbit of the vectors ${\vec r}_{\rm p}$, ${\vec v}_{\rm p}\vert _{\rm NR}$, ${\vec v}_{\rm p}\vert _{\Omega\rm F}$, and $\vec{F}$, which are respectively the position of the planet, its velocity in the non-rotating and rotating frames, and the force it feels from the disc, in normalised units. The force is measured for p=0.8 and is magnified by a factor $5\times 10^6$ to improve visibility. The vectors are drawn at the beginning of the orbit; a cross symbol is then drawn for each twentieth of the orbit. In the frame rotating at the velocity $\Omega _0$, the planet describes an epicycle centred on (x=1, y=0), in the clockwise direction starting at (x=a+e=1.15, y=0). The planet velocity marks out an epicycle, centred on (v_x=0, $v_y=a\Omega_0)$, in the anti-clockwise direction, starting at (v_x=0, v_y<1); the velocity in the rotating frame ${\vec v}_{\rm p}\vert _{\Omega\rm F}$ is plotted as a dashed line, describing a curve about zero in the clockwise direction. As could be expected, the force experienced by the planet is directed towards the inner disc (F_x<0), in the direction of the wake. The wake tends to rotate at a constant velocity. The wake leads the planet, and $F_{\rm y}$ is positive, between apoastron and periastron; this occurs because the planet is late with respect to the average rotation at speed $\Omega _0$ ($y_{\rm p}$is negative). Following periastron, the planet is ahead the root of wake, and $y_{\rm p}>0$, and F_y<0. This can be seen in the density maps of Fig. 6. This explains the variations of the sign of the power ${\vec v}_{\rm p}\vert _{\rm NR}\cdot\vec{F}$ observed in Fig. 4. The torque ${\vec r}_{\rm p}\wedge\vec{F}$ is always positive because the angle ( ${\vec r}_{\rm p}$, $\vec{F}$) remains slightly smaller than $\pi$, as can be seen in the figure.

2.3 Discussion

In this section, we consider one particular case, for fixed planet mass and disc parameters. The width, the depth, and the shape of the gap, opened by the planet would change if the disc viscosity, aspect ratio, or the planet mass vary, and the effect of the inner disc would consequently differ. The smoothing length for the potential, the tapering, and the choice of the frame can affect the measured force of the gas on the planet. As seen from the right panel of Fig. 4 a small change in the power, e.g. due to numerical aspects, can lead to a sign reversal of the average value, which would imply a reversed migration. Exploring all possible parameter space is beyond the scope of this paper. We want to demonstrate that the inner disc has a non negligible damping effect on eccentric giant planets, and we provide a qualitative explanation of this phenomenon.

The mass of the disc decreases, while evolving to an equilibrium state. At the end of the simulations the disc mass is between 30 and 50 percent of the planet mass in all cases.

We investigate whether eccentricity excitation can be induced by the planet in the disc via a tidal instability mechanism, which would operate through an inner 3:1 Lindblad resonance in the disc (Lubow 1991). If the planet is orbiting inside an outer disc, and for sufficiently large planet masses ( $M_{\rm p} \gtrsim 3~M_{\rm Jup}$), the outer disc will turn eccentric even in case of a planet on a circular orbit (Kley & Dirksen 2006; Papaloizou et al. 2001). The strength of such a mechanism scales inversely with planet mass such that for the smallest planet masses the timescale of eccentricity growth is several thousand planetary orbits. The gap created by the planet must be sufficiently wide for the instability to work; the outer 2:1 Lindblad resonance, which damps the eccentricity, must be cleared. Adding a small planetary eccentricity will reduce the necessary planet mass slightly (D'Angelo et al. 2006).

To check for this effect in our simulations, we studied models for e=0.15 for over 3500 orbits; we did not observe any eccentricity growth in the disk, even when the planet mass was increased to $5~M_{\rm Jup}$. We conclude that, for our planet masses, we do not expect to observe an eccentric inner disc. When we measure the torque and power of the inner disc force on the planet, our simulations have already reached a steady state; we therefore expect to observe an eccentricity damping of the planet by the almost circular disc. As pointed out above, we checked this behaviour using simulations where we release the planet from its fixed orbit and follow its orbital evolution self-consistently. We find exactly the results predicted for variations in the a and e of the planet.

Due to the small planet mass and the appreciable viscosity and pressure, the absence of an eccentric inner disc could be due to the lack of clearance of the eccentricity damping resonances in the disc, even for large planet eccentricities.

  
3 Application to the GJ 876 system

3.1 Characteristics of the system

Table 1: Parameters of the system GJ 876 as given by Laughlin et al. (2005) for GJ 876 b and GJ 876 c, and by http://exoplanets.eu/ for GJ 876 d.

GJ 876 is a $0.32~M_\odot$, M4 star that hosts 3 planets, the two largest of which are in a 2:1 Mean Motion Resonance (MMR); see Table 1. The resonant angles that characterise this resonance are $\theta_1=\lambda_c-2\lambda_b+\omega_c$, $\theta_2=\lambda_c-2\lambda_b+\omega_b$, and $\Delta\omega=\omega_b-\omega_c$, where $\lambda$ is the mean longitude and $\omega$ is the longitude of the periastron. All three angles librate about $0\hbox{$^\circ$ }$. Laughlin et al. (2005) estimated the amplitude of the librations to be: $\vert\theta_1\vert _{\rm max}=7\pm 1.8 \hbox{$^\circ$ }$, $\vert\theta_2\vert _{\rm max}=34\pm 12\hbox{$^\circ$ }$, and $\vert\Delta\omega\vert _{\rm max}=34\pm 12\hbox{$^\circ$ }$.

As for all hot Jupiters, the two giant planets in GJ 876 most likely formed further out in the protoplanetary disc (beyond the snow-line, where water is solid and can contribute to the formation of a massive core), and then migrated toward the central star. Migration is supported by the resonant motion: such a configuration can only be achieved by convergent motion of the planets. Thus, the outermost planet must have migrated inward more rapidly than the innermost one, and eventually captured it in resonance; or both planets may have been orbiting in a common gap, in which they approach each other under the repulsive action of the discs, until a resonant capture. They may then have migrated together to achieve their current configuration. Dissipation in the gas disc can also modify the amplitude of libration of the resonant angles.

The migration of two giant planets in MMR, in a gas disc, was studied first by Masset & Snellgrove (2001) for the case of Jupiter and Saturn, and in more detail by Morbidelli & Crida (2007). They showed that if the outermost planet is significantly lighter than the innermost one, the migration of both planets may be stopped, or reversed outward. In the case of GJ 876, the outer planet is substantially more massive than the inner one, and we would expect the inward migration of both planets. The observed small semi-major axes of GJ 876 b and GJ 876 c, are therefore not a surprise. Having been pushed inward by the outer planet (GJ 876 b), the innermost of the two resonant planets (GJ 876 c), should have had its eccentricity raised dramatically; this is according to the results of N-body simulations with artificial migration rate, and detailed hydrodynamical simulations, which were discussed in the Introduction. But we have seen that an inner gaseous disc has a damping effect on the eccentricity of a giant planet, at least for $e\geqslant 0.1$. Thus, we suggest that the presence of a gaseous disc inside the orbit of GJ 876 c could prevent its eccentricity from reaching unreasonably large values. To verify this hypothesis, we compute below numerical simulations of the GJ 876 system, embedded in a disc, using a hydro-code.

  
3.2 Code description and numerical parameters

3.2.1 Numerical scheme

To simulate the entire disc, we used a code developed by Crida et al. (2007), which was derived from FARGO by Masset (2000b,a). The region in which the planets orbit is modelled by a two-dimensional (hereafter 2D) polar grid, which is surrounded by a one-dimensional (hereafter 1D) grid extending over the entire disc; the disc is assumed to be axisymmetrical far from the planets. As pointed out by Crida et al. (2007), this 1D grid enables us to take account of the global disc evolution, which governs the type II migration of the giant planets. In addition, the inner disc evolution is self-consistently computed down to an arbitrarily small inner radius, which could not be reached by a 2D grid for numerical reasons. Thus, the planets can feel the influence of a realistic inner disc. Consequently, this code is well adapted to the problem that we wish to study.

In contrast to previous calculations in which we used a non-inertial frame centred on the star, the usage of an added 1D-grid requires that the frame is inertial and centred on the centre of mass of the system, which is considered to be the star, planet and disc.

The tapering function used here is f, as given in Eq. (5), with p=0.6.

The 2D grid spans the region where the planets orbit, $r\in [0.055;~0.655]$, with a resolution of N_s = 500 sectors in azimuth, and N_r=300 rings in radius; a ring width is thus $\delta r = 0.002$ AU, which applies also for the 1D grid. The outer edge of the 1D grid is fixed arbitrarily to r=10 AU, which is sufficiently distant. The inner edge is discussed below.

3.2.2 Disc parameters

The damping effect of the inner disc is proportional to its mass, in particular to the gas density and the amplitude of the wake close to the planet. The deeper the gap opened by the innermost planet, the smaller is the wake. The wider the gap, the further the disc lies from the planetary orbit. Consequently, a larger gap leads to a smaller damping. The shape of the gap is determined by the gas parameters $\nu$ and h (Crida et al. 2006). These parameters therefore play a crucial role in the damping. The gas viscosity is given by an $\alpha$-prescription, for which $\alpha=10^{-2}$. The chosen aspect ratio is h=0.07.

3.2.3 Radius of the inner edge of the disc

Crida & Morbidelli (2007) showed that the inner disc evolution is strongly dependent on the radius of the inner edge of the disc, and more precisely on the ratio between this radius and the radius of the planetary orbit. This radius is, in general, poorly constrained, and strongly model-dependent. The mass of the inner disc could be uncertain: if this radius is close to 0.13 AU, there would be no inner disc in GJ 876; if it were close to the stellar radius, a massive inner disc could be present.

Fortunately, in the case of GJ 876, a ``hot Neptune'' is present at 0.02 AU. This provides a strong constraint on the location of the inner boundary. Indeed, the migration of this planet stopped for some reason at the place where we observe it now.

This planet may have been caught in a planet trap (Masset et al. 2006), because it is insufficiently massive to open a gap and migrates in the type I regime. The planet trap could be the inner edge of the disc: at this location, the disc density increases rapidly from zero; this probably leads to a positive gradient in the vortensity and a strong corotation torque. In this case, the radius of the inner edge of the disc should be close to 0.02 AU.

	Figure 7: Density profiles at time 0 (dashed line), at the moment where the planets are released (solid line), and at the end of the simulation (after 250 years, dotted line).
Open with DEXTER

A second possibility is that GJ 876 d has migrated inwards in the gas disc, until it is within the empty cavity between the star and the inner edge of the disc. A few reasons can explain why this planet crossed the planet trap, which was created by the disc edge. The trap may have been too weak, due to the absence of a strong enough density gradient. The aspect ratio and viscosity of the disc may have been so low that the planet perturbed the disc profile and destroyed the trap. We note in addition that disc turbulence can create random torques, which can help jump over the trap. Once the planet has crossed the inner edge of the disc, it feels a negative torque from the disc through the outer Lindblad resonances (Goldreich & Tremaine 1980,1979). The planet continues to migrate inward until there is no more gas at the location of its 2:1 resonance (the outermost one). In this case, the inner edge of the disc must be located at the 2:1 resonance with GJ 876 d, i.e. at 0.033 AU. The planet would then remain, while the disc was evaporated by the star, from the inside outwards.

We focus on this second possibility. We claim that the inner disc cannot extend inwards further than 0.033 AU, otherwise it would draw GJ 876 d inward; it should have extended to precisely this radius, or GJ 876 d would not have migrated inward to its present position. The open inner edge of the 1D grid will therefore be located at r=0.033 in our simulation.

3.3 Results

The innermost planet, GJ 876 d, is not modelled in the simulation. The two largest planets are launched on circular orbits at r=0.36 and 0.21 AU respectively. At first, the planets influence each other and they influence the disc, but they are not affected by the disc. The planets shape a gap in the density distribution and the gas disc reaches the equilibrium state for this planetary configuration. This phase lasts for 75 years, which corresponds to approximately 200 orbits for the outer planet. In Fig. 7, we show the initial density profile, the profile after 75 years, and the final density profile. The mass of gas present in the 2D grid, after this first phase, is $2.69\times 10^{-2}$stellar mass ( ${\sim}1.7\times 10^{28}$ kg).

The planets are then released and allowed to move under the influence of the gas. The evolution of their orbital elements is shown in Fig. 8. As expected, the outer planet migrates inward, pushed by the outer disc (see the curve labelled a₂ in Fig. 8). However, the inner planet also migrates inward, although less rapidly. This is because it does not open a clean gap (see Fig. 7) and it occupies the inner edge of the gap opened by GJ 876 b; consequently the inner planet feels a strong negative corotation torque.

The three relevant resonant angles associated with the 2:1 resonance are shown in Fig. 9. After $\sim$110 years, the outer planet catches the inner planet in its 1:2 Mean Motion Resonance ($\theta _1$, $\theta _2$ and $\Delta \omega$ start librating about $0\hbox{$^\circ$ }$ with a small amplitude), and the pair of planets continue to migrate in this configuration. The eccentricities rise, as expected. But after a phase of eccentricity growth, a limit for e₁and e₂ is reached at about 150 years. The planets continue to migrate at the same rate, while their eccentricities remain constant. The value obtained for the eccentricities is close to the value provided by Laughlin et al. (2005). After $\sim$270 years, the planets reach their present semi-major axes, and their eccentricities oscillate within the ranges 0.019<e₂<0.032 and 0.21<e₁<0.25. The amplitudes of the libration of the resonant angles, are $\vert\theta_1\vert _{\rm max}\approx 18\hbox{$^\circ$ }$ in contrast to a value of about $8\hbox{$^\circ$ }$ at 220 years, $\vert\theta_2\vert _{\rm max}\approx 28\hbox{$^\circ$ }$, and $\vert\Delta\omega\vert _{\rm max}\approx 20\hbox{$^\circ$ }$. These amplitudes differ slightly from the values provided by Laughlin et al. (2005), but agreement for the semi-major axes and eccentricities is excellent.

	Figure 8: Semi-major axis (blue, left y-axis) and eccentricity (red, right y-axis) evolution of GJ 876 b (light colour, a1, e1) and GJ 876 c (dark, a2, e2). The horizontal lines correspond to the observed values, as given by Table 1.
Open with DEXTER

	Figure 9: Resonant angles associated with the 2:1 resonance $\theta _1$, $\theta _2$, $\Delta \omega$ as functions of time.
Open with DEXTER

3.4 Discussion

3.4.1 Role of the inner disc

We compute the same simulation without any action of the disc on the inner planet. The semi-major axis evolution of the two planets is almost unaffected. Migration is dominated by the outermost planet; this planet is pushed inward by the outer disc, which in turn causes the inner planet to move inward due to resonance locking. Consequently, the curves of a_i(t) overlap the curves of Fig. 8. In contrast, the behaviour of the eccentricity of the inner planet e₁ changes dramatically. It increases faster and continuously, as expected from N-body simulations. In fact, both eccentricities rise to high values that are incompatible with observations ( $e_1\sim 0.45$, e₂>0.1). The eccentricity evolution is displayed in Fig. 10 for the previous case (labelled ``ref''), and in the case where the action of the disc on the inner planet is switched off (labelled ``no inn.''). This convincingly demonstrates the important role of the inner disc in eccentricity damping for the inner planet, which in turn affects the eccentricity of the outer planet.

Figure 10: Eccentricity evolution of GJ 876 b (light colour, e2) and GJ 876 c (dark, e1) in three cases: (i) standard simulation (same as Fig. 8, curves labelled ``ref''); (ii) same as (i) but the inner planet is not affected by the disc (curves labelled ``no inn.''); (iii) same as (i) but the p parameter of the tapering function f is 0.8 (curves labelled ``p = 0.8''). The horizontal lines correspond to the observed values, as given by Table 1.

Open with DEXTER

3.4.2 Role of the tapering function

We should mention that in the standard simulation, if one takes p=0.8 instead of 0.6 in the tapering function f (see Eq. (5) and Fig. 2), the inner planet reaches a higher value of eccentricity (between 0.275 and 0.3), while the eccentricity of the outer planet saturates between 0.035 and 0.05. The eccentricity evolution in this case is displayed in Fig. 10 (curves labelled ``p = 0.8''). The eccentricities do not reach extremely high values, due to disc damping, but this is not in good agreement with observations (the horizontal lines in Fig. 10), in particular for the inner planet. This demonstrates that the tapering function can have an influence on the final eccentricity of the planets in numerical simulations (as could be expected from Fig. 4), and this should be taken into consideration.

3.4.3 Disc dispersal

To agree with the present state, the disc has to be dispersed when the planets reach their present semi-major axis (at time $t \sim 270$ yr). This is a common problem when modelling the extra-solar planets.

We applied a procedure also used by Morbidelli et al. (2007) and Thommes et al. (2007) to investigate the effect of removing gas from the disc. This procedure models the disappearance of gas from the planetary system, in a smooth way, as to not cause a sudden change in the potential. From time t=250 yr, the gas density decreases within each grid cell, on an exponential timescale of 27.5 years. The results are shown in Fig. 11. As the gas disappears, the planets remain in 2:1 MMR, while their migration speed decreases exponentially. It is indeed well known that when the planet is more massive than the disc, the inertia of the planet is the limiting factor in the type II migration regime. Since the migration speed and eccentricity damping are both proportional to gas density, the K factor is unaffected by this procedure. Consequently, the equilibrium values of the eccentricity are unaffected, and remain close to 0.03 and 0.22 respectively, while the semi-major axes converge to values of 0.21 and 0.13 AU respectively. At the end of the procedure, the gas has almost disappeared and the planetary system is similar to GJ 876.

	Figure 11: Semi-major axis and eccentricity evolution of GJ 876 b and GJ 876 c, with disc clearing from time t=250 yr on. The horizontal lines correspond to the observed values, as given by Table 1.
Open with DEXTER

The disc clearing process is complex and not well constrained, in particular in the vicinity of the star. In general, the gas density first slowly decreases while the disc accretes onto the central star and spreads outward. When the density is low, the photo-evaporation by the central star can play a significant role and erode the disc. The extreme UV photons ionise and heat the upper layer of the disc, gas leaves the potential well of the star, and the density decreases. This works more efficiently at radii greater than approximately 1 AU, because closer to the star, the gravity becomes too strong. Consequently, the region where the two giant planets of GJ 876 are orbiting should not be much affected; in particular, the inner disc should not disappear. The remnant disc inside $\sim$1 AU spreads viscously onto the star and outwards. The migration path of the planets should not be affected. In addition, the viscous evolution dominates the disc evolution until the density becomes very low. Photo-evaporation then happens on a timescale that is far shorter than the disc lifetime. The final phase of gas dispersal occurs when the disc mass density is too low to have a significant influence on the planets over such a short timescale.

We believe that the planetary configuration should not be perturbed significantly during this phase, and that the model of exponential damping of the density presented above, provides reliable results. In any case, observations show that these planets migrated toward their host star, until the gas density became too low; the planets stopped at 0.13 and 0.21 AU, when the disc disappeared. Our simulation demonstrates that when the planets stop migrating, they automatically have the correct eccentricities; our simulation is therefore consistent with the observed configuration. To our knowledge, this is the first time that an extra-solar system is reproduced in a fully hydro simulation taking into account all the protoplanetary disc and allowing for a significant migration of the planets with correct eccentricities.

  
4 Modelling an inner disc by $\textit {N}$-body calculations

Full hydrodynamical calculations (like the ones done in the previous section) require typically a large amount of computer time to study the effect of an inner disc on the evolution of resonant planetary systems. Fortunately, the effects of the outer and inner discs can be modelled approximately by gravitational N-body simulations using appropriately parameterised non-conservative forces. These forces can be derived using the migration rate $\dot a/a$ and the eccentricity damping rate $\dot e/e$ of a planet, embedded into the protoplanetary disc (see Beaugé et al. 2006; Lee & Peale 2002, for two different approaches). In place of the migration and damping rates we can also use the corresponding e-folding times defined as $\tau_a = -(\dot a/a)^{-1}$ and $\tau_e=-(\dot e/e)^{-1}$.

When studying the formation of a resonant system consisting of an inner and an outer giant planet, the outer planet is usually forced to migrate inward. When the ratio of the semi-major axes of the planets approaches a critical limit, a resonant capture between them can take place (for the conditions of a resonant capture into the 2:1 MMR see Kley & Sándor 2007). After the resonant capture, the two planets migrate inward as the outer planet continues to be affected by the negative tidal torques of the disc.

Following the hydrodynamical evolution of GJ 876 presented in the previous section, we provide further results for a three-body problem with dissipative forces. During migration of the giant planets, the presence of an inner disc is found to be consistent with the observed state of the resonant systems, in which two giant planets are engaged in a 2:1 MMR.

Since the eccentricity of the inner planet is excited by the resonant interaction, its orbit becomes more and more elongated penetrating into the inner disc. This represents a damping mechanism that acts against the eccentricity excitation, and may set a quasi-equilibrium between the processes keeping the eccentricity of the inner planet at a constant value (within certain limits) during the whole migration process. The effect of an inner disc can be investigated by using a repelling non-conservative force, acting on the inner planet, parameterised by a positive migration rate $\dot a/a$ and a negative $\dot e/e$. A ratio K of the above parameters can be defined, which is similar to the case of the inwardly-migrating outer planet. According to the definition of the e-folding times, $\tau_a$will have a negative sign.

We note, however, that when both the outward and inward migration of the outer and inner planets are considered simultaneously, the final state of the system cannot be characterised uniquely using only the K ratios to describe the inward and outward migration. The final values of the semi-major axes and the eccentricities depend directly on the migration rates (or on the corresponding e-folding times) $\dot a_1/a_1$, $\dot a_2/a_2$ and the eccentricity damping rates $\dot e_1/e_1$, $\dot e_2/e_2$ (where indices ``1'' and ``2'' represent the inner and outer planet, respectively), and not only on their ratios K₁ and K₂. At the end of the section, we characterise the system final state using migration parameters.

We repeat the three-body calculations for GJ 876 by Lee & Peale (2002) adding the effects of an inner disc; we then review the simulations for HD 73 526 by Sándor et al. (2007), and present new results for the modelling of the formation of both HD 82 943 and HD 128 311, with an inner disc.

4.1 GJ 876 and HD 73 526

Figure 12: Behaviour of the semi-major axes and the eccentricities of the resonant giant planets in the system GJ 876. The horizontal lines correspond to the observed values of the eccentricities. Top: only the outer disc is taken into account, and therefore only the orbital evolution of the outer planet is affected, with the e-folding times $\tau _{a_2}=2\times 10^4$ years, $\tau _{e_2}=2\times 10^2$ years, so K2=100. Bottom: same as above in the presence of outer and inner discs. In this case, $\tau _{a_1}=-2\times 10^4$, $\tau _{e_1}=2.5\times 10^3$, $\tau _{a_2}=2\times 10^4$, and $\tau _{e_2}=2.5\times 10^3$ years (giving K1=K2=8).

Open with DEXTER

GJ 876:

We repeated the simulations of Lee & Peale (2002) of the formation of GJ 876, using three-body calculations with dissipative forces. As for the study of Lee & Peale (2002), the outer giant alone is affected by an outer disc, and we set $K_2=\tau_{a_2}/\tau_{e_2}=100$. The evolution of the semi-major axes and eccentricities is shown in the top panel of Fig. 12. We assumed that $\tau _{a_2}=2\times 10^4$ and $\tau _{e_2}=2\times 10^2$ years, which implies that K₂=100. We note again that this ratio between the e-folding times is very high and may be physically unrealistic for massive planets.

To study the damping effect of the inner disc on the inner planet, while the giant planets revolve in a common gap, we performed further three-body simulations, with also $\tau _{a_1}=-2\times 10^4$ and $\tau _{e_1}=2.5\times 10^3$ years; we point out that the minus sign of $\tau _{a_1}$ stands for the outward migration. To model the damping effect of the outer disc on the outer planet we used the e-folding times $\tau _{a_2}=2\times 10^4$ and $\tau _{e_2}=2.5\times 10^3$ years. We note that the above migration parameters correspond to K₁=K₂=8. The result of our calculations is shown in the bottom of Fig. 12. These figures are very similar to those obtained by Lee & Peale (2002), even though we used different migration parameters.

HD 73 526:

In the case of HD 73 526, as already mentioned, the effect of an inner disc was studied in detail by Sándor et al. (2007). If only the outer planet is damped, it was shown that K₂ = 15-20 is required, which is a too high value according to the hydrodynamical simulations. Moreover, if only the outer planet is damped, the eccentricity of the inner planet continues to increase, which can result in the limit imposed by observations being exceeded ( $e_{\rm inn} \approx 0.3$). If, in addition, the inner planet is damped by an inner disc, the eccentricities will remain constant, and the outcome of the migration scenario can more easily reproduce the observations: the formation of HD 73 526 could be modelled successfully by using the e-folding times $\tau_{a_1}=-5\times 10^4$, $\tau_{a_2}=10^4$ years and $\tau_{e_1}=5\times 10^3$, $\tau_{e_2}=10^3$ (corresponding to K₁=K₂=10). The behaviour of the eccentricities in this case is shown in the top panel of Fig. 8 in the paper of Sándor et al. (2007).

We conclude that including the effects of an inner disc provides more physically-realistic model predictions of the formation of GJ 876 and HD 73 526, which are in good agreement with radial-velocity observations. In what follows, we show that inner discs may have been present in the systems HD 82 943 and HD 128 311, although a strong damping of the eccentricity of their inner giant planets may not be required.

4.2 HD 82 943 and HD 128 311

Figure 13: Behaviour of the semi-major axes and the eccentricities of the giant planets in the resonant system HD 82 943. The horizontal lines correspond to the observed values of the eccentricities. Top: only the motion of the outer planet is affected, with the e-folding times $\tau _{a_2}=2\times 10^4$ and $\tau _{e_2}=2.5\times 10^3$ years (K2=8). Bottom: an inner disc is also assumed, with $\tau _{a_1}=-2\times 10^5$, $\tau _{e_1}=5\times 10^4$ years e-folding times for the inner planet (K1=4), which allows the use of larger $\tau _{e_2}=2.5\times 10^3$ years and thus lower ratio K2 =4 for the outer planet.

Open with DEXTER

Figure 14: Behaviour of the semi-major axes and the eccentricities of the resonant giant planets in the resonant system HD 128 311. The horizontal lines correspond to the observed values of the eccentricities. Top: only the motion of the outer planet is damped by an outer disc, with $K_2=\tau _{a_2}/\tau _{e_2}=7$. Bottom: an inner disc is also assumed. The e-folding times for the inner planet are $\tau _{a_1}=2\times 10^4$, $\tau _{e_1}=10^4$ years (K1=2), which allows to use much lower ratio K2 =2 for the outer planet.

Open with DEXTER

HD 82 943:

Recently Lee et al. (2006) presented four sets of orbital solutions based on a best-fit double-Keplerian model for HD 82 943. Using their orbital elements as initial conditions, direct numerical integrations showed that three solutions exhibited ordered behaviour, while one solution (Fit I in the cited paper) was destabilised after a few thousand years. In one of the three dynamically-stable orbital solutions (Fit II), the variations if the eccentricities are negligible, and the resonance variables oscillate with small amplitudes, which indicates clearly that the system is deep inside a 2:1 MMR. Since this behaviour can be reached during an inward convergent migration of the giant planets, we investigate the formation of this system using the scenario of planetary migration.

We assume that only an outer disc is present and that the outer planet is affected by its damping. We derive the dynamical state calculated based on Fit II (e₁=0.422, e₂=0.14, a₁=0.74 AU, a₂=1.18 AU), by assuming the ratio K₂ = 8, using $\tau _{a_2}=2\times 10^4$ and $\tau _{e_2}=2.5\times 10^3$ years. This ratio K₂ does not appear to be too high, with values typically between 1 and 10. However, as in the case of HD 73 526, the eccentricities increase slightly during the entire migration process. This is not a problem if the migration of the planets is terminated gradually, about the observed values of their semi-major axes; this may lead, however, to the constraints placed by the observations being exceeded: see the top panel of Fig. 13. On the other hand, by assuming the presence of an inner disc, the eccentricities gradually reach their constant values, which do not appear to change further during the migration process. The latter behaviour is shown in the bottom panel of Fig. 13. During this simulation, we assumed for the inner planet that $\tau _{a_1}=-2\times 10^5$ and $\tau _{e_1}=5\times 10^4$ years, and for the outer planet that $\tau _{a_2}=2\times 10^4$ and $\tau _{e_2}=5\times 10^3$ years (giving K₁=K₂=4). This case appears to be a more appropriate formation scenario for HD 82 943 than the previous one, because the eccentricities appear to reach constant values.

We conclude that the present dynamical behaviour of the resonant system HD 82 943 (based on Fit II of Lee et al. 2006) can be explained in both ways: by assuming an outer disc alone or by assuming an inner and an outer disc. The latter case, however, appears to be more reasonable because it maintains a constant value of eccentricity, when the migration occurs over longer times.

HD 128 311:

Finally, we present our results for HD 128 311. The most recent orbital solution for this system is presented by Vogt et al. (2005), while possible formation scenarios are outlined by Sándor & Kley (2006). In the latter study, the formation of the resonant system is modelled by an inward convergent migration of the giant planets, followed by a sudden perturbation; only the effect of an outer disc was considered. Prior to the sudden perturbation, the outer giant planet migrates inward and captures the inner planet into a 2:1 resonance. It is important to find the final eccentricity values, produced by the migration process, because the perturbative events modify only the oscillations of the eccentricity. These values are $e_1\approx 0.45$ and $e_2 \approx 0.15$, obtained after a migration characterised by $K_2=\tau_{a_2}/\tau_{e_2}=5$. In the present study, we require a larger ratio, K₂=7, because in contrast to Sándor & Kley (2006), we do not stop migration when the planets reach their actual positions. The evolution of the semi-major axes and eccentricities are presented in the top panel of Fig. 14.

Assuming an inner disc, the same final state following the migration of HD 128 311 can be obtained. In our numerical experiments, we used $\tau _{a_1}=-2\times 10^4$, $\tau _{e_1}=10^4$ years for the inner planet, and $\tau _{a_2}=2\times 10^4$, $\tau_{e_2}=10^4$ years for the outer planet, which are equivalent to K₁=K₂=2. The evolution of the system is shown in the bottom panel of Fig. 14. Comparing the top and the bottom panels of Fig. 14, we conclude that the formation of HD 128 311 can be modelled by assuming the presence of an inner disc; in this case, however, it is not necessarily required.

   
4.3 Final state of a migration with an inner disc: numerical experiments

In this part we show the results of additional numerical simulations in which we studied how the different parameters of migration influence the final state of the resonant system. It was quite evident from the beginning that the characterisation based only on the ratios K₁ and K₂ would not yield unique results in the case when the inner planet is also affected by an inner disc.

Our calculations are based on the observed orbital and physical parameters of GJ 876, as given by Lee & Peale (2002). For the inward migration of the outer planet we fixed $\tau _{a_2}=4\times 10^4$, $\tau _{e_2}=5\times 10^3$ years (so, K₂=8), and we changed $\tau _{a_1}$ and $\tau _{e_1}$ in such a way to obtain the observed state of the system.

As starting values to our numerical simulations we used $\tau_{a_1} = -3\times 10^4$ years and found that $\tau_{e_1} = 3.7\times 10^3$years gives a correct result ( $K_1\approx 8$ here). Then we increased gradually the absolute value of $\tau _{a_1}$ corresponding to a weaker damping on a₁. To obtain the observed behaviour of the system, for larger $\tau _{a_1}$, we needed a larger $\tau _{e_1}$ meaning a weaker damping on e₁ too. To explore the mutual dependence of $\tau _{a_1}$and $\tau _{e_1}$, we performed a series of numerical experiments. Our results are shown in Fig. 15: we display on the x-axis the absolute value of $\tau _{a_1}$, while on the y-axis we display the eccentricity damping time $\tau _{e_1}$, on a logarithmic scale. The crosses indicate the values of $\tau _{a_1}$ and $\tau _{e_1}$ required to derive similar final values of the eccentricities e₁ and e₂ (0.26 and 0.035 respectively). The value of $\tau _{e_1}$ increases rapidly for small $\vert\tau_{a_1}\vert$; it is not, however, proportional to $\tau _{a_1}$: compare the data with the straight line in Fig. 15; it does level off and tends to a limiting value of $\tau_{e_1 \rm max}$, which in this case is slightly higher than $9\times 10^3$ years.

	Figure 15: Pairs of migration parameters $\tau _{a_1}$ and $\tau _{e_1}$ which result in the observed state of the system GJ 876 for fixed $\tau _{a_2}=4\times 10^4$ years and $\tau _{e_2}=5\times 10^3$ years. Crosses: empirically found with N-body simulations. Dashed curve: analytical expression Eq. (16) (see Sect. 4.4). Solid line: constant ratio K1=8, for comparison.
Open with DEXTER

We conclude that for a fixed pair of $\tau_{a_2}$ and $\tau_{e_2}$, there is no unique K₁ that determines the final state of the system. In contrast, there exists a $\tau_{e_1 \rm max}$, which determines the final state of the system if the inner disc does not affect the semi-major axis of the inner planet. In reality, however, the inner disc, which is rotating more rapidly, can transmit angular momentum to the inner planet, and energy as well, which increases the semi-major axis of the inner planet. In this case we require a smaller value of $\tau _{e_1}$, or equivalently, more effective damping on e₁.

We compare the two migration scenarios: (i) only an outer disc is present and there is no inner disc between the inner planet and the star; and (ii) both an outer and an inner disc are present and the planets orbit in a gap between the discs. We summarise our results. The final state of the system in case (i) can be described solely by the ratio $K_2=\tau_{a_2}/\tau_{e_2}$. In case (ii), the ratios K₁ and K₂ are insufficient to describe the final state of the system, which depends on additional parameters. In the most general case, when the inner disc forces the inner planet to migrate outward, we require three parameters, which can be K₁, K₂, and the ratio $\tau_{a_1}/\tau_{a_2}$. In the next subsection, we analyse theoretically this behaviour.

  
4.4 Final state of a migration with an inner disc: analytics

For the semi-major axis and eccentricity of a planet to decrease exponentially with damping timescales $\tau_a=-\dot{a}/a$ and $\tau_e=-\dot{e}/e$, its energy E and angular momentum H must vary with the following rates (through Eqs. (2) and (3) ):

  

$$\dot{E} E/\tau_a$$ = (6)

$$\dot{H} \frac{H}{2} \left(\frac{2e^2}{1-e^2}\frac{1}{\tau_e}-\frac{1}{\tau_a}\right)\cdot$$ = (7)

When two planets are considered, the variation of the energy of the system, i.e. pair of planets, $\dot{E}$ is the sum of the energy variations applied to both planets, such that

$$\dot{E} = \frac{E_1}{\tau_{a_1}}+\frac{E_2}{\tau_{a_2}}\cdot$$ (8)

If the two planets are in resonance, the ratio between their semi-major axes is kept constant, and consequently also the ratio between their energies. We define:

$$\varepsilon = E_2/E_1 = M_2a_1/M_1a_2.$$ (9)

Then,

$$E = E_1+E_2 = (1+\varepsilon)E_1,$$ (10)

$$\dot{E} = (1+\varepsilon) \dot{E_1},$$ (11)

we note that due to some energy exchange between the planets through the resonance, the variation of the energy of the first planet is not $E_1/\tau_{a_1}$, but $\dot{E_1}=\dot{E}/(1+\varepsilon)$, with $\dot{E}$ given by Eq. (8).

A similar situation holds for the angular momentum. For the angular momentum H of the system, i.e. pair of planets, (H=H₁+H₂), the total variation rate is given by (from Eq. (7)):

$$\dot{H} = \sum_{i=1,2}\frac{H_i}{2} \left(\frac{2{e_i}^2}{1-{e_i}^2}\frac{1}{\tau_{e_i}}-\frac{1}{\tau_{a_i}}\right)\cdot$$ (12)

When the planets are in resonance and their eccentricities are constant, i.e. in their final state, the ratio of their angular momentum is defined to be:

$$\eta = H_2/H_1 = \frac{M_2}{M_1}\sqrt{\frac{a_2(1-{e_2}^2)}{a_1(1-{e_1}^2)}},$$ (13)

which is also constant. Then, $H=H_1+H_2=(1+\eta)H_1$, and

$$\dot{H} = (1+\eta)\dot{H_1}.$$ (14)

If the eccentricities are constant ( $\dot{e_i}=0$), then, using Eqs. (3), (10), (11), and (8), we derive:

$$\dot{H_i} = -\frac{H_i}{2}\frac{\dot{E_i}}{E_i} = -\frac{H_i}... .../\tau_{a_1}+\varepsilon/\tau_{a_2}}{1+\varepsilon}\right)\cdot$$ (15)

Using Eqs. (12) and (15), Eq. (14) provides a relation between $\tau _{a_1}$, $\tau_{a_2}$, $\tau _{e_1}$, and $\tau_{e_2}$. In the problem that we study, the masses of the planets are known, as well as their eccentricities and the ratio between their semi-major axes. The question is: for a given effect on the outer planet, what should the effect on the inner planet be to be consistent with the observed eccentricities ? Below, we solve this equation in the unknown $\tau _{e_1}$ with the free parameter $\tau _{a_1}$, while $\tau_{a_2}$, $\tau_{e_2}$, $\varepsilon$ and $\eta$, are given constants ( $\varepsilon$ and $\eta$ being defined by Eqs. (9) and (13) respectively). Using Eqs. (12) and (15), Eq. (14) reads:

$$\sum_{i=1,2}\frac{H_i}{2} \left(\frac{2{e_i}^2}{1\!-\!{e_i}^2... ...frac{1/\tau_{a_1}+\varepsilon/\tau_{a_2}}{1+\varepsilon}\right)$$

$$\frac{2{e_1}^2}{1\!-\!{e_1}^2}\frac{1}{\tau_{e_1}} - \frac{1}... ...left(\frac{\varepsilon}{\tau_{a_2}}+\frac{1}{\tau_{a_1}}\right)$$

$$\frac{2{e_1}^2}{1\!-\!{e_1}^2}\frac{1}{\tau_{e_1}} = \frac{1}... ...a_1}}\left(1-\frac{1\!+\!\eta}{1\!+\!\varepsilon}\right)\ \cdot$$

Eventually, we obtain:

$$\frac{1}{\tau_{e_1}} = \underbrace{\frac{1\!-\!{e_1}^2}{2{e_1... ...{e_1}^2}\frac{\varepsilon-\eta}{1+\varepsilon}}_{K_{1,0}}\cdot$$ (16)

As one can see, the damping rate that should be applied to the inner planet ( $1/\tau_{e_1}$) is the sum of two terms.

The first one ( $1/\tau_{e_1 \rm max}$) is required to balance the action on the outer planet; when no force is applied to the inner planet, the energy loss rate of the outer planet is not the expected value of $E_2/\tau_{a_2}$, but $(E_2/\tau_{a_2}) \frac{\varepsilon}{1+\varepsilon}$, because the inner planet has to lose energy to preserve the resonant motion; the angular momentum loss rate is therefore overestimated by the expression Eq. (7), and both values of eccentricity rise. A damping of the eccentricity needs to be applied, in addition, to the inner planet as well.

The second term ( $K_{1,0}/\tau_{a_1}$) is proportional to ${\tau_{a_1}}^{\!-1}$; the coefficient K_1,0 is negative if e₁²>(1-(a₂/a₁)³)+(a₂/a₁)³e₂², which is always true in the case of a 2:1 MMR for e₂<0.866.

In the case studied in Sect. 4.3, one has (Lee & Peale 2002): $M_2=1.808~ M_{\rm Jup}=5.65\times 10^{-3}M_*$ , $M_1=0.5696~ M_{\rm Jup}=1.78\times 10^{-3}M_*$ , $e_1\approx 0.265$ , $e_2\approx 0.035$ , $\tau _{a_2}=4\times 10^4$ , $\tau _{e_2}=5\times 10^3$ years and the two planets are in 2:1 MMR, with a₂/a₁=1.602. Thus, $\eta=4.1641$ and $\varepsilon=1.9812$. This gives $\tau_{e_1 \rm max}=9~289$ years, and K_1,0=-4.8471. The dashed curve in Fig. 15 shows $\tau _{e_1}$ as a function of $\tau _{a_1}$ from Eq. (16). The fit is excellent.

This shows that Eq. (16) provides efficiently the parameters required to reproduce a system using N-body simulations with dissipative forces, or that one should have in the protoplanetary disc. In particular, the expression for $\tau_{e_1 \rm max}$ shows that if e₂ is relatively small, a large value of K₂ of the order of e₂^-2 is required for $\tau_{e_1 \rm max}\to \infty$. If the inner disc causes the inner planet to migrate outward with $\tau_{a_1}<0$, an additional damping on its eccentricity is required, which can be expressed as $\tau_{e_1 \rm add}=\tau_{a_1}/K_{1,0}$, with K_1,0<0, depending on all parameters of the system. If the inner disc attracts the planet sufficiently (which can happen as shown in Sect. 2), no eccentricity damping on the inner planet is required: in the above case, this occurs for $\tau_{a_1} = +45~025$ years.

Using this equation, reasonable parameters, i.e. not too large K₁and K₂, could be derived to explain how all known resonant, exoplanetary systems, were formed in the disc. Our numerical simulations have shown that this can be achieved for at least 4 such systems.

  
5 Conclusion

We address the problem of moderate eccentricity of extra-solar planets in resonance. The resonant configuration requires a convergent migration of the planets, but continued migration in resonance leads to unlimited eccentricity growth if no eccentricity damping mechanism is at work.

In Sect. 2, we showed that an inner disc has a non-negligible influence on a giant planet in an eccentric orbit, and modifies its orbital parameters. The strength of this effect varies with planetary eccentricity. For small eccentricities, $e \leq 0.05$, we find a small but positive $\dot{e}$, while for larger e, $\dot{e}$ becomes more and more negative. The induced change in semi-major axis remains positive, as expected for the Lindblad torques induced by an inner disc on a massive planet. The change in semi-major axis becomes negative only for larger eccentricities $e\geq 0.35$, leading to inward migration. The use of a constant K-factor between the e-folding timescales of the semi-major axis and the eccentricity is clearly an oversimplification.

The measured influence of the disc on the planet depends on the adopted tapering function, applied to exclude at least parts of the gas within the Hill sphere of the planet, and possibly other numerical effects, such as the smoothing length, the resolution, or boundary conditions. Determination of the precise, absolute magnitude of the influence is difficult.

An eccentricity damping should be applied to the inner planet to obtain realistic results, in particular if its eccentricity is larger than 0.1. Using a hybrid 2D-1D hydro-code that allows the simulation of the entire disc from its physical inner boundary to an arbitrarily large radius, we compute the long-term evolution of the GJ 876 system, for which the orbital elements are well known, and the radius of the inner edge of the disc can be estimated due to the presence of a third planet close to the star. We find that the inner disc does not disappear after the planets open a gap, and that it helps to dampen the eccentricity of the inner planet. With realistic disc parameters (viscosity and aspect ratio), we can to reproduce the observations.

Since hydro-simulations are time-consuming, we perform customised N-body simulations with non-conservative forces, which are added to simulate the effect of the outer and inner disc. We find that applying a significant damping to both planets, to reflect the influence of both the inner and outer discs, enables a few other exoplanetary systems to be reproduced well. In addition, N-body simulations show that when two planets are considered, the ratio between the eccentricity and semi-major axis damping applied to each of them (K₁ and K₂), are insufficient to determine the final state of the system. For a given damping in a and e applied on the outer planet, one can express analytically the eccentricity damping that should be applied to the inner planet to match the observed orbital configuration, as a function of its semi-major axis damping.

Given the satisfying fits of a few systems that we obtain, we claim that the problem of the low eccentricities of the resonant exoplanetary systems simply stems from the fact that the inner disc has not previously been taken properly into account, which is not reasonable. A pair of planets orbiting in a protoplanetary disc may well orbit in a common gap and enter a Mean Motion Resonance, but it does not necessarily open a complete cavity from the star to the outermost planet, so that the inner disc should influence the planets dynamics. Crida & Morbidelli (2007) suggested that the opening of such a cavity by a giant planet requires specific conditions; according to our hydro and N-body simulations, the low observed eccentricities of the exoplanets in resonance support the idea that the inner disc does in general not disappear.

Acknowledgements

A. Crida acknowledges the support through the German Research Foundation (DFG) grant KL 650/7. Zs. Sándor thanks the supports of the Hungarian Scientific Research Fund (OTKA) under the grant PD48424 and of the DFG under the grant 436 UNG 17/1/07.

Very fruitful discussions with F. Masset and C. Dullemond are gratefully acknowledged.

References

• Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006, MNRAS, 365, 1160 [NASA ADS] [CrossRef]
• Bryden, G., Rózyczka, M., Lin, D. N. C., & Bodenheimer, P. 2000, ApJ, 540, 1091 [NASA ADS] [CrossRef]
• Cresswell, P., Dirksen, G., Kley, W., & Nelson, R. P. 2007, A&A, 473, 329 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Crida, A., & Morbidelli, A. 2007, MNRAS, 377, 1324 [NASA ADS] [CrossRef] (In the text)
• Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 [NASA ADS] [CrossRef] (In the text)
• Crida, A., Morbidelli, A., & Masset, F. 2007, A&A, 461, 1173 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• D'Angelo, G., Lubow, S. H., & Bate, M. R. 2006, ApJ, 652, 1698 [NASA ADS] [CrossRef] (In the text)
• Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [NASA ADS] [CrossRef]
• Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [NASA ADS] [CrossRef]
• Kley, W. 2000, MNRAS, 313, L47 [NASA ADS] [CrossRef]
• Kley, W., & Dirksen, G. 2006, A&A, 447, 369 [NASA ADS] [CrossRef] [EDP Sciences]
• Kley, W., & Sándor, Z. 2007, in Extrasolar Planets: Formation, Detection and Dynamics, ed. R. Dvorak, 945 (In the text)
• Kley, W., Peitz, J., & Bryden, G. 2004, A&A, 414, 735 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Kley, W., Lee, M. H., Murray, N., & Peale, S. J. 2005, A&A, 437, 727 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Laughlin, G., Butler, R. P., Fischer, D. A., et al. 2005, ApJ, 622, 1182 [NASA ADS] [CrossRef] (In the text)
• Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [NASA ADS] [CrossRef] (In the text)
• Lee, M. H., Butler, R. P., Fischer, D. A., Marcy, G. W., & Vogt, S. S. 2006, ApJ, 641, 1178 [NASA ADS] [CrossRef] (In the text)
• Lin, D. N. C., & Papaloizou, J. 1979, MNRAS, 186, 799 [NASA ADS]
• Lubow, S. H. 1991, ApJ, 381, 259 [NASA ADS] [CrossRef] (In the text)
• Masset, F. 2000a, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences]
• Masset, F. 2000b, in Disks, Planetesimals, and Planets, ed. G. Garzón, C. Eiroa, D. de Winter, & T. J. Mahoney, ASP Conf. Ser., 219, 75
• Masset, F., & Snellgrove, M. 2001, MNRAS, 320, L55 [NASA ADS] [CrossRef] (In the text)
• Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478 [NASA ADS] [CrossRef] (In the text)
• Morbidelli, A., & Crida, A. 2007, Icarus, 191, 158 [NASA ADS] [CrossRef] (In the text)
• Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F., & Gomes, R. 2007, ArXiv e-prints, 706 (In the text)
• Papaloizou, J. C. B. 2003, Celestial Mechanics and Dynamical Astronomy, 87, 53 [NASA ADS] [CrossRef]
• Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 366, 263 [NASA ADS] [CrossRef] [EDP Sciences]
• Sándor, Z., & Kley, W. 2006, A&A, 451, L31 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Sándor, Z., Kley, W., & Klagyivik, P. 2007, A&A, 472, 981 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] (In the text)
• Snellgrove, M. D., Papaloizou, J. C. B., & Nelson, R. P. 2001, A&A, 374, 1092 [NASA ADS] [CrossRef] [EDP Sciences]
• Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388 [NASA ADS] [CrossRef] (In the text)
• Thommes, E. W., Bryden, G., Wu, Y., & Rasio, F. A. 2007, ArXiv e-prints, 706 (In the text)
• Vogt, S. S., Butler, R. P., Marcy, G. W., et al. 2005, ApJ, 632, 638 [NASA ADS] [CrossRef] (In the text)