© K. Baillié et al. 2019

1 Introduction

The huge diversity of observed exoplanets is a challenge for planetary formation scenarios, which must not only explain the trends in the distribution of the exoplanet orbital periods and masses but also retrieve the observational constraints of the solar system. The observations of protoplanetary disks also provide additional constraints necessary to model the environment and the conditions in which planets may form. This environment is the key element that determines which planets will fall onto their host star by spiraling inward by type-I migration in a few hundred thousand years (Goldreich & Tremaine 1979; Artymowicz 1993; Ward 1997) and which ones will avoid that by becoming trapped at the density gradient discontinuities (Masset et al. 2006; Paardekooper & Papaloizou 2009a), or at the opacity transitions (Menou & Goodman 2004). These latter transitions potentially result from the sublimation lines of the disk dust species (Baillié et al. 2015, 2016, hereafter referred to as BCP15 and BCP16). These planet traps also favor the accumulation of planetary embryos at certain radial distances from the star, where they may grow thanks to collisions (Morbidelli et al. 2008). Timescales are critical here since forming a planet by gas accretion on a solid core requires a few million years (Pollack et al. 1996), while the gas of the disk dissipates on a similar timescale (Font et al. 2004; Alexander & Armitage 2007, 2009; Owen et al. 2010), as confirmed by disk observations by Beckwith & Sargent (1996) and Hartmann et al. (1998).

Previous studies of planet migration mainly relied on simplified protoplanetary disk structures, mostly following power-law profiles as suggested by the classical minimum mass solar nebula (MMSN; Weidenschilling 1977; Hayashi 1981). Hasegawa & Pudritz (2011) and Paardekooper et al. (2011) used a similar framework to model the migration of planets in disks whose surface mass density and midplane temperature radial profiles followed power laws, while Bitsch & Kley (2011) and Bitsch et al. (2013) used density prescriptions to reconstruct (r, z)-temperature maps. Bitsch et al. (2013) allowed the geometry to be consistently calculated with the thermal structure.

In the present article, not only do we model the evolution of the disk structure (in density and temperature) by setting its geometry free (Baillié & Charnoz 2014b, hereafter BC14), but we also carry the initial state problem back from the usual MMSN to the primordial molecular cloud: we do so by modeling the formation of the disk itself by the gravitational collapse of the cloud, similarly to what Hueso & Guillot (2005) and Yang & Ciesla (2012) suggested, but this time accounting for the growth of the star and the evolution of its physical properties. The viscous evolution of the disk is described in BC14, BCP15 and BCP16 and takes into account the shadowed regions that are not irradiated directly by the star, the variations of the dust composition of the disk with the temperature, and the evolution of the dust-to-gas ratio. In addition, wenow take into consideration a proper model for the disk self-shadowing. We show that the disk tends to warm up during the collapse phase due to the stellar luminosity and that planet traps are carried away from the star, in particular at the sublimation lines (except for the silicates) and the heat transition barriers where the dominant heating process changes between viscous heating and stellar irradiation heating. We notice that a larger diversity of planets may be trapped than in the MMSN evolution models: in particular, these lines may hold Jupiter-like planets of a few hundred Earth masses located at a few astronomical units when the disk is so thick or viscous that this prevents gap opening. Finally, we highlight a migration mode dubbed “trapped migration” that allows planets to still migrate across the disk while remaining trapped, consistent with the results of Lyra et al. (2010).

Section 2 describes how we numerically model the star + disk system: from the growth of the star by the collapse of the molecular cloud that also feeds the disk, to the viscous spreading of the gas and dust. We highlight the improvements of the code detailed in BC14, BCP15 and BCP16 that account in particular for the early stellar evolution, the self-shadowing of the disk and the heating of the cloud. Based on the calculated evolution of the surface mass density and midplane temperature, Sect. 3 aims to estimate the influence of the disk structure on the protoplanet migration, and to determine the position of the planet traps. Section 4 investigates how naturally evolving disk structures may be related to the survival of growing planets. Finally, Sect. 5 summarizes our conclusions and details some perspectives that may provide a better understanding of the synthesis of the planet population.

2 Methods

2.1 Disk evolution model

Similar to BC14, BCP15 and BCP16, we model the protoplanetary disk as a viscous α-disk (Shakura & Sunyaev 1973) with a turbulent viscosity α_visc = 10⁻², as is usually taken for T Tauri-star disks without deadzones. Following Lynden-Bell & Pringle (1974) and Pringle (1981), we calculate the evolution of the disk surface mass density using the mass and angular momentum conservation. We assume that a cloud mass element d m_infall joins the disk at radius r over a time Δt, with an angular momentum dm_infall r² Ω(r). Due to the infall of material from the molecular cloud, the mass conservation reads $$\frac{\partial \Sigma (r,t)}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\Sigma (r,t)v_{\text{r}}(r,t)\right)=S(r,t),$$(1)

with Σ(r, t) being the surface mass density and v_r(r, t) being the radial velocity of the gas in the disk. The left part is similar to the expression derived by Pringle (1981) and the right part, S(r, t), is a source term that accounts for the infall of the molecular cloud gas onto the disk at radius r. This source term relates to the infall mass element by dm_infall = 2 π r Δr S(r, t) Δt.

The angular momentum conservation then reads $$\begin{array}{ll}r\frac{\partial }{\partial t}\hfill & \left(r^2\Omega \Sigma \right)+\frac{\partial }{\partial r}\left(r^2\Omega \cdot r\Sigma v_{\text{r}}\right)\hfill \\ \hfill & =\frac{\partial }{\partial r}\left(\nu \Sigma r^3\frac{\partial \Omega }{\partial r}\right)+rS(r,t)r^2\Omega ,\hfill \end{array}$$(2)

with ν being the viscosity. $$\begin{array}{l}{r}^3\Omega \frac{\partial \Sigma }{\partial t}+r\Sigma \frac{\partial }{\partial t}\left(r^2\Omega \right)+r^2\Omega \frac{\partial }{\partial r}\left(r\Sigma v_{\text{r}}\right)+r\Sigma v_{\text{r}}\frac{\partial }{\partial r}\left(r^2\Omega \right)\\ =\frac{\partial }{\partial r}\left(\nu \Sigma r^3\frac{\partial \Omega }{\partial r}\right)+rS(r,t)r^2\Omega .\end{array}$$(3)

While the term in $\frac{\partial }{\partial t}\left(r^2\Omega \right)$ can usually be forgotten when the star only accretes mass from the disk at limited mass accretion rates, we cannot neglect it in this study as the star is gaining mass from both the molecular cloud (during the collapse phase) and the disk. In addition, Eq. (1) allows to simplify the angular momentum equation: $$r\Sigma \frac{\partial }{\partial t}\left(r^2\Omega \right)+r\Sigma v_{\text{r}}\frac{\partial }{\partial r}\left(r^2\Omega \right)=\frac{\partial }{\partial r}\left(\nu \Sigma r^3\frac{\partial \Omega }{\partial r}\right).$$(4)

Assuming Ω ≈ Ω_Keplerian, we can write $$r\Sigma v_{\text{r}}=-3\sqrt{r}\frac{\partial }{\partial r}\left(\nu \Sigma \sqrt{r}\right)-r^2\Sigma \frac{{\dot{M}}_{*}}{M_{*}}.$$(5)

Finally, reporting this expression in Eq. (1), we find Eq. (6): $$\frac{\partial \Sigma (r,t)}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left(\sqrt{r}\frac{\partial }{\partial r}\left(\nu (r,t)\Sigma (r,t)\sqrt{r}\right)\right)+S(r,t)+S_2(r,t),$$(6)

where $$S_2(r,t)=\frac{1}{r}\frac{\partial }{\partial r}\left(r^2\Sigma (r,t)\frac{{\dot{M}}_{*}}{M_{*}}\right).$$(7)

Following Hueso & Guillot (2005) and Yang & Ciesla (2012) and their discussions about the chosen approximations, we model the cloud envelope as isothermal and spherically symmetric and assume that the presolar cloud rotates rigidly with a constant angular velocity and that the cloud material collapses below the centrifugal radius (defined as the point at which the maximal angular momentum in the shell is equal to the angular momentum in the disk): $$R_{\text{c}}(t)=10.5{\left(\frac{{\omega }_{\text{cd}}}{{10}^{-14}{\text{s}}^{-1}}\right)}^2{\left(\frac{T_{\text{cd}}}{15\text{K}}\right)}^{-4}{\left(\frac{M(t)}{1M_{\odot }}\right)}^3\text{AU},$$(8)

