© E. Lega et al. 2022

1 Introduction

Raw observations of extra-solar planets show two distinct populations of giant planets. The first is composed of objects orbiting very close to their host star (known as ‘hot-Jupiters’) and the second is characterised by planets with an orbital period larger than 100 days, the so-called warm and cold Jupiters.

From a theoretical standpoint, the observations of giant planets at distances larger than 0.5 au from their host stars (or periods larger than 100 days) are difficult to understand. Giant planets are expected to have started to form at the snowline, the place where the condensation of water creates an over-density of solid material and the rapid formation of the first planetesimals (Ida & Guillot 2006; Schoonenberg & Ormel 2017; Drazkowska & Alibert 2017; Drazkowska & Dullemond 2018). The snowline in young discs is expected to be at 5–7 au (Lee et al. 2021; van’t Hoff et al. 2020), thus the current location of the cold Jupiters implies that they migrated just a few au in the protoplanetary disc, at odds with classical theoretical expectations.

The migration of a planet changes regime as the planet grows in mass. When the planet has a mass smaller than 10–20 Earth masses, it migrates in the Type I regime. Then, as the planet grows more massive and starts opening a gap in the disc, the migration evolves towards the Type II regime, which is valid when a deep gap encapsulates the giant planet’s orbit.

This paper is part of a project focussing on the Type II regime. Type II migration is an important phase since Jupiter mass planets are expected to form in young discs on timescales of about one million years and this migration mode can therefore become the dominant one when acting over expected disc lifetimes of a few Myr (Johansen et al. 2019). For discs with traditional values of the viscosity parameter α¹ a giant planet needs to start Type-II migration at least 10–20 AU to remain beyond  1 AU from the parent star after a My of migration (Coleman & Nelson 2014; Bitsch et al. 2015; Voelkel et al. 2020; Robert et al. 2018). Our goal is to understand under which conditions Type-II migration can become slow enough to allow giant planets to start at the snowline (5–7 au) without overshooting the warm-Jupiter region after a ~ Myr of migration.

The classic view of the Type II regime is that migration has a speed comparable to the viscous radial motion of the gas, because the planet and gap have to migrate in concert (Ward 1997). This process has been revisited recently (Dürmann & Kley 2015; Robert et al. 2018). Results have confirmed that the migration rate is proportional to the disc’s viscosity for values of the α viscosity parameter larger than 10⁻⁴, though it isnot exactly equal to the unperturbed drift speed of the disc. For these values of α, Type II migration should reduce by a substantial fraction the orbital radius of a planet within the disc’s lifetime (Lega et al. 2021).

While it was originally thought that viscosity in discs might arise from turbulence generated by the magneto-rotational instability (MRI; Balbus & Hawley 1991), leading to α values largerthan 10⁻⁴, it was later realised that the ionisation of the gas near the midplane of the disc is too weak to sustain the MRI Gammie 1996; Stone et al. 1996. Even more recently, the inclusion of non-ideal MHD effects, such as ambipolar diffusion, led to the conclusion that the coupling between the magnetic field and the gas should not make the disc turbulent even at its surface (see Turner et al. 2014 for a review). Thus, discs are probably much less viscous than previously thought. Among the alternative mechanisms for generating turbulence, in purely hydrodynamic discs, the most promising is the vertical shear instability (VSI; Nelson et al. 2013; Stoll & Kley 2014). However, the VSI should not be active at the disc midplane within a few astronomical units from the star because the cooling rates are too slow; a recent study in discs with increasing physical realism (purely hydrodynamics) provides values of α with an upper limit of 10⁻⁴ within 5 au (Ziampras, in prep.). These results show the relevance of studying giant planet migration in very low viscosity discs.

In a recent paper (Lega et al. 2021), we have considered disc models with very low viscosities where α has values in the interval [0 : 10⁻⁴]. We have shown that for α ≤ 10⁻⁵ migration differs from classical Type-II migration in the sense that it is not proportional to the disc viscosity. In discs with a vanishingly small viscosity, gap formation by a giant planet leads to the formation of a vortex at the outer edge of the gap. The evolution of the system then depends on the importance of self-gravity in the disc (Zhu & Baruteau 2016). By comparing simulations with and without self-gravity, Lega et al. (2021) identified a region of parameter space, mainly characterised by disc mass and orbital distance, where self-gravity can be neglected. Here, the planet remains on an almost circular orbit and inward migration occurs only as long as the disc can refill the gap left behind by the migrating planet, either due to gas diffusion caused by the presence of the vortex or because of the inward migration of the vortex itself due to its interaction with the disc. We have called this type of migration ‘vortex-driven migration’. This migration is very slow and, in addition, cannot continue indefinitely because eventually the vortex dissolves. Thus, this result potentially provides a promising way to explain the limited migration that most giant planets seem to have experienced.

The Lega et al. (2021) study, as well as the present one, is limited to disc masses and orbital distances from the star for which self-gravity can be neglected. Further studies, including self-gravity, will be necessary to study the migration of planets forming farther out in the disc, at distances where self-gravity cannot be neglected.

