Free Access
Issue
A&A
Volume 564, April 2014
Article Number A135
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201323007
Published online 18 April 2014

© ESO, 2014

1. Introduction

In recent years, it has been shown that the migration of low mass planets (1020ME) can be changed from inwards to outwards in discs with heating and cooling (Paardekooper & Mellema 2006; Baruteau & Masset 2008; Paardekooper & Papaloizou 2008; Paardekooper & Mellema 2008; Kley & Crida 2008; Kley et al. 2009; Ayliffe & Bate 2010). The heating in these previous works is provided by viscous friction, while for the cooling a local cooling rate (in 2D simulations) or a radiative diffusion (in 3D simulations) is utilised. This change of migration has important consequences for the formation of planets, as it would provide for a zero-migration radius in the discs, where planets can survive (Lyra et al. 2010; Bitsch & Kley 2011). Additionally to that, these points of convergent migration in the disc can trap planetary cores in resonances (Cossou et al. 2013). But these cores could then break free from the resonance (due to turbulence), collide and form the cores of giant planets (Pierens et al. 2013).

Only recently, Bitsch et al. (2013, hereafter Paper I) have shown the importance of stellar irradiation on the disc structure and the migration of low mass planets in numerical simulations. In Paper I, we did show that the change of the disc structure from a shadowed to a flaring disc greatly influences the migration of embedded bodies. Roughly speaking, a local increase in the aspect ratio H/r of the disc, results in inward migration, while a local decrease in H/r leads to outward migration. As a result, planets can migrate outwards in shadowed regions of the disc.

In Paper I, we pointed out that the opacity in the disc plays an important role in order to determine the disc structure, as the opacity is relevant for the absorption of stellar irradiation and for the cooling of the disc. We only considered equilibrium discs that have no net-mass flux through the disc. In this work, we want to expand to discs with a constant mass flow through the disc. In principle a viscous accretion discs is only accreting in the inner parts of the disc, while it is viscously spreading in the outer parts of the disc (Lynden-Bell & Pringle 1974) and actually behaves as a decretion disc there. If strongly depends on r, the surface density profile varies, so that the high regions empty, and the low regions fill, which smooths the variations of on a quicker time scale for sharper variations. Thus, it is reasonable to assume that is almost independent of r at a given moment in the inner parts of the disc, on which we focus here. This accretion rate can actually be measured in observable discs.

By studying different rates for discs, we explore the disc structures as a function of disc evolution, where decreases with time. Additionally, we study the influence of metallicity on the disc structure, as the dust-to-gas ratio can change in time as the disc evolves. The discs dust-to-gas ratio could decrease if the radial drift gets rid of the dust particles faster than the gas falls onto the star, or if planetesimals form efficiently, binding the dust. But it could increase as the gas is removed and dust is left behind or is regenerated (e.g. due to collisional grinding); in the extreme case this leads to debris discs, where the dust-to-gas ratio is infinite.

Additionally, we investigate the influence of the disc structure on the migration of embedded planets. As 3D simulations of planets are computationally very expensive, we relate back to use torque formulae for prescribing the expected migration of planets in the disc. We therefore are able to lay out a complete history of migration during the evolution of protoplanetary discs, which would not be possible in 3D.

The paper is structured as follows. First, we give in Sect. 2 an overview over the numerical methods used in this work. We then discuss about the disc structure and planetary migration in discs that have a metallicity of 0.01 in Sect. 3. The implications of different rates of metallicity on the disc structure and migration is described in Sect. 4. We then point out what roles viscosity and gas surface density play for the migration of planets in discs with constant rates in Sect. 5. In Sect. 6, we discuss the implications of our results on the formation and migration of planets in accretion discs. Finally, we summarize in Sect. 7.

2. Methods

The protoplanetary disc is treated in this study as a three-dimensional (3D) non-self-gravitating gas whose motion is described by the Navier-Stokes equations. Without any perturbers like planets, the disc is an axisymmetric structure. We therefore can use only one grid cell in azimuthal direction, making the computational problem de-facto 2D in radial and vertical (r-z plane) direction, where we utilize spherical coordinates (r-θ). The colatitude θ is measured in such a way that θ = 90° is the mid-plane of the disc. A clear picture of the source of turbulence inside accretion discs is highly under debate (see e.g. Turner et al. 2014). Nevertheless, some sort of viscosity is needed to drive the accretion disc, so we treat the viscosity utilizing an α-prescription (Shakura & Sunyaev 1973). The dissipative effects can then be described via the standard viscous stress-tensor approach (e.g. Mihalas & Weibel Mihalas 1984). We also include the irradiation from the central star, which was described in detail in Paper I. For that purpose we modified and substantially extended an existing multi-dimensional hydrodynamical code FARGOCA, as it is presented in Lega et al. (2014).

The radiative energy associated with viscous heating and stellar irradiation is then diffused through the disc and emitted from its surfaces. To describe this process we utilize the flux-limited diffusion approximation (FLD, Levermore & Pomraning 1981), an approximation which allows the transition from the optically thick mid-plane to the thin regions near the disc’s surface.

The hydrodynamical equations solved in the code have already been described in detail (Kley et al. 2009) and the two-temperature approach for the stellar irradiation was described in detail in Paper I, so we refrain from quoting it here again. As in Paper I, we take R = 3 R and T = 5600 K, which gives the flux from the star: (1)where σ is the Stefan Boltzmann constant and r the distance to the star. Stellar heating is responsible for keeping the disc flared in the outer parts of the disc (Paper I).

2.1. Opacity

The opacity is a crucial parameter when using an energy equation with stellar heating. The Planck mean opacity κP is used in the coupled 2-energy equations (radiative energy and thermal energy), the Rosseland mean opacity κR is used for the radiative cooling and the stellar opacity κ is used for the absorption of stellar photons in the upper layers of the disc. For more information on how the different opacities are connected to the energy equation, see Paper I.

thumbnail Fig. 1

Opacity profiles of the opacities used in the code. The different bumps in the profile correspond to opacity transition, where the ice line is located at 170 K. The grey area marks the region where the opacity changes due to the melting of ice grains. The maximum temperature during the simulation is reached in the inner parts of the disc, where viscous heating is dominant and can reach up to 900 K for the 1 × 10-7 M/yr case. Opacities are plotted for a metallicity of 0.01.

Open with DEXTER

