OUTWARD MIGRATION OF JUPITER AND SATURN IN EVOLVED GASEOUS DISKS

Consider a gaseous disk orbiting a central star of mass M_s. In the framework of one-dimensional (1D) modeling, we assume azimuthal symmetry around the star and use vertically averaged quantities as a function of the radial distance r. For the current purposes, we assume that the evolution of the disk is driven by viscous torques, $\mathcal {T}_{\nu }$, and wind dispersal, $\dot{M}_{\mathrm{w}}$ at the disk's surface. The torque exerted on a disk ring of radius r, by material orbiting inside the ring, is $\mathcal {T}_{\nu }=-2\pi r^{3}\nu \Sigma \partial \Omega /\partial r$ (Lynden-Bell & Pringle 1974), where ν is the kinematic viscosity of the gas, Σ the surface density, and Ω the angular velocity. If Ω is identified as the Keplerian velocity (i.e., if effects of gas and magnetic pressure gradients are neglected), then $\mathcal {T}_{\nu }=3\pi \nu \Sigma \mathcal {H}$, where $\mathcal {H}=r^{2}\Omega$ is the specific angular momentum of the gas. Along with viscous diffusion, the disk is dispersed by a wind, whose origin is gas photoevaporation from the disk surface produced by photons emitted by the central star. Hence, we write $\dot{M}_{\mathrm{w}}=2\pi \int \dot{\Sigma }_{\mathrm{pe}} r dr$, where $\dot{\Sigma }_{\mathrm{pe}}$ is the mass per unit surface area and unit time removed from the disk.

The continuity equation for the disk requires that

$$\begin{equation} \frac{\partial }{\partial t}\Sigma + \frac{1}{r}\frac{\partial }{\partial r}(r\Sigma u_{r}) = -\dot{\Sigma }_{\mathrm{pe}}, \end{equation} \tag{ 1 }$$

where u_r is radial velocity of the gas. On the left-hand side, one can recognize the mass per unit time flowing through a circumference of radius r, $\mathcal {F}=2\pi r \Sigma u_{r}$. By using the relation $\mathcal {F}=-\partial \mathcal {T}_{\nu }/\partial \mathcal {H}$ (see Lynden-Bell & Pringle 1974) and since we assume Keplerian rotation ($\partial \mathcal {H}/\partial r = r\Omega /2$), Equation (1) becomes

$$\begin{equation} \pi r \frac{\partial }{\partial t}(\Sigma +\Sigma _{\mathrm{pe}}) - \frac{\partial }{\partial r}\left(\frac{1}{r\Omega }\frac{\partial \mathcal {T}_{\nu }}{\partial r}\right) = 0. \end{equation} \tag{ 2 }$$

To seek for numerical solutions of Equation (2), it is convenient to use $\mathcal {H}$ as an independent variable and $\mathcal {S}=\mathcal {H}^{3}\Sigma$ as a dependent variable and then solve

$$\begin{equation} \frac{\partial }{\partial t}(\mathcal {S}+\mathcal {S}_{\mathrm{pe}}) - \frac{3}{4}({\it G\mbox{$M_{s}$}})^2 \frac{\partial ^{2}}{\partial \mathcal {H}^{2}}\left(\frac{\mathcal {\nu \mathcal {S}}}{\mathcal {H}^2}\right) = 0. \end{equation} \tag{ 3 }$$

In the above equation, G is the gravitational constant and $\dot{\mathcal {S}}_{\mathrm{pe}}=\mathcal {H}^{3}\dot{\Sigma }_{\mathrm{pe}}$. Note that the quantity $\mathcal {S}/\mathcal {H}^{2}$ is the angular momentum per unit surface area. In writing Equation (1), we neglected the effects of the star's growth, which would introduce a term on the right-hand side of the order of $\Sigma \dot{M}_{s}/\mbox{$M_{s}$}$ (see Ruden & Pollack 1991). Given the initial values of $\dot{M}_{s}/\mbox{$M_{s}$}$ considered here (see Section 2.3) and the decline of $\dot{M}_{s}$ with time, this term would affect Σ only over a timescale of the order of 10⁷ years, or longer.

Photoevaporation involves contributions from far-ultraviolet (FUV), extreme-ultraviolet (EUV), and X-ray radiation emitted by the star (see Dullemond et al. 2007; Clarke 2011 and references therein). FUV radiation may be especially important in removing gas at large distances from the star, reducing the gas supply to the inner parts of the disk. However, a self-consistent calculation of FUV photoevaporation rates requires solving for the detailed vertical structure of the disk (e.g., Gorti et al. 2009; Gorti & Hollenbach 2009). Photoevaporation by EUV photons is more tractable since they ionize hydrogen at the very upper layers of the disk. Here we follow a simple approach and adopt the formulation of the EUV photoevaporation rate proposed by Dullemond et al. (2007):

$$\begin{equation} \frac{\dot{\Sigma }_{\mathrm{pe}}}{\dot{\Sigma }^{\mathrm{g}}_{\mathrm{pe}}} =\left\lbrace \begin{array}{@{}ll}\exp {\left[\frac{1}{2}\left(1-\frac{r_{\mathrm{g}}}{r}\right)\right]} \left(\frac{r_{\mathrm{g}}}{r}\right)^{2} & \mathrm{for}\ r\le r_{\mathrm{g}},\\ \ms\ms \left(\frac{r_{\mathrm{g}}}{r}\right)^{5/2} & \mathrm{for}\ r> r_{\mathrm{g}}. \end{array} \right. \end{equation} \tag{ 4 }$$

The radius r_g ≈ 10(M_s/M_☉) AU is the gravitational radius, beyond which gas at the disk surface is unbound (see, e.g., Armitage 2011 and references therein). The photoevaporation rate at r_g is

$$\begin{equation} \dot{\Sigma }^{\mathrm{g}}_{\mathrm{pe}}= 1.16\times 10^{-11}\sqrt{f_{41}}\left(\frac{1\,\mbox{AU}}{r_{\mathrm{g}}}\right)^{3/2}\, \left(\frac{\mbox{$M_{\odot }$}}{\mbox{AU}^{2}\,\mathrm{yr}}\right), \end{equation} \tag{ 5 }$$

where f₄₁ is the rate of EUV ionizing photons emitted by the star in units of 10⁴¹ s⁻¹. The total mass-loss rate due to photoevaporation is found by integrating Equation (4) over the entire disk according to the definition given above, hence $\dot{M}_{\mathrm{w}}= (0.55977\sqrt{e}+4\pi )\,r^{2}_{\mathrm{g}}\, \dot{\Sigma }^{\mathrm{g}}_{\mathrm{pe}}$ or

$$\begin{equation} \dot{M}_{\mathrm{w}}= 1.56\times 10^{-10}\sqrt{f_{41}\left(\frac{r_{\mathrm{g}}}{1\,\mbox{AU}}\right)}\, \mbox{$M_{\odot }$}\,\mathrm{yr}^{-1}. \end{equation} \tag{ 6 }$$