However, the discs considered in Lega et al. (2021), due to their small viscosities, have a negligible radial mass transport, incompatible with the mass accretion rates that are typically observed for young stars (~ 10⁻⁸ M_⊙ yr⁻¹, with a large, order of magnitude scatter around this value (Hartmann et al. 1998; Manara et al. 2016). The aim of the present paper is to test if a similar migration pattern holds in more realistic discs with radial transport of gas towards the star consistent with typical accretion rates of young stars. Magnetically driven disc winds have been proposed as a mechanism to remove angular momentum from thin ionised surface layers of low viscosity protoplanetary discs (Suzuki & Inutsuka 2009; Bai & Stone 2013; Turner et al. 2014; Gressel et al. 2015; Béthune et al. 2017) promoting in these layers a fast radial transport of gas towards the central star. However, using three dimensional non-ideal magneto-hydrodynamical simulations for planet migration studies would be prohibitive from the point of view of computational time; therefore in the present paper we use simple hydrodynamical simulations in which we mimic the effects of magnetically driven winds. We generate radial transport ofgas by applying a synthetic torque in a layer of finite thickness located either at the disc’s surface or close to its midplane. The thickness of the layer is characterised by its column density, which is denoted by Σ_A (indicating this is the “active layer”). A similar model has been introduced by McNally et al. (2020), who showed that the wind-driven accretion flow does not modify the migration in the specific case of low mass planets.

In this work, we extend that study to giant planet migration and explore the dependence of the results on the parameter Σ_A, looking for conditions that allow us to recover the slow vortex-driven migration mode unveiled in Lega et al. (2021). Additionally, we also consider the case where the torque due to the magnetic field induces a radial flow near the midplane. This is expected in the region of a disc within 10 au from the central star if the Hall effect is large and the magnetic field is aligned with the disc rotation axis (Bai & Stone 2013; Lesur et al. 2014; Béthune et al. 2017), or beyond 10 au if Ohmic resistivity is weak (Lesur 2021). The case of radial flow induced near the midplane has been studied in the case of small mass planets embedded in two dimensional inviscid discs (McNally et al. 2017). The authors have shown that the dynamical corotation torque can slow down planet migration or even reverse it such that the planet runs away outwards. This case has also been recently studied with a two dimensional model for the case of Saturn mass planets (Kimmig et al. 2020) showing that, in the case of very strong wind, the planet may undergo rapid type III outward migration. Clearly it is important to extend these studies to giant planet migration in three dimensional discs.

The paper is structured as follows: in Sect. 2 we describe our physical model; in Sect. 3 we explain the modifications made to the code to generate radial gas transport and we provide the setup of simulations in Sect. 4. The migration of Jupiter-mass planets is described in Sect. 5; the case of a migrating Saturn massplanet is considered in Sect. 6. Conclusions and discussion are provided in Sect. 7.

2 Physical model

The protoplanetary disc is treated as a three dimensional non self-gravitating gas whose motion is described by the Navier–Stokes equations. We use spherical coordinates (r, φ, θ) where r is the radial distance from the star, that is from the origin, φ is the azimuthal coordinate starting from the x-axis and θ the polar angle measured from the z-axis (the colatitude). The midplane of the disc is at the equator . We work in a coordinate system which rotates with angular velocity:

where M_⋆ is the mass of the central star, G the gravitational constant and r_p(0) is the initial distance to the star of a planet of mass m_p, assumed to be on a circular orbit. The gravitational influence of the planet on the disc is modelled as in Kley et al. (2009) using a cubic-potential of the form: (1)

where d is the distance from the disc element to the planet and ϵ is the softening length; our nominal value is ϵ = 0.4R_H where R_H is the Hill radius: . The function f is given by: f(x) = x⁴ − 2x³ + 2x.

To the usual Navier-Stokes equations (see for example Lega et al. 2014) we add an evolution equation for the internal energy per unit volume e = ρc_vT (where ρ and T are respectively the volume density and the temperature of the disc and c_v is the specific heat at constant volume): (2)

Here τ_c is the cooling time (τ_c = 2π∕Ω, and Ω is the local orbital frequency); T₀ is the initial temperature, defined as , with h₀ being the disc aspect ratio and μ the mean molecular weight (μ = 2.3 g mol⁻¹ for a standard solar mixture).

The system of equations is closed defining the pressure as: p = (γ − 1)e with γ = 1.4 the adiabatic index.

In short, we use an adiabatic EoS, with exponential damping of the temperature perturbations to the initial temperature profile. We use Eq. (2) instead of the simpler locally isothermal EoS or the pure adiabatic model in order to avoid disc instabilities like the vertical shear instability (Nelson et al. 2013).

3 Radial transport of gas

3.1 Generating an accretion flow with an imposed torque: layered discs

Instead of performing MHD simulations, which are prohibitive from the point of view of computer time, we mimic the removal of angular momentum from the surface layers of the disc by a magnetised wind, or the torque due to a horizontal magnetic field acting near the midplane, using a prescribed torque.

We remark that actual MHD winds would transport angular momentum towards the surface of the disc and eventually eject the gas, whereas in purely hydrodynamical models we can only consider the lower portion of the wind where the magnetic torque causes accretion. Moreover, disc winds in MHD simulations have a radially varying ‘accretion efficiency’, so that the radial mass accretion rate can vary in steady state because the vertical outflow rate would remove the excess mass. In the present model we do not have outflows, and therefore in order to set up a steady accretion disc we generate radial transport of gas that corresponds to a radially constant mass accretion rate, in a layer with vertically integrated density:

The underlying idea is that this layer has a high enough ionisation fraction to be strongly coupled to the magnetic field, such that angular momentum removal is efficient. We describe this as the “active layer”. Estimates of the ionisation rate as a function of Σ_A can be found in Armitage (2019) and references therein. We treat the column density of the active layer as a free parameter and, according to Armitage (2019) (see their Fig. 8), we test Σ_A = 0.1, 1 and 10 g cm⁻² as representative values of the column density of active layers with different ionisation rates. The mass flux towards the star in the active layer is denoted by Ṁ_A. The radial velocity of the gas depends on Ṁ_A through the relation: (3)

The accretion flow is obtained by applying a synthetic torque, γ_A, to the corresponding layers. Considering that the specific angular momentum is , the specific torque γ_A required to generate the radial velocity v_r is: (4)

In order to generate this torque we impose an acceleration in the azimuthal direction (φ) (5)

The acceleration in the azimuthal direction generates a radial velocity in the active layer of magnitude given by Eq. (3). We do not modify the other hydrodynamical variables; a mass flow will appear in the active layer as a consequence of the radial velocity generated by applying the synthetic torque.

In order to have a smooth transition between the active layer and the rest of the disc where no torque is applied, we multiply γ_A by a function with the Fermi function: (6)

The function g(z) is the ratio between the disc column density and Σ_A: g(z) = Σ(z)∕Σ_A; b is a parameter (b = 5 in our simulations) regulating the steepness of the transition of from 1 to 0. We remark that the Fermi filter satisfies: .

In order to simulate radial accretion in the vicinity of the midplane we apply a synthetic torque to the column of gas with z ≤ H = h₀r, while no torque is applied above. The transition between the torqued and untorqued regions is again implemented via the Fermi function with g(z) = z∕H.

3.2 Radial transport in classical viscous discs

We compare migration in a disc with an active layer to that in a classical viscous α-disc, and here we briefly review radial transport in a disc with the α viscosity prescription (Shakura & Sunyaev 1973). The viscosity ν is given by:

with c_s(r) the sound speed defined as and α is a constant. Angular momentum conservation in a steady-state disc links the radial transport (in the following Ṁ_D) to the disc surface density (Σ(r)) and the viscosity through: (7)

We assume α = 10⁻⁵, because this value gave rise to vortex-driven migration in the simulations of Lega et al. (2021) while making the simulations less sensitive to numerical resolution than those with α = 0. We note that we also assume α = 10⁻⁵ in the layered disc models.

4 Simulation setup

Our three dimensional simulations are done with the code FARGOCA (FARGO with Colatitude Added; Lega et al. 2014)². The code is based on the FARGO code (Masset 2000) extended to three dimensions. The fluid equations are solved using a second order upwind scheme with a time-explicit-implicit multi-step procedure. The code is parallelised using a hybrid combination of MPI between the nodes and OpenMP on shared memory multi-core processors. The code units are G = M_* = 1, and the unit of distance r₁ = 1 is arbitrary when expressed in au. The unit of time is therefore . For the simulations in this paper we adopt the Sun-Jupiter distance as the unit of length in au: r₁ = 5.2. When presenting simulation results distances are expressed in au and time in years. The mass of the planet is normalised with respect to the mass of the star and indicated byq ≡ m_p∕M_*.

We consider discs of initial aspect ratio h₀ = 0.05 and radial domain extending from r_min ≤ r ≤ r_max with r_min = 1 and r_max = 46 (au) as in Lega et al. (2021). In the vertical direction the discs have an opening angle of 12°, covering approximately 4.2 hydrostatic pressure scale heights from the midplane.

We do not consider inclined planets and therefore, using the symmetry of the disc, we can simulate only the half-disc above the mid-plane (multiplying the resulting force by a factor of two to obtain the total force exerted on the planet). Mirror boundary conditions are applied at the midplane as in Kley et al. (2009) and reflecting boundaries are applied at the disc surface. In the radial direction we use the classical prescription (de Val-Borro et al. 2006) of evanescent boundary conditions.

Discs have a viscosity α = 10⁻⁵ and a surface density with Σ₀ = 6.76 × 10⁻⁴ in code units (222 g cm⁻² at 5.2 au).

The exponent −1∕2 in the surface-density radial profile ensures that the viscosity-driven mass flow Ṁ_D is independentof r, given that the disc is not flared. The above values of viscosity and disc surface density provide a mass transport: Ṁ_D = 8 × 10⁻¹¹M_⊙ yr⁻¹ (from Eq. (7)).

The subscript D indicates that our discs are “dead” in the sense that a negligible amount of gas is transported due to viscosity.

The temperature profile is the same for all the discs: isothermal in the vertical direction with the midplane temperature scaling as 1∕r. The temperature is damped towards T₀(r) ∝ 1∕r according to Eq. (2).

Long run simulations in three dimensions are challenging in term of computational time, therefore we have used a moderate resolution of (N_r, N_θ, N_φ) = (576, 24, 360) corresponding to about six gridcells per scale height and four gridcells per planet’s Hill radius in the radial and azimuthal directions and height gridcells per planet’s Hill radius in the vertical direction. Convergence of results with respect to resolution has been shown for the same set of parameters (in classical viscous discs) in Lega et al. (2021).

In Table 1, we report the main simulations parameters. In the case of imposed radial flow in the vicinity of the midplane we consider also a control simulation with radial flow directed outwards: simulation set H_out. Although this simulation was mainly run to check the symmetry with respect to the H_in case, we remark that in some cases magnetised discs have both accreting and decreting regions (Béthune et al. 2017, Fig. 22); thus the planet can indeed experience a local radial outflow of the disc.

4.1 Setup of a disc with radial transport in thin surface layers

We show inthis section the new disc setup corresponding to simulation set S01. We assume that the mass flow carried in the active layer is Ṁ_A = 10⁻⁸ M_⊙ yr⁻¹ (each disc hemisphere carries Ṁ_A = 5 × 10⁻⁹ M_⊙ yr⁻¹), as this is the typical accretion rate observed in young stars (Hartmann et al. 1998; Manara et al. 2016).

Figure 1, top panel, shows in (r, θ) coordinates the radial velocity profile of the initial disc with no embedded planet. We clearly see that gas is transported through a thin surface layer. In Fig. 1, bottom panel, we show the total mass flow in M_⊙ yr⁻¹ transported across the disc. We notice that the transport in the active layers corresponds to 75% of the imposedflow and that the remaining part of the disc does not provide any radial transport (or negligible with respect to Ṁ_A as expected). We notice that most of the gas radial flow comes from a very thin layer covered by few grid-cells near the transition of from 1 to 0, where the product ρ(z)v_r has a sharp maximum because of the Gaussian shape of ρ(z). Hence, the integrated mass flux is underestimated unless the transition is finely resolved. Precisely, to reach 95% of the imposed flow in the active layers, a vertical resolution of 96 grid cells is required. We stress that, although the flow is underestimated, this has no impact on the radial velocity which is the key parameter in our study.

In order to run simulations over 10⁵ years we keep our modest vertical resolution. The discretisation effect on the flow transported to the star is the same for the three simulations sets with active layers(we have 80% of the nominal flow for S10) and this makes the three cases comparable. For the simulations with transport in the vicinity of the midplane the gas carried by the active layers corresponds to the expected value with no dependence on resolution.

We remark that, although we are considering a purely hydrodynamic model, we capture the essential result of a disc with basically no transport near the midplane and fast accretion at its surface. We notice for example that the accretion layer in Béthune et al. (2017) is radially localised while we impose here accretion through the whole disc and that the fastest inwards radial velocity in pure MHD simulations is close to c_s while it is about 10% of c_s in our model S01 and are respectively 1% and 0.1% of c_s in models S1 and S10.

Fig. 1 Disc with radial transport ṀA = 10−8 M⊙ yr−1 generated in the active layer with Σ(z) ≤ 0.1 g cm−2. Top panel: radial velocity normalised with respect to the sound speed; bottom panel: total mass flow as a function of the distance from the star. The flow transported by the active layer is 75% of the nominal flow MA = 10−8 M⊙ yr−1 (see text); the flow due to viscosity is clearly negligible with respect to ṀA, so that the total flow is transported by the active layer and the green and red curve are superposed (the minus sign in the plot indicates inflow).

4.2 Introduction of the planet in the simulation

We introduce the planet on a fixed orbit and grow its mass up to its final Jupiter-mass value q_J (or Saturn mass q_S, see Sect. 6) in T_growth = 2π × 600 according to: (8)

The simulation is then continued for an additional 200 orbits to stabilise the disc, before the planet is released in the migration runs at t = 9481 yr⁻¹. In this way we avoid the excitation of instabilities that would arise if the planet were initialised with its final mass (see also Hammer et al. 2017; Hallam & Paardekooper 2020).

We have reported in Fig. 2 the vertically integrated surface density at t = 9481 yr⁻¹ for simulation S01. The white line delimits the region above which the synthetic torque is applied. We remark that, as described in Sect. 3, the synthetic torque is applied everywhere in the disc on a column of density Σ_A, which means that in the gap region (in the vicinity of 5.2 au) the torque can be applied down to the midplane if the gas density in the gap is smaller than Σ_A (see Sect. 5.1).

One may question whether it is realistic that the gas is ionised at all scale heights in the gap, which is the condition to be active. The answer depends on the incidence of the ionising radiation: whereas irradiation from the central star is almost tangent to the disc and should not penetrate directly into the gap down to the midplane, scattered x-rays and cosmic rays can reach the midplane if the column density is thin enough (Kim & Turner 2020). Our model is certainly simplistic in this respect, given that we impose a thickness of the active layer without modelling self-consistently the value of this thickness, but it is worth noting that Kim & Turner (2020) report that gas in the gap for a Jupiter or Saturn mass planet is expected to be sufficiently ionised to couple strongly to magnetic fields.

In the present work, we also consider that gap formation does not affect the magnetic flux distribution. Few papers have investigated magnetic fields inside gaps (Carballido & Hyde 2017; Zhu et al. 2013) and this challenging and important issue certainly deserves future dedicated studies.

Fig. 2 Vertically integrated surface density, Σ(r, z) = ∫, in simulation set S01 at time t = 9481 yr before the planet is released in the migration phase. The white dashed line corresponds to the ΣA = 0.1 g cm−2. The synthetic torque is applied above this line up to the disc surface. We remark that the formation of the gap implies that the synthetic torque is applied deep down in the disc.

5 Migration of a Jupiter mass planet

We investigate the migration of the planet located initially at r = 5.2 au for the six simulations listed in Table 1.

The value of the Toomre Q parameter in a disc with initial density Σ₀ = 6.76 × 10⁻⁴ at distance r = 5.2 au from the star is much larger than 1 so that the effects of self-gravity on global disc evolution can be neglected (Lega et al. 2021).

We recall that a giant planet opening a gap in the disc can migrate only at the speed at which the disc can readjust itself in order for planet and gap to migrate together (Ward 1997). In viscous discs, the gas can displace radially on the viscous timescale, which sets the planet migration speed to be proportional to ν∕r. In low viscosity discs the growth of a giant planet leads to the formation of a vortex at the pressure bump located at the outer edge of the gap. The vortex clearly appears in the leftmost panels of Fig. 3 with the same qualitative characteristics for the six simulations sets. The vortex migrates inwards and spreads radially until it disappears (rightmost panels of Fig. 3). Vortex migration is due to its ability to exchange angular momentum with the disc, via the spiral density waves that it generates (Paardekooper et al. 2010). Vortex migration is inwards because the vortex is located at the outer edge of the gap, so that the disc inside the vortex’s orbit is strongly depleted and the outer disc dominates the interaction with the vortex. We have shown in Lega et al. (2021) that the vortex transfers angular momentum to the disc beyond its orbit, carving a secondary gap (second and third columns of Fig. 3). The secondary gap is also clearly shown in Fig. 4. Once this gap is formed, the vortex slows down, spreads radially and eventually dissipates completely, leaving the outer disc partially depleted near the planet’s orbit. It is also worth noticing in Fig. 4 that the gap depth and shape, in this first migration phase, is the same for all the simulation sets, an indication that the vortex dominates the dynamics despite the existence of an imposed advection.

We report in Fig. 5 the evolution of semi-major axes and eccentricities for the simulations with transport in surface and midplane layers, and for reference we also plot the case of a classical viscous disc (simulation set “NW”). An initial phase of rapid inward migration is observed for all simulations (up to t ~ 40 000 yr). Eccentricity is slightly excited in this phase and strongly damped later (Fig. 5, bottom panel) for all the simulations with the exception of simulation H_out³.

We have shown in Lega et al. (2021) that a planet’s migration is sustained because the migration of the vortex refills the gap left behind by the migrating planet. We have called this process ‘vortex-driven’ migration.

In this phase, the synthetic torque generating radial gas flow in the disc’s surface layer or midplane layer (up to one disc scale height) appears to have a negligible or moderate effect on the planet’s migration. We remark that the inward flow near the midplane (simulation H_in) slightly accelerates the planet’s inward migration while the outward flow of simulation H_out brakes the vortex-driven inward migration and contributes to an early reversal of the direction of migration. When the vortex has completely dissipated, the planet’s migration appears to depend on the properties of the layer that transports gas towards the star (Fig. 5, top panel). For Σ_A = 0.1 g cm⁻², we observe a short phase of outward migration (70 000 < t < 80 000 yr) and finally the planet slowly migrates inwards. The behaviour of this second phase is qualitatively similar to that observed for a planet embedded in a classical, non-layered disc with negligible radial flow (the cyan curve in Fig. 5, top panel). From Fig. 6, however, we see that for t > 80 000 yr the planet in simulation S01 is migratinginwards somewhat faster than in the classical disc case (simulation set NW). In the case of thicker active surface layers, more rapid inward migration is observed (simulations sets S1 and S10 in Figs. 5 and 6). In the cases H_in and H_out, the migration after the dissipation of the vortex is respectively inwards and outwards, like the imposed direction of the gas flow near the midplane (Figs. 5 and 6). We explain below the origin of these differences.

5.1 Radial gas flow through the planet’s gap

In order to understand the dependence of the migration speed on the properties of the active layers, we have computed the net radial gas flow near the region close to the gap (red curves in Fig. 7). The gas flow has important fluctuations over time; thus, to compute the net flow, we performed a temporal average over 1000 orbits in a time interval centred at about t = 75 000 yr. In this timeinterval the planet is very slowly migrating outwards for simulations NW, S01 and H_out, and is slowly migrating inwards for the other cases (see Fig. 5). The migration range is indicated in Fig. 7: a black and a red filled circles indicate respectively the planet–star distance at the beginning and at the end of the time interval considered for this average. The short migration range justifies the averaging procedure. For the stratified disc models the contributionfrom the active layers (green curve) as well as the one from the non-active or dead part of the disc (blue curve) are also plotted.

In the classical non-layered disc (Fig. 7, top left panel) we show the importance of flow fluctuations by reporting also the flow’s time average on five shorter windows of 200 orbits each, with time intervals given by: T_i = [(71 146 + (i − 1) × 200 × 11.85) : (71 146 + i × 200 × 11.85)] yr, with i = 1…5. In the outer disc (r > 6 au) the flow appears to fluctuate significantly over time; the average over 1000 orbits (red curve) appears effective in removing these fluctuations, remaining close to zero except near the 3 : 1 mean motion resonance.

In the inner disc the flow appears to be stationary and positive at the 1 : 2 mean motion resonance. From the observed gas flows, we expect a slowly pile-up of gas at these resonances (the location of the 1:2, 3:2, 2:1 resonances at t = 71 146 yr are indicated by the dotted vertical lines in Fig. 7). Scaling from our experience with more viscous discs, such resonant perturbations of the flow would require prohibitive long integration times to be averaged out. What is important to notice in Fig. 7, top left panel, is that there is no transport of gas across the gap as expected in low viscosity discs (see for example Robert et al. 2018). The shape of the gap at the beginning and at the end of the averaging time interval can be appreciated in Fig. 8.

In the case of a very thin active layer (Fig. 7, top right panel) we observe similar fluctuations of the radial flow at the resonant locations in the dead region of the disc (blue curve); instead the radial flow in the active layer (green curve) nicely corresponds to the flow generated by imposing the synthetic torque (Eq. (9)), which means that the advection of gas in the active layer is barely perturbed by the presence of the planet and its resonances. It is worth noticing from Fig. 8 that, for the simulation S01, the disc’s surface density at the bottom of the gap is of the order of Σ_A = 0.1 g cm⁻²; this means that all the gas in the deepest part of the gap is considered to be ‘active’. This confirms that the existence of the gap is not an obstacle to the flow of the active gas.

In the middle left panel of Fig. 7 (simulation S1), we clearly see that the radial flow transported by the active layers through the gap region is less than half the imposed value Ṁ_A = 10⁻⁸. This means that the planet acts as a partial barrier for the gas transported by the active layer; therefore gas transported by the active layer piles up at the outer gap’s edge. Remember that the planet can migrate only if the disc is able to refill the empty region left behind by the migrating planet; otherwise the planet “detaches” itself from the outer disc and, feeling a progressively weaker negative torque from the outer disc, its migration has to stop eventually. The pile up at the outer edge of the gap of the gas delivered by the active layer is the process that refills the disc and allows the planet to migrate.

Increasing further the thickness of the layer (Fig. 7, middle right panel, simulation set S10) we observe that there is no gas flow through the planet’s gap, which means that the planet is able to block the full flow of gas in the active layer and therefore the pile-up effect is even more pronounced. Consequently, inward migration is faster than in the case S1. Consistently, we also observe in Fig. 8 that the density of gas at the bottom of the gap is smaller in the S10 simulation than in the S1 and S01 simulations. In fact, the gas in the gap – which is active – is in steady state and its density is Σ_A f_gap where f_gap is the fraction of the flow of the active layer which is able to penetrate into the gap. The complement of the flow piles up at the outer edge of the gap, at a rate (1 − f_gap)Ṁ_A, driving the planet’s migration. We stress the fact that it is the gas pile up at the outer gap’s edge that sustains migration by readjusting the disc so that the planet and the gap can migrate together. Instead, the gas flowing through the gap has no impact on the migration; this flow has an impact only on the gap’s depth.

Precisely, we observe that the gap depth stabilises and is very similar in all the simulations with f_gap close to zero (middle and bottom panels of Fig. 8) including the reference case without a wind-driven accretion flow (NW) and the two simulations with inflow and outflow close to the midplane (H_in and H_out). As we observe in the bottom panels of Fig. 7 for the two simulations with radial flow imposed in a layer of thickness H at the midplane, the planet appears to be able to block the flow in the gap region. In the case of radial outflow we observe a large effect due to the 1:2 mean motion resonance. Associated with this effect, the surface density profile shows a kink around 3.5 au (Fig. 4). There is also a visible bump at 2 au, a relic of the gap-opening phase (Fig. 3) which, in this case, has not been eroded by inward migration. Both features in the surface density profile are responsible for the positive torque felt by the planet (and hence its outward migration), while the outward radial flow refills the gap left behind by the migrating plane, sustaining the planet migration. This positive torque contribution decreases when the planet migrates further out. Finally, as previously observed, the migration speed becomes similar in magnitude (and opposite in sign) to the case H_in (Fig. 6). In Fig. 7 bottom left panel, we observe that the contribution of the active layer close to the 1:2 resonance is positive, although the applied synthetic torque provides an inward radial flow as can be appreciated at distances larger than 8 au from the star. We investigate the vertical structure of the flow in the vicinity of the 2:1 resonance in the Appendix A.

Fig. 3 Contours of surface density (multiplied by r1∕2 to flatten the radial profile) for simulation sets NW, S01, S1, S10, Hin, Hout (from top tobottom) with a migrating Jupiter-mass planet. The panels correspond to different times, reported on top of each panel.

Fig. 4 Azimuthally averaged radial surface density profiles at two specific times (reported on top of the panels) within the rapid inward vortex-driven migration phase. No clear differences appear among the six simulations sets.

Fig. 5 Evolution of the semi-major axes and eccentricities for simulations in layered discs with different column densities of the active layers (ΣA). For comparison, data for a classical viscous disc are also plotted (label ‘No Wind’).

Fig. 6 Migration speed vs. time for the six simulation sets. A first phase of inward vortex-driven migration is common to all the cases. Atabout t = 50 000 yr the migration paths bifurcate: for the very thin active layer of simulation S01 the final speed of migration is very similar to the NW case; increasing the depth of the layer a regime of faster inward migration is observed. When the active layer is at the midplane the planet migrates inwards (Hin) and outwards (Hout) at opposite speeds in the final phase of the simulation.

Fig. 7 Total mass flow transported through the disc (red curve). A temporal average over 1000 outputs (spaced by Δt = 11.85 yr) has been performed in order to smooth fluctuations. For the non-stratified case (top left panel) we have split the averaging interval into five parts. We show the five averages to illustrate the time variations of the flow (see text). In the averaging time interval the planet is migrating outwards for simulations NW, S01 and Hout and is migrating inwards for the cases with thicker active layers: S1 and S10 and for simulation set Hin. The black and red filled circles at Ṁ = 0 indicate the planet–star distance respectively at the beginning and at the end of the averaging interval. For the stratified disc models the contribution from the active layers (green curve) as well as the one fromthe dead region of the disc (blue curve) are also plotted. Notice that for simulations Hin and Hout we apply thesynthetic torque in the region between the midplane and z = H (active layers); while for z > H no torque is applied (dead layers). The grey horizontal line at Ṁ = ∓ 10−8 M⊙ yr−1 indicates the expected flow transported by the active layers (the sign plus is for simulation set Hout). The vertical dotted grey lines provide the location of the main mean motion resonances computed using the planet position at t = 71 146 yr. The evanescent boundary condition is applied in the domain inside the vertical black dotted line at 1.25 au.

Fig. 8 Azimuthally averaged radial surface density profiles at specific times (reported on top of the panels). At the time of top and middle panels, the migration behaviour has clearly bifurcated with slow outward migration for discs NW and S01 and rapidinward migration for the simulation sets with thicker active layers as well as for the case Hin. Rapid outward migration is observed for the case Hout. Despite the differences in the migration regimes, the gap depth slowly evolves from the time corresponding to the middle panel to theone represented in the bottom panel: gaps are deeper in all the simulations sets in which the planet blocks the flow from the active layers.

5.2 Halting the accretion flow and planet migration rates

In this subsection, we consider the conditions under which the torque from the planet can halt the wind-driven accretion flow, and how the planet migration rate is expected to scale with the wind-driven accretion rate through the disc.

There are two limiting situations to consider. The first corresponds to a very rapid accretion flow in the active layer that can traverse the planet’s orbital location unimpeded by the torque from the planet. The active layer has little-to-no effect on the planet’s migration rate, which is therefore expected to be close to zero, as in the vortex-driven migration described in Lega et al. (2021). This case occurs only if Σ_A is sufficiently small because, for a given stellar accretion rate, the radial speed of the flow is inversely proportional to Σ_A. The other limit corresponds to the planet being able to completely halt the accretion flow, such that no gas crosses the planet’s orbit. In this case, as the planet migrates inwards and moves away from the original outer edge of the gap, the accretion flow fills in the gap left behind by the planet and migration proceeds on the time scale over which this occurs.

First, we compute when the transition should occur between the two limiting cases of unimpeded flow past the planet and no flow past the planet. The torque exerted by the planet on the disc (in our case the active layer) from a distance r > r_p to infinity can be found in Crida et al. (2006) and references therein: (9)

where Δ_r = r − r_p and q is the ratio between the mass of the planet and the mass of the central star. We consider in the following the torque at a distance from the planet equal to the Hill radius, meaning: . With this choice, the formula simplifies to: (10)

We now consider the specific torque exerted by the planet: and compare it to the specific torque γ_A (Eq. (4)) that is imposed in the active layer to generate the radial flow and associated accretion rate (Ṁ_A = 10⁻⁸ M_⊙ yr⁻¹ in this study). As we see from Eq. (4), γ_A is inversely proportional to Σ_A. Figure 9 shows γ_p and γ_A as a function of Σ_A, and demonstrates that the synthetic specific torque exerted on the active layer overcomes the specific torque received from the planet only in the case of the thinnest layer (simulation set S01) for our assumed accretion rate. This is consistent with the fact that the gas passes through the planet’s orbit at a rate that is unaffected by the planet as shown in Fig. 7, top right panel.

For the cases S10 and H_in, H_out (these last two cases have a column density of ~ 100 g cm⁻²), however, the torque exerted by the planet overcomes that imposed on the active layer, the accretion flow is completed halted by the planet and gas has to pile up at the outer edge of the gap, thus sustaining migration. The case S1 shows significant reduction in the flow past the planet, but not complete blockage of the accretion flow, and so represents an intermediate case between S01 and S10.

Equating the expression for γ_p with that for γ_A from Eq. (4) gives the column density for the active layer at the transition between impeded and unimpeded accretion flow: (11)

The previous equation can be rewritten as: (12)

The equivalence between the two specific torques gives the threshold between the two migration regimes: for γ_p < γ_A the gas should pass through the planet’s orbit without significant pile-up, and the planet is expected to experience only short range migration as in the vortex-driven migration regime. In the opposite limit, γ_p > γ_A, the flow of gas is inhibited by the presence of the planet and gas piles-up outside of the planet’s orbit. The inward migration speed, ṙ_p, is then controlled by the rate at which the accretion flow fills in the gap as the planet migrates inwards. Migration is always the outcome of the imbalance of the torques felt by the planet from the inner and the outer parts of the disc. Because the inner disc preserves its unperturbed density up to the gap’s edge, inward migration can be sustained only if the pile-up of gas outside of the gap restores the original unperturbed density Σ. This leads to the equation: (13)

which gives the migration rate (14)

For a planet orbiting at 5.2 au in a disc with Σ = 222 g cm⁻² and Ṁ_A = 10⁻⁸ M_⊙ yr⁻¹, Eq. (14) gives ṙ_p = 12.3 au yr⁻¹, which is in decent agreement with the migration rate observed in Fig. 6 for run S10 with Σ_A = 10 g cm⁻². We expect the migration rate to be slower when the accretion flow is only partially blocked because the rate at which the gap refills is decreased in this situation. This is supported by the fact that run S1 produces a slower migration rate than S10.

The picture presented above neglects second order effects that are likely to be important. One of these is that the migration speed is connected to the asymmetry of the gap associated with the whole disc, and not just to the ability of the planet to block the accretion flow in the active layer. Figure 8 (middle and bottom panels) shows the bottom of the gap is asymmetric with respect to the planet position in the cases S1 and S10, while it is more symmetric in simulation H_in because in this case γ_p ≫ γ_A. An effect may also arise from the fact that the accreting layer is near the disc surface in the run S10 whereas it is located near the midplane in run H_in. These factors likely explain why the final phase of planet inward migration is slower in the case H_in despite the fact that in runs S10 and H_in the same net flow is imposed on the disc and 100% of it is blocked by the torque from the planet.

We also notice that the column density for the active layer at the transition between impeded and unimpeded accretion flow scales with the planet–star distance as (Eq. (12)). As a consequence, a planet formed at about 5 au in a disc with ionised column density Σ_A slightly larger than at 5 au (or equivalently γ_p > γ_A), may experience fast inward migration until it reaches a distance to the star such that and the flow of gas is no more inhibited by the presence of the planet. The planet is then expected to experience short range migration and possibly stall in the warm-Jupiter region. As an example, in the case of an ionised column density of Σ_A = 0.5 g cm⁻², and Ṁ_A = 10⁻⁸ M_⊙ yr⁻¹, a Jupiter mass planet is expected to enter in the regime of short range migration at au.

Finally, in this section, we comment on how interactions between gap opening planets with classical viscous discs differ compared to interactions with laminar discs that have layered accretion. In the model of a layered disc we have considered, the torque per unit mass acting on the accreting material is set by the mass accretion rate and the column density of the active layer, and if these do not change then the torque has a constant and well defined value. A planet in this case can block the accretion flow by exerting the required torque on the active layer, as discussed above. In the case of a viscous disc, however, the formation of a gap generates radial gradients that increase the rate of diffusion of gas into the gap. This is illustrated by the equation for the radial velocity in a viscous disc (15)

which shows the radial velocity increases in magnitude in response to gradients developing in the flow. Thus, in the viscous case a gap is carved until an equilibrium is eventually reached between the torque exerted by the planet on the disc and the viscous torque (setting the gap’s radial profile); instead, in the case we consider in this paper no equilibrium between γ_p and γ_A occurs and the gas either piles up far enough from the planet’s orbit or passes through the planet’s orbit basically unimpeded. We recall that we consider the wind torque to be steady with a prescribed efficiency, neglecting possible inhomogeneities arising from a radial transport of magnetic flux in the gap’s vicinity. Therefore, in our model the radial flow due to the wind torque is not associated with an increase of α viscosity in the gap. A radial flux of angular momentum, inducing a radial flow like in an alpha viscous model, requires the modelling of the magnetic field and is expected to depend on its geometry (Lesur et al. 2013).

One consequence of the fact that the gas either piles-up or passes quasi unimpeded through the planet’s orbit is that we might expect the accretion of gas onto giant planets to occur at different rates, depending on whether or not they are embedded in viscous or layered discs, even when the global mass accretion rate through the disc is the same. This issue will be explored in a forthcoming paper (Nelson et al., in prep.).

Fig. 9 Specific torque γA and specific planetary torque exerted by a Jupiter mass planet on the outer disc (Eq. (10), see text) γp as a functionof ΣA. We notice that γA > γp only for the model with the thinnest layer (S01). The dot at ΣA = 100 g cm−2 refers to the case where the torque is applied to the layer around the midplane (Hin and Hout).

6 Migration of a Saturn mass planet

From the analysis of the migration of a Jupiter mass planet we have concluded that the equivalence between the two specific torques gives the threshold layer depth that separates the regime of short range migration from the one of fast inward migration. To confirm our findings, in this section we study the migration of a Saturn mass planet in discs with active layers that have different column densities.

We test the validity of Eq. (11) by simulating the migration of a Saturn-mass planet in discs with column layers of depth Σ_A = 0.3 g cm⁻² and Σ_A = 1 g cm⁻². These two values bracket the value of g cm⁻² obtained from Eq. (11) for a Saturn mass planet and Ṁ_A = 10⁻⁸ M_⊙ yr⁻¹. Additionally, we run a reference simulation for the disc without a wind-driven accretion layer. We use the same prescription as in Eq. (8) for the growth of the planet to its final mass value. We report in Table 2 the names and the main simulation parameters.

It is useful to start from the reference case NW_S for which we observe a phase of rapid inward migration (Fig. 10) followed by a phase of slower outward migration after vortex dissipation. The mechanism is the same as described for the Jupiter planet case, however the range of inward migration is slightly larger (for comparison we also plotted in Fig. 10 the data corresponding to simulation NW). By reducing themass of the planet the strength of the vortex appears to weaken. A weaker vortex opens a less deep secondary gap and the balance between the torques from the inner and the outer disc occurs later in time, when the planet is slightly closer to the star with respect to the more massive planet case. The planet migrates outwards in the final phase of the run with speed slightly reducing towards the end of the simulation.

The migration of the planet in the two layered discs confirms the expected behaviour: in both cases we observe “vortex-driven” migration and the migration paths bifurcate after vortex dissipation: in the case of the thin layer (simulation set S03_S) the planet’s migration slows down and practically stops towards the end of the simulation, while for the thicker layer, after a very short phase of outward migration, the planet migrates inwards again. No eccentricity excitation is observed in the three cases (not shown here).

We confirm therefore that the migration behaviour after vortex dissipation depends on the magnitude of the specific torque exerted by the planet at the outer edge of the gap with respect to the torque exerted by the wind.

Fig. 10 Evolution with time of semi-major axes for simulations in layered discs with different depths of the active layer ΣA) for the case of a Saturn mass planet. For comparison we also plotted (cyan dashed line) the data for simulation NW.

