© ESO, 2015

1. Introduction

The interaction of the growing protoplanet with the ambient disk leads to a change in the orbital elements of the planet. The most important for the overall evolution of the planet is the change in semi-major axis, i.e. the migration of the planet. The topic of planet-disk interaction and the orbital evolution of planets has been covered in a few recent reviews (Kley & Nelson 2012; Baruteau & Masset 2013; Baruteau et al. 2013), and here we present only a brief summary of the relevant issues. Depending on the mass of the planet, different regimes of migration are distinguished. Most important are the two limiting regimes of types I and II migration.

Type I migration occurs for low mass planets (with a mass less than about 50 M_⊕) that do not open a gap in the disk, and can be treated in the linear regime. Here the total torque is given by the effects of the spiral density waves (Lindblad torques) and the corotation torques generated by the gas flow in the coorbital horseshoe region. In contrast to the Lindblad torque that quite generally leads to inward migration (Ward 1997), the corotation torque suffers from saturation effects making it strongly dependent on on the magnitude of viscosity and thermal diffusion. These can significantly alter the migration speed and may even reverse the direction of the migration. Analytical formulae to account for these effects have been presented by Paardekooper et al. (2010) for adiabatic disks and Paardekooper et al. (2011) for disks with thermal diffusion.

Type II migration occurs for massive planets with masses comparable to M_J or larger. Because of angular momentum deposition in the disk, those planets open an annular gap in the protoplanetary disk at the location of the planet where the density is significantly reduced. In this case the migration speed is slowed down from the linear rate because of the reduced mass available near the planet. During the migration process the gap created by the planet will have to move with the planet through the disk. Hence, it is often assumed that in an equilibrium situation the gap moves exactly with the viscous accretion velocity of the disk, and that the planet is locked in the middle of the gap to maintain torque equilibrium (Ward 1982, 1997; Lin & Papaloizou 1986).

Direct numerical simulations of planet-disk interaction typically place the planet in a non-accreting disk and keep the planet at a fixed position to analyse the torque density distribution from which the migration speed can then be calculated, see Kley & Nelson (2012) for references. Calculations, where gap-opening Jupiter-type planets were allowed to move through the disk as given by the torques acting on them, have been performed for isothermal disks by Nelson et al. (2000) for two-dimensional (2D) simulations. They studied non-accreting and accreting planets where mass from inside the planet’s Roche lobe was added to the planet mass. Both cases lead to similar results indicating that a massive planet could migrate within 10⁵ years from about 5 au all the way to the center. In full 3D radiative simulations even faster migration rates were found, though only for planets with a maximum mass of about 0.6 M_J (Bitsch & Kley 2010). Numerical calculations of migrating massive planets were carried out by Edgar (2007, 2008). Using a constant kinematic viscosity he came to the result that there is no constant migration rate as predicted for the type II regime. In a related study Paardekooper (2014) analyzed the influence of planetary motion on the type I migration regime. He found that in non-isothermal disks the regime of outward migration can be larger than that estimated for fixed planets.

The migration of massive planets in evolving disks has been analysed by Crida & Morbidelli (2007) for 2D isothermal disks. In their models, the planets were placed in a disk with an initial Gaussian profile that evolves under a constant kinematic viscosity. Their results show that only planets that carve a deep gap in the disk experience genuine type II migration and move with the viscous speed of the disk. This is the case either for very massive planets or low viscosity disks, with a Reynolds number larger than about 5 × 10⁵. In the case of only partial cleared gaps, i.e. for less massive planets and larger viscosity, they found a reduced migration speed and in some cases even outwards migration, an effect due to the very effective action of the corotation torques. In addition, because of the positive density slope at the location of the planet in their simulations, the Lindblad torques are strongly reduced, which strengthens the effect.

Following this global study there have been recent models of planets embedded in disks with a given (constant) mass accretion rate. Bitsch et al. (2014) focussed on the regime of earth-mass planets in an irradiated accretion disk with a given Ṁ rate to determine the regimes of inward and outward migration. They modelled axisymmetric disks with vertical structure and used the formula by Paardekooper et al. (2011) to determine the migration properties of the embedded planets. Fung et al. (2014) investigated the structure of gaps created by stationary planet in an accreting disk, but they did not include planet migration. In their 2D isothermal studies they showed that despite the presence of a massive, gap opening planet, which was not allowed to accrete any material, the disks reached an equilibrium state with a constant mass accretion rate. Obviously, in this configuration the mass flow across the gap equals exactly the imposed disk accretion rate, ṁ. The timescale of type II migration was analysed with respect to the agreement with the orbital properties of the observed massive extrasolar planets by Hasegawa & Ida (2013). They argued that standard type II migration is too rapid and mechanisms need to be found to slow it down.

Recently, Duffell et al. (2014) looked at the migration of planets in accretion disks using an alternative, different approach. Instead of moving the planet according to the disk torques, they pulled a Jupiter mass planet through a zero ṁ non-accreting disk, measured the torques acting on the planet, and compared the resulting migration rate with the actual pull rate of the planet. In this way they obtained possible equilibrium solutions of the migration speed of a planet through a disk. In particular, they found that the possible equilibrium migration speed of the planet is independent of the viscous speed and can be lower and smaller than type II migration, depending on the local surface density of the disk.

