Planet Migration

Planet Migration Stephen H. Lubow Space Telescope Science Institute

Shigeru Ida Tokyo Institute of Technology Planet migration is the process by which a planet’s orbital radius changes in time. The main agent for causing planet migration is the gravitational interaction of a young planet with the disk from which it forms. We describe the migration rates resulting from these interactions based on a simple model for disk properties. These migration rates are higher than is reasonable for planet survival. We discuss some proposed models for which the migration rates are lower. The major uncertainties in migration rates are due to a lack of knowledge about the detailed physical properties of disks. We also describe some additional forms of migration.

1. INTRODUCTION

many giant planets with orbital radii substantially smaller than Jupiter’s added to the evidence (e.g., Butler et al, 2006). Their existence suggested that migration had occurred, since giant planet formation close to a star is not likely to occur (Bodenheimer et al, 2000; Ida & Lin, 2004). In Section 2 we describe calculations of migration rates due to disk-planet interactions. Section 2.2 describes the migration rate of a low mass planet interacting with a disk, modeled simply as a set of ballistic particles. This model assumes that disk structure is unaffected by the presence of the planet, the so-called Type I regime of planet migration. This section describes the essential physics of the Type I migration torque and provides an estimate of the migration rate. This estimate reveals a problem with the shortness of migration timescales. With such short timescales, it is difficult for cores of gas giants to form, which is apparently inconsistent with observationally inferred fraction of stars having giant planets (& 10%). This disparity is one of the most serious problems in the standard planet formation model that is based on core accretion. Sections 2.3 through 2.7 describe improvements to the migration rate calculation in Section 2.2, such as accounting for the effects of gas pressure. A reader interested in the the essential basics of migration could skip Sections 2.3-2.5. Sections 2.2-2.6 assume that the disk structure is unaffected by the presence of the planet. Section 2.7 describes effects of the back reaction of the disk due to the presence of a higher mass planet. The back reaction can lead to the formation of a gap in the disk, resulting in a different form of migration called Type II migration. Section 3 describes some of the outstanding questions in planet migration theory. Section 3 also describes some other forms of migration.

Planet formation occurs in disks of material that orbit stars (see Chapter 19). At early stages of evolution, the disks are largely gaseous and have masses of typically a few percent or more of the stellar masses. At later times, after several 106 y, the gas disperses and disks largely consist of solid material. Such disks last for timescales of ∼ 108 y. The interactions of a young planet with its surrounding disk affects the planet’s orbital energy and angular momentum. One consequence of such interactions is that a planet may move radially through the disk (toward or away from the central star), a process called migration that is the topic of this chapter. Although disks survive only a small fraction of the lifetime of the star and planet (∼ 1010 y), calculations of disk-planet interactions described in this chapter show them to be strong enough to have caused substantial migration. Early theoretical studies of planet migration such as Goldreich & Tremaine, 1980 understandably concentrated on the role of migration for the planets in our solar system. These studies suggested that there was substantial migration of Jupiter while it was immersed in the solar nebula. Migration should have occurred both at Jupiter’s early stages of formation when it was a solid core and at later stages when it was a fully developed gas giant planet interacting with the solar nebula. However, the evidence for migration was not apparent. The location of Jupiter is just outside the snow-line, where conditions for rapid core formation are most favorable. This situation seemed to suggest that Jupiter formed and remained near its current location. Therefore, a plausible conclusion was that migration did not play an important role for Jupiter and therefore perhaps for all planets. The situation for Jupiter currently remains somewhat of a puzzle. Evidence for the importance of migration changed dramatically in 1995 with the discovery of the first giant planet, 51 Peg b, that has an orbital period of only about 4 days (Mayor and Queloz, 1995). Subsequent discoveries of

2. Migration Rates Due To Disk-Planet Interactions 2.1 Aerodynamic Gas Drag Aerodynamic gas drag provides the dominant influence 1

somewhat greater than a few times RH . The solid line in Fig. 1 shows the path of a particle deflected by a planet whose mass is 10−6 Ms , Hill radius RH ≃ 0.007rp , and orbital separation r − rp ≃ 3.5RH . The planet and particle orbit counter-clockwise in the inertial frame with angular speeds Ωp and Ω(r), respectively. In the frame of the planet, the particle moves downward in the figure, since its angular speed is slower than that of the planet (Ω(r) < Ωp for r > rp ). To estimate the angular momentum change of a particle like that in Fig. 1, we consider a Cartesian coordinate system centered on the planet. The x axis lies along a line between the star and planet and points away from the star. The dashed line in Fig. 1 traces the path that the patricle would take in the absence of the planet, while the solid line shows the path in the presence of the planet. In both cases, the particle path is generally along the negative y direction. The two paths are nearly identical prior to the encounter with the planet (for y > 0). In that case, the particle velocity in the frame that corotates with the planet is then given by

on the orbital evolution of low mass solids embedded in a gaseous disk that orbits a star (Weidenschilling, 1977; Cuzzi & Weidenschilling, 2006). The gas is subject to a radial pressure force, in addition to the gravitational force of the star. For a disk with smooth structural variations in radius, this force induces slight departures from Keplerian speeds of order (H/r)2 Ωr, where H is the disk thickness and Ω is the angular speed of the disk at radius r from the star. For a thin disk (H ≪ r) whose pressure declines in radius, as is typically expected, the pressure force acts radially outward. The gas rotation rate required to achieve centrifugal balance is then slightly below the Keplerian rate. Sufficiently high mass solids are largely dynamically decoupled from the gas and orbit at nearly the Keplerian rate, but are subject to drag forces from the more slowly rotating gas. The drag leads to their orbit decay. Both the drag force and the inertia increase with the size of the solid, R. The drag force on the object increases with its area (∼ R2 ), while its inertia increases with its volume (∼ R3 ). In the high mass solid regime, the dominance of inertia over drag causes the orbital decay rate to decrease with object size as 1/R or with object mass as 1/M 1/3 . As we will see in the next section, the orbital migration rate due to gravitational interactions between an object and a disk increases with the mass of the object. There is then a cross-over mass above which the orbital changes due to disk gravitational forces dominates over those due to drag forces. For typical parameters, this value is much less than an Earth mass, 1M⊕ . Consequently, for the purposes of planet migration, we will ignore the effects of aerodynamic gas drag.

u ≃ r(Ω(r) − Ωp )ey ≃ x rp

ux ∼

As an initial description of gravitational disk-planet interactions, we consider the disk to consist of noninteracting particles, each having mass M much less than the mass of the planet Mp . The description is along the lines of Lin & Papaloizou, 1979. The particles are considered to encounter and pass by the planet that is on a fixed circular orbit of radius rp about the star of mass Ms . As a result of the interaction, the planet and disk exchange energy and angular momentum. This situation is a special case of the famous three-body problem in celestial mechanics. The tidal or Hill radius of the planet where planetary gravitational forces dominate over stellar and centrifugal forces is given by RH = rp

Mp 3Ms

1/3

,

(2)

Upon interaction with the planet, the particle is deflected slightly toward it. The particle of mass M experiences a force Fx = −GM Mp x/(x2 + y(t)2 )3/2 . This force dominantly acts over a time t when |y(t)| . |x| and is Fx ∼ −GM Mp /x2 . From equation (2), it follows that the encounter time ∆t ∼ |x/u| ∼ 1/Ωp , of order the orbital period of the planet, independent of x. As a result of the encounter, the particle acquires an x velocity

2.2 Ballistic Particle Model



3 dΩ e y ≃ − Ωp x e y . dr 2

GMp Fx ∆t ∼− 2 . M x Ωp

(3)

The particle is then deflected by an angle δ∼

ux ∼ u



Mp Ms



rp 3 ∼ x



RH x

3

(4)

after the encounter (see Fig. 1). To determine the change in angular momentum of the particle, we need to determine its change in velocity along the θ direction, that is the same as the y direction near the planet. To determine this velocity change, ∆uy , we ignore the effects of the star during the encounter and apply conservation of kinetic energy between the start and end of the encounter in the frame of the planet. The velocity magnitude u is then the same before and after the encounter, although the direction changes by angle δ.1 Since the particle in Fig. 1

(1)

where radius RH is measured from the center of the planet. Consider a particle that approaches the planet on a circular orbit about the star with orbital radius radius r sufficiently different from rp to allow it to freely pass by the planet with a small deflection. This condition requires that the closest approach between the particle and planet ∼ |r − rp | to be

1 We

are applying the so-called impulse approximation. The approximation involves the assumption that the duration of the interaction is much shorter than the orbital period and is only marginally satisfied here. Consequently, the expressions for deflection angle δ and the torque T cannot be determined with high accuracy in this approximation. They contain the proper dependences on various physical quantities. But the approximation

2

on Mp is a consequence of the deviations in torque between the perturbed and the unperturbed paths, since the angular momentum change along the unperturbed path is zero. These deviations involve the product of the linear dependence of the force on planet mass with the linear dependence of the path deflection on planet mass (equation (4)) and so are quadratic in Mp . As indicated by equation (7), the particle in Fig. 1 gains angular momentum. The reason is that the path deflection occurs mainly after the particle passes the planet, that is for y < 0. The planet then pulls the particle toward positive y, causing it to gain angular momentum. Just the opposite would happen for a particle with r < rp , for x < 0 in equation (7). The particle would approach the planet in the positive y direction, be deflected towards the planet for y > 0, and be pulled by the planet in the negative y direction, causing it to lose angular momentum. This process behaves like friction. The particle gains (loses) angular momentum if it moves slower (faster) than the planet. The angular momentum then flows outward as a result of the interactions. That is, for dΩ/dr < 0 as in the Keplerian case, a particle whose orbit lies interior to the planet gives angular momentum to the planet, since the planet has a lower angular speed than the particle. The planet in turn gives angular momentum to a particle whose orbit lies exterior to it. Fig. 4 shows the results of numerical tests of equation (7). It verifies the dependence of ∆J on x/rp and Mp /Ms . Departures of the expected dependences (solid lines) occur when x ≃ 3RH . At somewhat smaller values of x, the particle orbits do not pass smoothly by the planet. We now apply equation (7) to determine the torque on the planet for a set of particles that form a continuous disk with surface density Σ that we take to be constant in the region near the planet. The particle disk provides a flux of mass past the planet between x and x + dx

moves in the negative y direction, its pre-encounter y velocity is −u and its post-encounter y velocity is −u cos δ. We then have that the change of the velocity of the particle along the y direction is ∆uy = −u cos δ + u,

(5)

where u ∼ Ωp x is the velocity before the encounter (see equation (2)). We assume that the perturbation is weak, δ << 1, and obtain from equations (2), (4), and (5) 2    rp 5 Mp . (6) ∆uy ∼ u δ 2 ∼ rp Ωp Ms x It then follows that the change in angular momentum of the particle is given by  2   Mp rp 5 ∆J ∼ M rp ∆uy ∼ M rp2 Ωp . (7) Ms x Fig. 2 plots the particle orbital radius as a function of time for the particle plotted in Fig. 1, where t = 0 is the time of closest approach to the planet. Notice that after the encounter, the particle orbit acquires an eccentricity, as seen by the radial oscillations, and an increased angular momentum, as seen by the mean shift of the radius for the oscillations. Fig. 3 plots the torque on a particle due to the planet as a function of time. The solid curves are for particles that follow their actual paths (similar to the solid line in Fig. 1). The dashed curves are for particles that are made to follow the unperturbed paths (similar to the dashed line in Fig. 1). The net angular momentum change is the time-integrated torque. The torque of a particle along an unperturbed path (such as the dashed curve in Fig. 1) is antisymmetric in the time t. Consequently, there is no net change in angular momentum accumulated along this path. The change in angular momentum along the unperturbed path from t = −∞ to t = 0 is linear in the planet mass, since the force due to the planet is proportional to its mass. The departures from antisymmetry of the torque versus time result in the net angular momentum change. These departures are a consequence of the path deflection (solid curve in Fig. 1). For the case of the closer encounter plotted in the top panel, the departures from the unperturbed case are substantial. But for a slightly larger orbit (bottom panel), the torque along the perturbed path differs only slightly from the torque along the unperturbed path that integrates to zero. This behavior is a consequence of the steep decline of ∆J with x in equation (7). For x greater than a few RH , the angular momentum change acquired before the encounter is nearly completely canceled by the angular momentum change at later times t > 0. The net angular momentum change is given by equation (7) and is quadratic in the planet mass. The quadratic dependence

dM˙ ∼ uy Σ dx ∼ Σ Ωp x dx.

(8)

We evaluate the torque Tout on the planet due to disk material that extends outside the orbit of the planet from r = rp + ∆r to ∞ or x from ∆r to ∞, where ∆r > 0. We use the fact that the torque the planet exerts on the disk is equal and opposite to the torque the disk exerts on the planet. We then have Tout