7 Conclusions and discussion

The migration of giant planets in classical low viscosity discs has been recently studied by Lega et al. (2021). Interestingly, in discs with vanishing viscosity migration is driven by the migration of the vortex that forms at the outer edge of the gap. Under some conditions this vortex-driven migration can be very slow and then reverses, until it eventually stops. This result is very promising to explain the limited migration that most giant planets seem to have experienced, as is the case for the ‘cold-Jupiter’ population. However, classical low viscosity discs do not account for the gas accretion rates observed onto young stars.

In this paper, we have studied giant planet migration in low viscosity discs in which we have modelled the radial transport of gas in laminar accretion layers. This radial transport mimics the effect of angular momentum loss due to magnetised discs winds, or due to horizontal fields near the midplane being generated by the Hall effect, and account for the typically observed gas accretion rates of Ṁ = 10⁻⁸ M_⊙ yr⁻¹. We have modelled the effect of the disc wind by applying a synthetic torque in a surface layer of the disc characterised by a prescribed column density Σ_A. We have also considered a case with accretion focussed near the disc midplane to mimic transport properties of disc outer regions characterised by weak Ohmic diffusion or due to the Hall effect.

Migration appears to be characterised by two phases. In the first phase, planet migration is driven by the vortex and is directed inwards. This phase ends when the vortex disappears after having opened a secondary gap, as typically observed in vortex-driven migration (Lega et al. 2021). In this phase, the speed of migration is independent of the thickness of the active layer. In the second phase, migration is observed to depend on the thickness of the active layer.

