THREE-DIMENSIONAL DISK–PLANET TORQUES IN A LOCALLY ISOTHERMAL DISK

Important interactions occur between a young planet and its surrounding gaseous disk. Such interactions are generally strong and can lead to planet migration on relatively short timescales (Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Ward 1997, see also review by Lubow & Ida 2010). The shortness of the migration timescales may be a problem, since they can be shorter than planet formation timescales (≳10⁶ yr) predicted by the core nucleated accretion models of giant planet formation (e.g., Lissauer et al. 2009 and references therein).

There are distinct regimes of planet migration. The Type I regime applies to low-mass planets whose tidal torques are too weak to open a gap in the disk. In this case, the planet is fully embedded in the disk that is nearly unperturbed by the presence of the planet. Several calculations have been carried out to determine a Type I torque formula for disk–planet interactions, some of which are described below. The tidal torque on the planet involves opposing contributions from the material outside and inside the orbit of the planet. As a result, the net torque on the planet depends on both the values of the gas properties (such as density and temperature) and their radial gradients near the orbit of the planet. Several studies have investigated the possibility of stopping migration by trapping the planet in regions where the radial gradients of disk properties may be strong enough to counteract the general tendency to inward migration. Such disk regions may occur where the density and temperature vary rapidly in radius due to changes in disk opacity or turbulent viscosity (e.g., Menou & Goodman 2004, hereafter MG04; Matsumura et al. 2007). Disk temperature gradients are modified near a planet due to enhanced stellar exposure and shadowing effects caused by the variations in the vertical thickness of the disk. The modified temperature structure due to this effect may also act to somewhat slow migration (Jang-Condell & Sasselov 2005).

Most torque calculations to date have assumed that the disk is locally isothermal and the disk viscosity is not very low (α ≳ 0.001). Under the locally isothermal assumption, the gas temperature is a fixed imposed function of distance from the central star that is independent of height above the disk midplane. This approximation is formally justified in disk regions where the thermal timescale is shorter than the local orbital period. This approximation represents a limitation of our analysis, but it allows us to directly connect the results with those of previous studies that made the same approximation.

Some studies have suggested that important modifications to the migration rate can occur if the disk coorbital region responds nearly adiabatically (e.g., Baruteau & Masset 2008; Paardekooper & Papaloizou 2008) or if the disk turbulent viscosity is very low (Ward 1997; Rafikov 2002; Li et al. 2009). In the nearly adiabatic case, Type I outward migration appears possible if (1) the turbulent viscosity is sufficiently high to allow the corotation torques to be effective, (2) the disk thermal timescale near the planet is longer than the time for the gas to make U-turns in its horseshoe orbit near the planet (of order the planet's orbital period), (3) the disk thermal timescale is shorter than the gas libration timescale in the coorbital region, and (4) the disk radial entropy gradient is negative. The radial entropy gradient is zero if the surface density is ∝r^−3/2, as in the standard Minimum-Mass Solar Nebula (Hayashi 1981), and the temperature is ∝r⁻¹, as it would be the case in a disk with constant opening angle.

A weak turbulent viscosity regime (α ≲ 10⁻⁴) may occur in "dead zones" of disks where gas ionization levels may be too low to sustain magnetically driven turbulence below the disk's outer layers (e.g., Gammie 1996). Under such conditions, two-dimensional simulations suggest that the migration of a ≈10 M_E planet can be halted by the strongly nonlinear effects of waves launched at Lindblad resonances that result in shocks (Li et al. 2009; Yu et al. 2010). At sufficiently low levels of turbulent viscosity (α ≲ 10⁻⁵), vortices can form and collide with the planet, producing a jumpy and chaotic form of migration.

We restrict attention in this paper to locally isothermal disk regions with moderate levels of disk turbulent viscosity (5 × 10⁻⁴ ≲ α ≲ 0.05). For such cases, the well-known disk torque formula of Tanaka et al. (2002, hereafter TTW02) describes the situation of a disk that is spatially isothermal (the temperature is radially and vertically constant) by means of a linear analysis of the three-dimensional dynamical equations. Various three-dimensional nonlinear simulations have obtained good agreement with these migration rates (e.g., D'Angelo et al. 2003; Bate et al. 2003). However, there has not been a systematic analysis, based on three-dimensional simulations, to study the effects of temperature and density gradients on migration, which we provide here.

Others, such as MG04, generalized the three-dimensional analytic disk torque calculation of TTW02 to locally isothermal cases where radial disk temperature variations occur. The analysis of MG04 was done using two-dimensional resonant torque expressions, together with a softened potential to mimic the effects of the third dimension (perpendicular to the disk midplane). This calculation included only the effects of Lindblad torques and neglected the corotation torques. However, upon closer inspection of the three-dimensional Lindblad torque contribution of TTW02, we find that it disagrees with the results of MG04, applied to a fully isothermal disk, in important ways. This discrepancy suggests that the torque dependence on disk gradients and the trapping mechanisms in general require further examination. One goal of this paper is to determine the torque formula for a planet embedded in a three-dimensional, locally isothermal disk by means of numerical simulations that include effects of density and temperature gradients, as well as the co-orbital torque (see Section 2). This analysis is described in Section 3.

The distribution of torques per unit disk mass, as a function of radius, provides a description of where the disk contributions to the planet torque occur (e.g., D'Angelo & Lubow 2008, hereafter DL08). Such torque distributions are useful in the application to one-dimensional disk evolution models that account for planets, including cases where the planets open gaps in the disk, the so-called Type II regime of planetary migration (see, e.g., Ward 1997 and references therein). Current multi-dimensional hydrodynamical codes cannot practically follow the disk evolution over long times, such as disk viscous times, while one-dimensional codes can. The torque density distribution is typically estimated by means of the impulse approximation (Lin & Papaloizou 1986). But that distribution contains a long tail that is unphysical, particularly near the star (e.g., Veras & Armitage 2004). Another goal of the paper is to provide a more accurate formula for the torque density distribution that can be applied to such and similar problems, again based on three-dimensional simulations. These distributions are described in Section 4. Some applications of these formulae to migrating super-Earths are discussed in Section 5. Extensions of the torque density formula to situations of mildly varying disk gradients and gap-opening planets are presented in Sections 6 and 7, respectively. Section 8 contains the summary of the results.

2.1. Disk Model

The simulations are carried out along the lines of DL08. We use a spherical polar coordinate system {O; R, θ, ϕ}, whose origin O is located at the star–planet center of mass and which rotates at the same angular velocity as the planet orbits about the star. The planet's orbit is circular and has zero inclination relative to the disk midplane θ = π/2. Since the disk is assumed to be symmetric relative to its midplane, only the disk's northern hemisphere, θ ⩽ π/2, is explicitly simulated.

We assume that the disk's gas is locally isothermal and that the pressure, p, can be written as

$$\begin{equation} p(R, \theta, \phi) = \rho (R, \theta, \phi)\,c^{2}_{s}(r), \end{equation} \tag{ 1 }$$

where ρ(R, θ, ϕ) is the mass density. The sound speed of the gas, c_s(r), is assumed to be a function of cylindrical radius r = Rsin θ. From the equation of state for a perfect gas and Equation (1), it follows that the temperature T ∝ p/ρ = c²_s. Therefore, to determine the effects of temperature gradients on disk torques, we vary T in radius by varying the gas sound speed, which is taken to be a power law in r. The hydrostatic disk scale height is defined as

$$\begin{equation} H(r)=\frac{c_{s}}{\Omega _{\mathrm{K}}}, \end{equation} \tag{ 2 }$$

where Ω_K is the Keplerian angular velocity. Notice that, since the vertical motion and structure of the disk is directly computed from the three-dimensional Navier–Stokes equations, the actual scale height of the disk may deviate from the quantity H(r). For brevity, however, we will refer to H(r) simply as disk scale height.

The disk extends in azimuth over 2π radians around the star. In models with planet-to-star mass ratio M_p/M_s < 3 × 10⁻⁴, the radial disk region extends from 0.4a to 2.5a or from 0.5a to 1.55a, where a is the orbital radius of the planet. In models with planet-to-star mass ratio M_p/M_s ⩾ 3 × 10⁻⁴, the disk extends from 0.4a to 4.0a in radius. In the θ-direction, the disk domain contains from a few to many disk scale heights above the midplane, depending on both r and H(r).

The initial mass density, ρ, is independent of the azimuth, ϕ, and has a Gaussian profile in the meridional direction, θ. In the radial direction, ρ ∝ Σ/H and the initial surface density, Σ, is usually taken to be proportional to a power of r. In Section 6, we consider a more general case. The initial flow field has radial and meridional components u_R = −3 ν/(2R), where ν is the kinematic viscosity of the gas, and u_θ = 0. The initial azimuthal velocity is u_ϕ = r [Ω − Ω_K(a)]. The initial (unperturbed) rotation rate Ω accounts for the effects of pressure support produced by both density and temperature gradients. For a pressure supported disk

$$\begin{equation} \Omega ^{2}=\Omega ^{2}_{\mathrm{K}} +\frac{1}{r \rho }\left(\frac{\partial p}{\partial r}\right). \end{equation} \tag{ 3 }$$

Although Ω has a dependence on distance from the midplane (see, e.g., Equation (4) of TTW02), we adopt as an initial condition its value at the midplane that is given by

$$\begin{equation} \Omega ^{2}=\Omega ^{2}_{\mathrm{K}}\left[1- \left(\beta +\frac{\zeta }{2}+\frac{3}{2}\right) \left(\frac{H}{r}\right)^{2}\right]. \end{equation} \tag{ 4 }$$

In the above equation, we set β = −dln Σ/dln r and ζ = −dln T/dln r. Notice that, in an infinitesimally thin (i.e., two-dimensional) disk, Equation (3) involves the surface density and a vertically integrated pressure and thus Equation (4) becomes Ω² = Ω²_K[1 − (β + ζ)(H/r)²].

We typically impose a kinematic viscosity proportional to the inverse of the initial surface density distribution, so that the mass accretion rate through the disk (∝ν Σ) is nearly constant in radius and the disk can achieve a steady state. In such cases, the initial imposed value of the parameter β is therefore preserved in time. In the calculations discussed in Section 3, ν(a) = 10⁻⁵ a² Ω_K(a), which corresponds to a turbulent viscosity parameter α = 0.004 at the orbital radius of the planet when H(a)/a = 0.05. However, we also report on models with lower (as low as α ∼ 5 × 10⁻⁴) and higher kinematic viscosity, up to ν(a) ∼ 10⁻⁴ a² Ω_K(a) or α ∼ 0.05 (see also Section 4).

We consider several values of the planet-to-star mass ratio, M_p/M_s. For the Type I cases discussed in Section 3, we mainly concentrate on M_p/M_s = 3 × 10⁻⁶ (or 1 M_E planet if M_s = 1 M_☉). For Type II cases discussed in Section 7, we consider planets up to two Jupiter masses (M_p = 2 × 10⁻³ M_s). The gravitational potential of the planet, Φ_p, is smoothed over a fraction ε of the Hill radius of the planet R_H = a[M_p/(3 M_s)]^1/3 and is given by

$$\begin{equation} \Phi _{p}=-\frac{G\mbox{$M_{p}$}}{\sqrt{S^2+\varepsilon ^2 R_{\mathrm{H}}^2}}, \end{equation} \tag{ 5 }$$

where S is the distance from the planet. We generally apply a smoothing (or softening) length with parameter ε = 0.1. In the low-mass regime, another relevant length scale is the Bondi radius R_B = GM_p/c²_s (e.g., Masset et al. 2006; DL08), which can be shorter than R_H. In these situations, the smoothing length we apply is of order R_B or smaller.

The softening length plays a much weaker role in three dimensions than it does in two dimensions, since the vertical thickness of the disk effectively smooths the potential of the planet. We find that our results are insensitive to changes in values of the softening length. For cases with an M_p = 5 × 10⁻⁶ M_s planet, changes in the softening length by a factor of 10 produced a change in the net torque much less than 1%. For cases with an M_p = 10⁻³ M_s planet, changing ε by a factor of 2 resulted in a net torque variation of about 1%. Further details are given in the Appendix.

2.2. Numerical Method

3.1. Comparison of Three-dimensional and Softened Two-dimensional Torques in a Fully Isothermal Disk

Three-dimensional linear analytic calculations of the Type I migration rates were carried out by TTW02. They provide torque formulae for a disk where the gas sound speed is strictly constant, i.e., for an isothermal disk. The torque exerted on a planet, including contributions from Lindblad resonances and the corotation resonance, is given by

$$\begin{eqnarray} \mathcal {T}_{\rm{{TTW}}} &=&-\left(1.36 + 0.54\,\beta \right)\nonumber \\ & &\times \Sigma (a) \, \Omega ^{2}(a) \, a^{4} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2} \left(\frac{a}{H}\right)^{2}, \end{eqnarray} \tag{ 6 }$$

where the (azimuthally averaged) gas surface density Σ, the disk rotation rate Ω, and the disk vertical scale height H are evaluated at the orbital radius of the planet, a. We recall that β = −dln Σ/dln r. Equation (6) does not include the effects of temperature gradients, since the disk is assumed to be fully isothermal.

The corotation resonance contains the complication that it can be saturated (have zero strength). In the absence of irreversible processes, the flow in the coorbital region follows closed horseshoe and tadpole orbits in the corotating frame of the planet. The steady torque associated with such closed orbits is zero (saturated). Turbulent viscosity introduces an irreversible behavior that results in a nonzero torque. The effects of turbulent viscosity are stronger for lower mass planets (R_H < H), since the radial width of the horseshoe orbit region, $w\approx 2 a\sqrt{(\mbox{$M_{p}$}/\mbox{$M_{s}$})(a/H)}$ (see Masset et al. 2006), is narrower and the turbulence diffusion timescale across this region is therefore shorter (Ward 1992).

The tables of TTW02 describe the various contributions to the torque in the case that both corotation and Lindblad resonances contribute in a fully isothermal disk, based on linear theory. We infer from these tables that the three-dimensional torque contribution due to Lindblad resonances is

$$\begin{eqnarray} \mathcal {T}^{\mathrm{L}}_{\rm{{TTW}}} &=&-\left(2.34 - 0.10\,\beta \right)\nonumber \\ & &\times \Sigma (a) \, \Omega ^{2}(a) \, a^{4} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2} \left(\frac{a}{H}\right)^{2}. \end{eqnarray} \tag{ 7 }$$

Menou & Goodman (2004) determined the three-dimensional torque for a locally isothermal disk (i.e., the temperature depends on r) on the basis of two-dimensional linear calculations with a softened potential (of the type introduced in Equation (5)) to model the effects of the third (vertical) dimension. Only Lindblad resonances are included. The softening length was chosen to be H. As they note, it is not clear what value is most appropriate within a range of that order. From Figure 3 of MG04, we infer that the Lindblad-only torque is given by

$$\begin{eqnarray} \mathcal {T}^{\mathrm{L}}_{\rm{{SLC}}} &=&-\left(0.80 -0.77\,\beta + 1.12 \,\zeta \right)\nonumber \\ & &\times \Sigma (a) \, \Omega ^{2}(a) \, a^{4} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2} \left(\frac{a}{H}\right)^{2}. \end{eqnarray} \tag{ 8 }$$