Hence, it appears that the issue of type II migration is presently not resolved. Theoretically, it is clear that a planet can only be moved by the disk torques acting on it. On the other hand, disk material moves under the action of viscous torques and its local inflow speed is directly proportional to the viscous torque. If a planet migrated exactly in the type II regime, then the disk torque and viscous torque should be equal. In this paper, we make a new attempt at the problem of type II migration and study the migration of planets in 2D isothermal disks with a constant mass accretion rate. In a first step, we construct constant ṁ disks with stationary planets of different masses and then, in a second step, we move the planets according to the torques acting on it, measure their migration rate, and compare it to the type II migration rate.

Our setup is described in Sects. 2, and 3 we present our results for fixed planets in accreting disks. In Sect. 4 we move the planets according to the torques acting on them, and compare this in detail to type II migration in Sect. 5. The results are discussed in the final Sect. 6.

2. Setup

To study the migration of planets in disks, we assume that the disk is geometrically thin and simplify to a two-dimensional (2D) approximation. The disk is assumed to be locally isothermal and driven by an α-type viscosity. In addition we model explicit mass accretion through the disk. For our calculations we use the NIRVANA-code (Ziegler & Yorke 1997; Ziegler 1998) in a 2D setup. The star is located at the centre of a cylindrical coordinate system, which covers a radial range of 1.56 to 15.6 au corresponding to r = 0.3...3.0 in code units where the unit of length is given by r₀ = 5.2 au. In the following we will specify distances in code units unless specified otherwise. In the azimuthal direction we cover a complete annulus from 0 to 2π. In our standard resolution the domain is covered by 251 × 583 cells. The planet is placed on a circular orbit at a distance of 5.2 au corresponding to r = 1.0.

To speed up the initial relaxation we, first calculated the equilibrium of the disk, without the planet, with a reduced resolution of (N_r,N_φ) = (101,233). After 10⁶ time steps, corresponding to about 5000 orbits, these results were interpolated to our standard-resolution grid. Then the model could adapt to this new resolution for additional 10⁵ time steps (about 290 orbits) before the planet was released. The planet then could move freely in the disk, and change its semi-major axis and eccentricity. In addition, we performed another set of similar calculations with a lower resolution of (N_r,N_φ) = (135,405) and extending from r = 0.2 to 2.0. Both calculations give similar results, after 1000 orbits the positions of the planets differ by around 5%.

For our studies we varied the viscosity by changing the α-parameter, the accretion rate (and thus the disk surface density, see Eq. (4)below), and the planetary mass. The used parameter space with the highlighted standard model, is shown in Table 1. We limit our planet mass to a maximum of 2 M_Jup (or q = 0.002) because for larger m_p the outer disk becomes eccentric (Kley & Dirksen 2006) and migration properties change (D’Angelo et al. 2006).

2.1. Initial and boundary conditions

The goal of this work was to set up a disk with a steady accretion flow through the disk. Hence, matter has to be fed into the domain at the other boundary of the disk (at r_max), transported to the inner disk, leaving the domain at r_min. The local accretion rate through a disk is given by$\mathit{ṁ}\mathrm{=}\mathrm{-}\mathrm{2}\mathit{\pi r}\mathrm{\Sigma }{\mathit{u}}_{\mathit{r}}\mathit{.}$(1)It depends on the surface density, Σ, and the radial velocity u_r, of the gas. In equilibrium for a constant ṁ the viscous accretion velocity is given by${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\frac{\mathit{\nu }}{\mathit{r}}\mathit{,}$(2)where ν is the kinematic viscosity and r the radial distance to the central star. This result can be obtained for stationary accretion disks with constant accretion rate, derivations can be found in textbooks such as Frank et al. (1992) or Armitage (2010). In our case of an isothermal disk, we use the α-viscosity prescription, ν = αc_sH, (Shakura & Sunyaev 1973) with the disk scale height H = hr, where the relative disk thickness is constant, h = 0.05. With the isothermal sound speed, c_s = HΩ_K, this can be written as${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\frac{\mathit{\alpha }{\mathit{c}}_{\mathrm{s}}\mathit{H}}{\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{\alpha }{\mathit{h}}^{\mathrm{2}}\mathit{r}{\mathrm{\Omega }}_{\mathrm{K}}\mathit{,}$(3)where ${\mathrm{\Omega }}_{\mathrm{K}}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{/}{\mathit{r}}^{\mathrm{3}}}$ is the Keplerian orbital frequency with the stellar mass M_∗ and the gravitational constant G. With given ṁ and α this leads in equilibrium to the surface density$\mathrm{\Sigma }\mathrm{=}\frac{\mathrm{1}}{\mathrm{3}\mathit{\pi }}\frac{\mathit{ṁ}}{\mathit{\alpha }{\mathit{h}}^{\mathrm{2}}{\mathit{r}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{K}}}\mathrm{=}\frac{\mathit{ṁ}}{\mathrm{3}\mathit{\pi \alpha }{\mathit{h}}^{\mathrm{2}}\sqrt{\mathit{GM}}}\hspace{0.17em}{\mathit{r}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\mathit{r}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(4)From radial equilibrium (pressure gradient balanced by centrifugal force and gravity) the angular velocity is given by${\mathit{u}}_{\mathit{\phi }}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}\frac{{\mathit{h}}^{\mathrm{2}}}{\mathrm{4}}}\hspace{0.17em}\mathit{r}{\mathrm{\Omega }}_{\mathrm{K}}\mathit{.}$(5)Thus, for our initial conditions we use ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}$ from Eq. (3), Σ from Eq. (4)and u_φ from Eq. (5), with prescribed values of ṁ, α and h.

