On the migration of a system of protoplanets

Abstract

1 Introduction

It is generally assumed that planetary systems form in a differentially rotating gaseous disc. In the late stages of their formation, the protoplanets are still embedded in the protostellar disc and their orbital evolution is coupled to that of the disc. Gravitational interaction between the planets and the gaseous disc has basically two effects.

• (i)

  The torques by the planets acting on the disc tend to push away material from the orbital radius of the planet, and for sufficiently massive planets a gap is formed in the disc (Papaloizou & Lin 1984; Lin & Papaloizou 1993; Takeuchi, Miyama & Lin 1996). The dynamical process of gap formation has been studied through time-dependent hydrodynamical simulations for planets on circular orbits by Bryden et al. (1999) and Kley (1999, hereafter Paper I), and more recently by Lubow, Seibert & Artymowicz (1999).

  The results indicate that, even after the formation of a gap, the planet may still accrete material from the inner and outer disc and reach about 10 Jupiter masses (M_Jup). For very low disc viscosity and larger planetary masses the mass accumulation finally terminates (Bryden et al. 1999; Lubow et al. 1999).

• (ii)

  Additionally, the gravitational force exerted by the disc alters the orbital parameter (semimajor axis, eccentricity) of the the planet. Here, these forces typically induce some inward migration of the planet (Goldreich & Tremaine 1980) which is coupled to the viscous evolution of the disc (Lin & Papaloizou 1986). A two-dimensional computation of the migration of an embedded planet in an inviscid disc has been presented by Nelson & Benz (1999). Hence the present location of the observed planets (solar and extrasolar) may not be identical to the position at which they formed.

In particular, the migration scenario applies to some of the extrasolar planets (for a summary of their properties see Marcy, Cochran & Mayor 1999): the 51 Peg-type planets. They all have masses of the order of M_Jup, and orbit their central stars very closely, having orbital periods of only a few days. As massive planets, according to standard theory, have formed at distances of a few au from their stars, these planets must have migrated to their present position. The inward migration was eventually halted by tidal interaction with the star or through interaction with the stellar magnetosphere (Lin, Bodenheimer & Richardson 1996). The only extrasolar planetary system around a main-sequence star known so far (ν And) consists of one planet at 0.059 au on a nearly circular orbit and two planets at 0.83 and 2.5 au having larger eccentricities (0.18 and 0.41) (Butler et al. 1999).

In the case of the Solar system, the question arises of what prevented any further inward migration of Jupiter. As the net tidal torque on the planet is a delicate balance between the torque of the material located outside the planet and that of the material inside (e.g. Ward 1997), any perturbation in the density distribution may change this balance. The only evolutionary two-dimensional (r−ϕ) computations of migrating embedded planets in discs have been performed by Nelson et al. (2000), who studied the influence of various physical conditions on the rate of migration of the planet. In this Letter we consider the effect that an additional planet in the disc has on the migration rate.

We present the results of numerical calculations of a thin, non-self-gravitating, viscous disc with two embedded protoplanets. Initially the planets, each of mass 1 M_Jup, are on circular orbits at a=1a_Jup and 2a_Jup, respectively, where a_Jup is the semimajor axis of Jupiter's orbit around the sun. We take into account the back-reaction of the gravitational force of the disc on the orbital elements of the planets and star (see also Nelson et al. 2000). The models are run for about 3000 orbital periods of the inner planet, corresponding to 32 000 yr. In Section 2 a description of the model is given; the results are presented in Section 3 and our conclusions are given in Section 4.

2 The model

We consider a non-self-gravitating, thin accretion disc model for the protoplanetary disc located in the z=0 plane and rotating around the z-axis. Its evolution is described by the two-dimensional (r−ϕ) Navier-Stokes equations, which are given in detail in Paper I. The motion of the disc takes place in the gravitational field of the central star with mass M_s and the two embedded protoplanets with masses m₁ and m₂. In contrast to Paper I we use here a non-rotating frame, as both planets have to be moved through the grid. The gravitational potential is then given by

(1)

where G is the gravitational constant and r_s, r₁ and r₂ are the radius vectors to the star and the two planets, respectively. The quantities s₁ and s₂ are smoothing lengths which are 1/5 of the corresponding sizes of the Roche lobes. This smoothing of the potential allows the motion of the planets through the computational grid.

The motion of the star and the planets is determined first by their mutual gravitational interaction and secondly by the gravitational forces exerted on them by the disc. The acceleration of the star a_s, for example, is given by

(2)