The subscript emphasizes that this result derives from linear calculations that employ a softened gravitational potential. We recall that ζ = −dln T/dln r. As is apparent, there is substantial disagreement between Equations (7) and (8) for ζ = 0, where they should agree, in particular for the coefficient of β. Such gradient terms can play an important role in describing the possible trapping of a planet in regions of density and temperature radial variations.

The effects of the density gradient (the β term) on the torque enter in two ways. With a negative radial gradient of the surface density (β>0), the density interior to the orbit of the planet is higher than the density exterior to it. As a result, the inner Lindblad torques are favored over outer ones, giving rise to a outward (positive) torque contribution. On the other hand, the negative density gradient slows the disk rotation, due to outward pressure forces (see Equation (3)). This effect acts in the same sense as a drag term would in causing inward migration. The degree of cancellation of these two effects is critical in determining the net torque contributions by a density gradient. In particular, MG04 concluded that the cancellation is nearly complete in two dimensions (i.e., when an unsoftened potential is employed), but not in three dimensions. However, TTW02 found nearly exact cancellation in three dimensions (see Equation (7)).

There are at least three possible explanations for the discrepancy between the Lindblad-only torque of Equations (7) and (8) with ζ = 0. One possibility is that the softened potential applied to a two-dimensional calculation, Equation (8), does not adequately represent the full effects of three-dimensional torques, Equation (7). The analytic calculation of TTW02 suggests that a more complicated disk response may occur in a three-dimensional disk, especially in the presence of gradients, and that this response may not be expressible through the use of a softened two-dimensional potential. A second possibility is that the local two-dimensional torque formula used by MG04, from Ward (1997), is not sufficiently accurate to determine the effects of gradients on the torque. This may be the case since the two-dimensional (no softened potential) Lindblad-only torque formula of TTW02 does not agree with that of Ward (1997) for a fully isothermal disk. The former has a coefficient of β equal to ≈1.5, while the latter has a much smaller value (nearly zero). A third possibility is that the three-dimensional Lindblad torque of Equation (7) does not provide an accurate description of the Lindblad torque that would occur in the absence of a corotation torque. Equation (7) does describe the Lindblad torque contribution in the presence of a coorbital torque. The Lindblad and corotation resonances can overlap and may affect each other, even in the linear theory. However, for β = 3/2, the corotation torque nearly vanishes in the three-dimensional calculation and the effects of resonance overlap can be largely ignored. (For β = 3/2, there is a small, ≈1%, nonvanishing corotation torque contribution caused by three-dimensional effects.) For β = 3/2, the torque predicted by the three-dimensional isothermal analysis, Equation (7), has a dimensionless torque that evaluates to −2.2, while the corresponding result for the softened two-dimensional potential in an isothermal disk in Equation (8) evaluates to +0.36. This large difference in predicted torques (and directions) leads us to conclude that the resonance overlap, which should be absent in this comparison, is unlikely to explain the discrepancy between the calculations.