∼

Tout

=

∆J dM˙ dx, ∆r M dx  2  Mp rp 3 4 2 −CT Σrp Ωp , Ms ∆r −

Z

∞

(9) (10)

where CT is a dimensionless positive constant of order unity. The torque on the planet due to the disk interior to the orbit of the planet from x = −∞ to x = ∆r with ∆r < 0 evaluates to Tin = −Tout or  3 2  Mp rp 2 4 Tin = CT ΣΩp rp . (11) Ms |∆r|

does not lead to the correct dimensionless numerical coefficients of proportionality for these quantities. An exact treatment of the particle torque in the limit of weak perturbations exerted by a small mass planet is given in Goldreich & Tremaine, 1982.

3

assumption that particles pass by the planet with little deflection is invalid in this region. Using equation (15) with |∆r| ∼ RH , we obtain a torque 4/3  Mp 2 4 . (16) |T | ∼ ΣΩp rp Ms

The equations of motion for particles subject only to gravitational forces are time-reversible. If we time-reverse the particle-planet encounter in Fig. 1, we see that the eccentric orbit particle would approach the planet (both on clockwise orbits) and then lose angular momentum (apply t → −t in Fig. 2). What determines whether a disk particle gains or loses angular momentum? We have assumed that the particles always approach the perturber on circular orbits. The particles in the disk re-encounter the planet on the synodic period, of order rp /x planetary orbital periods. But we saw in Fig. 2 that the particles acquire eccentricity after the encounter. For this model to be physically consistent, we require that this eccentricity damp on the synodic period. The eccentricity damping produces an arrow of time for the angular momentum exchange process that favors circular orbits ahead of the encounter as shown in Fig. 1, resulting in a gain of angular momentum for particles with r > rp . Equations (10) and (11) have important consequences. The torques on the planet arising from the inner and outer disks are quite powerful and oppose each other. This does not mean that the net torque on the planet is zero because we have assumed perfect symmetry across r = rp . The symmetry is broken by higher order considerations, such as the radial density gradient. For this reason, migration torques are often referred to as differential torques. Since the torques are singular in ∆r, they are dominated by material that comes close to the planet. Consequently, the asymmetries occur through differences in physical quantities at radial distances ∆r from the planet. Consider for example the effect of the density variation in radius that we have ignored up to this point. If we expand the disk density in a Taylor expansion about the orbit of the planet, we have that dΣ , (12) Σout − Σin ≃ 2∆r dr where Σin and Σout are the surface densities at r = rp − ∆r and r = rp + ∆r, respectively. Consequently, for |dΣ/dr| ∼ Σp /rp , we expect that the sum of the inner and outer torques to be smaller than their individual values by an amount of order ∆r/rp . Similar considerations apply to variations in other quantities. That is, we have that the absolute value of the net torque T on the planet is approximately given by |T | = ∼ ∼

|Tin + Tout |, |∆r| , |Tin | rp 2 2   rp Mp 2 4 . ΣΩp rp Ms |∆r|

Another limit on ∆r comes about due to gas pressure. We have not yet described gas pressure effects, but will do so in Sections 2.4 and 2.5. One effect of gas pressure is to cause the disk to have a nonzero thickness H. The gas density is then smeared over distance H out of the orbit plane. Near the planet, the gas gravitational effects are smoothed over distance H. Distance ∆r is then in effect limited to H. The torque is then estimated as  2   Mp rp 2 |T | ∼ ΣΩ2p rp4 . (17) Ms H Which form of the torque applies (equation (16) or (17)) to a particular system depends on the importance of gas pressure. We expect equation (16) to be applicable in the case that H < RH and equation (17) to be applicable otherwise. For typically expected conditions in gaseous protostellar disks, it turns out that equation (17) is the relevant one for planets undergoing (nongap) Type I migration. The migration rate T /Jp with planet angular momentum Jp is then linear in planet mass, since T is quadratic while Jp is linear in planet mass. Therefore, the Type I migration rate increases with planet mass, as asserted in Section 2.1. This somewhat surprising result that more massive planets migrate faster is in turn a consequence of quadratic variation of ∆J with planet mass in equation (7). This quadratic dependence occurs because the possible linear dependence of ∆J on planet mass vanishes due to the antisymmetry of the torque as a function of time along the unperturbed particle path, as discussed in the second paragraph following equation (7). Based on equation (17) with typical parameters for the minimum mass solar nebula at rp = 5AU (Σ = 150 g/cm2, Ω = 1.8 × 10−8 s−1 , and H = 0.05rp ), we estimate the planetary migration timescale Jp /T for a planetary core of mass Mp = 10M⊕ embedded in a minimum mass solar nebula at 5AU as 4 × 105 y. This timescale is short compared to the disk lifetime, estimated as several 106 y or the Jupiter formation timescale of & 106 y in the core accretion model. The relative shortness of the migration timescale is a major issue for understanding planet formation. Since the migration is generally found to be inward, as we will see later, the timescale disparity suggests that a planetary core will fall into the central star before it develops into a gas giant planet. Research on planet migration has concentrated on including additional effects such as gas pressure and improving the migration rates by means of both analytic theory and multi-dimensional simulations. A more detailed analysis reveals that the torque in equation (17) does provide a reasonable estimate for migration rates in gaseous disks in the so-called Type I regime in

(13) (14) (15)

Equation (15) must break down for small ∆r, in order to yield a finite result. For values of |∆r| . 3RH , the particles become trapped in closed orbits in the so-called coorbital region (see Fig. 9). We exclude this region from current considerations, since the torque derivation we considered here does not apply in this region. In particular, the 4

which a planet does not open a gap in the disk. The model provides considerable physical insight, but is crude. In fact, it does not specify whether the migration is inward or outward (i.e., whether T is negative or positive). The density asymmetry about r = rp (see equation (12)) typically involves a higher density at smaller r, as in the case of the minimum mass solar nebula. This variation suggests that torques from the inner disk dominate, implying outward migration, as was thought to be the case in early studies. But this conclusion is incorrect. As we will see in Section 2.4, inward migration is typically favored. We have not included the effects of gas pressure in a proper way. A disk with gas pressure propagates density waves launched by the planet. The analysis has only considered effects of material that passes by the planet. In addition, there are effects from material that lies closer to the orbit of the planet. It is trapped in librating orbits of the coorbital region (see Fig. 9). This region can also provide torques. We have also assumed that the disk density is undisturbed by presence of the planet. Feedback effects of the disk disturbances and gaps in the disk can have an important influence on migration. Finally, there are other physical effects such as disk turbulence that should be considered.

purposes of describing waves in this section, we will ignore viscous forces, as well as gas self-gravity. We now consider the linear departures from the unperturbed state consisting of an axisymmetric disk that orbits the central star. The planet is taken to be on a circular orbit about the star. We express these perturbations as complex quantities for which it is implicit that the real part should be taken. We consider a Fourier decomposition of the gravitational potential in angle and time given by X Φ(r, θ, t) = Φm (r) exp [im(θ − Ωp t)], (22) m