where the integration is over the whole disc surface, and Σ denotes the surface density of the disc. The expressions for the two planets follow similarly. We work here in an accelerated coordinate frame where the origin is located in the centre of the (moving) star. Thus, in addition to the gravitational potential (1), the disc and planets feel the additional acceleration -a_s.

The mass accreted from the disc by the planets (see below) has some net angular momentum which in principle also changes the orbital parameter of the planets. However, this contribution is typically about an order of magnitude smaller than the tidal torque (Lin et al. 2000) and is neglected here.

As the details of the origin and magnitude of the viscosity in discs are still uncertain, we assume a Reynolds stress formulation (Paper I) with a constant kinematic viscosity. The temperature distribution of the disc is fixed throughout the computation and is given by the assumption of a constant ratio of the vertical thickness H and the radius r. Hence the fixed temperature profile is given by T(r)∝r⁻¹. We assume Hr=0.05, which is equivalent to a fixed Mach number of 20.

For numerical convenience we introduce dimensionless units, in which the unit of length, R₀, is given by the initial distance of the first planet to the star R₀r₁(t=0)=1a_Jup. The unit of time is obtained from the (initial) orbital angular frequency Ω₁ of the first planet, i.e. the orbital period of planet 1 is given by

(3)

The evolutionary time of the results of the calculations as given below will usually be stated in units of P₁. The unit of velocity is then given by v₀R₀t₀. The unit of the kinematic viscosity coefficient is given by ν₀R₀v₀. Here we take a typical dimensionless value of 10⁻⁵ corresponding to an effective α of 4×10⁻³.

2.1 The numerical method in brief

The normalized equations of motion are solved using an Eulerian finite difference scheme, where the computational domain r_min,r_max]×[ϕ_min,ϕ_max] is subdivided into N_r×N grid cells. For the typical runs we use N_r=128 and N=128, where the azimuthal spacing is equidistant, and the radial points have a closer spacing near the inner radius. The numerical method is based on a spatially second-order-accurate upwind scheme (monotonic transport), which uses a formally first-order time-stepping procedure. The methodology of the finite difference method for disc calculations is outlined in Kley (1989) and Paper I.

The N-body module of the program uses a fourth-order Runge-Kutta method for the integration of the equations of motion. It has been tested for long-term integrations studying the onset of instability in the three-body problem consisting of two closely spaced planets orbiting a star, as described by Gladman (1993). For the initial parameter used here, the error in the total energy after 1.2×10⁵ orbits (integration over 10⁶ yr) is less than 2×10⁻⁹.

2.2 Boundary and initial conditions

To cover fully the range of influence of the planet on the disc, we typically choose for the radial extent (in dimensionless units, where planet 1 is located initially at r=1) r_min=0.25, r_max=4.0. The azimuthal range covers a complete ring ϕ_min=0.0, ϕ_max=2π using periodic boundary conditions. To test the accuracy of the migration, a comparison calculation with r_max=8.0 and higher resolution N_r=256, N=256 has also been performed.

The outer radial boundary is closed to any mass flow v(r_max)=0, while at the inner boundary mass outflow is allowed, simulating accretion on to the central star. At the inner and outer boundaries the angular velocity is set to the value of the unperturbed Keplerian disc. Initially, the matter in the domain is distributed axially symmetrically with a radial surface density profile Σ∝r^−1/2.

Two planets, each with an initial mass of 1 M_Jup, are located at r₁=1.0, ϕ₁=π and r₂=2.0, ϕ₂=0. Thus they are not only spaced in radius, but are also positioned (in ϕ) in opposition to each other to minimize the initial disturbance. The radial velocity v is set to zero, and the angular velocity is set to the Keplerian value of the unperturbed disc.

Around the planets we then introduce an initial density reduction, the approximate extension of which is obtained from their masses and the magnitude of the viscosity. This initial lowering of the density is assumed to be axisymmetric; the radial profile Σ(r) of the initial distribution is displayed in Fig. 1 (solid line). The total mass in the disc depends on the physical extent of the computational domain. Here we assume a total disc mass within r_min=0.25 and r_max=4.00 of 0.01 M_⊙. The starting model is then evolved in time and the accretion rates on to the planets are monitored, where a given fraction of the mass in the inner half of the Roche lobe of the planets is assumed to accrete on to the planets at each time-step and is taken out of the computation and added to the masses of the planets. The Courant condition yields a time-step of 6.8×10⁻⁴P₁.

Figure 1.

The azimuthally averaged surface density for different times (stated in units of P₁). The solid line refers to the initial density distribution. The density inside planet 1 is reduced because of mass leaving through r_min.