3.2. Three-dimensional Torque in a Locally Isothermal Disk

In Figure 1, we plot the total torque on the planet as a function of density (left panel) and temperature (right panel) gradients (β and ζ, respectively) obtained from a subset of our simulations (filled circles, see figure caption for details) of an M_p = 3 × 10⁻⁶ M_s planet (1 M_E for a solar mass star) in a disk of fixed temperature and density at the orbital radius of the planet. In general, we consider disks with values of β in the range from −2 to 4 and values of ζ in the range from −2 to 2. But a few models had gradients outside of these ranges, as also indicated in Figure 1. The nonlinear simulation results we obtain fit well to the following torque equation:

$$\begin{eqnarray} \mathcal {T}&=& -\left(1.36 + 0.62\,\beta + 0.43\,\zeta \right)\nonumber \\ & &\times \Sigma (a) \, \Omega ^{2}(a) \, a^{4} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2} \left(\frac{a}{H}\right)^{2}. \end{eqnarray} \tag{ 9 }$$

The dependences on the planet-to-star mass ratio, M_p/M_s, and on disk scale height, H, at the planet's orbit (i.e., on the temperature at the orbit of the planet) are tested in Section 4 for fixed values of β and ζ. We find them to agree very well with the scalings in Equation (9) over a broad range of variations.

Zoom In Zoom Out Reset image size

If we subtract out the corotation contribution in the case of an isothermal disk, adopting the corotation torque expression given in TTW02, we obtain a Lindblad-only torque of

$$\begin{eqnarray} \mathcal {T}_{\rm{{ISO}}}^{\mathrm{L}}&=& -\left(2.34 - 0.02\,\beta \right)\nonumber \\ & &\times \Sigma (a) \, \Omega ^{2}(a) \, a^{4} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}} \right)^{2} \left(\frac{a}{H}\right)^{2}. \end{eqnarray} \tag{ 10 }$$

In this case, we infer a nearly exact cancellation between the two torque contributions by the density gradient, as described above. Not surprisingly, this result is close to the corresponding expression by TTW02, Equation (7).

The criterion for torque saturation of Ward (1992) is that the libration timescale τ_lib (over which the specific vorticity gradient is removed) is shorter than the viscous diffusion timescale across the coorbital zone τ_vd ≈ w²/ν (over which time the specific vorticity gradient is re-established). The relative velocity of the fluid on the librating orbit farthest from the planet's is ≈|3wΩ/4|, and thus τ_lib ≈ 16πa/(3wΩ), where Ω = Ω(a). The condition τ_lib ≲ τ_vd requires a viscosity such that

$$\begin{equation} \alpha \lesssim \left(\frac{a}{H}\right)^{7/2} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}} \right)^{3/2}. \end{equation} \tag{ 11 }$$

On this basis, we expect that the torques in Figure 1 involve unsaturated corotation torque contributions for the viscosity and sound speeds that we impose (α = 0.004 and H/a = 0.05 at r = a) on an M_p = 3 × 10⁻⁶ M_s planet, since the characteristic α required for saturation is only ≲10⁻⁴. In deriving the saturation condition above, we neglect a numerical factor of order unity multiplying the right-hand side of Equation (11) that depends on the precise value of w.

If the coorbital torque were partially saturated, it would introduce an additional dependence of the overall torque on planet mass. This dependence on planet mass would be inconsistent with the quadratic dependence predicted by linear theory. The torque we obtain in Equation (9) has a quadratic dependence on planet mass (see Section 4), as expected by linear theory. Therefore, our expectation that the torque is not saturated is consistent with its quadratic dependence on planet mass. We provide additional evidence for the lack of saturation below.

To test for saturation further, we run models with lower and higher viscosity, for a total variation of α of a factor 10. If saturation of the corotation resonance occurred, we would expect to observe a change in the total torque equal to 70% of the unsaturated torque, in a disk with zero surface density and temperature gradients (see Equations (6) and (7)). For such cases, however, we find only small differences (generally a few per cent) in total torques from those cases based on the standard viscosity that we apply. The results for somewhat lower and much higher viscosity are also plotted in Figure 1 (as empty squares and empty diamonds, respectively). The small differences in the figure (among different symbols) are not due to significant saturation effects, as indicated by the time behavior of the total torque illustrated in Figure 2, which is nearly constant over time, within a margin of 1% or less.

Zoom In Zoom Out Reset image size

If the coorbital torque were to saturate, it would generally decline (with some oscillations) from its unsaturated initial value to zero over several libration timescales of the coorbital region. Figure 2 plots three cases with turbulence parameter α = 0.002, represented as empty squares in the left panel of Figure 1. The torque is normalized to its average value between 100 and 200 orbits. Since τ_lib ≈ 170 orbital periods in these models, the torque evolution illustrated in Figure 2 spans a period of 12 libration timescales, or longer. There is no evidence of a substantial decline in the total torque and therefore no evidence of saturation.

In Section 4, we provide evidence that the torque behavior at even lower (α = 5 × 10⁻⁴) and higher viscosity (up to values of α = 0.05) is intrinsic to the torque density distribution, whose shape is subject to rather small changes with decreasing or increasing viscosity inside this range of the turbulence parameter, α. Therefore, we believe that our results apply to a disk where the corotation torques are unsaturated (i.e., fully effective).

The torques shown in Figure 1 disagree with the predictions, based on two-dimensional models, of Paardekooper & Papaloizou (2009) and Paardekooper et al. (2010). As we argue in Section 3.3, torques computed in a two-dimensional disk are susceptible to smoothing length issues and may not always reproduce three-dimensional torques. On the other hand, Masset & Casoli (2010) obtained a level of agreement with Equation (6) that is similar to ours and a similar variation (≲10%) of torques, as reported in our Figure 1, over the range of viscosities we consider and for zero disk gradients.

Equation (9) that we obtain agrees fairly well with the unsaturated torque expression in Equation (6) of TTW02, with ζ set to zero to match the radial isothermal assumption in the latter equation. Tanaka et al. (2002) provide the dependence of the torque contributions on the radial gradients of the surface density, disk thickness, and pressure in various tables. One might think that these tables can be used to infer a dependence of the linear torque on the radial temperature gradient. If these tables are used in this way, the resulting torque expression is in good agreement with Equation (9). However, it is not clear if the result is meaningful. The reason is as follows. Equation (11) of TTW02 is appropriate for a locally isothermal disk with radial temperature variations. But since there is a singularity in the term involving ∂Ω/∂z that they want to avoid, they set the radial temperature gradient to zero (so that ∂Ω/∂z = 0) in all subsequent equations, including the torque derivations. In doing so, they not only discard the ∂Ω/∂z term but also other terms, including a term that contributes to the coorbital torque. This coorbital torque term provides a contribution that is independent of the vortensity gradient (the usual coorbital torque contribution) and is due to the radial temperature gradient (see also Casoli & Masset 2009). Therefore, the effects of the temperature gradient on the linear three-dimensional torque cannot be readily deduced from TTW02.

The coefficient of the temperature gradient that we obtain in Equation (9) is significantly smaller than the corresponding coefficient in Equation (8), obtained by using a two-dimensional smoothed potential. Furthermore, the contrast between the results is even stronger when comparing the ratio of the gradient coefficients to the term that is independent of β and ζ in the respective equations (i.e., 1.12/0.8 versus 0.43/1.36). The slowing of migration (or trapping of planets) by means of the opacity variations, envisaged in MG04, relies on having a sufficiently steeply declining density profile and a sufficiently shallow temperature profile (see Equation (8)). However, the three-dimensional results (Equation (9) or (10)) suggest that this trapping mechanism does not occur under such conditions.

Equations (9) and (10) are derived under the requirement that surface density and temperature gradients (or, more precisely, the derivatives dln Σ/dln r and dln T/dln r) are roughly constant over a radial region of about 3H interior and exterior to the planet's orbit, where the torque is exerted (see Section 4). However, in Section 4 we will develop a formalism for computing torques that can be applied in the presence of mildly varying radial gradients (see Section 6) and of gap-opening planets (see Section 7).

3.3. A Two-dimensional Torque Formula for a Locally Isothermal Disk

Two-dimensional (r–ϕ) studies of disk–planet interactions often smooth the gravitational potential of the planet, Φ_p, over a length of order the disk thickness, H. This procedure is meant to account for the vertical stratification of disk material in proximity of the planet. However, as noted by MG04, it is unclear how well this approximation works since three-dimensional resonant torques depend on the vertical distribution of density in the disk (TTW02). Moreover, two-dimensional simulations of the Type I regime indicate that torques depend on the value of the softening length (e.g., Masset 2002; Nelson & Papaloizou 2004).