where M(t) is the total mass of the star + disk system at instant t, and ω_cd and T_cd are the rotational speed and temperature of the molecular cloud. The infall then happens in the inner parts of the disk, where r < R_c(t), and the source term is defined as $$S(r<R_{\text{c}}(t),t)=\frac{\dot{M}}{8\pi R_{\text{c}}^2}{\left(\frac{r}{R_{\text{c}}}\right)}^{-3/2}{\left[1-{\left(\frac{r}{R_{\text{c}}}\right)}^{1/2}\right]}^{-1/2},$$(9)

with Ṁ being the infall mass accretion rate from the cloud onto the star + disk system for which we use the expression of Shu (1977): $$\dot{M}=0.975\frac{c_{\text{S}}^3}{G},$$(10)

where $c_{\text{S}}=\sqrt{1.4\frac{k_{\text{B}}{T}_{\text{cd}}}{\overline{m}}}$ is the cloud isothermal sound speed, k_B is the Boltzmann constant and $\overline{m}=2.3m_{\text{H}}=3.8469{10}^{-27}$ kg, which is the mean molecular weight of the predominant H₂ gas in the disk.

In line with the works of Hueso & Guillot (2005) and Yang & Ciesla (2012), we consider that the isothermal cloud provides material to the disk at a constant rate throughout the entire collapse phase given by Eq. (10). Finally, we set the collapse phase to end when the total mass of the star + disk system reaches 1 M_⊙. The additional source term (Eq. (9)) is then canceled.

Among various attempts to model the collapse of the molecular cloud, Mellon & Li (2008) and Hennebelle & Fromang (2008) highlighted the limitations of the self-similar solution of Shu (1977) based on a purely hydrodynamical model: their MHD simulations showed that the collapsing material can drag magnetic field inward, intensifying its strength and leading to an outward transport of angular momentum. This phenomenon of “magnetic braking” detailed in Mellon & Li (2008) can slow down the gas rotation, possibly preventing the formation of a centrifugally supported disk. Comparing numerical models of Class-0 protostars with observations from IRAM Plateau de Bure, Maury et al. (2010) showed that observations are better explained by magnetized models of protostar formation than by a purely hydrodynamical simulation in the absence of turbulence. However, recent studies suggest a possible workaround to allow disk formation in a magnetic environment: while Dapp & Basu (2010) suggested a reduction of the magnetic field strength by Ohmic dissipation, Joos et al. (2012) investigated possible misalignments of the rotation axis and magnetic field direction, and Seifried et al. (2013) aimed at countering the magnetic braking by the disk turbulence. Li et al. (2014) provided an extensive review of the constraints and possible solutions for disk formation in the presence of a magnetic field, while Tobin et al. (2015) investigated the presence of such Class-0 disks in the Perseus region. Though our model is less realistic from a magnetic point of view, it still considers the magnetic field as being related to the turbulent viscosity α.

Masunaga et al. (1998) and Masunaga & Inutsuka (1999) showed that compressional heating might overtake radiative cooling and terminate the collapse phase. Though we believe that effect is of importance for the formation of the core, we chose to neglect that contribution for the disk, consistently with the models of Hueso & Guillot (2005) and Yang & Ciesla (2012). In addition, in our case, the termination of the collapse phase is determined by the initial mass of the cloud.

Similarly to BC14, BCP15 and BCP16, we divide the disk into consecutive annuli that are logarithmically distributed in radius between R_* and 1000 AU. We then model Eq. (6) numerically over a 1D grid of masses. The boundary conditions are set so that the disk cannot gain mass from the star at the inner edge. The mass accretion rate of the disk onto the star is therefore derived from the innermost mass flux.

For every radial position at every time, we model the thermodynamics and geometry of the disk by considering the heat equation described in Eqs. (15)–(18) from Calvet et al. (1991), accounting for stellar irradiation heating, radiative cooling, viscous heating defined as $F_v(r)=\frac{9}{4}\Sigma (r)\nu (r){\Omega }^2(r)$, and an external cloud envelope radiation heating $\sigma T_{\text{cd}}^4$, with σ being the Stefan–Boltzmann constant. Following the method thoroughly described in BC14, we solve this implicit equation numerically on the disk midplane temperature T(r). This allows us to derive the disk pressure scale height h_pr(r), its photosphere height H_ph(r), and its grazing angle α_gr(r) (the angle of incidence of the stellar irradiation upon the disk photosphere) defined as $${\alpha }_{\text{gr}}(r)=\text{arctan}\left(\frac{\text{d}{H}_{\text{ph}}}{\text{d}r}(r)\right)-\text{arctan}\left(\frac{H_{\text{ph}}(r)-0.4R_{*}}{r}\right).$$(11)

Though the self-shadowing of the disk was not taken into account in BCP16, we here consider that the irradiated parts of the disk cast shadows on the outer parts which have a lower incidence angle arctan(H_ph(r)∕r) than the shadowing inner parts. The value of the grazing angle determines whether a disk annulus is irradiated or not (in that case, we refer to “shadowed regions”). In addition, we consider the opacity model derived from Helling et al. (2000) and Semenov et al. (2003), and detailed in BCP16 (see their Fig. 1 and Table 1) to account for the composition of the dust species contained in the disk: water ice, volatile organics, refractory organics, troilite, olivine and pyroxene.

Since we aim at modeling the evolution of the disk over its lifetime and across wide radial scales, we are forced to neglect the heat diffusion. As we can see from the disks modeled by Bitsch et al. (2013, 2014, 2015a) that take into account the heat diffusion, we can expect the temperature structures to be smoothed radially. Therefore, the reader should keep in mind that our abrupt transitions (temperature plateau edges, heat transition barrier, optically thin Frontier) should probably be smoother in real disks.

2.2 Young star evolution model

Using the “PHYVE” code (Protoplanetary disk HYdrodynamical Viscous Evolution) thoroughly detailed in BC14, BCP15 and BCP16, we track the evolution of the mass distribution across the star + disk system: the viscous evolution governs the amount of material that is transferred from the disk to the star. In addition, the infall of mass onto the disk and onto the central star can be derived from Eqs. (6) and (9).

At every time step in the simulation, we interpolate the stellar radius R_*, luminosity ${\mathcal{L}}_{*}$ and temperature T_eff from tables of pre-calculated stellar evolutions modeled using the code CESAM detailed in Morel (1997), Morel & Lebreton (2008) and Piau et al. (2011), as well as Marques et al. (2013). These tables provide the radius, luminosity and temperatureof a constant-mass star as a function of its age and mass accretion rate. As these quantities are provided by the viscous evolution of the disk, we now have an empirical model for our star evolution, which is an interesting refinement compared to the fixed mass-accretion-rate star model used in Hueso & Guillot (2005). For the purpose of the present study, we make the approximation that there is no accretion luminosity associated with the stellar growth. This is indeed a necessary simplification in order to use the pre-calculated stellar evolution tables. However, we are confident that this approximation only marginally affects the disk evolution, which is the main focus of the present study. Indeed, this accretion luminosity would only be of the same order of magnitude as the stellar luminosity during the collapse phase. At this time, a larger luminosity would result in a hotter disk and therefore in a faster viscous spreading of the disk. Since this only takes place in the first few hundred thousand years, this will only temporarily accelerate the growth of the star, and is very unlikely to affect the later evolution of the disk and star.

2.3 Initial conditions

While previous works from BC14, BCP15 and BCP16 considered an MMSN around a classical T Tauri-type star as the initial condition of their simulations, in the present paper we relax that assumption and model the early evolution of the star and the disk while the collapse of the primordial molecular cloud supplies them with gas and dust. However, as we do not seek to model the star ignition, we assume that the star has already formed and grown up to 0.2 M_⊙ at the start of our simulation.

In addition, we chose the initial molecular cloud to be consistent with the initial clouds of Hueso & Guillot (2005) and Yang & Ciesla (2012). Following van Dishoeck et al. (1993), who estimated the temperature of the cloud falls between 10 and 20 K, and consistently with Yang & Ciesla (2012), we chose the temperature of the cloud to be T_cd = 15 K. We took its initial angular velocity ω_cd = 10⁻¹⁴ s⁻¹, in accordance with the observed velocity gradients in the clouds by Goodman et al. (1993), Barranco & Goodman (1998) and Lodato (2008) that provide a range from 10⁻¹⁵ to 10⁻¹³ s⁻¹.