In Fig. 1, the 3 different opacities used in the code are displayed. For the Rosseland mean opacity κR and for the Planck mean opacity κP, we use the same opacity law, as these opacities do not differ that much in the temperature region we are interested in (Paper I). However, compared to our previous work in Paper I, we changed the stellar opacity κ. In Paper I, we took the stellar opacity to be 0.1 of the Rosseland mean opacity, as indicated in Fig. 1. However, the stellar opacity depends on the temperature of the star T and not on the temperature of the gas, as dust grains absorb photons whatever their temperature is. As a consequence, κ is shown as a constant line in Fig. 1. The value of κ = 3.5 cm2/g corresponds to the typical value in a disc with a mixture of ice and silicate grains for T = 5600 K. The change of κ compared to Paper I will lead to a higher level of absorption of stellar photons in the upper layers of the disc, as the optical depth τ = κρΔr is much larger. Eventually, as the opacity was very low in the upper regions of the disc in the Paper I study, this will lead to flared discs for much lower gas densities than in Paper I. But at the same time, the main arguments about the disc structure (flared vs. non-flared) and for the migration still hold in the sense of the Paper I.

In the following, we define metallicity as the ratio of heavy elements to gas. In the condensed form, as ice or silicate grains, these heavy elements are, in our work, only of μm in size. We define ΣZ as the surface density of heavy elements in vapour or micrometer grains form and ΣG as the gas surface density. Thus, the metallicity Z is simply the ratio Z = ΣZ/ ΣG, assumed independent of r in the disc. This means that if grain growth would occur and the total amount of heavy elements (independent of size) stays the same, then the metallicity in our sense is still reduced.

2.2. Surface density

In a steady state accretion disc, the inward velocity of the gas is (2)where ν = αH2ΩK is the α-viscosity (Shakura & Sunyaev 1973) with H being the height of the disc and ΩK being the Keplerian orbital frequency. r denotes the orbital distance. With the radial velocity we can define an accretion rate as (3)where ΣG is the gas surface density. The slope of the surface density can be calculated using the viscosity (4)with h = H/r and where a = 9 / 7 as H/rr2 / 7 (Chiang & Goldreich 1997, Paper I). The proportionality, r2 / 7 describes the flaring index of the disc, which gives for a steady state accretion disc, where is constant for all radii ((d/dr) = 0) (5)where s denotes the power law index of the gas surface density in the disc. ΣG,0 is the value of the gas surface density at r0 = 1 AU. This leads to (6)We are using this slope of surface density for the initial conditions in all our simulations. Vertically, we initially set the density ρ of the disc to be in hydrostatical equilibrium (7)where ρ0(r) is the initial density.

For all simulations an adiabatic index of γ = 1.4 and a mean molecular weight of μ = 2.3 is set.

2.3. Boundary conditions

To get a disc with a given rate, we impose a -rate at the outer boundary, while we leave the inner boundary open. This way, the density structure inside the disc can change from the initial conditions and form bumps in the surface density that could not be created if a boundary was applied at the inner boundary as well, as the total mass inside the disc would stay constant and no re-filling could take place. The important quantities that have to be treated in a special way at the boundaries are the radial velocity vr and the gas surface density ΣG in the ghost rings. The details for the boundary conditions can be found in Appendix A.

2.4. Disc evolution

The initial conditions provided in Sect. 2.2 give an initial aspect ratio of the disc, (H/r)0, which has to be set to a higher value as the equilibrium H/r in order to get a flared disc, as a non flared disc cannot become flared any more (Dullemond 2002), see Fig. 2. In time then the viscosity would drop (as H drops and ν = αH2ΩK), until an equilibrium state is reached. Starting a simulation like that would imply that H and therefore Σ change a lot in time to guarantee a constant . So, if H drops by a factor of 23 as observed in Paper I, the surface density has to rise by a factor of 49. The time for the disc to adjust to that would be extremely long, as the disc has to be refilled from the outer boundary on a viscous timescale, which is essentially the life time of the disc.

We impose several steps in the simulations to save computation time. We first keep the viscosity constant in time, equal to its initial value , by using the initial h0 = (H/r)0 configuration to determine the viscosity during the whole first step. As the disc evolves, it compresses towards mid-plane and we reach an equilibrium state with a flared disc profile in the outer parts of the disc. We then fit this profile, which is labelled (H/r)1 in Fig. 2, with a new H/r value, where h1r2 / 7 (shown as (H/r)1,fit in Fig. 2). With this new estimate of H/r, we can compute a new α value for the viscosity by comparing it to (H/r)0. The change of H in the viscosity is compensated by the change of α, so that the disc will have the same viscosity as in the first step and therefore the rate is conserved as well.

The simulation is then restarted with the new α value (which is given by ) and the viscosity ν is now given by Eq. (4), with the local value of H to account for bumps and dips that might eventually form in the disc. For the outer parts, these changes of ν are not important as the disc structure is determined by stellar irradiation and not by viscous heating. In the inner parts, where viscous heating is the dominant heat source, the viscosity changes compared to the (H/r)0 case, leading to change of the disc structure. However, these changes are minimal. The advantage of this procedure is that the rate is the same as in the initial simulation with no need to refill ΣG as H drops.

thumbnail Fig. 2

Evolution of H/r in the disc for the different stages. The disc features = 2 × 10-7 M/yr. At the beginning, in the (H/r)0 stage, α0 = 0.001, which we change to α1 = 0.003 after the stage (H/r)1 is reached. The stage (H/r)2 shows the aspect ratio after the viscosity has change and (H/r)final shows the final stage after the surface density has been re-arranged.

Open with DEXTER

We then reach a new equilibrium state (H/r)2, where the inner parts of the disc changed compared to (H/r)1, as the viscosity changed. But, as the disc still needs to re-arrange the surface density on a viscous time scale, we relate back to a locally isothermal simulation, where we plug in the disc structure and sound speed as it is shown in the (H/r)2 disc. The disc is then evolved until a time independent mass flux through the disc is achieved. After that, we go back to the fully radiative energy equation with stellar irradiation to account for changes in the heating due to the re-arranged gas surface density ΣG and arrive at (H/r)final, which gives the final state of the disc. Please note that this whole process can only work as the outer disc is supported by stellar irradiation, keeping the input of mass at the outer boundary constant, not only in time, but also in vertical spacing as H stays constant.

Table 1

Simulation parameters for the used models.

In this way, much computation time is saved, as the computation time is mainly dominated by the stellar-irradiation algorithm that can take up to 90% of the whole computation time. Especially in the state where the surface-density has to be re-arranged, due to changes in viscosity this makes a big difference.

In principle during time the metallicity of small dust grains in the disc can change (e.g. due to grow of planetesimals), which would change the disc structure (see Sect. 4). This would imply that disc needs to re-arrange the gas inside, which happens on a viscous timescale. However, as we show below, this re-arranging of mass is much less pronounced in the outer parts of the disc, where the viscous timescale is much longer compared to the inner parts (τvisc = r2/ν = r2/ (αH2ΩK)). This is because the change of viscosity in the outer parts is much smaller compared to the change in the inner parts of the disc, as the H/r profile stays constant as it is dominated by stellar irradiation. This indicates that the steady-state of constant can be achieved in the inner parts on a timescale which is shorter than the lifetime of the disc, which is determined by the viscous timescale of the outer parts of the disc.

