© ESO, 2014

1. Introduction

Planet formation occurs in accretion discs around young stars. The evolution of accretion discs is still not understood completely, but it takes place on a Myr time-scale during which the accretion rate of gas onto the host star drops (Hartmann et al. 1998). The formation of gas giants has to take place during this time. In the most commonly favoured scenario (Pollack et al. 1996) a massive solid core of several Earth masses has to form first to accrete gas. However, the formation of cores via the accretion of embryos and planetesimals is not very efficient and therefore can not explain the growth of proto-planetary cores (Levison et al. 2010). Recently, the idea of pebble accretion onto planetary embryos offered a possible solution to this problem, because it occurs on a very short time-scale (Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012).

Regardless of how the cores may form, they will migrate through the disc because of interactions with the gas (Ward 1997). As long as discs had been assumed to be locally isothermal, cores were expected to migrate very fast towards the central star, which was a problem for understanding the formation of giant planets. It has been recently shown that the migration of the planetary cores in the disc depends on the thermodynamics in the disc (Paardekooper & Mellema 2006; Kley & Crida 2008; Baruteau & Masset 2008; Kley et al. 2009). For a review on planet migration see Baruteau et al. (2014). The thermodynamics inside the discs are determined by viscous heating, radiative cooling, and stellar illumination (Bitsch et al. 2013, hereafter Paper I). Stellar irradiation maintains a flared disc structure in the outer parts of the disc (Chiang & Goldreich 1997; Paper I). A consequence of the disc structure on the migration of planets is described in Paper I, where we find that outward migration is possible in the regions of the disc where H/r decreases, which effectively means that the regions of outward migration are smaller in a flared disc than a shadowed disc, where H/r decreases in the outer parts of the disc.

Bitsch et al. (2014, hereafter Paper II) studied the evolution of discs where the mass flow towards the star (Ṁ) is independent of radius. In these accretion discs, the mass flux Ṁ is defined as (1)where ν is the viscosity, following the α-viscosity approach (Shakura & Sunyaev 1973). H is the height of the disc, and Ω_K is the Keplerian rotation frequency. Σ_G denotes the gas surface density. We have found that during the evolution of the disc (decreasing Ṁ) the regions of outward migration shrink until for small Ṁ no outward migration is possible. Consequently, all planets migrate inwards, similarly to the isothermal disc scenario. This, again, may be problematic for understanding giant planet formation. Moreover, even when regions of outward migration exist, they typically concern only planets more massive than ≈5 Earth masses. Smaller planets migrate inwards as they grow, and therefore it is debated whether they can achieve this threshold mass quickly enough to be captured in the region of outward migration. Additionally, the inward migration of small-mass planets might be even faster than estimated in Paardekooper et al. (2011), see Lega et al. (2014). All this indicates that we still miss a mechanism for preventing inward migration of low-mass protoplanets.

The inward migration of small-mass planets can be stopped by a bump in vortensity (Masset et al. 2006), which is defined as (∇ × v)/Σ, where v is the velocity of the gas. The bump in vortensity (positive radial gradient) results in a positive contribution of the barotropic part of the corotation torque that acts on the planet which can overcompensate for the negative Lindblad torque and halt its inward migration.

Additionally, the formation of planetesimals and planetary embryos can be aided by a positive radial pressure gradient in a protoplanetary disc (Johansen & Youdin 2007). This stops the inward migration of aggregates that undergo gas drag (Brauer et al. 2008) and significantly facilitates planetesimal formation (Bai & Stone 2010a,b). This indicates that a bump in vortensity and pressure would be very useful for the formation of planetesimals and planetary embryos by the core accretion scenario, which raises the question of how such a vortensity bump and a pressure bump can appear in the disc.

The radial changes of the vortensity and pressure fields can be triggered by changes in viscosity, which can be caused by changes in the gas ionisation fraction or condensation at the ice line (see Sect. 2.2). In Paper II we only explored discs that have a unique value of α in time and space. In this paper, we investigate the effect of viscosity transitions on the structure of the disc, with special attention devoted to the radial gradients of vortensity and pressure. Our aim is to determine the conditions for the following appealing scenario for giant planet formation.

At a pressure bump in the disc, particles can accumulate and form planetesimals and planetary embryos. These embryos would then migrate inwards, therefore we assume that a bump in vortensity exists at the same location in the disc, which can trap the planetary embryos (Masset et al. 2006; Morbidelli et al. 2008). After they are formed, these embryos could then grow via pebble accretion (Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012) until they reach a mass of ≈10–20 M_Earth and gas accretion onto the core can start (Pollack et al. 1996). The inward migration of cores of a few Earth masses can also be stopped by the entropy related corotation torque (Paardekooper & Mellema 2006; Kley & Crida 2008; Baruteau & Masset 2008; Kley et al. 2009), so that at this stage of evolution a bump in vortensity is not necessarily needed any more. As the planet starts to grow beyond ≈50 M_Earth, it opens up a gap in the disc (Crida et al. 2006) and is then released into type-II-migration (Lin & Papaloizou 1986) and slowly moves towards the star as the disc is accreted. The final stopping point of the type-II migrating gas giant planet is then set by the inner edge of the disc (Crida & Morbidelli 2007) or the photo-evaporation radius of the disc (Alexander & Pascucci 2012), which explains the pile-up of large gas giants in the inner systems.

This paper is structured as follows: first, we give an overview of the methods used, especially of the viscosity (Sect. 2). Then we compare in great detail a Ṁ = 1 × 10^-7 M_⊙/yr disc with and without transition in viscosity in Sect. 3. In Sect. 4 we investigate the changes of the disc structure as the disc evolves to smaller Ṁ. Different variations of viscosity transitions are then discussed in Sect. 5. In Sect. 6 we discuss which viscosity transition is actually needed to realise our planet formation scenario, and we discuss how realistic it might be. Finally, we give a summary in Sect. 7.

2. Methods

2.1. General setup

The general simulation setup follows the descriptions of Papers I and II. The protoplanetary gas disc is treated as a three-dimensional (3D) non-self gravitating gas whose motion is described by the Navier-Stokes equations. Without any perturbers in the disc, the disc has an axisymmetric structure. We therefore can compute all quantities in 2D, in the r − θ-plane, where θ is the colatitude of the spherical coordinates centred on the star, such that θ = 90° in the midplane. We used a 386 × 66 grid cells for all simulations. The dissipative effects are described via the standard viscous stress tensor approach (e.g. Mihalas & Weibel Mihalas 1984). We also included the irradiation from the central star, which was described in detail in Paper I. For that purpose we modified and substantially extended the existing multi-dimensional hydrodynamical code FARGOCA, as presented in Lega et al. (2014) and Paper II.