Open in new tabDownload slide

3 Results

Starting from the initial configuration (Fig. 1), the planets exert torques on the adjacent disc material which tend to push mass away from the location of the planets. At the same time, the planets continuously accrete mass from their surroundings, and the mass contained initially between the two planets is quickly accreted by the planets and added to their mass. Finally one large gap remains in the region between r=1 and 2 (Fig. 1). The duration t_d for this ‘gap clearing’ process can be estimated by considering the time that it takes to transport material with a given angular momentum a certain distance Δ given the torque. Bryden et al. (1999) estimate this time to be of the order of

(4)

where r₀ is the location of the planet having mass ratio q that orbits the star with a period P. For a typical half-width of the gap Δ=0.2r₀ and a mass ratio q=0.001, one finds for the clearing time t_d≈320P. The time to clear the gap as inferred from Fig. 1 agrees very well with this estimate. As has been seen already in previous investigations, the gap formation in an otherwise unperturbed initial disc begins after only a few orbital periods (Bryden et al. 1999; Paper I). This initial clearing time is independent of the viscosity in the disc, and for the given value of q this time is much smaller than the typical viscous time-scale t.

Similarly to the one-planet calculations described in Bryden et al. (1999), Paper I and Lubow et al. (1999), each of the two planets creates a spiral wave pattern (trailing shocks) in the density of the disc. In the case of one disturber on a circular orbit, the pattern is stationary in the frame corotating with the planet. The presence of a second planet makes the spirals non-stationary, as is seen in the snapshots after 50, 100, 250 and 500 orbits of the inner planet that are displayed in Fig. 2. Near the outer boundary at r=4 the reflection of the spiral waves is visible. Using the higher resolution model with r_max=8.0 (Section 2.2), we have tested whether the numerical resolution or the wave reflection at r_max has any influence on the calculation of the net torques acting on the two planets and the accretion rates on to them. Owing to limited computational resources, the higher resolution model was run only for a few hundred orbital periods, and the largest difference (2.5 per cent) occurred in the mass m₃ of the outer planet. The difference in radial position (migration) is less than 1 per cent. We may conclude that our resolution was chosen to be sufficiently fine, and that the reflections at the outer boundary r_max=4 do not change our conclusions significantly.

Figure 2.

A grey-scale plot of the surface density at four different times in an ϕ−r diagram. The small ellipses indicate the positions of the planets and the sizes of their Roche lobes.

Open in new tabDownload slide

In previous calculations (Paper I) the equilibrium mass accretion rate from the outer part of the disc on to a 1 M_Jup planet for the same viscosity (ν=10⁻⁵) and distance from the star was found to be 4.35×10⁻⁵M_Jup yr⁻¹ for a fully developed gap. Here the accretion rate on to the planets is much higher in the beginning (≈5×10⁻⁴M_Jup yr⁻¹) as the initial gap was not as well cleared. Thus, during this gap clearing process, the masses of the individual planets grow rapidly at the onset of the calculations (Fig. 3). At t≈250 the mass within the the gap has been exhausted (see also Fig. 1) and the accretion rates on to the planets lower dramatically. At later times, after several hundred orbits (P₁), they settle to nearly constant values of about 10⁻⁵M_Jup yr⁻¹ for the outer planet and 2.2×10⁻⁶M_Jup yr⁻¹ for the inner planet (from Fig. 3). Since the mass inside planet 1 has left the computational domain and the initial mass between the two planets has been consumed by the two planets, this mass accretion rate on to planet 1 for later times represents the mass flow of material coming from radii larger than r₂ (beyond the outer planet). It is the material that has been flowing across the gap of the outer planet. Previously, this mass flow across a gap has been calculated to be about 1/7 of the mass accretion rate on to a planet (Paper I), and the present results are entirely consistent with that estimate.

Figure 3.

The increase of the masses of the planets in units of M_Jup.

Open in new tabDownload slide

The gravitational torques exerted by the disc lead to an additional acceleration of the planets, resulting in an expression similar to the acceleration of the star (equation 2). For one individual planet this force typically results in an inward migration of the planet on time-scales of the order of the viscous time of the disc (Lin & Papaloizou 1986). Through their detailed two-dimensional computations, Nelson et al. (2000) estimate this migration time to be of the order of 10⁴ initial orbits of the planet. Here, this inward migration (with a similar time-scale) is seen clearly for the outer planet in Fig. 4, where the time evolution of the semimajor axes of the two planets is plotted.