2.5. Numerical setup

We perform several simulations for different rates corresponding to different stages of the discs lifetime that are listed in Table 1. In the first sets of simulations, we assume that the dust is bound perfectly to the gas and as the gas accretes onto the star, the dust decreases at the same rate, keeping the dust to gas ratio constant. In the second set of simulations, we assume that the dust is not perfectly bound to the gas at the later stages of the discs evolution, resulting in an increasing and decreasing dust to gas ratio, depending on the favoured scenario for the evolution of dust particles.

All simulations are performed on a grid with 386 × 66 active cells in rθ direction, where the opening angle of the numerical simulation is 20°, meaning that 70°<θ< 90°. The radial domain extends from 1 AU to 50 AU.

3. -discs with constant gas to dust ratio

In this section we analyse the disc structure for discs with different rates, but with a constant metallicity of 0.01. The different rates correspond to different evolutionary states of the discs lifetime. We reduce the surface density of the gas, ΣG, to realize different states of the evolution, but keep the viscosity at the same value for all simulations. We then analyse the migration behaviour of planetary cores inside the discs by using the torque formula provided by Paardekooper et al. (2011).

3.1. Disc structure

In Fig. 3 the H/r, temperature and gas surface density profiles of discs with different are displayed. The H/r profiles show a flaring part in the outer disc and some bumps and dips in the inner part of the disc. These bumps and dips are caused by transitions in the opacity of the disc, for example at the ice line. This kind of structure was also observed for equilibrium discs with constant viscosity in Paper I. The height of the bumps and the depth of the dips in the disc is decreasing with decreasing rate (which means a decreasing ΣG). As decreases, so does the viscous heating, hence the temperature and H/r in the inner part of the disc. This also means that the shadowed region of the disc is becoming smaller as the discs density reduces. Additionally to that, the bumps of the disc structure move towards the central star as decreases. The bumps in the disc originate from the opacity transition, which always correspond to the same temperature region. Therefore less heating in the inner parts of the disc moves these bumps further inside.

As the gas surface density is reduced more and more, the bumps in the disc become smaller and smaller (as the disc heats less) until they finally disappear completely for small . For the 1 × 10-8 M/yr disc there is only a very small bump in the H/r profile visible. This bump vanishes for the 5 × 10-9 M/yr disc, which only shows an increase of H/r. In the outer parts of the disc, the structure is very similar for all . The heating in the outer disc is dominated by stellar irradiation, which is dependent on the optical depth (Δτ = ρGκΔr) of the disc. So, as long as the disc is optically thick (Δτ> 1) it can efficiently absorb stellar irradiation and be heated. The disc is then flared with H/rr2 / 7 as found for theoretical calculations in Chiang & Goldreich (1997).

thumbnail Fig. 3

H/r-profile (top), temperature (middle) and surface density profile (bottom) of discs with different rates. The profile is taken when the discs have reached their radially constant rate. The black dotted lines represent power-law fits to guide the eyes. The grey area in the temperature plot marks the temperature range where the opacity profile changes due to the melting of ice grains (see Fig. 1). Please note that we cut the displayed disc at 30 AU to enhance the details in the inner parts of the disc.

Open with DEXTER

In hydrostatic equilibrium, the temperature is related in mid-plane to H/r through (8)where μ is the mean molecular weight and the gas constant. The temperature profile (middle in Fig. 3) therefore reflects the fluctuations in the H/r profile by showing steeper gradients when H/r drops and less steep gradients when H/r increases. These changes are caused by the transition of opacity. The region of change in the opacity is marked in grey colour in the middle panel of Fig. 3.

The surface density is reduced by the appropriate factor for changing the rate. In the outer parts of the disc, this can well be observed (bottom in Fig. 3). However, in the inner parts of the disc this seems not to be exactly the case. For example between the 1 × 10-8 M/yr disc and the 5 × 10-8 M/yr disc, the difference in surface density at the inner boundary is only a factor of 2 instead of 5. This can be explained by the different viscosity in the inner parts of the disc, as ν = αH2Ω and H clearly varies for the different values. A smaller viscosity caused by a smaller H, results in a higher surface density as has to be constant, thus explaining the surface density profiles.

Additionally, there are wiggles in the ΣG profile at 5 AU for all discs. This also can be explained by changes in the viscosity due to changes in H (see top of Fig. 3) and that has to be constant throughout the disc. The changes of the gradients in the surface density ΣG and in the temperature T have important consequences for the migration inside these discs.

3.2. Migration maps

In order to estimate the torque acting on planets embedded in the discs with different , we use the formula by Paardekooper et al. (2011), which captures the effects of torque saturation in contrast to Paardekooper et al. (2010), where the torques are fully unsaturated. However, the formula by Paardekooper et al. (2011) might not be accurate for small mass planets. In fact, Lega et al. (2014) have found a new additional negative torque, which diminishes the total torque acting on embedded small mass planets. But these differences do not seem to be that large. So in absence of a formula that would capture this effect and to avoid large scale numerical simulations, we use the formula of Paardekooper et al. (2011), which gives a good approximation. The formula captures the torque caused by Lindblad resonances and horseshoe drag on low-mass planets embedded in gaseous discs in the presence of viscous heating and thermal diffusion. The formula does not include stellar heating, but stellar heating only changes the disc structure and not the mechanism responsible for outward migration (Paper I). This formula was also tested against 3D simulations in Bitsch & Kley (2011), who found a good agreement for 20 MEarth planets.

The formula of Paardekooper et al. (2011) is very complex, so we will not explain the whole formula here. The total torque acting on an embedded planet is a composition of its Lindblad torque and its corotation torque: (9)The Lindblad torque depends on the gradients of temperature Trβ and gas surface density ΣGrs. It is given in Paardekooper et al. (2011) by (10)where q is the mass ratio between planet and star, ΣP the gas surface density of the disc at the planet’s location and rP the distance of the planet to the host star. One can clearly see that a change in the gradient of temperature influences the Lindblad torque. The same applies to the corotation torque, which is strongly dependent on the gradient of entropy, Srξ, with ξ = β − (γ − 1.0)s. The largest contribution of the corotation torque arises from the entropy related horseshoe drag, which is given by (11)This indicates that we expect a change of the migration rate, when the temperature changes significantly. This is the case in the inner parts of the discs, where the aspect ratio H/r shows fluctuations (as H/r is proportional to T). In fact, in Paper I we found that outward migration seems only possible in regions of the disc where H/r drops. The torque formula that accounts for saturation effects (Paardekooper et al. 2011) will result in a smaller region of outward migration compared to the unsaturated torque formula (Paardekooper et al. 2010), which can clearly be seen in Fig. 4. Additionally the formula for the unsaturated torques is independent on the planetary mass, so that outward migration will be found for all planetary masses for the Paardekooper et al. (2010) formula, in contrast to the torque formula that accounts for saturation, as torque saturation depends on the planetary mass (Paardekooper et al. 2011).