By comparing the torque exerted by the planet to the prescribed torque applied in the active layers we have shown that gas advection modifies planet migration when the torque exerted by the planet overcomes the torque applied in the active layers. In this case, the gas transported by the active layer piles up at the outer gap’s edge. Considering that giant planet migration is possible only if the disc is able to refill the empty region left behind by the migrating planet, the gas pile up at the outer edge of the gap is the process that refills the gap.

On the contrary, when the torque applied on the active layers overcomes the torque exerted by the planet the gas transported by the active layers is not blocked by the planet and there is no pile up of gas to sustain fast inward migration. A slow inward migration is observed, similar to that seen in simulations of low viscosity discs without laminar accretion flows (Lega et al. 2021).

We remark that our result is different from the one of McNally et al. (2020), who found that planet migration was not affected by the presence of an accreting layer near the disc surface, but this is not surprising since the cited study concerned fully embedded low mass planets. For Jupiter-mass planets and a stellar accretion rate of 10⁻⁸ M_⊙ yr⁻¹, the separation between the two regimes, for a planet at 5.2 au occurs for a value of the vertically integrated column density of the active layer of about 0.2 g cm⁻², while in the case of a Saturn mass planet this happens at about 0.65 g cm⁻². We notice, however, that the column density giving the separation between the two regimes scales with the planet-star distance as so that a fast migrating giant planet should eventually enter in the regime of slow inwards migration once it is close enough to the star. The large number of giant planets in the warm-Jupiter region suggests that the column density of the active layer is of the order of 0.2–0.5 g cm⁻², which would stall a Jupiter-mass planet at ~4–1 au respectively.

