A&A 450, 833-853 (2006)
DOI: 10.1051/0004-6361:20054551

On the evolution of multiple protoplanets embedded in a protostellar disc

P. Cresswell - R. P. Nelson

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

Received 19 November 2005 / Accepted 17 January 2006

Abstract
Context. Theory predicts that low mass protoplanets in a laminar protostellar disc will migrate into the central star prior to disc dispersal. It is known that protoplanets on orbits with eccentricity $e \ga H/r$, where H is the disc scale height and r is the radius, can halt or reverse their migration.
Aims. We examine whether a system of interacting protoplanetary cores can excite and sustain significant eccentricity of the population, allowing some planetary cores to survive in the disc over its lifetime.
Methods. We employ two distinct numerical schemes: an N-body code, adapted to include migration and eccentricity damping due to the gas disc via analytic prescriptions, and a hydrodynamics code that explicitly evolves a 2D protoplanetary disc model with embedded protoplanets. The former allows us to study the long term evolution, the latter to model the systems with greater fidelity but for shorter times.
Results. After a brief period of chaotic interaction between the protoplanets that involves scattering, orbital exchange, collisions and the formation of co-orbital planets, we find that the system settles into a quiescent state of inward migration. Differential migration causes the protoplanets to form a series of mean motion resonances, such that a planet is often in resonance with both its interior and exterior neighbours. This helps prevent close encounters and leads to the protoplanetary swarm, or subgroups within it, migrating inward at a uniform rate.

In about $2 \%$ of runs a single planet is scattered onto a distant orbit with significant eccentricity, allowing it to survive in the disc for ${\sim } 10^6$ years. Over $20 \%$ of runs produce co-orbital planets that survive for the duration of the simulation, occupying mutual horseshoe or tadpole orbits.
Conclusions. Disc-induced damping overwhelms eccentricity growth through planet-planet interactions, such that a protoplanetary swarm migrates inward. We suggest co-orbital planets may be observed in future exoplanet searches.

Key words: stars: planetary systems: protoplanetary disks - planets and satellites: formation - stars: planetary systems: formation

   
1 Introduction

Of the 147 known extrasolar planetary systems, 17 are believed to be multiple-planet systems. As observational techniques improve and time-series of existing data lengthen, we can expect the number and diversity of multiple planet systems to increase substantially. Future missions such as KEPLER and COROT will extend the discovery-space of extrasolar planets further, covering the sub-Neptune mass range where current observations typically cannot detect planets (e.g. Bonfils et al. 2005; Marcy et al. 2005; Rivera et al. 2005; Vogt et al. 2005; Santos et al. 2004; Butler et al. 2003, 2004).

Current theories of planet formation suggest that a multistage process leads to the final assembly of planetary bodies within a protoplanetary disc. In simple terms, the process can be described as: (i) coagulation and settling of dust grains, leading eventually to the formation of 1 km-sized planetesimals; (ii) runaway growth of planetesimals leading to the formation of $\sim$100-1000 km-sized planetary embryos; (iii) oligarchic growth of these embryos, which accrete the surrounding planetesimals, eventually forming planetary cores of $\sim$10 $M_\oplus$ in the giant planet region beyond $\sim$3 AU, or Mars-mass bodies in the terrestrial planet zone. The 10-15 $M_\oplus$ rock and ice cores that form in the giant planet region are expected to accrete gaseous envelopes and become gas giants if they form prior to gas-disc dispersal, or remain as ice-giants if they form at the end of the gas-disc lifetime (e.g. Bodenheimer & Pollack 1986; Pollack et al. 1996). The Mars-mass embryos in the terrestrial zone are expected to eventually collide and accumulate, forming an inner terrestrial planet system (e.g. Chambers & Wetherill 1998).

There are a number of unsolved problems associated with this basic picture, in particular when it is applied to giant planet formation. Perhaps the most pressing is the fact that gravitational interaction between the gaseous protoplanetary disc and the massive planetary cores causes them to undergo rapid inward migration on a time scale of 10⁵ years for bodies in the 10 $M_\oplus$ range (e.g. Goldreich & Tremaine 1980; Ward 1986; Ward 1997; Tanaka et al. 2002). This process is normally referred to as type I migration, and the time scale for inward drift is an order of magnitude shorter than the time scale for gas accretion onto a forming gas giant planet (Papaloizou & Nelson 2005), making it difficult to understand how gas giant planets form at all before their cores accrete onto the central star.

There is a substantial body of work that has examined the interaction between multiple embryos, protoplanets, or fully formed planets within protoplanetary discs. Direct numerical simulations have been used to examine the interaction and growth of planetary embryos ("oligarchs'') during the oligarchic growth stage (e.g. Kokubo & Ida 2000; Thommes et al. 2003). Kominami et al. (2005) studied the formation of sub-Earth mass protoplanets from a planetesimal disc including the effects of type I migration, but found each protoplanet seed migrated inward before other bodies of comparable mass could form. The formation of resonant systems of giant planets has been examined by Kley (2000), Masset & Snellgrove (2001), Snellgrove et al. (2001), Lee & Peale (2001), Nelson & Papaloizou (2002), Kley et al. (2005), often with application to particular systems such as GJ876 (Marcy et al. 2001). More recently Papaloizou & Szuszkiewicz (2005) have examined the resonant capture between two protoplanets in the 1-20 Earth mass range. Thommes (2005) recently studied the effects of a giant planet on migrating Earth mass cores. He concluded that the formation of the first planet may be the most important phase of forming a planetary system, due to the giant's ability to capture the smaller cores into resonance, preventing their rapid migration. The problem of how to form the first planet still remains, however.

One issue that until now has not received significant attention is the question of how a swarm of protoplanetary cores, of Earth mass and above, embedded in a protoplanetary disc will evolve. The oligarchic growth picture of giant planet core formation predicts that a number of putative cores will form approximately coevally within the disc, separated by $\sim$8 mutual Hill radii (Kokubo & Ida 2000; Thommes et al. 2003). Differential type I migration may cause these bodies to undergo close encounters, inducing gravitational scattering and the pumping of significant eccentricities. It is known that type I migration can be substantially slowed or even reversed when a protoplanet attains an eccentricity $e \ga H/r$, where H/r is the ratio of disc scale height to radius (Papaloizou & Larwood 2000). In this paper we address the question of whether mutual interactions within a swarm of protoplanetary cores can excite and sustain significant eccentricities of the bodies, preventing substantial type I migration for at least some of the cores over the disc lifetime.

To address this problem we have performed numerical simulations using two distinct approaches. First we adapted an N-body code to include the effects of migration and eccentricity damping of the protoplanets due to the disc. This code allowed us to examine the long time scale evolution of planetary swarms embedded in protoplanetary discs. Second, we performed 2D hydrodynamical simulations of many-planet systems embedded in gaseous discs. These could not be run for as long, but they allowed us to check the results of the modified N-body simulations. The results of these simulations suggest that gravitational interaction within a planetary swarm is ineffective at maintaining significant eccentricities, due to the very strong eccentricity damping induced by the disc.

The plan of this paper is as follows. In Sect. 2 we give the equations of motion. In Sect. 3 we describe the two numerical schemes and calibrate the hydrodynamics code. In Sect. 4 we describe the initial conditions, and in Sect. 5 we discuss the results obtained using the modified N-body scheme. In Sect. 6 we present the hydrodynamic simulations, and compare their results with those obtained by the modified N-body code. We give our conclusions in Sect. 7.

   
2 Equations of motion

A protoplanetary disc undergoing near-Keplerian rotation has a ratio of scale height to radius $H/r \le 10^{-1}$. It is therefore reasonable to work with vertically averaged quantities and assume no vertical motion, reducing the problem to two dimensions. We work in cylindrical coordinates $(r,\theta,z)$, with the origin located at the central star. The continuity equation takes the form

$$\frac{\partial \Sigma}{\partial t} + \nabla \cdot \left ( \Sigma \mbox{\boldmath$v$ } \right ) =0,$$ (1)

while the momentum equation has radial and azimuthal components

$$\frac{\partial \left ( \Sigma v_r \right ) }{\partial t} + ... ...}{\partial r} -\Sigma \frac{\partial \Phi}{\partial r} + f_r,$$ (2)