Compared to the equilibrium discs with constant viscosity presented in Paper I that had seq = 0.5, we now have a much steeper surface density gradient with s ≈ 15 / 14. This reduces not only the entropy gradient in the disc and therefore the entropy related horseshoe drag, but also the barotropic part of the horseshoe drag, which is given by (12)by 50% compared to equilibrium discs with seq = 0.5. The effect on the Lindblad torque (Eq. (10)), however, is minimal. The expected total effect would be a reduced region of outward migration in the accreting discs compared to the equilibrium discs.

thumbnail Fig. 4

Torque acting on discs with different , with 5 × 10-8 M/yr (top), 1 × 10-8 M/yr (middle) and 5 × 10-9 M/yr (bottom) for the applied Paardekooper et al. (2011) formula. The black lines encircle the regions of outward migration for the Paardekooper et al. (2011) formula. The vertical solid yellow lines mark the outer edge of the zero-migration region for the unsaturated torques (Paardekooper et al. 2010), while the dashed-yellow line marks the inner edge of the zero migration region of the unsaturated torques. The vertical red lines indicate the ice line at 170 K.

Open with DEXTER

In Fig. 4 we present the migration maps for discs with = 5 × 10-8 M/yr, = 1 × 10-8 M/yr and = 5 × 10-9 M/yr, for which the disc structures are displayed in Fig. 3. For the highest disc model, we have two separated regions of outward migration, which are encircled by black lines in the figure. The first region is from the inner edge of the disc up to 3 AU, the second is from 5 AU <rP< 7 AU. These regions compare nicely to regions in the disc, where H/r drops. A drop of H/r corresponds to an increase of the temperature gradient resulting in an increased entropy gradient and hence a stronger entropy related corotation torque leading to outward migration. These changes in H/r are caused by opacity transitions, as was stated for equilibrium discs in Paper I. In fact, in the temperature profile shown in Fig. 3 (middle), the grey area marks the region where the opacity changes, which shows a shallower temperature gradient. For the = 5 × 10-8 M/yr disc this grey region (from 3 AU to 5 AU) corresponds nicely to the region of inward migration shown in Fig. 4.

Compared to the studies of equilibrium discs with constant viscosity (as in Paper I) where the minimum mass for outward migration was 5ME, it is now 8ME. This difference has its origin in the different disc structure. As the disc features a much steeper surface density gradient s ≈ 15 / 14 compared to the equilibrium disc seq = 0.5, it reduces the different positive contributions of the torque, so that very small mass planets do not migrate outwards any more. For small planetary masses, the horseshoe drag tends towards the linear corotation torque, resulting in inward migration. Infact, for small mass planets (<5 MEarth), new studies have shown an additional negative torque (Lega et al. 2014). The maximum mass of outward migration, 40ME seems to be the same compared to Paper I. However, one should be sceptical about the results for high mass planets, as the formula of Paardekooper et al. (2011) was derived in the linear regime for small mass planets and not for large planets that can even start to open up partial gaps.

We do not show the migration map for = 3 × 10-8 M/yr as it shows just a transition state between = 5 × 10-8 M/yr and = 1 × 10-8 M/yr, since they have no new features, as can also be guessed from the structures shown Fig. 3. However, in this case the inner and outer region of outward migration shrink and are shifted inwards. In fact when arriving at = 1 × 10-8 M/yr only one region of outward migration remains in our computational domain. This region is then much smaller in size and valid for a narrower range of planetary masses. But, it also implies that outward migration is possible for planets with a slightly smaller mass (5 MEarth).

The migration map for = 5 × 10-9 M/yr (bottom in Fig. 4) looks completely different as for the higher models. In fact there are no more regions of outward migration left. As stated before, outward migration is stronger if H/r drops, but for this case, no drop in H/r can be observed (see Fig. 3), leading to inward migration for all planetary masses and at all orbital distances. However, in the region where we observed outward migration for higher discs, the inward migration is slowest. In the region between 1.5 and 2 AU the inward migration is actually increased compared to the rest of the disc. This is supported by the temperature gradient of the disc, which actually seems to be even partly positive (see Fig. 3), which can even result in a negative corotation torque. But, keep in mind here that the actual migration rate is smaller for reduced disc masses, as Γ0 scales linearly with ΣG, see Eq. (10).

Interestingly, only the torque formula that accounts for saturation gives a negative torque for the whole disc, while the fully unsaturated torque formula (Paardekooper et al. 2010) still predicts outward migration. Infact in a disc, one can relate the total unsaturated torque, which consists of the Lindblad torque (Eq. (10)) and the barotropic (Eq. (12)) and entropy related (Eq. (11)) part of the corotation torque, to the flaring index of the disc. Assuming that ΣGν is independant of r as well as α (which gives s = 3 / 2 − β), and that γ = 1.4, one finds by using Eq. (8) (13)where b is the flaring index of the disc. This indicates that even the unsaturated torque is negative in discs with flaring indices larger than 0.1. By interpolating the trend of the H/r profile towards even smaller discs, we speculate that even the unsaturated torque formula will predict inward migration as the becomes smaller and smaller. However, the torques saturate in time, so that the formula for the unsaturated torques should not be used, especially for small discs (see discussion in Sect. 6).

4. -discs with different metallicity

In this section we explore the influence of different metallicity on the disc structure and the corresponding migration maps. A reduced metallicity in small grains could happen in time as the first planetesimals and embryos form in the disc and reduce the amount of heavy elements in the disc. Additionally, radial drift could wash out the dust grains, reducing the metallicity. A higher metallicity in small μm size dust grains in the late stages could be achieved by collisional grinding of planetesimals that produce small particles.

The change of metallicity has many implications, also for planet formation. For example, a higher metallicity helps the streaming instability to operate (Johansen & Youdin 2007) and the formation of planetesimals (Johansen et al. 2007), while a lower metallicity hinders the streaming instability to work.

4.1. Higher metallicity

A change in the metallicity of the disc translates to a change of the opacity of the disc. If the metallicity of a disc is doubled, so is the opacity. If the metallicity increases, the cooling rate of the disc decreases as can be seen from Eq. (15). In Fig. 5, it is shown that the aspect ratio H/r of discs increases with increasing metallicities, compared to the low metallicity cases (Fig. 3).