At the outer boundary gas enters the computational domain with ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}$, and Σ₀ is fixed according to the chosen ṁ. For the standard model with α = 0.003 and ṁ = 10^-7 M_⊙/ yr, the initial disk density at r = 1 is Σ₀ = 878 g cm^-3. At the inner boundary only the radial velocity is fixed according to Eq. (3), material then flows out with the local density at r_min. The azimuthal velocity u_φ is held fixed, according to Eq. (5), at both the inner and the outer boundary. In the azimuthal direction, we enforce periodic boundaries. For the gravitational potential, we use ϵ-smoothing with ϵ = 0.6H (Müller et al. 2012). To prevent density fluctuations due to wave reflection at the radial boundaries the radial and azimuthal velocity is damped towards its azimuthal mean value near the boundaries on a timescale of ${\mathit{\tau }}_{\mathrm{damp}}\mathrm{\approx }{\mathrm{\Omega }}_{\mathrm{K}}^{-1}$ with Ω_K taken at the inner and outer boundaries.

2.2. Accretion onto the planet

To study the influence of mass accretion by the planet on the migration speed, we performed additional simulations in the two extreme cases of non-accreting planets and maximally accreting planets. Both cases are of course unrealistic, but reality lies somewhere in between. To model maximal accretion, we implemented the accretion scheme devised by Kley (1999) where at each time step the gas density inside the Roche Lobe R_R ≈ a(q/ 3)^1/3 of the planet is reduced by a factor of 1 − f_redΔt with f_red = 1/2. In our simulations, the mass removed in this way is not added to the planet mass.

3. The disk structure for fixed planets

Fig. 1 Surface density of a fixed-orbit, non-accreting planet calculation for the standard model with (α,Mp/M∗,ṁ) = (0.003,0.001,10-7) after 2500 orbits. The planet is in the middle of a stable gap and the accumulation of mass near the planet is clearly visible. The spiral arms are damped at the inner and outer boundary so no reflections are seen.

Before starting the calculations with moving planets, we constructed equilibrium cases where the planets are not allowed to move and remain at their initial locations. This allows us to estimate the size and depth of a gap as a function of planet mass and the viscosity. The surface density of this calculation for the standard model after 2500 orbits is shown in Fig. 1. There is a very steady gap edge and stable spiral arms resulting from the perturbation of the planet.

Fig. 2 Azimuthally averaged surface density for a non-accreting (solid line) and an accreting (dashed line) planet with q = 0.001, for different values of the viscosity parameter α. The black line indicates the initial density profile at the beginning of the simulation. The red line corresponds to the standard model.

Fig. 3 Azimuthally averaged surface density for non-accreting (solid line) and accreting (dashed line) planets and α = 0.003, for different values of the planet mass. The red line corresponds to the standard model.

The azimuthally averaged surface density in the case of stationary planets is shown for different values of α (in Fig. 2) and q (in Fig. 3). As expected, the gaps are deeper for lower viscosity and higher planet mass. While a more massive planet causes a much wider gap, decreasing the viscosity primarily deepens the gap. In isothermal disks, for stationary planets, the surface density distribution is identical for different values of the accretion rate, and hence ṁ has no influence on the gap profiles that are therefore identical. The gap depth depends on the accretion onto the planet. For non-accreting planets the gap has a clearly visible density bump at the position of the planet. If accretion onto the planet is allowed, the bump disappears and the gap becomes deeper because gas is accreted onto the planet and removed from the simulations.

In the outer disk region, for r> 1.7, the surface density is increased over the case without a planet and the amount increases with planet mass. This is because the gas initially in the gap is pushed to the inner and outer regions of the disk, because of angular momentum transfer to the disk. Because, in the outer disk, the gas cannot leave the computational domain the surface density must be increased, the effect being stronger for more massive planets with a deeper gap. Hence, there occurs a jump in the surface density at the outer boundary. This is also clearly seen in Fig. 4 where no logarithmic scale is used. This jump in Σ at r_max has no influence on our results on the migration properties and it disappears if larger r_max are used.

Fig. 4 Azimuthally averaged surface density profile (top) and the migration rate profile (bottom) in case of (α,q,ṁ) = (0.003,0.001,10-7) as in Fig. 1. The planet is at r = 1.0. The profiles are taken at different times after the beginning of the calculation. The black line indicates the initial density profile at the beginning of the simulation before the planet was added.

The evolution towards the equilibrium state is shown in Fig. 4 where the radial dependence of the surface density and local accretion rate, ṁ(r), are displayed at different times. The migration rate is obtained by integrating the mass fluxes in the advection routine of the code over the φ-direction yielding a radial profile of accretion. This ensures consistency with the numerical calculations. In the bottom panel, it is clearly visible that the equilibrium is not yet reached after 3600 orbits. This is because the accretion rate is extremely sensitive to all kinds of perturbations. The equilibrium viscous inflow speed is given by u_r/c_s ≈ αh, which is about 10^-4 for our standard model. Hence even small pressure perturbations lead to large variations in the radial velocity, compared to the viscous speed. Thus, although the accretion rate profile is not flat, the changes in the surface density after 2700 and 3600 orbits are very small, and the gap is already fully developed. Obviously, in a steady-state solution the gas must be able to cross the gap. The inner and outer disk are connected and there is no significant pile up at the outer gap edge or depletion at the inner gap edge.

Fig. 5 Normalized torques acting on planet with fixed orbits at the end of the relaxation as a function of the planet mass. The data is taken from the low-resolution calculations. For isothermal disks the normalized torques are independent of the disk density.

Fig. 6 Planet migration tracks for different model parameter. Starting from the standard model (in red) with (α,ṁ,q) = (0.003,10-7,0.001), individual parameters have been varied while keeping two others fixed. We varied α (left panel), ṁ (middle panel), and q (right panel). The full and dotted lines correspond to non-accreting and accreting planets, respectively.