3 Results

3.1 Star and disk growth

Though we do not directly model the evolution of the star, the interpolation of the star characteristics over pre-calculated stellar evolutions allows us to follow the evolution of its mass (Fig. 1), temperature, radius and luminosity (Fig. 2). The evolutionary track on the HR diagram is shown in Fig. 3, together with some pre-main sequence tracks at fixed mass. The star approximately follows the evolution of an accreting protostar with Ṁ = 10⁻⁵ M_⊙ yr⁻¹ (as seen in, e.g., Palla & Stahler 1990); however, the luminosity of the star is higher than an accreting protostar at the low-mass end (M_* < 0.4 M_⊙): the modeled star begins higher in the Hayashi track, meaning that its luminosity is closer to that of an accreting protostar later, at higher mass. This difference in luminosity, though not physical, is only temporary; it is just an effect of the initial conditions, in order to have a realistic luminosity in the mass range of interest. We plan to extend this study to the more realistic case where the accreting protostar and the disk are modeled in a self-consistent way.

Starting with an initial star + disk system of total mass 0.2 M_⊙, we reach a total mass of 1 M_⊙ after 170 000 yr of evolution. At this date, we set Ṁ to 0, but note that the star has not yet reached its final mass: it is only 0.84 M_⊙, and the disk mass is 0.16 M_⊙, meaning thedisk-to-star mass ratio is around 0.19. We verify that the Toomre instability can only appear transiently in the outer disk around the end of the collapse phase, and that it cannot affect the accretion rate of the disk onto the central star.

In the first phase (the gravitational collapse), the disk and star masses grow linearly. The disk gains mass from its surface while, at its inner edge, it loses some material that falls onto the star. After the cloud has emptied, it can no longer provide material and the disk can only yield gas and dust to the star by its inner edge: the disk mass slowly decreases while the star tends asymptotically toward 1 M_⊙ (see the bend in the evolutionary track in Fig. 3). After 1.64 Myr, the disk-to-star mass ratio is equal to 0.087, which corresponds to the one that BC14, BCP15 and BCP16 used as their initial disk following the MMSN (Weidenschilling 1977; Hayashi 1981).

From Fig. 1 and Eq. (8), one can directly estimate the evolution of the centrifugal radius which tends to 10.5 AU at the end of the collapse phase. In conditions similar to Yang & Ciesla (2012; see their Fig. 2), the mass evolution of the central star and disk are consistent with their results. During early evolution (before 10 000 yr), when the centrifugal radius is not significantly larger than the stellar radius, the collapsing material falls directly onto the star. This explains why the disk seems to start growing a little later than the star. Subsequently, when the centrifugal radius becomes larger, the disk receives more material and grows until the cloud gas reservoir is empty.

As a consequence of the observed mass accretion rates, the star also evolves: its characteristics vary during the collapse phase before stabilizing at 170 000 yr, as can be seen in Fig. 2.

The star temperature increases from 3400 to 4400 K, and its radius decreases from 6 R_⊙ to 1 R_⊙ after 10 Myr. The luminosity first decreases during the first 30 000 yr, before increasing up to $5{\mathcal{L}}_{\odot }$ at the end of the collapse phase. After about 2 Myr of evolution, the star is very similar to the one used in BC14, BCP15 and BCP16, though with a slightly lower radius than the 3 R_⊙ taken in these previous works. The final star corresponds to a classical T-Tauri pre-main sequence star, slightly fainter than its luminosity at the zero-age main sequence (ZAMS).

Fig. 1 Time evolution of the star and disk masses. The gravitational collapse that feeds the disk and the star ends after 170 000 yr.

Fig. 2 Time evolution of the star effective temperature, radius and luminosity.

3.2 Disk evolution

Photo-evaporation due to the intense stellar irradiation will dissipate the gas of the disk after a few million years (Font et al. 2004; Alexander & Armitage 2007, 2009 and Owen et al. 2010). However, as we do not take that effect into consideration in the scope of the present paper, we are able to pursue our simulations over 10 Myr. Outputs of our code should be considered with great care for the last few million years however, given that the modeled gas may have dissipated at that time.

BCP16 not only modeled the evolution of the surface mass density, but also derived the thermodynamical profiles of such disks: from the pressure scale-height to the surface geometry and the midplane temperature. In particular, they described how the sublimation of the various dust species across the disk generates a temperature plateau at the sublimation temperature of each of these species (corresponding to the drops in opacities as a function of the midplane temperature) and triggers shadowing of the disk surface. As the disk evolves and cools down, these plateaus drift inward; the silicate sublimation region soon being confined at the inner edge of the disk. In the present paper, we refer to the mean locations of these sublimation plateaus as “sublimation lines”.

Figure 4 shows the evolution of the radial profiles of surface mass density. Early in the disk lifetime, this profile cannot be simply represented by a power law, as was the case for the MMSN disk models. We notice significant changes in the slopes of those profiles until 1 Myr, after which the density profiles become smoother. As the disk gains mass and forms the molecular cloud, its surface mass density grows, especially in the inner disk. The inner viscous heating tends to increase as well as the disk optical thickness, which in turn affects the disk midplane temperature. As a further consequence, the vertical extent of the disk (pressure scale height or photosphere height), and therefore its geometry, is also affected by the infall of cloud material. At the beginning of the simulation (black curve at 100 yr), the disk is concentrated inside 1 AU. It then spreads as it gains in mass: after 100 000 yr, it reaches several tens of astronomical units with a density at 1 AU, exceeding the initial density of the MMSN disk in BCP15 at the same location by almost an order of magnitude. After 1 Myr (several hundred thousand years after the end of the collapse), the disk evolves viscously in a similar manner to the MMSN: it spreads and its mass decreases since the star is accreting the material that the disk yields at its inner edge. We notice that the asymptotic trend in the surface mass density profile becomes shallower than the MMSN as the disk ages. These decreasing power-law indices are indeed expected since the mass flux becomes increasingly uniform (Lynden-Bell & Pringle 1974; Chiang & Goldreich 1997; Bitsch et al. 2014 and BC14).

Figure 5 indicates a viscous mass accretion rate of approximately a few times 10⁻⁶ M_⊙ yr⁻¹ after 100 000 yr, a few times 10⁻⁸ M_⊙ yr⁻¹ after 1 Myr, and an order of magnitude lower after 10 Myr. Thus, mass accretion rates are above the estimated ones for the MMSN in the first 100 000 yr. After the collapse phase, such a disk has a mass accretion rate similar to that of an MMSN that would be younger by a few million years. In other words, the MMSN corresponds to a phase that occurs a few million years after the start of the collapse.

Similarly, the evolution of the temperature radial profile (Fig. 6) shows that the disk is initially cold and will be heated for a few hundred thousand years while it gains mass from the infall of the molecular cloud. Later, when the collapse phase is finished, the disk follows a similar evolution to that of an MMSN (BCP15) and cools. In particular, we notice that the temperatureplateaux (where dust species sublimate) drift outward at first and then inward after a few hundred thousand years.

During the collapse phase, it is worth noticing that the outer regions are the most affected by the cloud envelope radiation term: this imposes a minimal temperature of T_cd = 15 K anywhere in the disk and in particular in the outer regions. This results in a lower sensitivity of the outer disk temperature to the variations inthe disk optical depth in the regions where the disk is the thinnest.

Given that the density and temperature structures are implicitly related, it is very difficult to estimate a precise error bar on these profiles. However, as BCP15 showed by studying the impact of the radial resolution, density and temperature structures are still sustained despite resolution changes. Therefore, the sublimation zones and planet traps can only marginally be affected for a given set of initial parameters.

Figure 7 details the various heating contributions after 1 Myr. Irradiation heating locally dominates between 0.1 and 0.15 AU before viscous heating takes over when the disk gains optical depth. This defines an inner heat transition barrier, located slightly below 0.15 AU after 1 Myr. This coincides with an immediately outer shadow region located around 0.3 AU, shortly after the silicate sublimation plateau. The part of the disk inner to 0.1 AU appears to be consistent with the inner disk structures described by Flock et al. (2016) using radiation hydrodynamical simulations: they expect an inner rim around the silicate sublimation front that would cast a shadow upon the few tenths of astronomical units immediately further out. Flock et al. (2017) confirmed this in 3D radiation nonideal magnetohydrodynamical simulations. Though these latter studies treat MHD more thoroughly than we do here (they rely on the results of MHD magnetorotational turbulence models while we only consider it implicitly in the turbulent viscosity parameter), the inner disk regions present consistent geometric and thermic structures between MHD and purely hydrodynamical simulations. In addition, the grazing angle radial profile shows that the disk is only marginally irradiated below 10 AU at 1 Myr: the grazing angle hardly passes 0.1 rad. Around 10 AU, the viscous heating rate drops and the grazing angle grows to heat up the disk by a few Kelvins.