The maximum of $\dot{\Sigma }_{\mathrm{pe}}$ occurs at r = r_g/4. Locally, gas is removed via photoevaporation and supplied by viscous diffusion, i.e., accretion through the disk, $\dot{M}=-\mathcal {F}$ (note that $\mathcal {F}$ is positive for an outward transfer of mass). Recalling the relations reported above, we can write $\partial \mathcal {T}_{\nu }/\partial \mathcal {H}= 3\pi \partial (\nu \Sigma \mathcal {H})/\partial \mathcal {H}$, and thus

$$\begin{equation} \dot{M}=3\pi \left[\nu \Sigma +2 r \frac{\partial }{\partial r}(\nu \Sigma)\right]. \end{equation} \tag{ 7 }$$

For νΣ nearly independent of r, i.e., in a stationary disk (Pringle 1981; see also Equation (1) with the right-hand side set to zero), $\dot{M}=3\pi \nu \Sigma$ is nearly constant throughout the disk. Therefore, if M_s = 1 M_☉, we expect gas depletion induced by photoevaporation to occur first around ∼3 AU.

In order to determine the thermal energy budget of the disk during its evolution, we assume that there is a balance among three terms: viscous heating, irradiation heating by the central star, and radiative cooling from the disk's surface. Viscous dissipation produces an energy flux equal to Q_ν = νΣ(r∂Ω/∂r)² (see, e.g., Mihalas & Weibel Mihalas 1999), which in the case of Keplerian rotation becomes

$$\begin{equation} Q_{\nu }=\frac{9}{4}\nu \Sigma \Omega ^{2}. \end{equation} \tag{ 8 }$$

Since Q_ν∝1/r³, for a disk with ∂(νΣ)/∂r ≈ 0, viscous dissipation becomes an ever less important source of heating as the distance from the star increases.

We follow the formulation of Hubeny (1990) for an irradiated disk and write the energy flux escaping from both sides of the disk surface as

$$\begin{equation} Q_{\mathrm{cool}}=2\sigma _{\mathrm{SB}}\,T^{4}\left(\frac{3}{8}\tau _{\mathrm{R}}+\frac{1}{2}+\frac{1}{4\tau _{\mathrm{P}}} \right)^{-1}, \end{equation} \tag{ 9 }$$

whereas the heating flux arising from stellar irradiation can be written as

$$\begin{equation} Q_{\mathrm{irr}}=2\sigma _{\mathrm{SB}}\,T^{4}_{\mathrm{irr}}\left(\frac{3}{8}\tau _{\mathrm{R}}+\frac{1}{2}+\frac{1}{4\tau _{\mathrm{P}}} \right)^{-1}. \end{equation} \tag{ 10 }$$