where m is a nonnegative integer, Φm is real, and Ωp is the orbital frequency of the planet. Physical quantities are expressed as X Σ(r, θ, t) = Σm (r) exp [i(m(θ − Ωp t)], (23) m

u(r, θ, t) =

X m

v(r, θ, t) =

X

(24)

vm (r)exp[im(θ − Ωp t)],

(25)

pm (r)exp[im(θ − Ωp t)].

(26)

m

p(r, θ, t) = 2.3 Waves in 2D Gas Disks

X m

We assume here that in lowest order the axisymmetric density component Σ0 (r) is equal to the unperturbed disk density. That is, the presence of the planet does not substantially Σ0 . We return to this point in Section 2.7. We also assume that in the lowest order the axisymmetric azimuthal velocity v0 is equal to the Keplerian speed ΩK (r)r. Pressure pm is assumed to involve an imposed axisymmetric locally isothermal (density independent) temperature variation. We then have that

In this section and the next section, we improve the particle disk description of the migration torque by considering a fluid disk, along the lines of Goldreich & Tremaine, 1979, 1980 and Meyer-Vernet & Sicardy, 1987 (see also Shu, 1991). A fluid disk introduces some additional physical effects. Unlike free particle orbits, disk fluid streamlines cannot cross each other. Tidal disturbances in a fluid disk generally result in the propagation of waves, often referred to as density waves. The disk is taken to be thin so that its thickness c ≪ r, (18) H= Ω

pm (r) = c2 (r)Σm (r).

(27)

Such a temperature variation arises in the case of an optically thin disk that is heated by a central star. This temperature distribution corresponds to a gas sound speed described as c(r). For each m > 1, the linearized forms of equations (19) (21) are

where c is the gas sound speed. We model the disk as two-dimensional and utilize a cylindrical coordinate system (r, θ) centered on the star. The disk has a density distribution Σ and velocity field (u, v). The equations for mass conservation, radial, and azimuthal force balance in the inertial frame are ∂Σ ∂(ruΣ) ∂(vΣ) + + = 0, ∂t r∂r r∂θ

um (r)exp[im(θ − Ωp t)],

d(rum Σ0 ) imvm Σ0 + = 0, rdr r

(28)

−iˆ ωum − 2Ωvm = −

dΦm 1 dpm − , Σ0 dr dr

(29)

−iˆ ωvm + 2Bum = −

impm imΦm − , rΣ0 r

(30)

−iˆ ω Σm +

(19)

∂u v ∂u v 2 1 ∂p ∂Φ ∂u +u + − =− + fvr − , (20) ∂t ∂r r ∂θ r Σ ∂r ∂r and

and

∂v ∂v v ∂v uv 1 ∂p ∂Φ +u + + =− + fvθ − , (21) ∂t ∂r r ∂θ r rΣ ∂θ r∂θ

where ω ˆ = m(Ωp − Ω) is the Doppler shifted frequency, B(r) = Ω+1/2r(dΩ/dr). For a Keplerian disk, B = Ω/4. We seek solutions to equations (28) - (30) by applying the WKB or tight-winding approximation. The approximation is that perturbed quantities vary much more rapidly in

where Φ is the gravitational potential due to the star and planet, p is the gas pressure, and fv is the viscous force per unit disk mass to model the effects of turbulence. For the 5

(OLR). The ILR occurs for r < rp and the OLR for r > rp . For a Keplerian disk, the turning point locations rm are determined by mΩp (37) ΩK (rm ) = m∓1 or 2/3  m∓1 . (38) rm = rp m

the radial direction than in the azimuthal direction. We regard azimuthal wavenumber m to be of order unity. The approximation is that |dym /dr| ≫ m |ym |/r, where ym is Σm , um , vm , or pm . From eqs (29) and (30) it follows that |um | ∼ |vm |. Under these approximations, equations (28) (30) simplify in lowest order to −iˆ ω Σ m + Σ0

dum = 0, dr

(31)

Fig. 6 plots rm as a function of m with lighter dots. For large m, the resonances approach the orbital radius of the planet, |rm − rp | ∼ 2rp /(3m). The sign of kr determines the direction of wave propagation. The radial group velocity of these waves is given by

c2 dΣm dΦm 2mΩΦm i(κ2 − ω ˆ 2) um = − − + , (32) ω ˆ Σ0 dr dr ω ˆr where equation (32) resulted from combining equations (29) and (30). If we ignore the driving terms involving Φm in the above equations, we determine the behavior of free waves having frequency mΩp . We express Z r ym (r) = ym exp [i kr (r′ )dr′ ], (33)

cg

=

(34)

where κ is the epicyclic oscillation frequency defined by κ2 = (1/r3 )(dr4 Ω2 /dr) = 4BΩ. For a Keplerian disk, κ = Ω. The dispersion relation determines radial wavenumber kr (r) as a function of frequency ω ˆ (r). The dispersion relation describes sound waves that are modified by the presence of rotation and shear. Dispersion relation (34) shows that the waves are propagating (kr is real) for |ˆ ω | > κ and evanescent otherwise. Equation (34) tells us that kr ∼ Ω/c ∼ H −1 for m of order unity.2 Therefore, the WKB approximation, kr rp ≫ m, then requires that H ≪ rp , as stated in equation (18). That is, the disk is thin and the gas is cold. For a given m, each radius where kr (rm ) = 0 (35)

dTd (r) dr

= =

(40)

−r iπr

Z

2π

Σ(r, θ)

0 X

∂Φ dθ, ∂θ

m Σm (r) Φm (r),

(41) (42)

m

=

is a radial wave turning point. The wave turning points coincide the locations of resonances, called Lindblad resonances. The turning points (or resonances) occur for ω ˆ = ∓κ,

(39)

Since ω ˆ is negative for an ILR and positive for an OLR, it then follows that for kr > 0, the waves propagate away from the orbit of the planet (cg < 0 at an ILR and cg > 0 at an OLR). Such waves also have the property that they are ”trailing”. That is, their wave fronts bend towards negative θ with increasing r (see Fig. 5). Along a path of constant wave phase kr ∆r + m∆θ = 0, so dθ/dr = −kr r/m < 0 for kr > 0. Leading waves have the opposite properties: they have kr < 0 and propagate towards the orbit of the planet. Up to this point we have disregarded the effects of the driving by the potential Φm . The key physical point is that the propagating waves generally do not interact with the potential because their wavelengths are short compared to the radial variations of the potential. To see this, consider the tidal torque per unit disk radius on the disk

where ym is um or Σm , kr (r) is the radial wavenumber, and ym on the right-hand side is taken to be constant. We assume |kr |r ≫ m (tight-winding approximation) and apply equation (33) in equations (31) and (32) with Φm = 0. We then obtain the following dispersion relation ω ˆ 2 − κ2 = kr2 c2 ,

dˆ ω , dkr kr c2 . ω ˆ

=

X dTm m

dr

(r),

(43)

where we applied equations (22) and (23) in going from equation (41) to (42) and used the fact that Φm is real. It is implicit that the real part of equation (42) should be taken. Each term in the sum involves the product Σm Φm . But in the wave propagating region, Σm (r) is rapidly varying in radius for a cold disk, since |kr |r ≫ 1, while Φm (r) is nonoscillatory and varies more slowly in r. Consequently, over a finite radial interval δr ≫ 2π/|kr | in the wave propagatR r+δr ing region, the net torque r (dTm /dr′ ) dr′ is in some sense small. However, near the Lindblad resonance, where kr is small (see eq (35)), the radial oscillations are not rapid and this integral develops a more substantial contribution.

(36)

where we follow the notation throughout that the upper sign denotes the condition for the inner Lindblad resonance (ILR) and the lower sign for the outer Lindblad resonance 2 Since

the radial wavelength is comparable to the disk thickness, 3D effects may be important. If the waves have a locally isothermal equation of state, the 2D assumption is valid. But if not, e.g., they are adiabatic or if the disk temperature varies with height from the orbit plane, then the waves have a 3D character (Lubow & Pringle, 1993; Lubow & Ogilvie, 1998). In such cases, the wave propagation and damping properties are quite different.

6

xAi(x) = 0 and Airy function Gi(x) satisfies Gi′′ (x) − xGi(x) = 1/π. These functions are described in Abramowicz & Stegun, 1972.

We then expect that torques develop near r = rm that we evaluate below. We determine the solution for equations (28) - (30) in the vicinity of the resonance located at radius rm , where equation (35) is satisfied, by applying a Taylor expansion about r = rm κ2 − ω ˆ 2 ≃ Dm x, (44) 2

2.3 Lindblad Torques To determine the torque on the disk for a given m, we apply the equation of mass conservation (28) and eliminate vm by equation (48), to express Σm near the resonance in terms of um as   i d(rum Σ0 ) mum κΣ0 . (53) ∓ Σm (r) = − ω ˆ rdr 2Ωr

2

where x = (r − rm )/rm and Dm = rd(κ − ω ˆ )/dr evaluated at r = rm . For a Keplerian disk, using Ω = κ and equation (37), we have Dm = ±3mΩ(rm )Ωp = ±3m2 Ω2p /(m ∓ 1). We reconsider equations (31) and (32) to find a solution for driven waves that this valid near a Lindblad resonance. We eliminate Σm in favor of um , again apply the approximation |dΣm /dr| ≫ m|Σm |/r together with equation (44) to obtain that 2

We then obtain from equations (41) - (43) that Z dTm Tm = (r)dr, dr Z = iπm rΣm (r)Φm (r)dr, Z Σ0 Ψ m um dr, = πm κ

2

c d um iκ(rm )Ψm (rm ) − um Dm x = , 2 dx2 rm rm

(45)

where

dΦm 2mΩ Φm ± r . (46) κ dr The left-hand side of equation (45) can readily be seen to be equivalent to the dispersion relation equation (34) by recognizing d/dx → ikr rm and Dm x → κ2 − ω ˆ 2 . The righthand side introduces the forcing. In addition, from equation (29) and (30), it follows from solving for um and vm near the resonance that Ψm =

vm

2iBum κ iκum . = ± 2Ω = ±

d log um (r) >0 dr

(47) (48)

Tm

w3 = and

c2 2 |D | rm m

Cm = −

πκΨm , rm Dm w

(56)

= −

π 2 mΣΨ2m . Dm

(58)

Fig. 7 shows that the torque contributions to the integral in equation (57) come from a radial region of order wrp ≪ rp . There is no dissipation in this model. One might expect the system to be time-reversible, so that the net torque on the planet is zero. But, we have imposed an arrow of time by the causality condition (49) that waves are only emitted at the resonance, implying that only trailing waves are present. In a dissipationless system of finite size, the trailing waves would reflect at the disk center and outer edge as leading waves that propagate back to the resonance (Fig. 5). These waves would cancel the torque from the trailing waves, resulting in a zero net torque. Consequently, the model implicitly assumes some dissipation somewhere in the disk prevents leading waves from returning to the resonance. Dissipation is needed for gas to produce a torque, as it was needed in the case of ballistic particles in Section 2.2. For a given m value there are ILR and OLR contributions from the gas interior and exterior to the orbit of the planet, respectively. We denote these contributions with subscripts

(49)

in the wave propagating region away from the resonance. The solution with these properties is  x  x um = Cm [Ai ± ± i Gi ± ], (50) w w with

(55)

where the integrals extend over all r. Equation (56) is obtained from equation (55) by applying the expression for Σm given by equation (53), integrating by parts, and dropping the surface term. We used the fact that ω ˆ 2 = κ2 near the resonance and applied Ψm given by equation (46). If we assume Σ0 , Ψm , Ω, and κ are nearly constant in the region near the resonance, we obtain from equations (50) and (56) that Z π 2 mΣΨ2m ∞ Tm = − Ai(±x/w)dx, (57) Dm w −∞

We seek solutions with the property that the waves propagate away from and not towards a resonance. That is, we impose a causality condition that waves are emitted from a resonance. These waves behave as trailing waves far from the resonance, as we saw in equation (40). That is, we require kr (r) = −i

(54)

(51)

(52)

where again the upper (lower) signs refer to inner (outer) Lindblad resonances. Airy function Ai(x) satisfies Ai′′ (x)− 7

not get arbitrarily close to the planet with increasing m, as indicated in Fig. 6. As a consequence, Ψm does not diverge as m goes to infinity, as seen in Fig. 8. Instead, it decays exponentially for m & rp /c and the total torque T is finite. 3D effects due to the disk thickness produce a similar drop-off. This means that m is effectively limited to

in and out. The torque on the planet T is minus the sum of the torques on the disk X T = −Td = − (Tm,in + Tm,out ). (59) m

Tm,in is Tm evaluated at the ILR and similarly Tm,out is Tm evaluated at the OLR. Since Dm > 0 at the ILRs, the torques given by equation (58) on the planet due to the inner disk, −Tm,in, are positive, i.e, outward. The direction changes to inward for the outer disk containing OLRs, since Dm < 0. This behavior is the same as we found for particle disks. Notice that the resonant torque T has a similar form to the particle torques in equations (10) and (11). Both depend on the square of the planet mass Mp (equation (58) involves Ψ2m that depends on Mp2 ). For a Keplerian disk, it can be shown that the correspondence between the particle torque and the gas torque on the planet is nearly exact. Notice that gas torque is independent of the gas sound speed in the Keplerian case, since Tm in equation (58) has that property. However, the resonance description for the gas disk torques is more precise than the particle description and allows us to include the nonKeplerian effects due to gas pressure. In this way, the sign and magnitude of the net torque can be determined. For large m, the resonances lie close to the planet (light dots in Figs. 6 and 8) and the torque T becomes singular. The singularity is mainly due to the Ψ2m factor in equation (58) (m/Dm is nearly independent of m). This singularity is related to the ones we encountered for Tout and Tin in equations (10) and (11), since m ∼ 1/∆r. Since the physical torque cannot be infinite, something must be wrong with our description. One issue is that we have approximated the disk as being two-dimensional. A 3D disk would have thickness H = c/Ω. Within a radial distance ∆r = r − rp of order H, we expect this description to break down, since the gas will be spread away from the planet in the direction perpendicular to the orbit plane. There is, however, another effect that yields a finite torque in the large m limit for a 2D disk. The dispersion relation (34) can be shown to have corrections due to azimuthal wave propagation. This ”extended” WKB approximation involves the replacement of kr2 → kr2 + kθ2 , where kθ2 = m2 c2 /r2 (Artymowicz, 1993a). We then have that ω ˆ 2 − κ2 = kr2 c2 +

m 2 c2 . r2

m . mcr =

c ∼ ∓H, rp

(62)

This critical value mcr is referred to as the torque cutoff. The m limit implies that the contributing resonances come from a region no closer to the planet than |∆r| ∼ H, as we applied in our torque estimate equation (17). 2.4 Differential Lindblad Torques In this section, we describe the net torque on the planet that results from the effects of the opposing torque contributions from the ILRs and OLRs, along the lines of Ward, 1997. A more rigorous and complete way to determine the torque involves considering the equations of motion in 3D and improving the validity of the equations near r = rp . In this section, we will describe the main qualitative features of the torque that remain valid when more accurate treatments are applied. We now describe the main reason why planets migrate inward. At the same distance d > 0 inside and outside the orbit of the planet, quantity Ψ2m (rp + d) outside the planet orbit is slightly larger than Ψ2m (rp − d) inside. This asymmetry is due to the circular geometry of the orbit. Quantity Ψm , defined in equation (46), contains a linear combination of radial and azimuthal derivatives of the gravitational potential Φm that are multiplied by radius r. It is the multiplication by r that causes Ψ2m (rp + d) > Ψ2m (rp − d). In addition, the resonance radius rm at the OLR is slightly closer to the planet orbital radius rp than is rm at the ILR. This effect further strengthens torques at OLRs. It can also −1 favor the Tm,out over be seen that the effects of factor Dm Tm,in . As a result, inward planet migration is the suggested outcome. These geometric effects are called curvature effects. They are of course not the only influences on the torque. In the remainder of this section, we will more fully describe various effects on migration. We consider a simplified disk model in which c(r) is taken to be constant in radius and the unperturbed density distribution is taken to be a power law in radius with constant β defined by

(60)

Recall that the site of wave excitation occurs where the wavelength is relatively large, equation (35). For large m, the wave turning point (where kr (rm ) = 0) in equation (60) becomes ω ˆ = ∓mc/rp , since |ˆ ω| ≫ κ. For a Keplerian disk, we have ∆rm = rm − rp ∼ ∓

Ωp rp r = . c H

Σ0 (r) ∝ r−β .

(63)

We consider some contributions to the gas disk rotation rate Ω2 = Ω2K + Ω2pr + Ω2pg ,

(61)

(64)

which are respectively due to the gravity of the star (Keplerian rotation), the axisymmetric component of the disk pressure force, and the axisymmetric component of planet

for large m. So unlike the case of the standard WKB approximation (equations (34) and (36)), the resonances do 8

gravity. The effects of disk self-gravity are ignored. The planet and disk are subject to nearly the same axisymmetric gravitational forces at the location of the planet. But, the disk feels the gravitational force of the planet that is not felt by the planet (no self-force). The disk feels pressure, while the planet does not. The angular frequency of the planet is the Keplerian rate Ωp = ΩK (rp ).

and

Ω2K (rm ) Ω2pg (rm ) Ω2pr (rm ) κ(r)

= rp + ǫrm,1 + ǫ2 rm,2 + . . . , = Ω2K + ǫ

dΩ2K

rm,1 + . . . , dr q r 3 p = ǫ2 Ω2 + . . . , 2πrm,1 p = −ǫ2 β Ω2p + . . . , = ΩK + ǫκ1 + . . . ,

(65)

(66) (67) (68) (69) (70)

where all functions of r on the right-hand side of the above are evaluated at the planet’s orbital radius rp . The scaled planet to star mass ratio q3 = Mp /(ǫ3 Ms ) is assumed to be of order unity. Substituting expansions (66) - (70) into equation (60), we obtain 2p rm,1 1 + (m/mcr )2 , (71) = ∓ rp 3 rm,2 rp

1 β m2 = − − − ∓ γ, 9 9m2cr 3

(74)

where f1 and f2 are dimensionless, order unity, positive functions of |rm,1 |/rp that contains a curvature term. Quantity ∆rm,2 is the mean shift of resonance locations for the ILR and OLR. The first term, -1/9, is a curvature term, while the next two terms involve pressure. Equation (73) shows that the first order antisymmetric shift rm,1 provides a net torque at the same level as the second order symmetric (mean) shift, ∆rm,2 . Since typically β > 0, we have that ∆rm,2 < 0 and both terms on the right-hand side of equation (73) contribute to inward migration. The effect of a negative ∆rm,2 can be understood as moving the OLR closer to the planet and the ILR further away. These shifts in turn cause the OLR torque to dominate and provide a negative torque contribution (see equation (10)). The effects of Ωpg are contained in the γ term. But since γ does not contribute to the mean shift ∆rm,2 , the planet gravity does not influence migration, provided that Mp ∼ ǫ3 Ms . We now consider the contribution to the differential torque due to the density factor that appears in equation (58). For a given m value, the net torque depends on the density change between the ILR and OLR that is related to β. This torque can be seen to be

The torque sum T in equation (59) is dominated by contributions having m ∼ mcr = Ωp rp /c. We consider an expansion of resonance radius rm and orbital frequency squared Ω2 (rm ) in powers of ǫ = 1/mcr ≪ 1 of the form rm

m2 β 1 − , ∆rm,2 = − − 9 9m2cr 3

Tm,Σ = CΣ

β |rm,1 | |Tm,in |, mcr rp

(75)

where constant CΣ is positive. For β > 0 as is typically the case, torque Tm,Σ is positive. It can be shown that torque Tm,Σ is be comparable to and opposite to Tm,cpg , but somewhat weaker. As a result, the differential (net) Lindblad torque is typically negative. Notice that by summing equations (73) plus (75) over m, we have that the differential torque X T ∼ Tm,cpg (76)

(72)

where quantity γ describes radius changes that are equal and opposite about r = rp . We see that the lowest order radius shifts rm,1 are equal and opposite. As we described above, the curvature asymmetry favors inward migration for resonances located at equal distances inside and outside radius rp . The torque contributions due to these resonance position and angular velocity shifts are obtained by applying equations (58) and (59).3 We first consider the differential torque contribution that does not involve the density factor Σ in equation (58). We denote this differential torque as Tm,cpg , since it involves curvature, pressure, and planet gravity. After considerable algebra, it follows that Tm,cpg due to the ILR and OLR for some 1 ≪ m . mcr is given by  2 Mp 4 2 Tm,cpg = −mcr Σrp Ωp Ms       |rm,1 | ∆rm,2 |rm,1 | + f2 (73) f1 rp rp rp

m

∼ mcr Tmcr ,cpg

(77)



2

Mp Ms  2  r 2 Mp p 4 2 Σrp Ωp , ∼ − H Ms

∼ −m2cr Σrp4 Ω2p

(78) (79)

which agrees with equation (17). A more careful evaluation of the sum can be performed to provide the detailed dependence of T on β. We will reconsider this issue in Section 2.6. 2.5 Coorbital Torques Thus far we have considered torques that arise from gas that passes the planet in the azimuthal direction (the y direction in Fig. 1). These torques can be described by Lindblad resonances, as we saw in Section 2.4. Gas that resides

3 The

derivation of equation (58) assumed that ω ˆ 2 = κ2 at the Lindlabd resonances. There are corrections to this condition and therefore the torque in the extended WKB approximation (equation (60)). For simplicity, we ignore such effects here.

9

Σ(ro ) at position o+. The ratio is given by

closer to the planet, |r − rp | . 3RH , in the so-called coorbital region, does not pass by the planet. Instead, it follows librating orbits in the corotating frame of the planet as seen in Fig. 9. We consider orbits that nearly fully circulate, the so-called horseshoe orbits. They approach the planet at both θ < θp and θ > θp , as seen in the figure. As the gas approaches the planet, it gets pulled by the planet from r > rp to r < rp or vice versa. This change in angular momentum of the gas in turn causes a torque to be exerted on the planet. Away from the region where the orbit transitions between different radii, the gas follows an approximately circular Keplerian path. For a disk with no dissipation and a planet on a fixed circular orbit, the net torque on the planet is zero. The reason is that the gas follows periodic orbits in the frame of the planet, so there is no change in angular momentum of the gas over a complete period of its motion. Whatever angular momentum is gained by a gas element as it changes from the inner radius ri to outer radius ro is lost when it later encounters the planet and shifts from ro to ri . So although gas on horseshoe orbits can come close to the planet, the symmetry limits the torque that it exerts on the planet. To analyze this situation more carefully, consider the labeled regions in Fig. 9. We describe the disk as consisting of a gas fluid, but will ignore the effects of gas pressure, and follow the approach of Ward, 1991. Gas loses angular momentum in going from position o+ to i+. This angular momentum is continuously gained by the planet. Similarly, the planet continuously loses angular momentum from the gas that passes from i− to o− . The gas spends a relatively short time in transition between ri and ro compared to the time it spends between encounters that is of order the synodic timescale ∼ rp /((ro − rp )Ωp ). Between encounters, the gas may experience the effects of turbulent viscosity. It acts to establish a characteristic density profile Σ0 (r) that would be present in the absence of the planet. As the gas re-encounters the planet, it does so with the turbulenceenforced background density. Such effects give rise to a net torque on the planet. The net torque on the planet due to a streamtube that extends on both the θ < θp and θ > θp sides of the planet is given by δTco = (M˙ + − M˙ − )

dr2 Ω w, dr

M˙ + Σ(ro ) B(ri ) ∼ , ˙ Σ(ri ) B(ro ) M−

where the ratio of B values comes into play because of a change in the area of fluid elements between ri and ro . We then obtain an expression for the torque on the planet due to the streamtube as ! ˙+ M − 1 Ωp rp w, (83) δTco ∼ M˙ M˙ −   ∆Σ ∆B 2 2 δr, (84) − ∼ ΣΩp rp w Σ B d log (Σ/B) ∼ ΣΩ2p w3 δr, (85) d log r where ∆Σ = Σo − Σi and ∆B = Bo − Bi . Integrating over all streamtubes to width w, we estimate the coorbital torque as d log (B/Σ) Tco ∼ −ΣΩ2p w4 , (86) d log r where we now interpret w as the width of the coorbital region that is a few times RH . Quantity B/Σ is sometimes called the vortensity, since B is related to the vorticity or curl of rΩeθ . The coorbital torque then depends on the gradient of the vortensity. A time-reversible system (no dissipation), with arbitary initial conditions would evolve towards a state in which the vortensity is constant in the coorbital region. This process occurs through ”phase mixing” on the libration timescale ∼ rp /(Ωp w) (see Fig. 10). The coorbital torque drops to zero or is said to be saturated. As discussed above, maintaining a density gradient and a nonzero torque could be accomplished by the effects of turbulent viscosity. We compare the coorbital torque with the net Lindblad torque for a particle (pressureless) disk as we discussed in Section 2.2. Taking w ∼ RH in equation (86), we see that the unsaturated coorbital torque (taking |d log (B/Σ)/d log r| ∼ 1 and nonzero) is comparable to the differential Lindblad torque in equation (16). It can be shown that this is also true for a disk where pressure effects are important, where H > RH . For such a disk, the coorbital torque is generally of the same order as the differential Lindblad torque, equation (17). The direction of the coorbital torque contribution depends on the sign of the vortensity gradient. For the minimum mass solar nebula model, Σ ∝ r−3/2 and B ∝ r−3/2 . Consequently, the coorbital torque is zero in that case. For smaller values of β < 3/2, the torque is positive. For the turbulent viscosity to prevent torque saturation, the viscous timescale across the coorbital region should be shorter than the libration timescale. This constraint can be translated into a condition on the minumum magnitude of the turbulent viscosity required for a nonzero corotation torque (see equation (89)). With pressure effects, the gas communicates disturbances with neighboring gas. But the gas is incapable of

(80)

where w = ro − ri is the radius change along the streamtube, M˙ − is the flux of mass passing the planet for θ < θp and analogously for M˙ + . Both mass fluxes are defined to be positive and are approximately given by M˙ ∼ M˙ ± ∼ Σrp (Ω(ri ) − Ωp )δr,

(82)

(81)

where δr is the radial width of the streamtube away from the planet, on the circular portion of the orbit. The ratio of the mass fluxes depends on the density Σ just ahead of the encounter, with density Σ(ri ) at position i− and density

10

launching a propagating wave because the region near rp is evanescent, as seen in Fig. 5. Instead, the acoustic disturbances remain trapped within a radial distance H of the planet’s orbit. The trapping of disturbances near the corotation radius rp limits the amount of angular momentum that a planet can gain. But in the case of a viscous disk, mass and angular momentum in the coorbital region can be transferred to the remainder of the disk by viscous torques. The corotation torque can then act without saturation. Although the coorbital torque is of the same order as the differential Lindblad torque, more detailed calculations as described in the next section show it typically does not dominate over Lindblad torques.

where α the dimensionless disk viscosity parameter in the standard α disk model (Ward, 1992). For disk parameters α = 0.004 and H/r = 0.05, this constraint implies that for planets of order 10 M⊕ or greater, the corotation torques should be saturated (small) and equation (87) should be applied. Nonlinear 3D hydrodynamical calculations have been carried out to test the migration rates, under similar disk conditions used to derive the analytic model. Figs. 11 and 12 show that the migration rates agree well with the expectations of the theory. We examine the comparison between simulations and theory in more detail by comparing torque distributions in the disk as a function of disk radius. We define the distribution of torque on the planet per unit disk mass as a function of radius as dT /dM (r) = 1/(2πrΣ(r)) dT /dr(r). Fig. 13 plots (Ms /Mp )2 dT /dM as a function of radial distance from the planet based on 3D simulations. The distributions show that the torque from the region interior (exterior) to the planet provides and positive (negative) torque, as predicted in equations (11) and (10) for a particle disk and as shown in Section 2.3 for a gas disk. Also, the integrated total torque is negative, as expected. The theory predicts that the torque density peak and trough occur at distance from the planet r − rp = ∆r ∼ ∓H, where the torque cut-off takes effect as indicated in equation (61). For the case plotted in Fig. 13 that adopts H = 0.05r, the predicted locations agrees well with the locations of the peaks and troughs in the figure. The theory in Section 2.4 and 2.5 also predicts that for a fixed gas sound speed, the shape of the scaled torque density distribution (Ms /Mp )2 dT /dM (r) is independent of planet mass. The reason is that the width of each contributing resonance is independent of planet mass, but depends on the sound speed (see equation (51)). In addition, the set of of contributing resonances (range of m values) is also independent of planet mass (see equation (62)). Furthermore, the overall torque scales with the square of the planet mass. As seen in Fig. 13, these expectations are well met for the two cases plotted.

2.6 Type I Migration Torques Thus far, our estimates of migration rates have assumed that the axisymmetric gas surface density is not perturbed by the presence of the planet. This regime is sometimes called the Type I case of migration. More detailed 3D linear analytic calculations of the Type I migration rates have been carried out by Tanaka et al, 2002. They assumed that the gas sound speed is constant in radius. For the case of saturated (zero) coorbital corotation torques, where only differential Lindblad torques are involved, the migration rate is given by 2    rp 2 Mp T = − (2.34 − 0.10 β) Σp Ω2p rp4 , (87) Ms H where β is given by equation (63). The torque on the planet resulting from the action of both Lindblad and (unsaturated) coorbital corotation torques is given by  2   Mp rp 2 2 4 . (88) T = − (1.36 + 0.54 β) Σp Ωp rp Ms H These migration rates are consistent with the estimate in equation (79). Notice that the Lindblad-only torque, equation (87), contains a small coefficient of β. The reason is that there is a near cancellation of the inward migration effects of the pressure in equation (73) (in ∆rm,2 ) with the outward migration effects of the density in equation (75). The differential Lindblad torque is then mainly due to the curvature effects. In going from equation (87) to equation (88), the effects of the coorbital torque reduce the inward migration rate for β < 1.5 and nearly vanish for β = 1.5.4 This behavior is consistent with equation (86) for which the coorbital torque is positive for β < 1.5 and is zero for β = 1.5. Based on scaling arguments, the condition on the turbulent viscosity for the coorbital torque to be effective (unsaturated), as discussed in Section 2.5, is given by α& 4 The



Mp Ms

3/2 

r 7/2 . H

2.7 Disk Response The various torque expressions (e.g., equations (17), (73), and (88)), contain a density term. This density is the lowest order density at the radius of the planet. Up to this point, we have assumed that this density is the unperturbed disk density, the disk density that occurs in the absence of a planet. If the planet mass is sufficiently small or the level of turbulent viscosity is sufficiently large, then the presence of the planet does not substantially modify the density distribution. But these conditions may be violated and the density can be affected. In particular, some our approximations, such as in going from equation (56) to equation (58), were reliant on setting the lowest order density distribution Σ0 to the unperturbed density distribution. Modifications to Σ0 caused by the planet need not be small and can produce a significant

(89)

small nonzero coorbital torque at β = 1.5 is due to 3D effects.

11

change in the migration torque. To analyze the disk response, we consider equations (19) and (21) and include the effects of shear viscosity with kinematic viscosity ν = αcH by taking ∂G , ∂r dΩ G = −2πr3 νΣ , dr fvθ = −

reminiscent of the case of ocean waves. They are generated by wind far from land, but undergo final decay when they break at the shore. Notice that the planet obtains its torque from gas near the Lindblad resonance, independent of where the waves damp and the disk gets torqued. Recall that it is necessary that the waves damp somewhere for there to be a nonzero Lindblad migration torque. c) Disk with strong dissipation and tidal forces In this case, the waves damp immediately by the strong dissipation and therefore Fw = 0. We have seen that diskplanet torques from an outer disk cause an inward torque on the planet and an outward torque on the disk. Similarly, the interactions cause an inward torque on the inner disk. As a result of these torques, material is pushed away from the orbit of the planet. In the case that the planet produces strong tidal force, we may expect that a gap is created near the orbit of the planet. Consequently, we expect that near the planet M˙ = 0 and Td = G in equation (93).