Viscous and stellar heating generates thermal energy in the disc, which diffuses radiatively through the disc and is emitted from its surfaces. To describe this process we used the flux-limited diffusion approximation (FLD, Levermore & Pomraning 1981), an approximation that allows the transition from the optically thick midplane to the thin regions near the disc’s surface; more details of this are given in Paper I. For comparison purposes with Paper II, we used a star with M_⋆ = 1 M_⊙ and R_⋆ = 3 R_⊙ and T_⋆ = 5600 K, which corresponds to L ≈ 9 L_⊙. The central star mainly influences the parts of the disc that are dominated by stellar irradiation. For the inner disc, we do not assume that this plays a significant role.

We adopted the opacity law of Bell & Lin (1994), with κ_⋆ = 3.5 cm²/ g as in Paper II. Throughout all the simulations presented here, we used a metallicity of 1% in μm size dust grains. We also apply the Ṁ boundary conditions described in the appendix of Paper II.

2.2. Prescription of viscosity transitions

Here we explore the effects of putative viscosity transitions at the ice line on the formation of planetesimals via the streaming instability as well as their capability of halting the migration of small planets. More precisely, we are interested in the viscosity gradients (that change the gradients in surface density and pressure) necessary to support our proposed scenario of planet formation. We modelled the transition in viscosity in a simple toy model, instead of using full-scale magneto-rotational-instability (MRI) simulations (Balbus & Hawley 1998). We made the general assumption that the α-parameter of our viscosity (Shakura & Sunyaev 1973) in the active layer (α_A) is much larger than in the midplane regions of the disc, which corresponds to a dead zone like structure (Balbus & Hawley 1998).

Pierens & Nelson (2010) have studied the influence of a dead zone on the accretion of gas onto giant cores in isothermal discs. They modelled the dead zone by an α transition in the vertical direction dependent on H in the disc. They found that the gas accretion process is not sufficiently altered by the presence of a dead zone, but that inward migration is somewhat reduced for giant planets. Here we instead focus on the migration of small-mass planets, and moreover, we link the transitions in α to the temperature and column density in the disc.

We took two main boundaries for the α transition into account. (i) An inner radial boundary that is dependent on the temperature (thermal ionisation); and (ii) a boundary that depends on the vertically integrated volume density (due to the absorption of X-ray and/or cosmic rays). In the midplane regions of the disc, we then explored the influence of an additional transition of viscosity at the ice line, where the different amount of condensed dust/ice grains there can quench the MRI at different degrees (Kretke & Lin 2007). The idea is not to give a priori quantitative estimates for these transition, but to understand the general behaviour of a disc that has viscosity transitions. Therefore the prescriptions also are slightly modified to test different parameters in Sect. 5.

The general assumption is that for T>T_crit = 1000 K the disc is fully ionised and therefore MRI active. This would correspond to the inner parts of the disc, and maybe to some super-heated parts in the atmosphere of the accretion disc. For T<T_crit the MRI activity has to be triggered by different sources. At the upper and lower layers of the disc, ionization originates from cosmic rays, stellar X-rays, ultraviolet photons and radionuclides. All these ionization sources have different penetration depths into the disc and therefore ionize the disc down to different levels.

Ultraviolet photons will only penetrate to a mass column of 0.1 g/cm² (Perez-Becker & Chiang 2011). X-rays can contribute down to Σ = 10 g/cm² (Bai & Goodman 2009), while cosmic rays can penetrate down to Σ = 100 g/cm² (Dzyurkevich et al. 2013). The different penetration depths for cosmic rays and X-rays suggests that an additional transition of α exists in the vertically integrated volume density.

Additionally, we imposed a transition of α at the ice line because there might be ice grains below the water condensation temperature that can then quench the MRI more effectively than in the region where only silicate grains are available (Kretke & Lin 2007). For the different α values we used we defined the following relation: (2)where we set α_A = 0.0264. This indicates that α_U = 0.00264 (α undead), which is still a significant, but low, alpha factor. In the α_D (α dead) region, the α factor is even smaller.

A zero viscosity in the dead zone can lead to outbursts in the disc triggered by gravitationally instabilities in the disc (Martin & Lubow 2014). Inside the dead zone, however, radionuclides can provide an additional source of ionization, and hydrodynamical instabilities might create turbulence, for instance, the vertical shear instability (Nelson et al. 2013) or the baroclinic instability (Klahr & Bodenheimer 2003), which will provide some viscosity. Therefore we assumed that α_D is not zero.

We set the following α regions in the disc: (3)

where → symbolizes a transition region for α. Σ_P(z) is the vertically integrated volume density from the top of the disc down to a height z, where the desired quantity of Σ_P is reached. This basically represents the penetration depth of cosmic and X-rays. We obtained the gas surface density by integrating both sides of the disc from top to midplane and from bottom to midplane. Therefore the gas surface density of the disc Σ_G is 2Σ_P(0). Transitions for the different α regions need to be smoothed out, because this is closer to real MHD simulations than a step function (Dzyurkevich et al. 2013).

Note here that if either T > 1000 K or Σ_P(z) < 10 g/cm², the α value corresponds to α_A, because the disc is either completely ionised by thermal ionisation or by X-rays and cosmic rays. The transition regions are smoothed in the following, linear way: and likewise for the other transition regions. Parts of the disc that are in these transition regions either follow the gradient of Σ_P or of T, depending on the steepness of these gradients (where the shallower gradient is favoured). Note here that the temperature range over which we imposed the α-transition (α_U → α_D) at the ice line is exactly the range over which the opacity changes in the opacity we used (Bell & Lin 1994). We therefore mark in our plots the ice line at 190 K, which is the middle of that transition range. The viscosity is given by (6)where Ω_K is the Keplerian frequency and c_s is the sound speed. Here we used the midplane sound speed c_s(z = 0), which provides a better control over the viscosity, which varies with height only because of changes in the α parameter. Additionally, the vertical changes in c_s are expected to be much smaller than the changes in α. The α map as a function of Σ_P and T is shown in Fig. 1 (top). We also display the α map of an Ṁ = 1 × 10^-7 M_⊙/yr disc (bottom panel) in equilibrium state (constant mass flux through each radial section of the disc) with α_A = 0.0264. The viscosity used for α_A in this work is generally ten times higher than in Paper II.