In the outer parts of the disc, the structure is quite similar to the one shown for the low metallicity discs (Fig. 3). This is not obvious, as the cooling rate reduces with increasing metallicity while the absorption of stellar irradiation increases. However, due to more grains in the disc, the absorption of stellar irradiation (which is ρκ) will happen for smaller densities, which implies that it happens at an higher altitude in the disc. Therefore, the heat that transfers down to mid-plane is actually reduced. But this effect is compensated by the reduced cooling, so the resulting disc structure in the outer parts is the same.

Clearly the aspect ratio is higher in the inner parts of the disc for the = 5 × 10-8 M/yr disc with a metallicity of 0.02 compared to the disc with a metallicity of 0.01 shown in Fig. 3. We recall that the inner part of the discs are dominated by viscous heating. The viscous heating does not depend on the metallicity, but only on the gas density. The higher metallicity reduces the cooling rate and therefore the balance between viscous heating and cooling turns to higher temperatures, hence to a higher h. This effect, of course, diminishes as the rate is reduced and therefore viscous heating as well.

thumbnail Fig. 5

H/r-profile of discs with different rates and metallicities. The profile is taken when the discs have reached their radially constant rate. The colour coding of the high metallicity discs match the rates displayed in Fig. 3. Please note that we cut the displayed disc at 30 AU to enhance the visibility in the inner parts of the disc.

Open with DEXTER

Interestingly, not only is the aspect ratio larger for higher metallicity discs, but also the bumps in the inner parts of the disc are more pronounced. For the = 5 × 10-9 M/yr disc, the bumps are actually visible for the first time. The increase of the bumps inside the disc structure is related to the reduced cooling rate in the high metallicity discs. This can have important implications on the migration of giant planet cores.

thumbnail Fig. 6

Torque acting on discs with different and metallicity, with 5 × 10-8 M/yr and 0.02 metallicity (top), 1 × 10-8 M/yr and 0.03 metallicity (middle) and 5 × 10-9 M/yr and 0.05 metallicity (bottom). The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K. The migration maps correspond to the disc profiles shown in Fig. 5.

Open with DEXTER

In Fig. 6, we present the migration maps for the same discs as in Fig. 4, but with different metallicities. The metallicities are 0.02 (top), 0.03 (middle) and 0.05 (bottom). The clear difference for the migration maps displayed here, with the ones in Fig. 4 is that the regions of outward migration are much larger and more extended. Additionally, the migration map for the 5 × 10-9 M/yr disc shows outward migration, which was not visible for the for disc with a metallicity of 0.01. This is clearly related to the change in the disc structure in the inner parts (see Fig. 5), where now a bump appears.

4.2. Lower metallicity

Due to the formation of planetesimals and planetary embryos, which require dust and ice grains, the amount of small solids can decrease, if there is no additional source of small grains. Discs with reduced metallicity, have a reduced opacity inside the disc, which works in the same way as described Sect. 4.1. The aspect ratio of discs with a metallicity of 0.001 are also displayed in Fig. 5.

Clearly, a lower metallicity leads to a higher cooling rate (Eq. (15)), which explains the drop of H/r in the inner parts of the disc. The outer parts of the disc remain flared and at about a constant level, as explained in Sect. 4.1. For the = 1 × 10-8 M/yr disc there is a dramatic change of the disc structure as no bumps in the discs profile are visible any more, which are still visible in the case of a metallicity of 0.01.

In Fig. 7 the migration maps for the low metallicity discs are displayed. In the = 5 × 10-8 M/yr case (Fig. 7, top), the regions of outward migration are much smaller compared to the 0.01 metallicity case (top of Fig. 4). In fact, the inner region of outward migration is dramatically reduced and the outer region is only valid to much smaller masses. Additionally, the outer region of outward migration is shifted inwards from 5 AU to 3.8 AU, which can be explained due to the higher cooling rate in the disc, which shifts the ice line further inwards.

thumbnail Fig. 7

Torque acting on discs with different and metallicity, with 5 × 10-8 M/yr and 0.001 metallicity (top) and 1 × 10-8 M/yr and 0.001 metallicity (bottom). The black lines encircle the regions of outward migration. The vertical red (top) and blue (bottom) lines indicate the ice line at 170 K. The migration maps correspond to the disc profiles shown in Fig. 5. Please note the different colour scale for the bottom plot.

Open with DEXTER

In the case of = 1 × 10-8 M/yr (Fig. 7, bottom), the changes are dramatic, as no region of outward migration exists any more. In fact, the region beyond the ice line now features a region of enhanced inward migration. In all other displayed migration maps, we observe a region of outward migration (or at least with slower inward migration). At this point in the disc, the aspect ratio profile undergoes a steep gradient (Fig. 5), which actually refers to a positive gradient in temperature, which leads to a negative entropy related corotation torque, resulting in an even larger negative total torque. The reason for that is probably related to the transition of viscous heating to stellar heating. As viscous heating is operating in the inner (mid-plane) parts of the disc, heat is produced there. Stellar irradiation heats the upper layers of the disc, which then transport heat down towards mid-plane. At 2.5 AU stellar heating is getting absorbed effectively, which heats the mid-plane, but viscous heating seems to diminish already at a smaller distance to the star. Therefore a higher temperature can be observed at 2.5 AU compared to the inner parts of the disc, which leads to a positive temperature gradient in the disc.

5. Influence of the interchange between ΣG and ν

As the mass flux through a disc is proportional to the gas surface density ΣG and the viscosity ν, a change in ΣG could be compensated by a change in ν to still have the same . We now want to construct a scenario where and the thermal structure are the same, but where ΣG and ν are not. For changing ν, we change α.

Let us look at the different heating and cooling sources of the disc. The viscous heating, which is the dominant heat source in the inner parts of the disc, is proportional to ΣGν (): (14)which implies that an interchange between ΣG and ν would not change the viscous heating. The radiative cooling is given by (15)where c is the speed of light, λ the flux-limiter (Kley 1989) and ER the radiation energy. In order to get the same cooling rate, a change in ρG has to be compensated by a change in κR. The same applies for the absorption of stellar photons, which scales with ρGκ. So, if the gas density is increased by a factor sf, the opacity, and therefore the metallicity, has to decrease by a factor sf. The ratio of ΣG/ ΣZ will change, but this will not change the thermal structure of the disc. This actually means that the surface density of heavy elements ΣZ stays constant.

Therefore we can use the data of the discs presented in the previous section and scale them by scaling factor sf to interchange ΣG and ν. When increasing ν by sf, the gas surface density ΣG has to decrease by sf. At the same time ΣZ has to be constant, but this means that the metallicity of the disc increases. So in scenario like that, the opacities κR, κP and κ increase by sf as well.

thumbnail Fig. 8

Torque acting on discs with = 1 × 10-8 M/yr for a rescaling of ΣGν by a factor of sf = 0.5 (top) and a factor of sf = 2 (bottom). This means that the disc displayed in top has a metallicity of 0.005, while the disc in bottom has a metallicity of 0.02. The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K.

