Giant Planet Migration in Viscous Power-Law Discs

TL;DR

This study systematically examines giant planet migration in viscous power-law discs, revealing that migration rate depends on disc surface density and is governed by angular momentum balance rather than viscous drift speed.

Contribution

It provides a new understanding of migration dynamics in power-law discs, challenging previous assumptions about viscous drift speed influence.

Findings

Migration rate is proportional to disc surface density.

Migration governed by angular momentum balance, not viscous drift.

Results are robust across different simulation parameters.

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.

Figures

Equations

$$q>40{\mathcal{R}}^{-1}$$

$$q>3{\left(\frac{h}{r}\right)}^3$$

$$\dot{a}=-\frac{3\nu }{2a}$$

$$\dot{a}=-3\frac{\nu \Sigma a}{M_{\text{p}}}$$

$$\Sigma (r)={\Sigma }_0{\left(\frac{r}{r_0}\right)}^{-\delta }$$

$$q_{\text{disc}}=\frac{{\Sigma }_0\pi a^2}{M_{\ast }}$$

$$\nu =\alpha {\left(\frac{h}{r}\right)}^2\Omega r^2$$

$$\varphi =-\frac{GM}{\sqrt{r^2+{ϵ}^2}}$$

Full text

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.111For the most up to date information, refer to http://exoplanet.eu/ and http://exoplanets.org/ 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 formation 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 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 circumstellar 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 $T_{\textrm{LR}}\propto\Sigma q^{2}$, where $\Sigma$ 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, $\nu$ (although the precise origin and exact behaviour of such a viscosity are still much debated), which leads to a viscous torque $T_{\nu}\propto\Sigma\nu$. Since the Lindblad torques scale with $q^{2}$, 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:

[TABLE]

where $\mathcal{R}=r^{2}\Omega/\nu$ 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

[TABLE]

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:

[TABLE]

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 $\Sigma$ 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 migration rate. They balance the angular momentum change of the planet, $\frac{1}{2}M_{\textrm{p}}\Omega_{\textrm{p}}a\dot{a}$ with the maximum viscous couple in the disc $\dot{J}=\frac{3}{2}\Sigma\nu\Omega_{\textrm{kep}}r^{2}$. Since we are using power law discs here, $\dot{J}$ will not have a maximum, so we use the nominal value at the planet’s orbit. This leads to the prediction

[TABLE]

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 systematic 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. Schäfer et al., 2004), but these have not been explored in detail. The work of Nelson et al. (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_{*}+M_{\textrm{p}}=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=M_{\textrm{p}}/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:

[TABLE]

and we take $r_{0}=1$ (the planet’s initial orbital radius). Four values for $\delta$ are considered: 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 $\Sigma_{0}$ through

[TABLE]

which provides a quick estimate of the disc’s mass within the planet’s orbit. This estimate is accurate for $\delta=0$ discs, but is an underestimate for larger $\delta$ values. We take four values for $q_{\textrm{disc}}$: $5\times 10^{-4}$, $10^{-3}$, $2\times 10^{-3}$, and $3\times 10^{-3}$. The total disc mass lies between $1.9\times 10^{-3}$ (for $q_{\textrm{disc}}=5\times 10^{-4}$ and $\delta=1.5$) and $0.018$ (for $q_{\textrm{disc}}=3\times 10^{-3}$ and $\delta=0$) in units of the stellar mass. By way of comparison, the Minimum Mass Solar Nebula (MMSN) requires at least ${5}\,{M_{J}}$ of gas in the vicinity of Jupiter’s orbit.222We 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 Thus, our lower mass discs are somewhat sub-Minimum.

The viscosity is taken to be uniform, and has values $\nu=10^{-4}$ and $10^{-5}$ in our units. With a uniform viscosity, $\delta=0.5$ yields a disc with an initially stationary surface density profile (cf equation 2.10 of Pringle, 1981). These viscosities may be related to the $\alpha$ prescription for viscosity, $\nu=\alpha c_{s}h$ using

[TABLE]

This implies that $\alpha$ varies with radius. With our aspect ratio, a viscosity of $\nu=10^{-5}$ gives $\alpha\approx 4\times 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:

[TABLE]

where $\epsilon=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 important. 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. However, this softening length is still substantially smaller than the expected gap size and the 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 $\nu=10^{-5}$ disc. The initial surface density profile had $\delta=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 $\nu=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.

Gaps in protoplanetary discs are known to vary smoothly with $q$ and $\nu$ (see, e.g. astro-ph/0608020), so what is taken to be a gap is somewhat arbitrary. We shall continue 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 $\nu=10^{-4}$. This viscosity means that for a Jupiter mass planet, the viscous gap opening criterion $q>40\mathcal{R}^{-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 of $q_{\textrm{disc}}$ (or equivalently, $\Sigma_{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 $\dot{a}\propto 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 $\dot{a}$ with $a$ falls as $\delta$ increases (that is, there is a pronounced curve in the migrations for $\delta=0$, whereas those for $\delta=1.5$ are almost straight lines). This is broadly consistent with equation 4, were the migration rate is proportional to $\dot{a}\propto a\Sigma\equiv a^{1-\delta}$. Complicating this is the viscous evolution of the disc itself, which is probably the reason why $\dot{a}$ is not strictly proportional to $a^{1-\delta}$ (which would predict accelerating migration for the $\delta=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.

5.2 Medium Viscosity

Here, we examine the results from runs with $\nu=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_{\textrm{disc}}$ – cf Equation 6) consistently undergo faster migration. The migration rate is roughly proportional to the disc surface density. Note also the variation of $\dot{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. 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.333However, they did not perform a parameter space search like we have done here 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.444Note that this finding in itself implicitly contradicted the usual assumption that the planet isolates the inner and outer discs 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 physical viscosity in the Navier-Stokes equations. In reality, the ‘viscosity’ in protoplanetary discs probably originates from MHD turbulence, and calculations have shown (Winters et al., 2003) that the gap structure obtained in an MHD calculation differs from that in a purely hydrodynamic one. In particular, the gaps tend to be wider and shallower. If material in the 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, $T_{z}(r)$ acting on the planet after 100, 500 and 1000 orbits. The disc viscosity was $\nu=10^{-5}$, and the surface density profile was initially flat ($\delta=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_{\textrm{disc}}\propto\Sigma_{0}$. These curves are plotted for a $\delta=0$, $\nu=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_{\textrm{disc}}$. 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, $\delta$, on the torque profiles. Again, the torque profiles are for a $\nu=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 $\delta$ 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 smoothing 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 calculation. 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 migration. 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 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 desirable. 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 consideration. 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.

Appendix 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 ($\nu=10^{-4}$) disc with $q_{\textrm{disc}}=0.002$ and $\delta=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 $\nu=10^{-5}$ case. The particular disc used in this comparison had $q_{\textrm{disc}}=0.002$ and $\delta=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.

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.

Appendix 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 $\nu=10^{-5}$, and $\delta=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.

Appendix 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 $\nu=10^{-5}$, $\delta=0$ discs. The four standard $q_{\textrm{disc}}$ 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 lines. The solid lines in figure 12 are equivalent to the $\delta=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 Origins 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å