Fig. 1 Top: α as a function of ΣP and T as specified by Eq. (3). Bottom: α parameter in an Ṁ = 1 × 10-7 M⊙/yr disc that is discussed in Sect. 3. Inside the black line α ≤ αU, inside the white line α = αD.

3. Viscosity transitions

In Paper II we focused on discs with constant Ṁ. These discs also featured a constant α parameter in the whole disc. As stated above, the α parameter of the viscosity now follows some transitions, so that α changes in the disc with r and z because the temperature and the volume density are also functions of r and z. The radial and vertical changes in α will lead to changes in volume density relative to the uniform α case, because Ṁ has to be independent of r in the steady state. We compare here a disc with the described transition in α with a disc that has constant α throughout the disc equal to the value of α_A.

3.1. Disc structure

In Fig. 2 (top) the midplane temperature profiles for Ṁ = 1 × 10^-7 M_⊙/yr discs with and without a transition in α are shown. Up to ≈20 AU the midplane temperature is smaller in the case with α transition than the case without α transition. Consequently, the H/r profile of the disc with α transition is smaller than the H/r profile of the disc without α transition (second panel in Fig. 2). Here we define H as the pressure scale height using the midplane sound speed c_s(z = 0) with H = c_s(z = 0)/Ω_K. The bumps in the H/r and T profiles at ≈4 AU are also caused by transitions in the opacity profile (see Paper II), which changes the heating and cooling properties of the disc.

In a 1D radial model of a disc with constant radial Ṁ, a reduction of α is compensated for by an increase of Σ_G of the same magnitude, as Σ_G ∝ Ṁ/ν. This would allow for the same amount of viscous heating in the 1D disc and hence the same temperature profile. However, in r − z-discs this is not the case. A reduction in α does not imply an increase of Σ_G by the same factor, because the active layer of the disc can carry most of the accretion rate (see Fig. 5), which is not possible in 1D models. Indeed, the surface density profile in the inner part of the disc with α transition (third from top in Fig. 2) only shows an increase by a factor of 3 in contrast with the reduction in α by a factor of 100 in the dead zone.

Fig. 2 Midplane temperature (top), H/r (second from top), integrated surface density ΣG (third from top) and midplane pressure (bottom) for the Ṁ = 1 × 10-7 M⊙/yr disc with and without an α transition.

In the outer parts of the disc, the surface density is lower for the disc with an α transition. This is caused by the higher H/r in the outer parts, which results in a higher viscosity and hence in a lower surface density for an equal Ṁ disc.

In the disc with an α transition, the surface density shows a significant local maximum at ≈4 AU, which is exactly where we set the transition of α_U to α_D. This is because, as Ṁ is independent of r, a change in viscosity is compensated for by a change in Σ_G, resulting in the high-density regions close to the location of the change in α. This local maximum in surface density at ≈4 AU is very appealing, because it might function as a planet trap (by changing the disc vortensity) and stop inward migration of low-mass planets that migrate in the type-I-migration regime (Masset et al. 2006). In contrast, the disc without an α transition shows a dip in surface density around ≈4 AU, which is caused by the transition in opacity at the ice line and the corresponding increase in viscosity caused by the increase in H/r (see Paper II).

The midplane pressure of the same discs is displayed in the bottom panel of Fig. 2. For the disc without an α transition, the pressure follows a power law without any derivations. No pressure bumps exist in this case. In the case with an α transition, however, small wiggles in the pressure are visible, especially at ≈4 AU, where α_U transitions to α_D.

Despite the bumps in surface density, no real pressure bump exists in the disc because the gradient of pressure is not inverted. In 2D discs the pressure is defined as (7)with R_gas being the gas constant and μ the mean molecular weight. Here s denotes the power-law index of the surface density profile and β the power-law index of the temperature profile. This indicates that when a jump in surface density can compensate for the gradient in T, a pressure bump exists. However, in 3D discs the pressure is related to the volume density ρ_G as (8)as in hydrostatic equilibrium ρ_G = Σ_G/, where H is the disc vertical thickness. Here f denotes the flaring index of the disc. By comparing P_3D with P_2D, we immediately see that for the same temperature profile a steeper positive radial surface density gradient is needed in 3D than in 2D to invert the radial pressure gradient, as long as the flaring index of the disc is not strongly negative, which generally is the case. In the models shown in Fig. 2, the local maximum of the surface density gradient is not large enough to result in a pressure bump.

Fig. 3 Temperature (top) and density (bottom) of the Ṁ = 1 × 10-7 M⊙/yr with α transition. The black and white lines mark the transitions in the α viscosity parameter as defined by Eq. (3) and as shown in Fig. 1. The grey line marks the disc’s pressure scale height H.

In Fig. 3 the (r,z) temperature (top) and density (bottom) distribution are shown. In the inner parts close to midplane, the disc is heated by viscosity. With increasing distance to the star, this heating decreases as the density decreases. The upper layers of the disc are heated by the star, therefore they are hotter than the midplane region. The regions of the disc far away from the star with low density are heated by heat diffusion downwards from the upper layers. Close to the α transition at ≈4 AU, the temperature is vertically constant for a very wide vertical range (up to ≈0.3 AU).

The volume density ρ_G near the midplane is nearly radially constant around the α transition at ≈4 AU, as shown in Fig. 3 (bottom). In this region the density is vertically constant up to ≈0.1 AU. To calculate the surface density, the density shown here was vertically integrated, which resulted in the local maximum in the surface density (middle panel in Fig. 2). But, this is not enough to generate a pressure bump because T decreases with r.

The radial velocity of the disc without an α transition is displayed in the top panel of Fig. 4. A positive value indicates an outward motion, while a negative value indicates an inward motion. Here, we see an outflow in the midplane regions of the disc and an inflow at the top layers of the disc, which was also observed by Urpin (1984), Kley & Lin (1992), Takeuchi & Lin (2002). This flow is caused by the α prescription for the viscosity (Fromang et al. 2011).