We perform two-dimensional simulations using settings and parameters described in Section 2 and applying a planet potential with a smoothing length equal to H, the same value of the smoothing length adopted by MG04. The potential is then $\Phi _{p}=-G M_{p}/\sqrt{S^{2}+H^{2}}$, where S is the distance from the planet and M_p = 3 × 10⁻⁶ M_s. Although the softening length is introduced as a proxy for the vertical disk thickness, it obviously affects the horizontal gradient of the gravitational potential as well. The torque density distribution of non-gap-opening planets peaks at a radial distance of about H from the planet's orbit (see Section 4). At those locations, the smoothed gravitational force, |∂Φ_p/∂S|, is $\sqrt{2}/4$ of that resulting from an unsoftened, point-mass potential.

The grid resolution in these calculations is Δr = a Δϕ = 1.5 × 10⁻³. We consider values of the disk gradients, β and ζ, in the range from −2 to 2. By fitting the results from the two-dimensional nonlinear models, we obtain a total torque that can be written as

$$\begin{eqnarray} \mathcal {T}_{\mathrm{2D}}&=& -\left(1.17 + 0.15\,\beta + 0.22\,\zeta \right)\nonumber \\ & &\times \Sigma (a) \, \Omega ^{2}(a) \, a^{4} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2} \left(\frac{a}{H}\right)^{2}. \end{eqnarray} \tag{ 12 }$$

We check the scaling of Equation (12) with disk thickness, H, by performing calculations for given values of β and ζ and varying the sound speed at r = a. The scaling with planet mass was already confirmed by previous two-dimensional studies (see, e.g., Nelson & Papaloizou 2004; Masset et al. 2006).

Clearly, there are important differences between our three-dimensional and two-dimensional results with smoothed potential, Equations (9) and (12), respectively. In fact, although there is not a wide discrepancy in total torque for the case of a spatially isothermal disk (ζ = 0) with constant surface density (β = 0), there are factors of 4 and 2 difference in the coefficients of parameters β and ζ, respectively.

We did not check how sensitive Equation (12) is to changes of the softening length and whether it is possible to reproduce more closely the three-dimensional formula (Equation (9)) with a different a value of this parameter. However, there are indications that the coefficients of β and ζ in Equation (12) are highly dependent on the softening length (S. Li & H. Li 2010, private communication).

The two-dimensional analytic linear torque of TTW02 produces a coefficient of β equal to 2.83 (see also the two-dimensional linear calculations by Korycansky & Pollack 1993), which is quite different from our value, although the torques are nearly identical for β = ζ = 0. However, TTW02 did not apply any smoothing in obtaining their result. Therefore, the likely explanation for this discrepancy is the smoothing length H that we adopt.

Waves carry torques until they are damped and deposit angular momentum in the disk. In order to gain some understanding of the radial distribution of torques, which is not provided by the integrated torque presented in Section 3, we analyze the azimuthally averaged torque density as a function of radius. This function measures the local strength of disk torques and should be insensitive to the details of wave damping that determine where the torque is exerted on the disk.

The torque distribution per unit disk mass as a function of radius, $d\mathcal {T}(r)/dM$, is defined such that

$$\begin{equation} \mathcal {T}=2 \pi \!\!\int ^{\infty }_{0}{\frac{d \mathcal {T}}{\it dM}(r) \,\Sigma (r)\, r\, {\it dr}}, \end{equation} \tag{ 13 }$$

where Σ(r) is the axisymmetric (i.e., azimuthally averaged) gas surface density. This quantity provides an important diagnostic for the nature of the disk–planet torques (DL08; Li et al. 2009). The theory of disk resonances (e.g., Meyer-Vernet & Sicardy 1987; Ward 1997) suggests that the torque distribution for a disk with smoothly varying properties should approximately follow

$$\begin{equation} \frac{d \mathcal {T}}{\it dM}(r)= \mathcal {F}\!\left(\frac{r-a}{H}, \beta, \zeta \right) \Omega ^{2}(a) \, a^{2}\! \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2}\! \left(\frac{a}{H}\right)^{4}\!, \end{equation} \tag{ 14 }$$

where $\mathcal {F}$ is a dimensionless function and H is evaluated at the planet's orbital radius, a. (There are some modifications to this expression due to resonance widths that we ignore.) Notice that, by using Equations (13) and (14) along with the approximation Σ(r) ≈ Σ(a) and performing the transformation r → r/H, one correctly finds from the integration that the total torque scales as (a/H)². In the past, function $\mathcal {F}$ has in effect been taken to be an inverse power law with distance from the planet that is modified close to the planet. The expression is based on the impulse approximation (Lin & Papaloizou 1986). We provide here a more accurate determination of $\mathcal {F}$.

In this section, we describe results for testing the universality of function $\mathcal {F}$ over a range of disk parameters by means of three-dimensional simulations. For this purpose, we fix surface density and temperature gradients and consider the case of a disk with β = 0.5 and ζ = 1. This situation corresponds to a steady-state viscous disk (accretion rate independent of radius) with constant aspect ratio, H/r, and constant viscosity parameter, α. Function $\mathcal {F}$ should then depend only on (r − a)/H, with these assumed values of β and ζ.