It is noteworthy that although X-rays and cosmic rays are important sources of ionisation in protoplanetary discs, and penetrate into columns of depth ~ 10 g cm⁻² and 100 g cm⁻², respectively, the strong magnetic coupling in the surface layers of protoplanetary discs that allows a magnetised wind to be launched is thought to depend on the ionsation of sulphur and carbon atoms by uv photons from the star (Perez-Becker & Chiang 2011). The column density of this ionised layer has been calculated to be in the range 0.01 to 0.1 g cm⁻². For this range of values we expect Jupiter mass planets stalling in the Warm Jupiter region.

Our model is certainly simplistic with respect to a full non-ideal MHD model. However it has its merits: on the one hand a full MHD model with an embedded planet evolving on long timescales is not yet accessible with available computational resources; on the other hand our model captures the essential interaction of the planet with a layer torqued by some external process, and hence provides a framework with which to understand future more complex simulations when they become possible. From the results we have obtained, we expect that giant planets embedded in low viscosity stratified discs may have a rich and still not fully explored variety of dynamical paths, depending on the properties of the active layers. Among these paths there is vortex-driven migration (Lega et al. 2021), which is extremely slow after a short-ranged phase, and therefore has the potential to explain the orbits of cold Jupiters (at a few au from the central star) even if these planets formed no farther than the snowline location (5–10 au).

