Toward a new paradigm for Type II migration
^{1}
Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Université Côte d’Azur, Bd de l’Observatoire,
CS 34229,
06304 Nice Cedex 4,
France
email: clement.robert@oca.eu
^{2}
Institut Universitaire de France,
103 Boulevard SaintMichel,
75005
Paris,
France
Received:
31
May
2018
Accepted:
15
July
2018
Context. Giant planets open gaps in their protoplanetary and subsequently suffer socalled type II migration. Schematically, planets are thought to be tightly locked within their surrounding disks, and forced to follow the viscous advection of gas onto the central star. This fundamental principle, however, has recently been questioned as migrating planets were shown to decouple from the gas’ radial drift.
Aims. In this framework, we question whether the traditionally used linear scaling of migration rate of a giant planet with the disk’s viscosity still holds. Additionally, we assess the role of orbitcrossing material as part of the decoupling mechanism.
Methods. We have performed 2D (r, θ) numerical simulations of pointmass planets embedded in locally isothermal αdisks in steadystate accretion, with various values of α. Arbitrary planetary accretion rates were used as a means to diminish or nullify orbitcrossing flows.
Results. We confirm that the migration rate of a gapopening planet is indeed proportional to the disk’s viscosity, but is not equal to the gas drift speed in the unperturbed disk. We show that the role of gapcrossing flows is in fact negligible.
Conclusions. From these observations, we propose a new paradigm for type II migration: a giant planet feels a torque from the disk that promotes its migration, while the gap profile relative to the planet is restored on a viscous timescale, thus limiting the planet migration rate to be proportional to the disk’s viscosity. Hence, in disks with low viscosity in the planet region, type II migration should still be very slow.
Key words: protoplanetary disks – planet–disk interactions – planets and satellites: formation
© ESO 2018
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Planetary migration is a key ingredient to understand the architecture of planetary systems. This radial displacement of planets is due to their gravitational interaction with the protoplanetary disk. These disks surround most young stars, and have a lifetime of a few million years. Planetary migration leads to significant changes in the semimajor axis of all planets (see Baruteau et al. 2014; for a recent review) and carves the structure of planetary systems.
Migration of planets has been extensively studied in recent decades. Small mass planets, for which the response of the disk can be considered linear, do not perturb the density profile of the disk, and are in a regime called type I migration. Giant planets, however, are massive enough to modify the disk radial density profile. They deplete the region around their orbit and create a gap (Lin & Papaloizou 1986a), separating the inner disk from the outer disk. Once a gap is open, the planet is repelled inward by the outer disk and outward by the inner disk. The position of the planet within the gap adjusts so that the torques from the inner and out disks cancel out. However, as the disk spreads viscously and the gas accretes onto the central star, the gap, as well as the embedded planet that carved it, is carried with it. This is the classical scheme of the socalled type II migration (Lin & Papaloizou 1986b), responsible for inward motion of giant planets. In this scheme, the planet does not migrate with respect to the gas, but together with the gas, and acts as a gasproof barrier between the parts of the disk.
This standard scheme of type II migration, where a planet follows exactly the viscous accretion speed of the gas has been questioned by several works. Quillen et al. (2004) note that if the inertia of the planet is much larger than that of the gas originally present in the gap, the disk has a hard time moving the planet, and the migration is slower than the viscous speed. Crida & Morbidelli (2007) add that the corotation torque, exerted on the planet by the gas still present in the gap, may play a role, especially in regions where the background density profile is steep (which promotes a high corotation torque). This could slightly decouple the planet from the gas evolution, but this process relies on a nonempty gap, hence it could be seen as a situation where perfect type II migration is not expected anyway.
Furthermore, Hasegawa & Ida (2013) remark that the assumption that a planet be locked in its gap had no solid physical ground. Lubow & D’Angelo (2006) and Duffell et al. (2014) show that, in simulations, gas is able to cross the gap during planetary migration. Dürmann & Kley (2015; DK15 hereafter) explicitly question the idea that the planet stays in equilibrium in the middle of the gap. They suggest that when the gas reaches an equilibrium gap profile, such that the torques from the viscosity, the pressure, and the planetary gravity balance on each side (Crida et al. 2006), the planet does not necessarily feel a zero torque. Hence, it moves to a different position inside the gap. The motion of the planet then forces the gas to readjust the equilibrium profile, by passing through the planet’s orbit from the inner part to the outer part of the disk. Because this transfer of gas is due to the planet–gas interaction and not due to the disk’s viscous spreading, the evolution of the planet becomes unlocked from the disk’s viscous evolution.
The authors of DK15 find that the migration speed is independent of the disk’s drift speed and mainly depends on parameters such as the mass of the disk or that of the planet. Furthermore, Dürmann & Kley (2017; DK17 hereafter) showed that planetary accretion, insome cases, is able to cut the gas flow across the planet’s orbit. In these cases the gap acts as a barrier between the inner and the outer disks and classical type II migration regime could be reestablished. However, the authors noticed that even in these cases the migration rate can differ from classical type II migration rate. Both studies (DK15 and DK17) considered a classical αviscosity disk and focused on the dependence of migration speed on parameters like the disk mass and the planetary mass for a fixed viscosity value.
These studies show that gas is crossing the gap and therefore the migration speed seems to be independent on the disk’s drift. Although the authors already provided evidence that migration speeds depend on viscosity, a direct comparison of this dependence with the fundamental assumptions of classical type II remains to be conducted. This is precisely the aim of the present paper.
In fact, though it is now admitted that giant planets’ migration does not follow classical type II migration, it is very important to check whether some scaling of migration rate with viscosity is still preserved. Precisely, the existence of giant planets at orbits larger than 1 au in semimajor axis (socalled warm Jupiters) is difficult to explain in the usual paradigm of viscously accreting disks unless considering low viscosity and confirming that migration speed scales with viscosity. Low viscosity is admitted in the central part of the disks, the socalled dead zone, where magnetorotational instability (Balbus & Hawley 1991) does not operate. Moreover, recent studies (Bai & Stone 2013; Bai 2016; Suzuki et al. 2016) have shown that the disk’s structure can be very different from that of a viscously accreting disk. Mass accretion onto the star could be ensured by magnetically driven winds providing angular momentum removal. These disks would have very low viscosity. It is beyond the purpose of this paper to model this kind of disks, however these studies motivate the investigation of the migration of giant planets in low viscosity disks. We will also investigate the role of gapcrossing flows on migration.
The paper is organized as follows. In Sect. 2, we present the physical model and the numerical scheme, while the setup of initial conditions is reported in appendices. In Sect. 3, we study the scaling of migration speed with viscosity by performing numerical simulations of Jupitermass planets in disks with controlled inflow rates and various viscosities, but same gas surface density. In Sect. 4, the accretion of gas by the planets is modeled by removing, with various efficiencies, the gas entering the Hill sphere; it allows us to find the influence of planetary accretion onmigration, and to quantify how important is cutting the gas flow across the gap. Finally, our findings are summarized in Sect. 5, where we propose a new, consistent paradigm to explain giant planet migration.
2 Physical setup
We consider a stationary accreting disk in which a planet is introduced. In this section, we present our background disk model the numerical code used and the prescription for planetary accretion.
2.1 Units and notations
We describe our 2D disk in polar coordinates (r, θ), centered osnto the host star. A subscript “_{0}” denotes values defined at a reference radius r_{0} = 1. Our time unit, hereafter called an orbit, is , where m_{*} is the mass of the central star. The planet mass is defined as m_{p} = qm_{*} = 10^{−3}m_{*}, and its initial semimajor axis as r_{p}(t = 0) = r_{0} = 1.
2.2 Accretion disk model
Because we are interested in comparing the radial drift of the embedded planet to that of the unperturbed disk, the easiest scheme is toset up a steadystate accreting disk. The local accretion rate is defined as (1)
hence, an accreting disk displays v_{r} < 0 and Ṁ > 0, where r is the distance to the central star, Σ is the gas surface density, and v_{r} is the gas radial velocity. We use a nonflared disk scaleheight H(r) = hr, h = 0.05 being the uniform aspect ratio, as well as anαviscosity model (Shakura & Sunyaev 1973), ν = αH^{2}Ω, and a power law density profile, . Within those prescriptions, we obtain the radial velocity for a uniform (hence steady) accretion rate as^{1} (2)
Combining Eqs. (1) and (2) gives , where . Therefore, the disk model is completely defined by fixing α and Ṁ. This leads to (3)
Unless specifically stated, all our simulations share this setup, following Eqs. (2) and (3). A summary of different values used in this study may be found in Table 1. In Sect. 2.7, we exhibit boundary conditions compatible with this physically stable initial state, and in Appendix B, we explain how they affect the disk in the presence of a planet.
Models parameters.
2.3 Hydrothermodynamics
The disk evolves along the Navier–Stokes equations:
where v is the gas velocity, ν is the αviscosity, Φ_{G} is the total gravitational potential yielded by the central star and the planet, and P is the pressure. Assuming a cooling time much shorter than the orbital period , the equation system is closed with a locally isothermal equation of state: (6)
where c_{s} is the sound speed and H(r) = hr is the disk scale height.
2.4 Numerical code
Our experiments were conducted using the 2D hydrodynamic grid code Fargo 2D (Masset 2000). This code solves Eqs. (4) and (5) using a finite difference multistep procedure. The fluid advection step is solved using a Van Leer method (Van Leer 1977). The Fargo algorithm (Masset 2000) is specifically suited for Keplerian rotation where the traditional Courant–Friedrichs–Lewy (CFL) condition (Courant et al. 1928) provides very small timesteps due to fast orbital motion at the inner boundary of the numerical domain. In the Fargo algorithm, the timestep is limited by the perturbed density arising fromdifferential rotation.
Our solver was set with a Courant parameter (Courant et al. 1928) of 0.4. The contributionΦ_{p} of the planet to the gravitational potential is smoothed as (7)
where d is the local distance to the planet, and is the planet’s Hill radius. As our model does not include the disk’s selfgravity, we exclude the material contained in the Hill region of a planet in the computation of the torque acting on it. To this end, we use the tapering function described by (Crida et al. 2008; Eq. (5)), with p = 0.6. This is done to avoid artificial “braking” in the planet’s migration due to the fact that this circumplanetary material does not feel any gravitational torque from the rest of the disk (see Crida et al. 2009).
2.5 Planetary accretion recipe
Accretion onto the planet is handled following the recipe of Kley (1999). At each time step, in every cell located within the Hill radius R_{H} of the planet, a fraction of the gas is removed. This fraction is given by K_{p} f(d) δt, where d is the distance to the planet, δt is the time step, and K_{p} is an arbitrary accretion efficiency parameter. K_{p} is typically ~1, and constrained such that 0 ≤ K_{p} δt < 1. In our simulations δt ≃ 10^{−3}. We use the smooth function f(d) proposed by Crida et al. (2016): (8)
For a uniform density around the planet, this function conserves the planetary accretion rate with respect to the original step function proposed by Kley (1999). It furthermore provides better accuracy when the gas is growing scarce in the planet’s vicinity, without the need for a higher resolution.
This mass taken away from the disk should fall onto the planet, with corresponding momentum. However, here we remove this gas (and its momentum) from the simulation, in order to compare the migration rates of planets of equal, fixed masses. By doing so, we highlight the influence of the gas flow on the migration of the planet. This is reasonable as we are interested in the migration rate of a planet given its mass and not in the planet’s mass accretion rate nor in the final position of an evolving planet.
2.6 Initial conditions for migration
We obtain our initial conditions after two stages, respectively coined introduction and relaxation, lasting T_{intro} = 1000 and T_{relax} = 4000 orbits, respectively. During introduction and relaxation, the planet is held at a constant semimajor axis r_{p} = 1.
During introduction, the planet mass is slowly increased from 0 to qm_{*} over T_{intro}, along a smooth function of time (9)
Classically, a gapopening planet quickly expels the gas from its horseshoe region. This gas tends to accumulate at the gap’s edges then to slowly spread at a viscous rate. To avoid this slow phase, we allow the planet to remove gas instead of scattering it, using the planetary accretion recipe described above. We use a smoothly decreasing planetary accretion efficiency (10)
with in our simulations. As illustrated by Crida & Bitsch (2017), this helps with the gap opening process, and overall saves computational time.
During relaxation, we let the system evolve to a neartosteady state. Convergence is considered achieved when both the mass flow Ṁ (r) and the gravitational torque Γ_{tot} exerted on the planet reach constant values with respect to time. See Appendix A. Hereafter, we call the time origin t = 0, the release date, from which the planet is allowed to migrate.
2.7 Boundary conditions
The choice of boundary conditions is a priori of nonnegligible importance, and can significantly affect the numerical steadystate the disk relaxes into. Our main concern in choosing appropriate boundaries is that the disk should be left unperturbed far from the planet’s orbital radius. In order to achieve this, we make the further distinction between the domain of interest, in other words, the radial vicinity of the planet, and the broad edges of the simulation domain. A smooth transition is used from the planet’s region of influence to the unperturbed initial state at edges of the simulation domain. This is done through a wavekillinglike algorithm (de ValBorro et al. 2006) where at each time step and in dedicated regions^{2}, perturbations in Σ, v_{r}, v_{θ} with respect to the initial state are damped out. This design choice is justified in more extensive details in Appendix B.
2.8 Models parameters
The simulation domain spans over r ∈ [0.1, 3.0]. We have chosen an inner boundary at r = 0.1, although this choice is very expensive in computational time. The advantage is the possibility to study migration of a planet initially located at r = 1 on a large orbital domain, going down in some cases to r = 0.4 with no impact of the inner boundary on migration. An arithmetic radialspacing of grid cells is used with a resolution of (n_{r}, n_{θ}) = (248, 628).
In this paper we aim at testing the scaling of the migration rate with the viscosity and therefore we consider different values of α in the interval [3 × 10^{−3}: 3 × 10^{−4}]. Table 1 describes in code units our models parametrization. Values are chosen so that Σ_{0} =5.3 × 10^{−4} code units in all simulations. This value is low enough that the coorbital mass deficit can not exceed m_{p}, to avoid the runawaytype III migration regime (Masset & Papaloizou 2003), but large enough that the disk is still able to push the planet efficiently^{3}. For m_{*} = m_{⊙} (hence m_{p} = m_{J}), and r_{0} = 5 au, model A1 physically translates into and Ṁ = 2 × 10^{−8} m_{⊙} yr^{−1}.
Figure 1 shows the initial density profiles at the end of the relaxation phase at t = 0 for simulations A1, A2, A5. As expected, gaps are deeper and wider at lower viscosities.
Fig. 1 Azimuthally averaged density profiles at t = 0, when the planet is being released. Black thin solid curve corresponds to the power law profile used at initialization. 

