A&A 450, 833-853 (2006)
DOI: 10.1051/0004-6361:20054551
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
,
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
of runs a single planet is scattered onto a distant orbit with significant eccentricity, allowing it to survive in the disc for
years. Over
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
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 100-1000 km-sized planetary embryos; (iii) oligarchic growth of these embryos, which accrete the surrounding planetesimals, eventually forming planetary cores of 10 in the giant planet region beyond 3 AU, or Mars-mass bodies in the terrestrial planet zone. The 10-15 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^{5} years for bodies in the 10 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 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 , 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.
A protoplanetary disc undergoing near-Keplerian rotation has
a ratio of scale height to radius
.
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
,
with the origin located at the central star.
The continuity equation takes the form
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:
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.
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
,
and adopt a locally isothermal equation of state such that the
pressure, P, and surface density,
are related by
The disc models are viscous, and we adopt the "alpha'' prescription for the kinematic viscosity (Shakura & Sunyaev 1973). We set .
The motion of protoplanets within the disc is integrated using a -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 .
When computing the time step size at each time level we first
calculate a "hydrodynamic time step'',
,
and an "N-body time step'',
.
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:
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
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 tend to pile up near the planet, standing-off at a distance of 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:
Figure 1: The inverse of migration rates vs. eccentricity for four planet masses. Top left panel - 1 , top right panel - 5 , lower left panel - 10 , lower right panel -15 . 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 , the inverse of the migration time scale, rather than in the following figures.
The code calibration simulations described here adopted a disc model with , with 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 , unless stated otherwise. The planet orbital evolution due to disc forces was calculated explicitly.
Figure 2: Eccentricity damping rates 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 e_{p}> 0.07. | |
Open with DEXTER |
Figure 1 shows the inverse of the migration time 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 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 at low e, and by a factor of 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 , 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 . Here we find that Eq. (13) predicts faster eccentricity damping than the hydrodynamic simulations, being about a factor of three faster for . 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.
We use softening when calculating the gravitational interaction between disc and protoplanets, as shown in Eqs. (5) and (7). Higher values of 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 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 , the planet Hill radius, and from , and . 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 5 on a fixed circular orbit. Four torque cutoffs are represented: (solid), (dotted), 1.0 (dashed) & (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 arises when , 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
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
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 .
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 order Runge-Kutta routine as used in NIRVANA. To best compare the two schemes we restrict the N-body code to the 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.
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:
We implement the following expressions in the N-body code as
accelerations experienced by the protoplanets due to the disc, using
values of
and
obtained from Eqs. (16)
and (17):
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_{1}=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,
where
usually took the value 4 or 5 (but other values were considered), and the mutual Hill radius is defined by
When setting the protoplanetary masses, our normal procedure was to define the mass of the innermost body (usually m_{1}=2 ), with subsequent bodies having m_{i+1} = m_{i}+2 . 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, , and a standard deviation , 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 and 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 , where was defined such that the disc contains 40 Jupiter masses within 40 AU. In NIRVANA simulations the disc had a viscous parameter of .
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.
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 , with the planet mass increasing by 2 as one moves out through the initial swarm, such that the outermost planet was 20 . This maximises convergent migration. The initial eccentricities were such that the mean value was and standard deviation . The variation of semi-major axes and eccentricities for this model are shown in Fig. 4 for a run time of 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 yr.
In Fig. 4, this scattering leads to a physical collision between the 8 and 10 bodies after yr, causing them to merge and form a new 18 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
yr into the run, the innermost and least massive of
the embryos scatters during an encounter
with the adjacent 4
protoplanet,
exchanging positions with it. Subsequent interactions cause this
2
embryo to be scattered inward due to close encounters with
the newly formed
18
protoplanet and the pre-exisiting 12
body.
It then collides with the 6
protoplanet, forming a new
8
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
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:
Figure 5: First order resonances between five of the bodies in Fig. 4. Each plot is labelled (P_{1}, P_{2})-p:q, where P_{1}/P_{2} 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 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 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 AU. The system, however, quickly returns to an orderly state of inward resonant migration.
For a qualitative change in the evolution to occur, the initial eccentricities needed to be around , 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 ( 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 ( compared to 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 ( as before). Close inspection of Fig. 6 shows that the innermost protoplanet, with mass m_{p}= 2 quickly interacts with its nearest neighbour, leading to them exchanging orbital locations. The 2 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 , respectively. Interaction between the 2 and 18 bodies leads to the 2 object being scattered into an orbit with a_{p}=45 AU and e=0.75. This orbit crosses that of the 20 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 protoplanets experience a collisional merger, resulting in the formation of an 18 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 ( ). 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 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 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.
Figure 7: As Fig. 4, but with initial separations of 4 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 body displacing the inner pair at 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 . 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 such configurations occur in % of the runs. This figure rises to % when is reduced to 4 . 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 , 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 (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 (green) and 16 (purple) protoplanets at 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 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_{4} and L_{5} points. This explains why the purple line in Fig. 7 becomes obscured by the green line.
Another co-orbital resonance forms between two 8 bodies, one of which is formed through collision at 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 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_{4}/L_{5} points at 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 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 , 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.
Separations of 6 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 to 10 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.
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 , with the standard deviation being 3-6 , usually subject to lower and upper mass cut-offs of 2 and 20 . 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 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 . The non-monotonic distribution of migration rates causes the population to fragment into smaller resonant groups. The 2.9 body (orange) is excited rather than captured by the neighbouring 10.6 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 and standard deviation of 6 . The masses of each body are listed in the caption of Fig. 8. The 2.9 (orange line) body is briefly trapped in the 8:7 MMR with the 10.6 (purple) planet, but the resonance breaks at yr. The lighter body is then scattered outward before being caught in a stable MMR with the outermost 6.6 (yellow) body.
We note that the mass ratio of the 2.9 and 10.6 bodies is 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 . For values of 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 . The protoplanets have reduced masses 1 , 1.5 , 2 ,..., 10.5 , 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 (green) protoplanet. The middle group of two planets migrates more slowly and is driven by the 14.7 (cyan) protoplanet. The outer group migrates the slowest and is driven by the 6.6 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.
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 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 yr, redistributing bodies throughout the disc. Without energy losses to the disc, migration is reversed at even moderate eccentricities. | |
Open with DEXTER |
It was found that decreasing the surface density (disc mass) to 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 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 0.2.
A run with 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 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.
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.
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 e_{p}>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 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 yr when the 6 (olive green) and 8 (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 . The resulting 14 body resonantly drives the 4 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 (orange line) and 12 (green) bodies which are each adjacent to the 10 (dark blue) protoplanet with which they are also in resonance. Another 3:2 resonance forms between the outermost 20 (yellow) and 16 (purple) bodies following the disturbance at 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 . 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 yr later. Rapid eccentricity damping is clearly visible early on in the 4 (light blue), 6 (olive green) and 2 (black) bodies; only the former acheives e_{p} > 0.1. | |
Open with DEXTER |
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 (P_{2}=2,4) both share the two resonances with exterior bodies (P_{1}=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 yr. Initially we observe that the 2 (black line) and 4 (light blue) planets scatter and move apart, which ultimately results in the formation of a 1:1 resonance between the scattered 4 protoplanet and the 8 body. The initial horseshoe orbits that form quickly evolve to become tadpole orbits at L_{4} and L_{5}. This system then remains stable over the full run time of the simulation.
A sole collision occurs later at yr between the 2 and 6 (olive green) protoplanets. This is a result of the 6 body being resonantly driven inward by the more rapidly migrating exterior planets after yr, breaking the 5:4 resonance between the 2 and 6 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.
We performed a hydrodynamic run with initial protoplanet separations reduced from 5 to 4 . 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 (dark blue line) body orbiting at 5.5 AU. Second, the 14 protoplanet (orange) undergoes a number of close encounters and orbital exchanges with interior planets, moving inward by AU in yr. Third, the 8 (red) body is repetitively scattered outward, and forms a 1:1 resonance with the 12 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 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 , 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 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 |
We performed a hydrodynamic run with initial mutual separations between the protoplanets of 4 , and initial eccentricities Gaussian distributed with mean and standard deviation . 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 (green line) protoplanet, which exchanges orbits with the 16 (purple) body whilst it is undergoing horsehsoe librations with the 14 (orange) protoplanet, resulting in the 12 and 16 bodies undergoing horseshoe orbits instead. This association is eventually broken by forcing from the approaching 18 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 ( ). 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 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 . 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 . | |
Open with DEXTER |
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 . (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 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 , cyan) and the lightest (2.4 , 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 ( ) 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 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.
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^{5} 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 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. ). 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 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 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 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:
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.