THE IMPORTANCE OF DISK STRUCTURE IN STALLING TYPE I MIGRATION

New high-precision radial velocity and transit surveys allow observers to detect low-mass and modest-size planets, particularly when those planets with short-period orbits. Evidence is emerging that these close-in super-Earth mass planets are significantly more ubiquitous than their higher mass cousins, the hot Jupiters (Mayor et al. 2009, 2011; Howard et al. 2010).

Most theoretical models of hot Jupiter formation assume that planets form far from their parent stars and migrate inward, either by planet–planet scattering (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997; Murray et al. 1998), secular interactions such as the Kozai resonance (Wu & Murray 2003), secular chaos (Nagasawa et al. 2008; Marchi et al. 2009; Wu & Lithwick 2011), or through tidal interactions with a gaseous circumstellar disk (Lin et al. 1996). This migration is required for hot Jupiters because it is difficult for them to form in situ. The gas disk, which must be present to supply material to these planets, damps the eccentricity of planetesimals and embryos. This means that the growth of a protoplanetary embryo, on a circular orbit, stalls once it acquires an isolation mass (M_iso); i.e., it has captured all solids within a feeding zone centered around its orbital radius r. The width of the feeding zone is approximately ∼10r_Hill where the Hill's radius for a planet of mass M_p around a star of mass M_* is r_Hill ≡ r[M_p/(3M_*)]^1/3 (Kokubo & Ida 2002). This leads to an isolation mass in a disk with the local surface density of dust and planetesimals (Σ_d) is M_iso ∼ (Σ_da²/M_*)^3/2M_*.

If we adopt a value for Σ_d inferred from the empirical minimum mass nebula (MMN) model (Hayashi 1981, see below), M_iso would be a fraction of Martian mass at 1AU and a few times Earth mass (M_⊕) at 5AU (Ida & Lin 2004). This maximum mass increases with the distance from the central star as there are more solids available at large r (Lissauer 1987). In the conventional core-accretion scenario, efficient gas accretion onto gas giants needs to be preceded by the emergence of cores with masses in excess of ∼10 M_⊕ (Pollack et al. 1996). Progenitor cores with such large masses can only be assembled beyond the snow line in an MMN.