Acknowledgements

We acknowledge support by DFG-ANR supported GEPARD project (ANR-18-CE92-0044 DFG: KL 650/31-1). We also acknowledge HPC resources from GENCI DARI n.A0100407233 and from ‘Mesocentre SIGAMM’ hosted by Observatoire de la Côte d’Azur. LE wish to thank Alain Miniussi for maintenance and re-factorisation of the code FARGOCA.RPN acknowledges support from STFC through the consolidated grants ST/M001202/1 and ST/P000592/1.

Appendix A Three dimensional investigation of the 1:2 mean motion resonance

The radial flow analysis that we provided in the paper is a one dimensional analysis that provides useful information. However, in order to have a better insight on the flow behaviour, and specifically to analyse the most persistent resonant effect that we have observed in all the simulations presented in the paper, we provide in this appendix a two dimensional (r − z) view of the radial flow in the vicinity of the 1:2 mean motion resonance. Fig.A.1 shows the azimuthally averaged radial flow for simulations NW, Hin and Hout. An important amount of outflow is observed mainly up toz = h₀r, while the flow is strongly inwards in a thin layer located above z > h₀r. When comparing simulation Hin and Hout to the reference simulation NW we see that the torque imposed up to a disc scale height respectively decreases and enhance the natural outflow at the 1:2 resonance, according to the direction of the imposed flow. The imposed inflow, in simulation Hin, appears not enough to overcome the natural outflow associated to the 1:2 resonance and the resulting total flow is still positiveas observed in Fig.7, bottom left panel. The imposed outflow of in simulation Hout, on the contrary, contributes in enhancing the gas motion at the resonant location (Fig.7, bottom right panel).