3.1. The torque acting on the planet

To estimate the expected migration rate of planets in disks with mass flow, we calculate the torques generated by the density perturbations in the disk as induced by the planet. In Fig. 5 we display the obtained torques as a function of planet mass M_P for different value of the viscosity. Here, we used the standard torque normalization (see e.g. Paardekooper et al. 2010)${\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{=}\mathrm{-}{\mathrm{\Sigma }}_{\mathrm{p}}{\mathrm{\Omega }}_{\mathrm{K}}^{\mathrm{2}}{\mathit{a}}_{\mathrm{p}}^{\mathrm{4}}\hspace{0.17em}{\left(\frac{\mathit{q}}{\mathit{h}}\right)}^{\mathrm{2}}\mathit{,}$(6)where Σ_p is the unperturbed surface density at the location of the planet and a_p its semi-major axis. This definition of Γ₀ takes the linear results into account for low mass planets that are subject to type I migration with $\mathrm{\Gamma }\mathrm{\propto }{\mathit{M}}_{\mathrm{p}}^{\mathrm{2}}$, massive planets open gaps in the disk that reduce the torque, and a different scaling is expected.

To strengthen the point that these torques are measured for fixed planets that are not allowed to migrate, we introduce the correction factor f_still and write for the torque acting on the planet$\mathrm{\Gamma }\mathrm{=}{\mathit{f}}_{\mathrm{still}}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathrm{0}}\mathit{.}$(7)As seen from Fig. 5, the total torque is indeed reduced for larger planet masses. For small planet masses, a flattening is expected as one approaches the type I regime. In our case, the planet masses are still too high, but this effect becomes marginally visible for the higher viscosity models, as expected. On the other hand, there is no visible dependence on variations in the disk mass or the accretion rate (which are coupled here), because in isothermal disks the disk mass is just a scaling factor that depends linearly on the specified ṁ value.

4. Migrating planets

After having analysed the equilibrium state for fixed planets, we release now the planets and evolve their orbital elements according to the action of the disk. We measured the migration tracks for many different parameters (see Table 1) for more than 1000 orbits for accreting and non-accreting planets. The migration tracks are shown in Fig. 6 for a variety of model parameters. All planets migrate inwards at a rate that decreases with time. At the beginning of the migration, shortly after the release of the planet, there is a short phase of rapid inward migration that then slows down. This transient phenomenon is caused by the unavoidable small mismatch of a disk structure with a stationary planet in comparison to the evolving case.

In cases where ${\mathit{M}}_{\mathrm{P}}\lesssim {\mathit{M}}_{\mathrm{D}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}$, with Σ_P being the undisturbed surface density at the position of the planet a_P, there is no slowing down, rather we observe very fast type III migration, such as for q ≤ 0.0005 (yellow and blue curve) in the right panel of Fig. 6. Type III migration sets in later for accreting planets.

4.1. Mass accretion onto the planet

Our implementation of the mass accretion onto the planet is very simple, and, therefore, we only used this method to obtain an idea what changes are expected in case of an accreting planet. Because the algorithm takes away the gas near the planet at a high rate, it serves as an upper limit of the possible accretion onto the planet. We found that the general behaviour of accreting planets is similar to the non-accreting cases. As seen in Fig. 6 accreting planets migrate more slowly in most cases. The reason is the reduced density in the vicinity of the planet, which leads to weaker torques. Although the contribution of the planet envelope inside 0.8 Hill-Radii is not considered when calculating the torque (see Crida et al. 2008), the region affected by the accretion is much bigger and thus important for the torques. This effect of accretion is similar to what has been seen in Nelson et al. (2000) who added the accreted material to the dynamical mass of the planet, however. This will lead to an even slower migration of the planet due to the increased inertia.

4.2. Gap profiles during the evolution

The gap profile has an important impact on the migration of a planet. Because the disk is depleted near the planet’s orbit, these regions cannot contribute to the torques on the planet. In Fig. 7, the change in global disk density is displayed for four different snapshots during the evolution of the planet in our standard model. At the inner edge of the gap, we observe a density pile up that initially increases with the planet moving inwards closer to the inner edge. But when the planet reaches r ≈ 0.7 at t = 530, it decreases again because the gap approaches r_min, and at the inner boundary an outflow condition with a given radial velocity is enforced. At the outer gap edge we see the gas lagging behind the movement of the gap. The gap moves Δr = 0.3 in 530 orbits. The viscous timescale τ_visc = Δr²/ν for the gas to cover this distance is about 2070 orbits. It is therefore obvious that the planet and the gap move faster than the viscous speed of the gas, which demonstrates that the gap must be dynamically created during the planet’s inward motion. The inward motion of the planet and the gap can only be maintained because gas can cross the gap faster than the viscous timescale.

To obtain a more detailed view on the local disk structure in the gap region, we display in Fig. 8 the gap at different times during the evolution of the planet in the disk. The gap does not change significantly during the migration through the disk if the width is rescaled with the local disk scale height H(a_p) at the actual position of the planet. Farther away from the planet the surface density varies because of the different positions in the disk.

Fig. 7 Global gap profiles at four different times during migration of the planet for the standard model with q = 0.001 and ṁ = 10-7. The black line indicates the unperturbed density profile without the planet.

Fig. 8 Gap profiles for the standard model taken at different positions of the planet in the disk during its migration path. The gap is rescaled using the local disk scale height, H, at the position of the planet. The curves correspond directly to those shown in Fig. 7.

Fig. 9 Specific torques for the same model as in Fig. 8 at the same times and respective positions.