Open with DEXTER 
3 Influence of viscosity on the migration rate
DK15 and DK17 demonstrated that the torque acting on a gapopening planet is primarily dependent on the disk’s mass rather than its accretion rate Ṁ(r) ∝ ν(r). Although they showed that viscosity still affected the torque, it remained to be clarify how viscosity’s role compares to the initial assumption of classical type II. A direct comparison with type II is the topic of this section.
“Classical” type II migration rate is given by , where v_{r} is the viscous speed of the unperturbed disk (Eq. (2)); i.e. the planet is assumed to migrate with the drift rate of the gas. We show in Fig. 2 the migration speed measured in simulations A1 to A5 plotted as a function of the semimajor axis of the planet. In all these runs, no planetary accretion is used. The migration speed was normalized in all simulations to so that scaling of the migration with v_{r} and/or withviscosity is apparent. Because planets do not migrate at identical speeds in different runs, note that a given semimajor axis corresponds to different dates; for instance, A1’s planet reaches r_{p} = 0.8 at t ≃ 1150, to be compared with t ≃ 8000 in case A5. Some observations can be drawn from Fig. 2:
As pointed out by Duffell et al. (2014), because the planet was first artificially maintained on a circular orbit, the corresponding gas distribution at release date is de facto inconsistent with a migrating perturber. Indeed, within this method, we obtain torques inducing migration timescales much shorter than the viscous spreading timescale of gap edges; this is a known effect in type II studies. It follows that the transitional stage immediately following release (r_{p} ≈ 1) ought tobe discarded from our analysis.
Following this transition, normalized migration tracks converge, demonstrating that steady migration rates scale linearly with the viscosity.
Because migration speeds scale with viscosity, convergence is reached on timescales . Therefore, very long integration times are required at the lower viscosities. In fact, properly comparing situations at different values of α requires snapshots with similar r_{p}, i.e. similar t∕τ_{ν} (and similar local disk mass, see below). In particular, the previously mentioned transitional stage corresponds to 0.9 ≲ r_{p} < 1 for every viscosity.
At odds with the classical speed expected for type II migration, the migration rate decreases with the semimajor axis and migration becomes slower than classical type II for r_{p} ≲ 0.54.
This last result is in agreement with DK15 who showed that migration rate decreases with decreasing local disk’s mass , Σ_{p} being the unperturbed surface density at planet position r_{p}. The top horizontal axis of Fig. 2 shows the local disk mass divided by the mass of the planet. The migration speed is equal to v_{r} when this mass ratio is ≃0.2, in remarkable agreement with DK15 (Fig. 15 therein).
Besides these observations, Fig. 2 reveals a puzzling fact for the question we are interested in the migration speed of a gap opening planet is proportional to the viscosity of the disk, but not equal to the radial drift of the gas. The fact that the planet migrates slower than the gas when the disk mass is low is not a surprise, but the reason for a faster migration remains unclear still (although already found by previous studies). To further inquire this possibility, we ran an additional simulation, B1, in a static unperturbed disk, i.e. where Eq. (3) is changed to s = 1 so that v_{r} = Ṁ = 0, Eq. (2), and Σ_{0} is arbitrary, hence kept to A1’s value. The migration speed found is displayed in Fig. 3 as the black curve, where it can be compared to the case s = 1∕2 (orange curve, same data as Fig. 2). Very clearly (and surprisingly), the unperturbed radial velocity of the gas has little influence on the migration speed of the planet.
What pushes the planet is not the radial drift of the gas onto the star. This supports DK15’s claim that the planet inside its gap is not necessarily at equilibrium with the gas, when the gas profile is itself at equilibrium with the planet. Hence, the planet can feel a torque, which drives its migration, even if the unperturbed disk has no radial drift. This is already a change of paradigm for type II migration. Why this torque should be proportional to the viscosity, however, is unclear at this point. This suggests that the picture may be more complicated than the one just described, as we will see later.
For the planet to migrate faster than the radial velocity of the unperturbed gas, one would naturally expect that it migrates with respectto the gas. In this picture, gas should cross the planet’s orbit from the inner to the outer disk to sustain migration. The following section is dedicated to studying the influence of planetary accretion on such mass exchanges and its resulting effects on migration.
Fig. 2 Normalized migration speed as a function of the semimajor axis for sets A1 to A5 with K_{p} = 0. The normalized local disk mass is indicated as a secondary graduation for the xaxis. 