The radial velocity of the disc with α-transition is displayed in the bottom panel of Fig. 4. We now observe several layers in the flow pattern. In the very top regions of the disc (where α=α_A) the flow is directed inwards. Farther down, the flow is directed outwards, while it then flows inwards again at the α_U transition. The flow farther inside the disc near the midplane is again directed outwards. Both discs show similar inflow behaviour in the top layers.

Fig. 4 Radial velocity of the disc without α transition (top) and with α transition (bottom). The black and white lines mark the transitions in the α viscosity parameter as defined by Eq. (3) and as shown in Fig. 1. In the bottom panel, the parts of the disc above the red line are fully active with αA. They are overplotted in the top plot as well to show where these transitions would be. The grey line represents the pressure scale-height H. Positive velocities mark an outward motion, while negative velocities indicate an inward motion.

Fig. 5 Ṁ rates for the different layers of the Ṁ = 1 × 10-7 M⊙/yr disc with an α transition. Here, a positive Ṁ rate indicates accretion, while a negative Ṁ indicates decretion (motion away from the star).

In Fig. 5 the Ṁ rates for the different α-layers of the disc with α transition are displayed. The layers are defined by their α values as indicated in the figure. Additionally, the total Ṁ rate is displayed, which is the sum of all the different layers. While the total Ṁ rate remains approximately constant in radius, the different Ṁ rates of the different layers vary strongly. Between 2–4 AU, the Ṁ rate in the active layer is highest (in inward direction) and in the same region, the layer with α<α_U has the largest outflow. At 4 AU the flux changes and each layer seems to carry the same amount of Ṁ. At 8 AU, α>α_U, so that only two layers in the disc exist. At that distance the flow returns to the state of an unlayered disc, with outflow in midplane and inflow in the top layers. We recall that the accretion rate of the active layer only depends on the radial velocity of the gas, because its density is fixed to 10 g/cm² at each r by construction.

In a simple view of a dead zone, with zero viscosity, the material flows inwards in the top layers of the disc and then dumps material into the dead zone, which then grows in density. Eventually, the dead zone becomes so massive that the disc becomes gravitationally unstable and produces outbursts of accretion (Martin & Lubow 2014). But this view is not supported by the flow picture we found in our simulations. Just above the region of low viscosity, we found outflow and in the very top layers inflow. This suggests that the mass flux in the active layer can adjust to the dead zone so that the mass flux is still constant in r and at the same time does not dump that much material in the dead zone. This consideration is also supported by the fact that we did not find an increase of Σ_G of a factor of 100 in the dead zone in our simulations, where α is reduced by a factor of 100. This basically suggests that the discs cannot become gravitationally unstable which disagrees with the hypothesis of Martin & Lubow (2014).

3.2. Planetesimal formation and embryo migration

The reduction of the gravitational force by the radially outwards-pointing force of the pressure gradient causes a difference Δv between the azimuthal mean gas flow and the Keplerian orbit. This is given by (9)where is the Keplerian velocity c_s/v_K = H/r (for the H/r profile, see the second panel from top in Fig. 2). η represents a measure of the gas pressure support (Nakagawa et al. 1986).

Because of the sub-Keplerian rotation of the gas, small solid particles drift towards the star. A local reduction in facilitates particle clumping in the streaming instability (Johansen et al. 2007; Johansen & Youdin 2007; Bai & Stone 2010a,b). Reduced would also help planetesimals to grow further by pebble accretion (Lambrechts & Johansen 2012).

The parameter is displayed in Fig. 6. In addition to the Ṁ = 1 × 10^-7 M_⊙/yr discs with and without α transition, we display the parameter in the standard minimum mass solar nebular (MMSN) case (Hayashi 1981). The MMSN model is constructed with simple power laws in surface density and temperature, which results in a power law for the parameter as well. The increasing with r means that it is much harder to form planetesimals via the streaming instability at larger distances.

Fig. 6 of the Ṁ = 1 × 10-7 M⊙/yr disc with and without an α transition. The black line marks the Δ parameter in the standard MMSN disc.

The disc without transition in α has a maximum at r ≈ 3.5 AU and a drop for larger r that is caused by the maximum in the H/r profile, which is originally caused by transitions in opacity (Paper II). That has the same value at 2 AU and at 5 AU implies that the formation of planetesimals at these locations are equally likely. But at 5 AU the disc is much colder than at 2 AU, which allows for the condensation of ice grains, meaning that more solid material is available at 5 AU. In the end, this might actually facilitate planetesimal formation more at 5 AU than at 2 AU.

In the case with an α transition is much reduced at 4 AU. This is the exact location of the α transition from α_U to α_D. This transition causes a wiggle in the midplane pressure (bottom panel in Fig. 2). The change of the pressure gradient then significantly reduces at that location, making planetesimals formation via the streaming instability much more likely there.

After they are formed, small planetesimal can migrate inwards. These planetesimals can then grow to planetary embryos by oligarchic growth (Kokubo & Ida 1998) or by pebble accretion (Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012). We now shift our attention to planet migration. To measure the migration rate, one should in principle extend our disc model to full 3D (with sufficient resolution in azimuth), introduce a planet in the disc and measure the torques. This would be very computationally expensive, however, in particular because for the small-mass planets that we are interested in here, a very high resolution would be required to resolve the planet’s horseshoe region. To estimate the smallest mass needed for outward migration of embryos triggered by the entropy related corotation torque, we therefore used the formula by Paardekooper et al. (2011), which captures the effects of torque saturation due to viscosity effects in contrast to Paardekooper et al. (2010), where the torques are fully unsaturated. The torque formula by Paardekooper et al. (2011) was observed to match with the torques in 3D simulations quite well (Bitsch & Kley 2011). However for small-mass planets, a new phenomenon has been observed that can result in a more negative torque that drives a faster inward migration than described by the formula of Paardekooper et al. (2011) (Lega et al. 2014). However, this phenomenon deserves more detailed studies. We therefore used the formula by Paardekooper et al. (2011) as an estimate for planet migration in gas discs. This was made in the same way as described in Papers I and II, so we do not state the torque formula explicitly. The total torque acting on an embedded planet is a composition of its Lindblad torque and its corotation torque: (10)The Lindblad and corotation torque depend on the local radial gradients of entropy S ∝ r^{− ξ}, with ξ = β − (γ − 1.0) s, and gas surface density Σ_G ∝ r^{− s}. β describes the radial gradient of the temperature profile, T ∝ r^−β. Very approximately, for large s a large ξ caused by a large β will lead to outward migration, while a flat radial entropy gradient will lead to inward migration.