In Fig. 9 the specific torques are displayed where we use the normalization of D’Angelo & Lubow (2010),${\left(\frac{\mathrm{d\Gamma }}{\mathrm{d}\mathit{m}}\right)}_{\mathrm{0}}\mathrm{=}{\mathit{q}}^{\mathrm{2}}{\mathit{h}}^{-4}{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{P}}^{\mathrm{2}}\mathit{.}$(8)The biggest contribution is generated within a 5H wide region inside and outside the planet. Apart from the initial adjustment of the gap to the moving planet (see Fig. 8) the specific torques reach an equilibrium profile after 100 orbits. As expected from very similar local gap profiles (in Fig. 8), the specific torques do not depend strongly on the position of the planet in the disk. This will be of importance later in explaining the overall migration properties of the planets (see Sect. 5).

As shown in Fig. 10 the gaps, of course, depend on the viscosity. For a fixed α-value, the gap structure near the planet does not change very much, even for different disk masses. Effects farther away from the planet ≳4H are results from the different positions in the disk as well as the movement of the gap through the disk with different velocities (see Fig. 6). At the inner gap edge the surface density is increased for the higher accretion rates and at the outer gap edge for the smaller accretion rate. This is because the inward migration rate of the planet is faster for a denser disk (and thus a higher accretion rate, see Fig. 6) and the gap edges cannot follow as fast as the planet moves. A higher migration rate therefore leads to a stronger gap deformation. The gap depth depends on α in the same manner as for the static planets in Fig. 2.

Fig. 10 Gaps of a q = 0.001 planet with different accretion rates 100 orbits after the release of its orbit. The gap is rescaled according to the disk scale height at the position of the planet, which varies for different values of α and ṁ. The gap profiles are taken from the low-resolution calculations.

Fig. 11 Azimuthally averaged local accretion rate, ṁ(r), at different times during the inward migration of the planet for the standard model (q = 0.001, ṁ = 10-7, α = 0.003). The plot is similar to Fig. 4, which has the same parameters, but now for moving planets. To compare the accretion rate near the planet the position is rescaled to units of the disk scale height at the position of the planet. The labels refer to the elapsed time (in initial orbits) after the release of the planet.

4.3. Flow across the gap

Above we have argued that the fast migration speed of the planet implies that mass has to cross the gap. Now we analyse this process in more detail. In Fig. 11 we display the local mass accretion rate through the disk, ṁ(r), similar to Fig. 4, but now for a migrating planet at different times after the release. At early times, when the planet is moving inwards very fast, the gas is flowing outwards instead of inwards. Even at later times, when the migration rate has slowed down, at the position of the planet the direction of the gas flow is still reversed. This effect is caused by the motion of the planet through the disk, which disturbs the local disk structure and leads locally, at the location of the planet, to a transfer of mass from the inner to the outer disk. The mass transfer is positive because the planet moves faster than the viscous accretion velocity ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}$. As shown above, the gap is bound to move with the planet, but the planet does not drive the disk structure as whole. In consequence, gas from the inner disk, not able to move inwards faster than ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}$, has to move outwards, which means crossing the gap in outward direction. For a planet moving more slowly than the viscous speed (for low disk densities) the radial flow remains always negative.

If the disk were separated and gas could not cross the gap, the inner disk would pile up and the region adjacent to the planet in the outer disk would be depleted because the gas could not follow the planet’s movement. The more massive inner disk would then yield higher positive torques and slow down or even reverse migration. Thus, seeing transport of the gas outwards is a clear sign that the planet is not bound to the viscous accretion velocity and moves independently.

5. Comparing to type II migration

Fig. 12 Migration tracks of planets embedded in disks for different planet masses, for given viscosity and disk accretion rate (α,ṁ) = (0.003,10-7). Coloured lines: theoretical tracks calculated by integrating the measured torques. The dashed black lines are the corresponding observed migration tracks in the simulations.

Fig. 13 Normalized torques acting on a Jupiter mass planet as a function of the local disk mass, with α = 0.003 and q = 0.001. Shown is the evolution of the torques during the migration process, each track starts at the rightmost point with the highest value of ${\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}$. The black dots mark the point when the planet has reached aP = 0.7. During the inward migration the local disk mass ${\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}$, where ΣP is the surface density in the undisturbed surface density profile at planet position ap, decreases. The models differ in their initial disk density Σ0 varied by the disk accretion rate ṁ, where ΣS corresponds to ṁ = 10-7. The torques are averaged over 50 orbits.

After having obtained the migration track of planets in disks with net mass flux we are now in a position to compare these to the type II migration speed of a viscously driven disk.

In any case, the migration of a planet can only be driven by the gravitational torques exerted on it by the disk material. To prove this statement, we plot in Fig. 12 the evolution of the planet’s semi-major as expected from the calculated total torque, Γ, during its motion through the disk (coloured lines). The migration track a(t) is obtained by integrating the formula$\mathit{ȧ}\mathrm{=}\frac{\mathrm{2}\mathrm{\Gamma }}{{\mathit{M}}_{\mathrm{p}}{\mathit{a}}_{\mathrm{P}}{\mathrm{\Omega }}_{\mathrm{K}}}\mathrm{·}$(9)As seen in Fig. 12, the observed migration (black dashed lines) agrees perfectly with the calculated migration tracks (coloured lines) throughout the whole simulation. Hence, the migration is indeed completely determined by the torques from the disk, as physically expected.