$$\frac{\partial \left ( \Sigma v_\theta \right ) }{\partial t... ...c{\Sigma}{r} \frac{\partial \Phi}{\partial \theta} + f_\theta.$$ (3)

In the above, $\Sigma$ denotes the surface density which is defined by

$$\Sigma=\int_{-\infty}^{\infty} \rho \mbox{\hspace{0.1cm}d}z,$$ (4)

while P is the vertically integrated pressure and f_r and $f_\theta$ are the vertically averaged viscous forces per unit volume in the radial and azimuthal directions, respectively. Full expressions for f_r and $f_\theta$are given in Nelson et al. (2000). The radial and azimuthal velocities are denoted v_r and $v_{\theta}$, and the corresponding angular velocity is given by $\Omega=v_\theta/r$. The gravitational potential $\Phi$ is given by

 

$$\Phi(\mbox{\boldmath$r$ }) = -\frac{GM_*}{r} - \sum_{p=1}^{N} \fr... ...rm {S}} \frac{\mbox{d$m($\boldmath$r$ }')}{r'^3} \mbox{\boldmath$r \cdot r$ }',$$ (5)

where M_* is the stellar mass, $\epsilon$ is a softening parameter, and the summations are over all protoplanets with masses m_p. The subscript "p'' denotes quantities evaluated at the location of the protoplanet. The latter two terms in Eq. (5) result from acceleration of the coordinate system due to the gravity of the protoplanets and the protostellar disc. The integral is performed over the surface area of the disc.

Each protoplanet experiences the gravitational acceleration from the central star, the other protoplanets and the protostellar disc. The equation of motion for each protoplanet is:

$$\frac{\mbox{d$^2$\boldmath$r$ }_p}{\mbox{d}t^2} = - \frac{GM... ...\vert^3} (\mbox{\boldmath$r$ }_p - \mbox{\boldmath$r$ }_{p'}) \mbox{\hspace{1.2cm}} - \sum_{p' = 1}^N \frac{G m_{p'}} {r_{p'}^3} \mbox{\boldmath$r$ }_{p'} - \nabla \Phi_{\rm {d}}$$ (6)

where

 

$$\Phi_{\rm {d}}(\mbox{\boldmath$r$ }) = - G \int_{\rm {S}} \frac{\... ... {S}} \frac{\mbox{d$m($\boldmath$r$ }')} {r^{'3}} \mbox{\boldmath$r \cdot r$ }'$$ (7)

is the gravitational potential of the disc. The integrals are again performed over the surface area of the disc. The final term in Eq. (7) is the indirect term due to the acceleration of the coordinate system by the disc; similar indirect terms from the protoplanets are contained in Eq. (6).

   
3 Numerical methods

Two distinct numerical schemes have been used: (i) a full hydrodynamic 2D disc model together with embedded planets; (ii) a much faster N-body code adapted to emulate the effects of orbital migration and eccentricity damping on the protoplanets due to the protoplanetary disc. We describe each in turn in the following sections. In Sect. 3.1.2 we present the results of test calculations that demonstrate the ability of the hydrodynamic code to produce results in reasonable agreement with previous results on migration and eccentricity damping.

   
3.1 Hydrodynamic scheme

We used a modified version of the grid-based code NIRVANA to conduct the hydrodynamic simulations (Ziegler & Yorke 1997). The code is based on the ZEUS-algorithm (e.g. Stone & Norman 1992), and uses the monotonic transport scheme (van Leer 1977) to calculate advection. The code has been applied to the study of disc-planet interactions in a number of previous publications, and further details may be found in Nelson et al. (2000).

We work in cylindrical polar coordinates $(r,\theta)$, and adopt a locally isothermal equation of state such that the pressure, P, and surface density, $\Sigma$ are related by

$$P = c_{\rm {s}}^2 \Sigma.$$ (8)

Here $c_{\rm s}$ is the isothermal sound speed, which is a fixed function of radius, and is computed according to $c_{\rm {s}}=\left(\frac{H}{r}\right) v_{\rm {K}}$, where $v_{\rm {K}}$ is the Keplerian velocity. The disc models all have aspect ratio H/r=0.05.

The disc models are viscous, and we adopt the "alpha'' prescription for the kinematic viscosity $\nu = \alpha c_{\rm {s}} H$ (Shakura & Sunyaev 1973). We set $\alpha = 5 \times 10^{-3}$.

The motion of protoplanets within the disc is integrated using a $5{\rm {th}}$-order Runge-Kutta scheme (Press et al. 1992). All bodies and the disc interact gravitationally, but disc self-gravity is not included as the discs considered here have a Toomre parameter $Q \ga 30$.

   
3.1.1 Time step control and boundary conditions

When computing the time step size at each time level we first calculate a "hydrodynamic time step'', $\Delta t_{\rm {h}}$, and an "N-body time step'', $\Delta t_{\rm {n}}$. The hydrodynamic time step is calculated according to the standard Courant-Friedrichs-Lewy criterion. The N-body time step, which is sensitive to the interplanetary forces, is determined according to:

$$\Delta t_{\rm {n}} = \frac{2 \pi}{400} \mathop{\min}_{i,j} \l... ...\vec r}_{ij} \vert^3}{G(m_i+m_j)}} \mbox{\hspace{2mm}} \right)$$ (9)

where ${\vec r}_{ij}$ is the distance between particles i and j. The time step adopted is then $\Delta t= \mathop{\min} (\Delta t_{\rm {h}}, ~ \Delta t_{\rm {n}})$Experimentation has shown that for the N-body system, this approach leads to energy conservation per time step of 1 part in 10⁹ for close encounters between protoplanets, and 1 part in 10¹⁴ for more distant encounters.

The boundary conditions we adopt are designed to prevent the reflection of density waves excited by the protoplanets from the radial boundaries of the computational domain. During the simulations we solve the equation

$$\frac{{\rm d}X}{{\rm d}t} = -\frac{X-X_0}{\tau} R(r),$$ (10)

where $X=X(r,\theta,t)$ is one of surface density $\Sigma$, radial velocity v_r or azimuthal velocity $v_{\theta}$. X₀ is the initial value of this quantity, and $\tau$ is the Keplerian orbital period at the inner/outer boundary. R(r) is a quadratic ramp function that ranges from 1 at the inner/outer boundary to 0 at a distance of 0.1/0.4 computational units inside the disc, respectively. (More information can be found at the website describing the comparison of hydrodynamic codes for which this function was introduced.) The inner boundary of the disc was set at 0.6. The outer boundary location was dependent on the number and separation of the protoplanets included, but was such that the distance between the outermost planet and the outer boundary radius was >0.4 so that the planet was not in the wave damping region.

   
3.1.2 Code calibration and resolution

The simulations must be of sufficient resolution to model the interaction between disc and protoplanets at Lindblad resonances. For sub-gap opening protoplanets on circular orbits, resonances associated with large azimuthal mode numbers $m \ge r/H$ tend to pile up near the planet, standing-off at a distance of ${\simeq} 2/3 H$ due to pressure support in the disk (e.g. Artymowicz 1993a). For protoplanets on eccentric orbits, eccentricity damping is mainly provided by interaction at co-orbital Lindblad resonances, which lie at the position of the planet's semi-majoraxis (Artymowicz 1993b). We have calibrated the hydrodynamic code by comparing results from simulations with analytic results in the literature.

In the following tests we demonstate the numerical convergence of NIRVANA in simulating orbital migration and eccentricity damping due to the surrounding disc, and establish which resolution is sufficient to model this with reasonable accuracy, while allowing long-term simulations to be performed on currently available computers. We define the migration time to be:

$$t_{\rm {m}} = \frac{r_p}{\dot r_p} = \frac{J}{2 \dot J}$$ (11)

where $\dot r_p$ is the radial velocity of the planet due to migration, J is the planet angular momentum, and $\dot J$ is the torque acting on the planet. We use the migration and damping rates given by Papaloizou & Larwood (2000) to compare with NIRVANA. The migration time is given by

$$t_{\rm {m}} = \frac{3.5 \times 10^5}{2} f_{\rm {s}}^{1.75} ... ...p} \right ) \left ( \frac{r_p}{1~\rm {~AU}} \right ) \rm {yr}$$ (12)

and the eccentricity damping time $t_{\rm {e}} \equiv e/{\dot e}$ is given by

$$t_{\rm {e}} = 2.5 \times 10^3 f_{\rm {s}}^{2.5} \left [ 1 +... ...p} \right ) \left ( \frac{r_p}{1\rm {~AU}} \right ) \rm {yr,}$$ (13)

where H/r_p is the disc aspect ratio, e is the planet eccentricity, and subscript "p'' indicates that quantities are evaluated at the planet location. $M_{\rm {J}}$ and $M_{\rm {d}}$ are the Jovian mass and disc mass contained within 5 AU. The factor $f_{\rm {s}} \equiv 2.5\epsilon/H$ results from the potential softening.

Figure 1: The inverse of migration rates vs. eccentricity for four planet masses. Top left panel - 1 $M_\oplus$, top right panel - 5 $M_\oplus$, lower left panel - 10 $M_\oplus$, lower right panel -15 $M_\oplus$. Three resolutions were considered: (150, 300) (dashed line), (300, 600) (dotted line) and (450, 900) (solid line). Also shown are migration rates predicted by Eq. (12), but multiplied by a factor of three (dot-dashed line).

Open with DEXTER

Equation (12) shows that the exchange of angular momentum from protoplanet to disc changes sign for eccentricities > 1.1H/r, indicating a possible reversal of migration. As a result of this discontinuity, we plot $1/t_{\rm {m}}$, the inverse of the migration time scale, rather than $t_{\rm {m}}$ in the following figures.

The code calibration simulations described here adopted a disc model with $\Sigma(r)= \Sigma_0 r^{-3/2}$, with $\Sigma_0$ defined such that the disc had 2 Jupiter masses within 5 AU. The planet was initially located at r=1, which for these test calculations was assumed to be equivalent to 1 AU. The gravitational softening parameter $\epsilon=0.5H$, unless stated otherwise. The planet orbital evolution due to disc forces was calculated explicitly.

	Figure 2: Eccentricity damping rates $t_{\rm {e}}$ vs. initial eccentricity. Protoplanet mass and grid resolution are the same as Fig. 1; Papaloizou & Larwood's prescription for e-damping is given as exact (i.e. no factor of 3 modification). Agreement is strong for ep> 0.07.
Open with DEXTER

Figure 1 shows the inverse of the migration time $t_{\rm {m}}$ versus eccentricity for a selection of planet masses and three different numerical resolutions. Also plotted for each mass is the inverse of migration time as obtained from Eq. (12) but with $t_{\rm {m}}$ multiplied by a factor of 3. We find Eq. (12) consistently predicts migration times that are faster than observed in the simulations by about a factor of three for low eccentricity, for the particular value of the gravitational softening parameter adopted (the role of softening is discussed in Sect. 3.1.3). The two higher resolution simulations clearly demonstrate convergence towards a single solution. There is consistent discrepancy between the hydrodynamic simulations and the migration times predicted by Eq. (12) (but multiplied by a factor of three) at higher eccentricities, with the simulations predicting shorter time scales for angular momentum exchange by a factor of two, approximately. In other words, Eq. (12) predicts faster migration than the simulations produce by a factor of ${\sim} 3$ at low e, and by a factor of ${\sim} 3/2$for high e.

An interesting feature of our test simulations is that we do not observe the semi-major axes of the planets to increase when the eccentricity is e > 1.1H/r, although we do see a reversal in the angular momentum exchange. This is because the loss of energy by the planet due to eccentricity damping is more than able to counterbalance the gain in angular momentum, preventing outward migration as measured by an increase in semi-major axis. For large values of $e \simeq 0.5$, we observe the migration to essentially stall until the eccentricity falls to smaller values, after which inward migration ensues. These two observations indicate that type I migration can be prevented or reversed by the excitation and maintenance of significant eccentricity.

Eccentricity damping rates are presented in Fig. 2. Overall, we find very good agreement between our simulations and the predictions of Eq. (13). The only region of serious discrepancy is for low eccentricities below $e \le 0.08$. Here we find that Eq. (13) predicts faster eccentricity damping than the hydrodynamic simulations, being about a factor of three faster for $e \le 0.05$. Comparing the hydrodynamic simulations with each other, we again find good agreement between the medium and high resolution simulations, with the low resolution simulation predicting slower eccentricity damping rates.

Due to the similarity of the results obtained using (300, 600) and (450, 900) grid cells, and the long run times required for the simulations, it was decided to adopt the former resolution. Note that these values are for a disc with a radial extent of 0.4-2.5 as used in single-body tests; models with differing distributions of protoplanets require larger discs, and the mesh was adjusted to keep the same resolution between models.

   
3.1.3 Influence of the softening parameter

