Free Access
Issue
A&A
Volume 575, March 2015
Article Number A28
Number of page(s) 17
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201424964
Published online 16 February 2015

© ESO, 2015

1. Introduction

Protoplanetary discs evolve over million year time scales during which the accretion rate onto the central star drops, from typically 5 × 10-6 M/yr around the youngest stars to 10-9 M/yr after 3–10 Myr (Hartmann et al. 1998). This reduction in can be interpreted as a reduction of the gas surface density ΣG. Ensembles of protoplanetary discs can be observed in star forming regions, like the Ophiuchus cluster (Andrews et al. 2010). These observations give constraints on the structure of protoplanetary discs (e.g. gradients in surface densities), although with large uncertainties. In the late stages of the discs, becomes so small that the disc is blown away by photo evaporation on a very short time scale (Alexander et al. 2014). The central star also evolves on these time scales (Baraffe et al. 1998). As the stellar luminosity changes, the amount of stellar heating received by the disc changes as well, which in turn affects the temperature and density structure.

During this evolution, planets form. Gas giant planets have to form, while the gas-rich protoplanetary disc is still present, and large parts of the growth of terrestrial planets and super Earths occur during the gaseous disc phase as well. Several important growth and formation stages rely on detailed knowledge of the protoplanetary disc structure:

Even though these mechanisms happen on different length scales and time scales, they are all strongly dependent on the underlying disc structure (temperature T, gas surface density ΣG, aspect ratio H/r), which makes the protoplanetary disc structure therefore a key parameter to understand the formation of planets and planetary cores. We now highlight some key processes, which will be discussed further in this paper.

Planetesimal formation and the formation of planetary embryos can be aided by reducing the radial pressure gradient in the protoplanetary disc. Bumps and dips in the surface density of the disc, which is where the pressure gradient changes, can locally reduce the pressure support of the gas. A reduction in the pressure support of the gas causes a reduction of the headwind acting on the pebble and would reduce their inward motion (Weidenschilling 1977a; Brauer et al. 2008; Birnstiel et al. 2012). Reduced pressure support also decreases the metallicity threshold for particle concentration through the streaming instability, significantly helping the formation of planetesimals (Bai & Stone 2010b).

Core accretion is the process where planetary embryos grow to the size of the cores of the giant planets (~10 MEarth). The accretion of planetesimals by embryos is slower than disc life times (Pollack et al. 1996; Rafikov 2005; Levison et al. 2010), but growth time scales can be drastically reduced by considering the accretion of pebbles (Johansen & Lacerda 2010; Ormel & Klahr 2010; Lambrechts & Johansen 2012; Morbidelli & Nesvorny 2012). In the latter scenario the cores can grow on time scales of 106 years at a radial distance of 5 AU from the host star (Lambrechts & Johansen 2014). Such fast growth depends on the point where a core grows large enough to enter a phase of rapid accretion (so-called Hill accretion), which occurs earlier in regions with lower scale heights.

Planetary migration describes the gravitational interactions of planets and planetary cores with the gas disc (Ward 1997). In a locally isothermal disc, cores are expected to migrate inwards on time scales that are much shorter than the discs’ lifetime (Tanaka et al. 2002), posing a problem for the formation of the cores of giant planets. However, the migration of cores cores depends on the thermodynamics in the disc, and even outward migration of planetary cores is possible (Paardekooper & Mellema 2006; Kley & Crida 2008; Baruteau & Masset 2008; Kley et al. 2009). The migration depends on the local radial gradient of entropy, which can drive outward migration (see Baruteau et al. 2014 for a review).

The structure of the protoplanetary disc that surrounded our own Sun can be approximated by the minimum mass solar nebula (MMSN) which comes from fitting the solid mass (dust and ice) of the existing planets in our solar system with a power law (Weidenschilling 1977b; Hayashi 1981). It is then often assumed that other protoplanetary discs have similar power law structures (Chiang & Laughlin 2013), however applying it to all extrasolar systems is troublesome (Raymond & Cossou 2014). Along with the MMSN model, Chiang & Youdin (2010) propose a model that features slightly different gradients in the disc (hereafter named CY2010). Both models, the MMSN and the CY2010 model, approximate the disc structure by uniform power laws in temperature T and gas surface density ΣG. The difference in the two models originates in updated condensate mass fractions of solar abundances (Lodders 2003) that lead to a higher estimate of the surface density. The different temperature profile is explained by different assumptions on the grazing angle of the disc that determines the absorption of stellar irradiation. In the CY2010 model that follows the calculations of Chiang & Goldreich (1997), T(r) ∝ r− 3 / 7, in contrast to the MMSN where T(r) ∝ r− 1 / 2 stemming from the assumption of an optically thin disc.

However, simulations including radiative cooling and viscous and stellar heating have shown that the disc structure features bumps and dips in the temperature profile (Bitsch et al. 2013, hereafter Paper I). Bitsch et al. (2014a, hereafter Paper II) highlighted the direct link in accretion discs between changes in opacity κ and the disc profile. The mass flux is defined for a steady state disc with constant at each radius r as (1)Here ΣG is the gas surface density and vr the radial velocity. Following the α-viscosity approach of Shakura & Sunyaev (1973) we can write this as (2)Here H is the height of the disc and ΩK the Keplerian rotation frequency. At the ice line, the opacity changes because of the melting and the sublimation of ice grains or water vapour. This change in opacity changes the cooling rate of the disc [D ∝ 1 / (ρκ)], hence changing the temperature in this region. A change in temperature is directly linked to a change in H, so that the viscosity changes. However, since the disc has a constant rate at all radii, a change in viscosity has to be compensated by a change in surface density, creating bumps and dips in the disc profile, which do not exist in the MMSN and CY2010 model.

We present detailed 2D (r,z) simulations of discs with constant , including stellar and viscous heating, as well as radiative cooling. We use the α approach to parametrize viscous heating, but the protoplanetary disc is not temporally evolved with the α approach. Instead α is used to break the degeneracy between column density ΣG and viscous accretion speed vr (Eq. (1)). Each value of is linked to an evolutionary time through observations (Hartmann et al. 1998). This then allows us to take the correct stellar luminosity from stellar evolution (Baraffe et al. 1998). In contrast to previous work (Bitsch et al. 2014b), we do not include transitions in the α parameter to mimic a dead zone, because we are primarily interested in investigating how the stellar luminosity and the disc’s metallicity affect the evolution of the disc structure in time. For the time evolution of the disc over several Myr (and therefore several orders of magnitudes in ), we present a semi-analytical formula that expresses all disc quantities as a function of and metallicity Z.

This paper is structured as follows. In Sect. 2 we present the numerical methods used and the link between and the time evolution. We point out important differences in the disc structure between our model and the MMSN and CY2010 model in Sect. 3. We show the influence of the temporal evolution of the star on the structure of discs with different in Sect. 4. We then discuss the influence of metallicity on the disc structure and evolution in Sect. 5. We then discuss the implications of the evolution of protoplanetary discs on planet formation in Sect. 6. We finally summarize in Sect. 7. In Appendix A we present the fits for the full time evolution model of protoplanetary discs.

2. Methods

2.1. Numerical set-up

We treat the protoplanetary disc as a three-dimensional (3D), non-self-gravitating gas, whose motion is described by the Navier-Stokes equations. We assume an axisymmetric disc structure because we do not include perturbers (e.g. planets) in our simulations. Therefore we use only one grid cell in the azimuthal direction, making the computational problem de facto 2D in the radial and vertical direction. We utilize spherical coordinates (r-θ) with 386 × 66 grid cells. The viscosity is treated in the α-approach (Shakura & Sunyaev 1973), where our nominal value is α = 0.0054. Here the viscosity is used as a heating parameter and not to evolve the disc viscously, because the viscous evolution of the disc is very slow. Instead we use an initial radial gas density profile from a 1D analytic model for each accretion rate . The vertical profile is computed from an analytic model of passive discs irradiated by the central star. The simulations are then started until they reach an equilibrium state between heating and cooling, which happens much faster than the viscous evolution of the disc. This final equilibrium state is different from the initial one, because the disc is not passive (viscous heating is included) and opacities depend on the temperature. This changes H/r in the disc, which in turn changes the local viscosity because of the α-prescription. Therefore we continue the simulations until a new radial density profile in the disc is achieved. This happens on a viscous time scale, but because the changes are only local variations relative to the initial profile, the new equilibrium state is achieved relatively fast, much faster than the global evolution of the disc and decay of the stellar accretion rate (see Sect. 2.2).

