Giant Planet Migration in Viscous Power-Law Discs
Richard .G. Edgar
Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627
rge21@pas.rochester.edu
ABSTRACT
Many extra-solar planets discovered over the past decade are gas giants in
tight orbits around their host stars. Due to the difficulties of forming these ‘hot
Jupiters’ in situ, they are generally assumed to have migrated to their present
orbits through interactions with their nascent discs. In this paper, we present
a systematic study of giant planet migration in power law discs. We find that
the planetary migration rate is proportional to the disc surface density. This is
inconsistent with the assumption that the migration rate is simply the viscous
drift speed of the disc. However, this result can be obtained by balancing the
angular momentum of the planet with the viscous torque in the disc. We have
verified that this result is not affected by adjusting the resolution of the grid, the
smoothing length used, or the time at which the planet is released to migrate.
Subject headings: hydrodynamics – planetary systems: protoplanetary disks –
planetary systems: formation
1. Introduction
The past decade has seen the discovery of over one hundred extra-solar planets. Most of
these planets (see, e.g. Marcy et al. 2005) are gas giants of Jupiter mass or greater, in orbits
close to their host stars.1 Such planets are often called ‘hot Jupiters.’ Although detection
methods (primarily radial velocity measurements) are strongly biased towards detecting such
planets, a sufficient number have been discovered to require detailed explanation.
The present orbits of the ‘hot Jupiters’ lie so close to their host stars that in situ forma-
tion is almost impossible. The problem is two-fold. Firstly, the expected disc temperature
is high, which suppresses gravitational capture of gas by a gas giant core. Secondly, there
1For the most up to date information, refer to http://exoplanet.eu/ and http://exoplanets.org/
http://arxiv.org/abs/0704.0448v1
http://exoplanet.eu/
http://exoplanets.org/
– 2 –
simply isn’t much material at such small orbital radii. We are therefore led to the conclusion
that the ‘hot Jupiters’ formed further from the star, and migrated to their present orbits.
In this paper, we present a study of giant planet migration in power law discs.
A planet in a circumstellar disc exerts a torque on the disc material at resonances
(Goldreich and Tremaine 1978, 1979, 1980). By Newton’s Third Law, there is an opposite
torque on the planet, which causes it to migrate. Planets of mass comparable to Jupiter (the
exact threshold depends on the disc conditions, see below) exert torques comparable to the
viscous torque in the disc. They can therefore open a gap in the disc. If the gap is perfectly
clean, then the planet will act as a ‘relay station’ between the inner and outer discs. We then
expect the planet’s orbital evolution to follow the viscous evolution of the disc – a process
known as Type II migration (Ward and Hahn 2000). We shall compare our computational
results to the predictions of Type II theory.
We summarise the theory of Type II migration in Section 2. We describe our numerical
method in Section 3, and demonstrate the formation of a gap in Section 4. Section 5 presents
our main results. In Section 6 we discuss the implications of our findings, before presenting
our conclusions in Section 7.
2. Type II Migration
In this section, we shall briefly review the theory of Type II migration. For a more
thorough analysis of the theory of planet–disc interactions, the reader should turn to one of
the many published reviews (e.g. Lin and Papaloizou 1993; Ward and Hahn 2000; Lin et al.
2000).
Planet–disc interactions are dominated by resonances. A planet embedded in a cir-
cumstellar disc excites waves at its Lindblad resonances (LR). These waves carry angular
momentum, and hence exert a torque on the disc material. The strength of the torque from
the Lindblad resonances is proportional to TLR ∝ Σq
2, where Σ is the local disc surface
density, and q is the planet-star mass ratio. These torques act to push material away from
the planet. At the same time, the disc gas is expected to have an intrinsic viscosity, ν (al-
though the precise origin and exact behaviour of such a viscosity are still much debated),
which leads to a viscous torque Tν ∝ Σν. Since the Lindblad torques scale with q2, we can
expect that as the planet accretes material, the Lindblad torques will eventually dominate
the viscous torque in the disc. Balancing the two torques leads to the so-called gap opening
criterion:
q > 40R−1 (1)
– 3 –
where R = r2Ω/ν is the Reynolds number of the disc (Bryden et al. 1999; Nelson et al.
2000). There is a similar requirement that the planet’s Hill sphere exceed the local scale
height of the disc, namely that
q > 3
where h is the disc scale height. See Lin and Papaloizou (1993) for a further discussion.
The tidal condition of Equation 2 leads to the gap width being at least twice the local scale
height. If this condition is not satisfied, then the edge of the gap will be Rayleigh unstable.
Consideration of the torque condition leads to the expectation that the gap will lie between
the m = 2 Lindblad resonances of the planet (that is, between r = 0.6 and r = 1.3). For
Jovian mass planets in circumstellar discs, these gap sizes are similar.
Once the planet has opened a gap, it is assumed to isolate the inner and outer discs
from each other. However, each part of the disc will still be undergoing viscous evolution.
According to Type II migration theory, the planet will be locked to the viscous evolution of
the disc (Ward and Hahn 2000). Setting the migration rate of the planet equal to the radial
drift velocity of the gas in a thin accretion disc (see e.g. Pringle 1981), we find:
ȧ = −
so long as the disc is sufficiently massive. Note that Equation 3 is independent of the disc
surface density. If the disc mass is too low, then the torques (also proportional to Σ above)
will not be able to force the planet to migrate at this rate. Syer and Clarke (1995) and
Ivanov et al. (1999) have examined this limit.
Ida and Lin (2004) suggest an alternative prescription for determining the Type II mi-
gration rate. They balance the angular momentum change of the planet, 1
MpΩpaȧ with the
maximum viscous couple in the disc J̇ = 3
r2. Since we are using power law discs
here, J̇ will not have a maximum, so we use the nominal value at the planet’s orbit. This
leads to the prediction
ȧ = −3
Since this equation was obtained by balancing angular momentum, we might expect it to be
valid for a full range of disc masses.
Surprisingly, giant planet migration does not seem to have been subjected to a system-
atic test (this is in sharp contrast to the theory of Type I migration, which applies to low mass
planets). Some curiosities in the behaviour of migrating giant planets have been seen (e.g.
Schäfer et al. 2004), but these have not been explored in detail. The work of Nelson et al.
– 4 –
(2000) provides the most complete set of runs to date, but the physical parameters were not
varied in a regular fashion.
In this paper, we shall present a series of numerical experiments, following giant planets
migrating in a variety of accretion discs. We will then compare our results to equations 3
and 4. Equation 3 makes particularly strong predictions about the expected migration rates
- namely that the migration rate depends solely on the disc viscosity. This is the first time
such a test has been performed.
3. Numerical Set up
We use the Fargo code of Masset (2000a,b) to perform our calculations. Fargo is a
simple 2d polar mesh code dedicated to disc planet interactions. It is based upon a standard,
Zeus-like (Stone and Norman 1992) hydrodynamic solver, but owes its name to the Fargo
algorithm upon which the azimuthal advection is based. This algorithm avoids the restrictive
timestep typically imposed by the rapidly rotating inner regions of the disc, by permitting
each annulus of cells to rotate at its local Keplerian velocity and stitching the results together
again at the end of the timestep. The use of the Fargo algorithm typically lifts the timestep
by an order of magnitude, and therefore speeds up the calculation accordingly. The mesh
centre lies at the central star, so indirect terms coming from the planets and the disc are
included in the potential calculation. We make use of an non-reflecting inner boundary, to
prevent reflected waves from interfering with the calculations. The pitch angle of the wake
is evaluated in the WKB approximation. The inner ring of active cells is then copied to the
ghost cells, with an azimuthal shift appropriate to the pitch angle. Material which flows off
the inner boundary is not added to the star (nor does the planet itself accrete). At the outer
boundary, mass was added, to compensate for the viscous evolution of the system.
We use units normalised such that G = M
+Mp = 1, while the planet’s initial orbital
radius is set at a = 1. References to times in terms of ‘orbits’ should be understood to
mean “orbital times at the planet’s initial radius.” The grid extends between r = 0.4 and
r = 2.5. Scaled to Jupiter’s orbit, this grid roughly covers the area between the asteroid belt
and Saturn’s orbit. We assume a constant aspect ratio disc, with h/r = 0.05. We set the
mass ratio q = Mp/M∗ to be 10
−3, approximately equal to the Jupiter-Sun value. In our
parameter space search, we varied the disc surface density profile and viscosity. The surface
density was initially a power law:
Σ(r) = Σ0
and we take r0 = 1 (the planet’s initial orbital radius). Four values for δ are considered:
– 5 –
0, 0.5, 1 and 1.5. These represent a reasonable range of alternatives, from a theoretically
simple constant surface density disc to the canonical Minimum Mass Solar Nebula. We set
Σ0 through
qdisc =
which provides a quick estimate of the disc’s mass within the planet’s orbit. This estimate
is accurate for δ = 0 discs, but is an underestimate for larger δ values. We take four values
for q
: 5×10−4, 10−3, 2×10−3, and 3×10−3. The total disc mass lies between 1.9×10−3
(for q
= 5× 10−4 and δ = 1.5) and 0.018 (for q
= 3× 10−3 and δ = 0) in units of the
stellar mass. By way of comparison, the Minimum Mass Solar Nebula (MMSN) requires at
least 5MJ of gas in the vicinity of Jupiter’s orbit.
2 Thus, our lower mass discs are somewhat
sub-Minimum.
The viscosity is taken to be uniform, and has values ν = 10−4 and 10−5 in our units.
With a uniform viscosity, δ = 0.5 yields a disc with an initially stationary surface density pro-
file (cf equation 2.10 of Pringle 1981). These viscosities may be related to the α prescription
for viscosity, ν = αcsh using
ν = α
Ωr2 (7)
This implies that α varies with radius. With our aspect ratio, a viscosity of ν = 10−5 gives
α ≈ 4 × 10−3 at the planet’s initial orbital radius. Note that for the highest viscosity, the
gap opening criterion of Equation 1 is not satisfied. The tidal condition of Equation 2 is
always satisfied in our numerical experiments.
The gravitational effect of the planet on the disc is smoothed at 0.6 of the disc thickness
at the planet’s orbital radius:
φ = −
r2 + ǫ2
where ǫ = 0.6h. There are two motivations for this, the first being the desire to avoid having
a singularity wandering around the grid. The second is physical. The 2d approximation
becomes poor close to the planet, where the vertical distribution of material becomes im-
portant. The actual distance of material from the planet ceases to be the well approximated
by the in-plane distance, which would lead to the gravitational effect being over-estimated.
Accordingly, we soften the potential over distances comparable to the disc scale height. How-
ever, this softening length is still substantially smaller than the expected gap size and the
2We calculate this by comparing Jupiter’s metal content to that of the Sun, and assuming that both
condensed from the same gas cloud. Scaled to the Solar System, our grid roughly covers the region between
the asteroid belt and Saturn
– 6 –
planet’s Hill sphere. When calculating the torque the disc exerts on the planet, material from
within the Hill sphere is subject to an exponential cut off, for similar reasons. At the start
of each run, the planet is introduced gradually (over about one orbit), and is not initially
permitted to migrate. This is done to minimise the effect of transients caused by the sudden
appearance of a planet in a smooth disc. We considered release times of approximately one
orbit, and 100 & 1000 orbits. The computational grid is covered by 128 radial and 384
azimuthal cells (all uniformly spaced).
So far as possible, this setup mirrors that used in the comparison project presented by
de Val-Borro et al. (2006). In that comparison, the Fargo code was seen to give similar
results to other codes used to study the disc–planet interaction problem.
4. Development of the Gap
Since we introduce the planet into an initially unperturbed disc, there is a period of
rapid evolution, as the planet clears a gap. In this section, we shall discuss the development
of this gap.
In Figure 1, we trace the evolution of the gap in a ν = 10−5 disc. The initial surface
density profile had δ = 0.5 (cf equation 5). There are no surprises in this plot, when
compared with the many other numerical calculations of gap formation. We see that the
gap is mostly cleared in the first 100 orbits (note that the y-axis is logarithmic). The gap
lies roughly between the m = 2 Lindblad resonances (located at r = 0.6 and r = 1.3). For a
q = 10−3 planet, this distance is also comparable to the co-rotation region. This plot draws
attention to the fact that the planet never completely clears the gap. Even after 1000 orbits,
the surface density in the gap is around 3% of its initial value. The gap edge is covered by
roughly ten grid cells.
If we increase the viscosity to ν = 10−4, the gap becomes far less pronounced. We show
the development of the gap in this case in Figure 2. Again, most of the depletion occurs
during the first 100 orbits, but the total depletion is far less. The density has only dropped
to around 50% of its initial value. This is not unexpected - according to Equation 1, a
Jupiter mass planet should not be able to open a gap in such a viscous disc. Of course, the
condition of Equation 2 is satisfied.
Figures 1 and 2 both show that most of the gap clearing occurs withing the first 100
orbits. We shall therefore use this as our canonical release time below. However, we shall
show the effect of varying the release time as well.
– 7 –
Gaps in protoplanetary discs are known to vary smoothly with q and ν (see, e.g.
astro-ph/0608020), so what is taken to be a gap is somewhat arbitrary. We shall con-
tinue with both viscosity values, but we must bear in mind that in the high viscosity case
the gap is quite shallow.
5. Results
We will now present the results of our numerical experiments of a Jupiter mass planet
migrating in power law discs, grouped by viscosity. Such planets are conventionally assumed
to undergo Type II migration. We have two predictions for the migration rate of giant
planets in Equations 3 and 4. We shall compare our results to these predictions.
5.1. High Viscosity
The higher viscosity runs had ν = 10−4. This viscosity means that for a Jupiter mass
planet, the viscous gap opening criterion q > 40R−1 of Equation 1 is not quite satisfied.
However, the tidal condition of Equation 2 is fulfilled.
We shall discuss the results from the runs where the planet was released after 100
orbits. In Appendix A, we demonstrate that the release time is not significant. The orbital
evolution of these planets is plotted in Figure 3. We cut the y-axis at 0.6, since at that point
the m = 2 ILR of the planet encounters the edge of the grid. The migration rate of the
planet therefore becomes unreliable. We can see that the migration rate is a strong function
(or equivalently, Σ0). This is in direct contradiction to the prediction of Equation 3.
Notice also how the migration rate varies with a. Equation 3 predicts that ȧ ∝ a−1, which
we do not see. We see that the migration rate generally falls with a, which is consistent
with equation 4. Figure 3 shows that the reduction of ȧ with a falls as δ increases (that
is, there is a pronounced curve in the migrations for δ = 0, whereas those for δ = 1.5 are
almost straight lines). This is broadly consistent with equation 4, were the migration rate is
proportional to ȧ ∝ aΣ ≡ a1−δ. Complicating this is the viscous evolution of the disc itself,
which is probably the reason why ȧ is not strictly proportional to a1−δ (which would predict
accelerating migration for the δ = 1.5 case).
We have seen that the migration rate of a Jupiter mass planet in these discs is strongly
affected by the disc surface density. However, since the gap is not particularly deep, whether
Type II behaviour should be expected is debatable.
– 8 –
5.2. Medium Viscosity
Here, we examine the results from runs with ν = 10−5. Again, we shall examine the
case where the planet was released after 100 orbits first.
In Figure 4, we show the orbital migration for planets embedded in a variety of discs.
We see that planets embedded in higher surface density discs (parameterised by q
Equation 6) consistently undergo faster migration. The migration rate is roughly propor-
tional to the disc surface density. Note also the variation of ȧ with a. It is similar to that
seen for figure 3 above. Eccentricities again remained low.
5.3. Summary of results
In this section, we have presented a series of giant planet migration runs. Such planets
have generally been thought to undergo Type II migration. One formulation of the theory
predicts that the migration rate depends solely on the disc viscosity (Equation 3). However,
we have found that the migration rates vary systematically with disc surface density. Higher
disc surface densities give faster migration, which is consistent with Equation 4. There is a
weaker variation with disc viscosity, which is inconsistent with both predictions.
In Appendix A, we show that our conclusions are not affected by varying the time at
which the planet was released to migrate. We demonstrate that the grid resolution does not
affect our results in Appendix B.
6. Discussion
In Section 5, we presented a series of runs designed to study how giant planets migrate.
This migration of such planets is thought to be controlled by the viscous torque within
the disc. Two different rates have been suggested. In the first (Equation 3), the planet is
locked to the viscous evolution of the disc, and the migration rate depends solely on the disc
viscosity. The second (Equation 4) computes an angular momentum balance between the
planet and disc. In this theory, the migration rate also depends on the disc surface density
and planet mass.
We have seen that the migration rate we obtain varies strongly with disc surface density,
indicating that Equation 3 in not appropriate. Although Equation 4 is more promising,
we do not recover the same variation with viscosity. The higher viscosity runs underwent
more rapid migration, but the difference in migration rates was not an order of magnitude.
– 9 –
However, this result is not as robust as the variation with surface density, since the high
viscosity runs did not satify both gap opening criteria. Although Figure 2 shows that the
high viscosity runs are definitely in the non-linear regime, the gap itself was not especially
clean.
We have shown (Appendix A) that our results are not simply ‘turn-on’ transients. The
migration rates are not significantly affected by the time at which the planet is released
from a fixed orbit. Our resolution tests (Appendix B) demonstrate that our results are not
significantly affected by a doubling of the grid resolution.
Our neglect of material within the Hill sphere when calculating the torque is a point
of concern. D’Angelo et al. (2003) noted that most of the torque in their calculations came
from within the Hill sphere.3 However, the theory of Type II migration takes no account of
this material either. It is a simple 2d theory, which assumes that the planet is merely acting
as a ‘relay station’ for the disc’s viscous torques. If the flow structure within the Hill sphere
is of critical importance, then we should not expect giant planet migration to be as simple
as Section 2 suggests.
The accretion behaviour of the planet could also affect migration. This is directly linked
to the previous point about flow within the Hill sphere. Kley (1999) showed that even in the
presence of a gap, a planet could continue to accrete material from the disc.4 Similar results
were reported by Lubow et al. (1999) and Kley et al. (2001). In this paper, we did not allow
the planet to accrete, and this caused material to build up around the planet. Since we
attentuated the torque from within the Hill sphere, this would not have affected our results
directly. However, if accretion were allowed, then the planet could gain an appreciable
amount of mass. This would both alter the gap structure, and make it more difficult for the
disc to move the planet (due to the planet’s increased inertia). Related to this issue is the
recent finding of Lubow and D’Angelo (2006) that the accretion rate through the gap could
be over 10% of the viscous accretion rate in the main disc, despite the drop in gas density.
Finally, there is the matter of viscosity. In our numerical experiments, we used a physi-
cal viscosity in the Navier-Stokes equations. In reality, the ‘viscosity’ in protoplanetary discs
probably originates fromMHD turbulence, and calculations have shown (Winters et al. 2003)
that the gap structure obtained in an MHD calculation differs from that in a purely hydro-
dynamic one. In particular, the gaps tend to be wider and shallower. If material in the
3However, they did not perform a parameter space search like we have done here
4Note that this finding in itself implicitly contradicted the usual assumption that the planet isolates the
inner and outer discs
– 10 –
corotation region is important to determining the migration rate, then this alteration in
gap structure will cause further changes to the migration rates. Nelson (2005) has already
demonstrated that a magnetic turbulence strongly affects the migration of a low mass planet.
Although our computations do not include MHD turbulence, the theory of Type II migration
neglects it too, so this cannot be the reason for the differences we have observed.
When might we expect migration to proceed according to equation 3? We believe this
might be possible for a planet of moderate mass, in a cold, very low viscosity disc, which
is more massive than the planet. Our reasoning is as follows: equation 3 is based on the
assumption that the planet completely isolates the inner and outer discs. This is easiest
to achieve in a very low viscosity disc (cf equation 1), which is also cold (cf equation 2).
We also require the disc to be more massive than the planet, to ensure sufficient angular
momentum reserves are available. The gap will lie roughly between the m = 2 Lindblad
resonances, and we would want these particular resonances to be responsible for most of the
gap clearing (i.e. the m = 2 resonances themselves must dominate the disc’s viscous torque).
This is because we would need disc material to be kept well away from the corotation region
of the planet. Masset and Papaloizou (2003) showed that corotation torques can give rise to
extremely rapid migration - known as ‘runaway’ or Type III migration. A planet less massive
than Jupiter will have its corotation region inside its m = 2 Lindblad resonances. However,
Masset and Papaloizou found that such planets tended to undergo runaway migration. This
reinforces the need for the disc itself to have a very low viscosity, so that the gap is as clean
as possible.
6.1. Origin of the Torque
We shall now discuss the radial origin of the torque.
In figure 5, we show the torque profiles, Tz(r) acting on the planet after 100, 500 and
1000 orbits. The disc viscosity was ν = 10−5, and the surface density profile was initially flat
(δ = 0). The planet was held on a fixed orbit for the entire calculation. In computing the
torques, the same exponential cut off used with Fargo was applied. We see that most of
the torque is generated within the range 0.8 < r < 1.2, and that the torques from the inner
and outer discs have opposite signs. At later times, the torque peaks on either side of the
planet lessen. This is due to the gap emptying further. The outer peak (which is pushing
the planet inwards) also broadens, while the inner peak does not. This ultimately ensures
that the planet migrates inwards.
Figure 6 shows how the torque felt by the planet scales with q
∝ Σ0. These curves
– 11 –
are plotted for a δ = 0, ν = 10−5 disc after 1000 orbits. As we might expect from section 5,
we see that the strength of the torques is directly proportional to the value of q
. This
leads to the migration rate varying strongly with the surface density of the disc.
Finally, in figure 7, we show the effect of varying the initial disc power law, δ, on the
torque profiles. Again, the torque profiles are for a ν = 10−5 disc, and are plotted after 1000
orbits of the (fixed) planet. We see that the torques are very similar, regardless of the initial
δ value, indicating that the perturbations induced by the planet are not dependent on the
background structure of the disc. Again, this is expected from section 5.
Figures 5, 6 and 7 all show that the torque generation peaks at radii of r = 0.9 and r =
1.1 (roughly 1.5 Hill radii from the planet). Comparing to figure 1, we see that these locations
lie deep within the gap, which is interesting for a number of reasons. The peaks are close to
the cutoff radius generally applied to obtain the numerical factor in equation 1, namely the
planet’s Hill sphere. They are also well within the corotation region, raising the possibility
that corotation torques are affecting the orbital evolution of the planet. Unfortunately, at
this resolution, the Hill radius is only covered by four or five grid cells, and the torque peak
only lies seven grid cells from the planet itself. The smoothing lengths are also comparable
to these distances.
Our resolution tests (Appendix B) show that our resolution is adequate for the smooth-
ing lengths used. However, with so much torque being generated close to the planet, it is
likely that the smoothing is significantly affecting the torque. Reduction of the smoothing
lengths is obviously desirable, but unfortunately not possible in a two dimensional calcula-
tion. As noted in section 3, we must smooth the planet’s gravity at about the local scale
height of the disc in order to make a 2d calculation valid. In Appendix C, we show the
effect of reducing the exclusion radius for the calculation of the migration torque. The effect
appears to be minimal, but these results should be treated with some caution, due to the
issue noted above. With the structure of the flow close to the planet so obviously critical to
determining the migration rate, a full 3d calculation would be required to determine a robust
migration rate. However, we would be most surprised if such calculations undermined our
main conclusion that migration rates of massive planets are proportional to the disc mass.
7. Conclusion
In this paper, we have performed a systematic test of giant (Jupiter mass) planet migra-
tion. The migration rates we obtained varied strongly with the initial disc surface density,
and less strongly with the disc viscosity. We have shown that the simplest theory of Type
– 12 –
II migration, where the planet is locked to the viscous evolution of the disc (Equation 3), is
incorrect. An alternative formulation, based on an angular momentum balance (Equation 4),
looks more promising. However, we have not tested this second theory fully. We verified
that our results were not simply ‘turn-on’ transients, or purely the effect of low resolution.
Neither doubling the grid resolution, nor allowing the planet to clear its gap for 1000 orbits
affected our central finding.
Separate confirmation of our results, using a different code would be highly desir-
able. Although we have no reason to believe that Fargo is misbehaving, the work of
de Val-Borro et al. (2006) underlines how codes can give varying results, even for the ‘same’
physical scenario (it is for this reason that ‘simulations’ should properly be referred to as
‘numerical experiments’). The issues of accretion and gravitational softening (both of the
planet’s effect on the disc, and the torque exerted on the planet) also merit closer consider-
ation. Indeed, if the gap shape and flow through the gap are critical for migration, one is
led to wonder if 2d calculations are sufficient. Two dimensional calculations have generally
been thought adequate for Jovian mass planets because the gap would keep material away
from the planet (where the 3d nature of the flow will become evident). If the flow within
the Hill sphere is important, then 2d calculations cease to be convincing.
A. Effect of Varying Release Time
In this Appendix, we shall demonstrate that varying the time at which the planet is
released does not affect our central conclusion. We start by showing that changing the
release time only affects the migration rates slightly. We then show that, since the effect is
consistent for all discs used, the conclusion that the planet migration rate varies with disc
surface density is robust.
In Figure 8, we show the effect of varying the release time on a planet in a high viscosity
(ν = 10−4) disc with q
= 0.002 and δ = 0.5. The effect is fairly minimal, and this plot
is typical. The explanation for this lies in Figure 2, which shows that only a minimal gap is
formed. Indeed, one can debate whether the surface density depression is a gap, since the
usual criterion (Equation 1) is not satisfied.
Figure 9 shows the effect of release time on the orbital migration of a planet in the
ν = 10−5 case. The particular disc used in this comparison had q
= 0.002 and δ = 0,
but the behaviour was generically the same for all cases. We see that the time at which
the planet is released does have an effect on the migration rate. However, the effect is not
especially dramatic.
– 13 –
We demonstrate that the release time does not affect our central conclusions in Figure 10.
This duplicates Figure 4, but with the release time increased to 1000 orbits (note that the x
axis has the zero point shifted, to improve the use of space). Although the exact migration
rates undoubtedly change, the main conclusion that these planets are not undergoing Type
II migration is unaffected. The migration rates continue to be affected by the disc surface
density.
B. The Effect of Resolution
In this appendix, we shall study the effect of increasing the grid resolution on our results.
We re-ran four of our numerical experiments, but with the grid resolution doubled. We
picked the set of four runs with ν = 10−5, and δ = 0.5. In Figure 11 we plot the results,
compared to the top right panel of Figure 4. In this plot, we can see that doubling the grid
resolution has a minimal effect on the migration rates obtained.
From this, we see that our conclusions are not simply an artifact of low resolution. Of
course, the smoothing we have used could be hiding some effects, but we would not expect
decreasing the smoothing length to change the planet migration to Type II behaviour. The
smoothing length, at 0.6h is already rather smaller than the gap width, so shrinking it further
is unlikely to enable the planet to make the gap cleaner. Furthermore, this smoothing length
is already as small as we can realistically make it. The flow close to the planet will really
have a 3d structure, not resolved in these calculations. By forcing all material to lie in the
disc plane (as required by a 2d calculation), we effectively make it closer to the planet –
significantly so within the gap. By softening the potential over a distance comparable to the
scale height, we approximate the true 3d strength of the planet’s gravity.
C. The Effect of Hill Sphere Exclusion
By excluding material within the planet’s Hill sphere when computing the migration
torque, we potentially reduced the migration torque substantially. Although this has a sound
physical motivation (cf section 3), it does potentially affect the migration rate of the planet.
Figure 12 shows that this has negligible effect on our results. This shows the migration
of eight planets, embedded in ν = 10−5, δ = 0 discs. The four standard q
values were
considered, and the planets were held for 1000 orbits before being released to migrate. The
only difference between each pair of curves is whether the exclusion radius for the torque
calculation was a full Hill radius, for the solid lines, or half the Hill radius, for the dotted
– 14 –
lines. The solid lines in figure 12 are equivalent to the δ = 0 (top left) panel of figure 4,
up to the difference in release time. In every case, the migration of each pair of planets is
almost identical.
In figure 13, we show the effect of the exclusion radius reduction on the torque profiles.
This figure should be directly compared to figure 5. We can see that the torque close to
the planet is increased, particularly for the first curve, plotted after 100 orbits. There is
also a stronger peak close to the planet for all the curves. Otherwise, the torque profiles are
remarkably similar to figure 5. This is as expected, given the results shown in figure 12.
The author acknowledges support from NSF grants AST-0406799, AST-0098442, AST-
0406823, and NASA grants ATP04-0000-0016 and NNG04GM12G (issued through the Ori-
gins of Solar Systems Program). I would like thank Frederic Masset for use of the Fargo
code. I am also very grateful to Eric Blackman and Alice Quillen, for reading early drafts
of this manuscript. Some of the computations presented here used the resources of HPC2N,
Ume̊a
REFERENCES
Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., and Papaloizou, J. C. B.: 1999, ApJ
514, 344
D’Angelo, G., Kley, W., and Henning, T.: 2003, ApJ 586, 540
de Val-Borro, M., Edgar, R. G., Artymowicz, P., Ciecielag, P., Cresswell, P., D’Angelo, G.,
Delgado-Donate, E. J., Dirksen, G., Fromang, S., Gawryszczak, A., Klahr, H., Kley,
W., Lyra, W., Masset, F., Mellema, G., Nelson, R. P., Paardekooper, S.-J., Peplinski,
A., Pierens, A., Plewa, T., Rice, K., Schäfer, C., and Speith, R.: 2006, MNRAS 370,
Goldreich, P. and Tremaine, S.: 1978, ApJ 222, 850
Goldreich, P. and Tremaine, S.: 1979, ApJ 233, 857
Goldreich, P. and Tremaine, S.: 1980, ApJ 241, 425
Ida, S. and Lin, D. N. C.: 2004, ApJ 604, 388
Ivanov, P. B., Papaloizou, J. C. B., and Polnarev, A. G.: 1999, MNRAS 307, 79
Kley, W.: 1999, MNRAS 303, 696
– 15 –
Kley, W., D’Angelo, G., and Henning, T.: 2001, ApJ 547, 457
Lin, D. N. C. and Papaloizou, J. C. B.: 1993, in E. H. Levy and J. I. Lunine (eds.), Protostars
and Planets III, pp 749–835
Lin, D. N. C., Papaloizou, J. C. B., Terquem, C., Bryden, G., and Ida, S.: 2000, Protostars
and Planets IV p. 1111
Lubow, S. H. and D’Angelo, G.: 2006, ApJ 641, 526
Lubow, S. H., Seibert, M., and Artymowicz, P.: 1999, ApJ 526, 1001
Marcy, G., Butler, R. P., Fischer, D., Vogt, S., Wright, J. T., Tinney, C. G., and Jones,
H. R. A.: 2005, Progress of Theoretical Physics Supplement 158, 24
Masset, F.: 2000a, A&AS 141, 165
Masset, F. S.: 2000b, in ASP Conf. Ser. 219: Disks, Planetesimals, and Planets, p. 75
Masset, F. S. and Papaloizou, J. C. B.: 2003, ApJ 588, 494
Nelson, R. P.: 2005, A&A 443, 1067
Nelson, R. P., Papaloizou, J. C. B., Masset, F., and Kley, W.: 2000, MNRAS 318, 18
Pringle, J. E.: 1981, ARA&A 19, 137
Schäfer, C., Speith, R., Hipp, M., and Kley, W.: 2004, A&A 418, 325
Stone, J. M. and Norman, M. L.: 1992, ApJS 80, 753
Syer, D. and Clarke, C. J.: 1995, MNRAS 277, 758
Ward, W. R. and Hahn, J. M.: 2000, Protostars and Planets IV p. 1135
Winters, W. F., Balbus, S. A., and Hawley, J. F.: 2003, ApJ 589, 543
This preprint was prepared with the AAS LATEX macros v5.2.
– 16 –
Fig. 1.— Development of the gap for a Jupiter-mass planet on a fixed circular orbit, em-
bedded in a ν = 10−5 disc, with δ = 0.5 (cf equation 5). The azimuthally averaged surface
density profile is show after 10, 100 and 1000 orbits. Note that the y-axis is logarithmic
– 17 –
Fig. 2.— Development of the gap for a Jupiter-mass planet on a fixed circular orbit, em-
bedded in a ν = 10−4 disc, with δ = 0.5 (cf equation 5). The azimuthally averaged surface
density profile is show after 10, 100 and 1000 orbits. For ease of comparison, the y-axis is
identical to that in Figure 1
– 18 –
Fig. 3.— Planetary migration in the ν = 10−4 case for disc power laws of δ = 0 (top left), 0.5
(top right), 1.0 (bottom left) and 1.5 (bottom right). The planets were released to migrate
after 100 orbits. The lines are marked by the value of q
(see Equation 6)
– 19 –
Fig. 4.— Planetary migration in the ν = 10−5 case for disc power laws of δ = 0 (top left), 0.5
(top right), 1.0 (bottom left) and 1.5 (bottom right). The planets were released to migrate
after 100 orbits. The lines are marked by the value of q
(see Equation 6)
– 20 –
Fig. 5.— Radial torque profile for a planet in a ν = 10−5 disc as a function of time. The
disc initially had a flat (δ = 0) surface density profile, and the planet was on a fixed orbit
for the entire time
– 21 –
Fig. 6.— Radial torque profiles for a planet in ν = 10−5 discs of differing surface density
after 1000 orbits. The disc had a δ = 0 initial surface density profile, and each line is marked
with its q
value. The planet’s orbit was fixed
– 22 –
Fig. 7.— Radial torque profiles for a planet in ν = 10−5 discs of differing δ, plotted after
1000 orbits. All discs had the same q
value, and the planet was held on a fixed orbit
– 23 –
Fig. 8.— The effect of release time on the ν = 10−4 case. The orbital evolution of a planet
in a sample disc is shown for three different release times
– 24 –
Fig. 9.— The effect of release time on the ν = 10−5 case. The orbital evolution of a planet
in a sample disc is shown for three different release times
– 25 –
Fig. 10.— Planetary migration in the ν = 10−5 case for planets released after 1000 orbits.
The disc power laws are δ = 0 (top left), 0.5 (top right), 1.0 (bottom left) and 1.5 (bottom
right). The lines are marked by the value of q
(see Equation 6). This plot should be
compared to Figure 4
– 26 –
Fig. 11.— Resolution test for a ν = 10−5 disc with δ = 0.5, where the planet is released
after 100 orbits. Four different q
values are compared (refer to the top right panel of
Figure 4) at two grid resolutions
– 27 –
Fig. 12.— Effect of reducing the exclusion radius to half the Hill radius when computing the
migration torque on the planet. This figure is based on the δ = 0 (top left) panel of Figure 4,
except that the planets were released after 1000 orbits. The standard four q
migrations
are plotted. The solid lines correspond to those in Figure 4 (up to the difference in release
time), while the dotted lines trace planets where the torque exclusion was only half the Hill
radius
– 28 –
Fig. 13.— Effect on the torque profiles of reducing exclusion radius to half the Hill radius.
This figure should be compared to figure 5
	Introduction
	Type II Migration
	Numerical Set up
	Development of the Gap
	Results
	High Viscosity
	Medium Viscosity
	Summary of results
	Discussion
	Origin of the Torque
	Conclusion
	Effect of Varying Release Time
	The Effect of Resolution
	The Effect of Hill Sphere Exclusion

  Many extra-solar planets discovered over the past decade are gas giants in