We use softening when calculating the gravitational interaction between disc and protoplanets, as shown in Eqs. (5) and (7). Higher values of $\epsilon$ will typically produce longer migration times for a given numerical resolution (Korycansky & Pollack 1993). The parameter can thus be used to adjust migration rates in 2D models, and an appropriately chosen value of $\epsilon$ will give migration rates that are consistent with those obtained for 3D disc models, where the vertical structure provides a natural softening of the disc potential (e.g. Tanaka et al. 2002; Nelson & Papaloizou 2004).

We have performed a number of test simulations to examine the influence of gravitational softening on migration and eccentricity damping rates, and to compare the obtained rates with values presented in the literature. We have also examined the effect of excluding material close to the planet from contributing to the force exerted by the disc on the planet. Here we excluded material from within a region of radius $R_{\rm {H}}$, the planet Hill radius, and from $\sqrt{2/3} R_{\rm {H}}$, $\sqrt{1.5} R_{\rm {H}}$ and $\sqrt{2} R_{\rm {H}}$. The rationale for this is simply that the simulations presented here are unable to accurately resolve material which enters the planet Hill sphere and becomes bound to the planet. If this bound material is allowed to contribute to the force on the planet, then its asymmtrical distribution due to lack of resolution can cause it to exert a large torque on the planet. The results shown in Fig. 3 indicate that the precise size of this "torque cut-off'' region does not strongly influence the results.

Figure 3: Migration time (computed from Eq. (11)) vs. softening for a body of $\approx$5 $M_\oplus$ on a fixed circular orbit. Four torque cutoffs are represented: $\sqrt {2.0}$ (solid), $\sqrt {1.5}$ (dotted), 1.0 (dashed) & $\sqrt {2/3}\times R_{\rm {H}}$ (dot-dashed); the latter two overlap. Also plotted is Papaloizou & Larwood's prescription for migration rates multiplied by a factor of 3 (triple dot-dashed); Tanaka et al.'s migration rates in 2D & 3D discs (lower and upper horizontal long-dashed lines); Korycansky & Pollack's results with softening parameter = 10-4 (shaded triangle). The peak at ${\approx } 0.1R_{\rm {H}}$ arises when $\epsilon \simeq \Delta r/2$, such that the softening length coincides with a radial cell boundary.

Open with DEXTER

Zero softening can give rise to numerical instabililty in hydrodynamic simulations where the planet moves through the computational grid, so we have first examined how simulations with small softening parameters agree with published results. Figure 3 shows migration times versus softening parameters assuming zero eccentricity. The migration time obtained for softening $\epsilon=10^{-4}$ by Korycansky & Pollack (1993), who solved numerically the linearised equations that describe type I migration, is shown by the shaded triangle at the bottom left of the figure. The simulations we performed with small softening parameters show reasonable convergence with this value. Tanaka et al. (2002) provide equations for the migration rate of planets in 2D and 3D discs, based on their linear calculations that neglected softening of the planet potential. These give

$$t_{\rm {m}}=\frac{1}{2.32 + 5.656 \beta} \left(\frac{M_*}{m_p... ...a_p r^2_p} \right) \left(\frac{H}{r_p}\right)^2 \Omega_p^{-1}$$ (14)

and

$$t_{\rm {m}}=\frac{1}{2.7 + 1.1 \beta} \left(\frac{M_*}{m_p} \... ...a_p r^2_p} \right) \left(\frac{H}{r_p}\right)^2 \Omega_p^{-1}$$ (15)

respectively for 2D and 3D discs, where $\beta$ is the power-law index for the surface density profile. The migration time given by Eq. (14) is shown by the lower horizontal dashed line, and that given for 3D discs by Eq. (15) is given by the upper horizontal dashed line. Additionally, the predictions obtained from multiplying Eq. (12) by a factor of three are shown by the dot-dot-dot-dashed line. It is clear from Fig. 3 that good agreement is obtained between the results for 3D discs presented by Tanaka et al. (2002), the hydrodynamic simulations, and the results obtained by Papaloizou & Larwood (2000) scaled by a factor of three, when the softening parameter $\epsilon \simeq 0.5 - 0.6 H$. All subsequent hydrodynamic models thus used $\epsilon=0.5H$. It is worth commenting that the non linear dependence of migration rate on the softening parameter given by Eq. (12) means that good agreement is obtained by this expression and Eq. (15) for $\epsilon \simeq H$.

We note that experiments conducted to examine the effect of softening on the eccentricity damping rates obtained using NIRVANA showed that the obtained rates were relatively insensitive to softening. Figure 2 shows the agreement obtained between the predictions of Papaloizou & Larwood (2000) and NIRVANA when $\epsilon=0.5H$.

   
3.2 Modified N-body scheme

The second method used to evolve a population of protoplanets is a simple N-body integrator with additional terms to model the effects of type I migration and eccentricity damping. The integrator used is the same $5{\rm {th}}$ order Runge-Kutta routine as used in NIRVANA. To best compare the two schemes we restrict the N-body code to the $(r,\theta)$ plane. We consider three dimensional simulations in a forthcoming paper. To prevent the time step becoming too low, embryos were removed from the simulation if they fell within 10 Solar radii of the origin.

   
3.2.1 Migration and eccentricity damping prescriptions

We wish to adopt migration and eccentricity damping prescriptions for use in the N-body code that enable us to follow the evolution of protoplanets that potentially achieve high eccentricities, and which also enable us to easily examine the evolution of protoplanet swarms in disc models with differing surface density profiles. The intention is to model such systems over long time periods, so that we can effectively extend the run times of the full hydrodynamic simulations. For the migration time we have adopted the prescription given by Tanaka et al. (2002), as written in Eq. (15), but modified to include the factor obtained by Papaloizou & Larwood (2000) to describe the evolution for large eccentricity:

$$t_{\rm {m}} = \frac{2}{2.7+1.1\beta} \left(\frac{M_*}{m_p} \r... ...t ( \frac{e_p r_p}{1.1 H} \right )^4} \right ) \Omega_p^{-1}.$$ (16)

The test calculations presented in Sect. 3.1.2 suggest that this prescription should give excellent agreement with the hydrodynamic simulations for small e, and good agreement for large e. Similarly, we wish also to adopt an eccentricity damping prescription that gives good agreement with the hydrodynamic simulations, and also allows us to calculate the evolution of protoplanets in disc models with different surface density profiles in a simple way. We thus adopt the formula obtained by Tanaka & Ward (2004) for the eccentricity evolution of low mass planets in 3D discs appropriate to small eccentricity. To model the evolution at high eccentricity we modify this formula by including the factor obtained by Papaloizou & Larwood (2000) which causes eccentricity damping to decrease at high eccentricity. The eccentricity damping time $t_{\rm {e}}$ is then given by

$$t_{\rm {e}} = \frac{Q_{\rm e}}{0.78} \left(\frac{M_*}{m_p}\ri... ... \left ( e_p \frac{r_p}{H} \right )^3 \right ] \Omega_p^{-1}.$$ (17)

Here $Q_{\rm e} \simeq 0.1$ is a normalisation factor that ensures that Eq. (17) gives good agreement with the eccentricity damping times obtained in the hydrodynamic simulations. To give maximum agreement with NIRVANA (see Fig. 2) we apply a lower bound to eccentricity damping, such that if a protoplanet has e_p < 0.08, the applied damping rate is what would be experienced if the body possessed an eccentricity of 0.08.

We implement the following expressions in the N-body code as accelerations experienced by the protoplanets due to the disc, using values of $t_{\rm {m}}$ and $t_{\rm {e}}$ obtained from Eqs. (16) and (17):

$${\vec a}_{\rm {m}}= - \frac{{\vec v}}{t_{\rm {m}}},$$ (18)

$${\vec a}_{\rm {e}} = -2 \frac{({\vec v} . {\vec r}) {\vec r}}{r^2 t_{\rm {e}}}\cdot$$ (19)

   
4 Initial conditions and model units

We first defined planetary initial conditions for the modified N-body code, and then certain models were rerun with NIRVANA using the same parameters. We first define the semi-major axis of the innermost body to be a₁=5 AU. Moving outward, successive planet semi-major axes were determined by choosing mutual separations to be a specified number of mutual Hill radii. Thus, $a_{i+1}=a_i+N_{\rm mH}R_{\rm mH}$where $N_{\rm mH}$ usually took the value 4 or 5 (but other values were considered), and the mutual Hill radius is defined by

$$R_{\rm {mH}}=\left( \frac{m_{\rm {i}}+m_{\rm {j}}}{3M_*} \right) ^ {\frac{1}{3}} \left( \frac{a_i + a_j}{2} \right) \cdot$$ (20)

We note that for two planets on initially circular orbits, rapid instability occurs if the separation $\Delta$ between planets is less than the critical value

$$\frac{\Delta_{\rm {crit}}}{R_{\rm {mH}}} = 2\sqrt{3} \approx 3.46$$ (21)

(Gladman 1993; Chambers et al. 1996). Simulations of oligarchic growth (e.g. Kokubo & Ida 2000; Thommes et al. 2003) suggest that the mutual separation of protoplanets is normally ${\simeq }8~R_{\rm mH}$. We adopted smaller spacings to maximise close encounters.

When setting the protoplanetary masses, our normal procedure was to define the mass of the innermost body (usually m₁=2  $M_\oplus$), with subsequent bodies having m_i+1 = m_i+2  $M_\oplus$. This rather artificial set up was chosen to maximise convergent migration in the planetary swarms, and hence to maximise interactions between the bodies. Different randomised mass distributions, however, were also considered, and are described in the relevant results sections below.

Planetary eccentricities were determined by defining a mean eccentricity for the planetary swarm, $\mu_{\rm e}$, and a standard deviation $\sigma_{\rm e}$, with eccentricities being randomly chosen according to a Gaussian distribution. Each body is given a random longitude of pericentre.

The distribution of semi-major axes, masses, and values of $\mu_{\rm e}$ and $\sigma_{\rm e}$ define a class of model. For each model class five realisations of the initial data were generated by modifying the random number seeds, giving rise to five different simulations. This allows us to gauge which trends in the results are due to changes in global parameters, and which are due to the stochastic nature of the problem.

The disc was initialised with scale height H/r=0.05 and $\Sigma(r) = \Sigma_0 r^{-0.5}$, where $\Sigma_0$ was defined such that the disc contains 40 Jupiter masses within 40 AU. In NIRVANA simulations the disc had a viscous parameter of $\alpha = 5 \times 10^{-3}$.

Computational units were adopted such that the central mass M_*=1corresponds to the Solar mass, the gravitational constant G=1, and the radius r=1 in the computational domain corresponds to 5 AU. The inner boundary of the computational domain for the hydrodynamic simulations was at r=0.6, and the outer boundary was set at a distance > 0.4 computational units beyond the outermost planet. We report our results in the physical units of years and astronomical units.

