Abstract

We examine the effect of a rapidly migrating protoplanet on a ring of planetesimals. The eccentricities of the planetesimals are usually increased by Δe∈ (0.01, 0.1), with the exact increase being proportional to the mass of the protoplanet, and inversely proportional to its migration rate. The eccentricity distribution is also substantially changed from a Rayleigh distribution. We discuss the possible implications for further planet formation, and suggest that the rapid passage of a protoplanet may not prevent the planetesimal disc from forming further planets.

1 INTRODUCTION

A planet in orbit around a star will disturb the orbits of test particles. This situation, known as the ‘restricted three body problem’ has been studied for centuries (see, e.g., Murray & Dermott 1999). Particles close to resonant locations have their orbital elements changed due to repeated interactions with the planet. These changes may be computed using various analytic techniques (ibid), and the analysis can be extended to slowly migrating planets. Recently, a rapid migration mode for protoplanets has been discovered (Masset & Papaloizou 2003). The planet can have its semimajor axis halved in only a hundred orbits or so, which is far too fast for the effect on planetesimal orbits to be calculated analytically.

In this paper, we attempt to quantify the expected eccentricity increase in a planetesimal ring, due to a rapidly migrating protoplanet. A simple theoretical calculation is described in Section 2. We describe a numerical model in Section 3 and the results obtained from it in Section 4. The possible implications for the formation of the Solar system are discussed in Section 5. We summarize our findings in Section 6.

2 THEORETICAL PREDICTIONS OF PLANETESIMAL EXCITATION

The orbits of the planetesimals can be changed in two ways:

  • (i)

    long term interactions with resonances (mean motion, secular and corotation);

  • (ii)

    close encounters with the planet.

Interactions between test particles and resonances have been studied for several centuries. However, for a rapidly migrating planet, these resonances are much less important, as they will be sweeping through the disc with the migrating planet. Hence, the resonance is unlikely to remain in the vicinity of a particle long enough to make significant changes. This is supported by the work of Tanaka & Ida (1999), who found that a slow migration rate led to ‘shepherding,’ (the planet gently brings the planetesimals with it), while fast migration leads to ‘predatory’ behaviour (the planet simply ploughs through the particle disc) – see also the work of Ward & Hahn (1995). In our calculations, migration is even faster than the ‘rapid’ migration of Tanaka & Ida (1999).

Hasegawa & Nakazawa (1990) found that the expected change in eccentricity of a single particle per collision was given by
1
where R1= 0.747, h= (mplan/3m*)1/3 is the Hill parameter and b is the impact parameter of the collision, in units of the semimajor axis. This is equation (38) of their paper, rewritten with conventional (rather than normalized) eccentricities. It is derived on the assumption that the stellar mass (m*) dominates, that the eccentricities and inclinations are small, and that the impact parameter is large enough to ensure that the particles cannot enter their mutual Hill sphere. Hasegawa & Nakazawa give another, similar, formula for the expected increase in inclination, but the coefficient is much smaller and we shall neglect inclination here (note also that the eccentricity increase is independent of the inclination). Equation (1) (with a different coefficient) may also be derived following the impulse approximation of Lin & Papaloizou (1979). Note the strong dependence on the impact parameter, curiously reminiscent of Rutherford scattering, although the link is not direct. This suggests that it is the closest encounter between the planet and the planetesimal that will be the most important.
However, what is the closest encounter that we should expect? Ever closer encounters will give larger changes (although note that equation 1 eventually ceases to be valid), but encounters with small impact parameters are rare. Both the planet and the planetesimals are on near-circular orbits, so an impact parameter of zero implies that the two bodies are on the same orbit. Such particles would have a (near) infinite synodic period. Consider a particle at semimajor axis a0, and the planet migrating from a0a to a0−Δa. We can expect an interaction with b≈Δa/a0 if the time for the planet to migrate between the two limits graphic is equal to the synodic period for orbits at a0 and a0a. Formally, this neglects the orbital inclination, but we expect the effect to be small in a disc. We have
2
3
4
where vkep is the circular Keplerian velocity, and Tkep is the Keplerian orbital period. Substituting into equation (1), we find
5