(90) (91)

and we will ignore fvr . We write the azimuthal velocity as v = Ω(r)r + v ′ ,

(92)

where velocity v ′ denotes the nonaxisymmetric departures from circular rotation that are assumed small compared to the axisymmetric circular velocity Ωr. Multiplying equation (21) by rΣ, applying equation (19), and integrating in θ, we obtain to high accuracy the torque equation ∂Td ∂Fw ∂G dr2 Ω = − − , M˙ dr ∂r ∂r ∂r

2.8 Type II migration (93)

Type II migration occurs for case c above when the region near the planet is strongly depleted of gas and a gap forms (Lin & Papaloizou, 1986). Gap formation occurs for Td = G in equation (93). Tidal forces on the disk interior or exterior to the planet are estimated by equation (10) with ∆r ∼ H, as we argued at the end of Section 2.3. For a given level of disk turbulent viscosity, the gap opening condition Td = G becomes a constraint on the planet mass. The condition is estimated as 1/2  3/2  H 40ν Mp & . (98) 2 Ms rp Ω rp

with M˙ = r ∂Td = −r ∂r

Z

2π

Z

0 2π

Σ

0

Fw = r2 Σ0

Z

Σudθ,

(94)

∂Φ dθ, ∂θ

(95)

uv ′ dθ,