Figure 7 shows a succession of shadowed regions between 0.3 and 10 AU, which starts at the outer edge of the silicate sublimation plateau, similarly to what was observed in BCP15. Drops in dust-to-gas ratio are consequently expected at the temperature plateaus as dust elements sublimate species by species according to their sublimation temperature. Consistently with our opacity model, all the dust is sublimated for temperatures above 1500 K, leading to a dust-to-gas ratio of zero, whereas for temperatures below 160 K, all the dust species are solid and the dust-to-gas ratio reaches a maximum value of 1%.

In addition, after 1 Myr, we notice that the temperature structures are shifted outward: the outer edge of the plateau related to the sublimation of the silicates is now around 0.25 AU in the present simulation while it was located around 0.15 AU in the MMSN case. Likewise, the outer edge of the water ice line is now around 3 AU instead of 2 AU in the MMSN case. Nonetheless, the outer heat transition barrier (coinciding with the outer edge of the outermost shadowed region) is still found around 10 AU.

We may now follow the evolution of the positions of the sublimation plateaus: the water ice line at 160 ± 2 K (Fig. 8 upper panel) starts initially below 1 AU and moves up to 12 AU at the end of the collapse phase before drifting inward below 2 AU after 3 Myr. Apart from at the beginning of the collapse phase (the first few thousand years), the water ice line is always further away from the star than it was in the MMSN simulations. In addition, its width regularly exceeds 1 AU. The observation is similar for the silicate sublimation line at (1500 ± 20 K) that can be located around 2 AU at the end of the infall before tending toward 0.1 AU after 2 Myr.

Fig. 3 Evolutionary track on the HR diagram of the protostar at the center of the disk (brown dashed line). Full lines indicate classical evolutionary tracks at constant mass and the dotted line indicates the ZAMS.

Fig. 4 Surface mass density radial profile evolution of a disk formed by the collapse of a molecular cloud. The dashed line shows the MMSN defined by Hayashi (1981).

Fig. 5 Evolution of the viscous mass flux radial profiles for a disk formed by the collapse of a molecular cloud.

Fig. 6 Midplane temperature radial profile evolution of a disk formed by the collapse of a molecular cloud.

Fig. 7 Mid-plane temperature radial profile (black) after 1 Myr of evolution with a self-consistently calculated geometry and a full continuous model of opacities. Disk-shadowed regions are displayed in gray. The ratio of the viscous heating contributionto the total heating (viscous heating rate) is presented in red, the grazing angle radial profile in yellow and the optical depth radial profile in blue.

3.3 Equivalent timeline

Due to the difference in the initial conditions between the present simulations and those from BCP16, we have to consider that the initial ages of the disks are not equivalent in these two simulations. Admitting that the disk age can be correlated to its mass accretion rate onto its central star, we should compare disks evolved from different initial conditions at equivalent mass accretion rates rather than elapsed time in the simulation. To that end, we focus on the mass flux profiles across the disks presented in Fig. 5 of the present paper and Fig. 4 of BCP15. This suggests that a disk evolved from the collapse of the molecular cloud for 500 000 yr has a mass flux (and therefore a similar age) that is similar to that of a disk evolved from an MMSN for 10 000 yr. Similarly, a collapse-formed disk after 3 Myr seems equivalent in mass flux to an MMSN disk after 100 000 yr; a collapse-formed disk after 4 Myr would be equivalent to an MMSN disk after 1 Myr as shown in Fig. 9; and a collapse-formed disk after 7 Myr would be equivalent to an MMSN disk after 5 Myr. The density and temperature profiles (Figs. 4 and 6 of the present simulation vs. Figs. 3 and 5 of BCP15) tend to validate these associations, along with the evolutions of the sublimation lines (Fig. 8 of this paper vs. Figs. 11 and 12 of BCP15). This suggests that a disk evolved from an MMSN could coincide with a collapse-formed disk that would be older by at least 2 Myr.

In the present paper simulations, it takes a few hundred thousand years for the disk to reach a profile similar in shape to the ones evolvedfrom an already formed MMSN. However, comparing this disk with an early MMSN-like disk is not as pertinent as comparing disks that are close to reaching their steady state. Such a comparison shows that the MMSN simulation seems to follow the collapse simulation by a couple of million years.

Fig. 8 Upper panel: time evolution of the water ice sublimation region (mid-plane radial location for which the temperature coincides with the water–ice condensation temperature 160 ± 2 K). Lower panel: time evolution of the silicates sublimation zone (mid-plane radial location for which the temperature coincides with the silicate sublimation temperature 1500 ± 20 K).

3.4 Type-I migration

A planet excites resonances (Goldreich & Tremaine 1979; Ward 1988; Artymowicz 1993; Jang-Condell & Sasselov 2005) across the disk: Lindblad resonances caused by the action of the spiral arms induced by the planet, and corotation resonances due to the horseshoe region around the planet. Discounting the back-reaction of the planet on the disk structure, we can derive the resonant torques that the planet exerts on the disk and then calculate the reaction torque exerted by the disk on the planet. Using the expressions derived in Appendices A and B for the Lindblad and corotation torques, we derive the total torque that a hypothetical planet of mass M_P, located at a radial distance r_P from the central star, would receive from a viscously evolving disk. $${\Gamma }_{\text{tot}}={\Gamma }_{\text{Lindblad}}+{\Gamma }_{\text{corotation}}.$$(12)

Our model assumes a constant and uniform turbulent viscosity α, and does not consider any deadzone or cavity inside the disk: we may neglect variations in the viscosity and apply the torque formulas detailed in Appendices A and B. In addition, these torque expressions are consistent with the ones used in the previous works we compare to.

Since the Lindblad torque is in principle dominated by the variations of the temperature compared to the ones of the density, a disk for which the density and temperature profiles are decreasing with the distance to the star will experience a negative Lindblad torque, leading to a possible inward migration. However, the corotation torque may counter the Lindblad torque and slow down or reverse the migration. In the temperature plateaus, the temperature variation is much lower than the density variation, which could lead to a simplification of the expressions provided in the appendix. However, the strong influence of the cut-off and saturation effects makes it difficult to predict the final sign of the total torque. The edges of the regions where the migration is outward define locations where planetary embryos converge (the so-called “planet traps”), and places from which embryos run away (“planet deserts”). BCP15 showed that such planet traps usually appear at the outer edge of the sublimation regions of the different dust elements, or at the heat transition barrier between a region dominated by viscous heating and a region dominated by irradiation heating. This barrier coincides with the outer edge of the most external shadowed region (region for which the photosphere is no longer directly irradiated by the star).

The planet mass is another key parameter for the estimation of torque amplitude: the normalization coefficient Γ₀ (r_P) depends on the planet mass through Eq. (A.5), and furthermore the saturation of the corotation torques is a function of the planet mass. This second effect is the only one due to the planet that is visible in the migration maps of Sect. 3.7 which are normalized by Γ₀(r_P). We note that the evolving temperature and density of the disk also affect these maps as the disk ages. BCP16 stated that low-mass planets (< 10 M_⊕) cannot have a positive corotation torque large enough to balance the negative Lindblad torque near the sublimation lines, and therefore remain in inward migration. For more massive planets located near the sublimation lines, the positive corotation torque increases, and induces a positive total torque. Further out in the disk, near the heat transition barrier, a locally positive temperature gradient induces an outward migration for planets more massive than a few tens of Earth masses. The cut-off at high viscosity introduced by Paardekooper et al. (2011; Eqs. (B.9) and (B.10)) reflects a drop of the corotation torque for massive planets, as illustrated by the results of BCP16 showing that planet traps appear to be correlated with sublimation lines only for planet masses between 10 and 100 M_⊕. However, additional traps related to the heat transitions can be found for planets not in this range.