Open with DEXTER

In Fig. 8 we present the migration maps of a = 1 × 10-8 M/yr disc, for two different cases for exchanging ΣG and ν. In the top, the viscosity is decreased by a factor of 2 and in the bottom plot it is increased by a factor of 2. This means that the top disc has 0.005 metallicity and the bottom one 0.02.

These two migration maps significantly differ from each other and from the one presented in Fig. 4. In the case of low viscosity (top in Fig. 8), no outward migration is detected any more. But in the case of high viscosity, a much larger region of outward migration is visible. However, outward migration in this case starts a much higher planetary masses (12 MEarth) compared to the normal viscosity case presented in Fig. 4, where outward migration starts for 5 MEarth.

The reason for this is the scaling of the thermal diffusion coefficient χ in the formula of Paardekooper et al. (2011), which is not linear in ρGκR. It reads: (16)where σ is the Boltzmann constant. Please note that a factor of 4 is missing in Paardekooper et al. (2011) according to Bitsch & Kley (2011). Additionally, the saturation parameters for the corotation torque and horseshoe drag are also dependent on the viscosity of the disc. These saturation parameters are responsible for the transition from horseshoe drag towards the linear corotation torque, which is much smaller than the horseshoe drag. This can then lead to inward migration, changing the migration maps.

This non-linearity causes different migration behaviour for discs with different ν to ΣG ratios, even if they feature the same value and the same thermal structure. In total, one could roughly say that a smaller viscosity allows outward migration for smaller planetary masses (unless it gets too low and no outward migration is possible any more), while a larger viscosity allows outward migration for larger planetary masses, where the minimum mass required for outward migration is also increasing. This is of crucial importance for N-Body simulations of planetary embryos in evolving discs, where just stating a value would not be enough.

6. Consequences on planet formation

The presented migration maps have important consequences for the resulting structures of planetary systems. Lyra et al. (2010) stated that planets formed in a disc, would migrate and then sit at the zero-migration radius, which is independent on the planetary mass in their case, because they use the torque formula of Paardekooper et al. (2010) that only accounts for the unsaturated torques. This picture is contradicted by the results presented here, as we use the torque formula with saturation (Paardekooper et al. 2011).

The planets in Lyra et al. (2010) consequently move inwards with the zero-torque radius in time as the disc disappears, which can also be observed from the migration maps presented here. However, they used a simple 1D model for their disc evolution and did not account for stellar heating, which clearly changes the disc structure and the migration pattern compared to fully radiative discs (Paper I). In particular Eq. (13) shows that the flaring index should be determined accurately for a good estimate of the unsaturated torque, but non stellar irradiated discs can not be flared in the outer parts (Paper I). The main difference between our model and the Lyra et al. (2010) model originate from the more complex disc structure presented here and from the different torque formula that accounts for saturation effects.

Lyra et al. (2010) also stated that in the late evolutionary stages of the disc, starting at 4.0 Myr, the surface density is so low that it can not transfer enough angular momentum to the planet for its orbital radius to evolve as fast as the equilibrium radius. This means that the planets detach from the evolution of the disc and are left behind. The beginning of their late stage of 4.0 Myr matches from the surface density at 1 AU quite well with our = 5 × 10-9 M/yr model with a metallicity of 0.01, where we do not observe a zero-torque radius any more for the torques that are prone to saturation (Paardekooper et al. 2011). This would indicate that planetary cores could start inward migration while the disc is disappearing and can not safely wait and detach themselves from the disc as proposed by Lyra et al. (2010). Therefore this could make the survival of planets in the outer parts of the disc in the late stages much harder. This problem becomes even more severe if the metallicity of the disc would drop in time (e.g. due to the formation of planetesimals out of dust grains) as can be seen from Fig. 7, where only inward migration is observed for = 1 × 10-8 M/yr.

The results of Lyra et al. (2010) could be interpreted in a way that the evolution of the disc naturally would account for the occurrence of Super-Earth and Mini-Neptune planets in the inner parts of the solar system. However, if the metallicity of the disc is large in the late stages of the evolution, the zero-migration radius for these planets is at 4 AU (see Fig. 6), which is much larger than in the model of Lyra et al. (2010), even though we used the torque formula accounting for saturation here. Additionally in the low metallicity case, the zero-torque radius for = 5 × 10-9 M/yr using the fully unsaturated torque is much further out in our simulations (5 AU) compared to 2 AU in Lyra et al. (2010), which is clearly a consequence of the different disc models.

From our simulations, we propose two different ways to form hot Neptunes. Either they form in the outer disc, where they get trapped in the outer region of outward migration and are then released at a late stage where they migrate inwards to the inner edge of the disc (low metallicity scenario). Or they assemble in the inner part of the disc from migrating embryos that were too small to be trapped in the regions of outward migration. Note that the second scenario can happen for low and high metallicities, as both scenarios require a minimal planetary mass to halt inward migration.

The migration maps for a metallicity of 0.01 (Fig. 4) can also be used as a guideline to discuss the evolution of giant planets during their formation. The core of a giant planet most likely forms beyond the ice line at the early stages of the disc. The reason for the core to actually grow to become a giant planet and not stay a hot Neptune, might be explained in the different formation time of the core itself. A faster formation of the core provides more time to accrete gas to become a gas giant, than a later formation of the core where the planet is stuck at Neptune size. Let us assume that a core is trapped at the outer zero-migration radius where it can grow until it reaches a mass of 40 MEarth. Thus, at that stage the planet is in a region between 4 and 5 AU, depending on which evolutionary stage the disc is in (the smaller is , the closer to the star is the equilibrium radius). When the 40 MEarth mass is exceeded, the planet is released and starts migrating towards the central star. A planet of that size is growing rapidly due to gas accretion (Pollack et al. 1996) and it starts to open a significant gap inside the disc changing its migration regime to that of type-II-migration. Due to this migration phase, the giant planet, once fully formed, is significantly closer to the star than the original location of its core, i.e. closer than ~4 AU. The formation of an inner cavity in the disc by photo-evaporation can eventually stop the inward migration of the planet, as suggested by Alexander & Pascucci (2012), preventing it from becoming a hot Jupiter.