The viscosity in protoplanetary discs can be driven by the magnetorotational instability (MRI), where ionized atoms and molecules cause turbulence through interactions with the magnetic field (Balbus & Hawley 1998). As ionization is more efficient in the upper layers of the disc (thanks to cosmic and X-rays), the midplane regions of the disc are not MRI active and therefore feature a much smaller viscosity (α parameter). In discs with a constant , a change in viscosity has to be compensated for by an equal change in surface density (Eq. (2)); however, in 3D simulations, much of the accretion flow can be transported through the active layer, so that the change in surface density is smaller than for 2D discs (Bitsch et al. 2014b). Additionally, hydrodynamical instabilities, such as the baroclinic instability (Klahr & Bodenheimer 2003) or the vertical shear instability (Nelson et al. 2013; Stoll & Kley 2014), can act as a source of turbulence in the weakly ionised regions of the disc. A realistic picture of the source of turbulence inside accretion discs is still being debated (see e.g. Turner et al. 2014). We therefore feel that it is legitimate to neglect the effects of a dead zone and assume a constant α throughout the disc.

The dissipative effects can then be described via the standard viscous stress-tensor approach (e.g. Mihalas & Weibel Mihalas 1984). We also include irradiation from the central star, as described in Papers I and II. For that purpose we use the multi-dimensional hydrodynamical code FARGOCA, as originally presented in Lega et al. (2014) and in Paper II. The radial extent of our simulations spans from 1 AU to 50 AU, which includes the range of the MMSN, which is defined from 0.4 to 36 AU. We apply the radial boundary conditions described in the Appendix of Paper II.

The radiative energy associated with viscous heating and stellar irradiation is transported 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 that 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. The flux F from the central star is given by (3)where R and T give the stellar radius and temperature and σ is the Stefan-Boltzmann constant. Stellar heating is responsible for keeping the disc flared in the outer parts (Paper I). The size and temperature of the star changes in time (see Sect. 2.2) and the corresponding values are displayed in Table 1, where t gives the age of the star. The stellar mass is fixed to 1 M. We describe the effects of time on the size and temperature of the star in Sect. 2.2 and take these effects into account in Sect. 4.

Table 1

Stellar parameters for different as time evolves.

The initial surface density profile follows ΣGr− 15 / 14, which follows from Eq. (2) for a flared disc with H/rr2 / 7. We model different values by changing the underlying value of the surface density, while we keep α constant. This does not imply that the same viscosity (ν = αH2Ω) is present at all the different stages, because H changes as the disc evolves (Paper II). In our simulations we set the adiabatic index γ = 1.4 and used a mean molecular weight of μ = 2.3.

We use the opacity profile of Bell & Lin (1994), which is derived for micrometer-sized grains. In fact dust growth can be quite fast (Zsom et al. 2010), depleting the micro meter-sized dust grains to some level. However, larger grains (starting from mm size) only make a very small contribution to the opacity at high temperatures. At lower temperatures (T< 15 K), the larger mm grains will dominate the opacity, but those temperatures are not relevant within 50 AU. Additionally, if these grains start to grow and form larger pebbles, they will make a zero contribution to the opacity.

We define ΣZ as the surface density of heavy elements that are μm size in condensed form and ΣG as the gas surface density. Thus, the metallicity Z is the ratio Z = ΣZ/ ΣG, assumed to be independent of r in the disc. We assume that the grains are perfectly coupled to the gas, meaning that the dust-to-gas ratio is the same at every location in the disc. In our simulations we assume metallicities from 0.1% to 3.0%. This means that if grain growth occurs and the total amount of heavy elements (independent of size) stays the same, then the metallicity in our sense is reduced (as in Paper II).

2.2. Time evolution of discs

Protoplanetary accretion discs evolve with time and reduce their . Observations can give a link between mass accretion and time t. Hartmann et al. (1998) find a correlation between the mass accretion rate and the star age (4)This correlation includes stars in the Taurus star cluster that span an range from = 5 × 10-6 M/yr to = 5 × 10-10 M/yr over a range of 10 Myr. If the disc is accreting viscously, the evolution of the disc is directly proportional to the viscosity and hence to α. However, Eq. (4) was derived without any parametrization in α. We therefore consider that α in our simulations is not the time evolution parameter of the surface density, but simply a parameter for viscous heating (Q+α) and for determining vr. The time evolution of the disc is parametrized by Eq. (4), where we use the values without the error bars for our time evolution.

This approach also implies that the evolution time of the disc has to be longer than the time to relax to a radially constant state. This is more critical in the early evolution of the disc (high ), because the disc evolves more rapidly at that point. But the time the simulations need to settle in a steady state (constant at all r) is about a factor of a few (5) shorter than the disc evolution time in Eq. (4) for high values and is much shorter for the small ranges, validating our approach.

A disc with = 1 × 10-9 M/yr can be cleared by photo evaporation quite easily, so that the remaining lifetime is very short, only a few thousand years (Alexander et al. 2014). Therefore we do not simulate discs with very low . Recent observations have reported that even objects as old as 10 Myr can have accretion rates up to = 1 × 10-8 M/yr (Ingleby et al. 2014), which is in contrast to Eq. (4). Such high accretion rates at that age cannot be explained by viscous evolution models, unless the disc was very massive in the beginning of the evolution. Even if more of these objects are observed, this will not change the validity of the approach we are taking here. It would just change the time evolution of the disc presented in Eq. (4), but the disc structure as a function of (presented in Sect. 4) would stay the same.

During the millions of years of the disc’s evolution, the star evolves as well (Baraffe et al. 1998). As the star changes temperature and size, its luminosity changes (Eq. (3)), influencing the amount of stellar heating received by the disc. The stellar evolution sequence used was calculated by Baraffe et al. (1998) (see Table 1), where we display the stellar age, temperature, density, luminosity, and the corresponding accretion rate from Hartmann et al. (1998). In Fig. 1 we plot the stellar luminosity and the rate of the disc as a function of time. As the stellar luminosity drops by a factor of 3, the rate decreases by two orders of magnitude. Summarizing our methods, we have a disc model that features the full 2D structure with realistic heating and cooling that is linked to the temporal evolution of the star and disc.

thumbnail Fig. 1

Time evolution of the stellar luminosity after Baraffe et al. (1998) and the evolution of after Hartmann et al. (1998). The luminosity of the star reduces by a factor of 3 in 10 Myr and the accretion rate reduces by over 2 orders of magnitude during the same time period. However, an accretion rate of = 1 × 10-9 M/yr is already reached after 5 Myr when the disc can be cleared by photo evaporation.

Open with DEXTER

The total mass flowing through the disc can be evaluated by integrating the accretion rate specified in Eq. (4) over time. This gives a minimum estimate of the total mass of the disc, as leftover material is blown away by photo evaporation as soon as < 1 × 10-9 M/yr. During the disc’s lifetime of 5 Myr, a total of 0.05 M of gas flows through the disc, marking the total mass of the disc.

We therefore present fits to our simulations in Appendix A that can then be easily used by other studies that need a simple, but accurate time evolution model of the accretion disc structure.

3. Disc structure

In this section we compare the structure of a simulated disc with the MMSN and the CY2010 nebula to point out crucial differences in the disc structure and their effect on the formation of planetesimals and planetary embryos. The simulations in this section feature a metallicity of 0.5% in μm-sized dust grains. This value allows the disc to contain more heavy elements that could represent pebbles, planetesimals, or planetary cores (that do not contribute to the opacity) without increasing the total amount of heavy elements (grains and larger particles) to very high values.

3.1. disc

In Fig. 2 the temperature (top), the aspect ratio H/r (middle), and the surface density ΣG (bottom) are displayed. The MMSN and the CY2010 model follow power laws in all disc quantities. These are quoted in Table 2. The simulated = 3.5 × 10-8 M/yr disc model features bumps and dips in all disc quantities. More specifically, the simulation features a bump in T at the ice line (at Tice ≈ 170 K, illustrated in Fig. 2). As ice grains sublimate at higher temperature, the opacity reduces for larger T and therefore the cooling rate of the disc [D ∝ 1 / (κρ)] increases. This reduces the gradient of the temperature for high T compared to T<Tice, creating an inflection in T. This also changes the scale height H/r of the disc, which in turn influences the viscosity. Since the disc features a radially constant , the disc adapts to this change in viscosity by changing ΣG, creating the flattening of ΣG at the ice line (Paper II).