5 Results of modified N-body simulations

We have performed more than 300 modified N-body simulations to examine the evolution of clusters of low mass protoplanets embedded in a protoplanetary disc. These models involved varying parameters such as initial planetary separation, initial eccentricities, initial mass distribution, disc mass etc. In spite of all of this parameter variation, a number of simple robust outcomes were observed. We present a few individual simulations below that characterise the range of outcomes observed.

5.1 The fiducial model 

Figure 4: Evolution of a 10 planet N-body model with the fiducial setup (see main text). Top: semi-major axes of the migrating embryos. A brief period of activity is followed by a long period of migration between bodies in first-order mean-motion resonances. Bottom: the embryos' eccentricities over the same time. Damping from the disc prevents eccentricities from growing outside periods of intense activity.

Open with DEXTER

We choose one particular model to act as the fiducial case, against which other models will be compared in subsequent sections. This model consisted of ten protoplanets, each separated by five mutual Hill radii. The innermost planet had mass m_p=2  $M_\oplus$, with the planet mass increasing by 2  $M_\oplus$ as one moves out through the initial swarm, such that the outermost planet was 20  $M_\oplus$. This maximises convergent migration. The initial eccentricities were such that the mean value was $\mu _{\rm e}=0.05$ and standard deviation $\sigma_{\rm e}=0.01$. The variation of semi-major axes and eccentricities for this model are shown in Fig. 4 for a run time of ${\simeq} 2 \times 10^{5}$ yr. The evolution of the system proceeds as follows.

At the beginning of the simulation all protoplanets begin to migrate inward. The differential migration caused by the variation in protoplanetary masses leads to a convergence of the orbits, which in turn leads to the formation of mean motion resonances (MMRs) between pairs of protoplanets. Due to their closer proximity, the inner pair of planets are the first to enter an MMR, followed by the third planet catching up with the second planet and entering an MMR with it. The fourth planet then catches up with the third planet, and enters a resonance with it, and so on. This series of "stacked'' MMRs causes modest excitation of eccentricities, and may push some planets through a sequence of mean motion resonances (e.g. 6:5, 7:6. 8:7) such that they undergo a close encounter with their interior neighbour. This leads to scattering of some of the protoplanets early on in the simulation after a few $\times 10^4$ yr.

In Fig. 4, this scattering leads to a physical collision between the 8 and 10  $M_\oplus$ bodies after $t \approx 1.5 \times 10^4$ yr, causing them to merge and form a new 18  $M_\oplus$ protoplanet. Thus, the collision causes one body to be removed from the disc, creating additional space between the embryos. Collisions between bodies have a stabilising effect in two ways: first because purely inelastic collisions tend to damp the eccentricity of the resulting body; second because the removal of a body reduces mutual interactions between embryos.

After about $1.5\times10^4$ yr into the run, the innermost and least massive of the embryos scatters during an encounter with the adjacent 4  $M_\oplus$ protoplanet, exchanging positions with it. Subsequent interactions cause this 2  $M_\oplus$ embryo to be scattered inward due to close encounters with the newly formed 18  $M_\oplus$ protoplanet and the pre-exisiting 12  $M_\oplus$ body. It then collides with the 6  $M_\oplus$ protoplanet, forming a new 8  $M_\oplus$ body. The subsequent evolution is then characterised by inward migration of the protoplanet swarm, with the system in this case forming a series of MMRs between each of the planetary pairs, with the system being driven towards convergence in large part by the most massive and rapidly migrating 20  $M_\oplus$ protoplanet located at the outer edge. The resonances that form in this particular case are (moving from innermost pair to outermost pair): 8:7, 6:5, 5:4, 6:5, 5:4, 5:4, 5:4. These are broadly consistent with the resonances obtained by Papaloizou & Szuszkiewicz (2005). Overall, our simulations show that the 5:4 and 6:5 resonances are most favoured, where the resonance into which bodies become captured depends on the masses of the bodies, the rate of convergence of their orbits, and their eccentricities. The existence of alternating 5:4 and 6:5 resonances also implies the existence of a 3:2 resonance. The libration of resonant angles associated with the above commensurabilities are illustrated in Fig. 5. The resonant angles for a p:q first order resonance are defined by:

 

$$\phi_1 p \lambda_1 - q \lambda_2 - \omega_1$$ =

$$\phi_2 p \lambda_1 - q \lambda_2 - \omega_2$$ = (22)

where $\lambda_1$ and $\omega_1$ are the mean longitude and longitude of pericentre of the outer planet, $\lambda_2$ and $\omega_2$ are the equivalent for the inner planet. The planets in Fig. 5 are labelled from 1-10, with 1 being the innermost planet at the beginning of the simulation, and 10 being the outermost.

Figure 5: First order resonances between five of the bodies in Fig. 4. Each plot is labelled (P1, P2)-p:q, where P1/P2 are the exterior/interior protoplanets respectively, and p:q gives the commensurability between them; left/right columns show the resonant argument associated with the exterior/interior body, respectively. Each plot takes the same colour as the corresponding exterior protoplanet in that resonance. Note that the y-axis differs by $\pi$ in the left and right columns.

Open with DEXTER

A common feature of almost all simulations conducted is a brief period of readjustment of the protoplanet swarm, which may lead to scattering of embryos and their mutual collision, followed by orderly inward migration in which the whole or most of the system forms "stacked'' MMRs. The formation of the MMRs, combined with the strong eccentricity damping provided by the disc and the occurrence of collisional mergers prevent the system from developing into a state in which substantial eccentricities are maintained. The formation of planetary cores and embryos coevally with neighbouring bodies of similar mass apparently does little to help prevent large scale type I migration. The eventual fate of these systems appears to be inward drift into the central star, in the absence of a stopping mechanism such as an inner magnetospheric cavity (Lin et al. 1996). A similar situation is observed to occur in all "realistic'' simulations performed. Occasionally, a burst of strong interaction between migrating protoplanets is observed to occur in simulations as they approach the central star. An example is shown in Fig. 4 at $t \simeq 1.75 \times 10^5$ yr. This arises when a protoplanet slips out of resonance and undergoes a close encounter with a neighbour, resulting in a brief period of scattering and at least one merger or, occasionally, one body thrown onto an orbit $a_p \leq 4$ AU. The system, however, quickly returns to an orderly state of inward resonant migration.

5.2 Effect of increasing initial eccentricities

We have performed numerous calculations for which the average initial protoplanet eccentricities ranged between $0.1 \le e \le 0.3$, but with semi-major axes identical to those in the fiducial model described in the previous section. For initial eccentricities at the lower end of this range, the evolution overall was found to be similar to that presented in Fig. 4, but with a modest increase in the amount of scattering occurring between the protoplanets during the early stages of the runs. Very quickly, however, the protoplanets settle into MMRs involving most of the swarm, and undergo inward resonant migration. We discuss such a model below in more detail.

For a qualitative change in the evolution to occur, the initial eccentricities needed to be around $e \simeq 0.3$, leading to numerous bodies being orbit crossing at the beginning of the simulation. These models then display substantially more scattering among the protoplanets, with some being flung out to large ($r \ge 40$ AU) radii. Such initial configurations, however, would appear to be very unrealistic so we do not discuss them further.

	Figure 6: As Fig. 4, but for a model with a higher initial eccentricity distribution ( $\mu _{\rm e}=0.1$ compared to $\mu _{\rm e}=0.05$ previously). Resonant migration dominates again once the excited body has been scattered beyond the system of protoplanets. The two plots to the right detail the evolution of the main population.
Open with DEXTER

Figure 6 depicts the evolution of a model identical to the fiducial case described in the previous section, except for a larger initial eccentricity distribution, with $\mu _{\rm e}=0.1$ ( $\sigma_{\rm e}=0.01$ as before). Close inspection of Fig. 6 shows that the innermost protoplanet, with mass m_p= 2 $M_\oplus$ quickly interacts with its nearest neighbour, leading to them exchanging orbital locations. The 2 $M_\oplus$ protoplanet undergoes similar interactions with the other protoplanets and is quickly propelled outward through the system of embryos until it encounters the outermost protoplanets with masses m_p=18 and 20  $M_\oplus$, respectively. Interaction between the 2 and 18 $M_\oplus$ bodies leads to the 2 $M_\oplus$ object being scattered into an orbit with a_p=45 AU and e=0.75. This orbit crosses that of the 20 $M_\oplus$ protoplanet, and subsequent interaction between these two causes the orbit of the scattered planet to be modified such that a_p=35 AU and e=0.67. During this period of activity, and as an indirect consequence of it, the 8 and 10 $M_\oplus$ protoplanets experience a collisional merger, resulting in the formation of an 18 $M_\oplus$ body.

Once the initial readjustment of the system is complete, consisting of orbit exchange, collisional mergers, and scattering, the system settles into a state of inward resonant migration. Here, each neighbouring pair of protoplanets forms a MMR as the swarm migrates inward in lockstep (excluding the outermost scattered body at 35 AU), at a rate determined by the total angular momentum content of the swarm and the total tidal torque acting on the bodies due to the disc. The final arrangement of the protoplanets in Fig. 6 is similar to those shown in Fig. 4, the resonances being 7:6, 5:4, 5:4, 5:4, 6:5, 6:5 and 4:3 working outward from the interior. Again the resonances are predominantly 5:4 and 6:5, with a closer resonance at the leading edge of the swarm where the smallest remaining body is pushed ahead of the group. Although resonant interaction between protoplanets leads to eccentricity excitation, the disc provides sufficient damping that the eccentricities remain small ( $e \le 0.03$). This fact, combined with the resonant configuration which helps prevent close encounters, means that no further scattering between protoplanets occurs, and the system evolves inward, eventually migrating into the central protostar.

Considering both the models depicted in Figs. 4 and 6, only the scattered body in the latter has the possibility to survive beyond the dispersal of the disc to remain a bound planet. Its large eccentricity substantially slows inward migration, while the reduced surface density at large radii naturally produces a slower migration rate and damps the eccentricity more slowly. If the disc were to disperse within ${\sim } 10^6$ yr, it would be left as the sole surviving planetary body. It is interesting, however, to speculate how such a body would evolve immediately after being scattered. One possibility is that after its eccentricity has damped, it could accrete rapidly much of the surrounding solid material that would be in the form of smaller bodies. Without the competition of neighbouring bodies of similar size, it could grow to become a massive planetary core with mass >15 $M_\oplus$ and subsequently accrete gas becoming a gas giant. Such a scenario has been discussed recently by Weidenschilling (2005).