(96)

2π

0

where M˙ is the mass flux through the disk, dTd /dr is the torque per unit radius on the disk, and Fw is the flux of angular momentum carried by waves. By integrating equation (19) in θ we have that ∂Σ0 1 ∂ M˙ =− . ∂t r ∂r

For disk parameters α = 0.004 and H/r = 0.05, the predicted gap opening at the orbit of Jupiter occurs for planets having a mass Mp & 0.2MJ . This prediction is in good agreement with the results of 3D numerical simulations (see Fig. 14). In addition to the above viscous condition, an auxiliary condition for gap opening has been suggested based on the stability of a gap. This condition is to preclude gaps for which steep density gradients would cause an instability that prevents gap opening. This condition, called the thermal condition, is given by rH & H (Lin & Papaloizou, 1986). The critical mass for gap opening by this condition is given by  3 Mp H &3 . (99) Ms rp

(97)

We briefly consider some specific cases of interest below. a) Unperturbed viscous disk Consider the case of a steady state viscous disk in which no planet is present, Fw = Td = 0. For a steady state, equation (97) requires that M˙ is constant. From equation (93), we have that M˙ Ωr2 = −G. For a Keplerian disk, we then recover the standard result that M˙ = −3πν(r)Σ0 (r). b) Disk with planet and little dissipation Consider the case that there is no dissipation in some region of space where waves are generated (i.e., near the resonance). In that region M˙ = 0, since there is reversibility. By equation (93) with G = 0 we have that Td = Fw . That is, the disk torque exerted near the planet is transferred by the wave flux Fw to some other region of the disk. Wherever the wave damps, there is an imbalance between Td and Fw , resulting in a nonzero mass flux M˙ that will generally result in a density change (equation (97)). This situation is

For H/r = 0.05, this condition requires a larger planet mass for gap opening than equation (98) for α . 0.01. Even if both the thermal and viscous conditions are satisfied, a substantial gas flow, M˙ ∼ −3πν(r)Σ0 (r), may occur in the presence of a gap in certain circumstances (Artymowicz & Lubow 1996).5 5 Near

12

˙ = 0, since the density the end of Section 2.7 we asserted that M

The migration rate of a planet embedded in a gap is quite different from the Type I (nongap) case that we have already considered. A planet that opens a gap in a massive disk, a disk whose mass is much greater than the planet’s mass, would be expected to move inward, pushed along with the disk accretion inflow. The planet simply communicates the viscous torques across the gap by means of tidal torques that balance them. The Type II migration timescale is then of order the disk viscous timescale  r 2 1 rp2 rp2 p , (100) ∼ ∼ tvis ∼ ν αcH H αΩp

shorter for a disk having a mass greater than the minimum mass solar nebula. For retention of habitable planets and giant planet cores, the short timescales for Type I migration is a serious problem. To be consistent with the ubiquity of extrasolar gas giants and formation of Jupiter and Saturn, some studies suggest that the Type I migration rates must be reduced by more than a factor of 10 (Alibert et al. 2005; Ida & Lin 2008). One major question is whether there are processes that could slow the migration. Several ideas have been proposed. The major uncertainty with them concerns our knowledge of the true state of the disk. Many young stars, such as the T Tauri stars, are surrounded by gaseous disks and have observational signatures of gas accretion. How the accretion operates on a global scale is not known. There are observational signatures of accretion onto T Tauri stars. Some form of turbulence is likely needed to produce the observationally inferred accretion rates. The nature of disk turbulent viscosity has an influence on the global disk structure and the disk response to the presence of a planet. Disk turbulence within planetforming regions of the disks in early stages of evolution is likely to be dominated by gravitational instability and later by magnetic instability. In the latter case, a major unknown is the level of disk ionization. Unless the gas is sufficiently ionized, magnetic instability will not occur, e.g., Salmeron & Wardle, 2008. We briefly discuss below a few of the several suggested mechanisms for slowing Type I migration. Low viscosity regions Disk viscosity suppresses the formation of the weak disk perturbations that are produced by a small mass planet, as suggested by equation (98). The perturbations are smoothed by viscous diffusion. But if the disk turbulent viscosity is sufficiently low, the disk density distribution can be affected by the presence of a low mass planet. We have seen that the ILR and OLR torques push material away from the orbit of the planet. If the planet is on a fixed orbit, the density perturbation is almost symmetric about r = rp . This near symmetry is broken by the migration of the planet. In the comoving frame of the planet, there is a radial steady-state gas flow past the planet. For an inwardly migrating planet, this flow causes a feedback effect that enhances the gas density, interior to the orbit, and lowers the density exterior to the orbit (Ward, 1997). The feedback then enhances the positive torques that arise in the inner disk and slows inward migration. This feedback grows with planet mass. Above some critical planet mass Mcr the steady-state radial flow of gas past the planet is not possible. As a result, the planet migration stops, and gap formation begins. Planet migration can be halted for a certain critical planet mass Mcr that depends on the gas sound speed, the turbulent viscosity, and the rate at which wave damping occurs. Shocks can provide wave damping, although the damping is not instantaneous. The waves launched by a low mass planet propagate some distance as they steepen and ultimately dissipate. The values of Mcr caused by shocks in 2D disks are typically ∼ 10M⊕

which is ∼ 105 years for α = 0.004, H = 0.05rp , and Ωp = 2π/12 y −1 . Therefore, the migration timescale can be much longer than the Type I migration timescale for the higher mass planet that open gaps, as is found in simulations (Mp > 0.1MJ in Figs. 11 and 12). However, this timescale is still shorter than the observationally inferred global disk depletion timescales ∼ 106 − 107 years. The actual migration rate may be somewhat smaller, in order to explain the abundant population of observed extrasolar giant planets beyond 1AU (Ida & Lin 2008). Such retardation is also required to retain terrestrial planets in habitable zones that are located inside the likely formation regions of gas giants that may be beyond the snow lines. In practice, the conditions for pure Type II migration are unlikely to be satisfied. The disk mass may not be very large compared to the planet mass and the disk gap may not be fully clear of material. However, simulations have shown that to within factors of a few, the migration timescale is of order the viscous timescale of the disk over a wide range of parameters, provided the tidal clearing is substantial and the disk mass is at least comparable to the planet mass (e.g., Fig. 12). This appears to be true even if there is a significant mass flow in the gap. The planet in Fig. 12 has mass ∼ 1MJ at a time 1000 orbits and is migrating at the disk viscous rate. At this time, the planet is accreting mass through a deep gap at a substantial rate, comparable to the accretion rate of the disk in the absence of a planet. Another way of understanding the Type II torques is to recognize that the distribution of tidal torques per unit disk mass, as seen in Fig. 13 still applies, even if the disk has a deep gap that is not completely clear of material. The disk density through the gap region adjusts so that the planet migration rate is compatible with the evolution of the diskplanet system. 3. OUTSTANDING QUESTIONS 3.1 Limiting Type I Migration As seen in Fig. 11, the timescales for Type I migration are short compared to disk lifetimes. They become even is low in a tidally produced gap. But, the mass flux need not be small, even if the density is small, provided that the radial velocity u increases sufficiently in the gap, as occurs certain situations.

13

(Rafikov, 2002; Li et al, 2008). Suppose a planet core could form in less than the migration timescale prior to reaching this mass, typically less than a few times 105 y. In this case, the overall migration timescale may be comparable to Type II migration timescale that is also long due to the small values of α (see equation (100)). Type I migration might then nearly stall before the planet migrates through the disk. Analytic calculations and simulations suggest that the feedback is strong, typically for α values up to of a few times 10−4 . For α & 0.001, Type I migration proceeds with little reduction (see Fig. 15). Disk Property Jumps We have seen in Section 2 that the torque on the planet depends on differences in disk properties across rp . We expect the density and temperature to generally smoothly decrease in radius with the outcome being inward migration (see equation (88)). However, this need not always be the case. It is possible that sudden changes in disk properties with radius could occur as a consequence of sudden changes in disk opacity, turbulent viscosity, or gravity variations (due to the planet itself). In such regions, it is possible that inner disk contributions could be enhanced or even be greater than the outer disk contributions. A planet could then experience slowed migration or even be trapped in such a region with no further inward migration (Menou & Goodman, 2004; Masset et al, 2006). Turbulent Fluctuations We have modeled the effects of disk turbulence by means of a turbulent viscosity in equation (91). But, in addition there are time-dependent small-scale density fluctuations that give rise to a fluctuating random torque on the planet. Unlike the Type I torque that acts continuously in the same direction, the fluctuating torque undergoes changes in direction on timescales characteristic of the turbulence that are short compared with the migration timescale. The fluctuating torque causes the planet to undergo a random walk. For the random walk to compete with Type I migration, the amplitude of the fluctuating torque must be much larger than the Type I torque. The change in angular momentum √ of the planet due to a random torque TR is given by ∼ N TR τ , where N is the number of fluctuations felt by the planet and τ is the characteristic timescale for the torque fluctuation. The angular momentum change caused by Type I torques Tin over time t is simply ∼ tTin . Since N = t/τ , it follows that √ for the random migration to dominate, we require TR > N Tin . If we take τ to be the orbital period of the planet and t to be the Type I migration timescale ∼ 105 y, the condition becomes TR > 300Tin. Torque TR depends linearly on the planet mass, while torque Tin increases quadratically with the planet mass and the migration time t decreases with the inverse of the planet mass. The random torque is then more important for lower mass planets. The nature of the random torque depends on the properties of the disk turbulence, in particular its power spectrum, that are generally not well understood. If there is power in the turbulence spectrum at low frequencies, then the fluctuating torques are more effective at counter-