For a 50-M⊕ planet in orbit around the Sun and migrating at 104 cm s−1 through a ring close to 3.5 au,equation (5) suggests 〈Δe21/2≈ 0.03. Several approximations were made in deriving equation (5). It assumes a single close encounter, and the estimated time-scale for this (based on the synodic period) is rather imprecise. The numerical factor could easily be incorrect. At this point, numerical simulations become useful in determining the evolution.

3 MODEL

To simplify the system as much as possible, we considered the interaction of a migrating planet with a narrow ring of test particles.

The bodies in our simulation belonged to one of three types:

  • (i)

    the central star, m*;

  • (ii)

    the migrating planet, mplan;

  • (iii)

    the test particles (planetesimals), mpart.

We evolved our system using Newtonian gravity, plus a torque to migrate the planet. To avoid the computational expense of a full n-body calculation, the test particles did not interact with each other. However, the gravitational interactions between the star and the planet, the star and the particles, and the planet and the particles, were all fully computed.

We neglected the effect of gas damping, as we were interested in the effects of rapid migration. This is reasonable, as our largest migration time-scales were still far shorter than the damping time-scales. Using the notation of Tanaka & Ida (1999), we typically had graphic, while graphic, even with enhanced gas density. This holds even if the gas damping time-scale is assumed to be due to wave excitation, as will be the case for the more massive planetesimals (Artymowicz 1993). Hence we did not expect the gas damping to be sufficient to make our migrating protoplanets ‘shepherd’ the planetesimals.

3.1 Orbital migration

Orbital migration of the planet is imposed by applying a torque to the orbit of the star and the planet. In the interest of simplicity, we apply a torque chosen to give a constant migration rate, graphic. The required torque may be computed using Kepler's laws:
6
where La2Ω is the orbital angular momentum and G is the required torque. The torque was applied as an extra force to the motion of the star and the planet. Migration was halted once the semimajor axis of the planet dropped below a preset value.

3.2 Initial conditions and integration

The star (1 M) and the planet were positioned in a circular orbit about their centre of mass. The initial orbital separation was usually 6 au and migration stopped at 0.5 au. We distributed the test particles in circular orbits around the centre of mass of the star–planet system. The particles were uniformly distributed in azimuth and radius (usually 3–4 au). Finally, we added small perturbations to their circular motion, following the prescription of Stewart & Ida (2000), to give a Rayleigh distribution of e and i (ibid). We tested a variety of planetary masses and migration rates.

We integrated the equations of motion using a Runge–Kutta integrator with adaptive step-size. To avoid catastrophically small time-steps, test particles that strayed too close to the star or the planet were eliminated. Similarly, particles that reached large radii were removed. At the end of each time-step, we computed the values of e, i and a for each particle using the Laplace–Runge–Lenz formalism.

4 RESULTS

4.1 General behaviour

When the planet was far from the planetesimal ring, waves in (e, a) space were observed propagating through the ring. These were caused by resonances sweeping through the ring – if the planet did not migrate, then sharp peaks in the eccentricity were observed, corresponding to resonant orbits. However, the increase in eccentricity prior to the encounter of the planet with the ring was fairly low. As the planet ploughed through the ring, the planetesimals were catapulted up lines of constant Jacobi energy, EJ, given by (Hayashi, Nakazawa & Adachi 1977):
7
where a′ is the particle semimajor axis in units of the semimajor axis of the planet (this energy is expressed in scaled units). The encounters with the planet were generally at a distance further than the L2 point (this is consistent with the assumptions made in deriving equation 1).

4.2 Numerical results

Several runs were performed, with the parameters and results summarized in Table 1. For all these runs, the ring of particles lay between 3 and 4 au initially. In this table ξ2=e2+i2, the conventional measure of deviation from circular, coplanar orbits. In practice, ξ was dominated by eccentricity, not inclination (cf. Section 2). We used a Kolmogorov–Smirnov test to investigate whether or not the final distribution of ξ was still drawn from a Rayleigh distribution. In every case, this hypothesis was strongly rejected. Qualitatively, we found that there was a tail of high-ξ particles. This made the rms values for ξ somewhat changeable, so we quote the median value as well. For the remainder of this discussion, references to the final ξ values refer to the median.

1

Numerical results. The quantity ξ is defined as ξ2=e2+i2. Recall that 104 cm s−1≈ 0.02 au yr−1.