Fig. 9 Comparison of the disk radial profiles (upper panel: surface mass density, middle panel: temperature and lower panel: mass flux) for the disk at 4 Myr in the present paper (black line) and a 1 million-year-old disk evolved from an MMSN in BCP16 (red). We notice that the new treatment of the disk self-shadowing induces the presence of temperature drops of up to 10 K in the region 10–100 AU that possibly generate narrow heat transition barriers prone to trap planets.

3.5 Type-II migration

If the planet is massive enough and/or the disk aspect ratio h∕r is low enough, the planet may open a gap in the gas of the protoplanetary disk. Crida et al. (2006) derived an empirical criterion for the gap-opening condition: $$\frac{3h_{\text{pr}}}{4R_{\text{Hill}}}+\frac{50}{q\mathcal{R}}\lesssim 1,$$(13)

with $R_{\text{Hill}}=r_{\text{P}}{\left(\frac{M_{\text{P}}}{M_{*}}\right)}^{1/3}$, which is the Hill radius, $\mathcal{R}=\frac{r_{\text{P}}^2{\Omega }_{\text{P}}}{{\nu }_{\text{P}}}$, the Reynolds number, and $q=\frac{M_{\text{planet}}}{M_{*}}$, which is the mass ratio of the planet to the star.

BCP16 found that in the case of the viscous evolution of an MMSN, gaps may be opened at the places where the dominant heating process switches from viscous heating to irradiation heating: for instance, this concerned planets more massive than 170 M_⊕ around 15 AU in the early disk (10 000 yr), or planets above 180 M_⊕ located around 7 or 9 AU after 1 Myr. This could be explained by the trough in aspect ratio that BCP15 observed in their Fig. 6.

In the following sections, we study how planet traps and deserts, and also gas gaps, are affected by the disk origin (molecular cloud collapse) and how it may impact the growth scenarios of planetary embryos.

3.6 Planetary trap evolution

Based on the simulated disk profiles, we derive the torques that hypothetical planets of given masses and radial distances to the star would experience, and we extract the planet trap and planet desert positions for each of the given times. Focusing on the planetary formation region, Fig. 10 details the locations of these traps and deserts for planet masses between 1 and 200 M_⊕.

As described in BCP16, the radial profile of the total torque at a given time shows a background of inward migration with possibly a few intervals of outward migration. Therefore, we expect each desert to be accompanied by a trap slightly further out, as the total torque returns to a background negative value outside of this outward migration range.

A low-mass planet (1 M_⊕) can only be trapped at the outer heat transition Frontier, either during the first 1000 yr or after the end of the collapse. In the second case, such a small planetary embryo would be trapped at distances from the star of between 8 and 40 AU.

For planets of 5 M_⊕ or more, we notice a second recurrent population of desert + trap located at the inner heat transition Frontier around 0.2 AU that appears very early in the simulation (during the collapse phase). Thus, a planetary embryo with the mass of a super-Earth could get trapped on the orbits of Mercury or Venus, in the very early instants of a disk formed by gravitational collapse. In addition, new traps begin to appear at the sublimation lines of the dust species after 2 Myr.

More massive planets (10 M_⊕) present traps at the same sublimation lines as previously noted, though they appear earlier, after 400 000 yr. Trapping planets at the sublimation lines may happen earlier for more massive planets. Planets of 100 M_⊕ may even be trapped at all the sublimation lines, even during the collapse phase. We also note that the inner heat transition Frontier can only trap embryos of up to 100 M_⊕.

Planets more massive than 200 M_⊕ can no longer be trapped below 1 AU: the outer heat transition barrier remains the only possibility for these planets to be trapped after the collapse phase.

Comparing Fig. 10 with Fig. 2 from BCP16, we notice that traps now appear later than they used to in the case of the MMSN evolution. However, they seem to survive longer: for instance, a planet of 10 M_⊕ can be trapped after 400 000 yr while it was possible as early as 20 000 yr in the MMSN disk simulations. Similarly, planets of 100 M_⊕ could remain trapped at the sublimation lines until 10 000 yr in an MMSN disk while they can now survive the first million years in the case of a collapse-formed disk. This is consistent with the time delay observed in Sect. 3.3 between the timelines of the two simulations. Nevertheless, the present simulations provide indications of the trapping possibilities for planetary embryos during the phase when the collapse-formed disk cannot be modeled by a stage of the evolution of an MMSN. Thus, an early trapping at the inner heat transition Frontier enables the transient survival of more massive planets until this trap eventually falls down onto the central star at the end of the disk simulation. Finally, both simulations agree on the possibility of trapping (very) massive planets at the outer heat transition barrier.

Fig. 10 Time evolution of the migration trap (blue “+”) and desert (red “x”) positions. The sublimation line positions and the heat transition radius are represented with the black dotted and dashed lines. Each subfigure shows the traps and deserts for agiven planet mass.

3.7 Migration maps

In the present section, we present migration maps (Figs. 11 and 12) that display the normalized total torque exerted by the disk on a putative planet as a function of its radial distance r_P to the central star and its mass M_P. Torques are normalized by Γ₀ (Eq. (A.5)). In this distance–mass representation, the blue background stands for a negative total torque synonymous of an inward migration of the planet. Red zones on the contrary show closed regions of outward migration. The contours of these regions are the zero-torque radii: inner edges of those zones are the planet deserts while outer edges are the planet traps. Figure 11 summarizes this by indicating the direction of migration with yellow arrows. In addition, the white area is the region where Eq. (13) is verified, that is, where a planet is massive enough to open a gap and enter type-II migration. Finally, the water ice sublimation line is symbolized as a yellow dotted line.

As detailed in BCP16, a trapped planet gaining mass remains trapped at least until it reaches the upper mass limit of that trap. Beyond that limit, an embryo would resume migrating inward until it reaches an inner trap or falls onto the star.

Figure 11 shows different possibilities for trapping planetary embryos after 200 000 yr. First, we notice three traps located at the icelines: the troilite trap is located at 1.9 AU for 38 M_⊕ < M_P < 145 M_⊕ and the refractory organics trap is at 3.4 AU for 35 M_⊕ < M_P < 200 M_⊕. Outward migration zones due to the volatile organics and water ice lines merge to allow planet trapping between 5 and 8.5 AU for planets between 30 and more than 250 M_⊕. Another trap appears to be located further away, where the irradiation heating becomes dominant (11–45 AU). A very narrow fifth zone of outward migration can be spotted with a trapping possibility slightly below 0.2 AU for planets between 2 and 60 M_⊕, due to the inner heat transition Frontier as can be seen in Fig. 10. In addition, planets more massive than 300 M_⊕ (though very unlikely to exist at this stage) may open a gap below 0.3 AU, as well as planets more massive than 200 M_⊕ between 40 and at least 100 AU. These gaps cannot be accessed by a trapped planet as was the case for most of the type-II migration regions obtained by BCP16 with the MMSN disk model. However, Crida & Bitsch (2017) suggests that planets more massive than a few tens of Earth masses may undergo a runaway growth that may help these planets to reach the gap-opening area fast enough to avoid falling onto the star.

We also notice that our migration maps present more regions of outward migration than the work of Bitsch et al. (2015a) for instance. This can be explained by the fact that these regions are correlated with the sublimation lines of dust elements and that our opacity model involves more sublimation transitions than the Bell & Lin (1994) model they used in their work. In addition, the absence of heat diffusion in our code enhances the temperature gradients, consequently helping to build traps.

Figure 12 shows the traps after 4 Myr, that is, more than 3.8 Myr after the end of the collapse phase. Two groups of traps can be identified. The ones that are inner to the water ice line (yellow dashed line) are associated with some of the sublimation lines, while those that are outer (around 9.5 AU for planets of a few tens of Earth masses, and between 10 and 40 AU) are associated with the heat transitions and the subsequent exit of the disk self-shadow, as seen in Fig. 7. Traps at the volatile organics line (1.2 AU) and water ice line (1.7 AU) are only efficient on planets less massive than 47 M_⊕. At 0.65 AU (refractory organics line), only planets between 3 and 30 M_⊕ can be trapped, while at 0.35 AU (troilite line), planets of 4–25 M_⊕ may be saved. Moreover, it is possible at this date for a planet to open a gap and enter type-II migration anywhere assuming it is massive enough. This figure can be directly compared to Fig. 5 of BCP16 (MMSN after 1 Myr), according to our derivation of the equivalent timeline (Sect. 3.3).