thumbnail Fig. 2

Mid-plane temperature T (top), aspect ratio H/r (middle) and integrated surface density ΣG (bottom) for a disc with = 3.5 × 10-8 M/yr and for the MMSN and the CY2010 model. The grey area marks in the temperature plot the region of the opacity transition at the ice line. The green line in the surface density plot indicates a fit for the surface density in the outer part of the disc. The main difference between the disc and the MMSN and CY2010 model is the higher temperature in the inner parts of the disc, which is caused by viscous heating. This increases H/r in the inner parts and then results in a reduced ΣG in the inner parts, owing to a change in viscosity that is compensated by a change in ΣG.

Open with DEXTER

In the inner parts of the disc, our models feature a much higher temperature than in the MMSN and the CY2010 models due to the inclusion of viscous heating. In the outer parts of the disc, where viscous heating becomes negligible and the temperature is solely determined by stellar irradiation, the temperatures of our simulations are comparable to the CY2010 model.

The aspect ratio H/r is related quite simply to the temperature (through the relation of hydrostatic equilibrium); (5)where is the gas constant, G the gravitational constant, and M the mass of the star. Therefore, the H/r plot (middle panel in Fig. 2) features the same properties as the temperature plot: i) a higher H/r in the inner disc for our simulations; and ii) about the same H/r in the outer disc compared to the CY2010 model. Strikingly, the aspect ratio of the MMSN is off by about 50% at 20 AU compared to our simulations. As the MMSN model does not feature viscous heating in the inner parts, H/r is smaller there, and because the radial change of H/r is roughly the same in both models in the outer parts of the disc, the MMSN has a much smaller H/r in the outer parts of the disc. The H/r diagram now features bumps and dips correlated to the wiggles in the temperature diagram. The dip in H/r starting beyond 3 AU represents a shadowed regions inside the disc. The stellar irradiation does not penetrate well into this region because it is absorbed by the bump in H/r at 3 AU. We emphasise that the drop in H/r is caused by the change of the cooling rate in the disc as the opacity changes and not by the heat transition from regions dominated by viscous heating to regions dominated by stellar heating (Paper II). If a constant opacity had been used, shadowed regions would not have appeared in the disc (Paper I). A drop in H/r normally implies outward migration for low mass planets (Papers I and II), a phenomenon that is not seen in the MMSN and CY2010 models.

The surface density profile (bottom panel in Fig. 2) of our simulations shows an inflection at the same location as where H/r shows a bump. Generally our simulations feature a lower surface density in the inner parts of the disc than in the MMSN and the CY2010 nebula. In the outer parts, on the other hand, the surface density of our simulations is much higher. This difference in surface density between our simulations and the other two models is caused by our underlying approximation (see Eq. (2)), which gives a shallower surface density slope for the outer disc. Andrews et al. (2010) find that the gradients of the surface density profile of accretion discs in the Ophiuchus star forming region are between 0.4 and 1.1. Our simulations feature a surface density gradient of 1 in the outer parts of the disc, but it is much shallower in the inner parts, matching the observations in contrast to the MMSN and CY2010 model.

3.2. Influence on planet formation

The streaming instability can lead to formation of planetesimals as a first step in forming planetary embryos and planets (Johansen et al. 2007; Johansen & Youdin 2007; Bai & Stone 2010a,b). The important quantity for triggering particle concentration by the streaming instability is Δ, which is the difference between the azimuthal mean gas flow and the Keplerian orbit divided by the sound speed. This difference is caused by the reduction of the effective gravitational force by the radially outwards pointing force of the radial pressure gradient. The parameter Δ is given by (6)where is the Keplerian velocity and cs/vK = H/r. Here η represents a measure of the gas pressure support (Nakagawa et al. 1986). The sub-Keplerian rotation of the gas makes small solid particles drift towards the star.

In Fig. 3 the radial pressure gradient ln(P) /ln(r) (top) and Δ (bottom) are displayed. The radial pressure gradient is constant for the MMSN and the CY2010 model, since these models are built upon strict power laws. The simulations, on the other hand, feature bumps and dips that result from the opacity transition. In the inner parts of the disc, the pressure gradient is much more shallow in the disc. We recall here that a shallower negative pressure gradient is preferred in order to form planetesimals.

Table 2

Parameters of the power laws used in the MMSN and the CY2010 models.

thumbnail Fig. 3

Radial pressure gradient (top) and pressure support parameter Δ = Δv/cs (bottom) for the disc with = 3.5 × 10-8 M/yr and for the MMSN and CY2010 models. The grey area indicates the radial range of the opacity transition at the ice line, as already indicated in Fig. 2. In contrast to the MMSN and CY2010 models that have a constant pressure gradient and therefore a steadily increasing Δ parameter, the model features a dip in the profile at 5 AU, which makes the formation of planetesimals in this region more likely. The reduced Δ parameter in the outer parts of the disc also makes it more likely to form planetesimals there compared to the MMSN and CY2010 model. The horizontal green lines at Δ = 0.025 and Δ = 0.05 mark the amount of heavy elements in the disc needed for the streaming instability to operate (1.5% and 2%, respectively), see Bai & Stone (2010b).

Open with DEXTER

The Δ parameter (bottom in Fig. 3) is nearly constant for our simulation in the inner parts of the disc, while it steadily decreases towards the star in the MMSN and CY2010 models. In the outer parts, Δ is up to 50% smaller in the model than in the MMSN, indicating significantly different conditions for forming planetesimals.

A reduced Δ parameter significantly helps the formation of large clumps via the streaming instability (Bai & Stone 2010b). Additionally, the formation of clumps is also dependent on the amount of heavy elements and particle sizes in the disc, where a larger amount of heavy elements strongly increases the clumping. The number of heavy elements is not restricted to the metallicity defined with μm dust grains, but also includes larger grains and pebbles that do not contribute to the opacity profile of the disc. For Δ = 0.025 (lower horizontal line in Fig. 3), a fraction of 1.5% in heavy elements is needed to form large clumps, while already for Δ = 0.05 (top horizontal green line in Fig. 3), a fraction of 2% in heavy elements is needed. For Δ = 0.1 a very high, probably not achievable, number of heavy elements is needed for the streaming instability to work. A reduction in Δ in the outer parts of the disc, as proposed by our model, therefore makes the formation of planetesimals much easier in the Kuiper belt.

Planetesimals can grow further by mutual collisions (Levison et al. 2010) or by the accretion of pebbles. In the latter case, core growth enters the fast Hill regime, when it reaches the “transition mass”, (7)where M is the stellar mass (Lambrechts & Johansen 2012). This corresponds to 0.03 MEarth at 5.2 AU in our simulated disc. A reduced disc scale height and Δ parameter thus help smaller embryos to reach this growth regime where cores are formed on time scales of 105 yr, even at wide orbits (50 AU), provided the pebble surface density is similar to MMSN estimates. Furthermore, a lower Δ makes accretion more efficient by increasing the proportion of pebbles that are accreted by a core versus those particles that drift past (Ormel & Kobayashi 2012).

Once planetary embryos have formed, these embryos are subject to gas-driven migration (for a review see Baruteau et al. 2014). The migration rate of planets can be determined by 3D simulations of protoplanetary discs (Kley et al. 2009; Bitsch & Kley 2011). However, these simulations are quite computationally expensive. Instead, one can calculate the torque acting on an embedded planet from the disc structure (Paardekooper et al. 2010, 2011). The torque formula of Paardekooper et al. (2011) takes torque saturation due to thermal diffusion into account and matches the 3D simulations of planets above 15 MEarth well (Bitsch & Kley 2011). However, for low mass planets (MP ≤ 5 MEarth), there is still a discrepancy between the torque formula and 3D simulations that actually show a more negative torque (Lega et al. 2014). Nevertheless, the predictions from the torque formula can already give first clues about the planetary migration history in evolving accretion discs.