Fig. A.1 Azimuthally averaged radial flow in the plane r − z at time t = 83004 yr for simulations sets NW, Hin and Hout. In simulations Hin and Hout a synthetic torque is applied below the disc scale height indicated in the figure by the white continous line. The dashed vertical line indicates the position of the 1:2 mean motion resonance considering planet position at time t = 83004 yr.

References

All Tables

All Figures

Fig. 1 Disc with radial transport ṀA = 10−8 M⊙ yr−1 generated in the active layer with Σ(z) ≤ 0.1 g cm−2. Top panel: radial velocity normalised with respect to the sound speed; bottom panel: total mass flow as a function of the distance from the star. The flow transported by the active layer is 75% of the nominal flow MA = 10−8 M⊙ yr−1 (see text); the flow due to viscosity is clearly negligible with respect to ṀA, so that the total flow is transported by the active layer and the green and red curve are superposed (the minus sign in the plot indicates inflow).

In the text

	Fig. 2 Vertically integrated surface density, Σ(r, z) = ∫, in simulation set S01 at time t = 9481 yr before the planet is released in the migration phase. The white dashed line corresponds to the ΣA = 0.1 g cm−2. The synthetic torque is applied above this line up to the disc surface. We remark that the formation of the gap implies that the synthetic torque is applied deep down in the disc.