Having shown that the torques determine the migration of the planet, we plot in Fig. 13 the normalized torques using the same normalization, Γ₀ as before in Fig. 5. Note that now the normalization factor Γ₀ is a function of the planet’s position, see Eq. (6). In Fig. 13 the torques are displayed during the evolution of the planet. The initial location refers to the rightmost points in the individual curves. They all start at the same height Γ/Γ₀ = 0.18 as given by the fixed planet resultis, see Fig. 5. During the inward migration the value of ${\mathrm{\Sigma }}_{\mathrm{p}}{\mathit{a}}_{\mathrm{p}}^{\mathrm{2}}$ decreases. The black dot in each curve corresponds to the time where the planet has reached a_p = 0.7.

Let us first compare the findings to the linear case, for small mass planets that do not open gaps. 2D simulations yield the following relation for the total Lindblad torque on planets (see Tanaka et al. 2002; Paardekooper et al. 2010)$\frac{{\mathrm{\Gamma }}_{\mathrm{L}}}{{\mathrm{\Gamma }}_{\mathrm{0}}}\mathrm{=}\mathrm{3.2}\mathrm{+}\mathrm{1.468}\hspace{0.17em}\mathit{p,}$(10)where p is the exponent of the surface density profile, with Σ(r) ∝ r^{− p}, and so p = 0.5 in our case. This yields Γ_L/ Γ₀ = 3.934. Expression (10)was derived for much less massive planets in the regime of type I migration and is not directly applicable to our calculations. The idea is to compare the migration speed of massive planets to the linear type I regime. As shown in Fig. 13 the torques are about 20 to 100 times smaller than the linear results, a consequence of the gap. Additionally, the results show a clear dependence on the density of the disk. In the linear case one would expect a constant ratio Γ/Γ₀ for all disk densities.

This fact is not apparent overall, but after the disk and the gap have adapted to the moving planet, the torques become constant for each individual curve, i.e. for each disk density. This implies that the normalized torques for a given planetary mass and α-value are completely defined by the disk density alone. Therefore, we introduce a slow down factor for migrating planets f_mig, which gives the factor by which the torques acting on the moving planet are smaller than those for the same planet on a fixed orbit. Because for gap-opening planets the fixed-orbit torques are again smaller by a factor f_still than the type I torques Γ₀ (see Fig. 5) we obtain the following expression for the torque of a moving planet$\mathrm{\Gamma }\mathrm{=}{\mathit{f}}_{\mathrm{mig}}\hspace{0.17em}{\mathit{f}}_{\mathrm{still}}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathrm{0}}\mathit{.}$(11)In Fig. 14 we show the slow down factor f_mig for various parameters. For small disk densities the slow down factor decreases with increasing disk density. For light planets with M_P<M_J this changes for Σ₀/ Σ_S between 0.1 and 0.5 and the slow down factor strongly increases to values much bigger than one, which marks the onset of type III migration. For the highest disk densities, the slow down factor for the lighter planets are in fact out of the plot range and reach values up to f_mig = 12. In case of smaller viscosity f_mig decreases, but for the same planetary mass the overall structure seems to remain the same.

Fig. 14 Slow down factor fmig depending on the disk density for different planetary masses and viscosities. Here, the surface density is given in units of ΣS, as defined in Fig. 13. The circles are results for α = 0.003 and the squares for α = 0.001.

Fig. 15 Migration rate of a Jupiter mass planet normalized by the viscous accretion velocity. Apart from the different y-axis the plot is similar to Fig. 13. The horizontal line corresponds to ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}\mathrm{=}\mathrm{1}$. The dashed line shows a fit of Duffell et al. (2014), but they used constant viscosity equal to α = 0.01 at r = 5.2 au, so their fit is normalized with ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}$ for that α.

In Fig. 15 we compare the obtained migration speed to the type II migration rate as given by the viscous inflow speed ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}$. In the case of classical type II migration, the curves should be independent from disk mass and identical to the viscous inflow speed, apart from the adjustment in the beginning. However, the results show a clear dependence with disk mass. Obviously, the actual migration of a planet is not coupled to the viscous evolution in the disk, as also pointed out recently by Duffell et al. (2014). Qualitatively, the relationship between disk mass and migration rate is similar to that found by Duffell et al. (2014) found by different means (see dashed curve in Fig. 15). In our case the inward migration is faster than theirs. For light disks with ${\mathrm{\Sigma }}_{\mathrm{D}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}\mathit{<}\mathrm{0.2}{\mathit{M}}_{\mathrm{J}}$ only, we find that migration becomes smaller than the viscous speed. Here, it is important to note that the viscous speed for our constant ṁ models scales directly as ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}\mathrm{\propto }{\mathit{a}}_{\mathrm{p}}\mathrm{\Omega }$ the quantity that Duffell et al. (2014) used in their plots.