Open with DEXTER 
Fig. 3 Migration speed as a function of the semimajor axis in two simulations with α = 0.003 and same surface density. Orange curve: case A1 (s = 1∕2); same curve as Fig. 2 (not normalized here). Black curve: case B1. 

Open with DEXTER 
4 Effect(s) of planetary accretion
Here we test how introducing planetary accretion into our model affects the gas flow through the planet’s orbit, and how the migration rate is changed in turn. Accretion’s efficiency is parameterized by the dimensionless number K_{p}, which we vary as K_{p} ∈{0, 0.2, 1.0, 5.0}, using run A2 as the reference case K_{p} = 0. The initial state in all runs is identical to that of A2. Unlike the timedependent accretion efficiency previously described, here we do not follow Eq. (10), and K_{p} is instantly switched to a nonzero value at t = 0, when the planet is released.
In orderto preserve comparability, we emphasize the importance of keeping the global structure of the disk as unperturbed as can be, despite the fact we are adding a sink point to the hydrodynamical model. Our boundary conditions ensure that the disk’s profile stays unchanged far away from the planet. Hence, we except all changes due to K_{p} to stay local tothe planet’s vicinity.
4.1 Material exchange between reservoirs
The inner and the outer disks constitute our two gas reservoirs of interest. In order to keep track of the planet’s relative radial displacement with respect to those reservoirs, we use passive, dimensionless, scalar tracers η_{i∕o}. Let us give a nonambiguous and partly arbitrary definition of their respective initial distributions. We should be concerned about avoiding confusion with material originated from the HorseShoe Region (HSR), whose width is given by Masset et al. (2006) as w_{HSR} ≈ 2.5 R_{H}. Hence, a precautionary choice is to consider only material initially distant of at least 3R_{H} from the planet’s orbit. Within this definition, we define η_{i,o}(r, θ, t = 0) = r∕r_{0}, for r ≷ r_{p} ± 2.5R_{H}, and η_{i∕o} = 0 elsewhere. Let us acknowledge that, as a tracer is advected along with the gas, its value in a given cell becomes the mass weighted average of the tracers that are found in the considered cell at the considered time. Therefore, at t > 0, we expect to find mixed material, displaying values of the tracer that do not correspond to the initial value but rather a weighted average of it.
Figures 4 and 5 display, for simulation sets A1 and A2 respectively, the evolution of those tracers, sampling over our parameter K_{p} from 0.0 to 5.0. Once more, following Duffell et al. (2014) and DK15, we find that, in the nonaccreting case, some gas is effectively transported from the inner disk to the outer disk as migration proceeds. Not only does the planet migrate with respect to the medium, it also actively ejects some material to larger orbits. This still holds in the weakly accreting case K_{p} = 0.2, although the transport efficiency is being slightly decreased. However, it is not so in the “standard” and “strong” accretion cases (K_{p} = 1, K_{p} = 5), where the outward flow from the inner disk is so efficiently blocked that the outer disk rushes into the planet’s vicinity.
Let us observe that while we successfully introduced accretion as a means to prevent gapcrossing flows, the procedure also profoundly modified the nature of the flow and added complexity to the picture. Indeed, not only did we prevent inner disk material to transfer into the outer disk, we also allowed the outer disk material to reach the planet’s feeding zone, hence causing depletion to happen in both halves of the disk. Additionally, we note that despite the origin of gas at a given radius being widely different depending on the accretion rate (see Fig. 5), at lower viscosity values (A2) the disk’s structure stays almost selfsimilar whatever the planetary accretion efficiency K_{p}, see Fig. 7. The net effect on migration is unclear at this point and will be discussed in Sect. 4.2. Furthermore, the transition from a naturally occurring gapcrossing flow to a twoway accretion flow with increasing K_{p} appears to be smooth. Indeed, we see in the weakly accreting case that, as less material flows from the inner disk to the outer disk, the latter immediately reacts and follows the planet more closely than in the nonaccreting scenario. From a quick comparison between Figs. 4 and 5’s respective “weak rate” panels, one can see that the stronger the viscosity, the quicker this compensation mechanism (normalizing time to the viscous timescale). This is confirmed in stronger accretion cases K_{p} = 1.0, 5.0. Hence, at least for high values of α, planetary accretion cannot cut the otherwise existing outward^{4} gapcrossing flow without causing a more rapid inward flow of the outer disk.
In short, we find that an accreting planet does not prevent gas from entering its HSR, even though it can nullify gapcrossing flows. Let us now discuss the implications of accretion on the migration speed.
Fig. 4 Evolution of passive tracers η_{i∕o} for simulation set A1, following material originated in the outer (top panels)/inner (bottom panels) disk, in polar coordinates. Left panel: initial state (white is 0). All but leftmost column: 1000 orbits after release, for varying values of K_{p}. We emphasize that color scales are different across rows. For the sake of readability, the disk is here displayed with angular coordinates such that θ_{p} = π in every frame. Filled lines show the planet’s feeding zone, 0.8R_{H} in radius; dashed lines are 3R_{H} large in radius, encompassing a somewhat broader region than the typical HSR. 