We test Equation (14) by considering simulations involving four different values for H/a and M_p/M_s and levels of viscosity in Figures 3–5, respectively. These figures show the results of $d\mathcal {T}(r)/{\it dM}$ obtained from simulations that are scaled such that the curves in the three figures should lie on top of each other, if the functional form of Equation (14) is correct. The figures demonstrate fairly good agreement with Equation (14) for variations of H/r. There are ∼30% variations (measured at the curve peaks) in the scaled values of $d\mathcal {T}(r)/{\it dM}$ (see Figure 3), while the unscaled torque densities differ by a factor of ∼20⁴. There is less than a 2% variation in the scaled $d\mathcal {T}(r)/{\it dM}$ for a factor of 20² change in torque density due to planet mass changes (see Figure 4). There is less than 5% change in the scaled $d\mathcal {T}(r)/{\it dM}$ for a factor of 100 change in α (see Figures 5 and 6). These results also add support to our claim in Section 3 that the coorbital disk torques are unsaturated. The reason is that there would be a decrease in the magnitude of net torque at about the 30% level at higher viscosity, if the coorbital torques were to switch from being saturated to being unsaturated (see Equations (9) and (10)). In particular, Figure 6 indicates that there is no significant change in total torque, over many libration timescales, since fractional variations stay at levels <0.5%. The slightly lower mass planet, M_p = 3 × 10⁻⁶ M_s (1 M_E), used in Section 3 (versus M_p = 5 × 10⁻⁶ M_s adopted here) is even less prone to saturation according to the criterion in Equation (11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The torque density distributions in, e.g., Figure 4 would produce a positive total torque, through Equation (13), for a sufficiently negative density gradient (β>0) since the multiplication by Σ would make the positive peak of the profile more positive and the negative peak less negative. However, as indicated in Equation (14), function $\mathcal {F}$ depends on the radial gradients of both density and temperature. For increasing values of |β| and |ζ|, the asymmetry between the positive and negative peaks becomes larger. As β and ζ increase, the positive peak becomes shallower whereas the negative peak deepens. This effect outweighs that produced by the surface density. This is illustrated in Figure 7, where $d\mathcal {T}(r)/dM$ is shown for various values of β (top left panel) and of ζ (bottom-left panel). In the right panels of the figure, we plot the cumulative torque, $\mathcal {T}(r)$, defined as the integral on the right-hand side of Equation (13) with integration limits from 0 to r (and Σ ∝ 1/r^β). As β and ζ increase, the magnitude of the total torque increases. The trend can persist for steeper surface densities, up to β ≈ 6, the largest radial gradient we consider. The total torque tends to decrease (i.e., it becomes less negative) for positive radial gradients of surface density and temperature, but relatively steep functions seem to be required to result in a positive torque (see bottom panels of Figures 7 and 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Analytic Approximations of Torque Density Distributions

Analytic expressions of function $\mathcal {F}$ in Equation (14) can be found by fitting data. For this purpose, we assume that $\mathcal {F}$ is a parametric function of β and ζ, of independent variable x = (r − a)/H(a), and has the following form:

$$\begin{eqnarray} \mathcal {F}(x,\beta,\zeta)&=&[p_{1}\,e^{-(x+p_{2})^{2}/p^{2}_{3}} +p_{4}\,e^{-(x-p_{5})^{2}/p^{2}_{6}}]\nonumber \\ & &\times \tanh {\left(p_{7}-p_{8}x\right)}, \end{eqnarray} \tag{ 15 }$$

and p_i = p_i(β, ζ), with i = 1, ..., 8, are the parameters resulting from the fit. This functional form is not based on a physical model of resonant torques, but rather the result of searching among simple combinations of functions whose parameters have direct physical interpretations. So, in Equation (15), parameters p₁ and p₄ are related to the amplitude of the positive and negative peaks, respectively, while parameters p₂ and p₅ are the distances of the peaks' positions from the planetary orbit. Parameters p₃ and p₆ regulate the width of the positive and negative portions of the function. The hyperbolic tangent affects the slope with which the profile transitions from positive to negative and the ratio p₇/p₈ gives the intercept of the function with the x-axis. Three examples of the fit are displayed in Figure 8, corresponding to gradients specified by pairs the (β, ζ) = (2, 2) in the top panels, (0.5, 1) in the middle panels, and (−2, − 2) in the bottom panels. In the figure, the rescaled torque density distribution $d\mathcal {T}(r)/dM$ (left panels) and the cumulative torque $\mathcal {T}(r)$ (right panels) from a simulations (solid dots) are compared against the fitting function (solid lines).

In Table 1, we give parameters p_i for a number of representations of function $\mathcal {F}$, corresponding to various parameter pairs (β, ζ). We find that expression (15) generally produces a very good or good fit to the distributions obtained from simulations. There are various ways to quantify how well Equation (15) fits the data, one of the most stringent ways is to compare cumulative torques generated by the fitting functions against those determined from simulations. Three such comparisons are plotted in the right panels of Figure 8, which have a level of agreement ranging from a few per cent (top and middle panels) to ≈10% (bottom panel). The approximations of function $\mathcal {F}$ by means of Equation (15), with parameters given in Table 1, produce asymptotic cumulative torques that deviate from simulation data by typically less than 10%.

Table 1. Fit Parameters for the Function $\mathcal {F}$ in Equation (15)

(β, ζ)	p1	p2	p3	p4	p5	p6	p7	p8
(−2, − 2)	0.0501524	0.849572	0.920263	0.0233445	1.20039	1.15790	0.145358	3.34263
(0, − 2)	0.0419354	0.984212	0.861833	0.0293424	1.03610	1.18903	−0.481275	4.38744
(1, − 2)	0.0379796	1.05644	0.826359	0.0328917	0.945701	1.21163	−0.936489	5.30156
(2, − 2)	0.0315764	1.15319	0.803355	0.0368841	0.760423	1.35752	−0.989595	4.54574
(4, − 2)	0.0202990	1.32162	0.713699	0.0464587	0.557390	1.48935	−1.11322	3.58864
(−2, − 1)	0.0466666	0.863280	0.966481	0.0264784	1.13846	1.12637	0.300767	3.29211
(0, − 1)	0.0382860	0.997974	0.910000	0.0328630	1.00086	1.15666	−0.171042	3.58431
(1, − 1)	0.0349835	1.05169	0.902117	0.0365165	0.922390	1.17324	−0.403456	3.31236
(2, − 1)	0.0309523	1.11777	0.889403	0.0405477	0.840493	1.22073	−0.593193	3.16310
(4, − 1)	0.0228084	1.34202	0.789789	0.0507717	0.629326	1.32866	−0.897774	3.07277
(−2, 0)	0.0421255	0.897070	1.01784	0.0304875	1.02268	1.14905	0.428988	2.89123
(0, 0)	0.0341462	1.04357	0.958885	0.0369814	0.933042	1.13253	0.0358698	3.23369
(1, 0)	0.0303427	1.11959	0.928609	0.0408822	0.861233	1.15412	−0.181823	3.07328
(2, 0)	0.0287186	1.21232	0.849933	0.0447468	0.802171	1.16904	−0.466765	3.59154
(4, 0)	0.0223466	1.32467	0.785512	0.0542731	0.646047	1.22339	−0.713550	2.69379
(−2, 1)	0.0342374	0.956221	1.11505	0.0352915	0.814786	1.25951	0.393368	2.10245
(0, 1)	0.0317793	1.06971	0.969339	0.0401954	0.931402	1.04558	0.221833	4.50245
(0.5, 1)	0.0297597	1.09770	0.938567	0.0421186	0.902328	1.03579	0.0981183	4.68108
(1.5, 1)	0.0270259	1.22111	0.891422	0.0467822	0.832732	1.06184	−0.226587	4.07528
(2, 1)	0.0255936	1.27569	0.867892	0.0489470	0.791497	1.06753	−0.358692	3.69761
(4, 1)	0.0178517	1.39810	0.912140	0.0595500	0.630116	1.16562	−0.487347	1.89963
(−2, 2)	0.0347398	0.864893	1.21887	0.0346082	1.00993	1.03900	0.823430	3.99184
(0, 2)	0.0274646	1.07401	1.04454	0.0435430	0.924349	0.975019	0.538068	5.37855
(1, 2)	0.0244322	1.20668	0.986380	0.0474291	0.866602	1.01086	0.323656	4.93965
(2, 2)	0.0211401	1.27928	0.987525	0.0518687	0.785682	1.05005	−0.134179	3.71919
(4, 2)	0.0184703	1.32365	1.04665	0.0666310	0.671786	1.01557	−0.446362	1.84684

Download table as:  ASCIITypeset image

Parameters in Table 1 can be used to generate an approximation of function $\mathcal {F}$ for a given pair (β, ζ). For example, one can proceed by selecting from the table four sets of parameters p_i corresponding to β₁, β₂, ζ₁, and ζ₂, such that β₁ ⩽ β < β₂ and ζ₁ ⩽ ζ < ζ₂. These parameters p_i can then be interpolated at (β, ζ), and the interpolated values can be used in Equation (15) to obtain function $\mathcal {F}$. We use bi-linear interpolations to generate torque density distributions $d\mathcal {T}(r)/dM$ every Δβ = Δζ = 0.2, which were then integrated according to Equation (13) to obtain the total torque as a function of the surface density and temperature gradients. By fitting these data, the coefficients of β and ζ in Equation (9) can be recovered with good agreement, within a margin of 10%.

4.2. A Migration Test Case

The formulae presented above in this section can be readily used to describe planet's migration, as indicated by the test case discussed here. We set up a three-dimensional disk model with a planet whose mass is M_p = 10⁻⁵ M_s (or 3 M_E if M_s = 1 M_☉). The disk is isothermal (ζ = 0) and has an initial surface density corresponding to a parameter β = 2. The slope of the azimuthally averaged Σ(r) is maintained throughout the disk evolution by imposing a kinematic viscosity ν ∝ r^β, so that the disk is in an approximate steady state.

The planet is allowed to migrate under the action of disk torques after 30 orbital periods. We adopt a treatment to approximately account for the axisymmetric effects of disk self-gravity by forcing the planet to rotate at (approximately) the same rate as the disk (see Appendix C of DL08). The grid rotates about its origin, O, at the same angular velocity as the planet (for details, see D'Angelo et al. 2005). Therefore, the planet moves only radially through the grid because of orbital migration. Gravitational forces exerted on the planet are computed via direct summation of elemental forces and the equations of motion of the planet are integrated by means of a Bulirsch–Stoer algorithm, which uses an adaptive time step control to achieve a relative accuracy of 10⁻⁵.

The evolution of the semimajor axis is illustrated in the top panel of Figure 9 as a thick solid line. We use coefficients p_i in Table 1, for the pair (β, ζ) = (2, 0), to obtain an approximation of function $\mathcal {F}$ (Equation (15)) and then calculate the torque density distribution $d\mathcal {T}(r)/{\it dM}$ (Equation (14)). The total torque $\mathcal {T}(a)$, as a function of the planet's semimajor axis a, is evaluated through the integral in Equation (13). Conservation of orbital angular momentum imposes that $da/dt=2\,\mathcal {T}(a)/(\mbox{$M_{p}$}\,a\,\Omega _{\mathrm{K}})$, which can be written more explicitly as

$$\begin{eqnarray} \frac{da}{dt}&=&4\pi \,\Omega _{\mathrm{K}}(a)\! \left(\frac{a}{\mbox{$M_{s}$}}\right)\! \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)\!\left(\frac{a}{H}\right)^{4}\nonumber \\ & &\times \int ^{r_{2}}_{r_{1}}{\!\! \mathcal {F}\!\left(\frac{r-a}{H}, \beta, \zeta \right) \!\Sigma (r)\, r\, dr}. \end{eqnarray} \tag{ 16 }$$

The integration limits in the above equation are such that r₂ − a ≫ H(a) and a − r₁ ≫ H(a).

Zoom In Zoom Out Reset image size

A numerical solution of Equation (16) is shown as a thin solid line in the top panel of Figure 9. The average relative difference between the orbital evolution obtained from the three-dimensional model and that simulated with the formulae presented above is at the 10⁻³ level for 10% variations of the planet's semimajor axis (see Figure 9, bottom panel).

Torque density distributions are often used to describe the tidal interaction between a disk and embedded bodies. In this section, we describe an application of the formulae provided in Section 4 to protoplanetary disks. Similar applications to other contexts involving astrophysical disks may employ an analogous strategy.

The synthesis of semimajor axis distributions of extrasolar planetary systems requires solving for the evolution of the protoplanetary disk and, simultaneously, for the gravitational interaction occurring between the planet(s) and the disk. Since these calculations involve evolution timescales equal to typical lifetimes of protoplanetary disks (millions of years) and length scales of hundreds of AU, they need to rely on one-dimensional models, as a consequence of the limitations of current computer capabilities.

The evolution of a thin, axisymmetric and viscous disk that interacts with one (or more) embedded planet(s) is described by the following equation (e.g., Lin & Papaloizou 1986):

$$\begin{equation} \pi r\frac{\partial }{\partial t}(\Sigma +\Sigma _{\mathrm{pe}})= \frac{\partial }{\partial r}\!\left\lbrace \frac{1}{r\,\Omega }\frac{\partial }{\partial r}\left[ \mathcal {T}_{\nu }(r)-\mathcal {T}(r)\right]\right\rbrace, \end{equation} \tag{ 17 }$$

where $\mathcal {T}_{\nu }(r)=-2\pi r^{3}\nu \Sigma \,d\Omega /dr =3\pi r^{2}\nu \Sigma \Omega$ (Lynden-Bell & Pringle 1974) is the disk viscous torque and $-\mathcal {T}(r)$ is the tidal torque exerted by the planet(s) at radius r. In Equation (17), it is assumed that the disk's rotation rate is unperturbed, i.e., Ω = Ω_K. Since one can write the second term in square brackets as $\partial \mathcal {T}/\partial r=2\pi r\Sigma \,\partial \mathcal {T}/\partial M$, Equation (17) becomes

$$\begin{equation} r\frac{\partial }{\partial t}(\Sigma +\Sigma _{\mathrm{pe}})= \frac{\partial }{\partial r}\!\left[ 3\sqrt{r}\frac{\partial }{\partial r}\!(\nu \Sigma \sqrt{r}) -\frac{2\,\Sigma }{\Omega }\frac{\partial \mathcal {T}}{\partial M}(r)\right]. \end{equation} \tag{ 18 }$$

In case of a multiple-planet system, the second term on the right-hand side involves a torque density distribution for each planet, and thus the summation $\partial \mathcal {T}/\partial M=\sum _{n}{\partial \mathcal {T}_{n}/\partial M}$.

The left-hand side of Equations (17) and (18) includes the term ∂Σ_pe/∂t, the mass loss flux from the disk. Here, we assume a simple model of gas removal due to photoevaporation by the central star that includes effects of UV photons

$$\begin{equation} \frac{\partial \Sigma _{\mathrm{pe}}}{\partial t}= 3.7\times 10^{-13} \left(\frac{10\,\mbox{AU}}{r}\right)^{5/2}\, \left(\frac{M_{\odot }}{\mbox{AU}^{2}\,\mathrm{yr}}\right), \end{equation} \tag{ 19 }$$

for r > 10 AU and ∂Σ_pe/∂t = 0 otherwise. Equation (19) applies to a disk surrounding a 1 M_☉ central star that emits ionizing photons at a rate of 10⁴¹ s⁻¹ (see, e.g., Martin et al. 2007 and references therein).

Equation (18) is solved by means of an implicit scheme, second-order accurate in time, over a grid that extends from 0.01 AU to 1850 AU. The grid is composed of 10,000 elements and uses a constant logarithmic spacing. Torque density distributions are obtained from Equations (14) and (15), at each time step, according to the interpolation procedure outlined in the last paragraph of Section 4.1. The parameter β required to construct function $\mathcal {F}$ is computed, every time step, as an average slope

$$\begin{equation} \bar{\beta }=-\frac{1}{(r_{2}-r_{1})}\int _{r_{1}}^{r_{2}} \!\frac{r}{\Sigma }\frac{d\Sigma }{dr}\,dr, \end{equation} \tag{ 20 }$$

in which the integration limits are r_2,1 = a ± 3H(a). The orbital radius of the planet is computed by numerically solving Equation (16).

The examples presented here are based on two disk configurations, both intended to reproduce the initial mass distribution of gas in a Minimum-Mass Solar Nebula. In the first, the initial disk surface density is taken from Davis (2005). At 1 AU, Σ ≈ 1100 g cm⁻² at time t = 0. Within 4 AU, $\Sigma \propto 1/\sqrt{r}$, followed by an exponential decay at larger radii. The disk contains a total mass of ∼0.02 M_☉ inside of about 40 AU. We assume that H/r is constant (i.e., the disk temperature T ∝ 1/r) and that the kinematic viscosity is ν = α H²Ω, where α is also a constant (hence $\nu \propto \sqrt{r}$). In the second configuration, Σ = 1700 g cm⁻²(1 AU/r)^3/2 at time t = 0 (Hayashi 1981). An exponential decay is applied beyond about 40 AU so that the disk mass within ∼70 AU is ∼0.02 M_☉. The aspect ratio of the disk is H/r ∝ r^1/4 (and thus $T\propto 1/\sqrt{r}$). In this case, we assume that $\alpha \propto \sqrt{r}$, hence the kinematic viscosity ν ∝ r^3/2.

Figure 10 shows the surface density at various times (see figure's caption for details). The top panels refer to the first disk's configuration, whereas the bottom panels refer to the second configuration. Each panel corresponds to a given value of α at 1 AU (but recall that α is constant in the first configuration and only varies with radius in the second), which produces initial accretion rates on the star between 10⁻⁸ M_☉ yr⁻¹ and 10⁻⁷ M_☉ yr⁻¹.

Zoom In Zoom Out Reset image size

The bottom-left panel of Figure 10 indicates that a density gap starts to form along the planet's orbit at late evolution times. This is due to the small viscosity and disk thickness at small radii. A condition for gap formation is that the symmetrized one-sided torque, i.e., the torque exerted by the planet on disk's gas interior and exterior of its orbit, exceeds the viscous torque (e.g., Lin & Papaloizou 1986). The one-sided torque is of order $(a/H)|\mathcal {T}|$. This can be found by using Equation (14) and approximating the integral $\int ^{\infty }_{a}\!|\mathcal {F}_{S}(x)|\Sigma (r) r dr$ in Equation (13) as $a H \Sigma (a)\!\int ^{\infty }_{0}\!|\mathcal {F}_{S}(x)| dx$, where x = (r −  a)/H(a) and $\mathcal {F}_{S}(x)=\left[\mathcal {F}(x)-\mathcal {F}(-x)\right]/2$ is a symmetrized form of function $\mathcal {F}(x)$. Therefore, the condition for gap formation becomes

$$\begin{equation} \xi \,\Sigma (a)\,\Omega ^{2}(a) a^{4}\! \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2}\left(\frac{a}{H}\right)^{3} \gtrsim \mathcal {T}_{\nu }(a), \end{equation} \tag{ 21 }$$

with $\xi =2\pi \!\!\int ^{\infty }_{0}\!|\mathcal {F}_{S}(x)| dx$. Recalling the definition of viscous torque given above, one obtains the condition

$$\begin{equation} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2}\gtrsim \frac{3\pi \alpha }{\xi } \left(\frac{H}{a}\right)^{5}, \end{equation} \tag{ 22 }$$

where the factor ξ is typically of order unity. In the second disk configuration, this condition is satisfied in disk regions interior of a few tenths of AU, but it is never fulfilled in the first disk configuration.

The planet evolution tracks illustrated in Figure 11 assume that a planetary embryo grows at an oligarchic rate, dM_p/dt ∝ M_p^2/3, from 0.1 M_E to 5 M_E over a ∼10⁶ yr period. All tracks start at 13 AU and are computed until disk gas within tens of AU has almost entirely dissipated and da/dt ≈ 0. The left and right panels of the figure refer, respectively, to the first and second disk configuration.

Zoom In Zoom Out Reset image size

Observations show that super-Earths, i.e., planets of several Earth masses, are also found within a few tenths of an AU from their stars, a likely result of orbital migration. This occurrence is predicted in the case of the second disk configuration, as indicated by the solid-line track in the right panel of Figure 11. The first disk configuration would also predict a similar outcome for a somewhat smaller turbulence parameter (α ≈ 0.001) or for different initial conditions of planet's mass and orbital radius (see Figure 12).

Zoom In Zoom Out Reset image size

The disk configurations introduced above are also used to produce the results illustrated in Figure 12. In both cases, a turbulent viscosity parameter α = 0.01, at 1 AU, is applied. The disk is allowed to evolve without any planet until time t₀. At time t₀, we assume that a planet of mass M_p is located at the orbital radius a₀. The disk–planet system is then allowed to evolve until negligible amounts of gas remain in the disk and da/dt drops virtually to zero. The orbital radius attained by the planet at the end of the calculation is indicated as a_∞, its asymptotic value. Figure 12 shows a_∞ as a function of a₀ for cases in which t₀ = 0 (left panels), 5 × 10⁵ yr (middle panels), and 10⁶ yr (right panels). Top and bottom panels refer, respectively, to the first and second disk's initial configuration. In each panel, different lines styles correspond to models with different planet masses (which is taken to be constant): M_p = 0.3 M_E (short-dashed line), 1 M_E (long-dashed line), 3 M_E (dot-dashed line), and 5 M_E (solid line). Clearly, initial conditions play a key role in the determination of semimajor axis distributions of planet populations.

As discussed in Section 1, a number of studies have looked into the possibility of stopping migration by trapping the planet at disk locations where the radial gradients of surface density and/or temperature may be strong and vary over short radial distances, as it may occur due to rapid changes in disk opacity or turbulent viscosity.

The results developed in Section 4 were derived under the condition that surface density and temperature slopes, β and ζ, are nearly constant over a radial region of about 3H(a) on either side of the planet's orbit, where most of the torque is produced. Here, we wish to show that the same formulae are applicable when disk radial gradients vary somewhat over this region. This is because torque density distributions respond smoothly to mildly varying gradients and can be approximated in terms of mean slopes.

We set up a three-dimensional model with a surface density Σ whose azimuthal average is illustrated in the left panel of Figure 13. The profile is such that Σ ∝ 1/r² at r ≪ 1 and constant at r ≫ 1, as also indicated by the long-dashed and short-dashed lines, respectively. This setup mimics the conditions of the opacity jumps in MG04 for a locally isothermal disk. Details of the slope β as a function of radius, in the transition region, are shown in the inset. This profile is maintained stable in time by imposing a kinematic viscosity ν ∝ 1/Σ. The temperature is constant throughout. In the right panel, rescaled torque density distributions $d\mathcal {T}(r)/dM$ are plotted when the planet (M_p = 10⁻⁵ M_s) is located at three different radial positions, indicated in the legend.

Zoom In Zoom Out Reset image size

The torque density distributions in Figure 13 exhibit a behavior reminiscent of that shown in the top left panel of Figure 7 where different, but radially constant, slopes β are imposed. This behavior is also seen by comparing the dot-dashed and short-dashed curves in the right panel with the thin solid curves marked by a filled square and diamond, respectively. These are torque densities per unit disk mass obtained from Equation (15) and parameters p_i (from Table 1) interpolated at slopes $\bar{\beta }\approx 0.2$ (filled square) and $\bar{\beta }\approx 1.8$ (filled diamond). Both slopes are obtained as mean values according to Equation (20). The torque density distribution produced when the planet is located at radius r = 1 (long-dashed line) has an intermediate behavior and can be approximated by using a mean slope $\bar{\beta }=1$, again resulting from Equation (20).

Therefore, if a disk has radial gradients that are mildly varying over a region of a few times H, one can still apply the expressions derived in Section 4 after having introduced appropriate mean slopes. One important aspect to appreciate is the small variation of the torque density distribution $d\mathcal {T}(r)/dM$ as the density gradient changes. Even though β varies by factors of order unity, the peaks of $d\mathcal {T}(r)/dM$ are only affected at ∼10% levels (see also Figure 7, top left panel). The results do not show any indication of a migration trapping.

Figure 14 shows, as a solid line, the migration track of an M_p = 10⁻⁵ M_s planet across the transition region of Figure 13. The migration track is obtained by integrating Equation (16) (Σ is that illustrated in Figure 13) and using Equation (20) to evaluate the mean slope $\bar{\beta }$ that characterizes function $\mathcal {F}$. For comparison purposes, the figure also shows the tracks that the same planet would follow if Σ was constant (short-dashed line) and Σ ∝ 1/r² (long-dashed line). The solid and short-dashed lines are similar until the density gradient $\bar{\beta }$ becomes greater than zero. Incidentally, the slope of the long-dashed line is similar to the other two initially since the reduced density value compensates for the larger density gradient (see Equation (9)).

Zoom In Zoom Out Reset image size

The rapid decay at small radii of the planet's semimajor axis seen in the solid and long-dashed curves of Figure 14 is what one would expect from, e.g., Equation (9) applied to a stationary disk: $|\mathcal {T}(a)|\propto a^{(\zeta -\beta)}$, hence

$$\begin{equation} \left|\frac{da}{dt}\right|\propto a^{(\zeta -\beta +1/2)}. \end{equation} \tag{ 23 }$$

The solid curve is therefore expected to steepen (|da/dt| diverges) as the semimajor axis decreases, whereas the short-dashed curve is expected to flatten (|da/dt| tends to zero) as a decreases.

Density transitions occurring at the edges of a "dead zone" may be quite abrupt and the slope β may change by large amounts over radial distances of order H (e.g., Matsumura et al. 2007). The analysis provided here may not apply to those sharp transitions, whose effects on torques would require a more detailed analysis.

The torque density distributions given in Section 4 were derived in a regime in which the mass density distribution along the planet's orbit remains basically unperturbed. In general, this condition ceases to hold when tidal torques exceed viscous torques (see Equation (22)), which is roughly expressed by the inequality (see Section 5)

$$\begin{equation} \left(\frac{\xi }{3\pi \alpha }\right)^{1/2}\left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right) \left(\frac{a}{H}\right)^{5/2}\gtrsim 1. \end{equation} \tag{ 24 }$$