From our results we obtained for moving planets that the normalized torque Γ/Γ₀ was constant during the evolution of a planet (see Fig. 13). This implies that the torque reduction factors f_mig and f_still do not depend on the local disk mass, and are only functions of the planet mass and the disk viscosity. Using Eq. (9)for the migration rate and the scaling for Γ₀, we obtain for each planet track$\mathit{ȧ}\mathrm{\propto }{\mathit{a}}_{\mathrm{P}}^{\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}{\mathit{f}}_{\mathrm{slow}}\hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{M}}_{\mathrm{P}}\mathit{,}$(12)where f_slow = f_migf_still. Using the relation ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}\mathrm{\propto }{\mathit{r}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$, as inferred from Eq. (3)for constant H/r, one finds$\frac{\mathit{ȧ}}{{\mathit{f}}_{\mathrm{slow}}{\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}{\mathit{M}}_{\mathrm{P}}}\mathrm{\propto }{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}{\mathrm{\Sigma }}_{\mathrm{P}}\mathit{.}$(13)This relation is plotted in Fig. 16 and indeed the data lie on a straight line, for all the models with different ṁ. The data for an accreting planet fall onto the same line. In case of different viscosity there is a shift due to the change in f_slow, but it follows the same trend.

Fig. 16 Normalized migration rates of planets with different masses between 0.2 and 2 MJ. The migration rate is scaled with the slow down factor and planet mass (see Eq. (13)). The slope of the black line is 1. The dotted line corresponds to α = 0.001. Accreting planets are indistinguishable from the non-accreting planets and therefore not shown here. Otherwise the plot is similar to Fig. 13.

6. Summary and conclusions

We have studied the migration of massive planets in locally isothermal disks with net mass flow, ṁ, through them. In our initial step, keeping the planet at a fixed location, we showed that despite the presence of the planet, the disk can nevertheless transport the full ṁ through the disk and the gap. We analysed the structure of the gap for different planet masses and viscosities. As expected, the depth of the gap increases with planet mass and is reduced with increasing viscosity. Our gap structure and depths for these models are in very good agreement with the empiric formulae presented in Fung et al. (2014). Additionally, the gap opening criterion by Crida et al. (2006) gives a suitable condition if there is a considerable gap depth, ≈0.1 of the ambient density. For large planet mass the normalized torque on the planet strongly decreases (see Fig. 5), in contrast to the linear rate where Γ/Γ₀ has no dependence on planet mass. This is a consequence of the increased gap width with a reduced mass that can drive the planet.

Upon releasing the planet we find rapid inward migration, often faster than the typical viscous inward drift of the disk material. As expected, the planet moves exactly according to the disk torques acting on it. A migration independent of the disk’s drift speed implies that material can cross the gap region during the migration process. We performed simulations without and with mass accretion onto the planet where, in the latter case, the accretion rate was near the maximum possible rate (Kley 1999). As expected, because of the reduced mass in the vicinity of the planet, the migration rate is reduced below the non-accreting case, but only by about 25%. Our main result is summarized in Fig. 15, which shows that the planet migration does not depend on the viscous inflow speed of the disk material. For small disk masses (M_D/M_J< 0.2) only, a Jupiter mass planet moves slower than the disk material. For higher disk masses (given by higher ṁ in our calculations), values and subsequently larger surface densities, Σ, the inward migration becomes faster as well.

An important result of our simulations is the finding that during the migration of the planet the normalized torques remain constant, a feature seen by Duffell et al. (2014). In principle, it might be possible to extract approximate analytical relationships for the slow down factors f_mig and f_still, which would then be useful to calculate theoretical migration tracks to be used, for example, in population synthesis models. To do this, more elaborate parameter studies will be necessary in the future.

Very qualitatively the shape of the curve is similar to that found by Duffell et al. (2014) by a complementary method, but we do not see signs of saturation. Our migration rate is typically faster than theirs up to a factor of 2, but note that their results have been obtained for a spatially constant viscosity, which implies a spatially variable α-value. Because of their special technique of pulling the planet through the disk and measuring the torque after reaching equilibrium, it is not know at which distance of the planet this occurs. Hence, the value of α is not known and the curves cannot directly be compared to each other as the torque depends on the value of α. Edgar (2007, 2008) also used constant viscosity. He also finds that the migration rate depends on the disk mass and the migration timescales are comparable to those we found.

Our findings are also in agreement with those of Nelson et al. (2000) who analysed the migration of a planet in a global disk and constant kinematic viscosity. They inferred a timescale for migration of a massive planet from its starting location at 5 au to 2.5 au of about 2500 orbits, but did not compare this to the viscous accretion rate.

For sufficiently low planetary mass, M_P ≲ M_D, but with M_P still large enough to open a partial gap, the planets can enter the very fast type III migration regime where the migration timescale becomes very short, of the order of about 100 dynamical times. The conditions necessary for type III migration (Masset & Papaloizou 2003) indicate that is may be relevant in the early phases of planet formation with larger disk masses, or in more massive self-gravitating disks.

For low disk masses M_D/M_J< 0.2, the migration rate becomes lower than type II migration. These small disk masses occur only during the end phases of the planet formation process when the accretion rate has already reduced significantly.

In this work, we studied only 2D disks and found that disk material can always cross the gap region as required by the migration speed of the planet. For massive, Jupiter type planets the averaged gap profile is identical for 2D and 3D disks (Kley et al. 2001) and one may expect very similar mass flow rates across the gap region, and hence similar migration accretion rates. The additional assumption of locally isothermal disks is not very restrictive as well because for massive planets that do open significant gaps, the dynamical behaviour is very similar to the isothermal case as has been shown by Bitsch & Kley (2010) for full 3D simulations of radiative and isothermal disks.

Nevertheless, it may be interesting to perform in the future 3D simulations of embedded, massive planets in radiative disks and analyse their migration properties. Those simulations are also important to understand better the mass growth of planets because the inclusion of the 3rd dimension may affect the mass accretion onto the planet, while it only mildly affects the migration.

Our results have consequence for population synthesis models in that a modification of the assumed type II migration speed is required for the models, based on the torques acting on the planets. Possibly, more elaborate models that cover a larger parameter space will allow the construction of suitable fit formulae for the migration of massive planets in the future, but this is beyond the scope of the present work.

Another point to consider is the angular momentum balance of the planet. We have taken out mass from the planet’s Roche lobe, but have not added this to the planet mass. In addition to the mass, the angular momentum of this material, which is gained by the planet, has to be considered. A fraction of this will change the orbital angular momentum, which will influence the migration of the planet. Future simulations need to consider this effect as well.

Acknowledgments

We thank Bertram Bitsch for very fruitful discussions and acknowledge the constructive and helpful comments from the referee to improve this paper. The calculations were performed on systems of bwGRiD, the computer grid of the Baden Württemberg state. The work of Christoph Dürmann is sponsored with a scholarship of the Cusanuswerk.

References

All Tables

All Figures

	Fig. 1 Surface density of a fixed-orbit, non-accreting planet calculation for the standard model with (α,Mp/M∗,ṁ) = (0.003,0.001,10-7) after 2500 orbits. The planet is in the middle of a stable gap and the accumulation of mass near the planet is clearly visible. The spiral arms are damped at the inner and outer boundary so no reflections are seen.
In the text

	Fig. 2 Azimuthally averaged surface density for a non-accreting (solid line) and an accreting (dashed line) planet with q = 0.001, for different values of the viscosity parameter α. The black line indicates the initial density profile at the beginning of the simulation. The red line corresponds to the standard model.
In the text

	Fig. 3 Azimuthally averaged surface density for non-accreting (solid line) and accreting (dashed line) planets and α = 0.003, for different values of the planet mass. The red line corresponds to the standard model.
In the text

	Fig. 4 Azimuthally averaged surface density profile (top) and the migration rate profile (bottom) in case of (α,q,ṁ) = (0.003,0.001,10-7) as in Fig. 1. The planet is at r = 1.0. The profiles are taken at different times after the beginning of the calculation. The black line indicates the initial density profile at the beginning of the simulation before the planet was added.
In the text

	Fig. 5 Normalized torques acting on planet with fixed orbits at the end of the relaxation as a function of the planet mass. The data is taken from the low-resolution calculations. For isothermal disks the normalized torques are independent of the disk density.
In the text

	Fig. 6 Planet migration tracks for different model parameter. Starting from the standard model (in red) with (α,ṁ,q) = (0.003,10-7,0.001), individual parameters have been varied while keeping two others fixed. We varied α (left panel), ṁ (middle panel), and q (right panel). The full and dotted lines correspond to non-accreting and accreting planets, respectively.
In the text

	Fig. 7 Global gap profiles at four different times during migration of the planet for the standard model with q = 0.001 and ṁ = 10-7. The black line indicates the unperturbed density profile without the planet.
In the text

	Fig. 8 Gap profiles for the standard model taken at different positions of the planet in the disk during its migration path. The gap is rescaled using the local disk scale height, H, at the position of the planet. The curves correspond directly to those shown in Fig. 7.
In the text

	Fig. 9 Specific torques for the same model as in Fig. 8 at the same times and respective positions.
In the text

	Fig. 10 Gaps of a q = 0.001 planet with different accretion rates 100 orbits after the release of its orbit. The gap is rescaled according to the disk scale height at the position of the planet, which varies for different values of α and ṁ. The gap profiles are taken from the low-resolution calculations.
In the text

Fig. 11 Azimuthally averaged local accretion rate, ṁ(r), at different times during the inward migration of the planet for the standard model (q = 0.001, ṁ = 10-7, α = 0.003). The plot is similar to Fig. 4, which has the same parameters, but now for moving planets. To compare the accretion rate near the planet the position is rescaled to units of the disk scale height at the position of the planet. The labels refer to the elapsed time (in initial orbits) after the release of the planet.

In the text

	Fig. 12 Migration tracks of planets embedded in disks for different planet masses, for given viscosity and disk accretion rate (α,ṁ) = (0.003,10-7). Coloured lines: theoretical tracks calculated by integrating the measured torques. The dashed black lines are the corresponding observed migration tracks in the simulations.
In the text

Fig. 13 Normalized torques acting on a Jupiter mass planet as a function of the local disk mass, with α = 0.003 and q = 0.001. Shown is the evolution of the torques during the migration process, each track starts at the rightmost point with the highest value of ${\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}$. The black dots mark the point when the planet has reached aP = 0.7. During the inward migration the local disk mass ${\mathrm{\Sigma }}_{\mathrm{P}}{\mathit{a}}_{\mathrm{P}}^{\mathrm{2}}$, where ΣP is the surface density in the undisturbed surface density profile at planet position ap, decreases. The models differ in their initial disk density Σ0 varied by the disk accretion rate ṁ, where ΣS corresponds to ṁ = 10-7. The torques are averaged over 50 orbits.

In the text

	Fig. 14 Slow down factor fmig depending on the disk density for different planetary masses and viscosities. Here, the surface density is given in units of ΣS, as defined in Fig. 13. The circles are results for α = 0.003 and the squares for α = 0.001.
In the text

Fig. 15 Migration rate of a Jupiter mass planet normalized by the viscous accretion velocity. Apart from the different y-axis the plot is similar to Fig. 13. The horizontal line corresponds to ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}\mathrm{=}\mathrm{1}$. The dashed line shows a fit of Duffell et al. (2014), but they used constant viscosity equal to α = 0.01 at r = 5.2 au, so their fit is normalized with ${\mathit{u}}_{\mathit{r}}^{\mathrm{visc}}$ for that α.

In the text

	Fig. 16 Normalized migration rates of planets with different masses between 0.2 and 2 MJ. The migration rate is scaled with the slow down factor and planet mass (see Eq. (13)). The slope of the black line is 1. The dotted line corresponds to α = 0.001. Accreting planets are indistinguishable from the non-accreting planets and therefore not shown here. Otherwise the plot is similar to Fig. 13.
In the text