In the above equations, σ_SB is the Stefan–Boltzmann constant, T the mid-plane temperature, and T_irr the irradiation temperature. Note that, for an irradiated disk, the constant in parenthesis on the right-hand side of Equation (9) is generally slightly different from that of a non-irradiated disk (compare with Equation (14) of D'Angelo et al. 2003). As in Menou & Goodman (2004), we set

$$\begin{equation} T^{4}_{\mathrm{irr}}=(1-\epsilon)\,T^{4}_{s} \left(\frac{R_{s}}{r}\right)^{2}W_{\mathrm{G}}, \end{equation} \tag{ 11 }$$

where is a measure of the disk's albedo, for which we adopt the value 1/2, and T_s and R_s are the effective temperature and radius of the star, respectively. This interpretation of the irradiation temperature, however, neglects the contribution of luminosity released by stellar accretion (e.g., Hartmann et al. 2011). In an actively accreting disk, quantity T⁴_s should be replaced with T⁴_* = T⁴_s + T⁴_acc, where T⁴_acc quantifies the luminosity due to accretion $L_{\mathrm{acc}}=G\mbox{$M_{s}$}\dot{M}_{s}/(2R_{s})$ (Pringle 1981) and thus

$$\begin{equation} T^{4}_{\mathrm{acc}}=\frac{1}{8\pi } \left(\frac{G\mbox{$M_{s}$}\dot{M}_{s}}{\sigma _{\mathrm{SB}}R^{3}_{s}}\right), \end{equation} \tag{ 12 }$$

where the accretion rate $\dot{M}_{s}$, computed as $-\partial \mathcal {T}_{\nu }/\partial \mathcal {H}$ (see Section 2.1) at the disk's inner radius, varies with time.

The quantity W_G in Equation (11) is a geometrical factor that accounts for the illumination of disk portions close to (first term) and far from (second term) the star (see Chiang & Goldreich 1997):

$$\begin{equation} W_{\mathrm{G}}=0.4 \left(\frac{R_{s}}{r}\right)+ \frac{H}{r}\left(\frac{d\ln {H}}{d\ln {r}}-1\right). \vspace*{4pt} \end{equation} \tag{ 13 }$$

The adiabatic scale height of the disk, $H=\sqrt{\gamma \,k_{\mathrm{B}}T/(\mu m_{\mathrm{H}})}/\Omega$, is derived from the requirement of vertical hydrostatic equilibrium. The adiabatic index, γ, is 1.4, the mean molecular wight, μ, is 2.39, k_B is the Boltzmann constant, and m_H the hydrogen mass.

If the second term on the right-had side of Equation (13) is negative, the disk is self-shadowed and that term should be dropped. A self-consistent calculation of this term from 1D, vertically averaged models may lead to numerical instabilities (see, e.g., Hueso & Guillot 2005). In fact, meaningful determinations of this term involve solving for the vertical thermal structure of the disk. Therefore, the last term in parenthesis on the right-hand side of Equation (13) is written as η and approximated to 2/7 (see, e.g., D'Alessio et al. 1998; Menou & Goodman 2004; Hueso & Guillot 2005; Rafikov & De Colle 2006).

The optical depths τ_R = κ_RΣ/2 and τ_P = κ_PΣ/2 in Equations (9) and (10) are based, respectively, on Rosseland (κ_R) and Planck (κ_P) mean opacities. Both κ_R and κ_P depend on T and the mass density ρ = Σ/(2H). We adopt grain opacities from Pollack et al. (1994), at temperatures below the vaporization temperatures of silicates, and gas opacities from Ferguson et al. (2005) for solar abundances, when all grain species have evaporated.

The thermal energy budget is given by

$$\begin{equation} Q_{\nu }+Q_{\mathrm{irr}}-Q_{\mathrm{cool}}=0. \vspace*{4pt} \end{equation} \tag{ 14 }$$

Note that if Q_ν ≪ Q_irr, a situation that may occur in an evolved disk, Equation (14) results in a gas temperature T = T_irr, that is,

$$\begin{equation} T=T_{*}\sqrt{\frac{R_{s}}{r}} \left[(1-\epsilon)\,W_{\mathrm{G}}\right]^{1/4}. \vspace*{4pt} \end{equation} \tag{ 15 }$$

The factor W_G is typically a weakly dependent function of T. If W_G is a constant, then T∝r^−1/2. If W_G∝R_s/r (e.g., at radii r ∼ R_s), then T∝r^−3/4. If W_G∝H/r (as we assume for r ≫ R_s), then W^1/4_G∝T^1/8, the temperature is T∝r^−3/7 (see also Chambers 2009), and the disk's aspect ratio is

$$\begin{equation} \left(\frac{H}{r}\right)^{7}=\eta \left(1-\epsilon \right) \left(\frac{\gamma \,k_{\mathrm{B}}T_{*}}{\mu m_{\mathrm{H}}}\right)^{4}\! \left(\frac{R_{s}}{G\mbox{$M_{s}$}}\right)^{4}\! \left({\frac{r}{R_{s}}}\right)^{2}. \end{equation} \tag{ 16 }$$

The choice of the parameter η may have some impact on the disk's thermal budget, yet Equation (16) suggests that this impact is low.

Equation (3) is evolved in time using an implicit numerical scheme, which avoids the sometimes prohibitively short time steps required by an explicit approach, especially when the inner disk radius extends very close to the star (see Bath & Pringle 1981). We use either a second-order Crank–Nicolson method (e.g., Press et al. 1992) or a fourth/fifth-order Dormand–Prince method with an adaptive step-size control based on the global accuracy of the solution (Hairer et al. 1993). In the latter case, the evaluation of derivatives (in the Runge–Kutta sequence) is performed by means of a backward Euler (implicit) method. A zero-torque boundary condition, $\mathcal {S}=0$, is applied at the disk's inner edge. At the outer edge, the applied boundary condition is such that $\partial \dot{M}/\partial \mathcal {H}= \partial ^{2}{\mathcal {T}_{\nu }}/\partial \mathcal {H}^{2}$ is constant. Figure 1 (left) shows a comparison between numerical (lines) and analytic (circles) solutions of Equation (3) (see the figure caption for details).

Zoom In Zoom Out Reset image size

Equation (14) is solved for the mid-plane temperature, T, at each radius, using a root-finding algorithm based on the Brent's method (Brent 1973). Convergence of the root-finding process is achieved within a tolerance of 10⁻³ K. An iterative procedure is implemented for each determination of T, so that the applied value of H and that corresponding to the converged temperature do not differ by more than 1%. In Figure 1 (right), the evolution of temperature is shown for the disk considered in the left panel. In this test, we set W_G = 0.05, so that the temperature evolves toward that in Equation (15), indicated as filled circles. The temperature profiles show major opacity transitions at T ≈ 160, 420, 680, and 1400 K, caused by vaporization of, respectively, water ice, refractory organics, troilite, and silicate grains (see Pollack et al. 1994). Note that heating via viscous dissipation is basically confined within ∼10 AU (see Section 2.2).

The solar nebula extends from r_in = 0.01 AU to r_out = 1850 AU and is discretized over 10000 grid points. The large outer radius is chosen to not interfere with viscous spreading of the disk. The numerical resolution is variable and such that Δr/r ≃ 1.2 × 10⁻³. In this study, we assume that M_s = 1 M_☉, T_s = 4280 K, and R_s = 2 R_☉ (Siess et al. 2000). The initial surface density distribution of the gas obeys the relation Σ = Σ⁰₁(r₁/r)^β, where β = 1/2, 1, or 3/2, within at least ∼10 AU. Farther away from the star, Σ is exponentially tapered. The extremes of the "slope" β bracket values derived for the solar nebula by Davis (2005) and by Weidenschilling (1977) and Hayashi (1981). The quantity Σ⁰₁, the surface density at r₁ = 1 AU, is such that the initial disk mass is M⁰_D ≃ 0.022, 0.044, or 0.088 M_s. (These will be regarded as nominal values. The total initial disk mass differs somewhat for the different values of β because of the tapering procedure.) The photoevaporation rate (Equation (4)) is specified by imposing f₄₁ in Equation (5). Here we use f₄₁ = 0.1, 1, 10, 100, and 1000 (Alexander et al. 2005). The kinematic viscosity is ν = ν₁(r/r₁)^β and ν₁ = 4 × 10⁻⁶, 8 × 10⁻⁶, and 1.6 × 10⁻⁵ r²₁ Ω₁, where Ω₁ is the rotation rate at r = r₁. As a reference, in a disk with constant aspect ratio H/r = 0.04, ν₁ = 8 × 10⁻⁶ r²₁ Ω₁ corresponds to a turbulence parameter (Shakura & Syunyaev 1973) α_t = 0.005. The initial accretion rate onto the star ranges from a few times 10⁻⁸ to a few times 10⁻⁷ M_☉ yr⁻¹. For comparison, the mass-loss rate in Equation (6) is between ∼10⁻¹⁰ and ∼10⁻⁸ M_☉ yr⁻¹.

The majority of disk models have an initial gas inventory of at least ∼0.02 M_☉ within a distance of 40 AU from the Sun, as required by a canonical minimum mass solar nebula (MMSN; e.g., Weidenschilling 1977; Hayashi 1981). This value is also consistent with the more recent MMSN model adopted by Chiang & Youdin (2010). Due to the steepness of the surface density, disk models with the lowest initial mass and parameter β = 3/2 have only 0.01 M_☉ worth of gas within 40 AU of the Sun.

Gas is removed via the combined action of accretion onto the star, $\dot{M}$ (Equation (7)), and photoevaporation, $\dot{M}_{\mathrm{w}}$ (Equation (6)). In particular, Equation (6) sets an upper limit to dispersal timescale, τ_D, ranging from ∼1.4 Myr for M⁰_D ≃ 0.022 M_s (when f₄₁ = 1000) to ∼560 Myr for M⁰_D ≃ 0.088 M_s (when f₄₁ = 0.1). For computational purposes, τ_D is defined as the time past which M_D ≲ 10⁻⁵ M_s.

The evolution of the disk mass for some selected cases is illustrated in Figure 2 for each reference viscosity (see the figure caption for details). A complete list of the disk lifetimes is reported in Table 1. The behavior of the disk mass as a function of time, for the different surface densities, can be qualitatively understood in terms of viscous evolution by means of the analytic solutions of Lynden-Bell & Pringle (1974, their Section 3.3): for equally massive disks, the more compact the disk is (i.e., the larger β), the more rapidly M_D reduces initially. By a somewhat conservative assumption, as discussed above, disks that survive beyond 20 Myr are discarded and will not be given any further consideration. This is the case, for example, for all models with a photoionizing rate characterized by f₄₁ ⩽ 1 and the flattest initial surface density (β = 1/2). Models of disks surviving less than 1 Myr will also be discarded based on considerations on planet formation timescales, as explained in Section 4.

Zoom In Zoom Out Reset image size

A quantity of primary importance for planetary migration is the average surface density around the planet's orbit. In Figure 3 (left panels), the evolution of Σ is shown for cases with different values of parameters β and f₄₁. As anticipated at the end of Section 2.1, once the accretion rate drops below some threshold, photoevaporation produces a gap in the surface density at a radial distance of a few AU. Then the disk inside the gap is removed by viscous diffusion on a timescale of the order of r²/(2πν) orbital periods. We shall see in the Sections 3.2 and 3.3 that the disk's aspect ratio has also a large impact on the rates and direction of migration. In the right panels of Figure 3, H/r is plotted at reference times for the same models as in the left panels. Once Σ becomes small enough and the viscous heating term Q_ν in Equation (14) becomes unimportant, the disk temperature, and hence H/r (see Section 2.2), is dictated only by the stellar irradiation temperature (Equation (11)).

Zoom In Zoom Out Reset image size

The surface density at 1 AU, Σ₁, versus time is illustrated in Figure 4 for selected models from Table 1. Thin and thick lines refer to smallest and largest values of ν₁, respectively. To obtain a better statistical characterization of the disk density and temperature, for each pair of parameters (M⁰_D, β) listed in Table 1, ν₁ and f₄₁ are varied randomly in the corresponding ranges indicated in the table, for a total of more than 2000 realizations. The histogram in Figure 5 (left) shows that after ∼1 Myr the value of Σ₁ is 150 g cm⁻², or less, in ∼80% of the models that may represent the solar nebula. The right panel of Figure 5 illustrates the occurrence frequency of the mid-plane temperature at 1 AU. As a reference, the equilibrium temperature established by stellar irradiation alone, i.e., Equation (15), is T₁ ≈ 100 K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The migration of a Jupiter–Saturn pair in a gaseous disk is evaluated by using a combination of two-dimensional (2D) and three-dimensional (3D) hydrodynamical calculations of tidal interactions between the planets and the disk. We adopt a reference frame {O; r, θ, ϕ} with origin, O, fixed on the star, radius ranging from 0.25 to 7 AU, and azimuth varying between 0 and 2π. In the 2D disk approximation, the co-latitude angle θ is equal to π/2, whereas it varies from θ_min to π/2 in a 3D disk. In the latter case, the disk opening angle, θ_min, is such that the disk's vertical extent locally comprises at least three pressure scale heights, H. Mirror symmetry with respect to the θ = π/2 plane is imposed on account of the planets orbiting in this plane of symmetry. The surface density Σ initially has a dimensionless gradient dln Σ/dln r = −1/2. We work in the assumption that the disk is locally isothermal, i.e., the temperature depends only on r, and that H/r is a constant. The gas pressure is therefore proportional either to Σ/r (2D) or to ρ/r (3D), where ρ is the mass density. It is further assumed that the kinematic viscosity of the disk, ν, is constant throughout the disk.

The coordinate system rotates about the axis perpendicular to the planets' orbital plane (θ = 0 axis) at a variable rotation speed, $\mathbf {\Omega }_{f}=\mathbf {\Omega }_{f}(t)$. Both $\mathbf {\Omega }_{f}$ and $\mathbf {\dot{\Omega }}_{f}$ are imposed by the requirement that the (relative) azimuthal position of the interior planet, ϕ₁, remains constant in time and the (relative) angular velocity, $\dot{\phi }_{1}$, is zero (for details about the procedure, see D'Angelo et al. 2005).

Naming $\mathbf {r}_{\mathrm{1}}$ and $\mathbf {r}_{\mathrm{2}}$ the vector positions of the interior (Jupiter) and exterior planet (Saturn), respectively, the gravitational potential in the disk is

$$\begin{eqnarray} \Phi &=& {-}\frac{{\it G\mbox{$M_{s}$}}}{r} -\frac{{\it G M_{\mathrm{1}}}}{\sqrt{|\mathbf {r}-\mathbf {r}_{\mathrm{1}}|^{2} +\varepsilon ^{2}_{\mathrm{1}}}} -\frac{{\it G M_{\mathrm{2}}}}{\sqrt{|\mathbf {r}-\mathbf {r}_{\mathrm{2}}|^{2} +\varepsilon ^{2}_{\mathrm{2}}}} \nonumber \\ & & +\frac{{\it G M_{\mathrm{1}}}}{r^{3}_{\mathrm{1}}}\mathbf {r}\cdot \mathbf {r}_{\mathrm{1}} +\frac{{\it G M_{\mathrm{2}}}}{r^{3}_{\mathrm{2}}}\mathbf {r}\cdot \mathbf {r}_{\mathrm{2}}, \end{eqnarray} \tag{ 17 }$$

which accounts for the contributions of non-inertial terms due to the reference frame being centered on the star (see Nelson et al. 2000). The potential softening lengths ε₁ and ε₂ are set equal to 1/4 (or 1/7 in some calculations) times the Hill radius, R_H, of the corresponding planet. It is worth stressing that the argument according to which ε should be a fraction of the disk's scale height, H, in the 2D geometry (e.g., Masset 2002; Müller et al. 2012) applies to fully embedded planets, when $\mbox{$R_\mathrm{H}$}<H$. If $\mbox{$R_\mathrm{H}$}\gtrsim H$, as in all our 2D calculations, the local disk scale height depends also on the gravity of the planet itself (e.g., D'Angelo et al. 2003). In such cases, one physical constrain on ε is that it should be smaller than the radius over which gas is effectively bound to (i.e., it rotates about) the planet (${\sim} \mbox{$R_\mathrm{H}$}/3$; see, e.g., D'Angelo et al. 2003).

The Navier–Stokes equations that characterize the disk evolution are solved by means of the finite-difference code described in D'Angelo et al. (2005 and references therein) with modifications detailed below. The disk is discretized in 678 × 16 × 700 grid zones, in r, θ, and ϕ, respectively (and 678 × 700 in 2D). Calculations were also performed at a higher resolution of 1353 × 28 × 2096. Comparisons of the evolution of the planets' orbital elements at these two resolutions yield good agreement. We apply the wave-damping boundary conditions of de Val-Borro et al. (2006) within r = 0.3 and beyond r = 6.65 AU, which are appropriate for planets far enough from the boundaries.

We implemented an orbital advection algorithm along the lines of the FARGO algorithm of Masset (2000; see also Kley et al. 2009). These types of algorithms exploit periodicity properties of the flow, as those naturally occurring in the azimuthal direction of a disk. (Note that these algorithms can also be applied to local disk simulations, if periodicity is imposed at the patch boundaries; see Gammie 2001.) As demonstrated by Masset (2000), when the highest velocity component is along the periodic direction, in a 2D disk such techniques can increase the time step limit required by the Courant–Friedrichs–Lewy condition (see, e.g., Stone & Norman 1992), relative to a standard advection scheme, by factors ∼10 or larger. In a 3D disk, the gain may critically depend on the numerical resolution in the vertical (θ) direction.⁵ The implementation requires care when handling the transport of quantities defined on staggered meshes. Kley et al. (2009) use split cells, apply the transport algorithm to each part, and then recombine the partial information to reconstruct the full transport. Unlike them, we define the auxiliary variables required in the procedure (see Masset 2000 for details) on the same staggered meshes as the transport quantities are defined, wherever they are necessary. This approach requires more copies of the standard auxiliary variables to be defined, and hence more memory storage, but offers the advantage that the advection of all hydrodynamical quantities can be performed in a single step. Contrary to the implementations of both Masset (2000) and Kley et al. (2009), the algorithm used here avoids any directional bias, maintaining the full symmetry of the advection scheme, by using a sequence that alternates the transport among directions (see Stone & Norman 1992).

The equations of motion of two planets orbiting in a disk around a star, written in a reference frame rotating at variable angular speed, are

$$\begin{eqnarray} \mathbf {\ddot{r}}_{\mathrm{1}}&=& {-}\frac{G(\mbox{$M_{s}$}+M_{\mathrm{1}})}{r^{3}_{\mathrm{1}}}\mathbf {r}_{\mathrm{1}} -\frac{{\it GM_{\mathrm{2}}}}{r^{3}_{\mathrm{2}}}\mathbf {r}_{\mathrm{2}} -\frac{{\it GM_{\mathrm{2}}}}{r^{3}_{\mathrm{12}}}\mathbf {r}_{\mathrm{12}} +\mathbf {\mathcal {A}}_{\mathrm{1}}-\mathbf {\mathcal {A}}_{s} \nonumber \\ & & -\mathbf {\Omega }_{f}\mathbf \;{\times }\;(\mathbf {\Omega }_{f}\mathbf \;{\times }\; \mathbf {r}_{\mathrm{1}}) -2\,\mathbf {\Omega }_{f}\mathbf \;{\times }\; \mathbf {\dot{r}}_{\mathrm{1}} -\mathbf {\dot{\Omega }}_{f}\mathbf \;{\times }\; \mathbf {r}_{\mathrm{1}} \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} & & \nonumber \\ \mathbf {\ddot{r}}_{\mathrm{2}}&=& {-}\frac{G(\mbox{$M_{s}$}+M_{\mathrm{2}})}{r^{3}_{\mathrm{2}}}\mathbf {r}_{\mathrm{2}} -\frac{{\it GM_{\mathrm{1}}}}{r^{3}_{\mathrm{1}}}\mathbf {r}_{\mathrm{1}} +\frac{{\it GM_{\mathrm{1}}}}{r^{3}_{\mathrm{12}}}\mathbf {r}_{\mathrm{12}} +\mathbf {\mathcal {A}}_{\mathrm{2}}-\mathbf {\mathcal {A}}_{s} \nonumber \\ & & -\mathbf {\Omega }_{f}\mathbf \;{\times }\; (\mathbf {\Omega }_{f}\mathbf \;{\times }\; \mathbf {r}_{\mathrm{2}}) -2\,\mathbf {\Omega }_{f}\mathbf \;{\times }\; \mathbf {\dot{r}}_{\mathrm{2}} -\mathbf {\dot{\Omega }}_{f}\mathbf \;{\times }\; \mathbf {r}_{\mathrm{2}}, \end{eqnarray} \tag{ 19 }$$