The same argument also suggests that it is difficult to form the rich population of known close-in (within ∼0.1AU) super-Earths (a few to 15 M_⊕) in situ. In order for there to be an adequate amount of building block material to assemble these planets, the magnitude of Σ_d must be more than an order of magnitude larger than that of the MMN. If this necessary condition is marginally satisfied then, in the presence of gas, instead of producing a single larger embryo, many embryos with much smaller isolation mass would emerge. Under these conditions some rocky Earth analogs or super-Earths could have formed in situ through collisional mergers after the gas disperses (as it can no longer damp embryos' mutually excited eccentricities). But giant impacts could still cause a substantial loss of volatiles in post-gas evolution, so this method may not lead to the emergence of the known low-density super-Earths.

In order for the magnitude of M_iso to reach a few M_⊕ at 0.1 AU, the local Σ_d must be increased by two orders of magnitude over the MMN. In disks with a solar composition, the corresponding gas surface density (Σ) would be sufficiently large that gravitational instability would induce a rapid redistribution of mass and angular momentum (Durisen et al. 2007; Gammie 2001; Boley et al. 2006). Protoplanetary embryo formation in such disk regions may be challenged because the disk temperature would be above the condensation temperature of most refractory grains and the Σ would decline substantially before the planetesimals can grow significantly. In gravitationally stable disks with a Σ distribution comparable to that of the MMN, the dust-to-gas ratio would need to be enhanced by more than an order of magnitude relative to its solar value (Ida & Lin 2005). This enhancement can be induced by planetesimal drift (e.g., Kretke et al. 2009; Hansen & Murray 2012) or by embryo migration.

Planetary embryos migrate due to their interaction with gas in their natal disks. Their tidal perturbation excites density waves in their Lindblad and corotation (or horseshoe) resonance regions. These waves carry angular momentum flux, which is determined by the surface density (Σ) and temperature (T) of the disk gas at the excitation location (Goldreich & Tremaine 1980). These waves dissipate and deposit angular momentum to the disk as they propagate away from the resonances (Papaloizou & Lin 1984). The disk responds with an evolving Σ and T distribution, which introduces a feedback effect (Takeuchi et al. 1996). In order to calculate the full nonlinear disk response, it is necessary to perform numerical hydrodynamic simulations of planets embedded in two- or three-dimensional disks. And mutual interactions between multiple planets embedded in a disk can lead to a wide varieties of diverse behaviors including outward migration and migration stalling (e.g., Masset & Snellgrove 2001; Podlewska-Gaca et al. 2012). Nevertheless, if we restrict our investigation to planets with modest masses, the disk structure is not strongly perturbed and their net linear Lindblad and corotation torque can be calculated by summing up the contribution of each distinctive resonance. In disks with a power-law Σ and T profiles it is possible to describe the net torque on the planets as a linear combination of the radial Σ and T gradient (Tanaka et al. 2002; Paardekooper et al. 2010).

In this paper, we are primarily interested in super-Earths' type I migration. The mass of these planets is generally too small to significantly perturb the structure of the disk. In such limit, super-Earths' migration rate can be deduced from the structure of unperturbed disk. Here, we apply various disk models to some recently obtained prescriptions for type I migration (Paardekooper et al. 2011, hereafter PBK11) and determine the direction and rate of migration for planets with a range of masses.

The simplest empirical prescription for the disk structure is the MMN model, which is based the assumption that all planets in the solar system attained their mass in situ and they accreted all the heavy elements in their neighborhood (Weidenschilling 1977b; Hayashi 1981). Through the augmentation of volatile elements for a solar elemental distribution, we can infer a Σ∝r^−3/2 distribution. With an additional assumption that protostellar disks are heated to their local equilibrium temperature, the disk temperature T∝r^−1/2. The MMN is a reasonable fiducial model, albeit the in situ formation assumption is inconsistent with the current concept of planetary mobility.

Alternatively, we can construct a set of dynamical models for the protostellar disks in which the Σ and T distribution are determined with physical principles. Due to uncertainties in the mechanism of angular momentum in disk it is common to use the α-viscosity model (Shakura & Sunyaev 1973) in which one assumes that the disk evolves due to a viscosity proportional to the sound speed (c_s) and disk scale height (h),

$$\begin{equation} \nu = \alpha c_s h. \end{equation} \tag{ 1 }$$

For the disk models we adopt here we first present results for disks with the common assumption that α is constant. We then investigate disks in which we assume the magnetorotational instability (MRI) excites turbulence and leads to angular momentum transfer in regions with sufficiently high ionization fraction (Hawley et al. 1995). In most planet-forming parts of the disk (ranging from a fraction to a few 10 or AUs), MRI-active regions are confined to the surface layers, which are exposed to ionizing stellar photons and cosmic rays (Gammie 1996). In order to include this effect, we generalize the quasi-steady-state accretion disk model to allow for a variable α parameter.

There are a number of groups who have presented sophisticated structure models of MRI-active disks (Sano et al. 2000; Terquem 2008; Kretke & Lin 2010). Due to the uncertainties on both MHD turbulence and dust grains' size distribution in the disk surface layer, we choose to use a model with fewer assumptions to simplify this issue.

A few studies have already begun investigating planet migration in MRI-active disks including the effects of a dead zone. Kuchner & Lecar (2002) and Masset et al. (2006b) investigated planet stalling at the inner edge of the dead zone. Matsumura et al. (2007, 2009) studied the growth and migration of Jovian planets in disks with dead zones. All of these studies indicate that planet migration is sensitive to the structure of the disk. In this paper, we present the evolution in simple disk models in order to discuss the generic disk profiles that can create "planet traps" and how these are likely to manifest themselves in physical disks.

As we have indicated above, the rate and direction of type I migration is determined by the competing contribution from differential Lindblad and corotation (or horseshoe) torque. In our computation of the condition for corotation torque saturation, we are mindful that the horseshoe region has a toroidal shape. Its height normal to the disk plane is comparable to its width in radial direction. For planets with mass less than a few M_⊕, the entire horseshoe region may be embedded in the dead zone. It is not clear whether there is any significant viscous flow across the horseshoe zone to prevent saturation. In these low-viscosity regions, embedded planets' long-term interaction may also modify the Σ distribution, trigger secondary instabilities, and halt their migration (Balmforth & Korycansky 2001; Li et al. 2005, 2009; Yu et al. 2010).

This paper is organized as follows. In Section 2, we briefly recapitulate relevant migration prescriptions. In Section 3, we construct three sets of disk models, which take into account surface irradiation and the presence of the "dead" midplane between MRI-active surface layers. In Section 4, we present the results on how planets will migrate in these models of static disks. In Section 5, we discuss potential complication introduced by layer structure in the dead zones and feedback. In Section 6, we discuss how this work connects to the ongoing effort of population synthesis modeling. In Section 7, we summarize our results and discuss their implications.

The type I migration rate is sensitive to the local Σ and T as well as to the radial gradients of these quantities. Following PBK11, we adopt the following general expression for the rate of type I migration,

$$\begin{equation} \frac{dr}{dt} = f(p,q,p_\nu ,p_\chi) \frac{M_p}{M_*}\frac{\Sigma r^2}{M_*} \left(\frac{r\Omega _{\rm K}}{c_s}\right)^2 r\Omega _{\rm K}, \end{equation} \tag{ 2 }$$

where M_p and M_* are the planet and stellar mass, respectively, Σ, T, c_s, and Ω_K all refer to the those quantities at the location of the planet. The magnitude and sign of the dimensionless efficiency factor f is a function of the dimensionless logarithmic surface density and temperature gradients, p(r) ≡ dlog Σ/dlog r and q(r) ≡ dlog T/dlog r, respectively. Note that the PBK11 α and β parameters correspond to −p and −q, respectively, in this paper.

In weakly perturbed disks, p and q are intrinsic properties of the disk alone. The efficiency factor f is also a function of the viscous and thermal saturation parameters

$$\begin{equation} p_\nu \equiv \frac{2}{3} \left({\rm Re\,} x_s^3\right)^{1/2} \end{equation} \tag{ 3 }$$

and

$$\begin{equation} p_\chi \equiv \left(\Omega _{\rm K}r^2 x_s^3 /2 \pi \chi\right)^{1/2}=3 p_\nu / 2\rm P_r^{1/2} \end{equation} \tag{ 4 }$$

where Re ≡ Ω_Kr²/2πν is the Reynolds number, x_s ≃ (M_pr/M_*h)^1/2 is the dimensionless width of the horseshoe orbit, h = c_sΩ⁻¹_K is the density scale height in the direction normal to the plane (i.e., the thickness) of the disk, and P_r = νχ⁻¹ is the Prandtl number. The magnitude of both p_ν and p_χ depends on the mass ratio M_p/M_* as well as disk parameters, such that

$$\begin{equation} p_\nu \propto \alpha ^{-1/2}(M_p/M_*)^{3/4}(h/r)^{-7/4} \end{equation} \tag{ 5 }$$

in a disk in which the viscosity is parameterized according to Equation (1).

Based on a three-dimensional linear analysis of wave excitation in a locally isothermal disk, Tanaka et al. (2002, hereafter TTW02) determined the migration rate to be

$$\begin{equation} f_{{\rm TTW}}(p,q) = 1.08(p+0.80q-2.52). \end{equation} \tag{ 6 }$$

In typical disks, differential Lindblad resonances generally lead to inward migration. Under some conditions, contribution from corotation resonances and horseshoe torques can lead to outward migration. Paardekooper et al. (2010, hereafter PBCK10) investigated disks in which there are fully unsaturated (i.e., with a maximum efficiency) horseshoe torques. From these upper torque limits, we can identify disk regions where planets have a potential to undergo outward migration. In a similar (to TTW02), locally isothermal limit, PBCK10 found that the combined contribution from linear Lindblad torques, linear-entropy-related corotation torque, and nonlinear-vortensity-related horseshoe drag (with a softening factor b/h = 0.4) leads to

$$\begin{equation} f_{\rm iso}(p,q) = -1.7-1.8q + 2p. \end{equation} \tag{ 7 }$$

PBCK10 also considered the effects of disk thermal structure. They found that in an inefficiently radiating disk with an adiabatic equation of state, a planet migrating under the fully unsaturated (i.e., with full strength) nonlinear horseshoe torques (with the adiabatic index γ = 1.4) yields

$$\begin{equation} f_{\rm ad}(p,q) = -1.22-5.64q+4.66p. \end{equation} \tag{ 8 }$$

With these migration laws, it becomes simple to identify the necessary conditions for stalling migration. The necessary criteria for outward migration as a function of p and q are listed in Table 1. For comparison and illustration purposes, we also included in this table the criteria relevant to the radial drift of particles due to aerodynamic drag (Weidenschilling 1977a).

The results in Table 1 indicate that with an adiabatic equation of state, outward migration in disks with negative surface density gradients (p < 0) is only possible if there is a strong negative temperature gradient (q < −0.22). Furthermore, in order for the torques to reach these fully unsaturated values, angular momentum must be extracted from the horseshoe region at the optimal rates, i.e., fast enough to avoid torque saturation, but slow enough to maintain a gradient across the horseshoe region. These requirements determine a set of sufficient conditions for stalling type I migration.

Here, we highlight migration rate's sensitive dependence on the disk structure. In the introduction, we have already indicated that the Σ and T distributions are determined by the not-fully-understood physics associated with the MRI turbulence and surface irradiation. In order to identify some robust features and assess range of uncertainties, we adopt three different disk models including (1) a relatively simple self-consistent constant α model (Shakura & Sunyaev 1973; Lin 1981; see Section 3.2) and (2) a model for the irradiated outer regions of the disk (see Section 3.3). We also consider (3) a disk model which takes into account the influence of a dead zone (see Section 3.4).

Before describing these three models, we orient the reader by describing some standard disk models and features common to most disk models. In Figure 1, we show where various disk models lie on the p–q plane. The Σ profile of the MMN corresponds to p = −1.5, and on this plot we place a star to indicate a geometrically flat (q = −1) MMN disk. It is interesting to note that in an MMN profile, the vortensity gradient is zero and the horseshoe toques do not contribute (see further discussion in Section 5).

Zoom In Zoom Out Reset image size

The observed spectral energy distributions of young stars indicate that at least the outer regions of most protostellar disks have a flaring structure. This feature means that h/r increases with r, an indication that these passive disks have slowly declining temperature gradients (q > −1). Chiang & Goldreich (1997, hereafter referred to as CG97) construct a model to match the observed features. In their model, there are three regions. Far from the star the disk is optically thin to both its own radiation and the stellar radiation, while in the inner region is optically thick to both. In between these two boundaries, there is a marginally opaque region in which the disk is optically thick to the stellar photons, while optically thin to the longer wavelength disk photons (see Section 3.3). The square symbols in Figure 1 indicate the temperature gradients of these three regions with an assumed MMN Σ profile. In this figure, the opaque and optically thin limits are virtually indistinguishable.

3.1. Steady-state Disk Models with Constant α

Instead of assuming a fixed Σ profile (such as that inferred empirically with the MMN scenario), it is also possible to deduce it from accretion disk models. In an viscous axisymmetric Keplerian disk, the mass accretion rate through each annulus is given by

$$\begin{equation} \dot{M} = 3\pi \nu \Sigma = 3\pi \alpha c_s^2 \Omega _{\rm K}^{-1}\Sigma. \end{equation} \tag{ 9 }$$

We assume that short-period super-Earths were formed within a few AUs prior to their type I migration. In these inner disk regions, gas attains a quasi-steady state; i.e., its mass accretion rate ($\dot{M}$) is independent of disk radius (Ruden & Lin 1986) albeit $\dot{M}$ does evolve over time. In such a state, Σ is inversely proportional to the effective viscosity ν.

It is generally assume that turbulence provides an effective angular momentum transfer mechanism and the magnitude of viscosity can be approximated with an ad hoc α prescription in which ν = αc_sH where c_s ≡ (R_gasT/μ)^1/2 (R_gas and μ are the gas constant and the mean molecular weight), H = c_s/Ω_K, and Ω_K ≡ (GM_*/r³)^1/2 are the sound speed, density scale height, and angular frequency of the disk gas, respectively (Shakura & Sunyaev 1973). If the efficiency factor α is constant with radius, this steady state would lead to a simple relationship between the surface density gradient and temperature gradient, such that p = −q − 3/2. The shaded rectangle in Figure 1 indicates the parameter space occupied by steady-state, constant α models.

In order to deduce a Σ distribution, it is necessary to independently or concurrently determine a T distribution. One possible approach is to assume an equilibrium T∝r^−1/2 distribution as in the MMN model and derive a Σ∝r⁻¹ distribution (Hartmann et al. 1998). Although this equilibrium T distribution is adequate for the optically thin outer regions of the disk, it is not appropriate for the inner disk regions, which is mostly heated by viscous dissipation. In standard accretion disk models, both T and Σ distribution can be obtained self-consistently under the assumption of thermal equilibrium, i.e., viscous dissipation is balanced by local radiative flux (Shakura & Sunyaev 1973; Lin 1981).

It is interesting to note that, although these standard disk models all lie in the lower region of this parameter space, the entire region indicated in this diagram satisfies the stability criteria: d/dr[(r²Ω)²] > 0 (Rayleigh's criterion) and dρ/dr > dρ/dz. Thus, the replacement of a constant α by a more general prescription does not generally introduce dynamically unrealizable models.

3.2. Description of Viscously Heated Disk Model

In order to determine temperature of a viscously heated disk, we balance the radiative losses at the surface (of temperature T_eff) with viscous heating and heating due to the background radiation, meaning

$$\begin{equation} 2\sigma _{\rm SB} T_{\rm eff}^4 = \frac{9}{4} \Sigma \nu \Omega _{\rm K}^2+2\sigma _{\rm SB}T_b^4, \end{equation} \tag{ 10 }$$

where σ_SB is the Stefan–Boltzman constant and T_b is the background temperature. The effective temperature can be related to the midplane temperature via the effective optical depth T⁴_effτ_eff = T⁴, where

$$\begin{equation} \tau _{\rm eff} \equiv \frac{3}{8}\tau + \frac{\sqrt{3}}{4} + \frac{1}{4\tau } \end{equation} \tag{ 11 }$$

and τ = κΣ/2 is the true optical depth (Hubeny 1990). For all of the disk models we parameterize the opacity as

$$\begin{equation} \kappa (T) = \kappa _0 T^\beta. \end{equation} \tag{ 12 }$$

3.3. Description of GL07 Model of Irradiated Disks

In the determination of the Σ and T distribution for the intermediate and outer disk regions, contribution from stellar irradiation needs to be taken into account. In an attempt to introduce a more comprehensive prescription, Garaud & Lin (2007) constructed an analytic disk model (hereafter referred to as GL07) based on the passive-disk model of CG97. However, instead of assuming an MMN surface density profile, GL07 calculated the structure of a disk evolving with a constant α-viscosity. Additionally, this model includes both stellar irradiation and viscous heating. This prescription results in five distinct regions. The outer three regions, the optically thin [v], marginally opaque [iv], and weakly optically thick [iii] regions are similar to those in CG97. In these regions the T distribution is independent of the Σ profile, although a steady-state solution can be preserved with a variable α prescription.

In the GL07 model, there are two regions close to the star that do not occur in the CG97 passive-disk model, a strongly opaque region [ii] and viscously heated region [i]. In region [ii], the optical depth to the disk's own radiation is high, causing the disk to be non-isothermal both in the disk plane and in the direction normal to it. In region [i], viscous dissipation dominates the thermal profile. In these regions q is a function of p. In the inner disk regions where accretion flow is in a quasi-steady state, both quantities (p and q) are well defined.

These five regions ([v] and [iii] are nearly indistinguishable) are indicated by triangles in Figure 1. The radial transitions between these regions evolve as a function of the mass accretion rate (see Figure 2). It is interesting to note that the optically thick region of GL07 lies remarkably close to the separatrix between outward and inward migration in Figure 1.

Zoom In Zoom Out Reset image size

3.4. Description of an MRI-active Disk Model

3.5. Inner Disk Truncation

A generic feature of protostellar disks is the presence of a cavity at their center, which is cleared by the magnetosphere of their central stars. The presence of this magnetospheric cavity quenches both Lindblad and horseshoe torques and therefore provides a natural post for hot Jupiters' type II migration to stall (Lin et al. 1996). Here, we consider the possibility that the termination of super-Earths' type I migration may have also been strongly influenced by this structure.

As an additional feature in our disk models, we consider the possibility of an inner boundary due to the disks' interaction with the stellar magnetosphere of their central stars. The location of this magnetospheric truncation radius is assumed to occur at

$$\begin{equation} r_{\rm mag}\simeq 0.5 \mu _*^{4/7} (2 {\it G M}_*)^{1/7} \dot{M}^{-2/7}, \end{equation} \tag{ 13 }$$

where μ_* is the stellar magnetic dipole moment (Königl 1991).

Many of the approximations integral to the disk model break down in this location. The actual Σ and T profiles in this disk region need a detailed analysis of the process of mass loading onto the stellar magnetic field as well as magnetic diffusivity and dissipation near the disk inner boundary. These studies are beyond the scope of this paper. For the purpose of studying the destiny of super-Earths, we impose a sharp positive surface density gradient near r_mag, where the surface density is reduced by the factor f_mag. In this paper we use

$$\begin{equation} f_{\rm mag} \equiv \exp \left(-\frac{1}{2}\left(\frac{r-r_{\rm mag}}{0.1 r_{\rm mag}}\right)^2\right), \end{equation} \tag{ 14 }$$

though we note that this is an arbitrary choice. And as noted by Tsang (2011), the interaction at the inner edge of the disk are likely more complicated than the simple picture presented here. As a planet approaches the inner edge of the disk, it is possible that reflections at the inner disk edge may cause the planets to have their migration halted at a few times the inner disk radius. But while the details are different, the unsaturated horseshoe torque near r_mag generally induce a rapid outward type I migration, creating a planet trap.

4.1. The Migration Rate of a Planet in a Steady-state Disk

In order to discuss the dominant and relevant physical effects, we illustrate the migration of a planet in a static (i.e., a steady-state) disk. In Figure 3, we show the relative migration rates (as compared to the TTW02 model) for a disk that is heated entirely by viscous dissipation (see Section 3.2). Such a structure is possible if the stellar irradiation is blocked by self-shadowing effects (Dullemond et al. 2002). We adopt the following model parameters: $\dot{M} = 10^{-8} \,M_\odot \,{\rm yr}^{-1}$, α = 10⁻³, T_b = 20K, and r_mag = 0.15 AU.

Zoom In Zoom Out Reset image size

From this standard viscous-disk model, we find distinct regions where the disk structure can allow outward migration. Other than the inner boundary of the disk, these locations include the regions far from the sublimation temperatures, i.e., zone 1 where the opacity is dominated by silicates (β = 1/2 in Equation (12)) and zone 2 where the opacity is dominated by ice grains (β = 2). In the regions interior to zone 1 (where silicate grains partially sublimate), between zones 1 and 2 (where ice grains partially sublimate), and outside zone 2 (where the disk become optically thin in the direction normal to the disk), fully unsaturated horseshoe torques cannot offset the differential Lindblad torque to induce planets to migrate inward.

4.2. Potential Location of Planet Traps

4.3. Migration During Planetary Growth

4.4. Sensitive Dependence on the Disk Models

In Figure 4, we show the landscape for a GL07 disk with $\dot{M}=10^{-8} \,M_\odot \,{\rm yr}^{-1}$, α = 10⁻³, β = 1. In the innermost region of the disk viscous heating dominates, and it is analogous to zone 1 in the fully viscous disk model (see discussions in the previous section). Consequently, planets with masses in the appropriate range undergo outward migration. Outside this region, stellar irradiation provides most of the heating for the disk and modifies the temperature distribution. Consequently, the horseshoe torque associated with the corotation resonance is weakened. There is a distinct strongly opaque region where the horseshoe torques reduce the impact of differential Lindblad resonance and significantly slow down but does not reverse planets' inward migration, even for those with optimum masses. The main reason for this reduction in the inward migration rate is due to the relatively steep decline of the midplane temperature with the radius. With the inner region being heated by viscous dissipation, the temperature gradient i.e., q(r) in GL07 is more negative than that for a totally passive, flaring disk in CG97.

Zoom In Zoom Out Reset image size

For these choices of parameters, the most massive embryos that can be trapped is ∼3–4 M_☉ in an unperturbed disk at a fraction of an AU. This mass range is well below that (M_c) required for the acquisition of massive gaseous envelope for gas giants through the onset of efficient gas accretion. The location is also too close to host stars to be compatible with that of most known gas giants. The upper mass limit of retainable embryos increases with the poorly determined and constrained magnitude of α. But, for a given $\dot{M}$, the surface density of the gas and presumably the heavy elements decreases with α. A larger value of $\dot{M}$ or a smaller value of α may also lead to the growth and retention of more massive cores. The torque applied by the cores can also lead to modifications in the Σ distribution and this feedback effect can also reduce the type I migration rate (see Section 5).

4.5. Disks with Dead Zones and Active Layers

In Figure 5, we show the relative migration rates for a planet in a disk with an MRI-active surface layer. For simplicity, we first consider a model in which the surface density Σ_A of the MRI-active layer is a constant (Σ_A = 50 g cm⁻²), similar to the seminal Gammie (1996) paper, and potentially there is a dead zone below it with surface density Σ_D. We take the viscosity in the active layers to be large with an α_A = 10⁻², while the viscosity in the dead zone is small α_D = 10⁻⁴. With this small but non-vanishing value for α_D, we approximate the propagation of turbulent disturbances near the disk surface into the midplane regions of the disk. The actual values of Σ_A and α_D are still poorly constrained (Fleming & Stone 2003; Turner et al. 2007).

Zoom In Zoom Out Reset image size

In regions where Σ_A is smaller than the total surface density, accretion flow in the active surface layer of the disk leads to a smaller effective magnitude for the vertically averaged α ≡ (α_AΣ_A + α_DΣ_D)/Σ. The main consequence of this modification is the expansion in both the lower and upper mass range for embryos that can be retained. The upper mass range exceeds the magnitude of M_c. The lower mass range extends well below the isolation mass for embryos even at relatively small radius. However, if we take into account the effect of vertical structure on the saturation parameter, the lower mass range would be limited to a fraction of the Earth mass. Although these tendencies maybe favorable to the formation of gas giants, the location of the planet-trapping radius is only a fraction of an AU. Larger values $\dot{M}$ are needed for the planet-trapping radius to become comparable with the typical orbital semimajor axis of known gas giants.

A robust feature in the above results is a tendency for super-Earths to migrate outward in the inner regions and inward in the outer regions of the disk. This property implies the presence of a planet-trapping radius where these planets may converge and accumulate. We also showed that the mass range of the potentially trapped planets and the location where their type I migration is stalled are model dependent and subject to uncertainties. In this section, we list several additional subtle effects that may also affect these quantities.

5.1. Disks with Variable Active Layers

In Figure 6, we show the relative migration rates for a planet in such a disk with a Σ_A that varies with the local dust-to-gas ratio. We use an approximate fit to the results of Ilgner & Nelson (2006, model 4) in which

$$\begin{equation} \Sigma _A \simeq \left(\frac{0.0025}{f_d}-0.018\right) \left(\frac{r}{\rm AU}\right)^{2.87}, \end{equation} \tag{ 15 }$$

where f_d is the dust-to-gas ratio for micron-sized particles. This prescription introduces a feature at the snow line where the dust-to-gas ratio decreases as a planet moves inward due to the sublimation of ice grains (Kretke & Lin 2007, 2010). In this confined transition region, embryos can be trapped at several AUs from their host stars. However, we note that for these parameters the mass of cores that can be trapped is relatively small, below M_c ∼ 10 M_⊕. Whether or not these cores are sufficiently massive to accrete atmospheres and become gas giants depends upon properties such as the planetesimal mass accretion rate and atmospheric opacity (Hubickyj et al. 2005; Rafikov 2006).

Zoom In Zoom Out Reset image size

5.2. Revisiting the Horseshoe Torques in the Dead Zone

In addition to the reduction in α and the corresponding modification in the Σ and T distribution, a layered disk structure may also affect the saturation of embedded planets' corotation torque. The horseshoe region where the corotation torque is applied has a dimensionless width in the disk plane x_s ≃ (M_pr/M_*h)^1/2 (Paardekooper & Papaloizou 2009). In the direction normal to the disk plane, the horseshoe region extend to a similar dimensionless height z_s ∼ x_s.

For a 10 M_⊕ planet (M_p/M_☉ ∼ 3 × 10⁻⁵) embedded in the midplane of a disk with a typical h/r ∼ 0.05, z_s ∼ 0.02. For these parameters, there is a substantial dead zone that extends over one disk scale height (so that the dimensionless height of the dead zone z_d > h/r) and the entire horseshoe region is engulfed within the dead zone. In this situation, the horseshoe region is entirely isolated from the active layer. Ideally, the continuous viscous diffusion of disk gas into the horseshoe region of relatively low mass planets may be quenched if their z_s < z_d. It is unclear whether the Reynold's stresses that provide a source of angular momentum transport in the dead zone (and therefore generate an effective α_D) may also redistribute vortensity and entropy between the horseshoe region and the surrounding disk with a similar efficiency. Without such redistribution, the corotation torque of low-mass planets with coplanar orbits would be saturated after a few libration periods.

In principle, the limited saturation for planets with h/r > z_s > z_d can be estimated from three-dimensional hydrodynamic simulations. But, such simulations are time consuming, especially for such low masses. Therefore, we utilize our approximation to highlight how including this effect may impact the migration rates. In Figure 7 we demonstrate the results if we assume that the effective viscosity is used in determining the Reynolds number (Re) in Equation (3) is negligible if z_s < z_d, meaning that, for the sole purpose of calculating Re, we set α = 10⁻⁸. Under this assumption, the horseshoe torque saturate for nearly the entire dead zone, resulting in inward migration for almost all masses and positions.

Zoom In Zoom Out Reset image size

As a comparison, in Figure 8 we assume that if z_s < z_d then the relevant effective viscosity is the dead zone viscosity α_D. In this case the horseshoe torques saturate at a much smaller mass so we still see stalling locations for small planetesimals, but we lose the ability to stall cores greater than around 3 M_⊕.

Zoom In Zoom Out Reset image size

In these estimates, we assume the cores have coplanar orbits. In N-body simulation of planetesimal dynamics and embryo formation, their mutual gravitational perturbation and collisions can lead to dispersions in the velocity (σ_z) normal to the disk plane. During the initial growth, σ_z of the most massive embryos is damped by dynamical friction. But during the advance Oligarchic growth stage, perturbation between comparable mass embryos may induce excitation σ_z to be a fraction of their surface escape speed, despite gas drag (Kokubo & Ida 1998). If σ_z > Ω_Kz_d, the embryos would encounter MRI-active regions during at least a fraction of each Keplerian orbital period and the reduction of the diffusion through the horseshoe regions would no longer be an issue.

5.3. Feedback and Transition to Type II Migration

In the above sections, we assume that planets with masses up to a few M_⊕ do not significantly alter the average disk Σ and T distribution. In our evaluation on the direction and rate of such planets' type I migration and the extent of saturation (i.e., the magnitude and sign of f), we calculate the Σ and T distribution neglecting any feedback effects. However, we note that previous research has demonstrated that, even with a modest mass (such as ∼5 M_⊕), planets may alter the disk structure, especially in the horseshoe region, and slow down the rate of migration (Masset et al. 2006a).

Additionally, the type I migration rate is modified once the interaction between the disk and a planet become nonlinear. This limit is reached when the planet's Hill sphere becomes comparable to the disk scale height, i.e., once

$$\begin{equation} \Theta _1 \equiv \frac{M_p}{M_\ast }\left(\frac{r}{h}\right)^3 > 1. \end{equation} \tag{ 16 }$$

Additionally, the planet begins to open a gap (Lin & Papaloizou 1993), when

$$\begin{equation} \Theta _2 \equiv \left(40 \alpha \left(\frac{h}{r}\right)^5\right)^{1/2} \frac{M_\ast }{M_p} < 1. \end{equation} \tag{ 17 }$$

Gap formation interrupts accretion flow through the disk and consequently triggers a transition from type I to type II migration (Ward 1997). This limit is more relevant to gas giants rather than super-Earths.

Though super-Earths do not generally have sufficient mass to open gaps and, their feedback can significantly affect type I migration process. In the low-viscosity (with α ≪ 1) limit, gas in the horseshoe region encounters weak shock when it flows close to a-few-Earth-mass planets' Roche lobe. Wave dissipation across the shock leads to angular momentum deposition and evolution in the Σ distribution in the horseshoe region (Li et al. 2009). On the horseshoe libration timescale, this effect can significantly reduce the torque magnitude and planets' type I migration rate (Yu et al. 2010).

Shock-induced dissipation along the streamlines also modifies the vortensity (the vorticity divided by the surface density) along each streamline. The development of extrema in the vortensity distribution across the streamlines promotes the growth of secondary or Rossby vortex instabilities (RVI), which can be suppressed by adequate viscosity. In the limit of negligible viscosity, the growth of RVI can lead to the emergence of large vortices (Li et al. 2005). The torque induced by the interaction between the planets and the vortices is chaotic. It can affect both the direction and rate of type I migration. Nevertheless, the onset of RVI significantly reduces the effective magnitude of type I migration rate (Yu et al. 2010). These feedback effects can extend the upper mass limit for the trapped planets.

In order to assess how these various disk models may impact the emergence and appearance of planetary system one may wish to include these prescriptions in population synthesis models. We note that without actually performing a detailed population synthesis model we cannot compare how precisely the models will be affected, but we can compare these models to existing formulations. In Ida & Lin (2008a) and Ida & Lin (2010) the authors assume that the gas follows the radial surface density gradient of the minimum mass solar nebula and the temperature follows a power law of T∝r^−1/2. As indicated in Figure 1, this disk structure would lead to inward migration under both the isothermal and adiabatic regimes. Therefore, given their disk model, the inclusion of the additional thermal and saturation effects discussed in this paper will not significantly change their result that the predicted type I migration rate is too fast to form the observed planet population. In Ida & Lin (2008b), the authors adjust the surface density to include a simple prescription for the dead zone. However, the same thermal profile is still used. In this case fully considering the possibly of torque saturation causes the snow-line "planet trap" to be only effective for a limited range of masses. The largest change for these models is that if they used a self-consistent temperature model they could create inner, viscously heated region. In these studies, the authors recognize that in this region they must artificially stop the planets migration, and the viscously heated inner region may be a natural way to achieve this effect.

In Mordasini et al. (2011) the authors used the viscous heating prescription of Lyra et al. (2010) and Sándor et al. (2011), which is very similar to the model described in Section 3.2 and Figure 3. This leads to various "concentration points" for planet migration. However, in this paper we argue that a viscously heated model may not be valid for large regions of the disk. If indeed the disk is passively heated from the star then the disks will not have planet traps far from the parent star unless there is dead zone or other change in viscosity. If this is the case, then the results may be more similar to those shown in Mordasini et al. (2012) in which the authors use a viscously heated model but do not account for the thermal effects. This means that they see inward migration throughout the disk without stalling zones.

In this paper, we discuss the migration of low-mass planets (super-Earths) in gaseous disks via type I migration including the impact of horseshoe torques as presented by Paardekooper et al. (2011). We find that planets with a few M_⊕ have a tendency to undergo outward migration through regions that are predominantly heated by viscous dissipation and inward through regions that are mainly heated by stellar irradiation. Since typical protostellar disks are mostly heated by viscous dissipation at small radii and by stellar irradiation at large radii, these planets tend to converge and stall at a fraction of 1 AU from their host stars. This result is robust and general for self-consistent protostellar disk structures in accordance with the conventional α models. Thus, we anticipate a high retention efficiency for super-Earths. This extrapolation is consistent with the observed ubiquity of super-Earths.

However, the location of the stalling radius is somewhat model dependent. Our numerical results indicate that its value is smaller in models that take into consideration stellar irradiation than in models that are primarily heated by viscous dissipation (e.g., Lyra et al. 2010). In the composite (viscous dissipation plus surface irradiation) models, it is difficult to sustain outward migration to large radii without any self-shadowed effects. While self-shadowing by a puffed up inner rim is possible (Ke et al. 2012), it is unclear whether this shadowing can extend beyond a few AUs.

In all of the models the planets migrate inward or outward until they reach a barrier such as the inner edge of the dead zone r_id, the magnetospheric truncation radius, or other stalling radii. However, even if a planet stalls at a given radius within the disk it may not determine its final location because, as a disk evolves, the locations of these barriers tend to change. For example, r_mag increases and r_id decreases with diminishing $\dot{M}$. This sensitive dependence on the disk structure may lead to large dispersion in the asymptotic period of super-Earths. In contrast to a 3 day pile-up in hot Jupiters' period distribution, super-Earth candidates, which are identified by the Kepler mission, show a nearly uniform logarithmic period distribution without any preferred orbital configuration. The observed spread-out period distribution of super-Earths does not contradict the expected large dispersion in both the disk structure and planets' stalling condition.

For planets that do reach the inner edge of the disk, the stalling locations for low-mass planets differ from those of more massive (Jovian mass) gap-opening planets. These high-mass planets terminate their type II migrating once they have reached a 2:1 resonance interior to the inner edge of the truncated disk or in the limit that the disk mass becomes smaller than the planet mass (Lin et al. 1996). Therefore, we expect that short-period high-mass planets (hot Jupiters), the product of type II migration, should be found on average interior lower mass planets, whose migration is stalled at the magnetospheric truncation radius. As we continue to detect low-mass planets on short-period orbits, they will become an important tool to constrain models of planet formation and migration. The final stopping location of short-period planets will reveal information about the structure of the inner disk region.

Additionally, as type I migration is proportional to the planetary mass, planets on the upper range of the type I migration regime are the most vulnerable. For these higher-mass embryos the horseshoe torque begin to saturate, meaning that they migrate inward at the full rate. However, these embryos are also the planets which are able to modify the Σ profiles of the gas disk in their vicinities. These feedback effect could quench type I migration, induce vortex formation through RVI, or lead to type II migration. These effect may lead to the falloff in the observed size distribution (at around 3–4 Earth radii).

High mass accretion rate promotes core retention at large radii. This correlation is primarily due to the fact that the aspect ratio (h/r) generally increases with $\dot{M}$. In this limit, the horseshoe torque does not saturate until the embryos have acquired a relatively large planet mass. For sufficiently large $\dot{M} ({>} 10^{-7} \,M_\odot$ yr⁻¹), stalled embryos can attained sufficiently large mass to initiate efficient gas accretion. A larger α can also prevent inward migration of these particularly vulnerable cores. However, the location of the transition from the viscously heated region to the strongly opaque region (i.e., the transition from inward to outward migration) depends only on $\dot{M}$ not on α. In general, the existence of a dead zone hinders the formation of cores at large radii because the low viscosity prevents the barrier for high-mass cores. However, the additional degree of freedom provided by the variable α parameter open the possibility of rapid outward migration for embryos with various masses. This effect may create planet traps, but only if the horseshoe torques can remain unsaturated.

A model in which planetary cores are retained at the "planet traps" begs one very important question. Once a planetary embryo acquires more than one Earth mass, it would migrate to one of these trapping locations on relatively short timescales, unless there is another mechanism to slow type I migration. Therefore, if, as expected, a system of multiple embryos emerge in a disk around a common host star, they would have a tendency to migrate to locations which are already occupied by others. This convergent dynamical evolution raises the question of how two (or more) embryos interact as they approach one another. This interaction needs to be explored in future studies, especially as it is clear that super-Earths appear to primarily be found in multi-planet systems.

The prescription generated in this paper can be applied to population synthesis models (Ida & Lin 2004, 2010). Through statistical comparisons between these models and the observed data, it is possible to highlight both robust physical processes that lead to common features and marginal effects that introduce diversity and dispersion.

We thank Dr. Clement Baruteau for useful conversation. This work is supported in part by NASA grants NNX07AI88G, NNX08AL41G, and NNX08AM84G as well as an NSF grant AST-0908807. K.K. appreciates funding from the NASA OSS program (PI: H. Levison).

While GL07 is presented for a disk model with a constant α parameter, it can be generalized to include radial dependence. The temperature models presented in GL07 are entirely local, and therefore independent of the radial variation of α. So for an arbitrary surface density profile one may simply use the prescriptions presented in GL07. However, in order to determine the surface density profile, assuming a steady-state disk,

$$\begin{equation} p=-q -\zeta - \frac{3}{2}. \end{equation} \tag{ A1 }$$

We define ζ ≡ dln α/ln r. In the optically thin, marginally opaque, and weakly opaque regions, the temperature profile is independent of the surface density profile so one can directly compute the surface density profile from this relationship. In the strongly opaque and viscously heated region, the temperature gradient depends on the surface density

$$\begin{equation} p_{\rm SO} = -\frac{15n+15-9\beta + \zeta (14n+14-2\beta)}{2(7n+15-\beta)}, \end{equation} \tag{ A2 }$$

where n is the polytropic index and the opacity varies according to κ(T) = κ_V(T/T_*)^β:

$$\begin{equation} q_V = \frac{9n-3\beta +15-4\zeta }{4n+8}-3 \end{equation} \tag{ A3 }$$

For the fiducial model in which n = 2 and β = 1,

$$\begin{equation} p_{\rm SO} = -\frac{9+10\zeta }{14}, \end{equation} \tag{ A4 }$$

and

$$\begin{equation} q_{\rm SO} =-\frac{6+2\zeta }{7}. \end{equation} \tag{ A5 }$$

In the viscously heated region

$$\begin{equation} p_{\rm V} = -\frac{3+6\zeta }{8} \end{equation} \tag{ A6 }$$

and

$$\begin{equation} q_{\rm V} = -\frac{9+2\zeta }{8}. \end{equation} \tag{ A7 }$$

We note that these are valid only under the same limit of GL07 that

$$\begin{equation} \frac{h}{r} \frac{d \ln (Z \rho _m)}{d \ln r} \ll \frac{z_s^2}{h^2}r\frac{d}{dr}\left(\frac{h}{r}\right). \end{equation} \tag{ A8 }$$