Fig. 11 Migration torque of a protoplanet with given radial distance to the central star rP and mass MP, in a protoplanetary disk after 200 000 yr of evolution. Black contours (zero-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 blackcontours while planetary deserts are at their inner edges. The yellow dotted line marks the water ice line and the white area verifies the criterion from Eq. (13).

Fig. 12 Migration map after 4 Myr of evolution. The legend is the same as in Fig. 12.

4 Discussion

4.1 Protoplanetary disk formation

Most of the previous works that studied the numerical evolution of protoplanetary disks would initiate their numerical simulations with a disk already formed and generally following the MMSN model. This model can be criticized by questioning the positions of the planets at the time of their formation compared to their present positions (Crida 2009); here, we address this discrepancy by considering the initial parameters of the molecular cloud: its temperature, angular velocity and mass.

An extensive analysis of the infall parameters (total mass, mass accretion rate, temperature of the cloud and turbulent viscosity of the disk) could be the object of a full follow-up study. However, even though the locations of the temperature plateaus and shadow regions are likely to be shifted by varying the initial conditions, these structures will still appear in the vicinity of the sublimation lines of the major dust components in the disk. In addition, as mentioned in Baillié & Charnoz (2014a), an increase in α_visc will mainly accelerate the disk viscous evolution while it may allow planetary embryos to desaturate their corotation torque for larger masses, therefore allowing more massive planets to become trapped in more turbulent disks.

In addition, we consider now a joint evolution model for the disk and the star. Thus, the collapse phase lasts for 170 000 yr before the disk begins to spread viscously, similarly to the MMSN simulations. However, we retrieve the observational asymptotic trend of the MMSN after the first million years. Though the date at which the collapse ends is not affected by the approximation that neglects the accretion luminosity of the star, the stellar growth could be underestimated due to the lower temperature of the disk that slows down its viscous spreading in the absence of accretion luminosity. We estimate that the MMSN is similar to a disk that would have formed by collapse of a molecular cloud between 2 and 3 Myr ago. Therefore, the first millions of years of evolution after formation of a protoplanetary disk are determined by the collapse and cannot be modeled starting from an MMSN. However, 2 Myr after the beginning of the collapse, the later disk evolution could be approximated starting from an MMSN.

4.2 Favorable conditions for CAI formation

As we see in Sect. 3, the inner disk is now hotter, over a longer period. The midplane temperature at 0.2 AU now exceeds the silicate sublimation temperature for at least the time of the collapse phase. Moreover, at the same time, the mass flux is directed outward in the inner part of the disk (below 10 AU) until the end of the collapse phase. In the MMSN simulations, these conditions could only be found in the very inner regions of the disk during the first 10 000 yr.

In the case of a disk built from collapsing material from the molecular cloud, Yang & Ciesla (2012) showed that refractory materials such as calcium–aluminum-rich inclusions (CAIs) could form during the collapse phase, with a maximum around the end of the infall when the sublimation lines are the furthest away from the star. They would then be transported further out in the disk where they could be preserved in primitive bodies. In line with the model of this latter study, the disk structures presented here validate the two specific conditions for this process to happen: reaching the dust sublimation temperature in the inner disk and being able to transport the formed materials outwards, toward colder regions where they can recondense. This must happen on longer timescales and over wider regions of the disk than in the MMSN case. This should certainly provide favorable conditions for the models of CAI formation described in Charnoz et al. (2015).

4.3 Planet migration and survival

Migration maps (Figs. 11 and 12) show that after the collapse phase, any planet more massive than 2 M_⊕ will encounter a planet trap during its inward migration. After 4 Myr, the embryos that may fall onto the star are the ones that would not have reached 2 M_⊕ or that would have grown to over a few tens of Earth masses very rapidly. Planet growth is one of the critical process that may force a planet out of the traps.

Based on their locations, traps naturally select which planets they allow to escape type-I inward spiraling, with respect to their mass: different types of planets may be saved.

• The inner heat transition barrier may temporarily save super-Earths and hot Neptunes around 0.2 AU (the trap however does not seem to survive the end of the collapse phase).

• Planets, possibly including those that are quite massive, may be trapped at several locations corresponding to the sublimation lines below 10 AU for any evolution time. In particular, we identify trapping possibilities for gas and ice giants or Neptune-like planets within the few inner astronomical units in the earliest phases.

• After a few million years, the inner traps may save some close-in super-Earths (we find up to four traps inner to 1.8 AU after 4 Myr).

• Toward the end of the collapse phase, planets may be trapped further away in the disk at the outer heat transition barrier, regardless of their mass: though they are not required to open a gap there, these planets may still be saved as this heat transition trap appears to be sustainable until the end of the simulation. However, massive planets that have not yet opened a gap may be harder to save after a few million years if they have not been trapped earlier, unless they undergo runaway growth as suggested by Crida & Bitsch (2017).

• Gap-opening giant planets can enter type-II migration and eventually escape falling onto the central star by type-I migration, consistently with the expectations of Bitsch et al. (2015b) that gas giants formed by pebble accretion are born between 10 and 50 AU, and then migrate inward. However, there is no longer a trap leading to a gap-opening zone inside 10 AU.

BCP16 already noticed that trapping a planet at the silicate sublimation line is not possible. Though this is also valid here, we may still trap planets close to the star at the inner heat transition barrier.

Toward the end of the collapse phase at around 100 000 yr, planets as massive as 200 M_⊕ may be trapped at the sublimation lines, while this maximal mass drops to 60 M_⊕ at the water ice line when the disk has evolved for 1 Myr. While Bitsch et al. (2015a) showed favorable conditions for trapping planets of up to 50 M_⊕, one has to be careful when trying to compare their work with Figs. 11 and 12 here: migration maps should be compared only at similar mass accretion rates. As the model of Bitsch et al. (2015a) assumes that the disk is in a steady state with a set mass accretion rate, their figures can only be compared with ours when our disk shows a similar accretion rate. Thus, our 1-million-yr-old disk may be compared with their model using a mass accretion rate of 3.5 × 10⁻⁸ M_⊙ yr⁻¹ (their Fig. 4). In addition, the maximal masses of the planets that can become trapped have to be compared for a similar outward region (i.e., a similar ice line). Thus, our planets trapped at the water ice line appear to be slightly more massive than theirs (around 37 M_⊕). BCP16 mentioned that part of the difference can be explained by the fact that their model uses a turbulent viscosity α_visc = 0.0054, about half of the one used in the present work: this affects the planet mass at which the corotation torque saturates, and therefore the maximal mass of the trapped planets, leading to more massive planets in the present work. In addition, one must remember that the stars and heat diffusion treatments are different between these two models.

Planets of a few Earth masses that are trapped before the end of the collapse phase at the sublimation lines may survive until the end of the disk life while enduring a “trapped migration”, outward at first and then inward after the end of the collapse. Some of the less massive planets will survive at the outer heat transition barrier and undergo an outward “trapped migration” before resuming an inward type-I migration. If such a planet gains mass in that latest migration, it may become trapped again at a sublimation line.

Observations of multi-exoplanet systems (Fig. 4a of Ogihara et al. 2015) reported periodicity ratio peaks around rational numbers such as 2, 3∕2, 4∕3, 5∕4 and 6∕5, but also wider significant bumps around numbers apparently not related to main resonances, such as 1.8 or 2.2 for instance.After 1 Myr, approximating that they are not planet-mass-dependent and that the modeled disk parameters are representative of the disk population, we report planet traps at 0.55, 0.94, 1.7 and 2.3 AU, meaning that the period ratios of two trapped planets on different orbits would be in the range from 1.57 to 8.55. The ratio of 1.57is obtained for the two outermost traps and is consistent with the observed peak slightly above 3∕2. The two innermost traps present a ratio of 2.23, also close to a peak in the observed distribution. Similarly, after 4 Myr, planet traps are located at 0.36, 0.62, 1.12 and 1.65 AU, providing periodicity ratios between 1.79 and 9.81, with 1.79 being due to the two outermost traps and 2.26 owing to the two innermost traps. Again, these two pairs of traps could explain some of the observed nonresonantbumps in the distribution of period ratios of adjacent pairs of planets. In addition, we notice that at both dates, the remaining adjacent pair of traps present a period ratio of 2.43 that could also echo the slightly less important bump in Fig. 4a of Ogihara et al. (2015).

Steffen & Hwang (2015) also observed a peak in the period ratio distribution of Kepler’s candidate multiplanet systems around 2.2. Assuming a purely viscously heated steady state inner disk following a power law (equating $T_{\text{m}}{(r)}^4$ with the viscous heat flux defined in Sect. 2.1), we can derive that the midplane temperature follows T_m(r) ∝ r^−5∕6. From there, we can derive the period ratios of the locations associated with the temperature of the troilite and refractory organics sublimation lines: we evaluate it around 2.33, very close to the observed 2.2 peak. As this 2.2–2.3 ratio is recurrent in our simulations between the two innermost snow-line traps, this could constitute evidence that planets are saved at the snow lines, as suggested by Shvartzvald et al. (2016) who claim that 55% of microlensed stars host a snow-line planet.

Though we cannot assume that planets will remain at the trap positions after photo-evaporation and late planet migration, there is a possibility that multiple planets trapped at the sublimation lines around the same star could present period ratios close to the observed peaks.

4.4 Trapped migration

As stated previously and detailed in Lyra et al. (2010), a trapped planet is likely to remain trapped unless its growth makes it more massive than the maximum mass of the trap. However, these traps do not remain at the same location: not only do they tend to be slightly further in for more massive planets at a given disk age, but they also drift outward as the disk ages during the collapse phase and inward after the infall is finished. Therefore, a trapped planet can still undergo a slow migration that we call trapped migration. Before 170 000 yr, this trapped migration may push planetary embryos quite far out in the disk before they resume an inward migration: a planet at the water ice line could be driven as far as 12 AU before turning back for instance. Similarly, a planet trapped at the heat transition barrier could reach up to 35 AU. It is worth noting that this trapped migration exists independently of the formation of the disk by the collapse of the molecular cloud. However, it can push planets outward for a few hundred thousand years in our present simulations while trapped planets would just drift inward in the MMSN case.

5 Conclusions and perspectives

5.1 Conclusions

While previous works relied on an initial disk density profile following the controversial MMSN model, we here address that debate by bringing the initial condition back to the parameters of the molecular cloud at the origin of the star and disk. Therefore, we have a better estimation of the disk age than with the MMSN model, since we now model the disk formation by the collapse of the initial molecular cloud instead of assuming it has already formed and reached a power-law density distribution. In that “collapse” scenario, the disk now presents higher midplane temperatures that move the sublimation lines outward.

Another consequence of that disk-building scenario is that a great diversity of planets can now escape falling onto the star by type-I inward migration: very massive planets may open a gap and enter type-II migration before the end of the collapse phase, while low-mass planets can be trapped even in the very inner parts of the disk at the inner heat transition barrier. Jupiter-like embryos may also be trapped between 4 and 9 AU at 200 000 yr, as the heat-transition barriers on one hand, and the sublimation lines associated with dust main components on the other, still act as planet traps (except for the silicate line).

In addition, planetary embryos may undergo a trapped migration: they are dragged by the trap that moves outward during the collapse phase and inward afterwards. Such planets may survive below 1.8 AU after 4 Myr of evolution.

Our modeling of the molecular cloud collapse phase actually confirms that such formed protoplanetary disks can be approximated by evolved “MMSN disks” only a few million years after the collapse phase is over which in our case happens around 170 000 yr. In addition, this provides indications on the disk structure and the consequences on the planetary embryos during that formation phase: we can then model a much younger disk than was possible with the power-law structure model inherited from the MMSN concept. In particular, the question of the formation of the first solids at the earliest times of the protoplanetary disk requires that disk formation be taken into account and cannot simply be modeled as the evolution of an MMSN. From these comparisons, we could also establish a correspondence between the timelines of disks, depending on whether they form by cloud collapse or are already formed: the initial states of these two types of simulations appear to be separated by approximately 2 Myr. In the long term, this will certainly allow us to better adjust the birth lines of disks, stars, and early solids in protoplanetary disks.

5.2 Perspectives

Our model is now able to take into consideration the stellar evolution, the viscous spreading of the disk and the type-I migration: we have access to the optimal conditions for preventing the fall of planetary embryos and we can model the physical and chemical characteristics of those conditions. We now intend to consider a growth model (similar to Alibert et al. 2005; Coleman &Nelson 2014; Bitsch et al. 2015b) for these embryos to estimate the influence of the planet mass accretion rate on the likelihood of trapping (and therefore the likelihood of a planet surviving a fatal fall onto its star).

In addition, our model cannot currently take into account the counter reaction of the planets onto the disk, or the interactions between multiple planets. These are limitations of our model that should be addressed in order to model the observed resonant planet configurations.

Though we use a detailed opacity model that accounts for the dust-to-gas ratio, we do not consider the dust flux as possibly independent from the gas flux. For instance, extracting dust radial profiles from Gonzalez et al. (2017) to update our disk thermal structure might be useful.

Now that we have an evolving star model, it would be very interesting to properly account for photo-evaporation based on evolving stellar characteristics. This would help determine the planet distribution at the end of the disk phase.

Finally, we expressed the limits of our stellar model: the use of pre-calculated tables of constant mass star evolution could be improved to account for the accretion luminosity, especially during the collapse phase. One way to properly do that could be to couple a disk-evolution code with a stellar-growth code that would use the disk accretion rate onto the star as an input. Observations of accretion luminosity could then provide initial parameters for modeling the radial profile of potential disks that have not yet been characterized.

Acknowledgements

The author also thanks Dominic Macdonald for valuable suggestions that improved the quality of the manuscript significantly. This work was supported by the Conseil Scientifique de l’Observatoire de Paris and the Centre National d’Études Spatiales. We also thank the referees for detailed and constructive comments which improved the quality of the paper.

Appendix A Lindblad torques

Assuming the disk to be adiabatic and at thermal equilibrium, the total Lindblad torque exerted by a 2D laminar disk in the absence of self-gravity on a circular planet can be estimated using the formulas detailed in Paardekooper & Papaloizou (2008). As thermal diffusion slightly affects the wave propagation velocity, Paardekooper et al. (2011) defined an effective index γ_eff to account for the fact that this velocity is comprised between the isothermal sound speed (maximum thermal diffusion) and the adiabatic sound speed (no thermal diffusion). Therefore, this effective index replaces the γ-index previously used in the formulas of Paardekooper & Papaloizou (2008): $${\gamma }_{\text{eff}}=\frac{2Q\gamma }{\gamma Q+\frac{1}{2}\sqrt{2\sqrt{{({\gamma }^2{Q}^2+1)}^2+16Q^2(\gamma -1)}+2{\gamma }^2{Q}^2-2}},$$(A.1)

where Q accounts forthe thermal diffusion: $$Q=\frac{2{\chi }_{\text{P}}{r}_{\text{P}}}{3h_{\text{pr}}^3(r_{\text{P}}){\Omega }_{\text{P}}},$$(A.2)

and χ_P is the thermal conductivity at the planet location, $${\chi }_{\text{P}}=\frac{16\gamma (\gamma -1)\sigma T_{\text{P}}^4}{3{\kappa }_{\text{P}}{\rho }_{\text{P}}^2{h}_{\text{pr}}^2(r_{\text{P}}){\Omega }_{\text{P}}},$$(A.3)

with ρ_P is the density, κ_P the Rosseland mean opacity and Ω_P = Ω(r_P) the Keplerian angular velocity at the planet position in the disk. The 16 factor is a correction introduced by Bitsch & Kley (2011).

In the isothermal case, γ_eff = 1, whereas in the adiabatic case, γ_eff = γ = 1.4. The Lindblad torque subsequently becomes $${\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{Lindblad}}}{{\Gamma }_0(r_{\text{P}})}=-{\left(2.5-1.7\frac{\partial \ln T}{\partial \ln r}+0.1\frac{\partial \ln \Sigma }{\partial \ln r}\right)}_{r_{\text{P}}},$$(A.4)