where $\mathbf {r}_{\mathrm{12}}=\mathbf {r}_{\mathrm{1}}-\mathbf {r}_{\mathrm{2}}$. Note that since the origin of the coordinate system is on the star, Equations (18) and (19) include the forces per unit mass exerted on the star by the planets. The gravitational acceleration terms imposed by the disk, $\mathbf {\mathcal {A}}_{\mathrm{1}}$, $\mathbf {\mathcal {A}}_{\mathrm{2}}$, and $\mathbf {\mathcal {A}}_{s}$, are defined by Equations (8) and (9) of D'Angelo et al. (2005) and updated every hydrodynamical time step, Δt. Equations (18) and (19) are integrated numerically over Δt by means of a high-order Gragg–Bulirsch–Stoer extrapolation algorithm with order and step-size control (Hairer et al. 1993). The algorithm chooses automatically a suitable order at each (internal) step, which basically depends on the required tolerance of the solution error. We set a relative tolerance of 10⁻⁹ and an absolute tolerance of 10⁻¹⁴.

The basic mechanism that may allow a pair of resonant-orbit planets to experience a positive torque exerted by a gaseous disk and migrate outward was first described by Masset & Snellgrove (2001). Labeling with subscripts 1 and 2 the inner and outer planet, respectively, in order for this mechanism to be active, the following conditions must be fulfilled.