Open with DEXTER 
Fig. 5 Same as Fig. 4 for simulation set A2. Snapshots are taken after a α_{1} ∕α_{2} longer period. 

Open with DEXTER 
Fig. 6 Radial position of the planet r_{p} VS time, forvarying accretion rates. Panels correspond to simulation sets A1 (left panel) and A2 (right panel). The plot includes analytical migration tracks Eq. (11) with varying time offsets as guidelines. 

Open with DEXTER 
Fig. 7 Density profiles snapshots at fixed time with varying accretion efficiency for simulation sets A1 and A2. Black thin lines indicate the initial states (same as Fig. 1). 

Open with DEXTER 
4.2 Impact of accretion on migration
It is clear from our previous observations that the density of the outer disk, hence the negative torque it yields on the planet, should be reduced by planetary accretion. DK17 already noted that planetary accretion can reduce migration speeds. Nonetheless, we stress that the positive contribution from the inner disk can also diminish as the disk is being forcedly depleted. Here we measure the effective migration rates against accretion efficiency and give further interpretation.
Because planetary accretion stops the gapcrossing gas flow, one may expect it to consolidate the classicaltype II migration scheme. We recall that in this scheme, the migration rate is given by , where v_{r} is the viscous speed of the unperturbed disk Eq. (2); i.e. the planet is assumed to migrate with the drift rate of the gas. Integrating Eq. (2), we get the analytical evolution of r_{p} (t) for a planet migrating in a classical type II fashion: (11)
Time evolution of r_{p} is shown in Fig. 6 for various accretion efficiencies (K_{p}), against this theoretical track, for simulation sets A1 and A2 (1 panel per value of α). We observe that planetary accretion has significant impact only in the highest viscosity case α = 0.003, where the migration rate can be reduced below the theoretical rate (cases K_{p} ≥ 1) after a few 1000 orbits in agreement with findings from DK17. However, changes in migration are barely noticeable in the second case where α = 0.001. We have checked that for lower viscosities (simulations A3, A4, and A5), the migration speed is hardly impacted by planetary accretion as well. Increasing the accretion efficiency of the planet (thus effectively cutting the flow of gas across the gap) does not allow the planet to decrease its migration speed down to v_{r} if α ≤ 0.001.
The role of accretion in the viscous case can be understood with corresponding density profiles Fig. 7. For the A1 runs, the local disk’s structure is being efficiently affected by the accretion, decreasing both inner and outer disk’s densities, hence their respective torques densities (as shown in Fig. 8). In contrast, the planet can hardly make a dent in the density profile in the lower viscosity cases. This is consistent with that a lower viscosity implies both a wider and deeper gap: as far less material is available in the vicinity of the separatrices of the HSR, the influence of planetary accretion on the dynamical evolution of the disk is hindered.
Except for high values of α, we can conclude from our study that planetary accretion is not of significant importance for the migration of giant planets. Actually, the fact that the curves in Fig. 2 overlap while there is no planetary accretion suggests that the latter is not important in setting the proportionality between the migration rate and the disk’s viscosity. The reason why planetary accretion plays a role in the migration speed in the viscous A1 case is not because it cuts the flow, but because it perturbs strongly the density profile and broadens the gap.
Fig. 8 Radial torque densities for simulation sets A1 and A2. 