with $q=\frac{M_{\text{planet}}}{M_{*}}$ which is the mass ratio of the planet to the star, $${\Gamma }_0(r_{\text{P}})={\left(\frac{q}{h}\right)}^2\Sigma (r_{\text{P}})r_{\text{P}}^4{\Omega }_{\text{P}}^2,$$(A.5)

and $h=\frac{h_{\text{pr}}(r_{\text{P}})}{r_{\text{P}}}$.

Appendix B Corotation torques

We estimate the corotation torque by considering the barotropic and entropic contributions separately, both of which contain linear and nonlinear parts. The barotropic was initially studied in Ward (1991) and Masset (2001) and detailed in Tanaka et al. (2002) as a torque arising from the density gradient, while the entropic corotation torque was expressed by Baruteau & Masset (2008) in the case of an adiabatic disk: the horseshoe region presents hotter material in the inner part than in the outer part, inducing, after a U-turn, an excess of mass leading the planet and a deficit of mass trailing behind the planet. This drives angular momentum exchange between the disk and the planet.

For low-enough viscosities (α_visc < 0.1), Paardekooper & Papaloizou (2009b) showed that the corotation torques are mostly non-linear, due 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.4 h_pr also used in Bitsch & Kley (2011) and Bitsch et al. (2014), the horseshoe drag torques read: $$\begin{array}{lll}{\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{hs},\text{entro}}}{{\Gamma }_0(r_{\text{P}})}\hfill & =\hfill & \frac{7.9}{{\gamma }_{\text{eff}}}{\left(-\frac{\partial \ln T}{\partial \ln r}+({\gamma }_{\text{eff}}-1)\frac{\partial \ln \Sigma }{\partial \ln r}\right)}_{r_{\text{P}}},\hfill \\ {\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{hs},\text{baro}}}{{\Gamma }_0(r_{\text{P}})}\hfill & =\hfill & 1.1{\left(\frac{\partial \ln \Sigma }{\partial \ln r}+\frac{3}{2}\right)}_{r_{\text{P}}}.\hfill \end{array}$$