• 1.  
  The planet-to-star mass ratios (q_i = M_i/M_s) must be such that q₁ > q₂.
• 2.  
  The separation of the semimajor axes Δa = a₂ − a₁ must be such that Δa = b(R_{H, 1} + R_{H, 2}), where b ≲ 4.5 (as we shall discuss below).
• 3.  
  q₂ must be large enough to open a gap, or partial gap, in the density distribution by tidal torques.

Condition (2) above implies that $\Delta a/a_{1}=(\sqrt [3]{q_{1}}+\sqrt [3]{q_{2}})/ (\sqrt [3] {3}/b - \sqrt [3]{q_{2}})$. However, Hill stability for close planets on circular orbits also imposes that $\Delta a/a_{1}\gtrsim 2.40\, \sqrt [3]{q_{1}+q_{2}}$ (this inequality strictly applies in the limit of vanishing masses and absence of gas; see Gladman 1993), whose right-hand side is about equal to 0.26 for q₁ = M₁/M_s = 9.8 × 10⁻⁴ and q₂ = M₂/M_s = 2.9 × 10⁻⁴, hence b ≳ 2 (or 2.2, adopting a more precise determination of the Hill stability criterion; see Gladman 1993, for details; see also Figure 6). Conditions (1) and (2) suggest that, for a Jupiter–Saturn pair, the first encountered first-order mean motion resonance in which the mechanism may be activated is the 3:2 (the second-order 5:3 commensurability being another possibility), as indicated in Figure 6. In principle, configurations external, but sufficiently near, to resonances may also promote outward migration, as we shall see in Section 3.3. It is worth stressing here that simple capture in a mean motion resonance does not imply outward migration (see, e.g., Zhang & Zhou 2010), as we shall see in Section 7.

Zoom In Zoom Out Reset image size

Figure 7 (top and middle) illustrates the surface density perturbed by resonant-orbit planets, derived from 3D calculations for disks of different thicknesses (see the figure caption). The bottom panels of the figure show the mass density in a vertical slice of the disk, while the planets are aligned with the star. The exterior planet opens a partial gap in the case of a thinner disk, but it does so to a lesser extent in the other case (see bottom panels).

Zoom In Zoom Out Reset image size

The occurrence of outward migration of resonant-orbit planets can be intuitively understood from Figure 8, which reports on the results obtained from the same calculations as in Figure 7, where Saturn is caught in a 2:3 mean motion resonance with Jupiter. The azimuthally averaged surface density in normalized units (long-dashed lines) indicates that Saturn has cleared a partial gap, whose inner part overlaps with the outer part of the gap opened by Jupiter. The short-dashed lines in the figure represent torque density distributions (D'Angelo & Lubow 2008, hereafter DL08) due to Jupiter, $d\mathcal {T}/dM$. These functions yield the total torque when integrated over the disk mass. The torque exerted on the interior planet peaks at a₁ ± R_{H, 1} and is mostly comprised in a radial region of average width ∼3.5 R_{H, 1} on either side of the orbit (note that R_{H, 1} > H in the case displayed in the left panel and R_{H, 1} ≈ H in the other case). But since gas depletion due to gap formation may extend somewhat beyond this distance (see long-dashed lines in Figure 8), one may allow for a maximum value of the factor b in condition (2) above between 4 and 5. It is also important to notice that the orbital eccentricity of a planet acts to widen and smooth out gap edges (see, e.g., Figure 2 of D'Angelo et al. 2006), which may also affect somewhat the factor b.