Most of the results of Table 1 bear out the qualitative behaviour suggested by equation (5), but differ quantitatively. The final median ξ values are close to those predicted (recall that ξ was dominated by e in our simulations). Increasing the migration time-scale or the mass of the migrating planet increases the final value of ξ. However, this does not quite happen in the linear manner predicted by equation (5). The detailed values are similar to those predicted, but not identical.

Fig. 1 shows part of the problem, plotting two sample particle trajectories in ea space. These particles have moved around significantly in ea space, and have not done this just once. The trajectories are characterized by large jumps (corresponding to those analysed above) interspersed with periods where the particles are almost stationary. Sometimes, the eccentricity is pumped up more slowly as well (presumably when the particle happened to be close to one of the resonances of the migrating planet). This is in sharp contrast to the ‘single encounter’ approximation of Section 2. Obtaining better predictions of the final ξ value is therefore problematic. Even if the resonant pumping is neglected, we cannot apply a simple random walk approach. First, the number of steps (large ‘jumps’) is also random – and not large. Worse, the steps are not truly random, but are constrained to have constant Jacobi energies (although the constant value changes for each step, as the planet will have migrated in the meantime).

1

Two sample particle trajectories in ea space.

Run 8 behaved slightly differently. However, even this is not surprising when compared with run 7. The two runs were identical apart from the initial ξ value, and the initial ξ for run 8 was higher than the final value for run 7. It is therefore not surprising that the final ξ value for run 8 is higher than that for run 7 – but note that the increase in ξ is similar for both runs.

We also ran a grid of 40 models, with mplan ranging between 10 and 400M⊕ and graphic in the range from 2 × 103 to 104 cm s−1. This covers the range expected by Masset & Papaloizou (2003) to undergo runaway migration, and a bit more on each end. All these runs started with average e and i values of 10−3. Fig. 2 plots a surface showing the resultant median ξ values. Planets more massive than 100M⊕ frequently managed to increase ξ to be greater than 0.1 (and reached 0.4 for a 400-M⊕ planet migrating at 2 × 103 cm s−1). Although not a perfect match, the behaviour predicted by equation (5) is seen.

2

Surface of median ξ values following migration.

5 IMPLICATIONS

One major problem with the migration of planets in their nascent discs is that protoplanets ‘could be too mobile for their own good’ (Ward & Hahn 2000). The rapid migration mode of Masset & Papaloizou (2003) makes the protoplanets even more mobile – it seems to be very easy to lose protoplanets into the Sun. Whilst this is obviously not very ‘good’ for the protoplanets in question, would such events make the formation of the Solar system impossible? Put another way, ‘Is the current Solar system the first Solar system?’

Considering the Q parameter of Toomre (1964) suggests that the disc around the young Sun may have been up to five times more massive than the minimum mass solar nebula (MMSN) model of Hayashi (1981), although Hayashi, Nakazawa & Nakagawa (1985) is likely to be a more accessible reference. This would give plenty of material for forming several generations of planets. The question is therefore whether a migrating planet will disrupt the disc sufficiently to prevent further planet formation within the nebula lifetime of a few Myr.

If the migration rate is low enough to permit shepherding (and hence depletion of planetesimals inward of the initial location of the protoplanet), Armitage (2003) has shown that it is unlikely that fresh material could diffuse into the inner portions of the disc. In this case, further planet formation would not be possible. We have considered a migration rate too high for shepherding, so we must examine whether our planetesimal disc is too hot to allow further planet formation.

In our simulations, Δe is typically in the range (0.01, 0.1), dependent on the mass of the protoplanet and its migration rate. Planetesimals are fairly fragile objects – bodies with mpart≈ 1022 g will suffer disruptive collisions if e≳ 0.01 (Kobayashi & Ida 2001). However, gas damping can recircularize the orbits of such bodies very rapidly – probably less than 104 yr [see fig. 1 of Tanaka & Ida (1999) and references therein]– especially if the disc is denser than the MMSN. The damping time-scale will reach a peak for mpart≈ 1025 g, but these bodies require e≳ 0.1 before suffering disruptive collisions. Still larger bodies will require even higher e values before they shatter, but these damp faster due to the excitation of density waves (Artymowicz 1993). Gravitational focusing will also fall as eccentricities increase, due to the higher relative velocities implied. This will cause a dramatic fall in the collision rate (cf. fig. A1 of Weidenschilling et al. 1997), so the likelihood of disruptive collisions is reduced. Finally, in a real planetesimal disc there will be a distribution of sizes. Wetherill & Stewart (1993) found that grinding moderately sized planetesimals into rubble helps the larger planetesimals accrete them, so some excitation may even be helpful.

