© ESO, 2016

1. Introduction

With the newest telescopes (either on the ground or in space), we now have access to a huge variety of observations of both protoplanetary disks and exoplanets. These observations combine to give a more global vision of planetary formation. It is very important, therefore, to understand how the disks select the planets that are to survive, but also how the planets reshape the disks. To achieve that, it is necessary to understand the disks’ structure, their evolution, and their interactions with the planets they host.

Observations by Beckwith & Sargent (1996) and Hartmann et al. (1998) estimated a typical disk lifetime of a few million years, consistent with the time required for gaseous planet formation involving gas accretion on solid cores (Pollack et al. 1996). However, Papaloizou & Terquem (1999) suggested that hot Jupiters may have formed at different locations than the ones where they are observed, involving some planetary migration. Long-known type I Lindblad inward migration was extensively described in Goldreich & Tremaine (1979), Artymowicz (1993) and Ward (1997), while Ward (1991) added an analytical description of the corotation torque. In particular, they estimated a migration timescale of a typical Earth-mass planet in a minimum mass solar nebula (MMSN; Weidenschilling 1977; Hayashi 1981) around 10⁵ yr. Korycansky & Pollack (1993), Alibert et al. (2005) and Ida & Lin (2008) note an inconsistency between the planet formation and migration timescales: planets should not have time to grow before they fall spiraling into their host star. Various attempts at slowing down the cores’ inward migration invoked a regular magnetic field (Terquem 2003), turbulent magnetism (Nelson & Papaloizou 2004), inner cavities (Kuchner & Lecar 2002), self-shadowing effects (Jang-Condell & Sasselov 2005) or strong positive density gradients (Masset et al. 2006). Paardekooper & Mellema (2008), Paardekooper & Papaloizou (2008), Baruteau & Masset (2008) derive analytical expressions from numerical simulations to model inward migration, while accounting for the disk profile and, in particular, its thermodynamical structure. Menou & Goodman (2004) estimate that opacity transitions could alter and even stop inward migration in steady state viscous α disks. Planets at these locations would be trapped, preventing their fall onto the star and favoring planet growth by creating a region where trapped planetary embryos would accumulate and collide. Paardekooper & Papaloizou (2009a) analyze the possibility to trap planets with positive surface mass density gradients, while Lyra et al. (2010), Bitsch & Kley (2011), Hasegawa & Pudritz (2011b) estimate the trap locations.

Most of the previous analysis of planet migration was enabled by modeling simplified disk structures using power-law prescriptions for the disk midplane density and temperature (Hasegawa & Pudritz 2011a; Paardekooper et al. 2011), density prescriptions with self-consistent 2D-temperature structure (Bitsch & Kley 2011; Bitsch et al. 2013) or steady-state accretion disk models (Bitsch et al. 2014). Though these disk models appear simpler on large scales, these approaches permitted the development of locally valid torque expressions, which are now widely used for modeling planetary migration. In this paper, we use the disk structure generated by the hydrodynamical numerical simulations described in Baillié & Charnoz (2014, hereafter BC14) and Baillié et al. (2015, hereafter BCP15). Coupling the disk dynamics, its geometry, and its thermal structure, this code models the long-term viscous evolution of α disks (Shakura & Sunyaev 1973) and is able to numerically retrieve some observational characteristics of protoplanetary disks. Based on Helling et al. (2000) and Semenov et al. (2003)’s opacity tables for the main elements that make up the disk dust, BCP15’s code relates the disk temperature to the disk dust composition and dust-to-gas ratio at that temperature. From the obtained disk structures in density and temperature, it is possible to derive the torques exerted by the disk on putative planets, that induce their inward or outward migration. These torques are particularly sensitive to the gradients of density and temperature. This leads logically to migration changes at the places where the gradients are most affected: at the sublimation lines of the main components of the dust and at the heat transition barrier that separates the inner disk under viscous heating domination and the outer disk that is dominated by stellar irradiation heating (as suggested by Hasegawa & Pudritz 2011b). BCP15 found that 10 M_⊕-planets could get trapped at these locations.

In this paper, we investigate how the planet mass affects its ability to get trapped and saved. We also analyze how growing planets may remain trapped as the disk evolves. We show that only intermediate-mass planets (a few tens of Earth masses) can get trapped at the sublimation lines (mainly the water ice, refractory organics and silicates lines) whereas the heat transition barrier might save more massive planets. However, the most massive planets at that heat frontier might open gaps and move from type I to type II migration if they become more massive than ~ 170 M_⊕. In addition, our disk simulations show the possibility for super-Earths to get trapped in several specific locations inside 2 au. The distribution of some of these traps could be related to the distribution of resonant pairs seen in exoplanet observations (Schneider et al. 2011; Wright et al. 2011). The super-Earth traps we found may also provide a plausible way to allow the in situ formation scenario of super-Earths inside 1 au, as described by Ogihara et al. (2015).

Section 2 details how Baillié et al. (2015)’s hydrodynamical code allows the viscous evolution of a protoplanetary disk to be tracked by accounting for the sublimation of the main components of the dust. Section 3 describes the interaction between the disk and potential planets, calculates the torques, and determines the locations of the planet traps and deserts. The influence of the planet mass on its ability to get trapped is detailed in Sect. 4. The impact of planetary growth and the possibility of opening a gap are investigated in Sect. 5, while the impact of the inner planetary traps on the close-in super-Earth formation scenarios is investigated in Sect. 6. Finally, Sect. 7 draws conclusions and details some perspectives to better model how planetary growth will affect the disk.

2. Methods

2.1. Disk evolution model