Our simulations strongly suggest that a group of protoplanetary cores that form coevally within a laminar protoplanetary disc will generally spiral into the central protostar, probably in resonance with some or all of the other bodies. Scattering may allow one or two relatively low mass, easily displaced bodies to survive, through relocation onto larger orbits, but survival still depends on the disc lifetime and the accretion history of the remaining planet. It seems clear that interactions between a cluster of protoplanets is unable to maintain sufficient eccentricity against disc damping such that type I migration can be reversed or substantially slowed.

5.3 Effect of decreasing initial protoplanet separations

Figure 7: As Fig. 4, but with initial separations $\Delta _0$ of 4 $R_{\rm {mH}}$ between the protoplanets. Resonant migration again dominates, but now with two pairs of co-orbital protoplanets. These two pairs lie in a 3:2 mean motion resonance, a consequence of the 14 $M_\oplus$ body displacing the inner pair at $3.1\times 10^4$ yr without breaking the commensurability. The bottom plot shows a detailed view of the active period following intial settling.

Open with DEXTER

We have performed numerous runs for which the initial mutual separation between protoplanets was reduced from 5 to 4 $R_{\rm {mH}}$. Overall, modestly decreasing the initial separation between the protoplanets has a similar effect to modestly increasing their eccentricities: a noticeable increase in the number of scattering events and orbital exchanges during the early stages of the evolution. Once again, however, resonant inward migration toward the central star ensues after this initial period of readjustment. One outcome of these runs, that is more common when the initial separations are smaller, is the formation of 1:1 resonances, in which two planets enter and maintain a configuration in which they are in mutual horseshoe or tadpole orbits. For initial separations $\Delta_0= 5$ $R_{\rm {mH}}$ such configurations occur in ${\simeq} 20$% of the runs. This figure rises to ${\simeq} 80$% when $\Delta _0$ is reduced to 4 $R_{\rm {mH}}$. Note these numbers apply to systems in which 1:1 resonances formed and survived for the duration of the run. Additional runs showed the formation of 1:1 resonances during the early stages of readjustment, but these sometimes did not survive.

Figure 7 shows the evolution of a model where the protoplanets were initially separated by 4 $R_{\rm {mH}}$ , but was otherwise identical to the fiducial run described in Sect. 5.1. A striking feature is that resonant migration dominates at late times, and the change in separations has not altered the ultimate fate of the bodies: absorption by the protostar. Closer inspection shows the 14  $M_\oplus$ (orange line) body migrating through the protoplanetary swarm through a sucession of orbital exchanges with interior protoplanets. These exchanges occur when two protoplanets enter each others' horseshoe regions, such that a close encounter causes the protoplanets to effectively swap orbits.

Another striking feature shown in Fig. 7 is the formation of 1:1 resonances by pairs of protoplanets. One example is the formation of such a configuration by the 12  $M_\oplus$ (green) and 16  $M_\oplus$ (purple) protoplanets at $\approx$ $2.7\times10^4$ yr; the smaller body is scattered into the horseshoe region of the larger, resulting in the horseshoe libration of these two protoplanets. Over a period of ${\la} 10^4$ yr, the libration amplitude decreases due to disc interaction, causing the planets to become more stably trapped in the 1:1 resonance. The horseshoe orbits eventually become tadpole orbits, with the protoplanets becoming "trojans'' that librate with small amplitude about the mutual L₄ and L₅ points. This explains why the purple line in Fig. 7 becomes obscured by the green line.

Another co-orbital resonance forms between two 8 $M_\oplus$ bodies, one of which is formed through collision at ${\simeq} 7.2 \times 10^3$ yr. Rather than following the gradual reduction of libration width described above however, in this instance the inner body is scattered almost directly into the path of the other, and the two begin small amplitude librations almost immediately. This arrangement is so stable that after the approaching 14 $M_\oplus$ body displaces the pair to larger radii, they immediately resume their co-orbital status. A small reduction in libration width is visible as the pair settle back into their mutual L₄/L₅ points at ${\sim} 3.1 \times 10^4$ yr.

Trojan planets have much the same effect as collisions in speeding the onset of resonant migration, effectively reducing the radial number density of bodies in the disc by having two bodies share very similar orbits. In this instance, the collisional merger (which removes the excited lowest mass protoplanet) and the formation of two co-orbitals, leaves only seven distinct orbits remaining, and resonant migration ensues once the 14 $M_\oplus$ body's inward journey through the population stops. Despite the initially reduced separations, the resonances are essentially unchanged from previous models, with (starting at the inner edge) 6:5, 5:4, 5:4, 6:5, 5:4 and 6:5 MMRs between each pair. Note the alternating 5:4/6:5 structure produces four additional 3:2 resonances, one of which is between the two co-orbital pairs. This run is typical of other simulations with $\Delta _0=4$ $R_{\rm {mH}}$, which produced very similar resonant structures to the previously considered models. The only significant difference between models with the two different separations is the profusion of long-lived trojan planets.

5.4 Effect of increasing initial separations

Simulations of the earlier stages of planet formation, including the oligarchic growth phase and beyond, indicate that massive planetary embryos tend to form with spacing in the range 7-10  $R_{\rm {mH}}$. It is not clear whether the abundance of 5:4 and 6:5 MMRs that arise in the simulations discussed so far are a result of the initial close proximity of the protoplanets we have adopted. We have performed numerous simulations with larger separations between the protoplanets, with values ranging from 6-10  $R_{\rm {mH}}$, to investigate the outcome in this more realistic situation.

Separations of 6 $R_{\rm {mH}}$ and greater allow differentially migrating protoplanets to pass through more distant first-order resonances such as 2:1 and 3:2. Nonetheless, neighbouring planets were never observed to become trapped in these. The 5:4, 6:5 and 7:6 resonances were still favoured because of the masses and convergence rates of the bodies. This is in basic agreement with Papaloizou & Szuzskiewicz (2005) who find formation of 2:1 and 3:2 resonances only when the planet masses are very similar, reducing the differential migration rate. Overall the runs with 6 $R_{\rm {mH}}$ to 10 $R_{\rm {mH}}$ separations evolved in a similar fashion to the fiducial run described in Sect. 5.1, with only a modest decrease in the amount of scattering, strong interactions and collisions observed.

5.5 Effect of changing mass distribution

The assumption that the protoplanetary masses increase monotonically as one moves outward in the disc, adopted by the runs discussed so far, is clearly open to question. Such a scenario is unlikely to occur in reality, and was adopted in our previous runs to maximise the convergence between neighbouring planets and increase the likelihood of strong interactions and scattering. In fact, it may be argued that if planet formation proceeds more rapidly in inner regions of the disk, where the dynamical time is shorter, then perhaps one should assume that protoplanet masses increase as one moves inward. This scenario, however, would lead to little interaction between the bodies as type I migration would cause the inner planets to runaway from the outer ones. Given that planetary accretion is a stochastic process, it is not unreasonable to examine randomised mass distributions.

We have performed numerous simulations in which the mass of the planet at each semi-major axis is chosen randomly from a Gaussian distribution with well defined mean and standard deviation. Typically we chose the mean to be in the range 5-10  $M_\oplus$, with the standard deviation being 3-6  $M_\oplus$, usually subject to lower and upper mass cut-offs of 2 and 20  $M_\oplus$. There are two significant differences that occur in these runs compared with those that have been described previously. The first occurs when two neighbouring bodies have significantly different masses, with the more massive protoplanet being the outer one. Rapid convergence of the orbits caused by strong differential type I migration in this case can lead to scattering rather than resonant trapping, and in general we observe more scattering and collisions in the runs with randomised mass distributions. A similar tendency for non gap forming planets with substantially differing masses to scatter was noted by Papaloizou & Szuszkiewicz (2005). Second, when the outermost planet is no longer the most massive and resonantly driving the rest of the swarm toward convergence, we find that inward resonant migration still occurs, but resonantly migrating planets now form distinct groups that migrate collectively at different rates.

Figure 8: As Fig. 4, but with a Gaussian mass distribution for the protoplanets. From the object at smallest initial radius, to the nearest 0.1 $M_\oplus$ the masses of the bodies are 13.8, 9.8, 17.6, 14.2, 8.5, 19.4, 2.9, 10.6, 14.7 & 6.6 $M_\oplus$. The non-monotonic distribution of migration rates causes the population to fragment into smaller resonant groups. The 2.9 $M_\oplus$ body (orange) is excited rather than captured by the neighbouring 10.6 $M_\oplus$ embryo (purple), but is subsequently able to form a stable resonance with another body of comparable mass.

Open with DEXTER

Figure 8 shows the results of one simulation in which the mass distribution has a mean of 10 $M_\oplus$ and standard deviation of 6  $M_\oplus$. The masses of each body are listed in the caption of Fig. 8. The 2.9  $M_\oplus$ (orange line) body is briefly trapped in the 8:7 MMR with the 10.6  $M_\oplus$ (purple) planet, but the resonance breaks at $t \simeq 2.5\times10^4$ yr. The lighter body is then scattered outward before being caught in a stable MMR with the outermost 6.6  $M_\oplus$ (yellow) body.

We note that the mass ratio of the 2.9 and 10.6 $M_\oplus$ bodies is $\approx$0.27. Examining our simulations as a whole, we find that the probability of scattering occurring between two neighbouring bodies, rather than prolonged resonant migration, is increased significantly when the two bodies possess a mass ratio $q_{\rm {m}} < 0.40{-}0.43$. For values of $q_{\rm {m}}$ below 0.4 scattering becomes almost certain. This is simply due to the enhanced differential migration rate that can cause planets to traverse the stable mean motion resonances and scatter after a close encounter. The lack of 9:8 MMRs suggests that the 8:7 MMR is the last stable resonance that the low mass planets can occupy without scattering, at least for the parameters that we study in this paper. The recent study of resonant capture by Papaloizou & Szuszkiewicz (2005) did occasionally find 9:8 and even 10:9 resonances forming, but the presence of additional planets in this work may render such configurations unstable.

Figure 9: As Fig. 4, but for a model with 20 protoplanets separated initially by 5 $R_{\rm {mH}}$. The protoplanets have reduced masses 1 $M_\oplus$, 1.5 $M_\oplus$, 2 $M_\oplus$,..., 10.5 $M_\oplus$, reflecting an earlier period of growth. The middle and bottom plots show the eccentricities of the initially inner and outer ten bodies, respectively. Collisions increase protoplanet masses and cause the formation of several smaller resonantly migrating groups.

Open with DEXTER