tight orbits around their host stars. Due to the difficulties of forming these
`hot Jupiters' in situ, they are generally assumed to have migrated to their
present orbits through interactions with their nascent discs. In this paper, we
present a systematic study of giant planet migration in power law discs. We
find that the planetary migration rate is proportional to the disc surface
density. This is inconsistent with the assumption that the migration rate is
simply the viscous drift speed of the disc. However, this result can be
obtained by balancing the angular momentum of the planet with the viscous
torque in the disc. We have verified that this result is not affected by
adjusting the resolution of the grid, the smoothing length used, or the time at
which the planet is released to migrate.

Introduction
The past decade has seen the discovery of over one hundred extra-solar planets. Most of
these planets (see, e.g. Marcy et al. 2005) are gas giants of Jupiter mass or greater, in orbits
close to their host stars.1 Such planets are often called ‘hot Jupiters.’ Although detection
methods (primarily radial velocity measurements) are strongly biased towards detecting such
planets, a sufficient number have been discovered to require detailed explanation.
The present orbits of the ‘hot Jupiters’ lie so close to their host stars that in situ forma-
tion is almost impossible. The problem is two-fold. Firstly, the expected disc temperature
is high, which suppresses gravitational capture of gas by a gas giant core. Secondly, there
1For the most up to date information, refer to http://exoplanet.eu/ and http://exoplanets.org/
http://arxiv.org/abs/0704.0448v1
http://exoplanet.eu/
http://exoplanets.org/
– 2 –
simply isn’t much material at such small orbital radii. We are therefore led to the conclusion
that the ‘hot Jupiters’ formed further from the star, and migrated to their present orbits.
In this paper, we present a study of giant planet migration in power law discs.
A planet in a circumstellar disc exerts a torque on the disc material at resonances
(Goldreich and Tremaine 1978, 1979, 1980). By Newton’s Third Law, there is an opposite
torque on the planet, which causes it to migrate. Planets of mass comparable to Jupiter (the
exact threshold depends on the disc conditions, see below) exert torques comparable to the
viscous torque in the disc. They can therefore open a gap in the disc. If the gap is perfectly
clean, then the planet will act as a ‘relay station’ between the inner and outer discs. We then
expect the planet’s orbital evolution to follow the viscous evolution of the disc – a process
known as Type II migration (Ward and Hahn 2000). We shall compare our computational
results to the predictions of Type II theory.
We summarise the theory of Type II migration in Section 2. We describe our numerical
method in Section 3, and demonstrate the formation of a gap in Section 4. Section 5 presents
our main results. In Section 6 we discuss the implications of our findings, before presenting
our conclusions in Section 7.
2. Type II Migration
In this section, we shall briefly review the theory of Type II migration. For a more
thorough analysis of the theory of planet–disc interactions, the reader should turn to one of
the many published reviews (e.g. Lin and Papaloizou 1993; Ward and Hahn 2000; Lin et al.
2000).
Planet–disc interactions are dominated by resonances. A planet embedded in a cir-
cumstellar disc excites waves at its Lindblad resonances (LR). These waves carry angular
momentum, and hence exert a torque on the disc material. The strength of the torque from
the Lindblad resonances is proportional to TLR ∝ Σq
2, where Σ is the local disc surface
density, and q is the planet-star mass ratio. These torques act to push material away from
the planet. At the same time, the disc gas is expected to have an intrinsic viscosity, ν (al-
though the precise origin and exact behaviour of such a viscosity are still much debated),
which leads to a viscous torque Tν ∝ Σν. Since the Lindblad torques scale with q2, we can
expect that as the planet accretes material, the Lindblad torques will eventually dominate
the viscous torque in the disc. Balancing the two torques leads to the so-called gap opening
criterion:
q > 40R−1 (1)
– 3 –
where R = r2Ω/ν is the Reynolds number of the disc (Bryden et al. 1999; Nelson et al.
2000). There is a similar requirement that the planet’s Hill sphere exceed the local scale
height of the disc, namely that
q > 3
where h is the disc scale height. See Lin and Papaloizou (1993) for a further discussion.
The tidal condition of Equation 2 leads to the gap width being at least twice the local scale
height. If this condition is not satisfied, then the edge of the gap will be Rayleigh unstable.
Consideration of the torque condition leads to the expectation that the gap will lie between
the m = 2 Lindblad resonances of the planet (that is, between r = 0.6 and r = 1.3). For
Jovian mass planets in circumstellar discs, these gap sizes are similar.
Once the planet has opened a gap, it is assumed to isolate the inner and outer discs
from each other. However, each part of the disc will still be undergoing viscous evolution.
According to Type II migration theory, the planet will be locked to the viscous evolution of
the disc (Ward and Hahn 2000). Setting the migration rate of the planet equal to the radial
drift velocity of the gas in a thin accretion disc (see e.g. Pringle 1981), we find:
ȧ = −
so long as the disc is sufficiently massive. Note that Equation 3 is independent of the disc
surface density. If the disc mass is too low, then the torques (also proportional to Σ above)
will not be able to force the planet to migrate at this rate. Syer and Clarke (1995) and
Ivanov et al. (1999) have examined this limit.
Ida and Lin (2004) suggest an alternative prescription for determining the Type II mi-
gration rate. They balance the angular momentum change of the planet, 1
MpΩpaȧ with the
maximum viscous couple in the disc J̇ = 3
r2. Since we are using power law discs
here, J̇ will not have a maximum, so we use the nominal value at the planet’s orbit. This
leads to the prediction
ȧ = −3
Since this equation was obtained by balancing angular momentum, we might expect it to be
valid for a full range of disc masses.
Surprisingly, giant planet migration does not seem to have been subjected to a system-
atic test (this is in sharp contrast to the theory of Type I migration, which applies to low mass
planets). Some curiosities in the behaviour of migrating giant planets have been seen (e.g.
Schäfer et al. 2004), but these have not been explored in detail. The work of Nelson et al.
– 4 –
(2000) provides the most complete set of runs to date, but the physical parameters were not
varied in a regular fashion.
In this paper, we shall present a series of numerical experiments, following giant planets
migrating in a variety of accretion discs. We will then compare our results to equations 3
and 4. Equation 3 makes particularly strong predictions about the expected migration rates
- namely that the migration rate depends solely on the disc viscosity. This is the first time
such a test has been performed.
3. Numerical Set up
We use the Fargo code of Masset (2000a,b) to perform our calculations. Fargo is a
simple 2d polar mesh code dedicated to disc planet interactions. It is based upon a standard,
Zeus-like (Stone and Norman 1992) hydrodynamic solver, but owes its name to the Fargo
algorithm upon which the azimuthal advection is based. This algorithm avoids the restrictive
timestep typically imposed by the rapidly rotating inner regions of the disc, by permitting
each annulus of cells to rotate at its local Keplerian velocity and stitching the results together
again at the end of the timestep. The use of the Fargo algorithm typically lifts the timestep
by an order of magnitude, and therefore speeds up the calculation accordingly. The mesh
centre lies at the central star, so indirect terms coming from the planets and the disc are
included in the potential calculation. We make use of an non-reflecting inner boundary, to
prevent reflected waves from interfering with the calculations. The pitch angle of the wake
is evaluated in the WKB approximation. The inner ring of active cells is then copied to the
ghost cells, with an azimuthal shift appropriate to the pitch angle. Material which flows off
the inner boundary is not added to the star (nor does the planet itself accrete). At the outer
boundary, mass was added, to compensate for the viscous evolution of the system.
We use units normalised such that G = M
+Mp = 1, while the planet’s initial orbital
radius is set at a = 1. References to times in terms of ‘orbits’ should be understood to
mean “orbital times at the planet’s initial radius.” The grid extends between r = 0.4 and
r = 2.5. Scaled to Jupiter’s orbit, this grid roughly covers the area between the asteroid belt
and Saturn’s orbit. We assume a constant aspect ratio disc, with h/r = 0.05. We set the
mass ratio q = Mp/M∗ to be 10
−3, approximately equal to the Jupiter-Sun value. In our
parameter space search, we varied the disc surface density profile and viscosity. The surface
density was initially a power law:
Σ(r) = Σ0
and we take r0 = 1 (the planet’s initial orbital radius). Four values for δ are considered:
– 5 –
0, 0.5, 1 and 1.5. These represent a reasonable range of alternatives, from a theoretically
simple constant surface density disc to the canonical Minimum Mass Solar Nebula. We set
Σ0 through
qdisc =
which provides a quick estimate of the disc’s mass within the planet’s orbit. This estimate
is accurate for δ = 0 discs, but is an underestimate for larger δ values. We take four values
for q
: 5×10−4, 10−3, 2×10−3, and 3×10−3. The total disc mass lies between 1.9×10−3
(for q
= 5× 10−4 and δ = 1.5) and 0.018 (for q
= 3× 10−3 and δ = 0) in units of the
stellar mass. By way of comparison, the Minimum Mass Solar Nebula (MMSN) requires at
least 5MJ of gas in the vicinity of Jupiter’s orbit.
2 Thus, our lower mass discs are somewhat
sub-Minimum.
The viscosity is taken to be uniform, and has values ν = 10−4 and 10−5 in our units.
With a uniform viscosity, δ = 0.5 yields a disc with an initially stationary surface density pro-
file (cf equation 2.10 of Pringle 1981). These viscosities may be related to the α prescription
for viscosity, ν = αcsh using
ν = α
Ωr2 (7)
This implies that α varies with radius. With our aspect ratio, a viscosity of ν = 10−5 gives
α ≈ 4 × 10−3 at the planet’s initial orbital radius. Note that for the highest viscosity, the
gap opening criterion of Equation 1 is not satisfied. The tidal condition of Equation 2 is
always satisfied in our numerical experiments.
The gravitational effect of the planet on the disc is smoothed at 0.6 of the disc thickness
at the planet’s orbital radius:
φ = −
r2 + ǫ2
where ǫ = 0.6h. There are two motivations for this, the first being the desire to avoid having
a singularity wandering around the grid. The second is physical. The 2d approximation
becomes poor close to the planet, where the vertical distribution of material becomes im-
portant. The actual distance of material from the planet ceases to be the well approximated
by the in-plane distance, which would lead to the gravitational effect being over-estimated.
Accordingly, we soften the potential over distances comparable to the disc scale height. How-
ever, this softening length is still substantially smaller than the expected gap size and the
2We calculate this by comparing Jupiter’s metal content to that of the Sun, and assuming that both
condensed from the same gas cloud. Scaled to the Solar System, our grid roughly covers the region between
the asteroid belt and Saturn
– 6 –
planet’s Hill sphere. When calculating the torque the disc exerts on the planet, material from
within the Hill sphere is subject to an exponential cut off, for similar reasons. At the start
of each run, the planet is introduced gradually (over about one orbit), and is not initially
permitted to migrate. This is done to minimise the effect of transients caused by the sudden
appearance of a planet in a smooth disc. We considered release times of approximately one
orbit, and 100 & 1000 orbits. The computational grid is covered by 128 radial and 384
azimuthal cells (all uniformly spaced).
So far as possible, this setup mirrors that used in the comparison project presented by
de Val-Borro et al. (2006). In that comparison, the Fargo code was seen to give similar
results to other codes used to study the disc–planet interaction problem.
4. Development of the Gap
Since we introduce the planet into an initially unperturbed disc, there is a period of
rapid evolution, as the planet clears a gap. In this section, we shall discuss the development
of this gap.
In Figure 1, we trace the evolution of the gap in a ν = 10−5 disc. The initial surface
density profile had δ = 0.5 (cf equation 5). There are no surprises in this plot, when
compared with the many other numerical calculations of gap formation. We see that the
gap is mostly cleared in the first 100 orbits (note that the y-axis is logarithmic). The gap
lies roughly between the m = 2 Lindblad resonances (located at r = 0.6 and r = 1.3). For a
q = 10−3 planet, this distance is also comparable to the co-rotation region. This plot draws
attention to the fact that the planet never completely clears the gap. Even after 1000 orbits,
the surface density in the gap is around 3% of its initial value. The gap edge is covered by
roughly ten grid cells.
If we increase the viscosity to ν = 10−4, the gap becomes far less pronounced. We show
the development of the gap in this case in Figure 2. Again, most of the depletion occurs
during the first 100 orbits, but the total depletion is far less. The density has only dropped
to around 50% of its initial value. This is not unexpected - according to Equation 1, a
Jupiter mass planet should not be able to open a gap in such a viscous disc. Of course, the
condition of Equation 2 is satisfied.
Figures 1 and 2 both show that most of the gap clearing occurs withing the first 100
orbits. We shall therefore use this as our canonical release time below. However, we shall
show the effect of varying the release time as well.
– 7 –
Gaps in protoplanetary discs are known to vary smoothly with q and ν (see, e.g.
astro-ph/0608020), so what is taken to be a gap is somewhat arbitrary. We shall con-
tinue with both viscosity values, but we must bear in mind that in the high viscosity case
the gap is quite shallow.
5. Results
We will now present the results of our numerical experiments of a Jupiter mass planet
migrating in power law discs, grouped by viscosity. Such planets are conventionally assumed
to undergo Type II migration. We have two predictions for the migration rate of giant
planets in Equations 3 and 4. We shall compare our results to these predictions.
5.1. High Viscosity
The higher viscosity runs had ν = 10−4. This viscosity means that for a Jupiter mass
planet, the viscous gap opening criterion q > 40R−1 of Equation 1 is not quite satisfied.
However, the tidal condition of Equation 2 is fulfilled.
We shall discuss the results from the runs where the planet was released after 100
orbits. In Appendix A, we demonstrate that the release time is not significant. The orbital
evolution of these planets is plotted in Figure 3. We cut the y-axis at 0.6, since at that point
the m = 2 ILR of the planet encounters the edge of the grid. The migration rate of the
planet therefore becomes unreliable. We can see that the migration rate is a strong function
(or equivalently, Σ0). This is in direct contradiction to the prediction of Equation 3.
Notice also how the migration rate varies with a. Equation 3 predicts that ȧ ∝ a−1, which
we do not see. We see that the migration rate generally falls with a, which is consistent
with equation 4. Figure 3 shows that the reduction of ȧ with a falls as δ increases (that
is, there is a pronounced curve in the migrations for δ = 0, whereas those for δ = 1.5 are
almost straight lines). This is broadly consistent with equation 4, were the migration rate is
proportional to ȧ ∝ aΣ ≡ a1−δ. Complicating this is the viscous evolution of the disc itself,
which is probably the reason why ȧ is not strictly proportional to a1−δ (which would predict
accelerating migration for the δ = 1.5 case).
We have seen that the migration rate of a Jupiter mass planet in these discs is strongly
affected by the disc surface density. However, since the gap is not particularly deep, whether
Type II behaviour should be expected is debatable.
– 8 –
5.2. Medium Viscosity
Here, we examine the results from runs with ν = 10−5. Again, we shall examine the
case where the planet was released after 100 orbits first.
In Figure 4, we show the orbital migration for planets embedded in a variety of discs.
We see that planets embedded in higher surface density discs (parameterised by q
Equation 6) consistently undergo faster migration. The migration rate is roughly propor-
tional to the disc surface density. Note also the variation of ȧ with a. It is similar to that
seen for figure 3 above. Eccentricities again remained low.
5.3. Summary of results
In this section, we have presented a series of giant planet migration runs. Such planets
have generally been thought to undergo Type II migration. One formulation of the theory
predicts that the migration rate depends solely on the disc viscosity (Equation 3). However,
we have found that the migration rates vary systematically with disc surface density. Higher
disc surface densities give faster migration, which is consistent with Equation 4. There is a
weaker variation with disc viscosity, which is inconsistent with both predictions.
In Appendix A, we show that our conclusions are not affected by varying the time at
which the planet was released to migrate. We demonstrate that the grid resolution does not
affect our results in Appendix B.
6. Discussion
In Section 5, we presented a series of runs designed to study how giant planets migrate.
This migration of such planets is thought to be controlled by the viscous torque within
the disc. Two different rates have been suggested. In the first (Equation 3), the planet is
locked to the viscous evolution of the disc, and the migration rate depends solely on the disc
viscosity. The second (Equation 4) computes an angular momentum balance between the
planet and disc. In this theory, the migration rate also depends on the disc surface density
and planet mass.
We have seen that the migration rate we obtain varies strongly with disc surface density,
indicating that Equation 3 in not appropriate. Although Equation 4 is more promising,
we do not recover the same variation with viscosity. The higher viscosity runs underwent
more rapid migration, but the difference in migration rates was not an order of magnitude.
– 9 –
However, this result is not as robust as the variation with surface density, since the high
viscosity runs did not satify both gap opening criteria. Although Figure 2 shows that the
high viscosity runs are definitely in the non-linear regime, the gap itself was not especially
clean.
We have shown (Appendix A) that our results are not simply ‘turn-on’ transients. The
migration rates are not significantly affected by the time at which the planet is released
from a fixed orbit. Our resolution tests (Appendix B) demonstrate that our results are not
significantly affected by a doubling of the grid resolution.
Our neglect of material within the Hill sphere when calculating the torque is a point
of concern. D’Angelo et al. (2003) noted that most of the torque in their calculations came
from within the Hill sphere.3 However, the theory of Type II migration takes no account of
this material either. It is a simple 2d theory, which assumes that the planet is merely acting
as a ‘relay station’ for the disc’s viscous torques. If the flow structure within the Hill sphere
is of critical importance, then we should not expect giant planet migration to be as simple
as Section 2 suggests.
The accretion behaviour of the planet could also affect migration. This is directly linked
to the previous point about flow within the Hill sphere. Kley (1999) showed that even in the
presence of a gap, a planet could continue to accrete material from the disc.4 Similar results
were reported by Lubow et al. (1999) and Kley et al. (2001). In this paper, we did not allow
the planet to accrete, and this caused material to build up around the planet. Since we
attentuated the torque from within the Hill sphere, this would not have affected our results
directly. However, if accretion were allowed, then the planet could gain an appreciable
amount of mass. This would both alter the gap structure, and make it more difficult for the
disc to move the planet (due to the planet’s increased inertia). Related to this issue is the
recent finding of Lubow and D’Angelo (2006) that the accretion rate through the gap could
be over 10% of the viscous accretion rate in the main disc, despite the drop in gas density.
Finally, there is the matter of viscosity. In our numerical experiments, we used a physi-
cal viscosity in the Navier-Stokes equations. In reality, the ‘viscosity’ in protoplanetary discs
probably originates fromMHD turbulence, and calculations have shown (Winters et al. 2003)
that the gap structure obtained in an MHD calculation differs from that in a purely hydro-
dynamic one. In particular, the gaps tend to be wider and shallower. If material in the
3However, they did not perform a parameter space search like we have done here
4Note that this finding in itself implicitly contradicted the usual assumption that the planet isolates the
inner and outer discs
– 10 –
corotation region is important to determining the migration rate, then this alteration in
gap structure will cause further changes to the migration rates. Nelson (2005) has already
demonstrated that a magnetic turbulence strongly affects the migration of a low mass planet.
Although our computations do not include MHD turbulence, the theory of Type II migration
neglects it too, so this cannot be the reason for the differences we have observed.
When might we expect migration to proceed according to equation 3? We believe this
might be possible for a planet of moderate mass, in a cold, very low viscosity disc, which
is more massive than the planet. Our reasoning is as follows: equation 3 is based on the
assumption that the planet completely isolates the inner and outer discs. This is easiest
to achieve in a very low viscosity disc (cf equation 1), which is also cold (cf equation 2).
We also require the disc to be more massive than the planet, to ensure sufficient angular
momentum reserves are available. The gap will lie roughly between the m = 2 Lindblad
resonances, and we would want these particular resonances to be responsible for most of the
gap clearing (i.e. the m = 2 resonances themselves must dominate the disc’s viscous torque).
This is because we would need disc material to be kept well away from the corotation region
of the planet. Masset and Papaloizou (2003) showed that corotation torques can give rise to
extremely rapid migration - known as ‘runaway’ or Type III migration. A planet less massive
than Jupiter will have its corotation region inside its m = 2 Lindblad resonances. However,
Masset and Papaloizou found that such planets tended to undergo runaway migration. This
reinforces the need for the disc itself to have a very low viscosity, so that the gap is as clean
as possible.
6.1. Origin of the Torque
We shall now discuss the radial origin of the torque.
In figure 5, we show the torque profiles, Tz(r) acting on the planet after 100, 500 and
1000 orbits. The disc viscosity was ν = 10−5, and the surface density profile was initially flat
(δ = 0). The planet was held on a fixed orbit for the entire calculation. In computing the
torques, the same exponential cut off used with Fargo was applied. We see that most of
the torque is generated within the range 0.8 < r < 1.2, and that the torques from the inner
and outer discs have opposite signs. At later times, the torque peaks on either side of the
planet lessen. This is due to the gap emptying further. The outer peak (which is pushing
the planet inwards) also broadens, while the inner peak does not. This ultimately ensures
that the planet migrates inwards.
Figure 6 shows how the torque felt by the planet scales with q
∝ Σ0. These curves
– 11 –
are plotted for a δ = 0, ν = 10−5 disc after 1000 orbits. As we might expect from section 5,
we see that the strength of the torques is directly proportional to the value of q
. This
leads to the migration rate varying strongly with the surface density of the disc.
Finally, in figure 7, we show the effect of varying the initial disc power law, δ, on the
torque profiles. Again, the torque profiles are for a ν = 10−5 disc, and are plotted after 1000
orbits of the (fixed) planet. We see that the torques are very similar, regardless of the initial
δ value, indicating that the perturbations induced by the planet are not dependent on the
background structure of the disc. Again, this is expected from section 5.
Figures 5, 6 and 7 all show that the torque generation peaks at radii of r = 0.9 and r =
1.1 (roughly 1.5 Hill radii from the planet). Comparing to figure 1, we see that these locations
lie deep within the gap, which is interesting for a number of reasons. The peaks are close to
the cutoff radius generally applied to obtain the numerical factor in equation 1, namely the
planet’s Hill sphere. They are also well within the corotation region, raising the possibility
that corotation torques are affecting the orbital evolution of the planet. Unfortunately, at
this resolution, the Hill radius is only covered by four or five grid cells, and the torque peak
only lies seven grid cells from the planet itself. The smoothing lengths are also comparable
to these distances.
Our resolution tests (Appendix B) show that our resolution is adequate for the smooth-
ing lengths used. However, with so much torque being generated close to the planet, it is
likely that the smoothing is significantly affecting the torque. Reduction of the smoothing
lengths is obviously desirable, but unfortunately not possible in a two dimensional calcula-
tion. As noted in section 3, we must smooth the planet’s gravity at about the local scale
height of the disc in order to make a 2d calculation valid. In Appendix C, we show the
effect of reducing the exclusion radius for the calculation of the migration torque. The effect
appears to be minimal, but these results should be treated with some caution, due to the
issue noted above. With the structure of the flow close to the planet so obviously critical to
determining the migration rate, a full 3d calculation would be required to determine a robust
migration rate. However, we would be most surprised if such calculations undermined our
main conclusion that migration rates of massive planets are proportional to the disc mass.
7. Conclusion
In this paper, we have performed a systematic test of giant (Jupiter mass) planet migra-
tion. The migration rates we obtained varied strongly with the initial disc surface density,
and less strongly with the disc viscosity. We have shown that the simplest theory of Type
– 12 –
II migration, where the planet is locked to the viscous evolution of the disc (Equation 3), is
incorrect. An alternative formulation, based on an angular momentum balance (Equation 4),
looks more promising. However, we have not tested this second theory fully. We verified
that our results were not simply ‘turn-on’ transients, or purely the effect of low resolution.
Neither doubling the grid resolution, nor allowing the planet to clear its gap for 1000 orbits
affected our central finding.
Separate confirmation of our results, using a different code would be highly desir-
able. Although we have no reason to believe that Fargo is misbehaving, the work of
de Val-Borro et al. (2006) underlines how codes can give varying results, even for the ‘same’
physical scenario (it is for this reason that ‘simulations’ should properly be referred to as
‘numerical experiments’). The issues of accretion and gravitational softening (both of the
planet’s effect on the disc, and the torque exerted on the planet) also merit closer consider-
ation. Indeed, if the gap shape and flow through the gap are critical for migration, one is
led to wonder if 2d calculations are sufficient. Two dimensional calculations have generally
been thought adequate for Jovian mass planets because the gap would keep material away
from the planet (where the 3d nature of the flow will become evident). If the flow within
the Hill sphere is important, then 2d calculations cease to be convincing.
A. Effect of Varying Release Time
In this Appendix, we shall demonstrate that varying the time at which the planet is
released does not affect our central conclusion. We start by showing that changing the
release time only affects the migration rates slightly. We then show that, since the effect is
consistent for all discs used, the conclusion that the planet migration rate varies with disc
surface density is robust.
In Figure 8, we show the effect of varying the release time on a planet in a high viscosity
(ν = 10−4) disc with q
= 0.002 and δ = 0.5. The effect is fairly minimal, and this plot
is typical. The explanation for this lies in Figure 2, which shows that only a minimal gap is
formed. Indeed, one can debate whether the surface density depression is a gap, since the
usual criterion (Equation 1) is not satisfied.
Figure 9 shows the effect of release time on the orbital migration of a planet in the
ν = 10−5 case. The particular disc used in this comparison had q
= 0.002 and δ = 0,
but the behaviour was generically the same for all cases. We see that the time at which
the planet is released does have an effect on the migration rate. However, the effect is not
especially dramatic.
– 13 –
We demonstrate that the release time does not affect our central conclusions in Figure 10.
This duplicates Figure 4, but with the release time increased to 1000 orbits (note that the x
axis has the zero point shifted, to improve the use of space). Although the exact migration
rates undoubtedly change, the main conclusion that these planets are not undergoing Type
II migration is unaffected. The migration rates continue to be affected by the disc surface
density.
B. The Effect of Resolution
In this appendix, we shall study the effect of increasing the grid resolution on our results.
We re-ran four of our numerical experiments, but with the grid resolution doubled. We
picked the set of four runs with ν = 10−5, and δ = 0.5. In Figure 11 we plot the results,
compared to the top right panel of Figure 4. In this plot, we can see that doubling the grid
resolution has a minimal effect on the migration rates obtained.
From this, we see that our conclusions are not simply an artifact of low resolution. Of
course, the smoothing we have used could be hiding some effects, but we would not expect
decreasing the smoothing length to change the planet migration to Type II behaviour. The
smoothing length, at 0.6h is already rather smaller than the gap width, so shrinking it further
is unlikely to enable the planet to make the gap cleaner. Furthermore, this smoothing length
is already as small as we can realistically make it. The flow close to the planet will really
have a 3d structure, not resolved in these calculations. By forcing all material to lie in the
disc plane (as required by a 2d calculation), we effectively make it closer to the planet –
significantly so within the gap. By softening the potential over a distance comparable to the
scale height, we approximate the true 3d strength of the planet’s gravity.
C. The Effect of Hill Sphere Exclusion
By excluding material within the planet’s Hill sphere when computing the migration
torque, we potentially reduced the migration torque substantially. Although this has a sound
physical motivation (cf section 3), it does potentially affect the migration rate of the planet.
Figure 12 shows that this has negligible effect on our results. This shows the migration
of eight planets, embedded in ν = 10−5, δ = 0 discs. The four standard q
values were
considered, and the planets were held for 1000 orbits before being released to migrate. The
only difference between each pair of curves is whether the exclusion radius for the torque
calculation was a full Hill radius, for the solid lines, or half the Hill radius, for the dotted
– 14 –
lines. The solid lines in figure 12 are equivalent to the δ = 0 (top left) panel of figure 4,
up to the difference in release time. In every case, the migration of each pair of planets is
almost identical.
In figure 13, we show the effect of the exclusion radius reduction on the torque profiles.
This figure should be directly compared to figure 5. We can see that the torque close to
the planet is increased, particularly for the first curve, plotted after 100 orbits. There is
also a stronger peak close to the planet for all the curves. Otherwise, the torque profiles are
remarkably similar to figure 5. This is as expected, given the results shown in figure 12.
The author acknowledges support from NSF grants AST-0406799, AST-0098442, AST-
0406823, and NASA grants ATP04-0000-0016 and NNG04GM12G (issued through the Ori-
gins of Solar Systems Program). I would like thank Frederic Masset for use of the Fargo
code. I am also very grateful to Eric Blackman and Alice Quillen, for reading early drafts
of this manuscript. Some of the computations presented here used the resources of HPC2N,
Ume̊a
REFERENCES
Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., and Papaloizou, J. C. B.: 1999, ApJ
514, 344
D’Angelo, G., Kley, W., and Henning, T.: 2003, ApJ 586, 540
de Val-Borro, M., Edgar, R. G., Artymowicz, P., Ciecielag, P., Cresswell, P., D’Angelo, G.,
Delgado-Donate, E. J., Dirksen, G., Fromang, S., Gawryszczak, A., Klahr, H., Kley,
W., Lyra, W., Masset, F., Mellema, G., Nelson, R. P., Paardekooper, S.-J., Peplinski,
A., Pierens, A., Plewa, T., Rice, K., Schäfer, C., and Speith, R.: 2006, MNRAS 370,
Goldreich, P. and Tremaine, S.: 1978, ApJ 222, 850
Goldreich, P. and Tremaine, S.: 1979, ApJ 233, 857
Goldreich, P. and Tremaine, S.: 1980, ApJ 241, 425
Ida, S. and Lin, D. N. C.: 2004, ApJ 604, 388
Ivanov, P. B., Papaloizou, J. C. B., and Polnarev, A. G.: 1999, MNRAS 307, 79
Kley, W.: 1999, MNRAS 303, 696
– 15 –
Kley, W., D’Angelo, G., and Henning, T.: 2001, ApJ 547, 457
Lin, D. N. C. and Papaloizou, J. C. B.: 1993, in E. H. Levy and J. I. Lunine (eds.), Protostars
and Planets III, pp 749–835
Lin, D. N. C., Papaloizou, J. C. B., Terquem, C., Bryden, G., and Ida, S.: 2000, Protostars
and Planets IV p. 1111
Lubow, S. H. and D’Angelo, G.: 2006, ApJ 641, 526
Lubow, S. H., Seibert, M., and Artymowicz, P.: 1999, ApJ 526, 1001
Marcy, G., Butler, R. P., Fischer, D., Vogt, S., Wright, J. T., Tinney, C. G., and Jones,
H. R. A.: 2005, Progress of Theoretical Physics Supplement 158, 24
Masset, F.: 2000a, A&AS 141, 165
Masset, F. S.: 2000b, in ASP Conf. Ser. 219: Disks, Planetesimals, and Planets, p. 75
Masset, F. S. and Papaloizou, J. C. B.: 2003, ApJ 588, 494
Nelson, R. P.: 2005, A&A 443, 1067
Nelson, R. P., Papaloizou, J. C. B., Masset, F., and Kley, W.: 2000, MNRAS 318, 18
Pringle, J. E.: 1981, ARA&A 19, 137
Schäfer, C., Speith, R., Hipp, M., and Kley, W.: 2004, A&A 418, 325
Stone, J. M. and Norman, M. L.: 1992, ApJS 80, 753
Syer, D. and Clarke, C. J.: 1995, MNRAS 277, 758
Ward, W. R. and Hahn, J. M.: 2000, Protostars and Planets IV p. 1135
Winters, W. F., Balbus, S. A., and Hawley, J. F.: 2003, ApJ 589, 543
This preprint was prepared with the AAS LATEX macros v5.2.
– 16 –
Fig. 1.— Development of the gap for a Jupiter-mass planet on a fixed circular orbit, em-
bedded in a ν = 10−5 disc, with δ = 0.5 (cf equation 5). The azimuthally averaged surface
density profile is show after 10, 100 and 1000 orbits. Note that the y-axis is logarithmic
– 17 –
Fig. 2.— Development of the gap for a Jupiter-mass planet on a fixed circular orbit, em-
bedded in a ν = 10−4 disc, with δ = 0.5 (cf equation 5). The azimuthally averaged surface
density profile is show after 10, 100 and 1000 orbits. For ease of comparison, the y-axis is
identical to that in Figure 1
– 18 –
Fig. 3.— Planetary migration in the ν = 10−4 case for disc power laws of δ = 0 (top left), 0.5
(top right), 1.0 (bottom left) and 1.5 (bottom right). The planets were released to migrate
after 100 orbits. The lines are marked by the value of q
(see Equation 6)
– 19 –
Fig. 4.— Planetary migration in the ν = 10−5 case for disc power laws of δ = 0 (top left), 0.5
(top right), 1.0 (bottom left) and 1.5 (bottom right). The planets were released to migrate
after 100 orbits. The lines are marked by the value of q
(see Equation 6)
– 20 –
Fig. 5.— Radial torque profile for a planet in a ν = 10−5 disc as a function of time. The
disc initially had a flat (δ = 0) surface density profile, and the planet was on a fixed orbit
for the entire time
– 21 –
Fig. 6.— Radial torque profiles for a planet in ν = 10−5 discs of differing surface density
after 1000 orbits. The disc had a δ = 0 initial surface density profile, and each line is marked
with its q
value. The planet’s orbit was fixed
– 22 –
Fig. 7.— Radial torque profiles for a planet in ν = 10−5 discs of differing δ, plotted after
1000 orbits. All discs had the same q
value, and the planet was held on a fixed orbit
– 23 –
Fig. 8.— The effect of release time on the ν = 10−4 case. The orbital evolution of a planet
in a sample disc is shown for three different release times
– 24 –
Fig. 9.— The effect of release time on the ν = 10−5 case. The orbital evolution of a planet
in a sample disc is shown for three different release times
– 25 –
Fig. 10.— Planetary migration in the ν = 10−5 case for planets released after 1000 orbits.
The disc power laws are δ = 0 (top left), 0.5 (top right), 1.0 (bottom left) and 1.5 (bottom
right). The lines are marked by the value of q
(see Equation 6). This plot should be
compared to Figure 4
– 26 –
Fig. 11.— Resolution test for a ν = 10−5 disc with δ = 0.5, where the planet is released
after 100 orbits. Four different q
values are compared (refer to the top right panel of
Figure 4) at two grid resolutions
– 27 –
Fig. 12.— Effect of reducing the exclusion radius to half the Hill radius when computing the
migration torque on the planet. This figure is based on the δ = 0 (top left) panel of Figure 4,
except that the planets were released after 1000 orbits. The standard four q
migrations
are plotted. The solid lines correspond to those in Figure 4 (up to the difference in release
time), while the dotted lines trace planets where the torque exclusion was only half the Hill
radius
– 28 –
Fig. 13.— Effect on the torque profiles of reducing exclusion radius to half the Hill radius.
This figure should be compared to figure 5
	Introduction
	Type II Migration
	Numerical Set up
	Development of the Gap
	Results
	High Viscosity
	Medium Viscosity
	Summary of results
	Discussion
	Origin of the Torque
	Conclusion
	Effect of Varying Release Time
	The Effect of Resolution
	The Effect of Hill Sphere Exclusion