Fig. 7 Migration map for discs without (top) and with (bottom) an α transition. The encircled regions in black mark regions in which outward migration is possible. The vertical red line marks the ice line at 190 K in both plots. The two plots feature slightly different radial extensions and planetary masses.

The migration maps for discs without (top) and with (bottom) an α transition are displayed in Fig. 7. In the case without an α transition, two regions of outward migration are visible, one around 2 AU and one between 4 and 7 AU. For these regions of outward migration a planetary mass of at least ≈17 M_Earth is needed, which is very high. The reason for this lies in the relatively high viscosity with α = 0.0264 (Paper II).

In the case of an α transition this changes dramatically. There are still two regions of outward migration, but they are much more narrow in radial distance than for the case without an α transition.

More precisely, now the region of outward migration is only a small band between 3.5 to 4.0 AU. This is because farther away the viscosity is too low to keep the entropy-driven corotation torque unsaturated. The outer migration region is now due to the vortensity-driven corotation toque, which is positive and locally exceeds the Lindblad torque because of the positive radial gradient of the surface density. This is the planet-trap mechanism first studied in Masset et al. (2006)¹. The planet trap mechanism also changes the smallest and largest masses for outward migration, relative to the entropy-driven corotation case. The smallest mass for outward migration is now ≈0.5 M_Earth, which is much smaller than before. On the other hand, the largest mass for undergoing outward migration is reduced to ≈23 M_Earth.

Planetesimal formation for discs with and without the α transition, the streaming instability is more likely to operate in a region of outward migration. However, this may not be enough to explain the formation of a giant planet at the snowline (i.e. at the location of the α-transition). In fact, only objects more massive than 0.5 M_Earth avoid inward migration in this region. Thus, if the streaming instability and the subsequent pebble accretion (Lambrechts & Johansen 2012) do not form a 0.5 M_Earth object sufficiently fast, the object will have the time to migrate into the inner system. Therefore, a steeper viscosity gradient than assumed here may be needed to allow a giant planet to form at the snow line, because a steeper surface density gradient (capable of stopping the inward migration of smaller-mass planets) and possibly even a pressure bump (preventing the inward migration of planetesimals) is needed.

4. Disc evolution

We follow here the prescription for the transition in α from Eq. (3) and decrease the disc’s Ṁ by decreasing the Σ_G value to follow the evolution of the disc through different accretional stages (same approach as in Paper II). We expect that with the decrease of the surface density of the disc, the radial gas-flow pattern eventually returns to the pattern of an unlayered disc. This is because the vertical integrated surface density from the disc’s surface down to the midplane, which cannot exceed Σ_G/ 2, will not exceed the 100 g/cm² limit to allow a full transition from α_U to α_D (see Eq. (3)).

Fig. 8 Surface density (top) and midplane temperature (middle) for discs featuring different Ṁ rates with α transitions. The bottom plot features of the streaming instability for the same discs.

The surface density (top) and midplane temperature (middle) are displayed in Fig. 8. The ratio between the surface densities in the Ṁ = 1 × 10^-7 M_⊙/yr and the Ṁ = 1 × 10^-8 M_⊙/yr is higher than a factor of 10 (the ratio of the accretion rates) in the inner parts of the disc. This is because the disc with the smaller Ṁ has a larger vertically averaged α (because the column density of the gas is lower). Consequently, the full surface density is reduced more than proportionally to the accretion rate. The outer parts of both discs are fully active, so that the difference in Σ_G is ≈10, i.e. proportional to the change of the Ṁ rate. This means that the change in total surface density Σ_G is not linear in the change of Ṁ, because the change of viscosity is not linear in the first place.

A local maximum in the surface density (positive surface density gradient) is visible only for the Ṁ = 1 × 10^-7 M_⊙/yr disc. For the Ṁ = 5 × 10^-8 M_⊙/yr disc a dip in the surface density profile at ≈2–3 AU is visible, which is at the outer edge of the α_U layer. In fact, inside of ≈2.5 AU, this reduction in α causes the surface density to increase. For lower Ṁ values no such dips are visible, because the highest surface density of these discs at any radius is lower than the condition for reducing α_A fully towards α_U.

is displayed in the bottom panel of Fig. 8. The dip in that was visible for the Ṁ = 1 × 10^-7 M_⊙/yr disc vanishes for smaller Ṁ. This is because the transition in α form α_U to α_D disappears, because the disc’s surface density is decreased. The big bump of at ≈2 AU in the Ṁ = 5 × 10^-8 M_⊙/yr disc is caused by the very steep negative surface density gradient in this region, which results in a steeper negative pressure gradient (Eq. (8)), which increases (Eq. (9)). In the later stages of the disc evolution, becomes much smaller in the inner regions of the disc, which in principle facilitates planetesimal formation by the streaming instability. However, at such late stages the lifetime of the disc is quite short, so that there may not be enough time to form the cores and then accrete gas to form giant planets.

Fig. 9 Migration map for discs with an α transition. The discs feature Ṁ = 5 × 10-8 M⊙/yr (top) and Ṁ = 1 × 10-8 M⊙/yr (bottom). The regions encircled in black mark regions in which outward migration is possible. The vertical red line in the top plot marks the ice line at 190 K. The bottom plot has a different colour scale.

The migration maps for discs with Ṁ = 5 × 10^-8 M_⊙/yr (top) and Ṁ = 1 × 10^-8 M_⊙/yr (bottom) discs are displayed in Fig. 9. A large region of outward migration exists around ≈3–4.5 AU in the Ṁ = 5 × 10^-8 M_⊙/yr disc. This region is quite similar to the region of outward migration for the Ṁ = 1 × 10^-7 M_⊙/yr disc without an α transition shown in the top panel of Fig. 7. This is because the disc with Ṁ = 5 × 10^-8 M_⊙/yr has a surface density lower than 100 g/cm² beyond 3 AU and therefore is close to being fully ionised. The region of outward migration here is related to the entropy-driven corotation torque and not to the vortensity-driven corotation torque.

The Ṁ = 1 × 10^-8 M_⊙/yr disc features for 2 < r < 4 AU a region of strong inward migration, which is caused by the positive temperature gradient in the disc. In the very inner parts of the disc (r < 1.0 AU) a region of outward migration exists, because of the strong negative temperature gradient, as shown in the top panel of Fig. 8. This migration map is also similar to the migration map of a disc without an α transition. This means that the region of outward migration moves inwards in time as the disc evolves to smaller Ṁ.