The subsequent evolution involves resonant inward migration. Figure 8 shows that the planets have divided into three distinct groups, migrating inward at different rates. The innermost group is being driven largely by the 19.4  $M_\oplus$ (green) protoplanet. The middle group of two planets migrates more slowly and is driven by the 14.7 $M_\oplus$ (cyan) protoplanet. The outer group migrates the slowest and is driven by the 6.6 $M_\oplus$ body. We note that scattering and collisional mergers in the previous runs described, where the initial mass distributions are monotonically increasing with orbital radius, do lead to periods of evolution during which groups of protoplanets migrate at different rates. In these cases, however, the presence of the most massive and rapidly migrating planet at the outer edge of the swarm ensures that - unless collisions have created a massive body in the middle of the swarm - these resonant groups are eventually driven to form a single resonantly migrating group.

5.6 Effect of changing protoplanet number

We have performed simulations in which we varied the number of protoplanets, running cases with five and twenty bodies rather than the default number of ten. In general, decreasing the number down to five results in less scattering and fewer collisions, with the system quickly settling into a mode of resonant inward migration.

Increasing the number of bodies to twenty results in more scattering and collisions during the early phase of readjustment, but the system again settles into a prolonged period of resonant, inward migration. Such a run is shown in Fig. 9. Runs with larger numbers of bodies indicate that stable resonant migration cannot occur for more than ${\simeq} 10{-} 12$ protoplanets in resonance. Instead, collisions and the creation of more massive bodies break up large resonant groupings into smaller clusters of protoplanets, each with a similar or fewer number of bodies to the models presented in other sections.

	Figure 10: As Fig. 4, but for a model with reduced (gaseous) disc mass. Eccentricity damping (and migration rates) are reduced by a factor of 5, but collisions reduce the total number of bodies while eccentricity damping is still too strong for mutual interactions to excite the population.
Open with DEXTER

	Figure 11: As Fig. 4, but for a model with eccentricity damping (Eq. (17)) reduced by a factor of 50. Scattering of only ten bodies is now self-sustaining over ${\sim } 10^6$ yr, redistributing bodies throughout the disc. Without energy losses to the disc, migration is reversed at even moderate eccentricities.
Open with DEXTER

5.7 Effect of decreasing disc mass

The main obstacle to sustaining the high eccentricities needed to halt migration in these models is two-fold: strong eccentricity damping by the disc, and inelastic collisions. In each of Figs. 4 and 6-9, no more than one body ever achieves - let alone sustains - an eccentricity greater than 0.1 (barring instantaneous increases - the large spikes - due to close encounters). With this seemingly the main obstacle to maintaining an excited population of embryos, a series of simulations were run with lower (gaseous) disc masses. If planet formation occurs late on during the life of the gas disc, after it has largely accreted onto the central star or is being dispersed through UV photoevaporation, then the gas-solids ratio could be significantly lower than assumed in our previous runs, decreasing damping and migration rates.

It was found that decreasing the surface density (disc mass) to $0.5\times$ the default model did produce more scattering and collisions in about twice as many models as before. The scattering, however, usually involved only two or three excited bodies, as found in previous models. Increasing initial eccentricities created sufficent early mixing so that the active period of the population was prolonged considerably in some models, though this extra activity, prolonged or otherwise, led to extra collisions and thus a greater number of isolated migrating bodies. Only when the disc mass was reduced to $0.2\times$ the mass of the default disc did considerable activity spontaneously develop from a low-eccentricity population; but once again this inevitably led to collisions and resonant migration. Such a run is shown in Fig. 10, and illustrates the potential importance of collisions in maintaining a dynamically cold population of protoplanets. We note that 2D runs overestimate the collision frequency, so this result may change in 3 dimensions. No evolution involving sustained scattering to high eccentricities and stalled migration was observed for disk mass reduction factors of $\ge$0.2.

5.8 Effect of reducing e-damping efficiency

To examine the role of eccentricity damping further in determining the evolution of a swarm of planetary embryos, we performed numerical experiments for which the eccentricity damping rates were artificially decreased by a factor $f_{\rm {d}}$ in Eq. (17). Ultimately we found that $f_{\rm {d}} \la 0.2$ was required before scattering occurred in quantities substantially above those seen for the unaltered rate, in agreement with the reduced disc mass models; otherwise resonant migration dominated completely. Below this value, eccentricities were frequently sustained at $e_p \sim 0.1$. Decreasing damping rates by $f_{\rm {d}} \la 0.1$ we found sustained eccentricities 0.1<e_p<0.3 to be common throughout the population, and multiple embryos were frequently moved onto larger orbits or ejected completely. This was also sufficent for prolonged orbit crossing even when separations between remaining embryos have been increased by the loss of other embryos from the swarm.

A run with $f_{\rm {d}}=1/50$ is shown in Fig. 11, showing prolonged excitation of the protoplanetary swarm, and long periods of outward migration arising because large eccentricities are maintained over long times. Because the disc is no longer effective at damping eccentricities, migration is able to reverse. Such a situation leads to lifetimes of $t \sim 10^6$ yr for bodies in the disc. This confirms that the main action of the disc is to prevent eccentricity growth among a population of protoplanets, leading to an absence of scattering and the ineffectiveness of mutual perturbations at raising eccentricities.

6 Results of hydrodynamic simulations

The results from the modified N-body simulations presented in previous sections show that a number of interesting phenomena arise in the evolution of a swarm of low mass protoplanets: a modest degree of scattering between protoplanets; orbital exchange between protoplanets; collisional mergers; large-scale resonant migration involving numerous bodies; formation of 1:1 resonances leading to "tadpole planets''; eventual migration of the swarm into the central protostar. We now present results of hydrodynamic calculations that directly simulate the evolution of similar protoplanetary swarms. Given that the prescriptions used in the modified N-body simulations, to describe the disc-planet interaction, were normalised to the results of hydrodynamic simulations, it is not surprising we find a similar range of phenomena emerging from the full hydrodynamic simulations. Inevitably, the hydrodynamic simulations are less conclusive about the long-term behaviour of the protoplanet swarms due to restrictions on the run times currently possible. Indeed, we are forced to halt most simulations at the point when the innermost protoplanet enters the damping region close to the inner boundary of the computational domain. This is located at a radius of 3.5 AU in the plots shown in the following sections.

   
6.1 A fiducial hydrodynamic simulation

The initial conditions of this simulation are designed to be as close as possible to those for the fiducial modified N-body simulation described in Sect. 5.1. The initial interplanetary separation was 5 $R_{\rm {mH}}$, and the initial eccentricies were Gaussian distributed with mean $\mu _{\rm e}=0.05$ and standard deviation $\sigma_{\rm e}=0.01$. The masses are distributed by setting the innermost protoplanet mass to 2  $M_\oplus$, and increasing subsequent masses as one moves out through the swarm by 2  $M_\oplus$, such that the outermost body has 20  $M_\oplus$. The evolution of the semi-major axis and eccentricity distribution is shown in Fig. 12, covering a run time of ${\simeq} 3.5 \times 10^4$ yr, which should be contrasted with the run times of > 10⁵ yr possible with the modified N-body runs.

	Figure 12: Evolution of a model of ten bodies using the NIRVANA hydrodynamic code, similar to Fig. 4. Collisions remove smaller, easily perturbed bodies as differential migration produces a successively more crowded system, until stable resonances form and the group migrate inwards at a single rate. As in the majority of N-body integrations, no body achieves and sustains ep>0.1.
Open with DEXTER

The evolution of the semi-major axes in the upper panel of Fig. 12 shows that the outcome is ultimately determined by resonant inward migration. The first $5 \times 10^3$ yr are characterised by the protoplanets undergoing differential inward migration, leading to convergence of their orbits and formation of mean motion resonances between planetary pairs. The system remains relatively quiescent until after ${\sim} 1.5 \times 10^4$yr when the 6  $M_\oplus$ (olive green) and 8  $M_\oplus$ (red) bodies collide during a period when the inner protoplanets undergo a burst of stronger interaction leading to a momentary increase of eccentricities up to $e_p \simeq 0.05{-} 0.08$. The resulting 14 $M_\oplus$ body resonantly drives the 4 $M_\oplus$ body inward, ultimately leading to a collisional merger between the two innermost bodies.

At the end of the run the system consists of 8 protoplanets undergoing resonant inward migration. The two innermost bodies form a separate resonance and are moving away from the larger group. The main MMRs in the population are found to be 5:4 and 6:5 commensurabilities, as found in the N-body runs, giving rise to a 3:2 resonance between the 14  $M_\oplus$ (orange line) and 12  $M_\oplus$ (green) bodies which are each adjacent to the 10  $M_\oplus$ (dark blue) protoplanet with which they are also in resonance. Another 3:2 resonance forms between the outermost 20 $M_\oplus$ (yellow) and 16 $M_\oplus$ (purple) bodies following the disturbance at ${\approx} 1.7 \times 10^4$ yr. Prior to collisions the innermost protoplanets were typically in 7:6 MMRs. These numbers are in broad agreement with the simulations performed by Papaloizou & Szuszkiewicz (2005) who find that low mass protoplanets tend to form commensurabilities in the range 4:3-8:7, depending on the rate of convergent migration. As found in the N-body results, the eccentricities typically remain below e < 0.03 except during close encounters, after which the disc quickly damps the eccentricty down to a background value of $e \le 0.03$. Thus, the combination of mean motion resonances and strong disc damping appears to prevent the system maintaining a dynamically excited state in which high eccentricities can prevent rapid inward type I migration. Although it cannot be claimed with certainty that the interesting phase of evolution is over in Fig. 12, the N-body runs suggest that this system will undergo large-scale inward resonant migration resulting in the protoplanets entering the central protostar.

Figure 13: As Fig. 12, but with higher initial eccentricities, similar to Fig. 6. All surviving bodies end the run in 5:4 and/or 6:5 resonances with adjacent protoplanets, while a co-orbital system forms and remains stable until the end of integrations, over ${>}2.93\times 10^4$ yr later. Rapid eccentricity damping is clearly visible early on in the 4 (light blue), 6 (olive green) and 2  $M_\oplus$ (black) bodies; only the former acheives ep > 0.1.

Open with DEXTER

6.2 A hydrodynamic run with larger initial eccentricity

We now present the results of a simulation which is very similar to that presented in the previous section, except that the eccentricities are Gaussian distributed with mean $\mu _{\rm e}=0.1$ and standard deviation $\sigma_{\rm e}=0.01$. The results are shown in Fig. 13, and one is again struck by the inward migration of the protoplanet swarm. Interestingly, the evolution in this case displays less scattering than the lower eccentricity run described in Sect. 6.1, and results in fewer collisional mergers. Inspection of the lower panel in Fig. 13 shows that the initial eccentricities are quickly damped by the disc.

	Figure 14: First order resonances between the protoplanets in Fig. 13. Each plot is labelled and coloured as in Fig. 5. Note that the trojan planets (P2=2,4) both share the two resonances with exterior bodies (P1=5,6), as shown in the lower four pairs of plots.