Zoom In Zoom Out Reset image size

The solid lines in Figure 8 represent the cumulative torque, which is defined as

$$\begin{equation} \mathcal {T}_{\mathrm{CM}}(r)= 2\pi \int _{0}^{r}\frac{d\mathcal {T}}{dM}\langle \Sigma \rangle r^{\prime } dr^{\prime }. \end{equation} \tag{ 20 }$$

Looking at $\mathcal {T}_{\mathrm{CM}}$ in the left panel of Figure 8, it appears to be clear that the positive torque exerted by the disk interior of Jupiter's orbit is larger than that exerted by the disk exterior of the orbit, principally because Saturn has lowered the density there. The right panel illustrates the situation for a thicker disk, in which gas depletion operated by the exterior planet is not sufficient to reverse the sign of the torque. Figure 9 allows for a comparison of the cumulative torque for three different values of the disk thickness: H/r = 0.04 (top), 0.05 (center), and 0.07 (bottom) (see the figure caption for further details). A closer inspection of $d\mathcal {T}/dM$ and $\mathcal {T}_{\mathrm{CM}}$ in Figure 8 (left) and of the cumulative torque in Figure 9 (top) indicates that the total (positive) torque is basically driven by Lindblad resonances and that corotation torques are unimportant ($\mathcal {T}_{\mathrm{CM}}$ does not vary significantly over the radial width of the corotation region), as also argued by Morbidelli & Crida (2007).

Zoom In Zoom Out Reset image size

It thus appears from Figures 8 and 9 that the discriminant factor for outward, stalled, or inward migration is the depth (and width) of the outer planet's gap. If the outer planet opened a very deep and sufficiently wide gap, the inner planet would be subjected only to a one-sided Lindblad (positive) torque exerted by the interior disk that, to within a factor of order unity, can be written as (see Lubow & Ida 2010 and references therein)

$$\begin{equation} \mathcal {T}_{\mathrm{OS}}\sim a^{4}\,\Omega ^{2}\,\Sigma \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2} \left(\frac{a}{\widetilde{\Delta }}\right)^{3}, \end{equation} \tag{ 21 }$$

where $\widetilde{\Delta }=\max {(H,R_{\mathrm{H}})}$. The concept of one-sided torque is very useful to evaluate the presence of a tidally induced gap, i.e., condition (3) above. To first-order approximation, a density depletion begins to form when the one-sided torque exceeds the viscous torque (see Section 2.1), $\mathcal {T}_{\mathrm{OS}}\gtrsim \mathcal {T}_{\nu }$, which yields a simple order-of-magnitude condition $q^{2}\gtrsim 3\pi \alpha _{\mathrm{t}}(H/a)^{2}(\widetilde{\Delta }/a)^{3}$ (see also Papaloizou & Lin 1984; Ward & Hahn 2000), or

$$\begin{equation} g=\frac{q}{\sqrt{3\pi \alpha _{\mathrm{t}}}}\left(\frac{a}{H}\right) \left(\frac{a}{\widetilde{\Delta }}\right)^{3/2} \gtrsim 1. \end{equation} \tag{ 22 }$$

This conditions should be regarded as a measure of how much the density along the planet's orbit is depleted, hence it can be considered a condition for gas depletion. A condition for tidal truncation (gap formation) is then g ≫ 1 (see also Lin & Papaloizou 1986b). If predictions from the inequality (22) are compared with results from direct 3D calculations (DL08, Figures 6 and 8), one finds that g ≈ 1 corresponds to a ∼20% density depletion (relative to the unperturbed state, i.e., with no planet), and g ≈ 2.7 to ∼60% depletion. If applied to Saturn in the disks of Figure 8, g ≈ 3.4 for a density depletion of roughly 75% (left panels) and again g ≈ 1 for a density drop of ∼20% (right panels).

In order to derive migration rates, we shall assume that the orbits of the interior and exterior planets, Jupiter and Saturn, respectively, are in the 3:2 mean motion resonance and that this resonance is maintained during migration. Results from a calculation that support this assumption are plotted in Figure 10. All calculations resulting in outward migration behave similarly, but there are also instances in which the resonance is broken (see below). Since the ratio a₂/a₁ is supposed to be a constant, we concentrate on the migration rate of the interior planet, Jupiter. We seek an expression for $\dot{a}$ of the form

$$\begin{equation} \frac{d a}{d t}=\dot{a}_{\mathrm{ref}}\, k_{1}(\Sigma)\,k_{2}(a)\,k_{3}(g), \end{equation} \tag{ 23 }$$

where $\dot{a}_{\mathrm{ref}}$ is a reference migration speed, and k₁, k₂, and k₃ are dimensionless functions. To derive such expression, we employ results from both 2D and 3D calculations. In this section, the planets are fully formed since the beginning of the calculations and non-accreting.

Zoom In Zoom Out Reset image size

In Figure 11 (top-left panel), the semimajor axis evolution of the reference model is shown for both planets (see the figure caption for details). The disk has an aspect ratio H/r = 0.04 and turbulence parameter α_t = 0.005. At r = 1 AU, the surface density is Σ₁ = 50 g cm⁻² at time t = 0. This value is chosen from the disk evolution calculations discussed in Section 2.4, which show that ≳ 50% of disks have Σ₁ ≲ 50 g cm⁻² after ∼1 Myr. Following Pierens & Raymond (2011), we set a₁ = 1.5 AU and a₂ = 2 AU at t = 0 (see the caption of Figure 11 for further details), slightly outside the 3:2 mean motion resonance (see Figure 10). In the top-right panel of the figure, a comparison between 2D and 3D calculation results is presented. The 3D migration rate is about 25% smaller than the 2D one, which correction we apply to all 2D results. Note that this comparison is carried out at the same numerical resolution in the r – ϕ plane (see Section 3.1). The parameter that may produce the largest differences between 2D and 3D outcomes is the disk thickness. In this case, however, we rely only on 3D calculations to approximate Equation (23). Similar plots are shown in the bottom panels of Figure 11, but for the orbital eccentricity. For the duration of the evolution we consider, the eccentricity of Jupiter does not exceed ∼0.03 in any of the models discussed in this section. The orbital eccentricity of the exterior planet grows larger, to values e₂ ∼ 0.1 (see also Pierens & Raymond 2011).

Zoom In Zoom Out Reset image size