The time the disc spends at high Ṁ rates (Ṁ > 5 × 10^-8 M_⊙/yr) is expected to be very short (a few 100 kyr, see Hartmann et al. 1998). But in the evolution from Ṁ = 1 × 10^-8 M_⊙/yr to 5 × 10^-8 M_⊙/yr the lowest mass for outward planet migration increases from 2.5 to 10 Earth masses. This means that the formation of a giant planet in this region can occur only if the growth rate of the core is fast enough. Otherwise the core eventually finds itself below the lower mass boundary of the outward migration region and starts to migrate towards the inner disc. Similarly, when the disc reaches Ṁ ≈ 1 × 10^-8 M_⊙/yr, the outward migration region disappears for all planetary masses. Thus, it is necessary that by this time, the core has reached a large enough mass to accrete a massive envelope and transit to a slow type-II migration mode. We now examine whether such a fast growth is possible. Considering pebble accretion, the growth rates of planetesimals can be estimated (Eq. (42) Lambrechts & Johansen (2012)) (11)where t_d defines the growth timescale, ρ_P the pebble density, and M₀ the initial seed mass that accretes the pebbles. Assuming that a fraction of 10% of solids is in pebbles and an initial seed mass of 10^-4M_E, the growth time-scale is about 600 kyr, which is similar to the time-scale of the disc evolution according to the observations of Hartmann et al. (1998). This indicates that it may be possible for the cores to grow at the same rate as the outward migration region shifts upward in mass while the stellar accretion rate is reduced (see top panel in Fig. 9).

5. Different α transitions

We explore here different prescriptions for the change of α as a function of Σ_G and T.

5.1. Change of the cosmic ray penetration depth

The penetration of cosmic and X-rays into the disc is of crucial importance for the drive of the MRI. It is still debated how deep cosmic rays penetrate the disc. We therefore changed the transition from α_A to α_D to model a reduced penetration depth of cosmic and X-rays in the following way: (12)This implies that α_D can still be reached at lower surface density values, hence it can be reached at lower Ṁ values. Accordingly, we expect to observe the features related to the presence of a dead zone (surface density bump, planet trap, etc.) for discs with smaller Ṁ with this prescription of α (Eq. (12)) and not with that of Eq. (3).

Fig. 10 Density (top) and migration map (bottom) of the Ṁ = 5 × 10-8 M⊙/yr disc with a lower penetration depth for cosmic rays. In the top plot, the black and white lines mark the transitions in the α viscosity parameter as defined by Eq. (12), the grey line marks the disc’s height H. In the bottom plot, the regions encircled in black enclose the regions of outward migration, while the vertical red line marks the ice line at 190 K.

In Fig. 10 the (r,z) density structure of a Ṁ = 5 × 10^-8 M_⊙/yr disc with the α transition defined in Eq. (12) is shown in the top panel. Clearly, the discs features a local maximum in density at ≈3.3 AU. This then translates into a local maximum in the vertically integrated surface density (not displayed). No real pressure bump (positive pressure gradient) is visible in this case either, only a flattening of the pressure gradient, because the negative gradient in T overcompensates for the positive gradient in ρ_G (see Eq. (8)).

The migration map (bottom panel in Fig. 10) shows a very similar pattern as that of the Ṁ = 1 × 10^-7 M_⊙/yr disc with the α transition described by Eq. (3) (Fig. 7). More precisely, the stopping point of inward migration is shifted slightly inwards to ≈3 AU. This is expected because the surface density gradients, the α distribution, and the temperature gradients are similar. For smaller Ṁ, the prescription of Eq. (12) is not enough to create a low viscosity region. This means that low Ṁ discs are always nearly fully ionised, and their structure is similar to that described above.

5.2. Change of the value of α_U

In this section the ratio between the different α is changed, that is, we modify Eq. (2). These three independent modifications are (13)where we kept the transition parameters for temperature and surface density as in Eq. (3). This means that the change of viscosity at the ice line is reduced, because we aim to keep the same total reduction of α from α_A to α_D. For the last test case, there is only one transition in α because the transition α_U → α_D does not exist any more.

Fig. 11 Surface density for Ṁ = 1 × 10-7 M⊙/yr discs with different ratios of αA/αU (red, blue, and black) and one (light blue) with αU → αD for 185 K <T< 195 K.

In Fig. 11, the surface density for the four different α ratios for a disc with Ṁ = 1 × 10^-7 M_⊙/yr is displayed. In the outer parts of the disc, the surface density remains the same for all α cases because the outer parts of the disc are fully ionised and hence always have α_A, which is the same in all simulations. In the inner parts of the disc, the surface density increases as α_U becomes smaller, because all discs have the same Ṁ.

The local maximum in surface density (positive surface density gradient) at the snow line also decreases when α_U decreases, because the contrast α_U to α_D decreases. In the case of α_A = 50α_U the positive surface density gradient at the ice line has completely vanished and only a flat plateau is visible. This might be enough to stop inward migration.

Fig. 12 Migration map for the disc with Ṁ = 1 × 10-7 M⊙/yr that features αA = 100αU. The regions encircled in black mark regions in which outward migration is possible. The vertical red line in the plot marks the ice line at 190 K.

The migration map for α_A = 10α_U is discussed in Sect. 3. We do not display the migration maps for α_A = 20α_U and α_A = 50α_U here because the migration map looks quite similar to that presented in Fig. 7 (bottom). The only difference is that the region of outward migration is smaller (e.g. outward migration is still possible only between 3 and 10 M_Earth). In Fig. 12 the migration map for α_A = 100α_U is displayed. Interestingly, this migration map still features a small region of outward migration, even though that there is only a plateau in the surface density and not a bump.

This is because the barotropic parts of the corotation torque still generate a positive (although weak) contribution, and the steep negative radial temperature gradient in the region produces a positive entropy related corotation torque. The viscosity, although low, is still high enough to prevent torque saturation (this is not the case in the example discussed in Appendix A).

5.3. Change of the transition α_U to α_D