Following the model described in BC14 and BCP15, we consider a viscous α disk (Shakura & Sunyaev 1973) and we use similar terminology to that used in those papers. We set the turbulent viscosity to α_visc = 10^-2, as is commonly taken for disks around T Tauri stars in the absence of deadzones (Fromang & Papaloizou 2006). Using the mass and angular momentum conservation, Eq. (1) from Lynden-Bell & Pringle (1974) describes the time evolution of the surface mass density. $\frac{\mathit{\partial }\mathrm{\Sigma }\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\frac{\mathrm{3}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial r}}\left(\sqrt{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{\left(}\mathit{\nu }\mathrm{\right(}\mathit{r,t}\mathrm{)}\hspace{0.17em}\mathrm{\Sigma }\mathrm{(}\mathit{r,t}\mathrm{)}\hspace{0.17em}{\sqrt{\mathit{r}}}^{\mathrm{)}}\right)\mathit{,}$(1)where Σ is the surface mass density and ν is the viscosity.

Similarly to BC14 and BCP15, we define annuli of masses that are logarithmically distributed in radii between R_∗ and 1000 au. We then apply Eq. (1) to the 1D grid in surface mass density. We set the boundary conditions so that the flux at the inner edge cannot be directed outward (the disk cannot gain mass from the star). The inner mass flux is therefore the mass accretion rate of the disk onto the star.

At each time step and for each annulus, the disk midplane temperature, T_m(r), is calculated by accounting for stellar irradiation heating, radiative cooling, and viscous heating of the disk material. The iterative process detailed in BC14, by which the disk radial temperature profile is determined, estimates the vertical structure of the disk, the shape of the photosphere, and the composition of the disk midplane at the same time.

The angle at which the star sees the photosphere at a given radial location r is the grazing angle, α_gr(r). Equation (2) describes how the photosphere height H_ph(r) is related to this angle, which actually controls the amount of energy that the star provides to the disk photosphere. ${\mathit{\alpha }}_{\mathrm{gr}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\arctan \left(\frac{\mathrm{d}{\mathit{H}}_{\mathrm{ph}}}{\mathrm{d}\mathit{r}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right)\mathrm{-}\arctan \left(\frac{{\mathit{H}}_{\mathrm{ph}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}\mathrm{0.4}{\mathit{R}}_{\mathrm{\ast }}}{\mathit{r}}\right)\mathrm{·}$(2)Equation (18) from Calvet et al. (1991) describes the resulting heating in the midplane while the viscous heating contribution depends on the surface mass density and the viscosity (i.e., the midplane temperature): ${\mathit{F}}_{\mathit{v}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\left(\mathit{R}\frac{\mathrm{d\Omega }}{\mathrm{d}\mathit{r}}\right)}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{9}}{\mathrm{4}}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}$(3)Therefore, we solve the energy equation by considering the irradiation heating as an implicit function of the midplane temperature. Assuming the vertical distribution of material to follow the Gaussian law that results from a vertical isothermal distribution in hydrostatic equilibrium, we can use Eq. (A.9) from Dullemond et al. (2001) to calculate the ratio χ of the photosphere height H_ph(r) to the pressure scale height h_pr(r). This approximation is reasonable below a few pressure scale heights, where most of the disk mass is located. The irradiation term directly depends on the opacity of the medium, which in turn depends on the physical states of the various components of the disk material, and therefore depends on the disk midplane temperature as explained in Sect. 2.2. Then, for a given initial grazing angle, a temperature guess sets the disk local composition, its opacities, its pressure scale height, and therefore its photosphere height (through χ). We can then derive $\frac{\mathrm{d}{\mathit{H}}_{\mathrm{ph}}}{\mathrm{d}\mathit{r}}$ and use Eq. (2) to aim at a convergence on the grazing angle after a few iterations. At radial locations where the algorithm converges towards a solution involving a positive grazing angle, the disk photosphere is considered irradiated. Conversely, when the converged solution requires a negative grazing angle, we consider the disk local photosphere as shadowed by inner regions and therefore not directly irradiated by the star: we then use a simpler algorithm that cancels the contribution of the stellar irradiation. Thus, the geometrical structure (grazing angle, photosphere, and pressure heights), the temperature, and the disk composition are determined jointly by iterating numerically: the algorithm is thoroughly described in BC14.

2.2. Disk composition

BC14 used a very simple model of opacities: either all the dust particles were sublimated or none, depending on whether the local temperature was higher than 1500 K (considered as the sublimation temperature of the silicate dust) or not. BCP15 used a refined model of opacities based on more dust species with optical constants measured in lab experiments, to account for the variations of the dust composition as a function of the local temperature. The computation of the Rosseland mean opacities follows the procedure described by Helling et al. (2000) and Semenov et al. (2003). BCP15 considered the dust grains to be composed of olivine silicate, iron, pyroxene, troilite, refractory and volatile organics, and water ice, with the initial abundances described in Table 1. Sublimation temperatures are provided by Pollack et al. (1994), corresponding to gas densities of about 10^-10 g/cm³. At a given temperature, their opacity model considers a mixture of dust grains with rescaled abundances at that temperature to account for the sublimated components. These opacities vary by several orders of magnitude, in particular around the sublimation temperature of the dust’s main components.

The Rosseland (χ_R) and Planck (κ_P) opacity variations with temperature are presented in Fig. 1 together with the Planck mean opacities at stellar effective temperature T_∗ = 4000 K in extinction (${\mathit{\chi }}_{\mathrm{P}}^{\mathrm{\ast }}$) and absorption (${\mathit{\kappa }}_{\mathrm{P}}^{\mathrm{\ast }}$).

Fig. 1 Mean-opacity variations with local temperature. Black: Rosseland mean opacity in extinction. Red: Planck mean opacity in absorption. Green: Planck mean opacity in extinction at stellar irradiation temperature. Blue: Planck mean opacity in absorption at stellar irradiation temperature. The corresponding dust-to-gas ratio is displayed with a dashed line.

3. Results

To better compare our results with previous studies, we decide to model the evolution of a MMSN protoplanetary disk around a classical T Tauri-type young star. Our initial surface mass density profile is given by Weidenschilling (1977) with the scaling from Hayashi (1981): $\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{17}\hspace{0.17em}\mathrm{000}{\left(\frac{\mathit{r}}{\mathrm{1}\hspace{0.17em}\mathrm{au}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{kg}\hspace{0.17em}{\mathrm{m}}^{-2}\mathit{.}$(4)The central star is a classical T Tauri-type young star with constant M_∗ = 1 M_⊙, R_∗ = 3 R_⊙, T_∗ = 4000 K and ${\mathrm{ℒ}}_{\mathrm{\ast }}\mathrm{=}\mathrm{4}\mathit{\pi }\hspace{0.17em}{\mathit{R}}_{\mathrm{\ast }}^{\mathrm{2}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{B}}\hspace{0.17em}{\mathit{T}}_{\mathrm{\ast }}^{\mathrm{4}}$ along the simulation.

This choice of initial state is discussed in BC14 and is also validated by Vorobyov & Basu (2007) who show that the MMSN density profile could be interpreted as an intermediate stage along the evolution of a protoplanetary disk under self-regulated gravitational accretion.

3.1. Time evolution

Though the gas of this type of disk is believed to photo-evaporate by stellar irradiation in a few million years (Font et al. 2004; Alexander & Armitage 2007, 2009; Owen et al. 2010), BCP15 followed the evolution of this type of disk over 10 million years.

In addition to following the evolution of the surface mass density, BCP15 also derived the associated midplane temperature profiles and noticed temperature plateaux at the sublimation temperatures of the main dust components (referenced in Table 1). These plateaux drift inward as the star accretes material form the disk and the disk cools down with time. The mean locations of these sublimation plateaux are referred as “sublimation lines” and are represented with the black dotted and dashed lines in Fig. 2.

Fig. 2 Time evolution of the migration traps (blue +) and deserts (red ×) positions. The sublimation line positions and the heat transition radius are represented in black. Each subfigure shows the traps and deserts for a given planet mass.

3.2. Type I migration

The results reported in BCP15 clearly show that radial gradients of surface mass density and temperature that were expected to be mostly negative could strongly vary and even become positive locally when the disk evolves.

Following their methods, we derive the total torques that a viscously evolving disk would give to a putative planet of mass M_P, which is located at a radial distance r_P from the central star. We assume that the planet has already formed. Goldreich & Tremaine (1979), Ward (1988), Artymowicz (1993), Jang-Condell & Sasselov (2005) early described that a planet excites resonances in the disk: Lindblad resonances that were caused by the action of the spiral arms were induced by the planet and corotation resonances owing to the horseshoe region around the planet. Since we know the various components of the associated resonant torques (Sects. 3.2.1), assuming that the disk structure is not affected by the planet, we can derive the reaction torque exerted by the disk on the planet.

BCP15 also noticed that putative planets could experience total torques that could change sign when encountering these density and temperature gradient irregularities, which results in planets being possibly trapped at specific locations called “planet traps” and usually located at the outer edge of the sublimation regions of the different dust elements, or at the outer edge of the shadowed regions (regions for which the photosphere is no longer directly irradiated by the star). They also provided evidence that some regions of the disk called “planet deserts” can be totally depleted in planets. Though this was previously discussed in imposed-density or steady-static disks (Hasegawa & Pudritz 2011a; Paardekooper et al. 2011; Bitsch & Kley 2011; Bitsch et al. 2013, 2014), BCP15 was the first attempt at modeling these traps in viscously evolving disks.

In this paper, we study how these traps and deserts are affected by the mass of the planet. To that purpose, we need to fully account not only for the Lindblad torque but also for the detailed corotation torque, including saturation effects.

3.2.1. Lindblad and corotation torques

A planet excites Lindblad resonances of multiple orders in a disk (Goldreich & Tremaine 1979, 1980). Assuming thermal equilibrium and an adiabatic disk, Paardekooper & Papaloizou (2008) provide a formula for the total Lindblad torque exerted by a 2D laminar disk in the absence of self-gravity on a circular planet. Owing to thermal diffusion, waves propagate at a velocity between the isothermal sound speed (maximum thermal diffusion) and the adiabatic sound speed (no thermal diffusion). Paardekooper et al. (2011) then generalized the previous formula to account for this difference, therefore replacing the adiabatic index γ by an effective index γ_eff: ${\mathit{\gamma }}_{\mathrm{eff}}\frac{{\mathrm{\Gamma }}_{\mathrm{Lindblad}}}{{\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}}\mathrm{=}\mathrm{-}{\left(\mathrm{2.5}\hspace{0.17em}\mathrm{-}\mathrm{1.7}\frac{\mathit{\partial }\ln \mathit{T}}{\mathit{\partial }\ln \mathit{r}}\mathrm{+}\mathrm{0.1}\frac{\mathit{\partial }\ln \mathrm{\Sigma }}{\mathit{\partial }\ln \mathit{r}}\right)}_{{\mathrm{r}}_{\mathrm{P}}}\mathit{,}$(5)with $\mathit{q}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{planet}}}{{\mathit{M}}_{\mathrm{\ast }}}$ the mass ratio of the planet to the star, ${\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}\mathrm{=}{\left(\frac{\mathit{q}}{\mathit{h}}\right)}^{\mathrm{2}}\hspace{0.17em}\mathrm{\Sigma }\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}\hspace{0.17em}{\mathit{r}}_{\mathrm{P}}^{\mathrm{4}}\hspace{0.17em}{\mathrm{\Omega }}_{\mathrm{P}}^{\mathrm{2}}\mathit{,}$(6)$\mathit{h}\mathrm{=}\frac{{\mathit{h}}_{\mathrm{pr}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}}{{\mathit{r}}_{\mathrm{P}}}$, and Ω_P = Ω(r_P) the Keplerian angular velocity at the planet position in the disk.

In the isothermal case, γ_eff = 1, whereas in the adiabatic case, γ_eff = γ = 1.4. Following Paardekooper et al. (2011), the effective index is defined by ${\mathit{\gamma }}_{\mathrm{eff}}\mathrm{=}\frac{\mathrm{2}\mathit{Q\gamma }}{\mathit{\gamma Q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\mathrm{2}\sqrt{\mathrm{(}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{Q}}^{\mathrm{2}}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{16}{\mathit{Q}}^{\mathrm{2}}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}}\mathrm{+}\mathrm{2}{\mathit{\gamma }}^{\mathrm{2}}{\mathit{Q}}^{\mathrm{2}}\mathrm{-}\mathrm{2}}}\mathit{,}$(7)where Q accounts for the thermal diffusion: $\mathit{Q}\mathrm{=}\frac{\mathrm{2}{\mathit{\chi }}_{\mathrm{P}}{\mathit{r}}_{\mathrm{P}}}{\mathrm{3}{\mathit{h}}_{\mathrm{pr}}^{\mathrm{3}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}{\mathrm{\Omega }}_{\mathrm{P}}}$(8)and χ_P is the thermal conductivity at the planet location ${\mathit{\chi }}_{\mathrm{p}}\mathrm{=}\frac{\mathrm{16}\mathit{\gamma }\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\mathit{\sigma }{\mathit{T}}_{\mathrm{P}}^{\mathrm{4}}}{\mathrm{3}{\mathit{\kappa }}_{\mathrm{P}}{\mathit{\rho }}_{\mathrm{P}}^{\mathrm{2}}{\mathit{h}}_{\mathrm{pr}}^{\mathrm{2}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}{\mathrm{\Omega }}_{\mathrm{P}}}\mathit{,}$(9)with σ the Stefan-Boltzmann constant, ρ_P the density, κ_P the Rosseland mean opacity at the planet location, and the 16 factor is a correction introduced by Bitsch & Kley (2011) of a previous 4 factor from Paardekooper et al. (2011).