Open with DEXTER 
5 Conclusion
To summarize, we have found in Sect. 3 that although type II migration speed is proportional to the gas’ viscosity, it is not driven by the radial inward drift of the gas. In particular, we confirmed that the giant planet can actually migrate faster than the gas drifts, and even in a stationary disk. We concluded that what drives type II migration isthe imbalance between the torques felt by the planet from the inner and outer disks, as pointed out by DK15. However, the width and shape of the gap is not directly linked to the viscosity, especially at low ν where the pressure effects are dominant (Crida et al. 2006); hence, we do not expect this torque imbalance to be proportional to ν, in contrast with the observed migration speed. In Sect. 4, we have seen that gapcrossing flows are actually negligible at low viscosity, and that cutting this small gas flow with planetary accretion hardly impacts the migration speed. Thus, the planet migrates faster than the disk drift, even when no gas is exchanged between the inner and outer disks. Gapcrossing flows cannot be responsible for the observed fast migration, in contrast with the of case type III migration (Masset & Papaloizou 2003).
These two results allow us to draw a new, consistent picture of type II migration. As a giant planet forms, it opens a gap by perturbing the gas profile with the gravitational torque it exerts. The gas reaches a new equilibrium profile on each side of the gap. Nonetheless, the planet inside its gap feels a nonzero torque, because the inner and the outer torques have no reason to balance out (as recently studied by Kanagawa et al. 2018). Thus, the planet has to migrate inward. As it does so, some gas may cross the gap from the separatrix of the HSR, although this is not enough to restore the initial gap profile in the frame of the planet if the viscosity is low ad the gap is wide (regardless of whether the planet accretes or not). Therefore, the density distribution has to adapt to the new position of the planet^{5}, and this is done over a viscous time. Once the gas is again at equilibrium with the planet, the planet is not in an equilibrium inside the gap anymore, and we are back to the initial situation.
In this scheme, the planet may well migrate faster than the gas drifts, because it is pushed by a torque that has no connection with the drift of an unperturbed disk. But because the gapcrossing flow is negligible, the planet must migrate at a rate proportional to viscosity, otherwise it would pile gas up in the inner disk and leave a depleted outer disk behind, eventually halting its migration.
Although the final result (migration speed of gapopening planets is proportional to the viscosity) is in line with the standard picture of type II migration, this new scheme is conceptually revolutionary in our understanding of this phenomenon, and allows us to reconcile all the puzzling observations that have been made recently, questioning the standard picture. Additionally, it confirms that even if some gas may cross the gap, a giant planet in a lowviscosity disk should migrate slowly. In this frame, the abundance of warm Jupiters, who did not migrate all the way toward their star, may suggest that most protoplanetary disks have a low effective viscosity in the planetforming region.
Acknowledgements
We acknowledge support by the French ANR, project number ANR13BS05000301, projet MOJO (Modeling the Origin of JOvian planets). C.M.T.R., H.M., and A.C. acknowledge funding from ANR grant No. ANR16CE310013 (PlanetFormingDisks). HPC resources from GENCI [IDRIS] (Grant 2017, [i2017047233]) and from “Mesocentre SIGAMM”, hosted by Observatoire de la Côte d’Azur were used. C.M.T.R. and E.L. wish to thank Alain Miniussi for his valuable help in maintaining and developing the code base used in this work. Figures in this paper were produced with matplotlib (Hunter 2007). The authors wish to thank the anonymous referee for their valuable observations and help in making this paper clearer.
Appendix A Disk relaxation: equilibrium of torques on a fixedorbit planet
Here we give details on the relaxation of the disk toward initial conditions for planetary migration. In Fig. A.1, the total gravitational torque Γ_{tot} that the disk exerts on the planet is plotted alongside with its inner and outer components.
Fig. A.1 Evolution of the gravitational torque Γ_{tot} acting on the planet during introduction and relaxation ( here denotes time during those stages). Γ_{i} is the contribution of the inner disk (up to r = r_{p} = 1), Γ_{o} of the outer disk, and Γ_{tot} = Γ_{i} + Γ_{o}. A vertical line indicates the end of the planet’s introduction and starting date of the disk’s relaxation. 