For standard disk parameters H/a ≈ 0.05, α of a few times 10⁻³, and $\sqrt{\xi }$ of order unity, Equation (24) requires a planet-to-star mass ratio M_p/M_s ≳ 10⁻⁴. Gap formation can induce average variations in the surface density, relative to the unperturbed state, of two orders of magnitudes when the planet mass grows from one Earth mass to one Jupiter mass. However, the corresponding variation in the strength of the peaks of the rescaled $d\mathcal {T}(r)/dM$ is only a factor between two and three, as shown in Figure 15. This figure displays torque density distributions produced by planets of various masses, for cases in which R_H < H (top panel) and R_H > H (bottom panel). Notice that the radial positions of the peaks are located at r − a ≈ ±R_H, when R_H > H (see also DL08).

Zoom In Zoom Out Reset image size

The thin solid curve with filled circles in the top panel of Figure 15 is function $\mathcal {F}$ (Equation (15)) of variable x = (r − a)/H(a) and parameters p_i given in Table 1, for disk gradients corresponding to β = 0.5 and ζ = 1. The same function, but divided by 2 and of variable x = (r − a)/R_H, is plotted in the bottom panel. The peak strength of the rescaled $d\mathcal {T}(r)/dM$ shows only small changes for R_H > H. (For R_H ≫ H, the torque density scalings can be different from those in Equation (14)). The factor two is given by the left-hand side of Equation (24) for R_H ≈ H and a more accurate determination of $\sqrt{\xi }$ (see definition in Section 5).