However, in the case of our own solar system, this is not sufficient. A mechanism is needed to move the giant planet (Jupiter) outwards and bring it to 5 AU. This consideration supports the Grand Tack scenario (Walsh et al. 2011), in which the formation of Saturn eventually forces Jupiter to reverse its migration direction and move outwards. In the Grand Tack scenario, Saturn, once formed, migrates faster than Jupiter and cached it in resonance inside their common gap, which causes the migration reversal (Masset & Snellgrove 2001; Morbidelli & Crida 2007; Pierens & Nelson 2008). A potential problem, often invoked with this scenario, is that a priori Saturn and Jupiter should have followed exactly the same growth-migration histories, so that Saturn could never catch up with Jupiter. Our work now shows that this assumption is actually not true. If Saturn forms after Jupiter, will have decreased in the disc, and the migration map will differ. In fact, as times goes by, the maximum mass for being still trapped in a zero-torque zone decreases. So, Saturn should start migration inwards at a smaller mass than Jupiter did. Possibly, Jupiter never had inward type-I-migration, while Saturn had, and may even have experienced fast type-III-migration, which leads Saturn to catch up with Jupiter.

Of course this is dependent on the ratio between the migration speed of the giant planets and the actual accretional evolution of the disc. In addition, opening of a gap by Jupiter would change the thermal structure of the disc, hence changing the migration map for Saturn. A detailed study of this process is beyond the scope of this paper, but it is worth noting that our results qualitatively support the concept of the Grand Tack scenario.

It also seems that it is easier to keep the first cores in discs with higher metallicity, as the minimum mass needed for outward migration is reduced compared to lower metallicity discs. This would potentially make it easier to keep the first cores that form giant planets. As potentially more cores can be kept inside the disc, it is more likely that these cores can form giant planets. Therefore, high metallicity discs would imply a higher probability to form gas giant planets, which is supported by observations.

Hasegawa & Pudritz (2011) stated that among the ice line and dead zone, the heat transition between viscous and stellar heating can act as a planet trap. In a disc, stellar heating is dominant in the outer parts of the disc and therefore responsible for the flaring part of the disc (see Chiang & Goldreich 1997, Paper I). In the inner parts viscous heating is dominant. The different heat sources result in a change of the gradient of temperature between these two different regions, which can lead to a region of outward migration, according to Hasegawa & Pudritz (2011). Additionally they stated that the ice-line (located at rice) can function as a planet trap.

However, in the simulations presented here, we see a different picture. Here, the ice line functions as a region of divergent migration, meaning that planets interior to the ice line (r<rice) migrate inwards, while planets located outside (r > rice) the ice line migrate outwards. Additionally, one could argue that the disc is locally radial isothermal at the ice line, as melting and resublimating of ice grains keep the temperature constant for a certain distance around the actual ice line (Lyra et al. 2010; Hasegawa & Pudritz 2011), which is not explicitly taken into account in our simulations, but it is slightly reflected by the smoothed transition between the different branches of opacity (Fig. 1). As the disc is locally radial isothermal, no temperature gradient exists around the ice line, so that outward migration is not possible there and therefore the ice line can not act as a point of convergent migration.

The outward migration beyond the ice line then stops in a region of the disc, where H/r reaches a local minimum, which is a convergent migration radius in the disc. A minimum of H/r can only exist, if H/r increases again, which, in this case, is achieved by stellar irradiation. In that sense, a change form inward to outward migration is caused by stellar irradiation, because it flares the disc. But at the same time the change of the opacity at the ice line, provides a change in the disc structure as the cooling properties of the disc change. In this sense, the divergent migration point is related to the opacity transition at the ice line and the stopping point of outward migration is related to stellar heating, as the disc becomes flared (Paper I).

The turbulence inside the disc is thought to be driven by the MRI, which would then also account for a non-turbulent inner part of the disc, the so called dead zone. The dead zone basically is a region of low viscosity inside the disc, which can have important consequences on the disc structure and on the migration of embedded planets. The study of the influence of the dead zone on disc structures is planned in future work. Additionally to that, we plan to use the presented code in 3D to test the influence of embedded planets on the disc structure of stellar irradiated discs and if the migration behaviour is still similar to the one predicted by the torque formulae. With the presented code these simulations seem to be possible in a reasonable amount of computation time (several weeks).

7. Summary

We presented here hydrodynamical simulations of stellar irradiated accretion discs, featuring viscous heating, radiative cooling and stellar irradiation. These discs feature different radially constant accretion fluxes = 3πνΣG through them that represent different stages of the lifetime of an accretion discs, as the rate is reducing as the disc becomes older. We especially focused on the thermal structure of the disc. From the resulting disc structure, we computed the torque acting on embedded planets by using a torque formula (Paardekooper et al. 2011) in order to find sweet spots for planetary growth.

At the ice line, the opacity profile (Fig. 1) changes due to the melting and condensation of ice grains. This transition in opacity changes the stellar heating and radiative cooling rates, resulting in bumps and dips in the inner disc. These bumps and dips in the aspect ratio H/r of the disc, translate into changes in viscosity, ν = αH2ΩK. This therefore causes changes in the gas surface density ΣG as the product of viscosity and surface density is constant (the -rate). The resulting disc structure (density and temperature) is therefore not uniform and does not follow a simple power law. Additionally, the bumps and dips are smaller and closer to the central star when the rate reduces. In the final stages, the disc is flared with no bumps in the inner part.

The thermal structure of the disc depends on the gas surface density ΣG and the viscosity ν in such away that they are interchangeable, as long as the surface density of heavy elements ΣZ is constant. In other words, for discs with a given ΣG × ν product, the same disc structure is reached as long as ΣZ is the same. Additionally, we pointed out that a non linear scaling in the torque calculations exists, which indicate that an interchange between ΣG and ν leads to a different migration behaviour, even when the thermal disc structure is the same. In fact, a lower disc viscosity results in a possible outward migration for lower mass planets than a higher viscosity disc. We can therefore conclude that the thermal structure is a two parameter space in ΣGν and ΣZ, while the migration is a three parameter space in ΣGν, ν and ΣZ.

As small mass cores can migrate, collide and form bigger objects at convergent migration distances (Pierens et al. 2013), the location of these is of high importance. We calculate the migration by applying the (Paardekooper et al. 2011) torque formula onto our disc structures. We should point out here, that in all the disc scenarios, a minimum planetary mass is needed in order to find outward migration. This mass lies (dependent on the underlying disc structure) in the range of 3 to 5 MEarth. This is caused by a transition from the horseshoe drag towards the linear corotation torque in viscous discs, which reduces the total torque acting on an embedded planet.

We find a nice match between the regions in the disc where the opacity law changes (and therefore the disc structure) and where the migration changes. In fact, we confirm here the results of Paper I, where we stated that outward migration is more likely in regions where H/r drops. Of course, this means that for smaller bumps and dips in the H/r profile, a smaller region of outward migration is the consequence. For a metallicity of 0.01 outward migration is no longer supported in the late stages of the disc (low ).

As the disc accretes onto the star, the metallicity of the disc can change. Our simulations show that a disc with metallicity of 0.001 actually gives rise to only inward migration, compared to discs with a metallicity of 0.01. A larger metallicity, on the other hand, allows outward migration for discs with an even smaller rate. So, if the metallicity of the disc increases in the late stages, it is easier to keep planets in the disc with gas driven migration.