Open with DEXTER 
While the onesided torques decrease in absolute value as the gap opens, the total torque remains nearly constant. We measure its final value as only a fraction (~0.12) of the differential Lindblad torque: (A.1)
where Σ_{p} is the unperturbed surface density at planetary semimajor axis, r_{p}. This is expected for a gapopening planet. For lack of an existing normalization factor for type II migration, let us derive one here. The viscous torque exerted by inner material r < R onto its complementary part r > R in an unperturbed disk may be expressed as (A.2)
Hence, the expected value the total torque in classical type II can be estimated as the differential viscous torque that a gas annulus occupying the gap would feel as (A.3)
where w_{HSR} is the width of the horseshoe region, estimated as 2.5 R_{H}. We observe that Γ_{tot} converges toward ~4 Γ_{II}. This discrepancy has been tackled by Dürmann & Kley (2015) who showed that in fact, the final value of this torque is primarily determined by the disk mass has little to do with the disk accretion rate Ṁ, or equivalently its viscous torque T_{ν} ∝ ν ∝Ṁ. Furthermore, they showed that eventhough the inner and outer disks reach stationary density profiles, the planet is not in an equilibrium position between them. Indeed, a stationarystate disk embedding a fixedorbit planet does not yield the expected
type II torque but rather a significantly stronger one.
Appendix B Design of boundary conditions
Here we provide more insight on how the boundary conditions were selected and acknowledge the consequences of this choice. An obvious choice for boundaries would consist in simply fixing the values of velocities and surface density to their initial state. In principle, this design would allow for conservation of both the total mass in the simulation domain and the initial flow Ṁ throughout the relaxation. However, this does not hold within the staggered grid scheme used in Fargo. Indeed, radial velocities are defined at lower edge of cells while density is centerdefined. While this is a convenient design for numerical integration of hydrodynamics equations, it makes it impossible to keep the density, velocities and the net mass flux entering the simulation domain constant all at once. This is because the net mass entering the simulation at the outer edge is effectively defined as a combination of fixed and free numbers. Figure B.1 illustrates this point, and define the two equally problematic possible designs. A solution to this issue would be to correct the imposed value for Σ_{n} (case 1) or v_{r,n+1} (case 2) at each time step to achieve the desired value in Ṁ. In practice this proved less than efficient, as Ṁ(r) is very sensitive to even the most subtle spatial or time variation, convergence was not achieved in practical times.
Fig. B.1 Two possible designs for fixed boundary conditions in a 1D grid where radial velocities are defined at inner cell edge (red arrows) and density is centerdefined (blue circles). Boldface is used to indicate fixed quantities. The net mass inflow is determined as a combination of quantities inside of a green dashed ellipse, which always contain both fixed and free values. The last cell harbouring free quantities is indexed by n. 