In Paardekooper et al. (2011) the total torque acting on an embedded planet is a composition of its Lindblad torque and its corotation torque (8)The Lindblad and corotation torques depend on the local radial gradients of surface density ΣGrs, temperature Trβ, and entropy Srξ, with ξ = β − (γ − 1.0)s. Very roughly said, for not too negative ΣG gradients, a radially strong negative gradient in entropy, caused by a large negative gradient in temperature (large β), will lead to outward migration, while a shallow gradient in entropy will not lead to outward migration.

thumbnail Fig. 4

Torque acting on discs with = 3.5 × 10-8 M/yr. The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K. The region of outward migration is exactly correlated to the region in the disc, where the temperature gradient is the steepest and can trap planets between 5 and 30 Earth masses.

Open with DEXTER

In Fig. 4 the migration maps for the = 3.5 × 10-8 M/yr disc is displayed. The torque parameter Γ0 is defined as (9)where q is the planet-to-star mass ratio, h = H/r, r the semi-major axis and ΩP the Keplerian rotation. The actual speed of inward migration changes with planetary mass not only because the torque is proportional to the planetary mass squared, but also because the mass changes the saturation effects of the corotation torque (Paardekooper et al. 2011).

The regions of outward migration correspond to the regions in the disc where H/r decreases with r. The MMSN and the CY2010 model feature a flared disc for all r (Fig. 2). This imposes a very shallow temperature gradient leading to a shallow entropy gradient, which is not enough to produce a corotation torque that can overcompensate the Lindblad torque. Planets, regardless of mass, therefore migrate inwards towards the star in both the MMSN and CY2010 models. We do not display these migration maps here. This lack of zones of outward migration makes the formation of giant planet cores much harder.

4. Influence of the time evolution of disc and star

The protoplanetary disc gets accreted onto the star over millions of years, but the star also evolves on these time scales (Baraffe et al. 1998), contracting and changing its size and luminosity. The luminosity in turn determines the stellar heating that is absorbed by the disc, and thus influences the disc structure. By using Eq. (4), we can link different stages to different times. These times then give us the stellar evolution time and thus the stellar radius and temperature (Table 1). The stellar temperature stays roughly constant with time, while the stellar radius becomes significantly smaller as time evolves. This means that for discs with lower , the stellar heating will decrease as well. This will influence the outer regions of the disc, which are dominated by stellar heating. An increase of a factor of 2 in stellar luminosity results in a change of H/r of up to 20% and of up to 45% in temperature in the parts dominated by stellar irradiation. In this section all shown simulations feature a metallicity of 0.5%, as in Sect. 3.

In Fig. 5 the mid-plane temperature (top), H/r (middle), and the surface density ΣG (bottom) are displayed for different values of ranging from = 1 × 10-7 M/yr to = 4.375 × 10-9 M/yr. The temperature in the inner parts of the disc drops as decreases, because the viscous heating decreases as ΣG decreases. In the outer parts of the disc, the temperature also drops, because the stellar irradiation decreases in time. However, this drop in temperature is not as large as in the inner parts of the disc. In the late stages of the disc evolution, the disc becomes so cold that ice grains will exist throughout the disc so that no opacity transition is visible any more. For low there is a temperature inversion at approximately 5 AU. This inversion is caused by the different vertical heights of the absorptions of stellar photons compared to the vertical height of the cooling through diffusion.

thumbnail Fig. 5

Mid-plane temperature (top), H/r (middle), and surface density ΣG (bottom) for discs with different values of around an evolving star. The grey area in the temperature plot marks the transition range in temperature for the opacity law at the ice line. The black lines mark the fits discussed in Appendix A. For lower rates, the inner parts of the disc are colder, as viscous heating is reduced. For very low , the temperature is below the ice condensation temperature throughout the disc’s mid plane, and therefore the bump in the inner part of the protoplanetary disc vanishes.

Open with DEXTER

As evolves, the shadowed regions of the disc (seen by a reduction in H/r as a function of r) shrink. For very low , no bumps in H/r exist any more, because the temperature of the disc is so low that the opacity transition is no longer inside the computed domain, but further inside at r< 1 AU. This will have important consequences for planet migration, since outward migration is only possible when radially strong negative temperature gradients exist.

The wiggles in the temperature profile directly translate to dips in the surface density profile, because a change in the viscosity of the disc must be directly compensated for by a change in the disc’s surface density (Eq. (2)) to maintain a constant . In the very late stages of the disc evolution, when the accretion rate and the surface density are very low ( ≤ 1 × 10-9 M/yr), the disc will experience rapid photo-evaporation (Alexander et al. 2014).

thumbnail Fig. 6

Δ parameter for the disc simulations with evolving . The black lines mark the fits discussed in Appendix A. The Δ parameter describes the triggering of particle concentrations in the streaming instability, which can lead to planetesimal formation. In the regions of lower Δ at 5 AU, the formation of planetesimals is more likely, compared to the inner regions of the disc (2 AU) where Δ is higher.

Open with DEXTER

In Fig. 6 the Δ parameter (Eq. (6)) for the discs with different is displayed. In the inner parts of the disc, Δ only slightly reduces as shrinks, as long as is still high enough that the inner parts of the disc are dominated by viscous heating. As soon as the disc starts to become dominated by stellar heating in the inner parts, Δ drops by a significant factor. This is because Δ is proportional to H/r, which shows exactly that behaviour as well. This reduction of Δ for low helps the formation of planetesimals significantly, since lower metallicity is needed to achieve clumping (Bai & Stone 2010b), indicating that in the very late stages of the disc evolution, planetesimal formation becomes easier. In the outer parts of the disc, Δ is very high throughout the different stages and changes only slightly, following the reduction in H/r as the stellar luminosity decreases. This indicates that for the formation of planetesimals, very high metallicity is needed in the outer disc.

The growth of planetary embryos via pebble accretion is significantly accelerated, when the embryo reaches the Hill accretion regime, defined in Eq. (7). Through all stages, the minimum of H/r and Δ coincides, reducing the transition mass towards the Hill accretion regime and making planet formation in these locations easier. Additionally, as H/r and Δ drop in the late stages of evolution, pebble accretion is much more efficient because the transition mass is reduced, making core growth more efficient in the late stages of the disc, provided that enough pebbles are still available in the disc.

thumbnail Fig. 7

Gravitational torque acting on planets in discs with = 7 × 10-8 M/yr (top) and = 8.75 × 10-9 M/yr (bottom). The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K. The white line represents the zero-torque line of the fits presented in Appendix A. As drops, the regions of outward migration shrink so that only planets with lower mass can be saved from inward type-I-migration. Additionally, the orbital distance at which outward migration acts becomes smaller with decreasing . This is caused by the shallower gradient in temperature for lower discs.

Open with DEXTER

In Fig. 7 the migration map for the = 7 × 10-8 M/yr (top) and = 8.75 × 10-9 M/yr disc are displayed. The migration map for = 3.5 × 10-8 M/yr can be found in Fig. 4. As the disc evolves and decreases, the area of outward migration shrinks and moves inwards. This is caused by the disc becoming colder, and therefore the region of opacity transition at the ice line moves inwards as well. This implies that the strong radial negative gradients in temperature also move inwards, shifting the regions of outward migration to smaller radii (Paper II). In the late stages of the disc evolution, not only is the region of outward migration shifted towards the inner regions of the disc, but it also seems that outward migration is only supported for lower mass planets than for earlier disc evolution stages (high ). This was already observed in Paper II.

thumbnail Fig. 8

Mid-plane temperature for discs with various and metallicity. The metallicity is 0.1% in the top panel, 1.0% in the second from top panel, 2.0% in the third from top panel, and 3.0% in the bottom panel. An increasing metallicity reduces the cooling rate of the disc, making the disc hotter.

Open with DEXTER

5. Influence of metallicity

As the disc evolves in time, the micro-metre dust grains can grow and form larger particles, which (above mm size) do not contribute to the opacity of the disc any more and the opacity of the disc decreases, as the number of micro-metre dust grains reduces. In the later stages of the disc, the number of small dust grains can be replenished because of destructive collisions between planetesimals and larger objects, which would increase the opacity again. We therefore extend our model to range over a variety of metallicities, namely from 0.1% to 3.0%.