acting the Type I torques. The reason is that the effective N value is smaller. Some simulations and analytic models suggest that turbulent fluctuations arising from a magnetic instability (the magneto-rotational instability Balbus & Hawley, 1991) are important for migration of lower mass planets (Nelson, 2005; Johnson et al., 2006; see Fig. 18). However, if turbulent fluctutations are important for migration, then the eccentricities of planetesimals are pumped up so highly that collisions between them may result in destruction rather than accretion (Ida et al., 2008). Therefore, although the turbulent fluctuations may inhibit the infall of planetary cores into the central star by migration, they tend to inhibit the build up of the cores necessary for giant planet formation in the core accretion model. 3.2 Other Forms of Migration Kozai Migration A planet that orbits a star in a binary star system can periodically undergo a temporary large increase in its orbital eccentricity through the a process known as the Kozai effect. Similar Kozai cycles occur in multi-giant planet systems. The basic idea behind Kozai migration is that the increased eccentricity brings the planet closer to the star where it loses orbital energy through tidal dissipation. In the process, the planet’s semi-major axis is reduced and inward migration occurs (Wu & Murray, 2003). We describe this in more detail below. Consider a planet in a low eccentricity orbit that is well interior to the binary orbit and is initially highly inclined with respect to it. The orbital plane of the planet can be shown to undergo tilt oscillations on timescale of ∼ Pb2 /Pp , where Pb is the binary orbital period, Pp is the planet’s orbital period, and by assumption Pb ≫ Pp . Under such conditions, it can be shown that the component of the planet’s angular momentum perpendicular to the binary orbit plane q (the z-component) is approximately conserved, Jz = Mp GMs ap (1 − e2p ) cos I, where ap and ep are respectively the semi-major axis and eccentricity of the planet’s orbit and I is the inclination of the orbit with respect to the plane of the binary. The conservation of Jz is easily seen in the case that the binary orbit is circular and the companion star is of low mass compared to the mass of the star about which the planet orbits. On such long timescales ≫ Pb , the companion star can be considered to be a continuous ring that provides a static potential. In that case, the azimuthal symmetry of the binary potential guarantees that Jz is conserved. By assumption, we have cos I ≪ 1 and ep ≪ 1 in the initial state of the system. As the planet’s orbital plane evolves and passes into alignment with the binary orbital planet, cos I ∼ 1, conservation of Jz requires ep ∼ 1. In other words, Jz is initially small because of the high inclination of the orbit. When the inclination drops, the orbit must become more eccentric (radial), in order to maintain the same small Jz value. The process then periodically trades high inclination for high eccentricity.

14

During the times of increased eccentricity, the planet may undergo a close encounter with the central star at periastron distance ap (1 − ep ). During the encounter, the tidal dissipation involving the star and planet results in an energy loss in the orbit of the planet and therefore a decrease in ap . This process then results in inward planet migration. The energy loss may occur over several oscillations of the orbit plane. Another requirement for the Kozai process to operate is that the system must be fairly clean of other bodies. The presence of the other object could induce a precession that washes out the Kozai effect. Given the special requirements needed for this process to operate, it is not considered to be the most common form of migration. However, there is good evidence that it does operate in some systems. The Kozai effect can also occur in two-planet systems, where the outer planet plays the role of the binary companion. The process can be robust due to the proximity of outer planet (Nagasawa et al. 2008). Runaway Coorbital Migration In Section 2.5 we saw that the coorbital torque will be saturated (reduced to zero) unless some irreversibility is introduced, such as turbulent viscosity. However, planet migration itself introduces irreversibility and could therefore act to prevent torque saturation. The coorbital torque model presented in Section 2.5 ignored the effects of planet migration on the horseshoe orbits near the planet (Fig. 9). The coorbital torque for a migrating planet could then depend on the rate of migration. Under some conditions, the coorbital torque could in turn cause faster migration and in turn a stronger torque, resulting in an instability and a fast mode of migration (Masset & Papaloizou, 2003). The resulting migration is sometimes referred to as Type III migration. To see how this might operate in more detail, we consider the evolution of gas trapped in the coorbital region (Artymowicz, 2004; Ogilvie & Lubow, 2006). For a sufficiently fast migrating planet, the topology of the streamlines changes with open streamlines flowing past the planet and closed streamlines containing trapped gas (see Fig. 16). The leading side of the planet contains trapped gas acquired at larger radii, while the gas on the trailing side is ambient material at the local disk density. The density contrast between material on the trailing and leading sides of the planet gives rise to a potentially strong torque. A major question centers around the conditions required for this form of migration to be effective. If the planet mass is very small, the process appears ineffective. Once a more massive planet forms a gap, there is an insufficient amount of gas in the coorbital region to cause a substantial torque. Some simulations suggest that this type of migration requires a somewhat massive planet that is not allowed grow in mass to be immersed in a disk and allowed to migrate before gap opening is complete, artificially bypassing the gap opening that would occur naturally for a growing planet (e.g, Zhang et al, 2008; D’Angelo & Lubow 2008). These simulations suggest that this form of migration may not typically arise, due the the special conditions required.

Migration Driven By Nonisothermal Effects in the Coorbital Region Many studies of disk planet interactions simplify the disk temperature structure to be locally isothermal. The locally isothermal assumption, frequently applied in numerical simulations and as we applied in equation (27), means that the temperature structure is prescribed and is independent of disk density. The behavior in the isothermal limit tends to suggest that coorbital torques do not typically dominate migration (e.g., equation (88)). The nonisothermal regime has been recently explored in simulations by Paardekooper & Mellema, 2006 who find that outward migration due to coorbital torques may occur in certain regimes. The conditions required for this possible effect is an active area of investigation. Migration in a Planetesimal Disk After the gaseous disk is cleared from the vicinity of the star, after about 107 y, there remains a disk of solid material in the form of low mass planetesimals. This disk is of much lower mass than the original gaseous disk. But the disk is believed to have caused some migration in the early solar system with important consequences (Hahn & Malhotra, 1999; Tsiganis et al, 2005). There is strong evidence that Neptune migrated outward due to the presence of Kuiper belt objects that are resonantly trapped exterior, but not interior, to its orbit. The detailed dynamics of a planetesimal disk are somewhat different from the case of a gaseous disk, as considered in Section 2. The planetesimals behave as a nearly collisionless system of particles. Jupiter is much more massive than the other planets and can easily absorb angular momentum changes in Neptune. As Neptune scatters planetesimals inward and outward, it undergoes angular momentum changes. It is the presence of Jupiter that breaks the symmetry in Neptune’s angular momentum changes. Once an inward scattered planetesimal reaches the orbit of Jupiter, it gets flung out with considerable energy and does not interact again with Neptune. As a result of the loss of inward scattered particles, Neptune gains angular momentum and migrates outward, while Jupiter loses angular momentum and migrates slightly inward. A similar process occurs in gaseous decretion disks of binary star systems (Pringle, 1991). The circumbinary disk gains angular momentum at the expense of the binary. The binary orbit contracts as the disk outwardly expands. A gapopening planet embedded in a circumbinary disk (or under some conditions, a disk that surrounds a star and massive inner planet) would undergo a form of Type II migration that could carry the planet outward (Martin et al, 2007). In the solar system case, the Sun-Jupiter system plays the role of the binary. Viscous torques are the agent for transferring the angular momentum from the binary outward in the gaseous circumbinary disk, while particle torques play the somewhat analogous role in the planetesimal disk. In a planetesimal disk, another process can operate to cause migration. This process is similar to the runaway coorbital migration (Type III migration) described above, 15

but applied to a collisionless system of particles (Ida et al. 2000). Interactions between the planestimals and the planet in the planet’s coorbital zone can give rise to a migration instability. Migration of Eccentric Orbit Planets The analysis in Section 2 assumed that planets reside on circular orbits. This assumption is not unreasonable, since there are strong damping effects on eccentricity for a planet that does not open a gap in the disk (Artymowicz, 1993b). Some eccentricity may be continuously produced by turbulent fluctuations in the gas, as described above, or by interactions with other planets. In general, eccentricity damping is faster than Type I migration. For planets that open a gap, it is possible that they reside on eccentric orbits in the presence of the gaseous disk. In fact, one model for the observed orbital eccentricities of extra-solar planets attributes the excitation of eccentricities to disk-planet interactions (Goldreich & Sari, 2003; Ogilvie & Lubow, 2003). A planet on a sufficiently eccentric orbit embedded in a circular disk will orbit more slowly at apoastron than the exterior gas with which it tidally interacts. Similarly, a planet can orbit more rapidly than the tidally interacting gas at periastron. These angular velocity differences can change the nature of the ”friction” between the planet and the disk discussed in Section 2.2. For example, at apoastron the more slowly orbiting planet could gain angular momentum from the more rapidly rotating gas. This situation is then just the opposite of the case in Section 2.2. Furthermore, since the planet spends more time at apoastron than periastron, the effects at apoastron could dominate over effects at periastron. It is then possible that outward migration could occur for eccentric orbit planets undergoing Type I migration, assuming such planets could maintain their eccentricities (Papaloizou, 2002). In the case of a planet that opens a gap, simulations suggest that outward torques dominate as the planet gains eccentricity from disk-planet interactions (D’Angelo et al, 2006). If giant planet orbits evolve this way, then their orbital distribution might favor their presence at larger radii, beyond the snowline where they may form (see Chapter 19). The situation is complicated by the fact that the gaseous disk generally gains eccentricity from the planet by a tidal instability (Lubow, 1991). For the outward torque to be effective, there needs to be a sufficient difference in the magnitude and/or orientation between the planet and disk eccentricities, so that the planet moves slower than nearby disk gas at apoastron. Multiplanet Migration Thus far, we have only considered single planet systems. Of the more than 200 planetary systems detected to date by Doppler techniques, over 20 reside in multi-planet system (Butler et al, 2006). About 5 systems have been found to have orbits that lie in mutual resonance, typically the 2:1 resonance. The resonant configurations are likely to be the result of convergent migration, migration in which the separation of the orbital radii decreases in time. This process occurs as the outer planet migrates inward faster than the

inner planet. Planets can become locked into resonant configurations and migrate together, maintaining the planetary orbital frequency ratio of the resonance. The locking can be thought of as a result of trapping the planets within a well of finite depth. Just which resonance the planets become locked into depends on their eccentricities and the relative rate of migration that would occur if they migrated independently. As planets that are initially well-separated come closer together, they lock into the first resonance that provides a deep enough potential to trap them against the effects of their convergence. We discuss below the consequences of resonant migration. To maintain a circular orbit, a migrating planet must experience energy and angular momentum changes that satisfy dE/dt = Ωp (t)dJ/dt, where Ωp is the angular speed of the planet. As the planets migrate together, their mutual interactions cause deviations from this relation. As a result, their energies and angular momenta evolve in a way that is incompatible with maintaining a circular orbit. Orbital eccentricities, as well as mutual inclinations, can develop (Lee & Peale, 2002; Yu & Tremaine, 2001; Thommes & Lissauer, 2003). To see how this process operates in more detail, consider the case that the mass of the inner planet Mi is much less than the mass of the outer planet planet Mo . We assume the planets undergo convergent migration. As the planets migrate together locked in a resonance, the inner planet undergoes energy and angular momentum changes as result of its interaction with the outer planet. The much more massive outer planet is unperturbed by the small mass inner planet. We ignore the effects of the disk interactions on the inner planet compared with the effects of the outer planet. We assume that the outer planet migrates inward due to its interaction with the disk and maintains a circular orbit. We consider a cylindrical coordinate system as in Section 2. The energy of the inner planet is given by Mi vi2 + Mi Φs (ri ) + Mi Φo (ri , ro (t), θo (t) − θi ), 2 (101) where Φs is the potential due to the star. Potential Φo (ri , ro (t), θo (t) − θi ) is due to the outer planet. It contains an explicit time dependence due to the position of the outer planet at (ro (t), θ0 (t)). Taking the time derivative of the energy of the inner planet, we obtain   ∂Φo dvi dEi = Mi vi  + ∇i (Φs + Φo ) + Mi , (102) dt dt ∂t Ei =

where the gradient ∇i is taken in the inner planet coordinates (ri , θi ). The first term on the right-hand side of equation (102) vanishes as a consequence of the equation of motion of the planet. We then have that dEi dt

16

∂Φo , ∂t ∂Φo dro ∂Φo − M i Ωo , = Mi dt ∂ro ∂θi = Mi

(103) (104)

where Ωo = dθo /dt and we have used the fact that Φo depends on θo (t) − θi . We recognize that −Mi ∂Φo ∂θi is the torque on the inner planet that is equal to dJi /dt. The ratio of the first term to the second term in equation (104) is then easily shown to be ∼ Mo /Ms ≪ 1. We then have to high accuracy that dJi dEi = Ωo (t) . (105) dt dt The inner planet is assumed to initially be on an approximately circular orbit, but locked in resonance, at t = 0. Using the standard Keplerian relations that p Ei p = −GMs Mi /(2ai ), Ωo = GMs /a3o , and Ji = Mi GMs ai (1 − e2i ), with semi-major axis ai , it is straightforward to show from equation (105) that s p ai (0) 2 1 − ei (t) = λ + (1 − λ) , (106) ai (t)