Open with DEXTER 
Although a steadystate is characterized by uniformity in accretion rate Ṁ (r), there is in fact no strong argument in favor of enforcing its initial value at the grid’s edges. Indeed, DK15 already clarified that the migration rate is not bound to the unperturbed Ṁ. Furthermore, as previously observed, enforcing a constant mass flow through both radial edges would cause the total mass, dictated solely by s and the radial limits of the grid, (r_{min}, r_{max}), to be constant. In this case, one would obtain qualitatively different results only by enlarging the simulation domain, corresponding to different physical scenarios. For instance, strong pressure bumps may be found at the edges of the gap for a narrow enough grid, corresponding to a fast planetformation scenario where a jovian mass would agglomerate over the course of a few 100 orbits.
For these reasons, we chose to relax the assumption that the final uniform value in Ṁ should be exactly equal to that of the initial, analytical state described in Sect. 2.2, and instead favor the fact that the disk should stay close to unperturbed far from the planet. As explained in Sect. 2.7, we extend the notion of “boundaries” to broad radial domains where perturbations with respect to the initial state are consistently damped out. This artificial process is applied to velocities (both radial and azimuthal) as well as to surface density. A direct consequence of this is that mass is no longer being conserved in wavekilling regions, so that desired uniformity in Ṁ (r) is only relevant outside those regions. Additionally, this damping process prevents reflections of sound waves, such as spiral wakes, as is its original purpose (de ValBorro et al. 2006).
Figure B.2 shows the evolution of the radial flow profile obtained in the case A1. Although not perfectly converged, at , we settled for what we considered an acceptable level of uniformity to save computational time.
Fig. B.2 Relaxation of the azimuthaly averaged radial flow, case A1. denotes time during introduction and relaxation stages (). Limits of the domain of interest (excluding wavekilling zones) are shown as dashed lines. Ṁ_{init} denotes the uniform value in the analytic state of the unperturbed disk described in Sect. 2. 628 snapshots covering oneorbit were averaged to obtain each line in this figure. 