The temperature for discs with metallicities of 0.1%, 1.0%, 2.0%, and 3.0% (from top to bottom) are displayed in Fig. 8. The H/r and surface density profiles of those discs are shown in Appendix A. Clearly, discs with less metallicity are colder than their high-metallicity counterparts. This is caused by the increased cooling for low-metallicity discs, because D ∝ 1 / (ρκR). The general trends of the disc structure, however, hold for all metallicities, namely, that all discs feature a bump in temperature at the opacity transition at the ice line.

However, discs with < 1 × 10-8 M/yr and a metallicity of 0.1% no longer show a bump in temperature, and they follow the perfect flaring disc, because they are so cold that the transition in opacity at the ice line is no longer present in the outer parts of the disc. Additionally, the temperature inversion observed in Sect. 4 occurs at higher rates for low-metallicity discs than for discs with higher metallicity.

Discs that have the same value, but different metallicity do not have the same surface density profile (see Fig. A.3). As the opacity increases, the disc becomes hotter because it cools less efficiently, and therefore H increases, which then results in higher viscosity, reducing ΣG as is constant in r (Eq. (A.1)). An increase in metallicity therefore does not result in an increase of the same factor in temperature. The opposite applies for a lower metallicity, which results in a cooler disc (smaller H), decreasing the viscosity and therefore increasing ΣG.

The increase in temperature is much greater in the inner parts of the disc than in the outer parts. The reason for that lies in the fact that the inner disc is dominated by viscous heating, which does not depend on the opacity, and only the change in the cooling rate depends on the opacity and therefore changes the temperature in the inner disc. In the outer disc, however, the cooling and the heating both depend on the opacity. A reduced opacity increases the cooling, making the disc thinner (H/r decreases), but at the same time, the absorption of stellar photons is also reduced in the upper layers because fewer dust grains are available. This means that stellar irradiation can be efficient at much smaller heights from the mid plane, heating the disc at a smaller vertical height z, shifting the heated region closer to the mid plane, and making the disc hotter. This effect counterbalances the increased cooling, so that the disc stays roughly at the same temperature in the parts that are dominated by stellar irradiation for all metallicities, as long as the changes in metallicity are not too big.

The higher temperature in the inner disc, especially for lower , caused by higher metallicity has important implications for planetesimal formation, planetary migration, and the evolution of the ice line in time. These are discussed in Sect. 6.

6. Discussion

Planets will form during the evolution of the protoplanetary accretion disc. We discuss here the general implications of the disc structure on the formation of planetesimals and embryos (Sect. 6.1) and giant planets (Sect. 6.2). Additionally, we focus on the inward motion of the ice line as the disc evolves in time (Sect. 6.3).

6.1. Planetesimal and embryo formation

The first building blocks of planetary embryos are planetesimals, which can be formed by the streaming instability (Johansen & Youdin 2007). The streaming instability requires not only a small Δ parameter (Fig. 6), but also an increased amount of heavy elements (Bai & Stone 2010b). Both of these requirements are satisfied at the ice line, where an increased amount of heavy elements is likely to be present owing to the condensation of ice grains and pebbles (Ros & Johansen 2013). Additionally this region features a low Δ value, because it is in the shadowed region of the disc, the minimum of H/r. This feature is independent of the overall metallicity of the disc (Sect. 5), but because the disc is hotter for higher metallicity discs, the Δ parameter is also larger, as H/r increases. This implies that for higher metallicity discs, the streaming instability might not be able to be operated, because H/r is greater. However, the formation of the first planetesimals and planetary embryos is still most likely at the location of the ice line, since the disc features a minimum of Δ in that region. In the inner parts of the disc, the Δ parameter is also smaller. But since the disc is very hot there, the icy particles have evaporated and fewer metal grains are thus available for triggering the streaming instability.

If indeed the first planetesimals form at the ice line, then embryo formation is most likely to occur here as well. These embryos can then grow via collisions between each other and through pebble accretion (Lambrechts & Johansen 2012). The growth rate versus the migration rate is now key to determining what kind of planet emerges. If the planet grows rapidly to several Earth masses, then it can overcome the inward type-I-migration, be trapped in a region of outward migration (Lambrechts & Johansen 2014), and eventually become the core of a giant planet at a large orbital radius. If the planet does not grow fast enough to be trapped in a region of outward migration, the planet migrates inwards and can end up as either a giant planet in the inner system (if the core can grow fast enough in the inner parts of the disc) or “just” be a hot super-Earth or mini-Neptune (Cossou et al. 2014).

As the disc evolves in time, the Δ parameter becomes smaller in the inner disc, meaning that the formation of planetesimals might be more likely at later times. However, the farther the disc has evolved in time, the less time is left for giant planets to form, because those have to accrete gas from the surrounding disc. Planetesimal formation might also be triggered indirectly by grain growth. As the metallicity of μm dust grains drops, e.g. due to grain growth, H/r decreases, because of a drop in opacity, which increases the cooling of the disc. This then leads to a lower Δ value, making it easier for the streaming instability to operate. These two effects would imply that the formation of planetesimals is easier in the late stages of the disc evolution, thus potentially explaining the abundance of small planets compared to gas giants (Fressin et al. 2013), which require an earlier core formation in order to accrete an gaseous envelope.

In the outer parts of the disc (r> 10 AU), the Δ parameter of the discs, regardless of the discs metallicity, is up to 50% smaller than in the MMSN and CY2010 model. A reduction in Δ by this factor significantly helps the clumping of particles in the streaming instability that can then lead to the formation of planetesimals, making the formation of planetesimals much more likely in our presented model. This makes it much easier to form Kuiper belt objects and the cores of Neptune and Uranus in our model. During time, as ΣG reduces, the Δ parameter even reduces slightly in the outer disc, making the formation of planetesimals more likely in the later stages. Such late formation of planetesimals in the outer disc could explain why Neptune and Uranus did not grow to become gas giants (Lambrechts et al. 2014; Lambrechts & Johansen 2014).

A full analysis of the effect of the disc structure on planetary formation and migration will require introducing detailed prescriptions for planetary growth. This will be the topic of a future paper.

6.2. Giant planet formation

For solar-like stars, the occurrence of a giant planet (MP> 0.1 MJup) is related to the metallicity of the host star. In particular, higher metallicity implies a higher occurrence rate of giant planets (Santos et al. 2004; Fischer & Valenti 2005). This implies that the formation of planets in systems with higher metallicity might be systematically different from systems with lower metallicity.

In the core accretion scenario, a core forms first, which then accretes gas to become a giant planet. The core itself migrates through the disc as it grows. Cossou et al. (2014) show that planets that are not trapped inside regions of outward migration are more likely to become super-Earths rather than giant planets, making zones of outward migration important for the final structure of planetary systems.

The mid-plane temperature profiles for discs with different and metallicity are shown in Fig. 8. They imply that higher metallicity results in a hotter disc that can keep the shadowed regions longer (see also Fig. A.2), which are able to support outward migration at a few AU for a longer time compared to lower metallicity discs. This is illustrated in Fig. 9, which shows the migration regions where outward migration is possible for discs with different metallicity and for two values; = 1.0 × 10-7 M/yr and = 8.75 × 10-9 M/yr.

thumbnail Fig. 9

Migration contours for discs with different metallicities. Top panel: = 1.0 × 10-7 M/yr; bottom panel: = 8.75 × 10-9 M/yr. Planets migrate outwards in the regions of the disc that are enclosed by the solid lines, where the colours mark different metallicities.

Open with DEXTER

For the = 1.0 × 10-7 M/yr discs, there are three clear trends visible for higher metallicities:

  • the regions of outward migration are shifted farther away from the central star, because the disc is hotter and therefore the opacity transition at the ice line is farther away;

  • the regions of outward migration are larger in radial extent, because the shadowed regions in the disc are larger, which enlarges the region where a steep radial temperature gradient exits, in turn enlarging the entropy driven corotation torque;

  • outward migration is only possible for higher minimal masses, which is caused by the changes in the disc structure that overcompensate for the reduced cooling time caused by the higher metallicity, which would actually reduce the minimal mass required for outward migration.

In principle one could also infer from Fig. 9 that outward migration is also possible for higher planetary masses for increasing metallicity. However, as the torques acting on the planets are calculated by a torque formula that was derived in the linear regime for low mass planets (5 MEarth), one cannot extend it towards masses of a few 10 Earth masses. Additionally planets that are that massive start to open up gaps in the disc, which means that the planet transitions into type-II-migration.