Numerical simulations provide an important tool for analyzing planet migration. They can provide important insights in cases where nonlinear and time-dependent effects are difficult to analyze by analytic methods. Some powerful grid-based hydrodynamics codes (such as the Zeus code Stone & Norman, 1992) have been adapted to the study of disk-planet interactions. In addition, particle codes based on the Smoothed Particle Hydrodynamics (SPH) (Monaghan, 1992) have sometimes been employed. A systematic comparison between many of the codes has been carried out by de Val-Borro et al, 2006. We discuss a few basic points. For planets that open a gap, grid-based codes offer an advantage over particlebased codes. The reason is that the resolution of grid based codes is determined by the grid spacing, while the resolution of particle-based codes is determined by the particle density. If a planet opens an imperfect gap, the particle density and resolution near the planet is low. Poor resolution near the planet can give rise to artificial torques. Higher resolution occurs where the particle density is higher, but this occurs in regions that interact less strongly with the planet. There have been variable resolution techniques developed for grid-based codes in which the highest resolution is provided in regions near the planet where it is needed, such as nested grid methods (D’Angelo et al, 2002). Such techniques need to provide a means of joining the regions of high and low resolution without introducing artifacts (such as wave reflections) or lowering the overall accuracy of the scheme. Grid-based codes that simulate disk-planet interactions typically employ numerical devices to improve convergence. For example, the gravitational potential of the planet is often replaced by one that does not diverge near the planet. The potential is limited by introducing a smoothing length, a distance within which the potential does not increase near the planet. Another limitation is that the simulated domain of the disk is typically limited to a region much smaller than the full extent of the disk. Techniques have been developed to ensure that reflections from the boundaries do not occur, e.g., by introducing enhanced wave dissipation near the boundary or approximate outgoing wave boundary conditions. The time-steps of codes are limited by the Courant condition. Short time-steps often result from the region near the inner boundary of the computational domain (smaller radii) where the disk rotation is fastest. As a result, it is difficult to extend the disk very close to the central star. The FARGO scheme (Masset, 2000) is very useful in overcoming this limitation. However, the method is difficult to apply to a variable grid spacing code. Convergence is a major issue with these simulations. Ideally, one should demonstrate that the results of the simulations are sufficiently insensitive to the locations of the boundaries, the size of the smoothing length, the size of the time steps, and the grid resolution. Even the direction of migration can be affected by the size of potential smoothing length in certain cases. In a 2D simulation, a finite smooth-

where λ = Ωi /Ωo > 1 is the ratio of orbital frequencies that is a constant for resonantly locked planets. The result implies that eccentricity formally goes to unity where ai (t) = ai (0)



λ−1 λ

2

.

(107)

For planets locked in a 2:1 resonance, we have that λ = 2 and the eccentricity goes to unity for ai (t) = ai (0)/4. The analysis has assumed the orbits remain coplanar. Numerical simulations have shown that before reaching a radial orbit, the inner planet becomes strongly inclined relative to the outer one, if it starts with a small nonzero initial inclination. If the inner planet manages to miss striking the star as it passes ei = 1, further inward migration can cause the inner planet’s orbit to flip over and change the sense of its orbital motion to be counterrotating. Qualitatively similar effects occur for planetary systems with nonextreme mass ratios. Eccentricity is generated in both planets (see Fig. 17). Excitation of inclination requires that Mo > Mi /2. Planetary system GJ876 is a well-studied case in which the planets are in a 2:1 resonance. If the system’s measured eccentricities are due to resonant migration, then according to theory (Lee & Peale, 2002), the system migrated inward by less than 10%. Such a small amount of migration is hard to understand. It is possible that eccentricity damping through disk-planet interactions could have limited the eccentricities to the observed levels as the planets underwent further migration. But the required damping rate is quite high. This level of damping is more than an order of magnitude higher than is found in hydrodynamic simulations of this system (Kley et al, 2005). Once the disk dissipates, migration ceases and eccentricity growth by this mechanism is terminated. The disk might have dissipated after the planets achieved convergent migration and underwent a small amount of further migration, but the timing seems somewhat unlikely. 3.3 Validity of Numerical Simulations 17

ing length ∼ H provides a means of simulating the reduced effects of planet gravity on a disk of finite thickness. But in the 2D case, the limit of zero smoothing length is unphysical. In practice, testing for convergence is computationally expensive, but can be done for a subset of the models of interest, perhaps over a limited time range. Demonstrating convergence is important for providing reliable results.

4. FUTURE PROSPECTS The theory of planet migration is interwined with the theory of planet formation (Chapter 19). The timescale for a planet to grow within a disk is a key element in understanding whether planet migration is a major obstacle to planet formation. We have pointed out that the formation timescales for gas giant planets in the core accretion model do present a problem for the simplest planet migration theories. However, alternative migration models, some of which are described in Section 3, may be appropriate. Future prospects for resolving this issue depend on advances in the theory of planet formation. As we have seen in Section 2, planet migration in disks is a consequence of the action of Lindblad and corotational resonances. The Lindblad resonances are better understood. Their linear and nonlinear properties have been analyzed in more detail. The corotational resonances are somewhat more delicate and less well understood. Progress on the theory planet migration will likely involve further investigations of the role of corotational resonances. Improvements to computer capabilities should allow longer simulations to be carried out with higher resolution. Progress will be made by also including more physical effects. For example, most multi-dimensional simulations have made only the simplest assumptions about the thermal properties of the disk. Some calculations have suggested that higher mass eccentric orbit planets could migrate outwards due to their interactions with a gaseous disk. But we do not know whether such planets have acquired their eccentricities at this early stage. Observations of young planets would be quite valuable in understanding this issue. A major uncertainty in the theory of planet migration is the physical state of the disk. We have seen in Section 3 that low mass planet migration behaves very differently depending on the level of disk turbulence and the structural properties of the disk. For example, in weakly turbulent disks, feedback effects can limit Type I migration, while in highly turbulent disks torque fluctuations may play a role in modifying Type I migration. The presence of rapid radial density variations in the disk can substantially alter migration, since it depends on the competition between torques involving material just inside and outside the orbit of the planet. The better determination of disk properties will likely rely on some combination of improved theory and observations. It is unlikely that theory alone will be able to make much

Fig. 1.— Path of a particle that passes by a planet of mass Mp = 10−6 Ms . The coordinates are in units of the orbital radius of the planet rp . The planet lies at the origin, while the star lies at (−1, 0). The dashed line follows the path that is undisturbed by the planet with x = 0.025rp , while the solid line follows the path resulting from the interaction with the planet. In the frame of the planet, the particle moves in the negative y direction.

1.030 1.028 1.026 1.024 1.022

-3

-2

-1

1

2

3

Fig. 2.— The solid line plots the particle distance from the star r/rp as a function of time in units of planet orbit periods for the particle that follows the perturbed path in Fig. 1. The particle passes the planet at time t = 0. Immediately after passage by the planet, the particle is deflected toward smaller radii, toward the planet, and acquires an eccentricity, as indicated by the radial oscillations. The dashed line plots the mean radius of these oscillations. Since the mean radius of the oscillations is larger than the initial orbital radius (dashed versus solid line at t < −1), the particle gained energy and angular momentum, as a consequence of its interaction with the planet.

18

Fig. 3.— Torque on the particle normalized by M Ω2p rp2 as a function of time in units of planet orbit periods along the unperturbed (dashed lines) and perturbed paths (solid lines). The planet mass Mp = 10−6 Ms . Top panel is for a particle having x = 0.02rp , the case in Fig. 1. The bottom panel is for particle for a particle with x = 0.03rp .

Fig. 5.— Schematic of acoustic (pressure) wave propagation in a gas disk. Leading waves (dashed wavefronts) propagate towards corotation (CR), towards the orbit of the planet. Trailing waves (solid wavefronts) propagate away from corotation (CR), away from the orbit of the planet. The region between the inner Lindblad resonance (ILR) and outer Lindblad resonance (OLR) is evanescent (nonpropagating) and propagating elsewhere.

Fig. 6.— Radius of a Lindblad resonance in units of rp as a function of azimuthal wave number m. The two sets of lighter dots follow the standard WKB dispersion relation for density waves with pressure, equation (34), in a Keplerian disk with radii given by equation (38). The two sets of darker dots follow the extended WKB dispersion relation (60) that accounts for azimuthal effects with c = 0.1Ωr. The dashed line is the location of the planet. Points below (above) the dashed line are for inner (outer) Lindblad resonances. In the standard WKB approximation, the resonances get closer to the planet with increasing m. In the extended WKB approximation, the resonances maintain a fixed separation from the planet with increasing m.

Fig. 4.— Numerical test of equation (7) based on orbit integrations. Top panel: Log-log plot of 104 ∆J/(M rp2 Ωp ) as a function of x/rp . The lower set of points is for a planet of mass Mp = 10−6 Ms . The upper set is for Mp = 10−3 Ms . The solid lines are for ∆J ∝ x−5 that pass through the respective rightmost points. Bottom panel: Log-log plot of 105 ∆J/(M rp2 Ωp ) as a function of 105 Mp /Ms . The points are the results of numerical simulations for a fixed value of x = 0.2rp . The solid line is for ∆J ∝ (Mp /Ms )2 that passes through the left-most point.

19

Fig. 7.— Cumulative outer Lindblad resonance torque in a gas Fig. 9.— Coorbital streamlines near a Saturn mass planet Mp =

disk starting at some radius well inside the resonant radius rm to some radius r. The abscissa is the scaled radial distance from the resonance y = (r − rm )/(wrm ), where w ≪ 1 is the dimensionless resonance width. The cumulative torque is normalized by its value at large y. A wave launched at the resonance propagates for y > 0 and is evanescent for y < 0. The torque accumulates within a distance of ∼ wrm from rm .

3 × 10−4 Ms located at the origin in the corotating frame of the planet. Angle θ increases to the left. Solid streamlines are the horseshoe orbits. Locations i∓ and o∓ label the positions near the encounter with the planet.

Fig. 10.— Total angular momentum for a set of 60 particles, each of mass M , on horseshoe orbits in units of M rp2 Ωp as a function of time in units of the planet orbit period in a star-planet system with Mp = 3 × 10−6 Ms . The particles start at t = 0 distributed between r = rp + RH /60 to r = rp + RH (with RH ≃ 7 × 10−3 rp ) along the star-planet axis 180 degrees from the planet. The changes in angular momentum cause a torque to be exerted on the planet. The torque oscillates and declines in time as angular momentum becomes more constant in time because the particles undergo phase mixing on the libration timescale ∼ 150 planet periods.

Fig. 8.— Ψ2m /(Mp2 Ω4p rp4 ), defined in equation (46), versus m. The vertical axis is plotted with logarithmic spacing. The lighter dots are evaluated for a Keplerian disk at radii of outer Lindblad resonances given by the standard WKB approximation, while the darker dots are evaluated at radii of outer Lindblad resonances given by the extended WKB approximation with c = 0.1Ωp rp . These radii are plotted in Fig. 6. For m & Ωp rp /c = 10, quantity Ψ2m declines in the extended WKB approximation.

20

Fig. 12.— Migration of a planet undergoing growth via gas accretion. The disk parameters are similar to those in Fig 11. Top: The vertical axis is the orbital radius in units of the initial orbital radius 5.2AU . The horizontal axis is time in units of the initial orbital period 12y. The solid curve is the result of 3D hydrodynamical simulations. The lower and upper short dashed curves are based on equations (87) and (88) respectively, applied to a planet of variable mass. The long dashed curve corresonds migration on the disk viscous timescale. Bottom: Average disk density near the planet relative to the initial value as a function of time. The density is averaged over a band of radial width 2H centered on the orbit of the planet and is normalized by its initial value. The first solid circle marks the time when Mp = 16.7M⊕ , subsequent circles occur at integer multiples of 33M⊕ . The planet initially follows the predictions of Type I migration theory in the top panel while there is no substantial gap in the disk (no drop in the curve on the bottom figure). After gap opening, the planet follows Type II migration. Obtained from D’Angelo & Lubow, 2008.

Fig. 11.— Migration timescales versus planet mass for a planet embedded in a 3D disk of mass ∼ 0.02Ms with density Σ ∝ r −1/2 and H/r = 0.05. The planets are on fixed circular orbits and have fixed masses. The dots with error bars denote results of 3D numerical simulations with the same disk parameters and α = 0.004 (Bate et al, 2003). The dashed line plots equation (88) based on linear theory (Tanaka et al, 2002). Above about 0.1MJ , the planet opens a gap in the disk, Type I theory becomes invalid, and Type II migration occurs.

21

Fig. 13.— Scaled torque per unit disk mass on the planet as a function of (r − rp )/rp based on 3D simulations. The vertical scale is in units of GMs (Mp /Ms )2 /rp . The solid and longdashed curves are for 1 M⊕ and 10 M⊕ mass planets, respectively. The disk parameters are H/r = 0.05 and α = 4 × 10−3 for both the cases. According to linear theory, these two curves should overlap. Torque distributions are averaged over one orbital period. Figure based on D’Angelo & Lubow, 2008.

Fig. 15.— Influence of disk viscosity parameter α on the migration of a planet with mass 10M⊕ in a disk with H/r = 0.035 and disk of mass ∼ 0.1M⊙ based on 2D simulations. The vertical axis is the orbital radius in units of the initial orbital radius. The horizontal axis is the time in units of the initial planet orbital period. For α = 8 × 10−4 the migration follows the Type I rate. At the lower values, the migration halts due to a feedback effect. Figure based on Li et al, 2008.