Open with DEXTER

The scattering activity that does occur in this run is completed within $6 \times 10^3$ yr. Initially we observe that the 2  $M_\oplus$ (black line) and 4  $M_\oplus$ (light blue) planets scatter and move apart, which ultimately results in the formation of a 1:1 resonance between the scattered 4 $M_\oplus$ protoplanet and the 8 $M_\oplus$ body. The initial horseshoe orbits that form quickly evolve to become tadpole orbits at L₄ and L₅. This system then remains stable over the full run time of the simulation.

A sole collision occurs later at ${\simeq} 2.4 \times 10^4$ yr between the 2 and 6  $M_\oplus$ (olive green) protoplanets. This is a result of the 6 $M_\oplus$ body being resonantly driven inward by the more rapidly migrating exterior planets after ${\simeq} 1.3 \times 10^4$yr, breaking the 5:4 resonance between the 2 and 6 $M_\oplus$ bodies.

Once the initial phase of scattering is complete, the system settles into a state of resonant inward migration. Examples of the resonances that are formed are shown in Fig. 14, which displays the two resonant angles associated with each first order resonance (defined by Eq. (22)). The planets in this plot are labelled from 1-10, with 1 being the innermost planet at the beginning of the simulation and 10 being the outermost. Note that the existence of alternating 5:4 and 6:5 resonances implies the existence of a 3:2 resonance between pairs of planets that are not adjacent to one another, but are separated by another planet. Examples of this are shown in Fig. 14. One particularly interesting feature of the resonances displayed in Fig. 14 is those involving the co-orbital planetary system, which consists of protoplanets 2 and 4. It can be seen that protoplanet 5 forms a 5:4 resonance with the tadpole planets 2 and 4. Protoplanet 6 forms a 6:5 resonance with planet 5. This results in 3:2 resonances forming between protoplanet 6 and tadpole protoplanets 2 and 4.

N-body simulation results like those presented in earlier sections strongly suggest that the long-term evolution of the protoplanets shown in Fig. 13 will be resonant inward migration into the protostar.

6.3 A hydrodynamic run with smaller initial separations

We performed a hydrodynamic run with initial protoplanet separations reduced from 5 to 4 $R_{\rm {mH}}$. In other respects this model set up is the same as the fiducial hydrodynamic run described in Sect. 6.1. The evolution of the semi-major axes and eccentricities are shown in Fig. 15. Comparing this figure with Figs. 12 and 13 shows that the system immediately shows greater interaction between the protoplanets. A number of features show close resemblance to phenomena observed in the modified N-body runs. First, the innermost planet undergoes a number of close interactions with exterior protoplanets, resulting in repetitive orbit exchange and eventually collisonal merger. Although not obvious from Fig. 15, this collision occurs with the 10  $M_\oplus$ (dark blue line) body orbiting at $\sim$5.5 AU. Second, the 14 $M_\oplus$ protoplanet (orange) undergoes a number of close encounters and orbital exchanges with interior planets, moving inward by ${\simeq} 2$ AU in ${\simeq} 2 \times 10^3$ yr. Third, the 8 $M_\oplus$ (red) body is repetitively scattered outward, and forms a 1:1 resonance with the 12 $M_\oplus$ protoplanet, which quickly evolves from a horseshoe to a tadpole orbit. These configurations are usually highly stable, as shown by the N-body simulations and the hydrodynamic run shown in Fig. 13. In this case, however, the tadpole planets are disrupted by a close encounter which returns the two to horseshoe orbits and eventually leads to the dissolution of co-orbital motion.

Once the initial period of orbital readjustment and collisions has come to an apparent end, the system reverts to the usual process of resonant inward migration. In this case integrations end with three separate groups on diverging orbits: one consists of the outer five protoplanets and another comprises the inner two, with the lone 14 $M_\oplus$ body between them. The remaining resonances are again mostly 5:4 and 6:5 in nature.

Figure 15: As Fig. 12, but with reduced initial separations $\Delta _0=4$ $R_{\rm {mH}}$, similar to Fig. 7. The excitation and subsequent collision of the smallest body is repeated from some N-body runs (cf. Fig. 4), and a co-orbital system forms for ${>}1.7\times 10^4$ yr. Long-term stable resonances are predominantly 5:4 and 6:5. More massive bodies progress inwards by displacing smaller protoplanets to larger radii.

Open with DEXTER

6.4 A hydrodynamic run with smaller separations and larger eccentricities

We performed a hydrodynamic run with initial mutual separations between the protoplanets of 4  $R_{\rm {mH}}$, and initial eccentricities Gaussian distributed with mean $\mu _{\rm e}=0.1$ and standard deviation $\sigma_{\rm e}=0.01$. The evolution of semi-major axes and eccentricities are shown in Fig. 16. The evolution once again follows a now familiar path of initial close encounters, gravitational scattering, and collisional mergers. Of note is the 12  $M_\oplus$ (green line) protoplanet, which exchanges orbits with the 16  $M_\oplus$ (purple) body whilst it is undergoing horsehsoe librations with the 14  $M_\oplus$ (orange) protoplanet, resulting in the 12 and 16 $M_\oplus$ bodies undergoing horseshoe orbits instead. This association is eventually broken by forcing from the approaching 18 $M_\oplus$ body.

In this particular example this early phase of evolution has not reached completion before the innermost planet reaches the inner boundary of the computational domain. The integration was continued for a short time to allow the outer seven bodies to relax, indicating a possible evolution and sequence of resonances for the population of protoplanets, and halted shortly afterwards.

	Figure 16: As Fig. 15, but with the increased initial eccentricities of Fig. 6 ( $\mu _{\rm e}=0.1$). Compare this and Fig. 15 with Figs. 12 and 13; in both pairs of models, the run with higher intial eccentricities does little to raise the number of encounters beyond the first 3000 years.
Open with DEXTER

Figure 17: As Fig. 15, but with with a randomised, Gaussian protoplanet mass distribution, similar to Fig. 8. Working radially outwards through the initial distribution, to the nearest 0.1 $M_\oplus$ the bodies have masses 9.5, 14.1, 8.8, 5.4, 4.9, 11.2, 2.4, 5.9, 18.0 and 10.7 $M_\oplus$. As in the N-body case, the non-sequential arrangement of masses causes the population to split into several resonant groups, while scattering only occurs between bodies with mass ratios greater than ${\sim } 2.5$.

Open with DEXTER

6.5 A hydrodynamic run with a randomised mass distribution

We performed a single run for which the initial mass of each protoplanet was drawn from a Gaussian random distribution, with mean $\mu_{\rm m}=10$ $M_\oplus$ and standard deviation $\sigma_{\rm m}=6$ $M_\oplus$. The initial mutual separation was 5 $R_{\rm {mH}}$ and the initial eccentricities were Gaussian distributed with mean $\mu _{\rm e}=0.05$ and standard deviation $\sigma_{\rm e}=0.01$. The evolution of semi-major axes and eccentricities are shown in Fig. 17.

The results presented in Fig. 17 are in basic agreement with the equivalent N-body runs that we performed. The effect of randomising the protoplanet masses is to increase slightly the degree of scattering and strong interactions. This happens in part due to the more rapid convergence of orbits when two objects of significantly differing mass and migration rate are formed adjacent to one another, with the exterior object being the heavier one. This can lead to the formation of a closer resonance, or allow the planets to pass through successive resonances without being trapped, leading to gravitational scattering. Another outcome of randomising masses is that the global migration of the protoplanetary swarm, which occurs once the initial period of orbital readjustment and collisions is complete, tend to occur in distinctly migrating groups which form separate resonant populations within the swarm.

The existence of three separate resonantly migrating groups of protoplanets can be observed in Fig. 17. The two interior most planets form one group, with the exterior planet of the two being 14  $M_\oplus$. (These bodies quickly move into the damping region, but their subsequent evolution does not affect the other protoplanets due to the large separation between them, and the integration was allowed to continue. This pair are not plotted after the inner body enters the damping zone.) The next group consists of six bodies, two of which are in a 1:1 horseshoe resonance. This group can be observed to be migrating more slowly than the interior resonant pair of planets. After passing through the 6:5 and 7:6 MMRs, the outer two planets form a stable 8:7 resonance at ${\sim} 2 \times 10^4$ yr. Their moderate masses (see text accompanying Fig. 17) mean they will fall behind the other bodies as resonant migration drives the large central group rapidly inwards.

	Figure 18: A series of surface density plots of the simulation shown in Fig. 12. The time in years of each snapshot is shown above each plot. The protoplanets are represented by white circles, with size proportional to mass.
Open with DEXTER

The 1:1 resonance that forms in this case is an interesting example as it takes longer to evolve from horsehoe to tadpole orbits than co-orbital planets do in other simulations. We note that this pair of planets are the heaviest (18.0  $M_\oplus$, cyan) and the lightest (2.4  $M_\oplus$, orange) in the simulation, and as such have the greatest difference in their migration rates. The strong differential migration may explain the extended period of time spent in the horsehoe configuration, as the transformation to tadpole orbits occurs when the migration rate of the pair is slowed as they join the inner resonant group.

We remark that, given the chaotic nature of these systems, we would not expect to see quantitative agreement between the hydrodynamic and modified N-body simulations. Nonetheless, qualitative agreement between the different models is obtained when we consider the simulations as a whole. Common to all models is the strong propensity for the protoplanet swarms to settle into a state of resonant inward migration. Only very rarely ( ${\sim} 2 \%$) do we see a protoplanet get ejected to a much larger semi-major axis such that its lifetime in the disc is substantially increased. This has not been observed in any of the hydrodynamic runs, which may be a consequence of not being able to run for long enough or because its occurrence is too rare. During the first few $\times 10^3$ yr, most of the simulations show that the initial orbital configuration readjusts, leading to gravitational scattering, orbital exchange, collisional mergers, and formation of 1:1 resonances. Very quickly, however, the system settles into a prolonged period of inward migration. Strong eccentricity damping and the existence of numerous mean motion resonances prevent further close encounters and gravitational scattering from occurring, and ensures that the probable end result is migration of the protoplanets into the central protostar.

7 Conclusions