The regions of outward migration then shrink and move closer towards the star as the disc evolves in time and reduces. For the = 8.75 × 10-9 M/yr disc, the discs with high metallicity (>2.0%) maintain two large regions of outward migration at a few AU, which are able to trap planets of up to 30 MEarth. The lower metallicity discs, on the other hand, only feature one region of outward migration in the inner region of the disc, which is only able to trap cores up to 15 MEarth.

Additionally, as the regions of outward migration last longer in the outer disc for high metallicity, the growing core can stay farther out, being released at a larger orbital distance into type-II-migration when it becomes large enough. This core would then migrate inwards in the type-II regime and finally be stopped when the disc dissipates (Alexander & Pascucci 2012), which might explain the pile-up of Jupiter-sized planets around 1 AU. On the other hand, in a low metallicity disc, the planetary core would be trapped closer to the star, and when it is released in type-II-migration as a gas giant it is already closer to the star and might therefore become a hot Jupiter.

6.3. Ice line

As the disc evolves in time and drops, the inner regions of the disc become colder, moving the ice line closer to the star. In the evolution of all our disc models, the ice line (at 170 K) moves to 1 AU at 2 Myr. However, this result is troublesome considering evidence from meteors and asteroids in our own solar system. Ordinary and enstatite chondrites contain very little water and must have formed on the warm side of the ice line. The parent bodies formed 2–3 Myr after CAIs (which mark the zero age of the solar system). However, at this time the ice line in our nominal models had moved to approximately 1 AU. This is potentially in conflict with the dominance of S-type asteroids, believed to be the source of ordinary chondrites in the inner regions of the asteroid belt. We propose two different ideas that could potentially add to the solution of this problem.

The metallicity is the key parameter for the temperature profile of the disc. Our simulations show that a higher metallicity in μm dust increases the temperature of the inner disc. If the metallicity is even higher (e.g. 5%) the ice line would have been farther outside. Therefore an intrinsic higher metallicity of the disc could help this problem. Alternatively, a higher metallicity in μm-sized dust grains can be created by dust-polluted ice balls that drift across the ice line and then release their dust grains as the ice melts, increasing the metallicity in μm-sized dust in the inner parts of the disc (Sirono 2011).

The time evolution of the accretion disc contains large error bars in time (Hartmann et al. 1998), which could simply imply that the solar system was one of the slower evolving discs, keeping a higher rate at later times, which would keep the ice line farther out and thus the inner system dry at later stages of the disc evolution. For example, in Hartmann et al. (1998) one observed disc still has ≈ 2.0 × 10-8 M/yr at 3 Myr, which has a high enough temperature to keep the ice line at 2 AU for Z ≥ 0.005. Ingleby et al. (2014) report discs with high accretion rates ( = 1 × 10-8 M/yr) that are 10 Myr old, which would then have an ice line farther out and thus avoid the mentioned problem.

However, it is unlikely that these simple scenarios are solely responsible for keeping the snow line far out during the lifetime of the disc, and more complicated scenarios can play a role. This includes local heating from shocks in disc with dead zones (Martin & Livio 2013) or the interplay among the radial motion of the gas, drifting icy particles, and growing planets, which will be the subject of a future paper.

7. Summary

In this work we have compared the power law assumptions of the minimum mass solar nebula and the Chiang & Youdin (2010) models with simulations of protoplanetary discs that feature realistic radiative cooling, viscous, and stellar heating. The modelled disc structures show bumps and dips that are caused by transitions in the opacity, because a change in opacity changes the cooling rate of the disc (Bitsch et al. 2014a). These features can act as sweet spots for forming planetesimals via the streaming instability and for stopping the inward migration of planetary cores of a few Earth masses. Regions of low pressure support also enhance the growth of planetary cores via pebble accretion. These attributes are lacking in the MMSN and CY2010 models, making the formation of planets in these models much harder. Additionally, as the disc evolves in time and the accretion rate decreases, the radial gradients of ΣG, H/r and T in the disc change. This temporal evolution is not taken into account in the MMSN and CY2010 models either.

During the time span of a few Myr, the star evolves along with the disc. The star changes its size and temperature and therefore its luminosity. This changes the amount of stellar heating received by the disc and strongly changes the parts of the disc where stellar heating is the dominant heat source. For higher this mostly affects the outer parts of the disc, while for lower , as viscous heating becomes less and less important, the whole disc structure is affected.

We present a simple fit of our simulated discs over all accretion rates and for different metallicities. The fit consists of three parts that correspond to three different regions in the disc, which are dominated by different heat sources. The inner disc is dominated by viscous heating, while the outer disc is dominated by stellar irradiation. In between, a transition region in the disc exists where viscous heating starts to become less important, but is in the shadow of direct stellar illumination at the same time.

The different values can be linked to a time evolution of the disc obtained from observations (Hartmann et al. 1998). With this simple relation between time and accretion rate, our presented fit can easily be used to calculate the disc structure at any given evolution time of the star-disc system. This can then be used as input for planet formation models.

A Fortran script producing T, Σ, and H/r as functions of Z, and time is available upon request.

Acknowledgments

B.B., A.J., and M.L. thank the Knut and Alice Wallenberg Foundation for their financial support. A.J. was also supported by the Swedish Research Council (grant 2010-3710) and the European Research Council (ERC Starting Grant 278675-PEBBLE2PLANET). A.M. is thankful to the Agence Nationale pour la Recherche under grant ANR-13-BS05-0003-01 (MOJO). The computations were done on the “Mesocentre SIGAMM” machine, hosted by the Observatoire de la Côte d’Azur. We thank M. Havel for discussing with us on the early evolution model of the Sun.

References

Online material

Appendix A: A simple model for evolving accretion discs

In this Appendix we provide a fit for our simulated disc models for all the different stages. We provide here only a fit for the temperature T, because for discs with constant , all other quantities can be derived, from the value of and the α parameter. This is done by using the relation and the hydrostatic equilibrium equation (A.1)where is the gas constant. The temperature given in the following fit is the mid-plane temperature of the disc, which does not represent the vertical structure of the disc that can feature a super-heated upper layer (e.g. Paper I). With the fit of the temperature, H can be determined and then finally ΣG.

We identify three different regimes for the temperature gradients in the disc. The inner disc is dominated by viscous heating and features a higher temperature than the water condensation temperature (Fig. A.1), while the outer disc is dominated by stellar irradiation. In between, the disc is in the shadow of stellar irradiation. This can be seen more clearly in the middle panel of Fig. 5, where H/r drops with increasing r. We fit each of these three regimes with individual power laws, which are indicated in Fig. A.1 as well.

thumbnail Fig. A.1

Radial temperature profile of the disc with = 3.5 × 10-8 M/yr. The different background colours mark different fit-regimes. The grey area marks the inner part of the disc that is dominated by viscous heating, while the golden part indicates the outer parts of the disc that is dominated by stellar heating. The white middle region indicates the shadowed region of the disc.

Open with DEXTER

The fit for the temperature for the outer part of the disc is Tr− 4 / 7 which is steeper than in the MMSN and CY2010 models that are based on analytical arguments (e.g. Chiang & Goldreich 1997). The reason for this originates in the shadowed parts of the disc. The disc is cooler in the shadow and is heated by diffusion from the hotter inner part (heated viscously) and by the hotter outer part (heated from the star). This leads to a steeper temperature gradient in the outer parts, where the disc just emerges from the shadow, since it has to (partially) compensate for the lack of heating in the shadowed region. This compensation for heat influences the gradient in temperature to quite large distances. But we suspect that the gradient in temperature will return back to Tr− 3 / 7 in the very outer parts of the disc, because the influence of the shadowed region reduces with larger orbital distances. This is also true for discs in the later stage evolution, since those discs no longer feature a prominent bump at the ice line that can shade parts of the disc from stellar irradiation.

The changes in the temperature gradient originate in the changes of opacity, which changes the cooling rate of the disc and hence the temperature (Paper II). The changes in opacity occur at the ice line, where ice grains melt and sublimate, and at the silicate line (1000 K). As decreases, the viscous heating in the inner disc decreases as well, making the disc colder and moving the ice-line bump closer towards the central star. If the temperature of the disc is always below the water condensation temperature, the disc will therefore not feature the ice line bump any more (see < 8.75 × 10-8 M/yr in Fig. 5), and thus the inner part of the fit, which requires a higher temperature than the water condensation temperature, is not needed any more. Consequently, the fit of the temperature profile will only consist of fits of the middle and outer parts of the disc (indicated in white and yellow in Fig. A.1).