The transition from α_U to α_D in Eq. (3) was implemented across a range of 60 K, because this corresponds to the temperature at which the opacity transition is smoothed at the ice line. However, 60 K can be considered a broad temperature range just for the sublimation and condensation of ice grains. We now change this range to 10 K. This is the only modification to Eq. (3). The corresponding surface density is shown in Fig. 11 as the light-blue line.

A much steeper surface density gradient at the location of the ice line is clearly visible. This change of surface density will result in a planet trap, similar to Fig. 7, but with a higher intensity. This also traps planets with masses smaller than 0.5 M_Earth (Fig. 13). Note that an inverse pressure gradient now exists (Fig. 14). An inverse pressure gradient in a protoplanetary disc stops the inward migration of small bodies that undergo gas drag (Brauer et al. 2008), promotes the onset of the streaming instability (Johansen & Youdin 2007), and significantly facilitates planetesimals formation (Bai & Stone 2010a,b), as also discussed in Sect. 3.2.

Fig. 13 Migration map for the disc with Ṁ = 1 × 10-7 M⊙/yr with steep gradient in the α transition. The encircled regions in black mark regions in which outward migration is possible. The vertical red line in the plot marks the ice line at 190 K.

Fig. 14 Midplane pressure for the Ṁ = 1 × 10-7 M⊙/yr discs with different transition ranges in the temperature for αU to αD.

The maximum steepness a density gradient can have is constrained by the disc height H because a density bump would be unstable to the Rayleigh instability if its width were smaller than H (Yang & Menou 2010). Here the transition occurs on a length scale of 2H.

6. Gradients at the α transition

We have presented here detailed simulations of discs with transitions in the α parameter that resulted in changes in the disc structure. If the α transition is large enough, small mass planets can be trapped. The smallest mass of trapped planets decreases with increasing value of the radial positive surface density gradient. At the same location of the disc, of the streaming instability (Eq. (9)) is greatly reduced, which facilitates planetesimal formation. Again, an inversion of the pressure gradient is only visible if the disc features a very steep positive surface density gradient that can overcompensate for the negative temperature gradient (Eq. (8)).

We now investigate which changes in α, or more precisely, in viscosity are needed to create a pressure bump in the disc. In Fig. 15 the pressure gradient, dln(P)/dln(r), is displayed as a function of the viscosity gradient, dln(ν)/dln(r). In addition to the simulation data, we provide a rough fit through the data points to estimate the viscosity gradient necessary to cause a pressure bump in the disc.

Fig. 15 Pressure gradient at the α transition as a function of the viscosity gradient. The data are taken from simulations of the discs with Ṁ = 1 × 10-7 M⊙/yr and different α transitions.

The change in viscosity required to observe a pressure bump is quite large, dln(ν)/dln(r) ≈ −28. The reduction in viscosity needs to be this high because of the vertical structure of the disc. As stated above, in a 1D (radial) disc, a reduction in viscosity will be equally compensated for by an increase in surface density. If the vertical structure is taken into account, however, this is different because most of the accretion rate can be carried by the active layer (see Fig. 5), which results in a lower density increase in the midplane regions of the disc. Clearly, a pressure bump is much harder to achieve in a disc with a layered vertical structure than in a 1D disc structure.

The only simulation in which we observed a bump in the midplane pressure featured a jump of α by a factor of 10 over a temperature range of 10 K, which corresponds to 0.13 AU in this simulation. We justified this jump in α by more particles being available beyond the ice line (r>r_ice) that can quench the MRI (Kretke & Lin 2007). It seems realistic that the temperature range over which water vapour condensates into ice grains is only 10 K. But a change of α of a factor of 10 caused by more ice grains is doubtful. A smaller reduction of α seems more reasonable. However, if α is only changed by a factor of 2, no bump in the midplane pressure can be observed. This indicates the reduction in α needed to create a pressure bump is ≈5.

7. Summary

We have investigated the influence of viscosity transitions on the structure and migration rate of planets in accretion discs, which feature the same mass flow Ṁ through every radial section of the disc. The structure of the disc was calculated by using 2D hydrodynamical simulations that feature viscous and stellar heating as well as radiative cooling. The viscosity transitions were implemented by reducing the α parameter of the viscosity prescription. When the viscosity is reduced in parts of the disc, the disc’s flow adapts to keep the vertically integrated Ṁ independent of radius.

For high Ṁ rates, the gas surface density of the disc Σ_G is high enough to shield the inner parts of the disc from cosmic and X-rays, resulting in a region of reduced α, the dead zone. This region of reduced viscosity is denser than the active layers, which creates a local maximum in surface density. Because we imposed a radial α transition at the ice line as well, we created a planet trap at ≈4 AU. This planet trap is very efficient in trapping small-mass planets (M_P > 0.5 M_Earth), in contrast to discs without a transition in α, where a much larger planetary mass is needed to stop inward migration by the positive corotation torque (Fig. 7).

As the disc evolves to lower Ṁ rates, the surface density decreases, so that the disc is unable to completely shield itself from cosmic and X-rays. This implies that as the disc evolves, the regions of low α shrink and the disc finally reaches a state without reduced α. In this case the only possibility to stop planet migration is the positive entropy related corotation torque (when it exists).

The migration of small solids (pebbles) is instead related to gas drag and is sensitive to how sub-Keplerian the gas-disc is. A pressure bump is needed to stop the inward migration of pebbles. To create a pressure bump a positive radial volume density gradient is needed that is strong enough to compensate for the negative temperature gradient. This can only be achieved if the viscosity transition is very steep. More precisely, our simulations show that a reduction of about a factor of 5 of α is needed over a temperature range of ≈10 K.

We explained the reduction in α at the ice line by more grains that can block the MRI (Kretke & Lin 2007). This view is under debate and the magnitude of this change is not known.

Nevertheless, we can conclude here that quenching the MRI at the ice line by ice grains (Kretke & Lin 2007) can be an effective trap for Earth-mass planets. This does not necessarily imply a pressure bump in the disc structure, however small dust and ice particles can be trapped there. It implies a significantly reduced pressure gradient, which in turn can significantly facilitate planetesimal formation by the streaming instability (Bai & Stone 2010a,b).

Acknowledgments

B. Bitsch and A. Morbidelli have been partly sponsored through the Helmholtz Alliance Planetary Evolution and Life. We thank the Agence Nationale pour la Recherche under grant ANR-13-BS05-0003-01 (MOJO). We also thank A. Johansen for useful discussions. In addition we thank an anonymous referee for his/her input. The computations were made done on the Mesocentre SIGAMM machine, hosted by the Observatoire de la Côte d’Azur.