Based on this discussion, it seems that the rapid passage of a protoplanet will not prevent further planet formation. Coagulation will be suppressed for a while, but this should be short compared with the disc lifetime. However, the arguments in the preceding paragraph are not rigorous, and a longer term calculation of the ‘end’ states would be required to give a firm answer. This is particularly true just as the planetesimal disc re-achieves equilibrium.

6 CONCLUSIONS

Combining our results (including some varying a0 not explicitly listed here), we conclude that a rapidly migrating protoplanet will increase the eccentricity of a planetesimal ring by roughly
8
where a0 is the initial semimajor axis of the planetesimal ring and graphic is the migration rate of the protoplanet. During the passage through the ring, the distribution of eccentricities will also be substantially changed from a Rayleigh distribution. In particular, there will be a tail of high eccentricity bodies.

Simple arguments suggest that this effect should not prevent subsequent planet formation by the surviving planetesimals (very few are ejected or accreted). The resultant eccentricities should not be dangerously high for at least some of the planetesimals, collisions are less likely due to reduced gravitational focusing, and damping time-scales are short. However, problems could arise shortly before the planetesimal disc achieves equilibrium once more, and further simulations (spanning longer time-scales and including all gravitational interactions) are needed to assess the effect in detail.

Acknowledgments

RGE acknowledges financial support provided through the European Community's Human Potential Programme under contract HPRN-CT-2002-00308, PLANETS. Some of the calculations used the resources of the High Performance Computing Centre North (HPC2N) and the National Supercomputing Centre, Linköping. The KS test was performed using a publically available routine written by W. Landsman. The authors would like to acknowledge helpful discussions with Willi Kley, which took place while attending the KITP workshop on planet formation.

References

Armitage
P. J.
,
2003
,
ApJ
,
582
,
L47
DOI:

Artymowicz
P.
,
1993
,
ApJ
,
419
,
166
DOI:

Hasegawa
M.
Nakazawa
K.
,
1990
,
A&A
,
227
,
619

Hayashi
C.
,
1981
,
Prog. Theor. Phys. Suppl
.,
70
,
35

Hayashi
C.
Nakazawa
K.
Adachi
I.
,
1977
,
PASJ
,
29
,
163

Hayashi
C.
Nakazawa
K.
Nakagawa
Y.
,
1985
, in
Protostars and Planets II
.
Univ. Arizona Press
Oxford, UK
, Tucson , pp.
1100
1153

Kobayashi
H.
Ida
S.
,
2001
,
Icarus
,
153
,
416
DOI:

Lin
D. N. C.
Papaloizou
J.
,
1979
,
MNRAS
,
186
,
799

Masset
F. S.
Papaloizou
J. C. B.
,
2003
,
ApJ
,
588
,
494
DOI:

Murray
C. D.
Dermott
S. F.
,
1999
,
Solar System Dynamics
.
Cambridge Univ. Press
Oxford, UK
, Cambridge

Stewart
G. R.
Ida
S.
,
2000
,
Icarus
,
143
,
28
DOI:

Tanaka
H.
Ida
S.
,
1999
,
Icarus
,
139
,
350
DOI:

Toomre
A.
,
1964
,
ApJ
,
139
,
1217
DOI:

Ward
W. R.
Hahn
J. M.
,
1995
,
ApJ
,
440
,
L25
DOI:

Ward
W. R.
Hahn
J. M.
,
2000
,
Protostars and Planets IV
.
Univ. Arizona Press
Oxford, UK
, Tucson , p.
1135

Weidenschilling
S. J.
Spaute
D.
Davis
D. R.
Marzari
F.
Ohtsuki
K.
,
1997
,
Icarus
,
128
,
429
DOI:

Wetherill
G. W.
Stewart
G. R.
,
1993
,
Icarus
,
106
,
190

Footnotes

1

Note that 104 cm s−1≈ 0.2 au yr−1, which is not unreasonable for the rapid migration discussed by Masset & Papaloizou (2003).

2

For the short time-scales we are studying this is reasonable, as the close proximity of the (relatively) massive planet will have a far greater effect.

3

We also checked that the KS test did allow the initial distribution to be Rayleigh.