As mentioned above, Figures 10 and 11 indicate that outward migration can be activated also if orbital configurations are external, but somewhat near, the 3:2 commensurability. This effect is related to the gap widths and the extent to which density perturbations compound. In this context, orbital eccentricity may play some important role, since it affects the shape of a gap (D'Angelo et al. 2006).

In Section 3.2, we argued that the torque exerted on Jupiter would tend to the one-sided Lindblad torque $\mathcal {T}_{\mathrm{OS}}$ (Equation (21)), if Saturn (i.e., the exterior planet) carved a very deep and wide gap in the disk. In the opposite limit of very large disk thickness, H, and/or viscosity parameter, α_t, neither Jupiter nor Saturn would be capable of depleting the disk significantly and therefore it is expected that the torque exercised upon the interior planet will be of type I (Ward 1986; Lin & Papaloizou 1986a)

$$\begin{equation} \mathcal {T}_{\mathrm{I}}\sim -a^{4}\,\Omega ^{2}\,\Sigma \left(\frac{a}{H}\right)^{2} \left(\frac{\mbox{$M_{p}$}}{\mbox{$M_{s}$}}\right)^{2}, \end{equation} \tag{ 24 }$$

where, again, we neglect a factor (typically) of the order of unity in front of the right-hand side. Since both $\mathcal {T}_{\mathrm{OS}}$ and $\mathcal {T}_{\mathrm{I}}$ are linear in Σ, and since in the limit of zero orbital eccentricity

$$\begin{equation} \frac{da}{dt}=\frac{2\mathcal {T}}{a\Omega \mbox{$M_{p}$}}, \end{equation} \tag{ 25 }$$

one obvious guess is to approximate k₁ as a linear function. In the left panel of Figure 12, the migration velocity $\dot{a}$ of the inner planet from calculations (symbols), normalized to the velocity from the reference model (Figure 11), is plotted against the value of the normalized Σ. Here the value of the surface density is that at a distance of 5.5 R_{H, 1} from the (inner) planet's orbit and interior to it. A function proportional to Σ (solid line) appears a reasonable approximation of k₁ over the range of densities shown in the figure. Hence, we will assume that k₁(Σ) = Σ/Σ_ref, where the density is sampled as stated above.

Zoom In Zoom Out Reset image size

We note in passing that if the interior, and hence the exterior, planet is subjected to a torque of the type given in Equation (24), the resonance may be broken since the inner, more massive planet may drift inward at larger speed than the outer, less massive planet. This is indeed observed in some calculations of relatively thick disks, as illustrated in Figure 13 (see the figure caption for further details).

Zoom In Zoom Out Reset image size

The form of function k₂ can be guessed following a similar line of argument. In their natural units of a²Ω² (times a mass), both torques $\mathcal {T}_{\mathrm{OS}}$ (Equation (21)) and $\mathcal {T}_{\mathrm{I}}$ (Equation (24)) scale as a²Σ (the dependence on H/a is considered later), which can be seen as a measure of the local disk mass. Therefore, if Σ was constant, $\dot{a}\propto a^{2}$ in units of aΩ (see Equation (25)) and thus $\dot{a}$ would scale as a^3/2. In the right panel of Figure 12, the migration speed of the inner planet from calculations (symbols), normalized to the reference migration speed, is plotted as a function of a, normalized to the semimajor axis in the reference model, a_ref. Also in this case, the approximation seems satisfactory (solid line) and so we shall assume that k₂(a) = (a/a_ref)^3/2.

In Section 3.2, it was anticipated that the depth and width of the density depletion produced by the exterior planet play a fundamental role in determining magnitude and direction of the interior planet's migration, and hence of the pair as a whole. Quantity g (Equation (22)) can be used as a proxy to discriminate among the various situations, i.e., different combinations of H/r and α_t (and q) that may affect Σ. On account of the compact orbital configuration, since q₁ > q₂ we have that g₁ > g₂ (unless H/r and/or α_t are rapidly varying with disk radius). If g₂ ≫ 1, then g₁ ≫ 1, and the interior planet is likely subjected to a torque whose limit is the one-sided Lindblad torque in Equation (21) and the two planets migrate outward. Otherwise, if g₁ ≪ 1, then g₂ ≪ 1, and migration will be dictated by a type I torque (Equation (24)) and be directed inward.

The migration velocity of the inner planet should depend on both g₁ and g₂. However, in the present context g₂ should have the larger impact of the two and therefore, for the sake of simplicity, we assume that function k₃ depends only on g₂. In particular, referring to the value of g₂ in the reference model (Figure 11) as g_ref and indicating g₂ simply as g, k₃ will be approximated as a function of g/g_ref. Since it is unclear how the agreement between 2D and 3D calculations varies as a function of the disk thickness (it should likely worsen as H increases), we only use 3D calculations to find an approximation to function k₃. Figure 14 shows the ratio $\dot{a}/\dot{a}_{\mathrm{ref}}$ obtained from calculations for various values of the ratio g/g_ref. The thick horizontal line in the plot indicates the type I migration speed that would apply to the inner planet if H/r = 0.1, corresponding to the value used for leftmost data point on the graph. For reference, the normalized migration speed corresponding to the one-sided torque in Equation (21) would be ∼90. The broken line is a linear interpolation of the data, which will be used as a representation of k₃.

Zoom In Zoom Out Reset image size

Summarizing the results of Section 3.3, we write the migration speed of the interior planet as

$$\begin{equation} \frac{d a}{d t}=\dot{a}_{\mathrm{ref}} \left(\frac{\Sigma }{\Sigma _{\mathrm{ref}}}\right) \left(\frac{a}{a_{\mathrm{ref}}}\right)^{3/2} k_{3}(g/g_{\mathrm{ref}}), \end{equation} \tag{ 26 }$$

where the dimensionless function k₃ is obtained via linear interpolation of the numerical data in Figure 14, i.e., the thick solid line in the figure. The reference values in Equation (26) are taken from the reference model, corrected for 3D effects, as discussed in Section 3.3: $\dot{a}_{\mathrm{ref}}= 2.7\times 10^{-6}\,\mbox{AU}\,\mathrm{yr}^{-1}$ for Σ_ref = 42 g cm⁻² and a_ref = 1.57 AU. Recall that both Σ and Σ_ref are evaluated at a distance of 5.5 R_{H, 1} interior of the inner planet's orbit. The argument of function k₃ is given by

$$\begin{equation} \frac{g}{g_{\mathrm{ref}}}=\sqrt{\frac{\alpha _{\mathrm{t,ref}}}{\alpha _{\mathrm{t}}}} \left(\frac{H_{\mathrm{ref}}}{H}\right) \left(\frac{R_{\mathrm{H}}}{\widetilde{\Delta }}\right)^{3/2}, \end{equation} \tag{ 27 }$$

where $\alpha _{\mathrm{t,ref}}=0.005\sqrt{1\,\mbox{AU}/a_{\mathrm{ref}}}$, (H/a)_ref = 0.04, and $\widetilde{\Delta }=\max {(H,R_{\mathrm{H}})}$. Recall that here parameters g and g_ref are computed from Equation (22) are applied to the exterior planet. For a mass ratio q different from that of the reference model, the right-hand side of Equation (27) should be multiplied by q/q_ref. Equation (26), without the correction due to 3D effects, is also in reasonable agreement with the results presented by Pierens & Raymond (2011) in their Figure 21, for a disk of 0.4 M_J within 1.5 AU (Σ ∼ 500 g cm⁻² at 1 AU). As explained above, the exterior planet's orbit may not always be resonant with that of the interior planet, when migration is inward (see Figure 13). Nonetheless, for the outer planet's orbit we set a₂ = (3/2)^2/3a₁.

Equation (26) predicts stalling points (where k₃, and hence $\dot{a}$, is ≈0) at g₀ ≈ 0.8 g_ref, which will be regarded as a nominal value. But, notice that since the approximation to k₃ is sampled at a limited number of points, its zero could be located at a somewhat different abscissa, yet it is located between 0.7 g_ref and g_ref. In principle, there could be multiple stalling points in a disk (if disk properties are not monotonic), which would represent locations of stable equilibrium since they are convergent radii for the planets' semimajor axes. The argument, however, is based on the assumption that the 3:2 commensurability is preserved even for inward migration. This is not true in general (as illustrated, for example, in Figure 13), but it further requires that the inward migration of the exterior planet be faster than that of the interior planet and that the condition for gap overlap (see condition (2) in Section 3.2) be satisfied. When these conditions are not met, the pair of planets can migrate past a stalling point, toward the star. This may happens, for example, if disk temperature and viscosity are low enough so that both planets open up a gap and drift according to type II migration before capture into resonance occurs.

Some constraints on the range of outward migration predicted by Equation (26) can be derived for the disk models discussed in Section 2.4. The disk thickness, H, is affected by internal (viscous) heating typically over the first few million years of evolution (see Figure 3, left). If we neglect that source of heating, H can be approximated by using Equation (16) and then Equation (27) can provide a rough measure of a disk's radial range over which the planets can drift away from the star, as shown in Figure 15. The addition of internal heating would raise the value of H, hence reducing the ratio g/g_ref, which suggests that outward migration may not be activated in the warm interiors of a young disk. The radial region over which a curve extends above or, to some degree, inside the shaded area is favorable to outward migration. The intersection of a curve with the horizontal line in the shaded area gives the nominal radius of the stalling point of each planet (see above). According to the plot, outward migration of a Jupiter-mass planet locked in the 3:2 mean motion resonance with a Saturn-mass planet cannot proceed beyond ∼7 AU (top x-axis). The numerical experiments discussed in Section 4 are broadly consistent with this prediction.

Zoom In Zoom Out Reset image size

The validity of Equation (26) for planet-to-star mass ratios different from those adopted here (q₁ = M₁/M_s = 9.8 × 10⁻⁴ and q₂ = M₂/M_s = 2.9 × 10⁻⁴) was not investigated. Parameter g (Equation (22)) is ∝q₂ if $\mbox{$R_\mathrm{H}$}<H$ and ${\propto} \sqrt{q_{2}}$ otherwise. Therefore, a faster outward migration might be expected as q₂ increases, provided that the ratio q₁/q₂ stays roughly constant. Condition (2) of Section 3.2, along with the requirement of Hill stability, suggests that, as q₂ increases beyond ∼0.001, the 2:1, rather than the 3:2, commensurability may be available to activate outward migration for a ratio q₁/q₂ ≈ 3, in accord with the findings of Crida et al. (2009). This is schematically illustrated in Figure 16 (see the figure caption for further details). As the ratio q₁/q₂ approaches 1, the shaded area in the graph shifts upward.

Zoom In Zoom Out Reset image size

We shall consider here the possibility that Saturn forms within the 1:2, but outside the 2:3, commensurability with Jupiter, while both planets are beyond several AU from the Sun. If during the course of its evolution Saturn remains inside the 1:2 mean motion resonance with Jupiter, then capture in the 2:3 orbital resonance is still possible.

In order to maintain such a compact orbital configuration throughout the evolution of the two planets, the migration rates must be very similar over time. In fact, if the constraint 1.31 < a₂/a₁ < 1.59 has to be preserved, a change of this ratio of at most 20% over the formation timescales essentially implies that a₂/a₁ is roughly constant. Therefore, by taking the time derivative, one has

$$\begin{equation} \frac{\dot{a}_{2}}{a_{2}}\sim \frac{\dot{a}_{1}}{a_{1}}. \end{equation} \tag{ 35 }$$

But since the interior planet has to grow faster than the exterior planet does, their migration rates are bound to differ, at some point in time at least. Thus, the condition in Equation (35) is unlikely to be (always) satisfied and the difference a₂ − a₁ will either increase or decrease.

If Saturn forms within the 1:2 orbital resonance and Jupiter has still to acquire most of its mass, there will be a phase when the migration rate of Jupiter significantly exceeds that of Saturn (presumably, around the time of runaway gas accretion; see DL08). Hence, it is most likely that Saturn is left behind and probably becomes trapped in the 1:2 commensurability. Instead, if Jupiter has already acquired the bulk of its mass, the opposite is likely to occur and most probably Saturn becomes locked in the 2:3 orbital resonance while it is still growing.

We consider this last situation in some detail by performing 3D calculations similar to those of the reference model of Section 3.3, but applying a lower mass to the exterior planet: q = M₂/M_s = 10⁻⁴ and 2 × 10⁻⁴. Equation (27), appropriately modified for mass ratios q ≠ q_ref (see the discussion after Equation (27)), suggests that if q ≈ 10⁻⁴ migration of the pair is inward, but if q ≈ 2 × 10⁻⁴ then migration is outward. Direct calculations agree with these predictions, as indicated by the cumulative torques shown in the top-left panel of Figure 22 for three different values of M₂ (see the caption for details). The density depletion due to the tidal torques of the exterior planet amounts to ∼35% for M₂ = 10⁻⁴ M_s and to ∼60% for M₂ = 2 × 10⁻⁴ M_s (top-right panel), in accord with the expectations of Equation (22). The bottom panels of Figure 22 show vertical distributions of the mass density (see the figure caption).

Zoom In Zoom Out Reset image size

Therefore, if Saturn grows while locked in the 2:3 orbital resonance with Jupiter, there is a mass smaller than Saturn's final mass for which migration stalls and then reverses, as suggested by $\mathcal {T}_{\mathrm{CM}}$ as function of M₂ in Figure 22 (top-left panel). This situation would prevent Jupiter from reaching the inner disk regions, unless Saturn achieved the mass for migration reversal when Jupiter is already there. For H/r ≈ 0.04, such mass is between 10⁻⁴ M_s and 2 × 10⁻⁴ M_s (and presumably closer to 10⁻⁴ M_s for H/r ≈ 0.03), which implies that a substantial fraction of Saturn's envelope is acquired in the inner regions of the solar nebula.

This circumstance, however, would likely be at odds with the elemental abundances of some species measured in Saturn's atmosphere (see Hersant et al. 2008). The abundances relative to hydrogen of elements such as C, N, S, As, and P are a few to several times as high as the solar abundances (see Lodders 2003; Asplund et al. 2009 and references therein). In fact, the presence in large amount of these elements is believed to have arisen from accretion of gas (e.g., Guillot & Hueso 2006) and/or solids (e.g., Hersant et al. 2008) in a cold disk environment.