Tanaka et al. (2002) initially describe a torque arising from the density gradient, called “barotropic” corotation torque. Paardekooper & Mellema (2006) show that the total torque exerted on a planet by the disk increases with the disk opacity and describe how that torque excess is related to the fact that the coorbital region trailing the planet was hotter and underdense. (Baruteau & Masset 2008) then show that, in an adiabatic disk, these regions (a colder denser region leading the planet, and a hotter underdense one trailing the planet) could form by entropy advection. The positive resulting “entropic” corotation torque can then become larger than the Lindblad torque.

These two contributions are known to have both linear and non-linear contributions. In the usual range of viscosity (α_visc< 0.1), Paardekooper & Papaloizou (2009b) state that the corotation torques are mostly non-linear, owing to the horseshoe drag caused by the interaction of the planet with the gas in its vicinity (Ward 1991). As the horseshoe region is closed, it contains a limited amount of angular momentum and therefore is prone to saturation, which cancels the horseshoe contributions to the corotation torque.

Paardekooper et al. (2010) described the fully unsaturated horseshoe drag expressions for both entropic and barotropic (or vortensity) terms. Using the gravitational softening b = 0.4h_pr also used in Bitsch & Kley (2011) and Bitsch et al. (2014), the horseshoe drag torques read: $\begin{array}{ccc}{\mathit{\gamma }}_{\mathrm{eff}}\frac{{\mathrm{\Gamma }}_{\mathrm{hs}\mathit{,}\mathrm{entro}}}{{\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}}& \mathrm{=}& \frac{\mathrm{7.9}}{{\mathit{\gamma }}_{\mathrm{eff}}}\hspace{0.17em}{\left(\mathrm{-}\frac{\mathit{\partial }\ln \mathit{T}}{\mathit{\partial }\ln \mathit{r}}\mathrm{+}\mathrm{(}{\mathit{\gamma }}_{\mathrm{eff}}\mathrm{-}\mathrm{1}\mathrm{)}\frac{\mathit{\partial }\ln \mathrm{\Sigma }}{\mathit{\partial }\ln \mathit{r}}\right)}_{{\mathrm{r}}_{\mathrm{P}}}\\ {\mathit{\gamma }}_{\mathrm{eff}}\frac{{\mathrm{\Gamma }}_{\mathrm{hs}\mathit{,}\mathrm{baro}}}{{\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}}& \mathrm{=}& \mathrm{1.1}{\left(\frac{\mathit{\partial }\ln \mathrm{\Sigma }}{\mathit{\partial }\ln \mathit{r}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}_{{\mathrm{r}}_{\mathrm{P}}}\mathrm{·}\end{array}$In the general case (including possibly saturation), the total corotation torque is the sum of the barotropic and entropic contributions: ${\mathrm{\Gamma }}_{\mathrm{corotation}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{c}\mathit{,}\mathrm{baro}}\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{c}\mathit{,}\mathrm{entro}}\mathit{,}$(12)each of these contributions including a combination of non-linear (Eqs. (10), (11)) and linear parts (Eqs. (13), (14)). The fully unsaturated linear expressions are $\begin{array}{ccc}{\mathit{\gamma }}_{\mathrm{eff}}\frac{{\mathrm{\Gamma }}_{\mathrm{lin}\mathit{,}\mathrm{entro}}}{{\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}}& \mathrm{=}& \left(\mathrm{2.2}\mathrm{-}\frac{\mathrm{1.4}}{{\mathit{\gamma }}_{\mathrm{eff}}}\right)\hspace{0.17em}{\left(\mathrm{-}\frac{\mathit{\partial }\ln \mathit{T}}{\mathit{\partial }\ln \mathit{r}}\mathrm{+}\mathrm{(}{\mathit{\gamma }}_{\mathrm{eff}}\mathrm{-}\mathrm{1}\mathrm{)}\frac{\mathit{\partial }\ln \mathrm{\Sigma }}{\mathit{\partial }\ln \mathit{r}}\right)}_{{\mathrm{r}}_{\mathrm{P}}}\\ {\mathit{\gamma }}_{\mathrm{eff}}\frac{{\mathrm{\Gamma }}_{\mathrm{lin}\mathit{,}\mathrm{baro}}}{{\mathrm{\Gamma }}_{\mathrm{0}}\mathrm{\left(}{\mathit{r}}_{\mathrm{P}}\mathrm{\right)}}& \mathrm{=}& \mathrm{0.7}{\left(\frac{\mathit{\partial }\ln \mathrm{\Sigma }}{\mathit{\partial }\ln \mathit{r}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\right)}_{{\mathrm{r}}_{\mathrm{P}}}\mathrm{·}\end{array}$Thus, Paardekooper et al. (2011) express the barotropic and entropic torques that account for the saturation effects, as follows: $\begin{array}{ccc}{\mathrm{\Gamma }}_{\mathrm{c}\mathit{,}\mathrm{entro}}& \mathrm{=}& \mathit{F}\mathrm{\left(}{\mathit{p}}_{\mathit{\nu }}\mathrm{\right)}\mathit{F}\mathrm{\left(}{\mathit{p}}_{\mathit{\chi }}\mathrm{\right)}\sqrt{\mathit{G}\mathrm{\left(}{\mathit{p}}_{\mathit{\nu }}\mathrm{\right)}\mathit{G}\mathrm{\left(}{\mathit{p}}_{\mathit{\chi }}\mathrm{\right)}}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathrm{hs}\mathit{,}\mathrm{entro}}\\ & & \mathrm{+}\hspace{0.17em}\sqrt{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{K}\mathrm{(}{\mathit{p}}_{\mathit{\nu }}\mathrm{\left)}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{K}\mathrm{(}{\mathit{p}}_{\mathit{\chi }}\mathrm{\left)}\mathrm{\right)}}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathrm{lin}\mathit{,}\mathrm{entro}}\\ {\mathrm{\Gamma }}_{\mathrm{c}\mathit{,}\mathrm{baro}}& \mathrm{=}& \mathit{F}\mathrm{\left(}{\mathit{p}}_{\mathit{\nu }}\mathrm{\right)}\mathit{G}\mathrm{\left(}{\mathit{p}}_{\mathit{\nu }}\mathrm{\right)}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathrm{hs}\mathit{,}\mathrm{baro}}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{K}\mathrm{(}{\mathit{p}}_{\mathit{\nu }}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathrm{\Gamma }}_{\mathrm{lin}\mathit{,}\mathrm{baro}}\mathit{,}\end{array}$where the function F(p) governs saturation: $\mathit{F}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{(}\mathit{p}\mathit{/}\mathrm{1.3}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(17)and the functions G(p) and K(p) govern the cut-off at high viscosity: $\begin{array}{ccc}\mathit{G}\mathrm{\left(}\mathit{p}\mathrm{\right)}& \mathrm{=}& \{\begin{array}{c}\\ & & \\ & & \end{array}\\ \mathit{K}\mathrm{\left(}\mathit{p}\mathrm{\right)}& \mathrm{=}& \{\begin{array}{c}\\ & & \\ & & \end{array}\end{array}$with p_ν the saturation parameter related to viscosity and p_χ the saturation parameter associated with thermal diffusion: $\begin{array}{ccc}{\mathit{p}}_{\mathit{\nu }}& \mathrm{=}& \frac{\mathrm{2}}{\mathrm{3}}\hspace{0.17em}\sqrt{\frac{{\mathit{r}}_{\mathrm{P}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{P}}{\mathit{x}}_{\mathrm{s}}^{\mathrm{3}}}{\mathrm{2}\mathit{\pi }{\mathit{\nu }}_{\mathrm{P}}}}\\ {\mathit{p}}_{\mathit{\chi }}& \mathrm{=}& \sqrt{\frac{{\mathit{r}}_{\mathrm{P}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{P}}{\mathit{x}}_{\mathrm{s}}^{\mathrm{3}}}{\mathrm{2}\mathit{\pi }{\mathit{\chi }}_{\mathrm{P}}}}\mathit{,}\end{array}$ν_P the kinematic viscosity at the planet location, χ_P the thermal conductivity at the planet position, and x_s the half width of the horseshoe ${\mathit{x}}_{\mathrm{s}}\mathrm{=}\frac{\mathrm{1.1}}{{\mathit{\gamma }}_{\mathrm{eff}}^{\mathrm{1}\mathit{/}\mathrm{4}}}\hspace{0.17em}{\left(\frac{\mathrm{0.4}}{\mathit{ϵ}\mathit{/}\mathit{h}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\sqrt{\frac{\mathit{q}}{\mathit{h}}}$(22)where the smoothing length is ϵ/h = b/h_pr = 0.4.

The various contributions of the corotation torque are strongly sensitive to the temperature and surface mass density gradients. Both the corotation and Lindblad torques also scale with ${\mathit{M}}_{\mathrm{P}}^{\mathrm{2}}$ through Γ₀. The corotation torque also varies with the mass of the planet through the half-width of the horseshoe region.

Paardekooper & Papaloizou (2009a) show that the planet mass for which the corotation torque saturates is an increasing function of the viscosity. Actually, their Fig. 14 anticipates this mass to grow as ${\mathit{\nu }}_{\mathrm{visc}}^{\mathrm{2}\mathit{/}\mathrm{3}}$. Whereas Bitsch et al. (2015a) use a viscosity of α_visc = 0.0054, ours is α_visc = 0.01. Therefore, we should expect a slightly larger planet mass where the corotation torque saturates and cannot cancel the Lindblad torque in our case (around 1.5 times greater than in the Bitsch et al. 2015a, case).

3.3. Planetary trap evolution

In this paper, we successively model a variety of planet masses from 1 M_⊕ to 200 M_⊕ in interaction with an evolving protoplanetary disk, which is located between 0.1 and 30 au, to focus on the planetary formation region. The total torque exerted by the disk over the planet is given by ${\mathrm{\Gamma }}_{\mathrm{tot}}\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{Lindblad}}\mathrm{+}{\mathrm{\Gamma }}_{\mathrm{corotation}}\mathit{.}$(23)This total torque is positive when the disk yields angular momentum to the planet, making it migrate outward. A negative torque however results in the disk gaining angular moment from the planet, making the planet migrate inward. Lyra et al. (2010) described the zero-torque interface as an “equilibrium radius”. There are two different kinds: places of convergence called “traps” where the inner torque is positive and the outer is negative, and places of divergence, the planetary “deserts” (Hasegawa & Pudritz 2012 and BCP15), where the inner torque is negative and the outer positive. Planetary embryos accumulate in traps, while deserts are expected to be depleted in planets.

Assuming fixed surface mass density radial profiles, Hasegawa & Pudritz (2011b) estimated planet trap locations in steady-state disks based on Lindblad torques alone. Bitsch et al. (2014) proposed that Neptune like planets could get trapped in the regions of outward migration owing to shadowed regions, while BCP15 then studied dynamically evolving disks and showed that the traps are mostly located at the outer edge of the temperature plateaux, and the outer edge of the shadowed regions.

Figure 2 shows the evolution of the traps and deserts along the disk dynamical evolution for various planet masses from 1 M_⊕ to 200 M_⊕. In a typical MMSN-disk, density and temperature gradients are such that the total torque induces inward migration. Therefore, a change of sign of the total torque is generally accompanied by a second one to maintain inward migration in the inner and outer disk regions. Thus, there is always an even number of zero-torque lines with, for each couple, an inner desert and an outer trap.

For a low mass-planet (1 M_⊕), we notice a single “double inversion” located around 10−12 au, at the location of the heat transition barrier between an inner region that is dominated by viscous heating and an outer region being dominated by irradiation heating. BCP15 show that this barrier affects the temperature, density and geometrical profiles of the disk, and that it is a potential planet trap location. However, 2D-disk models from Bitsch et al. (2014, 2015a) do not show evidence of such a strong barrier in the temperature profile and, consequently, do not find any trap in the migration maps around that location. Although their vertical disk profile is initially passive, their simulations involve a proper treatment of the viscous heating. And, although their model accounts for a consistent geometry (including the shadowing effects), we believe that our approach is more self-consistent in relating the disk structure to its evolution: the disk temperature provides the viscosity necessary for calculating the temperature as well as the mass accretion rates across the disk that will, in turn, affect the density profile that constrains back the temperature profile. However, their simulations take into account the heat diffusion, which is not considered in the code presented in this paper and should be in the scope of a future work. This effect should definitely smooth the midplane temperature discontinuities at the transition. In a future paper, we will investigate the question of the sustainability of the heat-transition barrier trap with radial heat diffusion.

Between 5 M_⊕ and 20 M_⊕, we observe several “double inversions”. The outer one is now located a few au inside the heat transition barrier and is now slightly different to the ones observed for M_P = 1 M_⊕ and M_P ≥ 50 M_⊕. BCP15 describe in detail the importance of the temperature plateaux with regard to the positions of traps and deserts. These plateaux, owing to the sublimation of one of the dust’s major components, present the most brutal temperature gradient variations, along with the heat-transition barrier that separates the inner disk, which is dominated by viscous heating, and the outer one, which is dominated by irradiation heating. As the torque expressions show, temperature gradients are expected to have a stronger influence on the total torque than the density gradients. Therefore, total torque inversions are correlated to the edges of the temperature plateaux. As BCP15 state, we should refer to sublimation zones rather than lines: in this paper, we assume that the sublimation line just marks the middle of the so-called sublimation zone. Therefore, at the inner edge of this zone, where the temperature gradient is negative inside and almost zero outside, we expect a transition from a positive torque to a negative torque, resulting in a convergence location for planetary embryos (trap). Similarly, a divergence location (planetary desert) is expected at the outer edge of the temperature plateau, where the temperature will resume decreasing. Thus, traps and deserts are expected to be related to sublimation zones. This is consistent with the observation of traps located at the sublimation lines of the water ice, volatile organics, troilite, and between refractory organics and troilite; while deserts are found 1 au outside the water, volatile organics, and refractory organics ice lines, just inside the troilite line and at the silicates’ sublimation line. However, we notice that planets more massive than 10 M_⊕ are less inclined to be trapped at the volatile organics sublimation line.

In addition, intermediate mass planets (5 M_⊕ to 50 M_⊕) can be trapped around 0.3 au. These traps seem to be transient since they appear early in the disk evolution and disappear in less than 10 000 yr, a long time before the steady state is reached.

For planets more massive than 20 M_⊕, the silicate sublimation line no longer generates deserts after 10 000 yr. For more massive planets, traps and deserts associated with the sublimation lines are only possible in the early ages of the disk and progressively disappear completely for masses greater than 100 M_⊕. In this case, the only remaining trap and desert are located at the heat transition barrier.

Based on the relative abundances of the dust elements (Table 1) and on the opacity jumps associated with the sublimation of these elements (Fig. 1), we can expect the sublimation of the water ice (the major component), refractory organics, and silicates to affect the medium opacity in the most powerful way, while the opacity will be less sensitive to the troilite sublimation and even less to the volatile organics sublimation. The huge variation in the opacity within a very small range of temperature is reflected in the width of the corresponding temperature plateau, and therefore in the radial separation between the associated trap and desert. The amount of planetary embryos that are expected to accumulate at a given trap grows with the area of convergence towards that trap in the disk, i.e., with that radial separation. The most effective traps should therefore be the water ice trap, the refractory organics trap, and the silicates trap.

4. Discussion

4.1. Influence of the planet mass

Planetary growth influences the planet dynamics. In this section, we investigate the possibility of trapped planets remaining trapped while they grow. It is therefore necessary to understand how a change in the planetary mass affects the total torque experienced by the planet.

As stated in Eq. (6), the planet mass affects the torque amplitudes through the normalization coefficient Γ₀(r_P). It also changes the saturation of the corotation torques and, therefore, the total torque value. Section 4.3 details migration maps, which display the total torque amplitude as a function of the radial distance to the star and the planet mass. As these migration maps are normalized by Γ₀(r_P), the corotation torque saturation remains the only effect from the planet that may affect these maps; the dynamical evolution of the disk itself also obviously affects them as the temperature and density profile vary while the disk ages.

The total torque sign is affected by the variations in temperature gradients (more important than the density-gradient variations) at the sublimation temperature of the dust main components. However, these gradients remain negative around the sublimation regions, which is not necessarily the case around the heat transition barrier between viscous heating and irradiation heating domination, where the temperature gradient may become positive. In the case of a disk that is evolved from a typical MMSN, BCP15 show that, for a 10 M_⊕-planet in the inner disk (around the sublimation regions), a negative Lindblad torque was compensated for by a positive corotation torque so that the total torque could become positive itself; whereas near the heat transition, it was the temperature gradient that became positive, therefore reversing the Lindblad and the total torques. For lower mass planets, the mass is too low for the corotation torque to grow and balance the Lindblad torque.

Equations (15)−(22) show that the viscosity parameter p_ν increases with the planet mass following ${\mathit{p}}_{\mathit{\nu }}\mathrm{\propto }{\mathit{M}}_{\mathrm{P}}^{\mathrm{3}\mathit{/}\mathrm{4}}$. Paardekooper et al. (2011) mentioned that for low viscosity parameters (i.e., small planet masses), the corotation torque could possibly be lower than half its maximum value (corresponding to a stronger viscosity). However, as the planet mass increases, the corotation torque increases as well and is now able to compensate for the Lindblad torque. In addition, we notice from these equations that the corotation torque is more sensitive to the density gradient than the Lindblad torque is. As the planet mass gets significantly higher, the corotation torque drops: this cut-off at high viscosity is reflected by the G and K functions tending to 0.

This is consistent with what we observe in Fig. 2: for low-mass planets (e.g., M_P = 1 M_⊕) or too massive planets (e.g., M_P ≥ 100 M_⊕), the corotation torque does not allow the sign of the total torque to reverse and thus traps can only be found at the heat transition where the temperature gradient inversion alone is sufficient to reverse the Lindblad torque just outside the shadow region that is located around this transition.

As the torque exerted by the disk on a planet located in a trap is null, it is very likely that a trapped planet will remain trapped if the trap survives and the planet mass does not vary much. In that case, Fig. 2 shows that the planets slowly migrate towards the interior of the disk. However, for planetary masses M_P ≤ 50 M_⊕ located beyond 1 au, there can be several inner traps that would mean an escaping planet is trapped again.

4.2. Type II migration

The planet mass also favors the opening of a gas gap in the protoplanetary disk. In their Eq. (15), Crida et al. (2006) derived a criterion for opening a gap: $\frac{\mathrm{3}{\mathit{h}}_{\mathrm{pr}}}{\mathrm{4}{\mathit{R}}_{\mathrm{Hill}}}\mathrm{+}\frac{\mathrm{50}}{\mathit{q}\mathrm{ℛ}}\lesssim \mathrm{1}\mathit{,}$(24)where R_Hill is the Hill radius $\left({\mathit{R}}_{\mathrm{Hill}}\mathrm{=}{\mathit{r}}_{\mathrm{P}}{{}^{\mathrm{\right(}}\frac{{\mathit{M}}_{\mathrm{P}}}{{\mathit{M}}_{\mathrm{\ast }}}^{\mathrm{\left)}}}^{\mathrm{1}\mathit{/}\mathrm{3}}\right)$ and the Reynolds number is $\mathrm{ℛ}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{P}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathrm{P}}}{{\mathit{\nu }}_{\mathrm{P}}}$.

The first term varies according to $\mathrm{(}\mathit{h}\mathit{/}\mathit{r}\mathrm{)}{\mathit{M}}_{\mathrm{P}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}$ and the second one in $\mathrm{(}\mathit{h}\mathit{/}\mathit{r}{\mathrm{)}}^{\mathrm{2}}{\mathit{M}}_{\mathrm{P}}^{-1}$. Therefore, either an increase in the planet mass or a decrease in the local aspect ratio h/r may favor a gap opening. The aspect ratio profile from BCP15, Fig. 6 showed that the aspect ratio minima are located around 15 au in young disks and may get closer to almost 7 au when the disk ages. Thus, these are places where gaps are most likely to appear in our model if there were planets located in those locations.

In this paper, we consider the action of the disk over a putative planet without taking into account the back-reaction of the planet on the disk. In addition, we chose to consider a simplistic bimodal model for the disk-planet interaction: either the planet and the disk satisfy the criterion from Eq. (24) for opening a gap and the planet changes from type I to type II migration, or we consider it to be still embedded in a gas disk, not opening any gap, experiencing the torques defined by Paardekooper et al. (2011), and following type I migration. These torque formulas have been obtained for low-mass planets (for which the horseshoe width has to be lower than the disk scale-height). In our case, the aspect ratio profiles suggest that these torque formulas are valid for a planet up to a few tens of Earth-masses. Though our study extends beyond this mass limit, we consider these expressions as providing a first approximation of the torques experienced by a massive planet in our bimodal migration model. This could be improved by considering the partial opening of gaps in the disk by planets that are more massive than a few tens of Earth masses, by introducing a transition torque similar to Dittkrist et al. (2014) for example. Additionally, modeling the gap opening realistically would require taking into account the planet back-reaction on the disk, which could be done in the scope of a future work by coupling our hydrodynamical code with a planetary growth code.

4.3. Migration maps

In this section, we focus on the torque’s dependency on the planet mass and its distance to the star. To allow the torques to be compared, we display normalized torques by Γ₀ (Eq. (6)).

Figures 3−7 show the disk density and temperature profiles in parallel with the torque amplitudes at different radial and planet mass scales. In this distance-mass representation, we notice closed zones of outward migration in a continuum of inward migration. The contours of these zones correspond to the location where the total torque cancels. The outer borders are the planetary traps, while the inner ones are the planet deserts. As seen in the previous section, sublimation lines can be associated with the locations of the planet traps.

Fig. 3 Upper panel: midplane temperature (black), surface-mass density (blue), grazing angle (yellow) and viscous heating rate (red) radial profiles after 10 000 yr of evolution. Shadowed regions are displayed in gray and sublimation lines are shown in dotted lines. Lower panels: migration torque of a protoplanet with given radial distance to the central star rP and mass MP, in a protoplanetary disk after 10 000 yr of evolution. The white area verifies the criterion from Eq. (24). Black contours (0-torque contour) delimit the outward migration conditions while the rest of the migration map shows inward migration. Planetary traps are located at the outer edges of the black contours while planetary deserts are at the inner edges. The yellow dotted line marks the water ice line.

We note that there is an actual possibility for a 20 M_⊕-planet to get trapped at 12.5 au. And we also note that the outward migration zones correspond to the shadowed regions of the disk as indicated in gray in the upper panel of Figs. 3−7.

These outward migration zones associated with the ice lines are wider for lower-mass planets: as the planet mass increases, the trap is found further inside the disk, while the deserts remain at the same location. While we can see from Figs. 3−7 that these outward migration regions are actually shadowed, BC15’s Fig. 6 shows that the aspect ratio is locally decreasing across those regions. Equation (22) states that a decrease in the aspect ratio has a similar effect as a mass increase on the corotation torque saturation. Therefore the corotation torque saturates for lower planet mass in the outer region of the outward migration zone. However, for a disk of a given age, the maximal mass of the planets that can be trapped in the different traps is similar.

The most exterior outward migration zone is also the one that may trap the most massive planets. Though it is not associated with a sublimation line, it is related to the temperature and density profile irregularities that are generated by the heat transition region, as is visible in the upper panel. Furthermore, we note that below 50 M_⊕, the outermost trap is no longer related to the heat transition barrier, unlike for more massive planets: as we consider a lower planet mass, the corotation torque is no longer sustained and its saturation prevents the compensation of the Lindblad torque.

In addition, we note that there are zones in the migration maps (white areas) where the planets are massive enough and the aspect ratio low enough for a gap to open: the planet is now in type II migration. After 10 000 yr of evolution, planets more massive than 170 M_⊕ are able to open gaps around 15 au. Moreover, it appears that any planet located outside the most inner desert can escape type I inward migration: they migrate toward a trap and therefore avoid falling onto the central star. Though we consider the influence of the disk on putative planets of various masses, we do not actually follow their growth, migration, nor back-reaction over the disk.

Fig. 4 Disk profiles and migration maps for a disk after 100 000 yr of evolution. Legend is the same as in Fig. 3.

After 100 000 yr, it is no longer possible to trap planets more massive than 60 M_⊕ at the sublimation lines. That maximal mass is comparable with that of 38 M_⊕ from Bitsch et al. (2015a), Fig. 4, which was obtained for an accretion rate of 3.5 × 10^-8 M_⊙/ yr, and which is compatible with a disk that is slightly younger than 100 000 yr in our case. This maximal mass ratio is consistent with the 1.5 ratio expected in Sect. 3.2.1. We note that the outward migration zone beyond the water ice line for the low-mass planets is now split into two different outward migration zones: we explain the separation between the two zones by the presence of a shadowed region in the disk, which induces a change in the dominant heating processes in this place. In addition, the type II migration area disappears as the aspect ratio is higher at 100 000 yr than it was after 10 000 yr. It would require a more massive planet (M_P ≥ 200 M_⊕) to open a gap at this location.

Fig. 5 Disk profiles and migration maps for a disk after 1 million years of evolution. Legend is the same as in Fig. 3.

BCP15 indicate that after 1 million years of evolution, the minimal value of the aspect ratio was reached much further into the disk (around 7 au) than it was previously. This generates the white areas in Figs. 5−7 where potential planets would open gaps in the disk.

After that date, it becomes very difficult for a planet that is more massive than 40 M_⊕ to get trapped at a sublimation line. This upper limit decreases with time. BCP15 shows that, at a given radius, the pressure scale height decreases with time. The half-width of the horseshoe x_s then increases, leading to an increase in the viscosity parameter p_ν, which in turn induces a stronger saturation. Therefore, the planet mass at which the corotation torque cannot compensate the Lindblad torque decreases.

We also note the apparition of a few new outward migration zones that are not related to the sublimation lines: the increase in the number of shadowed regions leads to an alternation of dominant heating processes, which generate numerous consecutive outward migration zones for the most massive planets. This then helps to trap these massive planets (M_P ≥ 50 M_⊕) at locations further away from the sublimation lines.

Fig. 6 Disk profiles and migration maps for a disk after 5 million years of evolution. Legend is the same as in Fig. 3.

Fig. 7 Disk profiles and migration maps for a disk after 10 million years of evolution. Legend is the same as in Fig. 3.

These phenomena are even more obvious after 10 million years (Fig. 7). With time, the maximal mass of a planet that can be trapped decreases. This makes the trapping process more and more difficult as the planet mass grows over time.

5. Trapping growing planets

From Figs. 3 to 7, we can infer that a trapped planet will remain trapped as it grows. Indeed, for a 60 M_⊕-planet trapped at 10 au after 10 000 yr, if it were to grow, this planet will enter the inward migration zone that would bring it back inside the trap located just inside its initial position. Similarly, a 100 M_⊕-planet that is initially located at 15 au would enter the outward migration region if it grew. It would then migrate outward back to the trap, just outside its previous location. Thus, trapped planets will remain trapped and will follow the trap migration as they gain mass. In most cases, this involves a slow inward migration. If the planets reach the upper mass limit of the outward migration region where they are, gaining mass may free them and they will resume migrating inward. They may be trapped again when they cross the next inner outward migration region if their mass is still lower than the mass limit of that zone. Otherwise, they may not be able to escape the fatal inward migration. The planet growth speed is then the key parameter when estimating whether it will become massive enough to open a gap quickly or if it will not grow fast enough to avoid resuming a fatal inward migration. Since massive planets are already inclined to enter the runaway gas accretion phase, it is highly unlikely that they will remain trapped in the inner traps with a mass not exceeding 100 M_⊕. These types of planets would be saved if they do not pass that mass limit before photo-evaporation can dissipate enough gas to slow the planet type I migration.

Though a gap can be opened by a massive enough planet anywhere in the disk, we can infer a plausible scenario for gap-opening with a trapped minimal mass planet: a massive planet that follows a trap is likely to open a gap and enter type II migration at the radial distance at which the trap verifies the type II migration criterion from Eq. (24) (i.e., at the junction between the black trap contour and the white area in Figs. 3−7). After 10 000 yr, this could only happen around 15.8 au for planets more massive than 180 M_⊕. Similarly, after 1 million years, this can happen at 7.7 au for planets around 180 M_⊕. Even though it is very likely that the gas will be dissipated by photo-evaporation at longer evolution times, gaps could appear around 5 au for planets slightly more massive than 170 M_⊕.

Thus, we can list several types of planets that may escape type I migration: those that will grow very quickly and open a gap (with masses higher than 200 M_⊕), and those that will grow slowly and stay trapped until the gas dissipates and possibly avoid spiraling inward onto the star. The former may open gaps without necessarily experiencing type II migration rates: Dürmann & Kley (2015) show that in a disk where mass exceeds approximately 0.2 M_J, a Jupiter mass planet would migrate inward faster than if it was migrating at the viscous rate of the disk, and that the actual migration rate would depend on the disk surface mass density. Therefore, planets forming in the inner disk are still prone to fall onto the star: this is consistent with the results of planetary growth by pebble accretion from Bitsch et al. (2015b), which states that gas giants mainly form beyond 10 au (and up to 50 au if they form very early in the disk evolution) before migrating inward over several au. Therefore, investigating the survival of the giant planets requires a more complex model than the type I migration maps we present in this paper. The latter can be found at any distance from the star, although they can exceed 100 M_⊕ only if they get trapped in the outermost traps (those not associated with the sublimation lines). In addition, it seems possible to find planets trapped below 1 au.

6. Super-Earths

Exoplanet observations mentioned in Schneider et al. (2011) and Wright et al. (2011) show the existence of a population of planets with masses ranging between a few Earth masses and a few tens of Earth masses, and located inside 1 au from the star. A certain proportion of these super-Earths appear to be in mutual resonances (see Fig. 4a from Ogihara et al. 2015). Hansen & Murray (2012) and Hansen & Murray (2013) aimed at modeling the in situ formation of hot-Neptunes and super-Earths inside 1 au, using N-body numerical simulations. However, their simulations required a much more massive solid disk than the MMSN solid part. Although adding gas to such a massive disk would probably be gravitationally unstable, a more realistic disk model would consider a smaller quantity of gas by enriching the inner disk in solids by inward migration of dust grains, pebbles, and planetesimals. Their model however did not manage to reproduce the distribution of period ratios between adjacent planets.

Performing similar but more complete simulations, Ogihara et al. (2015) modeled super-Earths in situ formation starting with planetary embryos and planetesimals. They included planet-disk interactions, although with a simplistic disk model that assumes the temperature profile of an optically thin disk. They show that the interaction between the planets and the gas disk, and the resulting migration, should not be neglected. However, given their simplistic temperature profile, their disk cannot account for shadowing effects, temperature bumps, sublimation lines and, therefore, planetary traps. They claimed that the in situ formation of close-in super-Earths would require canceling type I migration inside 1 au.

Our disk structure shows that it is possible for planets of a few Earth masses to a few tens of Earth masses to get trapped at various locations inside 2 au. Thus, assuming that these planets formed either in situ (yet beyond 0.2 au) or in the outer disk, they may survive if they do not accrete too much mass. After 1 million years of evolution, we found traps located at 0.3 (for planets up to 21 M_⊕), 0.6 (up to 30 M_⊕), 1.3(up to 49 M_⊕), and 1.9 au (up to 40 M_⊕), which translates into periodicity ratios of 2.83 between the inner most and the second trap, and 1.77 between the third one and the fourth one (the others ranging between 3 and 16). After 5 million years, the inner traps are located at 0.3, 0.5, 0.9, and 1.45 au (for planets up to 20,25,44, and 34 M_⊕, respectively), which translates into periodicity ratios of 2.15 between the inner most and the second trap, 2.41 between the second and the third, and 2.04 between the third and the fourth.

Thus, periodicity ratios lower than 2 are consistent with the exoplanet observations that showed peaks at ratios of 2 and 5:3. Yet, we do not find any other correlation in periodicity ratios. Since our model only involves the viscous evolution of the disk and no real modeling of the individual planet dynamics, the mentioned correlation most certainly emerged by coincidence. However, the trap locations suggest that planets that would form further out in the disk and migrate towards these traps of periodicity ratios around 2, could be counted as in resonance.

Our simulations show that there is a possibility of retaining super-Earths (up to a few tens of M_⊕) that formed beyond 0.2 au in the simulations of Ogihara et al. (2015) by trapping them in our traps below 1 au. Though the second, third, and fourth traps can hold planets slightly more massive, the innermost trap can save planets up to 20 M_⊕, which is consistent with the fact that the vast majority of detected super-Earths that are below 20 M_⊕.

7. Conclusions and perspectives

7.1. Conclusions

Based on the viscous evolution of a MMSN, we have studied the possibility of trapping planets of various masses and preventing them from falling onto their host star. We have estimated the total torque exerted by the disk on a planet of a given mass at a given distance from the star. The mass of the planet mainly influences the saturation of the corotation torque, thereby affecting the positions of the planet traps and deserts across the disk.

For low mass (M_P< 5 M_⊕) or massive planets (M_P ≥ 100 M_⊕), there appears to be only one possibility of long-term trapping, located at the location of the heat transition barrier around 10 au (varying between 8 and 20 au with time). Only intermediate-mass planets (10 M_⊕ ≤ M_P ≤ 50 M_⊕) can get trapped at lower distance from the star: in these cases the traps usually coincide with the sublimation line of one of the major dust components (preferentially water ice, refractory organics, and silicates).

Trapped intermediate-mass planets are likely to remain trapped as they grow, until they reach the maximum mass of the outward migration zone where they are located. Then they might resume their inward migration, maybe get trapped to the next inner trap, and eventually fall onto the star. More massive trapped planets will also remain trapped in their respective traps slightly outer in the disk. If their mass and the local aspect ratio allow it, they might open a gap and then change from type I to type II migration. We estimate that massive planets are more likely to open gaps at certain trap locations, correlated with the edge of the shadowed regions: planets more massive than 170 M_⊕ can open a gap around 15.8 au after 10 000 yr of viscous evolution, or around 5.5 and 9 au after 5 million years.

We have identified several super-Earth traps in the disk, some of which are located around 5 and 7.5 au, though some others, associated with the sublimation lines, are located inside 2 au. These innermost traps would allow the survival of close-in super-Earths up to a few tens of Earth masses, such as the ones reported in exoplanet observations. They may also prevent the inward migration of super-Earths formed in situ below 1 au in the numerical simulations of Hansen & Murray (2012), Hansen & Murray (2013) and Ogihara et al. (2015). In addition, beyond the identification of resonance pairs in the planet periodicities, it would be very interesting to correlate these periodicity ratios with the trap locations.

Finally, observing young disks that possibly host planets and gaps (e.g., HL Tau), and coupling these observations with exoplanet detection and analysis would be very interesting to validate the giant planet-trapping at the heat transition barrier.

7.2. Perspectives

By coupling the disk evolution with planetary growth, we may be able to model planet dynamics, and then estimate at what distance the planet is more likely to get trapped and possibly open a gap. Knowing how the planet mass grows will determine its dynamics and its chances of survival. Modeling the planetary growth consistently with the disk evolution will therefore be the object of a future study.

In addition, it becomes essential to take into account type II migration since gap openings are going have a very strong effect on the gas density distribution in the disk.

Moreover, the presence of multiple planets located at the same trap line (outer edge of the same outward migration zone, though with various masses) will alter the disk by inducing planet-planet resonant interactions, as well as disk-planet resonances, as related in Cossou et al. (2013): planets can be trapped in chains of mean motion resonances. This requires a strong coupling between the hydrodynamical models of the disk evolution and the N-body models for planet interactions. It would be interesting to study the effect of this refinement on the counterreaction of the planet on the disk.

Although we used a detailed opacity model that accounts for the main dust components with a sublimation temperature above that of water ice, it might be interesting to add this opacity refinement for species with lower sublimation temperature to account for, for example, the presence of CO, CO₂, CH₄, and NH₃ ices.

Finally, recent studies by Fouchet et al. (2010) and Gonzalez et al. (2015), using bi-fluid SPH simulations, tend to show that the dust-to-gas ratio cannot be considered as constant and uniform in presence of planets. Adapting the opacity model to reflect dust-to-gas ratio variations would be the object of a future study.

Acknowledgments

We thank Dominic Macdonald for valuable suggestions that improved the quality of the manuscript significantly. We also thank the referee for detailed and constructive comments that improved the quality of the paper. This work was supported by IDEX Sorbonne Paris Cité and the Conseil Scientifique de l’Observatoire de Paris. We acknowledge the financial support from the UnivEarthS Labex program of Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02) and from the Institut Universitaire de France.

References

All Tables

All Figures

	Fig. 1 Mean-opacity variations with local temperature. Black: Rosseland mean opacity in extinction. Red: Planck mean opacity in absorption. Green: Planck mean opacity in extinction at stellar irradiation temperature. Blue: Planck mean opacity in absorption at stellar irradiation temperature. The corresponding dust-to-gas ratio is displayed with a dashed line.
In the text

	Fig. 2 Time evolution of the migration traps (blue +) and deserts (red ×) positions. The sublimation line positions and the heat transition radius are represented in black. Each subfigure shows the traps and deserts for a given planet mass.
In the text

Fig. 3 Upper panel: midplane temperature (black), surface-mass density (blue), grazing angle (yellow) and viscous heating rate (red) radial profiles after 10 000 yr of evolution. Shadowed regions are displayed in gray and sublimation lines are shown in dotted lines. Lower panels: migration torque of a protoplanet with given radial distance to the central star rP and mass MP, in a protoplanetary disk after 10 000 yr of evolution. The white area verifies the criterion from Eq. (24). Black contours (0-torque contour) delimit the outward migration conditions while the rest of the migration map shows inward migration. Planetary traps are located at the outer edges of the black contours while planetary deserts are at the inner edges. The yellow dotted line marks the water ice line.

In the text

	Fig. 4 Disk profiles and migration maps for a disk after 100 000 yr of evolution. Legend is the same as in Fig. 3.
In the text

	Fig. 5 Disk profiles and migration maps for a disk after 1 million years of evolution. Legend is the same as in Fig. 3.
In the text

	Fig. 6 Disk profiles and migration maps for a disk after 5 million years of evolution. Legend is the same as in Fig. 3.
In the text

	Fig. 7 Disk profiles and migration maps for a disk after 10 million years of evolution. Legend is the same as in Fig. 3.
In the text