References

Appendix A: Interchange between Σ and ν

In Paper II we mainly focused on discs that have a ten times higher surface density and a ten times smaller α than the discs presented here, which gives the same Ṁ rate. Because the penetration of cosmic rays into the disc is dependent on the discs surface density, this is a crucial quantity for developing a low viscosity region through the shielding of cosmic rays. Here, we take a disc with Ṁ = 1 × 10^-8 M_⊙/yr, but with the viscosity parameter α_A = 0.003, which is ten times lower than presented in Sect. 3.

Fig. A.1 Surface density for Ṁ = 1 × 10-8 M⊙/yr disc with αA = 0.003, which is ≈10 lower than in Fig. 2. The Ṁ = 1 × 10-7 M⊙/yr disc with larger αA is overplotted for comparison.

In Fig. A.1 the surface density of the Ṁ = 1 × 10^-8 M_⊙/yr disc with α_A = 0.003 is displayed. The prescription of the α transition follows the one given in Eq. (3). As in Fig. 2 (middle), a bump in the surface density profile is visible.

Fig. A.2 Migration map for Ṁ = 1 × 10-8 M⊙/yr disc with αA = 0.003, which is ≈10 lower than in Fig. 9. The vertical blue line marks the ice line at T = 190 K.

However, the migration map (Fig. A.2) does not show a positive torque that acts on the planet at any radial distance. The negative torque is smallest in the region close to the surface density bump at ≈2.5 AU. This indicates that in principle the bump in the surface density should work as a planet trap, but as the α parameter in that region of the disc is very small (α_D = 3 × 10^-5), the torques do saturate and cannot sustain outward migration any more (Masset et al. (2006)). This result shows that the parameter range of the α transitions that allow the formation of a planet trap is very narrow.

All Figures

	Fig. 1 Top: α as a function of ΣP and T as specified by Eq. (3). Bottom: α parameter in an Ṁ = 1 × 10-7 M⊙/yr disc that is discussed in Sect. 3. Inside the black line α ≤ αU, inside the white line α = αD.
In the text

	Fig. 2 Midplane temperature (top), H/r (second from top), integrated surface density ΣG (third from top) and midplane pressure (bottom) for the Ṁ = 1 × 10-7 M⊙/yr disc with and without an α transition.
In the text

	Fig. 3 Temperature (top) and density (bottom) of the Ṁ = 1 × 10-7 M⊙/yr with α transition. The black and white lines mark the transitions in the α viscosity parameter as defined by Eq. (3) and as shown in Fig. 1. The grey line marks the disc’s pressure scale height H.
In the text

Fig. 4 Radial velocity of the disc without α transition (top) and with α transition (bottom). The black and white lines mark the transitions in the α viscosity parameter as defined by Eq. (3) and as shown in Fig. 1. In the bottom panel, the parts of the disc above the red line are fully active with αA. They are overplotted in the top plot as well to show where these transitions would be. The grey line represents the pressure scale-height H. Positive velocities mark an outward motion, while negative velocities indicate an inward motion.

In the text

	Fig. 5 Ṁ rates for the different layers of the Ṁ = 1 × 10-7 M⊙/yr disc with an α transition. Here, a positive Ṁ rate indicates accretion, while a negative Ṁ indicates decretion (motion away from the star).
In the text

	Fig. 6 of the Ṁ = 1 × 10-7 M⊙/yr disc with and without an α transition. The black line marks the Δ parameter in the standard MMSN disc.
In the text

	Fig. 7 Migration map for discs without (top) and with (bottom) an α transition. The encircled regions in black mark regions in which outward migration is possible. The vertical red line marks the ice line at 190 K in both plots. The two plots feature slightly different radial extensions and planetary masses.
In the text

	Fig. 8 Surface density (top) and midplane temperature (middle) for discs featuring different Ṁ rates with α transitions. The bottom plot features of the streaming instability for the same discs.
In the text

	Fig. 9 Migration map for discs with an α transition. The discs feature Ṁ = 5 × 10-8 M⊙/yr (top) and Ṁ = 1 × 10-8 M⊙/yr (bottom). The regions encircled in black mark regions in which outward migration is possible. The vertical red line in the top plot marks the ice line at 190 K. The bottom plot has a different colour scale.
In the text

Fig. 10 Density (top) and migration map (bottom) of the Ṁ = 5 × 10-8 M⊙/yr disc with a lower penetration depth for cosmic rays. In the top plot, the black and white lines mark the transitions in the α viscosity parameter as defined by Eq. (12), the grey line marks the disc’s height H. In the bottom plot, the regions encircled in black enclose the regions of outward migration, while the vertical red line marks the ice line at 190 K.

In the text

	Fig. 11 Surface density for Ṁ = 1 × 10-7 M⊙/yr discs with different ratios of αA/αU (red, blue, and black) and one (light blue) with αU → αD for 185 K <T< 195 K.
In the text

	Fig. 12 Migration map for the disc with Ṁ = 1 × 10-7 M⊙/yr that features αA = 100αU. The regions encircled in black mark regions in which outward migration is possible. The vertical red line in the plot marks the ice line at 190 K.
In the text

	Fig. 13 Migration map for the disc with Ṁ = 1 × 10-7 M⊙/yr with steep gradient in the α transition. The encircled regions in black mark regions in which outward migration is possible. The vertical red line in the plot marks the ice line at 190 K.
In the text

	Fig. 14 Midplane pressure for the Ṁ = 1 × 10-7 M⊙/yr discs with different transition ranges in the temperature for αU to αD.
In the text

	Fig. 15 Pressure gradient at the α transition as a function of the viscosity gradient. The data are taken from simulations of the discs with Ṁ = 1 × 10-7 M⊙/yr and different α transitions.
In the text

	Fig. A.1 Surface density for Ṁ = 1 × 10-8 M⊙/yr disc with αA = 0.003, which is ≈10 lower than in Fig. 2. The Ṁ = 1 × 10-7 M⊙/yr disc with larger αA is overplotted for comparison.
In the text

	Fig. A.2 Migration map for Ṁ = 1 × 10-8 M⊙/yr disc with αA = 0.003, which is ≈10 lower than in Fig. 9. The vertical blue line marks the ice line at T = 190 K.
In the text