As long as the temperature is below the sublimation temperature of silicate grains (T< 1000 K), the fit can be applied to regions that are inside 1 AU. The reason for that is that there are no additional bumps in the opacity profile for 170 K <T< 1000 K that could change the disc structure by changing the cooling rate of the disc (Paper I).

As decreases and the viscous heating in the inner parts of the disc decreases, the disc’s heating is dominated more and more by stellar irradiation. This transition in the dominant heating process occurs for ≈ 1.75 × 10-8 M/yr. We therefore observe a huge drop in the disc’s temperature (top panel in Figs. 5 and 8) around these rates. For even lower rates, the changes in temperature seem to be smaller again, which is caused by the fact that the disc’s heating is now entirely dominated by stellar irradiation. We therefore introduce different fitting regimes for different values.

Since the heating and cooling of the disc is directly proportional to the value, our fit is a function of . We express each of the three individual power law fits in the following form (A.2)where A is a reference temperature expressed in K and B the factor by which the scaling with is expressed through μR = log (/ref), where ref is the reference value. The factor ψij expresses the transition from regime i to regime j. The factor l corresponds to the power law index of the fit, shown in Fig. A.1.

In the inner regions of the disc for high , the disc is dominated by viscous heating, and higher metallicity reduces the cooling, making the disc hotter, therefore increasing the temperature and the gradient of the temperature, which is then a function of metallicity. For < 1.75 × 10-8 M/yr, the disc starts to become dominated by stellar irradiation, so that the temperature gradient in the inner parts of the disc flattens. At the same rates, the shadowed regions of the disc shrink, making the gradient of temperature in the outer parts of the disc flatter as well, because the outer parts of the disc no longer have to compensate for such a large unheated shadowed region.

The fit for the whole disc with ≥ 1.75 × 10-8 M/yr can therefore be expressed as (A.3)For < 1.75 × 10-8 M/yr and ≥ 8.75 × 10-9 M/yr, the fit is expressed as (A.4)where the gradient of the inner part of the disc relaxes back to r− 6 / 7 at = 8.75 × 10-9 M/yr, and the gradient of the outer part of the disc becomes less steep as a function of Z. For < 8.75 × 10-9 M/yr the fit is expressed as (A.5)These formula are used for both high and low metallicity discs.

The middle part contains two ψ functions since the second regime has to be smoothed into the first and into the third regime. For all different fitting regimes the values of A, B, and ψ change, but keep the same shape.

Along with the changes in , the disc changes as different metallicities provide different cooling rates and therefore different temperatures. These changes in temperature then reflect back on changes in H/r and ΣG. We therefore also introduce different regimes for the metallicity, where we distinguish between high metallicity (Z ≥ 0.005) and low metallicity (Z< 0.005). The factors A, B, and ψ are therefore functions of and Z. The scaling with Z is expressed through χ = log  (100Z/ 0.5) / log  (2).

Appendix A.1: Low metallicity discs

In this section the fit for discs with Z< 0.005 is presented. The first regime of is for > 3.5 × 10-8 M/yr, which corresponds to the early evolution of the disc. We define ref = 3.5 × 10-8 M/yr, and the fit reads as (where r is in units of AU) (A.6)The next regime spans the interval 1.75 × 10-8 M/ yr < ≤ 3.5 × 10-8 M/yr and ref = 3.5 × 10-8 M/yr, and the fit reads as(A.7)The next regime spans the interval 8.75 × 10-9 M/ yr < ≤ 1.75 × 10-8 M/yr and marks the regime where T undergoes a large change in the inner parts of the disc, because viscous heating is replaced by stellar heating as the dominant heat source. Here, ref = 1.75 × 10-8 M/yr (which changes μR), which results in some reduction factors R that need to be taken into account: (A.8)The fit then reads as (A.9)For even lower values, the disc is dominated by stellar heating. We therefore have to introduce another regime that is ≤ 8.75 × 10-9 M/yr. The reduction factors are now (A.10)Here ref = 8.75 × 10-9 M/yr and the fit reads as(A.11)For accretion rates of 1 × 10-10 M/ yr << 1 × 10-8 M/yr, photo evaporation of the disc can be dominant and clear the disc on a very short time scale (Alexander et al. 2014). This implies that going to accretion rates of < 1 × 10-9 M/yr might not be relevant for the evolution of the disc. Therefore we recommend not using these fits for < 1 × 10-9 M/yr because the disc will be photo-evaporated away in a very short time for these parameters.

In Fig. 5 the mid-plane temperature (top), H/r (middle), and the surface density (bottom) of our simulations and the fits from Eqs. (A.6) to (A.11) with Z = 0.005 are displayed (shown in black lines). The general quantities of the discs are captured well by the provided fits. In Fig. 6 the Δ parameter (Eq. (6)) of the simulations with Z = 0.005 and the one derived from the fit is displayed. The general trend towards a decreasing Δ in the inner parts of the disc with decreasing is reproduced well. The pressure in the disc is calculated through its proportionality with the disc’s density (A.12)where cV is the specific heat at constant volume. In hydrostatic equilibrium the density is calculated through . The fit in Δ seems to be a bit off for < 4.375 × 10-9 M/yr. The reason for that lies in the different quantities that are needed to compute Δ. Small errors in independent variables (e.g. T, ΣG) translate into larger errors in the final quantity (here Δ).

In Fig. 7 the torque maps for the discs with = 7 × 10-8 M/yr (top) and = 8.75 × 10-9 M/yr (bottom) with Z = 0.005 are displayed. The white lines indicate the zero-torque lines resulting from the fitted disc structures presented by the black lines in Fig. 5. The overall agreement in the torque acting on the planet from the fit compared to the simulation is quite good, but there seem to be some variations on the “edges” of the region of outward migration. This is caused by the fit having small deviations from the simulations for all quantities of the disc. But in the migration formula of Paardekooper et al. (2011), all these quantities and their gradients play a role. This is especially important for the torque saturation at higher planetary masses. However, the general trend, especially for low mass planets, is captured well.

Appendix A.2: High metallicity discs

In this section fits for discs with high metallicity (Z ≥ 0.005) are presented. The fits follow the basic principles as for the low metallicity discs, but the parameters for A, B, and ψ differ, because the dependence on Z and therefore on χ changes. The first regime of for > 3.5 × 10-8 M/yr is represented by the following fit (where r is in units of AU) (A.13)The next regime spans the interval 1.75 × 10-8 M/ yr < ≤ 3.5 × 10-8 M/yr and ref = 3.5 × 10-8 M/yr and the fit reads as (A.14)The next regime spans the interval 8.75 × 10-9 M/ yr < ≤ 1.75 × 10-8 M/yr and marks the regime where T undergoes a large change in the inner parts of the disc, because viscous heating is replaced by stellar heating as dominant heat source. Here, ref = 1.75 × 10-8 M/yr (which changes μR), which results in some reduction factors R that need to be taken into account: (A.15)The fit then reads

thumbnail Fig. A.2

H/r for discs with various and metallicity. The metallicity is 0.1% in the top panel, 1.0% in the second from top panel, 2.0% in the third from top panel, and 3.0% in the bottom panel. An increasing metallicity reduces the cooling rate of the disc, making the disc hotter, and therefore increases H/r. The black lines in the plots indicate the fit presented in Appendix A.

Open with DEXTER

as(A.16)For even lower values, the disc is dominated by stellar heating. We therefore have to introduce another regime that is ≤ 8.75 × 10-9 M/yr. The reduction factors are now (A.17)Here ref = 8.75 × 10-9 M/yr and the fit reads as (A.18)The black lines in Fig. 8 represent the fits obtained with the formulas of Eqs. (A.13) to (A.18) for discs with different Z. In Figs. A.2 and A.3, the surface density and H/r of the same discs and the corresponding fits are displayed.

thumbnail Fig. A.3

Surface density ΣG for discs with various and metallicity. The metallicity is 0.1% in the top panel, 1.0% in the second from top panel, 2.0% in the third from top panel, and 3.0% in the bottom panel. An increasing metallicity reduces the cooling rate of the disc, making the disc hotter, increasing the viscosity, and therefore decreasing ΣG compared to discs with lower metallicity. The black lines in the plots indicate the fits presented in Appendix A.