For a gap-opening planet, the details of the torque density distribution well within the gap region (|r − a| ≲ R_H) are not very important for computing torques in the framework of one-dimensional disk models (of the type discussed in Section 5) since those portions of $d\mathcal {T}(r)/dM$ are multiplied by a surface density that is relatively small (because of mass depletion). For reference purposes, in the bottom panel of Figure 15 we also plot as a thin solid curve with empty circles the torque density per unit disk mass that results from the impulse approximation, as implemented by Veras & Armitage (2004).

We obtained an expression, Equation (9), for the disk–planet torque in three dimensions for a locally isothermal disk that contains radial density and temperature variations, based on numerical simulations. This expression was derived for density and temperature gradients that are fairly constant over the torque producing region of the disk, which extends roughly over a radial distance of 3H inside and outside the orbit of the planet. (see, e.g., Figures 3–5). In the case of a fully isothermal disk, the torque expression agrees well with previous analytic, linear results obtained by TTW02 (Equation (6)). The analysis presented by TTW02 cannot be used to deduce the dependence of the torque on the temperature gradient (see discussion in Section 3.2), which we provided in Equation (9). We tested the torque expression for a range of disk densities, temperatures, viscosities, and planet masses and found good agreement. The results agree well with the expected scalings of linear theory with planet mass and gas sound speed (see Figures 3 and 4 for which linear theory predicts the curves overlap in each figure). Although these torque scalings agree with linear theory, we cannot claim that Equation (9) with a nonzero temperature gradient agrees with linear theory because, as noted above, there is no expression available for the linear torque in a three-dimensional locally isothermal disk. These results suggest that substantially slowing or halting planet migration by means of smoothly varying but large negative density gradients and small temperature gradients in a locally isothermal disk is difficult. In addition, considerably reducing or stopping migration in a locally isothermal disk by means of opacity variations (e.g., MG04) or shadowing effects (e.g., Jang-Condell & Sasselov 2005) appears unlikely (see Section 6).

We also obtained an expression for the torque density per unit disk mass (Equations (14) and (15)). We tested the range of applicability of this expression for a range of disk temperatures, planet masses, and viscosities and found good agreement with the expected scaling predicted by linear theory. This expression can be applied to one-dimensional simulations of disk–planet interactions, such as, but not limited to, the cases of super-Earth evolution discussed in Section 5. The torque density expression commonly used is based on the impulse approximation and contains a long power-law tail that extends far from the planet. Ad hoc modifications have been made to make it more physically reasonable (e.g., Veras & Armitage 2004). The expression we obtained is more accurate and overcomes that difficulty. The formalism presented in Section 4 can also be employed for computing torques in the presence of mildly varying disk gradients (see Figure 13) and of gap-opening planets (see Figure 15). Moreover, applications to other contexts involving tidal interactions between astrophysical disks and embedded objects may also be possible.

The Lindblad-only torque formula of MG04, Equation (8), based on the torque density of Ward (1997) modified for a planet potential smoothed over a distance H, results in a positive torque (outward migration) for a steep enough and radially decreasing surface density. We notice that, in the fully isothermal case, there is disagreement between the two-dimensional (no softened potential) formula of  TTW02 (for Lindblad torques only) and the results of Ward (1997). At the end of Section 3.1, we discussed some possible reasons for the discrepancy. But, we cannot be certain at this point what the reason is. In addition, we found that corotation torques play a significant role and, therefore, the use of a Lindblad-only torque formula may not be appropriate in the regimes that we investigated.