Figure 4.

The evolution of the radial distance between the star and the planets. For times longer than t=1000, the two planets appear to be in a 2:1 resonance.

Open in new tabDownload slide

The inner planet, on the other hand, initially for t<200 moves slightly inwards, but then the semimajor axis increases and, showing no clear sign of migration anymore, settles to a constant mean value of 1.02. This different behaviour of the two planets is caused by the contributions of various parts of the disc to the direction and magnitude of the migrations (see also Lubow et al. 1999). The parts of the disc lying beyond the orbital distance of a planet cause, for the given planetary masses, an inward migration, and the inner parts of the disc cause an outward migration. As the inner planet has neither inner nor outer disc remaining, it does not migrate anymore, in contrast to the outer planet which moves inward, driven by the torques created by the outer disc material.

This decrease of the semimajor axis of the outer planet reduces the orbital distance between the two planets. From three-body simulations (a star with two planets) and analytical considerations (Gladman 1993; Chambers, Wetherill & Boss 1996) it is known that, when the orbital distance of two planets lies below the critical value of where

(5)

is the mutual Hill radius of the planet, the orbits of the planets are no longer stable. In the calculations this effect is seen in the strong increase of the eccentricity of the inner planet. At t=2500 its eccentricity has grown to about e₁=0.1, while the eccentricity of the outer planet remains approximately constant at a level of e₂=0.03.

We should remark here that, in the pure three-body problem without any disc and the same initial conditions (r₁=1.0, m₁=1.0; r₂=2.0, m₂=1.0) for the three objects, the semimajor axes of the planets stay constant as this system is definitely Hill-stable (Gladman 1993). However, if one takes as initial conditions for the pure three-body system the parameters for the planets as obtained from the disc evolution at t=2500 (r₁=1.0, m₁=2.3; r₂=1.5, m₂=3.2) then the evolution becomes chaotic on time-scales of hundreds of orbits, and the eccentricity grows up to e=0.6 for both planets within 4000 orbits.

4 Conclusions

We have presented calculations of the long-term evolution of two embedded planets in a protoplanetary disc that cover several thousand orbital periods. The planets were initially located at 1a_Jup and 2a_Jup from the central star with initial masses of 1 M_Jup each. The gravitational interaction with the gaseous disc, having a total mass of 0.01 M_⊙, leads to an inward migration of the outer planet, while the semimajor axis of the inner planet remains approximately constant and even increases slightly. At the same time, the ongoing accretion on to the planets increases their masses continuously until, at the end of the simulation (at t=2500), the outer planet has reached a mass of about 3.2 M_Jup and the inner planet has reached a mass of about 2.3 M_Jup.

This increase in mass and the decreasing distance between the planets eventually causes the orbits to become unstable, resulting in a strong increase of the eccentricities on time-scales of a few hundred orbits.

From the computations we may draw three major conclusions.

• (i)

  The inward migration of planets immersed in an accretion disc may be halted by the presence of additional protoplanets located, for example, at larger radii. After having grown sufficiently in mass, they disturb the density distribution significantly by forming, for example, one large joint gap. This in turn reduces the net gravitational torque acting on the inner planet. Thus the migration of the inner planet is halted, and its semimajor axis remains nearly constant.

• (ii)

  When disc depletion occurs sufficiently rapidly to prevent a large inward migration of the outer planet(s), a planetary system with massive planets at a distance of several au remains. This scenario may explain why, for example, in the Solar system the outer planets (in particular Jupiter) have not migrated any closer to the Sun.

• (iii)

  If the initial mass contained in the disc is sufficiently large then the inward migration of the outer planet(s) will continue until some of them reach very close spatial separations. This will lead to unstable orbits, resulting in a strong increase of the eccentricities. Orbits may cross, and then planets either may be driven to highly eccentric orbits or may leave the system altogether (see e.g. Weidenschilling & Marzari 1996). This may then explain the high eccentricities in some of the observed extrasolar planets, in particular the planetary system of ν Andromedae.

As planetary systems containing more than two planets display different stability properties (Chambers et al. 1996), it will be interesting to study the evolution of multiple embryos in the protoplanetary nebula.

Acknowledgments

I thank Dr F. Meyer for lively discussion on this topic. Thanks also to an unknown referee who helped to improve the quality of the presentation. Computational resources of the Max Planck Institute for Astronomy in Heidelberg were available for this project, and are gratefully acknowledged. This work was supported by the Max-Planck-Gesellschaft, Grant No. 02160-361-TG74.

References

Author notes