In the text

	Fig. 3 Contours of surface density (multiplied by r1∕2 to flatten the radial profile) for simulation sets NW, S01, S1, S10, Hin, Hout (from top tobottom) with a migrating Jupiter-mass planet. The panels correspond to different times, reported on top of each panel.
In the text

	Fig. 4 Azimuthally averaged radial surface density profiles at two specific times (reported on top of the panels) within the rapid inward vortex-driven migration phase. No clear differences appear among the six simulations sets.
In the text

	Fig. 5 Evolution of the semi-major axes and eccentricities for simulations in layered discs with different column densities of the active layers (ΣA). For comparison, data for a classical viscous disc are also plotted (label ‘No Wind’).
In the text

Fig. 6 Migration speed vs. time for the six simulation sets. A first phase of inward vortex-driven migration is common to all the cases. Atabout t = 50 000 yr the migration paths bifurcate: for the very thin active layer of simulation S01 the final speed of migration is very similar to the NW case; increasing the depth of the layer a regime of faster inward migration is observed. When the active layer is at the midplane the planet migrates inwards (Hin) and outwards (Hout) at opposite speeds in the final phase of the simulation.

In the text

Fig. 7 Total mass flow transported through the disc (red curve). A temporal average over 1000 outputs (spaced by Δt = 11.85 yr) has been performed in order to smooth fluctuations. For the non-stratified case (top left panel) we have split the averaging interval into five parts. We show the five averages to illustrate the time variations of the flow (see text). In the averaging time interval the planet is migrating outwards for simulations NW, S01 and Hout and is migrating inwards for the cases with thicker active layers: S1 and S10 and for simulation set Hin. The black and red filled circles at Ṁ = 0 indicate the planet–star distance respectively at the beginning and at the end of the averaging interval. For the stratified disc models the contribution from the active layers (green curve) as well as the one fromthe dead region of the disc (blue curve) are also plotted. Notice that for simulations Hin and Hout we apply thesynthetic torque in the region between the midplane and z = H (active layers); while for z > H no torque is applied (dead layers). The grey horizontal line at Ṁ = ∓ 10−8 M⊙ yr−1 indicates the expected flow transported by the active layers (the sign plus is for simulation set Hout). The vertical dotted grey lines provide the location of the main mean motion resonances computed using the planet position at t = 71 146 yr. The evanescent boundary condition is applied in the domain inside the vertical black dotted line at 1.25 au.

In the text

Fig. 8 Azimuthally averaged radial surface density profiles at specific times (reported on top of the panels). At the time of top and middle panels, the migration behaviour has clearly bifurcated with slow outward migration for discs NW and S01 and rapidinward migration for the simulation sets with thicker active layers as well as for the case Hin. Rapid outward migration is observed for the case Hout. Despite the differences in the migration regimes, the gap depth slowly evolves from the time corresponding to the middle panel to theone represented in the bottom panel: gaps are deeper in all the simulations sets in which the planet blocks the flow from the active layers.

In the text

	Fig. 9 Specific torque γA and specific planetary torque exerted by a Jupiter mass planet on the outer disc (Eq. (10), see text) γp as a functionof ΣA. We notice that γA > γp only for the model with the thinnest layer (S01). The dot at ΣA = 100 g cm−2 refers to the case where the torque is applied to the layer around the midplane (Hin and Hout).
In the text

	Fig. 10 Evolution with time of semi-major axes for simulations in layered discs with different depths of the active layer ΣA) for the case of a Saturn mass planet. For comparison we also plotted (cyan dashed line) the data for simulation NW.
In the text

	Fig. A.1 Azimuthally averaged radial flow in the plane r − z at time t = 83004 yr for simulations sets NW, Hin and Hout. In simulations Hin and Hout a synthetic torque is applied below the disc scale height indicated in the figure by the white continous line. The dashed vertical line indicates the position of the 1:2 mean motion resonance considering planet position at time t = 83004 yr.
In the text