For zero radial density and temperature gradients in the disk, the torques resulting from nonlinear calculations of a two-dimensional disk, with a planet potential smoothed over a distance H, are in fair agreement with the three-dimensional results. However, there is substantial quantitative disagreement in the effects of these gradients on the torques (compare Equations (9) and (12)). One possible source of this disagreement may be due to the effects of the softened two-dimensional potential in modifying the torques in the presence of disk gradients.

The analysis presented here is limited by the range of the parameter space we considered, involving a turbulent viscosity 5 × 10⁻⁴ ≲ α ≲ 0.05. In addition, it is limited by the assumption that the disk thermodynamics behaves in a locally isothermal manner. The results may not hold under conditions very different from those we assumed. For example, if the gas behaves somewhat adiabatically, the torque can be modified and may even result in outward migration (e.g., Baruteau & Masset 2008; Paardekooper & Papaloizou 2008; Kley et al. 2009). Under such conditions, nonlinear effects may be more important than they are in the locally isothermal case.

In weakly turbulent disk regions (α ≲ 10⁻⁴), such as in "dead zones," we expect the disk response to a planet to be dominated by the highly nonlinear effects of waves launched at Lindblad resonances resulting in a complex behavior, including the formation of vortices (Li et al. 2009; Yu et al. 2010). Such effects produce strong changes in the underlying gas distribution near a low-mass (∼5 M_E) planet that in turn affect its migration rate. Type I planet migration can cease and a vortex dominated, jumpy form of migration can prevail. The linear torque formulae presented above are inapplicable in such a regime. For α ≳ 5 × 10⁻⁴, the viscosity regime we considered, the nonlinear feedback effects of shocks on migration subside due to turbulent diffusion that washes out these sharp features in the gas density (Li et al. 2009). However, these simulations were performed in two dimensions and it would be useful to understand how three-dimensional effects could alter this form of migration. In particular, the formation of vortices at gap edges may operate differently in three dimensions.

In highly turbulent disk regions (α ≳ 0.1), there may be other types of changes in the properties of the flow that can potentially affect disk–planet interactions. Preliminary results suggest that the shape of the torque density distribution may also be altered. Therefore, in such disk regions, the torque acting on low-mass planets may not be represented by the formulae presented here.

Paardekooper & Papaloizou (2009) performed linear and nonlinear two-dimensional calculations of torques acting on low-mass planets in fully isothermal disks and found that nonlinearities can develop, depending on the level of viscosity. For an M_p ∼ 10⁻⁵ M_s planet embedded in a disk with H/a = 0.05, they found agreement between linear and nonlinear calculations only at very high viscosity, beyond α ∼ 0.1. Those two-dimensional results, however, are highly sensitive to the value of the smoothing length applied in the planet's potential. As discussed above, in the fully isothermal case, our three-dimensional torques are in good agreement, including the effects of density gradients, with the linear theory of TTW02 at much smaller viscosity levels, α ∼ 0.001. In the more general locally isothermal case, we found that the overall scalings of the torque with planet mass and gas sound speed agree with the predictions of linear theory. Nonlinear feedback effects due to partial saturation of the corotation torque would in general give rise to different scalings with planet mass (e.g., Ogilvie & Lubow 2003). Such effects are not expected to be present in the parameter space we investigated (see Equation (11)), as we also verified in simulations (see Figures 2 and 6). Other nonlinear effects are possible, especially in other regimes as mentioned above, but we did not see any evidence for them to have a significant impact in these three-dimensional simulations.

The analysis provided here may not be applicable to cases where substantially large variations of density and temperature gradients occur within the disk zone where torques are produced, i.e., a few times H or less in radius. Such situations may arise near the edges of "dead zones" (e.g., Matsumura et al. 2007), although there are limits on the slope of these edges imposed by disk instabilities (e.g., Yang & Menou 2010). It would be desirable to explore such situations along the lines of the approach taken in this paper.

We benefitted from discussions with Frédéric Masset and Kristen Menou at the KITP program entitled "The Theory and Observation of Exoplanets." We are grateful to Gordon Ogilvie and Jim Pringle for helpful discussions and useful feedback on this work. We thank Hui Li and Shengtai Li for running some two-dimensional models and for informative discussions. We acknowledge support from NASA Origins of Solar Systems Program grants NNX08AH82G (G.D.) and NNX07AI72G (S.L. and G.D.). G.D. was also supported in part by the National Science Foundation under grant No. NSF PHY05-51164. S.L. acknowledges visitor support from the IoA and Newton Institute at Cambridge University and many beneficial discussions at the Newton Institute program "Dynamics of Discs and Planets." Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

In numerical simulations of the kind discussed in this paper, there is often much emphasis placed on the form of the gravitational potential of the planet. In a context in which disk self-gravity is neglected and the "effective" radius of the planet is much smaller than the smallest length scale resolved numerically, this potential is that of a dimensionless mass M_p, i.e., Φ_p = −G M_p/S, where S is the distance from said mass. Since this potential is singular at the position of the planet and the planet lies on the computational grid, the singularity could generate numerical instabilities because the source terms of the momenta equations could involve the calculation of a diverging gravitational force. (Notice that there is no problem with employing a point-mass potential for the gravity field generated by the star, when it is located outside of the computational grid.) Therefore, the singularity is usually removed by introducing a modified potential as in Equation (5). In this case, the potential is softened over a given length, which is either a fraction of the planet's Hill radius, R_H, or on the order of the local (hydrostatic) disk scale height, H(a). The second option is meaningful only in a two-dimensional (r–ϕ) geometry, where the smoothing length is used, in effect, as a proxy for the third (vertical) dimension of the disk. We choose the first option and adopt a softening length that is a tenth of R_H.

The deviation of the softened potential in Equation (5) from the point-mass potential, normalized to the point-mass potential, is ≈(ε²/2)(R_H/S)², which is ⩽12% for S/R_H ⩾ 2 ε and ⩽2% for S/R_H ⩾ 5 ε. We check that our results are largely insensitive to the choice of ε = 0.1 by performing calculations for an M_p = 5 × 10⁻⁶ M_s planet in a disk, whose viscosity parameter is α = 0.005, and using softening lengths (εR_H) with ε = 0.05 and 0.01. The resulting net torques agree at about the 10⁻³ level, implying that they are basically converged with respect to the parameter ε. We perform a similar experiment with an M_p = 10⁻³ M_s planet, using ε = 0.05, and found discrepancies at the 10⁻² level. Torque density distributions are plotted in Figure 16 (top panels) for the M_p = 5 × 10⁻⁶ M_s planet (left) and the M_p = 10⁻³ M_s planet (right), with values of ε indicated in the legends. Given the level of agreement, profiles corresponding to different softening length cannot be distinguished from one another, or barely so. The differences of $d\mathcal {T}(r)/dM$, taking cases with closest values of ε (see the figure's caption for details), are shown beneath the torque density distributions.

Zoom In Zoom Out Reset image size

The outcome illustrated in the top panels of Figure 16 can be explained with the property of the relative distance between the planet and the locations where momenta and gravitational potential are computed. The planet position is fixed on the (rotating) grid at a "corner" of a volume element of the grid (where eight such elements are in contact). Momenta are computed at the geometrical center of each of the six faces of a volume element, whereas the gravitational potential is computed at the geometrical center of the volume element. Therefore, there is always a buffer distance of $\sqrt{3}L/2$ (neglecting curvature), where L is the average grid spacing, which prevents Φ_p from diverging, even if ε = 0. This argument also applies when the planet migrates since, by virtue of the time-varying rotation rate of the reference frame (see Section 4.2), it only moves along the radius, R, but not along the azimuth, ϕ. Hence, there still exists a minimum buffer distance of $\sqrt{2}L/2$ between the planet position and the center of a volume element. (Notice that, on a two-dimensional grid, these buffer distances become shorter and equal to $\sqrt{2}L/2$ and L/2, respectively.)

Numerical resolution is another potential issue that may have an impact on simulation results. In the bottom panels of Figure 16, we test for convergence of torque density distributions with respect to numerical resolution. As mentioned in Section 2.2, the average linear resolution in the disk region of interest (|R − a| ⩽ 4 H) applied in this study ranges from 2 × 10⁻³a to 5 × 10⁻³a. These grid systems contain typically from 10 to over 40 millions of grid elements. The torque density distributions shown in Figure 16 (bottom panels) are calculated using average linear resolutions of 2 × 10⁻³ and 10⁻³ (left) and of 4 × 10⁻³ and 2 × 10⁻³ (right). The higher resolution cases in these plots are based on grid systems that have nearly 110 (left) and 90 (right) million grid elements. The profiles lie on top of each other and, again, can be hardly distinguished. The differences are displayed beneath the plots of $d\mathcal {T}(r)/{\it dM}$, showing agreement at a ∼10⁻² level or better.