The migration map is therefore not fixed in time, but evolves as and the metallicity change. This if of crucial importance for the migration scenarios. This influences the formation of hot Neptunes that can halt their migration in these planet traps, because it allows for two formation scenarios. Either they form in the outer disc, get trapped at the convergent migration radius and are then released towards the inner disc in the late stages, or they assemble from smaller bodies that migrate towards the inner disc, that then grow and stay there. Additionally to that, the evolution of the migration maps give also indications about the final location of giant planets. Especially in the solar system case our migration maps indicate that Jupiter would be formed well inside 4 AU and thus needs a mechanism that transports it outwards again, which is provided by the Grand Tack scenario (Walsh et al. 2011).

Acknowledgments

B. Bitsch and A. Morbidelli have been sponsored through the Helmholtz Alliance Planetary Evolution and Life. The Nice group is thankful to ANR for supporting the MOJO project. The computations have been done on the Mesocentre SIGAMM machine, hosted by the Observatoire de la Côte d’Azur. We thank W. Kley and C. Dullemond for useful discussions. Additionally we want to thank an anonymous referee for comments that helped to improve the manuscript.

References

Appendix A: Boundary conditions

Appendix A.1: Inner boundaries

We impose open inner boundaries, which means the velocity can be direct inwards towards the star and material can leave the disc. For the radial velocity at the boundary between active and ghost domain follows: (A.1)where vr,A is the radial velocity in the first active grid cell. A negative velocity implies an inward flux (towards the central star). The polar velocity vθ is copied from the active domain and the azimuthal velocity is set via the following equation: (A.2)where rA and rG are the distances of the active and ghost cells to the central star. Note that rA>rG.

The density in the ghost cells is copied from the active cells to the ghost cells to mimic the outflow condition. This applies to all other cell centred quantities (ϵ, T, ER) as well.

Appendix A.2: Outer boundaries

From the starting configuration with a given H, the disc evolves into its final state with a different H. A change of H also implies a change of the viscosity, ν = αcsH. But the viscosity is vertically constant, as is predicted by MHD simulations (Flaig et al. 2012). We impose the same H/r in the ghost and active cell, where we can now calculate the height in the ghost cell HG and the sound speed cs,G. The sound speed in mid-plane is used to calculate the viscosity in the ghost cells, where then an averaged viscosity between active and ghost cells gives the velocity (as the velocity is defined on the edges of the grid cells): (A.3)where cs is the sound speed in the ghost cells at mid-plane height. The minus sign defines an inward (to the star) velocity.

The polar velocity vθ is copied from the active domain and the azimuthal velocity is set via the following equation: (A.4)where rA and rG are the distances of the active and ghost cells to the central star. Note that rA<rG.

The density in the ghost cells has to be provided by . Taking into account that ρ varies with height, we get (A.5)The vertical density distribution (hydrostatic equilibrium) is given by (A.6)with H = csΩ. ρ0 can be determined by putting Eq. (A.6) in Eq. (A.5) and solving as a function of . The resulting ρ0 will then be used to calculate ρ(z) in the ghost cells. This way, the density can change in time as the viscosity changes until the equilibrium state is reached, but the disc will always have the same inflow. The other cell centred quantities are copied from the active ring, as it was done for the inner boundaries.

All Tables

Table 1

Simulation parameters for the used models.

All Figures

thumbnail Fig. 1

Opacity profiles of the opacities used in the code. The different bumps in the profile correspond to opacity transition, where the ice line is located at 170 K. The grey area marks the region where the opacity changes due to the melting of ice grains. The maximum temperature during the simulation is reached in the inner parts of the disc, where viscous heating is dominant and can reach up to 900 K for the 1 × 10-7 M/yr case. Opacities are plotted for a metallicity of 0.01.

Open with DEXTER
In the text
thumbnail Fig. 2

Evolution of H/r in the disc for the different stages. The disc features = 2 × 10-7 M/yr. At the beginning, in the (H/r)0 stage, α0 = 0.001, which we change to α1 = 0.003 after the stage (H/r)1 is reached. The stage (H/r)2 shows the aspect ratio after the viscosity has change and (H/r)final shows the final stage after the surface density has been re-arranged.

Open with DEXTER
In the text
thumbnail Fig. 3

H/r-profile (top), temperature (middle) and surface density profile (bottom) of discs with different rates. The profile is taken when the discs have reached their radially constant rate. The black dotted lines represent power-law fits to guide the eyes. The grey area in the temperature plot marks the temperature range where the opacity profile changes due to the melting of ice grains (see Fig. 1). Please note that we cut the displayed disc at 30 AU to enhance the details in the inner parts of the disc.

Open with DEXTER
In the text
thumbnail Fig. 4

Torque acting on discs with different , with 5 × 10-8 M/yr (top), 1 × 10-8 M/yr (middle) and 5 × 10-9 M/yr (bottom) for the applied Paardekooper et al. (2011) formula. The black lines encircle the regions of outward migration for the Paardekooper et al. (2011) formula. The vertical solid yellow lines mark the outer edge of the zero-migration region for the unsaturated torques (Paardekooper et al. 2010), while the dashed-yellow line marks the inner edge of the zero migration region of the unsaturated torques. The vertical red lines indicate the ice line at 170 K.

Open with DEXTER
In the text
thumbnail Fig. 5

H/r-profile of discs with different rates and metallicities. The profile is taken when the discs have reached their radially constant rate. The colour coding of the high metallicity discs match the rates displayed in Fig. 3. Please note that we cut the displayed disc at 30 AU to enhance the visibility in the inner parts of the disc.

Open with DEXTER
In the text
thumbnail Fig. 6

Torque acting on discs with different and metallicity, with 5 × 10-8 M/yr and 0.02 metallicity (top), 1 × 10-8 M/yr and 0.03 metallicity (middle) and 5 × 10-9 M/yr and 0.05 metallicity (bottom). The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K. The migration maps correspond to the disc profiles shown in Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 7

Torque acting on discs with different and metallicity, with 5 × 10-8 M/yr and 0.001 metallicity (top) and 1 × 10-8 M/yr and 0.001 metallicity (bottom). The black lines encircle the regions of outward migration. The vertical red (top) and blue (bottom) lines indicate the ice line at 170 K. The migration maps correspond to the disc profiles shown in Fig. 5. Please note the different colour scale for the bottom plot.

Open with DEXTER
In the text
thumbnail Fig. 8

Torque acting on discs with = 1 × 10-8 M/yr for a rescaling of ΣGν by a factor of sf = 0.5 (top) and a factor of sf = 2 (bottom). This means that the disc displayed in top has a metallicity of 0.005, while the disc in bottom has a metallicity of 0.02. The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.