We have presented the results of simulations, performed using two different numerical schemes, that examine the evolution of swarms of low mass protoplanets embedded in a protoplanetary disc. We employed a modified N-body scheme that included analytic prescriptions for modelling migration torques and eccentricity damping due to disc interaction, and calculated the evolution of the systems for times exceeding 10⁵ yr. We also performed hydrodynamic simulations of similar systems to verify the modified N-body results. These runs were unable to follow the evolution for times exceeding $4 \times 10^4$ yr, due to the high computational cost of running long-term hydrodynamic simulations. The results produced by the two simulation methods, however, were found to be very similar during the earlier phases of the evolution, suggesting that the long-term behaviour predicted by the modified N-body code is reliable.

A primary motivation for performing this study was to examine whether interactions between low mass protoplanets embedded in a protoplanetary disc could lead to the growth and maintenance of significant eccentricies (i.e. $e_p \ge 0.1$). It is known that type I migration, associated with low mass protoplanets, can be slowed or even reversed when eccentric orbits are maintained (Papaloizou & Larwood 2000), and an outstanding question was whether gravitational interactions within a protoplanet swarm could allow at least some planets to survive in the disc for extended time periods. The results of our simulations indicate that this is not possible, chiefly because the disc-induced eccentricity damping is too strong (by more than a factor of ten).

The simulations gave rise to a number of dynamical phenomena that occurred commonly in the runs, which followed similar evolutionary paths overall despite variations in the initial conditions. A typical evolutionary sequence can be described as follows: the protoplanets begin to migrate inward differentially due to their different masses and migration rates. Convergence of the orbits leads to mean motion resonances being established between pairs of planets, with resonances involving numerous bodies being common. Close encounters between a few bodies during this early stage leads to gravitational scattering, orbital exchange between adjacent bodies, and occasional collisional mergers. Pairs of planets were often found to form 1:1 resonances, beginning with horseshoe orbits which quickly evolved into stable tadpole orbits, that usually survived for the duration of the simulation. Once this period of initial readjustment is over, typically within $5 \times 10^4$ yr after the beginning of the simulation, the system settles into a state in which the protoplanets are trapped in mean motion resonances involving either the whole swarm of bodies, or groups of protoplanets within the swarm. These resonances were found to typically be 4:3, 5:4, 6:5 and 7:6. Those bodies that form resonant groups migrate inward in lockstep, eventually being accreted by the central protostar in the absence of a stopping mechanism such as a magnetospheric cavity.

There are some exceptions to this basic picture. Occasionally (in ${\sim} 2 \%$ of cases) a lone protoplanet will be ejected to large radius in the protoplanetary disc, causing its lifetime within the disc to be lengthened substantially. There are simulations which show short bursts of scattering and strong interactions between the protoplanets as they approach the protostar within ${\sim} 1$ AU, but these systems quickly settle down into a state of orderly inward migration again. Overwhelmingly, however, the systems show large scale migration of the swarm into the star. The reasons for this are basically threefold:

(i)

Eccentricity damping by the disc is simply too strong for substantial eccentricities to be generated and maintained by either close encounters or the majority of eccentricity-driving mean motion resonances.

(ii)

The formation of mean motion resonances, and the resulting tendency for protoplanets to migrate inward in lockstep, means that close encounters between bodies occur rarely.

(iii)

Inelastic collisions between bodies provide a source of energy dissipation that reduces the tendency for strong scattering to occur.

These results show that if multiple protoplanets form approximately coevally in the giant planet zone between 5-20 AU within a laminar disc with simple structure due to oligarchic growth (e.g. Kokubo & Ida 2000; Thommes et al. 2003), then the long-term evolution of that system will be collective inward migration and eventual accretion by the protostar. This will occur on a time scale much shorter than that required for the cores of putative giant planets to accrete their gaseous envelopes (e.g. Pollack et al. 1996; Papaloizou & Nelson 2005). If protoplanetary discs do indeed have simple density and temperature structure, and are non turbulent near their midplane due to the low ionisation fraction there (e.g. Gammie 1996; Fromang et al. 2002; Ilgner & Nelson 2005) then the question of how giant planets form via the core accretion model remains. Recent work examining the interaction of low mass protoplanets with fully turbulent discs shows that a process of stochastic migration ensues, which may help to overcome type I migration (Nelson & Papaloizou 2004; Nelson 2005). It remains a possibility that the combination of turbulent fluctuations and planet-planet interactions may be even more effective in overcoming type I migration.

The frequency with which co-orbital planets in horseshoe and tadpole orbits form in the simulations suggests that such configurations may be observed, if they are able to survive against type I migration. Whether such bodies could accrete massive gaseous envelopes and survive as co-orbital giants is not known. We note that Laughlin & Chambers (2002) have suggested alternative scenarios for forming co-orbital planets. One involves formation in situ. The other involves direct physical collision, leading to a binary planet that subsequently separates and forms a trojan pair. We find co-orbitals form naturally from existing bodies on initially distinct orbits due to close encounters during chaotic phases of evolution.

The models that we have presented here can be significantly improved. Most important is the inclusion of the 3rd dimension in the simulations, as this will reduce the collision rate between protoplanets, and may thus prolong the time during which the protoplanets interact strongly. Preliminary 3D simulations, however, indicate that the picture presented in this paper does not alter radically. These three dimensional simulations will be the subject of a forthcoming paper.

Acknowledgements

The simulations reported here were performed using the UK Astrophysical Fluids Facility (UKAFF) and the QMUL HPC facility purchased under the SRIF initiative. P.C. is supported by a PPARC Ph.D. studentship.

References

• Artymowicz, P. 1993a, ApJ, 419, 166 [NASA ADS] [CrossRef] (In the text)
• Artymowicz, P. 1993b, ApJ, 419, 155 [NASA ADS] [CrossRef] (In the text)
• Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 [NASA ADS] [CrossRef] (In the text)
• Bonfils, et al. 2005, A&A, 443, L15 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Butler, R. P., Marcy, G. W., Vogt, S. S., et al. 2003, ApJ, 582, 455 [NASA ADS] [CrossRef] (In the text)
• Butler, R. P., Vogt, S. S., Marcy, G. W., et al. 2004, ApJ, 617, 580 [NASA ADS] [CrossRef] (In the text)
• Chambers, J. E., & Wetherill, G. W. 1998, Icarus, 136, 304 [NASA ADS] [CrossRef] (In the text)
• Chambers, J. E., Wetherill, G. W., & Boss, A. 1996, Icarus, 313, 261 [NASA ADS] [CrossRef] (In the text)
• Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18 [NASA ADS] [CrossRef]
• Gammie, C. F. 1996, ApJ, 457, 355 [NASA ADS] [CrossRef] (In the text)
• Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [NASA ADS] [CrossRef] (In the text)
• Gladman, B. 1993, Icarus, 313, 247 [NASA ADS] [CrossRef] (In the text)
• Ilgner, I., & Nelson, R. P. 2006, A&A, 445, 205 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kley, W. 2000, MNRAS, 313, 47 [NASA ADS] [CrossRef] (In the text)
• Kley, W., Lee, M. H., Murray, N., & Peale, S. J. 2005, A&A, 437, 727 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Kokubo, E., & Ida, S. 2000, Icarus, 143, 15 [NASA ADS] [CrossRef] (In the text)
• Kominami, J., Tanaka, H., & Ida, S. 2005, Icarus, 178, 540 [NASA ADS] [CrossRef] (In the text)
• Korycansky, D. G., & Pollack, J. B. 1993, Icarus, 102, 150 [NASA ADS] [CrossRef] (In the text)
• Laughlin, G., & Chambers, J. E. 2002, AJ, 124, 592 [NASA ADS] [CrossRef] (In the text)
• Lee, M. H., & Peale, S. J. 2001, BAAS, 33, 1198 [NASA ADS] (In the text)
• Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606 [NASA ADS] [CrossRef] (In the text)
• Marcy, G. W., Butler, R. P., Fischer, D., et al. 2001, ApJ, 556, 296 [NASA ADS] [CrossRef] (In the text)
• Marcy, G. W., Butler, R. P., Vogt, S. S., et al. 2005, ApJ, 619, 570 [NASA ADS] [CrossRef] (In the text)
• Masset, F., & Snellgrove, M. 2001, MNRAS, 320, 55 [NASA ADS] [CrossRef] (In the text)
• Nelson, R. P. 2005, A&A, 443, 1067 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Nelson, R. P., & Papaloizou, J. C. B. 2002, MNRAS, 333, 26 [NASA ADS] [CrossRef] (In the text)
• Nelson, R. P., & Papaloizou, J. C. B. 2004, MNRAS, 350, 849 [NASA ADS] [CrossRef] (In the text)
• Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 [NASA ADS] [CrossRef] (In the text)
• Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823 [NASA ADS] [CrossRef] (In the text)
• Papaloizou, J. C. B., & Nelson, R. P. 2005, A&A, 433, 247 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Papaloizou, J. C. B., & Szuszkiewicz, E. 2005, MNRAS, 363, 153 [NASA ADS] [CrossRef] (In the text)
• Pollack, J. B., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] (In the text)
• Press, W., et al. 1992, Numerical Recipes in FORTRAN (Cambridge: Cambridge Univ. Press), 710 (In the text)
• Rivera, E. J., Lissauer, J. J., Butler, R. P., et al. 2005, ApJ, 634, 625 [NASA ADS] [CrossRef] (In the text)
• Santos, N. C., Bouchy, F., Mayor, M., et al. 2004, A&A, 426, L19 [EDP Sciences] [NASA ADS] [CrossRef] (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 [EDP Sciences] [NASA ADS] (In the text)
• Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 75 [NASA ADS] (In the text)
• Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388 [NASA ADS] [CrossRef] (In the text)
• Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 [NASA ADS] [CrossRef] (In the text)
• Thommes, E. W. 2005, ApJ, 626, 1033 [NASA ADS] [CrossRef] (In the text)
• Thommes, E. W., Duncan, M. J., & Levison, H. F. 2003, Icarus, 161, 431 [NASA ADS] [CrossRef] (In the text)
• Van Leer, B. 1977, Comput. Phys., 23, 276 [NASA ADS] (In the text)
• Vogt, S. S., Butler, R. P., Marcy, G. W., et al. 2005, ApJ, 632, 638 [NASA ADS] [CrossRef] (In the text)
• Ward, W. R. 1986, Icarus, 67, 164 [NASA ADS] [CrossRef] (In the text)
• Ward, W. R. 1997, Icarus, 126, 261 [NASA ADS] [CrossRef] (In the text)
• Weidenschilling, S. J. 2005, Space Sci. Rev., 116, 53 [NASA ADS] [CrossRef] (In the text)
• Ziegler, U., & Yorke, H. 1997, Comput. Phys. Commun., 101, 54 [NASA ADS] (In the text)