Fig. 16.— Coorbital streamlines near an inwardly migrating Sat-

Fig. 14.— Azimuthally averaged disk surface density normal-

urn mass planet Mp = 3 × 10−4 Ms located at the origin in the comoving frame of the planet. The inward migration rate is 0.002Ωp rp . Angle θ increases to the left. The heavy streamlines are open and pass by the planet. The light streamlines are closed and contain trapped material on the leading side of the planet (streamlines based on Ogilvie & Lubow, 2006). The trapped material is retained from regions further from star acquired at earlier time. The open streamlines carry ambient disk material. In contrast to nonmigrating case in Fig. 9, the asymmetry and the density differences between the trapped (retained) and open (ambient) gas gives rise to a coorbital torque.

ized by the unperturbed value at r = 1 as a function of r in units of rp = 5.2AU . The disk is simulated in 3D with parameters H/r = 0.05 and α = 0.004. The density profiles are for planets with masses of 1 (long-dashed), 0.3 (dot-dashed), 0.1 (dotted), 0.03 (short-dashed), and 0.01 (thin solid) MJ . Only planets with masses Mp & 0.1MJ produce significant perturbations. The thick solid line is based on a 2D simulation of a 1MJ planet by Lubow et al (1999). Obtained from Bate et al, 2003.

22

progress along these lines because the density structure depends on the disk turbulence which is difficult to accurately model from first principles. The resolution of disk properties in the inner parts of protostellar disks is important for such purposes. New telescopes such as ALMA and JWST may be quite valuable in making such determinations. Acknowledgments. This work was partially supported by NASA grant NNX07AI72G to SL. We thank Gennaro D’Angelo and Jim Pringle for carefully reading a draft and suggesting improvements.

REFERENCES Alibert, Y., Mousis, O., Mordasini, C. and Benz, W. (2005) New Jupiter and Saturn Formation Models Meet Observations. Astrophys. J., 626, L57-L60. Artymowicz, P. (1993a) On the Wave Excitation and a Generalized Torque Formula for Lindblad Resonances Excited By External Potential. Astrophys. J., 419, 155-165. Artymowicz, P. (1993b) Disk-Satellite Interaction via Density Waves and the Eccentricity Evolution of Bodies Embedded in Disks. Astrophys. J., 419, 166-180. Artymowicz, P. and Lubow, S. H. (1996) Mass flow through gaps in circumbinary disks. Astrophys. J., 467, L77. Artymowicz, P. (2004) Migration Type III. KITP Conference on Planet Formation, http://online.kitp.ucsb.edu/online/planetf c04 Balbus, S. A., and Hawley, J. F. (1991) A powerful local shear instability in weakly magnetized disks. I. Linear analysis. Astrophys. J., 376, 214-222. Bate, M. R., Lubow, S. H., Ogilvie, G. I., and Miller, K. A. (2003) Three-dimensional calculations of high- and low-mass planets embedded in protoplanetary discs. Mon. Not. Roy. Astr. Soc., 341, 213-229. Bodenheimer, P., Hubickyj, O., and Lissauer, J. J. (2000) Models of the in Situ Formation of Detected Extrasolar Giant Planets. Icarus, 143, 2. Butler, R. P., et al. (2006) Catalog of Nearby Exoplanet. Astrophys. J., 646, 505-522. Cuzzi, J. N., and Weidenschilling, S. J. (2006), Particle-Gas Dynamics and Primary Accretion, in Meteorites and the Early Solar System II, (Lauretta & McSween Jr. eds.), Univ. Arizona Press, Tucson. D’Angelo, G., Lubow, S. H., and Bate, M. R. (1996) Evolution of Giant Planets in Eccentric Disks. Astrophys. J., 652, 16981714. D’Angelo, G., Henning, T., and Kley, W. (2002) Nested-grid calculations of disk-planet interaction. Astron, Astrophys., 385, 647-670. D’Angelo, G., and Lubow, S. H. (2008) Evolution of Migrating Planets Undergoing Gas Accretion. Astrophys. J., 685, 560583. Goldreich, P., and Tremaine, S. (1979) The excitation of density waves at the Lindblad and corotation resonances by an external potential. Astrophys. J., 233, 857-871. Goldreich, P., and Tremaine, S. (1980) Disk-satellite interactions. Astrophys. J., 233, 425-441. Goldreich, P., and Tremaine, S. (1982) The Dynamics of Planetary Rings. Ann. Rev. Astron. Astrophys., 20, 249-283.

Fig. 17.— Evolution of orbital eccentricity (top) and inclination (bottom) of a two planet system in which an inward migration of 10−5 AU y −1 is imposed on the outer one. The horizontal axis is in units of years. Both planets have Mp = 1MJ . The orbital eccentricities of the planets grow as a result of resonant migration in the 2:1 resonance. The eccentricity of the inner planet grows faster than the outer one. A small initial mutual inclination grows to large values as the eccentricities increase. From Tommes & Lissauer, 2003.

Fig. 18.— Orbital radius as a function of time in orbits for 3M⊕ planets embedded in a gaseous disk based on simulations. The dotted lines plot the migration of planets in a laminar disk, a disk without turbulent fluctuations. The solid lines plot the migration of planets in a disk with turbulent fluctuations due to the MHD turbulence. The random motions are due to the fluctuating torques. Obtained from Nelson, 2005.

23

tion, and the Kozai mechanism. Astrophys. J., 678, 498-508. Nelson, R. P. (2005) On the orbital evolution of low mass protoplanets in turbulent, magnetised disks. Astron. Astrophys. 443, 1067-1085. Ogilvie, G. I., and Lubow, S. H. (2003) Saturation of the Corotation Resonance in a Gaseous Disk. Astrophys. J., 587, 398406. Ogilvie, G. I., and Lubow, S. H. (2006) The effect of planetary migration on the corotation resonance. Mon. Not. Roy. Astr. Soc., 370, 784-798. Paardekooper, S.-J., and Mellema, G. (2006) Halting type I planet migration in non-isothermal disks. Astron. Astrophys. 459, L17-L20. Papaloizou, J. C. B. (2002) Global m = 1 modes and migration of protoplanetary cores in eccentric protoplanetary discs. Astron. Astrophys. 388, 615-631. Pringle, J. E. (1991) The properties of external accretion discs. Mon. Not. Roy. Astr. Soc., 248, 754-759. Rafikov, R. R. (2002) Planet Migration and Gap Formation by Tidally Induced Shocks. Astrophys. J., 572, 566-579. Stone, J. M., and Norman, M. L. (2002) ZEUS-2D: A radiation magnetohydrodynamics code for astrophysical flows in two space dimensions. I - The hydrodynamic algorithms and tests. Astrophys. J. Supp., 80, 753-790. de Val-Borro, M., et al. (2006) A comparative study of disc-planet interaction. Mon. Not. Roy. Astr. Soc., 370, 529-558. Weidenschilling, S. J. (1977) Aerodynamics of solid bodies in the solar nebula. Mon. Not. Roy. Astr. Soc., 180, 57-70. Wu, Y., and Murray, N. (2003) Planet Migration and Binary Companions: The Case of HD 80606b. Astrophys. J., 589, 605-614. Salmeron, R., & Wardle, M. (2008) Magnetorotational instability in protoplanetary discs: the effect of dust grains. Mon. Not. Roy. Astr. Soc., 388, 1223-1238. Shu, F. (1991) Physics of Astrophysics, Vol. II: Gas Dynamics, University Science Books, NY

Goldreich, P., and Sari, R. (2003) Eccentricity Evolution For Planets in Gaseous Disks. Astrophys. J., 585, 10241037. Hahn, J. M., & Malhotra, R. (1999) Orbital Evolution of Planets Embedded in a Planetesimal Disk. Astron. J., 117, 3041-3053. Ida, S., Bryden, G., Lin, D. N. C. and Tanaka, H. (2000) Orbital migration of Neptune and orbital distribution of transNeptunian objects. Astrophys. J., 534, 428-445. Ida, S., and Lin, D. N. C. (2004) Toward a deterministic model of planetary formation. I. A desert in the mass and semimajor axis distributions of extrasolar planets. Astrophys. J., 604, 388-413. Ida, S., Guillot, T., and Morbidelli, A. (2008) Accretion and destruction of planetesimals in turbulent disks. Astrophys. J., 686, 1292-1301. Ida, S., and Lin, D. N. C. (2008) Toward a deterministic model of planetary formation. IV. Effects of type I migration. Astrophys. J., 673, 487-501. Johnson, E. T., Goodman, J., and Menou, K. (2006) Diffusive Migration of Low-Mass Protoplanets in Turbulent Disks. Astrophys. J., 647, 1413-1425. Kley, W., Lee, M. H., Murray, N., and Peale, S. J. (2005) Modeling the resonant planetary system GJ 876. Astron. Astrophys. 437, 727-742. Lee, M. H., and Peale, S. J. (2002) Dynamics and Origin of the 2:1 Orbital Resonances of the GJ 876 Planets. Astrophys. J., 567, 596-609. Li, H., Lubow, S. H., Li, S., and Lin, D. N. C. (2008) Type I Planet Migration in Nearly Laminar Disks. Astrophys. J.,, arXiv:0812.0586. Lin, D. N. C., and Papaloizou, J. (1979) Tidal torques on accretion discs in binary systems with extreme mass ratios. Mon. Not. Roy. Astr. Soc., 186, 799-812. Lin, D. N. C., and Papaloizou, J. C. B. (1986) On the tidal interaction between protoplanets and the protoplanetary disk. III Orbital migration of protoplanets. Astrophys. J., 309, 846-857. Lubow, S. H. (1991) A model for tidally driven eccentric instabilities in fluid disks. Astrophys. J., 381, 259-267. Lubow, S. H. and Pringle, J. E. (1993) Wave Propagation in Accretion Disks. Astrophys. J., 409, 360-371. Lubow, S. H. and Ogilvie, G. I. (1998) Three Dimensional Waves Generated at Lindblad Resonances in Thermally Stratified Disks. Astrophys. J., 504, 983-995. Martin, R. G., Lubow, S. H., Pringle, J. E., & Wyatt, M. C. (2007) Planetary migration to large radii. Mon. Not. Roy. Astr. Soc., 378, 1589-1600. Masset, F. S. (2000) FARGO: A fast eulerian transport algorithm for differentially rotating disks. Astron. Astrophys. Supp., 141, 165-173. Masset, F. S., and Papaloizou, J. C.B. (2003) Runaway migration and the formation of hot Jupiters . Astrophys. J., 588, 494508. Masset, F. S., Morbidelli, A., Crida, A., and Ferreira, J. (2006) Disk surface density transitions as protoplanet traps. Astrophys. J., 642, 478-487. Mayor, M., and Queloz, D. (1995) A Jupiter-mass companion to a solar-type star. Nature, 378, 355-359. Menou, K., and Goodman, J. (2004) Low-Mass Protoplanet Migration in T Tauri Disks. Astrophys. J., 606, 520-531. Meyer-Vernet, N., and Sicardy, B. (1987), On the physics of resonant disk-satellite interaction. Icarus, 69, 157-175 Monaghan, J. J. (1992) Smoothed particle hydrodynamics. Ann. Rev. Astron. Astroph., 30, 543-574. Nagasawa, M., Ida, S., and Bessho, T. (2008) Formation of hot planets by a combination of planet scattering, tidal circulariza-

Tanaka, H., Takeuchi, T., & Ward, W. R. (2002) ThreeDimensional Interaction between a Planet and an Isothermal Gaseous Disk. I. Corotation and Lindblad Torques and Planet Migration. Astrophys. J., 565, 1257-1274. Thommes, E. W., and Lissauer, J. J. (2003) Resonant Inclination Excitation of Migrating Giant Planets. Astrophys. J., 597, 566-580. Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. (2005) Origin of the orbital architecture of the giant planets of the Solar System. Nature, 435, 459-461. Ward, W. R. (1991) Horeshoe Orbit Drag. Lunar and Planetary Institute Conference Abstracts, 22, 1463 Ward, W. R. (1992) Coorbital Corotation Torque. Lunar and Planetary Institute Conference Abstracts, 23, 1491 Ward, W. R. (1997) Protoplanet Migration by Nebula Tides. Icarus, 126, 261-281. Yu, Q., and Tremaine, S. (2001) Resonant Capture by Inward-migrating Planets Astron. J., 121, 1736-1740. Zhang, H., Yuan, C., Lin, D. N. C., and Yen, D. C. C. (2008), On the Orbital Evolution of a Jovian Planet Embedded in a Self-Gravitating Disk. Astrophys. J., 676, 639-650.

24