In the general case (including possibly saturation), the total corotation torque is the sum of the barotropic and entropic contributions: $${\Gamma }_{\text{corotation}}={\Gamma }_{\text{c},\text{baro}}+{\Gamma }_{\text{c},\text{entro}},$$(B.3)

each of these contributions including a combination of nonlinear (Eqs. (B.1)–(B.2)) and linear parts (Eqs. (B.4)–(B.5)). The fully unsaturated linear expressions are $$\begin{array}{lll}{\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{lin},\text{entro}}}{{\Gamma }_0(r_{\text{P}})}\hfill & =\hfill & \left(2.2-\frac{1.4}{{\gamma }_{\text{eff}}}\right){\left(-\frac{\partial \ln T}{\partial \ln r}+({\gamma }_{\text{eff}}-1)\frac{\partial \ln \Sigma }{\partial \ln r}\right)}_{r_{\text{P}}},\hfill \\ {\gamma }_{\text{eff}}\frac{{\Gamma }_{\text{lin},\text{baro}}}{{\Gamma }_0(r_{\text{P}})}\hfill & =\hfill & 0.7{\left(\frac{\partial \ln \Sigma }{\partial \ln r}+\frac{3}{2}\right)}_{r_{\text{P}}}.\hfill \end{array}$$

Thus, Paardekooper et al. (2011) expressed the barotropic and entropic torques accounting for the saturation effects as follows. $$\begin{array}{lll}{\Gamma }_{\text{c},\text{entro}}\hfill & =\hfill & F(p_{\nu })F(p_{\chi })\sqrt{G(p_{\nu })G(p_{\chi })}{\Gamma }_{\text{hs},\text{entro}}\hfill \\ \hfill & \hfill & +\sqrt{(1-K(p_{\nu }))(1-K(p_{\chi }))}{\Gamma }_{\text{lin},\text{entro}},\hfill \\ {\Gamma }_{\text{c},\text{baro}}\hfill & =\hfill & F(p_{\nu })G(p_{\nu }){\Gamma }_{\text{hs},\text{baro}}+(1-K(p_{\nu })){\Gamma }_{\text{lin},\text{baro}},\hfill \end{array}$$

where the function F(p) governs the saturation $$F(p)=\frac{1}{1+{(p/1.3)}^2},$$(B.8)

and the functions G(p) and K(p) govern the cut-off at high viscosity: $$G(p)=\{\begin{array}{ccc}\frac{16}{25}{\left(\frac{45\pi }{8}\right)}^{3/4}{p}^{3/2}& \text{for}& p<\sqrt{\frac{8}{45\pi }},\\ 1-\frac{9}{25}{\left(\frac{8}{45\pi }\right)}^{4/3}{p}^{-8/3}& \text{for}& p\ge \sqrt{\frac{8}{45\pi }},\end{array}$$(B.9) $$K(p)=\{\begin{array}{ccc}\frac{16}{25}{\left(\frac{45\pi }{28}\right)}^{3/4}{p}^{3/2}& \text{for}& p<\sqrt{\frac{28}{45\pi }},\\ 1-\frac{9}{25}{\left(\frac{28}{45\pi }\right)}^{4/3}{p}^{-8/3}& \text{for}& p\ge \sqrt{\frac{28}{45\pi }},\end{array}$$(B.10)

with p_ν being the saturation parameter related to viscosity and p_χ the saturation parameter associated with thermal diffusion: $$\begin{array}{lll}{p}_{\nu }\hfill & =\hfill & \frac{2}{3}\sqrt{\frac{r_{\text{P}}^2{\Omega }_{\text{P}}{x}_{\text{s}}^3}{2\pi {\nu }_{\text{P}}}},\hfill \\ p_{\chi }\hfill & =\hfill & \sqrt{\frac{r_{\text{P}}^2{\Omega }_{\text{P}}{x}_{\text{s}}^3}{2\pi {\chi }_{\text{P}}}},\hfill \end{array}$$

where ν_P and χ_P are the kinematic viscosity and thermal conductivity at the planet position, and x_s is the half width of the horseshoe. $$x_{\text{s}}=\frac{1.1}{{\gamma }_{\text{eff}}^{1/4}}{\left(\frac{0.4}{\varepsilon /h}\right)}^{1/4}\sqrt{\frac{q}{h}},$$(B.13)

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 scale also with $M_{\text{P}}^2$ through Γ₀. The corotation torque also varies with the mass of the planet through the half-width of the horseshoe region; the latter determining the saturation and cut-off coefficients.

References

All Figures

	Fig. 1 Time evolution of the star and disk masses. The gravitational collapse that feeds the disk and the star ends after 170 000 yr.
In the text

	Fig. 2 Time evolution of the star effective temperature, radius and luminosity.
In the text

	Fig. 3 Evolutionary track on the HR diagram of the protostar at the center of the disk (brown dashed line). Full lines indicate classical evolutionary tracks at constant mass and the dotted line indicates the ZAMS.
In the text

	Fig. 4 Surface mass density radial profile evolution of a disk formed by the collapse of a molecular cloud. The dashed line shows the MMSN defined by Hayashi (1981).
In the text

	Fig. 5 Evolution of the viscous mass flux radial profiles for a disk formed by the collapse of a molecular cloud.
In the text

	Fig. 6 Midplane temperature radial profile evolution of a disk formed by the collapse of a molecular cloud.
In the text

Fig. 7 Mid-plane temperature radial profile (black) after 1 Myr of evolution with a self-consistently calculated geometry and a full continuous model of opacities. Disk-shadowed regions are displayed in gray. The ratio of the viscous heating contributionto the total heating (viscous heating rate) is presented in red, the grazing angle radial profile in yellow and the optical depth radial profile in blue.

In the text

	Fig. 8 Upper panel: time evolution of the water ice sublimation region (mid-plane radial location for which the temperature coincides with the water–ice condensation temperature 160 ± 2 K). Lower panel: time evolution of the silicates sublimation zone (mid-plane radial location for which the temperature coincides with the silicate sublimation temperature 1500 ± 20 K).
In the text

Fig. 9 Comparison of the disk radial profiles (upper panel: surface mass density, middle panel: temperature and lower panel: mass flux) for the disk at 4 Myr in the present paper (black line) and a 1 million-year-old disk evolved from an MMSN in BCP16 (red). We notice that the new treatment of the disk self-shadowing induces the presence of temperature drops of up to 10 K in the region 10–100 AU that possibly generate narrow heat transition barriers prone to trap planets.

In the text

	Fig. 10 Time evolution of the migration trap (blue “+”) and desert (red “x”) positions. The sublimation line positions and the heat transition radius are represented with the black dotted and dashed lines. Each subfigure shows the traps and deserts for agiven planet mass.
In the text

Fig. 11 Migration torque of a protoplanet with given radial distance to the central star rP and mass MP, in a protoplanetary disk after 200 000 yr of evolution. Black contours (zero-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 blackcontours while planetary deserts are at their inner edges. The yellow dotted line marks the water ice line and the white area verifies the criterion from Eq. (13).

In the text

	Fig. 12 Migration map after 4 Myr of evolution. The legend is the same as in Fig. 12.
In the text