Open with DEXTER 
References
 Bai, X.N. 2016, ApJ, 821, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 769, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Baruteau, C., Crida, A., Paardekooper, S.J., et al. 2014, Protostars and Planets VI, 667 [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Ann., 100, 32 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Crida, A., & Bitsch, B. 2017, Icarus, 285, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., & Morbidelli, A. 2007, MNRAS, 377, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Sándor, Z., & Kley, W. 2008, A&A, 483, 325 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crida, A., Baruteau, C., Kley, W., & Masset, F. 2009, A&A, 502, 679 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crida, A., Bitsch, B., & Raibaldi, A. 2016, in SF2A2016: Proc. of the Annual Meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, J. Richard, L. Cambrésy, et al. 473 [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Duffell, P. C., Haiman, Z., MacFadyen, A. I., D’Orazio, D. J., & Farris, B. D. 2014, ApJ, 792, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Dürmann, C., & Kley, W. 2015, A&A, 574, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dürmann, C., & Kley, W. 2017, A&A, 598, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hasegawa, Y., & Ida, S. 2013, ApJ, 774, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W. 1999, MNRAS, 303, 696 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1986a, ApJ, 307, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1986b, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H., & D’Angelo, G. 2006, ApJ, 641, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Masset, F. S., & Papaloizou, J. C. B. 2003, ApJ, 588, 494 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. S., D’Angelo, G., & Kley, W. 2006, ApJ, 652, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Quillen, A. C., Blackman, E. G., Frank, A., & Varnière, P. 2004, ApJ, 612, L137 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van Leer B. 1977, J. Comput. Phys., 23, 276 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Azimuthally averaged density profiles at t = 0, when the planet is being released. Black thin solid curve corresponds to the power law profile used at initialization. 

Open with DEXTER  
In the text 
Fig. 2 Normalized migration speed as a function of the semimajor axis for sets A1 to A5 with K_{p} = 0. The normalized local disk mass is indicated as a secondary graduation for the xaxis. 

Open with DEXTER  
In the text 
Fig. 3 Migration speed as a function of the semimajor axis in two simulations with α = 0.003 and same surface density. Orange curve: case A1 (s = 1∕2); same curve as Fig. 2 (not normalized here). Black curve: case B1. 

Open with DEXTER  
In the text 
Fig. 4 Evolution of passive tracers η_{i∕o} for simulation set A1, following material originated in the outer (top panels)/inner (bottom panels) disk, in polar coordinates. Left panel: initial state (white is 0). All but leftmost column: 1000 orbits after release, for varying values of K_{p}. We emphasize that color scales are different across rows. For the sake of readability, the disk is here displayed with angular coordinates such that θ_{p} = π in every frame. Filled lines show the planet’s feeding zone, 0.8R_{H} in radius; dashed lines are 3R_{H} large in radius, encompassing a somewhat broader region than the typical HSR. 

Open with DEXTER  
In the text 
Fig. 5 Same as Fig. 4 for simulation set A2. Snapshots are taken after a α_{1} ∕α_{2} longer period. 

Open with DEXTER  
In the text 
Fig. 6 Radial position of the planet r_{p} VS time, forvarying accretion rates. Panels correspond to simulation sets A1 (left panel) and A2 (right panel). The plot includes analytical migration tracks Eq. (11) with varying time offsets as guidelines. 

Open with DEXTER  
In the text 
Fig. 7 Density profiles snapshots at fixed time with varying accretion efficiency for simulation sets A1 and A2. Black thin lines indicate the initial states (same as Fig. 1). 

Open with DEXTER  
In the text 
Fig. 8 Radial torque densities for simulation sets A1 and A2. 

Open with DEXTER  
In the text 
Fig. A.1 Evolution of the gravitational torque Γ_{tot} acting on the planet during introduction and relaxation ( here denotes time during those stages). Γ_{i} is the contribution of the inner disk (up to r = r_{p} = 1), Γ_{o} of the outer disk, and Γ_{tot} = Γ_{i} + Γ_{o}. A vertical line indicates the end of the planet’s introduction and starting date of the disk’s relaxation. 

Open with DEXTER  
In the text 
Fig. B.1 Two possible designs for fixed boundary conditions in a 1D grid where radial velocities are defined at inner cell edge (red arrows) and density is centerdefined (blue circles). Boldface is used to indicate fixed quantities. The net mass inflow is determined as a combination of quantities inside of a green dashed ellipse, which always contain both fixed and free values. The last cell harbouring free quantities is indexed by n. 

Open with DEXTER  
In the text 
Fig. B.2 Relaxation of the azimuthaly averaged radial flow, case A1. denotes time during introduction and relaxation stages (). Limits of the domain of interest (excluding wavekilling zones) are shown as dashed lines. Ṁ_{init} denotes the uniform value in the analytic state of the unperturbed disk described in Sect. 2. 628 snapshots covering oneorbit were averaged to obtain each line in this figure. 

Open with DEXTER  
In the text 