Open with DEXTER

thumbnail Fig. A.4

Mid-plane temperature (top), H/r (middle), and surface density ΣG (bottom) for discs at different time evolutions. The values for T, H/r, and ΣG have been calculated from Eqs. (A.6) to (A.11) with a metallicity of 0.5%. As time evolves, the surface density reduces and with it the inner region that is dominated by viscous heating, so that the disc is dominated by stellar irradiation at the late stages of the evolution, after a few Myr. In this late stage, the aspect ratio is nearly constant in the inner parts of the disc and T and ΣG follow simple power laws.

Open with DEXTER

Appendix A.3: Time evolution of the disc

To link these fits with time evolution, we relate back to Eq. (4), which links time with a corresponding value. For each value, we can calculate the disc structure through Eqs. (A.6) to (A.11) for discs with Z ≤ 0.005 and through Eqs. (A.13) to (A.18) for discs with Z> 0.5. We can therefore calculate what the disc structure looks like at different times. This is presented in Fig. A.4, where temperature (top), H/r (middle), and surface density (bottom) of discs with Z = 0.005 and different ages are displayed.

The age of the disc, tdisc = 0 is not defined in Eq. (4). Here we take the age tdisc = 0 to be 100 kyr after the formation of the disc, because during the first stages, the disc is still undergoing massive changes in its structure (Williams & Cieza 2011). We therefore recommend using this definition of tdisc = 0 age because it avoids discontinuities in the disc evolution. The age of the star, presented in Table 1, follows a different time scale, meaning age tdisc = 0 for the disc corresponds to a stellar age of t = 100 kyr. This means we add 100 kyr in Eq. (4) to define the new tdisc = 0 age for the disc in the time of the star. Equation (4) then looks (without the measurement errors): (A.19)This equation describes the age of the disc tdisc in our set-up (indicated in Fig. A.4), while Eq. (4) describes the age of the star t as stated in Table 1.

All Tables

Table 1

Stellar parameters for different as time evolves.

Table 2

Parameters of the power laws used in the MMSN and the CY2010 models.

All Figures

thumbnail Fig. 1

Time evolution of the stellar luminosity after Baraffe et al. (1998) and the evolution of after Hartmann et al. (1998). The luminosity of the star reduces by a factor of 3 in 10 Myr and the accretion rate reduces by over 2 orders of magnitude during the same time period. However, an accretion rate of = 1 × 10-9 M/yr is already reached after 5 Myr when the disc can be cleared by photo evaporation.

Open with DEXTER
In the text
thumbnail Fig. 2

Mid-plane temperature T (top), aspect ratio H/r (middle) and integrated surface density ΣG (bottom) for a disc with = 3.5 × 10-8 M/yr and for the MMSN and the CY2010 model. The grey area marks in the temperature plot the region of the opacity transition at the ice line. The green line in the surface density plot indicates a fit for the surface density in the outer part of the disc. The main difference between the disc and the MMSN and CY2010 model is the higher temperature in the inner parts of the disc, which is caused by viscous heating. This increases H/r in the inner parts and then results in a reduced ΣG in the inner parts, owing to a change in viscosity that is compensated by a change in ΣG.

Open with DEXTER
In the text
thumbnail Fig. 3

Radial pressure gradient (top) and pressure support parameter Δ = Δv/cs (bottom) for the disc with = 3.5 × 10-8 M/yr and for the MMSN and CY2010 models. The grey area indicates the radial range of the opacity transition at the ice line, as already indicated in Fig. 2. In contrast to the MMSN and CY2010 models that have a constant pressure gradient and therefore a steadily increasing Δ parameter, the model features a dip in the profile at 5 AU, which makes the formation of planetesimals in this region more likely. The reduced Δ parameter in the outer parts of the disc also makes it more likely to form planetesimals there compared to the MMSN and CY2010 model. The horizontal green lines at Δ = 0.025 and Δ = 0.05 mark the amount of heavy elements in the disc needed for the streaming instability to operate (1.5% and 2%, respectively), see Bai & Stone (2010b).

Open with DEXTER
In the text
thumbnail Fig. 4

Torque acting on discs with = 3.5 × 10-8 M/yr. The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K. The region of outward migration is exactly correlated to the region in the disc, where the temperature gradient is the steepest and can trap planets between 5 and 30 Earth masses.

Open with DEXTER
In the text
thumbnail Fig. 5

Mid-plane temperature (top), H/r (middle), and surface density ΣG (bottom) for discs with different values of around an evolving star. The grey area in the temperature plot marks the transition range in temperature for the opacity law at the ice line. The black lines mark the fits discussed in Appendix A. For lower rates, the inner parts of the disc are colder, as viscous heating is reduced. For very low , the temperature is below the ice condensation temperature throughout the disc’s mid plane, and therefore the bump in the inner part of the protoplanetary disc vanishes.

Open with DEXTER
In the text
thumbnail Fig. 6

Δ parameter for the disc simulations with evolving . The black lines mark the fits discussed in Appendix A. The Δ parameter describes the triggering of particle concentrations in the streaming instability, which can lead to planetesimal formation. In the regions of lower Δ at 5 AU, the formation of planetesimals is more likely, compared to the inner regions of the disc (2 AU) where Δ is higher.

Open with DEXTER
In the text
thumbnail Fig. 7

Gravitational torque acting on planets in discs with = 7 × 10-8 M/yr (top) and = 8.75 × 10-9 M/yr (bottom). The black lines encircle the regions of outward migration. The vertical red lines indicate the ice line at 170 K. The white line represents the zero-torque line of the fits presented in Appendix A. As drops, the regions of outward migration shrink so that only planets with lower mass can be saved from inward type-I-migration. Additionally, the orbital distance at which outward migration acts becomes smaller with decreasing . This is caused by the shallower gradient in temperature for lower discs.

Open with DEXTER
In the text
thumbnail Fig. 8

Mid-plane temperature for discs with various and metallicity. The metallicity is 0.1% in the top panel, 1.0% in the second from top panel, 2.0% in the third from top panel, and 3.0% in the bottom panel. An increasing metallicity reduces the cooling rate of the disc, making the disc hotter.

Open with DEXTER
In the text
thumbnail Fig. 9

Migration contours for discs with different metallicities. Top panel: = 1.0 × 10-7 M/yr; bottom panel: = 8.75 × 10-9 M/yr. Planets migrate outwards in the regions of the disc that are enclosed by the solid lines, where the colours mark different metallicities.

Open with DEXTER
In the text
thumbnail Fig. A.1

Radial temperature profile of the disc with = 3.5 × 10-8 M/yr. The different background colours mark different fit-regimes. The grey area marks the inner part of the disc that is dominated by viscous heating, while the golden part indicates the outer parts of the disc that is dominated by stellar heating. The white middle region indicates the shadowed region of the disc.

Open with DEXTER
In the text
thumbnail Fig. A.2

H/r for discs with various and metallicity. The metallicity is 0.1% in the top panel, 1.0% in the second from top panel, 2.0% in the third from top panel, and 3.0% in the bottom panel. An increasing metallicity reduces the cooling rate of the disc, making the disc hotter, and therefore increases H/r. The black lines in the plots indicate the fit presented in Appendix A.

Open with DEXTER
In the text
thumbnail Fig. A.3

Surface density ΣG for discs with various and metallicity. The metallicity is 0.1% in the top panel, 1.0% in the second from top panel, 2.0% in the third from top panel, and 3.0% in the bottom panel. An increasing metallicity reduces the cooling rate of the disc, making the disc hotter, increasing the viscosity, and therefore decreasing ΣG compared to discs with lower metallicity. The black lines in the plots indicate the fits presented in Appendix A.

Open with DEXTER
In the text
thumbnail Fig. A.4

Mid-plane temperature (top), H/r (middle), and surface density ΣG (bottom) for discs at different time evolutions. The values for T, H/r, and ΣG have been calculated from Eqs. (A.6) to (A.11) with a metallicity of 0.5%. As time evolves, the surface density reduces and with it the inner region that is dominated by viscous heating, so that the disc is dominated by stellar irradiation at the late stages of the evolution, after a few Myr. In this late stage, the aspect ratio is nearly constant in the inner parts of the disc and T and ΣG follow simple power laws.

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.