EDP Sciences
Free Access
Issue
A&A
Volume 549, January 2013
Article Number A124
Number of page(s) 14
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201220159
Published online 10 January 2013

© ESO, 2013

1. Introduction

In accretion discs around young stars planetary embryos are born and grow through collisions to form bigger and bigger objects. The largest protoplanets constitute the cores of giant planets. Embedded planets interact with the gas in the disc and can move through the disc. This process depends on the disc’s physical properties (Ward 1997), in particular the local density and temperature and the gradients thereof.

In a real disc, not only viscous heating and radiative cooling determine the disc’s thermal structure, but also stellar irradiation. In the outer parts of the disc the heating from the star can keep the disc flared (Chiang & Goldreich 1997), in contrast to the disc profile without stellar irradiation. The effect of stellar irradiation on the disc structure has largely been investigated utilizing a 1D + 1D numerical approach (e.g. Bell et al. 1997; Dullemond et al. 2001) with the goal of fitting the spectral energy distributions (SEDs) of observed discs

The migration of a low-mass planet embedded in a fully radiative gaseous disc can be significantly different from migration in an isothermal or purely adiabatic disc (Paardekooper & Mellema 2006, 2008; Paardekooper & Papaloizou 2008; Kley & Crida 2008; Kley et al. 2009; Ayliffe & Bate 2010). While all authors agree that radiation transport can slow the rate of inward migration, there is still a lack of consensus whether the direction of migration can be outward and, if so, in which part of the disc (Bitsch & Kley 2011). Part of this confusion may be due to the sensitivity of the direction and magnitude of migration on local disc parameters (Paardekooper et al. 2010, 2011; Masset & Casoli 2010), including, for example, the radial disc temperature gradient (Ayliffe & Bate 2011).

In all previous works on planetary migration, the disc structure is determined by viscous heating and cooling. The equilibrium disc structure is thus determined by the disc mass and the magnitude of viscosity (Bitsch & Kley 2011). Including stellar irradiation can lead to a flared disc profile in the outer parts of the disc (Dullemond & Dominik 2004). The resulting change in the disc structure in the outer parts can have a huge effect on the migration of embedded planets. As the aspect ratio of the disc H/r changes, so does the temperature profile, which in turn determines the strength of the entropy gradient. The entropy gradient determines the magnitude of the corotation torque (Baruteau & Masset 2008), which constitutes a significant fraction of the total torque, and hence determines the strength and direction of migration.

In this paper we investigate the structure of the disc due to changes in opacity and the effects of stellar irradiation using a 2D (r − z) disc model. In Sect. 2 we describe the energy equation including stellar heating and the opacities used. We then compare test studies with analytical calculations of the disc structure in the case of constant opacity in Sect. 3, followed by studies with non constant opacity in Sect. 4. The differences between discs including stellar heating and those with only viscous heating are described in Sect. 4.3. We then apply torque formulae to estimate the torque acting on embedded planets for these disc structures in Sect. 5. In Sect. 6 the summary and conclusions are presented.

2. Methods

The protoplanetary disc is treated in this study as a two-dimensional (2D) non-self-gravitating gas whose motion is described by the Navier-Stokes equations. In this paper we focus on the axisymmetric structure of the disc and compute a vertical slice (in the r − z plane) utilizing 2D spherical coordinates (r − θ). Turbulence in discs is thought to be driven by magneto-hydrodynamical instabilities (Balbus & Hawley 1998) but here we treat viscosity utilizing either a constant kinematic viscosity or an α-viscosity (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 will be described in detail in Sect. 2.1. For that purpose we modified and substantially extended an existing multi-dimensional hydrodynamical code Nirvana (Ziegler & Yorke 1997; Kley et al. 2001).

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 that allows us to 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), so we focus here on our newly implemented two-temperature approach in Sect. 2.1 and how to solve the coupled equations in Appendix B.

2.1. Energy equation

In the two-temperature approach (featuring the radiative energy density ER and the thermal energy density ϵ) the evolution equations for the thermal and radiation energy read (Kley 1989; Dobbs-Dixon et al. 2010; Commercon et al. 2011): (1)Here, ER is independently evolved from the thermal component with B(T) = 4σT4 (σ being the StefanBoltzmann constant and T the temperature of the gas). F denotes the radiative flux. From Eq. (1), one can calculate the coupling time-scale τcoup = 1/(ρκc) with the density ρ, the opacity κ, and c the speed of light. The coupling time-scale defines the coupling between B(T)/c and ER. If ρκ is high, the coupling time is small, indicating that ER relaxes quickly towards B(T)/c, meaning ER ≈ 4σT4/c (optically thick limit). If ρκ is low, the coupling time is large, and ER and B(T)/c can decouple (optically thin limit). B(T) − cER represents the exchange of energy between the thermal and radiative components through emission and absorption of low energy photons. κP is the Planck opacity (see Sect. 2.3). The thermal energy ϵ is given by ϵ = ρcvT (cv being the specific heat at constant volume which we hold fixed and uniform), and u = (ur,uθ,uϕ) the gas velocity. P is the gas pressure, Q+ the viscous dissipation function and S the contribution from the stellar heating.

The radiative flux, using FLD, can be written as (2)where κR is the Rosseland mean opacity and λ the flux limiter. Please note that κP and κR are different opacities. More details concerning the opacities used here are given in Sect. 2.3. Using FLD allows us to solve for stable accretion disc models that cover several vertical pressure scale heights. We use here the FLD approach described in Levermore & Pomraning (1981) with the flux-limiter of Kley (1989) given by (3)where (4)These are the same specifications as in Kley et al. (2009), who used a one-temperature approach with Er = arT4, where ar = 4σ/c is the radiation constant.

For our purposes the physical extent of the star is small with respect to the disc extension so we approximate it as a point source. In spherical coordinates, the stellar radiation is thus propagated along the radial direction only. With this method, the code simultaneously treats the stellar heating component via ray-tracing and subsequent re-radiation via FLD. The stellar heating density S (i.e. the energy deposited per unit time and volume), received by a grid cell i of width Δr is given by: (5)where R is the stellar radius and T the stellar surface temperature. The factor eτi expresses how much stellar irradiation has been absorbed in the grid before it hits the radial grid cell i, with τi being the optical depth integrated radially towards grid cell i. The factor 1−eρκΔr describes how much stellar irradiation is absorbed in the actual grid cell, with κ being the optical opacity, calculated using the stellar spectrum (Dobbs-Dixon et al. 2010). For more information regarding the opacities, please see Sect. 2.3. In the optically thin limit (ρκΔr < 1.0), we make the following approximation: (6)More details concerning the stellar heating function are given in Appendix A.

The coupled energy equations (Eq. (1)) have to be solved simultaneously and we follow here the approach outlined in Commercon et al. (2011). The advection (u·∇)ϵ and the compressional heating (− P∇·u) terms are solved in separate routines. Appendix B describes our approach to simultaneously solving the coupled energy equations (Eq. (1)).

2.2. Physical setup

2.2.1. Computational parameters

For computational reasons, the inner disc (0.0416 AU  < r < 1.04 AU) and outer disc (1.04 AU  < r < 50.96 AU) are treated in two different sets of simulations. Both sets of simulations are 2D discs in the r-θ direction (radial and vertical) with 384 × 32 (inner disc) or 384 × 64 (outer disc) active grid cells. The opening angle used by the grid, i.e. the range of θ, differs for the inner and outer disc simulations. In the inner disc 83° < θ < 90°, while in the outer disc 70° < θ < 90°, where θ is the colatitude measured from the z-axis and θ = 90° is the mid-plane.

The reason for introducing two sets of simulations for the inner and outer disc is numerical. The time step for the inner disc, due to the small inner radius, is very small, significantly increasing the computation time. It is therefore not feasible to simulate the whole disc in one simulation. Additionally, the opening angle of the numerical disc needs to increase to larger values at larger radial distances as the initial disc structure is flared (H/r increasing with r). Setting an opening angle too large (e.g. 20°) for the inner disc results in a collapse of the time step due to very low densities in the upper regions of the disc. Test runs using a density floor showed the same behaviour. Therefore, the two sets of simulations utilize different opening angles. The simulations then can be attached to each other at the respective boundaries of the numerical grids, resulting in a continuous profile. The transfer of the stellar heating is described more precisely in Sect. 2.2.2.

We assume the upper and lower parts of the disc are identical and therefore apply symmetric boundary conditions at the disc’s mid-plane, i.e. at θmax = 90°. At the top of the disc (at θmin) we set the radiation energy to ER = aRT4 with T = 3.0 K (the temperature of the interstellar medium). In this way the disc will always be cooled by the upper boundary and all the heating created in the disc (viscous and stellar) can be transported outwards, as the boundary radiation energy is always lower than in mid-plane. This type of cooling at the top of the disc has been used in Kley et al. (2009), and is also adopted in SPH simulations of fully radiative discs (Ayliffe & Bate 2011). In the radial direction we use reflecting boundary conditions for all disc parameters as described in Kley et al. (2009).

We start the simulations with an initial surface density profile of Σ = Σ0(r/1 AU)−1/2, where Σ0 = 1000 g/cm2 or Σ0 = 3000 g/cm2, which is approximately the value used in Weidenschilling (1977) at 1 AU. The initial aspect ratio H/r ∝ r2/7 indicates a flared disc, following the “flaring disc principle” (Dullemond 2002): if the disc can flare, it will. Self-shadowed discs are discs that cannot be made to flare, which may occur for discs that are initialized with constant H/r.

We use either a constant kinematic viscosity of ν = 1015 cm2 s-1 or an α prescription with different values for α, following Shakura & Sunyaev (1973) with , where cs is the mid-plane sound speed and ΩK the Keplerian orbital frequency. In case of α-viscosity we adopt a vertically constant viscosity utilizing the mid-plane values. MHD simulations indicate that the turbulent stresses do not strongly change with height (Flaig et al. 2012). For discs with constant viscosity, the initial surface density profile is the equilibrium profile, which cancels viscous transport. For α-discs the profile will evolve until a new equilibrium profile is achieved under the influence of reflecting boundary conditions. In this way, the mass inside the disc is conserved. This case is what we call equilibrium disc.

2.2.2. Stellar heating

The star at the centre of our grid has M = M, T = 5656 K and R = 3.0   R, corresponding to a typical protostar. In our simulations the inner discs starts at r = 0.0416 AU, close to the corotation radius of the star.

Let us first assume that the inner edge of the disc is sharply cut off and it cools through that edge as a blackbody. If all of the stellar irradiation is absorbed at this edge, the received energy per time is (7)where 2πrmin2H is the surface of the inner face of the disc. rmin denotes the radius of the inner edge of the disc. The blackbody cooling of the surface of the inner edge is given by (8)with TD being the disc’s temperature. In equilibrium (heating is equal to cooling) we get (9)Using vertical hydrostatic equilibrium, (10)with μ being the mean molecular weight and ℛ the gas constant, we obtain (11)Equation (11) indicates that the aspect ratio of the inner rim scales with the grid truncation radius (rmin) as ; decreasing rmin results in a decreased H/r. Therefore, regardless of our choice of rmin, provided that the chosen value of rmin is smaller than the radius of the innermost bump in H/r due to opacity transitions (see Sect. 4), we are guaranteed that the disc interior to this point will not have a larger aspect ratio than that at the chosen truncation radius. This has important consequences for the irradiation received by this first cell, as it guarantees this first cell will not be in the shadow of a higher H/r region interior to it. Though the simulation naturally captures self-shadowing in the simulated regions of the disc, Eq. (11) shows this is not relevant for the region interior to rmin. Therefore, if rmin is larger than the physical truncation radius of the disc, we can mimic the absorption of stellar irradiation due to the inner disc by simply reducing the stellar irradiation by a factor of eτinner before it hits the first active grid cell. In our case τinner corresponds to the optical depth of the inner two ghost cells of the numerical grid, that are located inside of rmin. These ghost cells have the same properties (e.g. ρ, κ, Δr) as the first active cell. By choosing the values of the ghost cells for eτinner, we follow the radial profile of the disc, making τinner consistent with the rest of the simulation.

As some of the stellar irradiation is absorbed by the inner disc, the outer disc will not receive stellar irradiation in the optically thick mid-plane region. From the simulations of the inner disc, we can measure to what angle from mid-plane in the disc the stellar irradiation is absorbed, ϑabs, by looking at the distribution of stellar heating in the inner disc. This angle ϑabs varies with the disc mass and viscosity. In the outer disc simulation, the disc will only receive stellar heating for an opening angle ϑo larger than ϑabs (keep in mind that ϑ = θ − 90° = 0 represent the mid-plane of the disc): (12)For the simulations of the outer disc, we reduce the stellar irradiation on the top layers by eτinner as we did for the simulations of the inner disc. This has the effect that the transition from zero stellar irradiation to stellar irradiation is smoothed out a little bit. This, of course, implies that the simulations of the inner disc have to be done before the simulations of the outer disc to obtain the correct ϑabs.

In principle one could radially integrate the stellar flux from the inner boundary of the inner simulation to the outer boundary to see what remains. This would then act as the starting stellar heating of the simulations of the outer grid. However, as the simulations of inner and outer disc have different opening angles, this could only be done up to the opening angle of the inner disc. Also, as the inner parts are (radially) optically thick, most of the stellar irradiation is absorbed in the first cells, resulting in a zero stellar irradiation at the outer boundary of the inner disc simulations. Only the very top layers of the inner disc simulations still receive stellar irradiation at their outer boundary, which can be mimicked by the introduction of ϑabs. In Appendix A.3 we present test simulations showing that the inner and outer disc simulations match perfectly with our simple prescription of ϑabs.

2.3. Opacities

In the previous sections we have introduced three different opacities, the Planck mean opacity κP, the optical opacity κ and the Rosseland mean opacity κR. In principle the opacities can be calculated as in Dobbs-Dixon et al. (2010). The Rosseland and Planck mean opacities are defined in the usual manner, so the local Planck mean opacity for the low-energy photon group is given by (13)while the optical opacity is given by the high-energy photon group (14)The frequency dependent opacity is given by κν and the subscript ns indicates scattering processes, which are neglected when calculating wavelength dependent opacities. In principle, the spectrum of impinging radiation can differ from a blackbody, but for our purposes we set Jν(T) = Bν(T). Finally, the Rosseland mean for the low-energy group is given by (15)The wavelength-dependent opacities include the effect of scattering (subscript s). Dobbs-Dixon et al. (2010) state that the ratio between κ and κP is about a factor of 10. They calculated the opacities directly from wavelength-dependent opacities using gaseous, solar-composition opacities. The net result of their work is that the photospheres for the stellar irradiation and the cooling are physically disconnected. The upper layers of the disc (where stellar irradiation plays the largest role) is largely depleted of dust, particularly at a later stage of evolution, when planets are forming. Thus the opacity due to grains is reduced (as is done in Ayliffe & Bate 2009). In our first set of simulations, we leave the Rosseland and Planck opacities constant at 1.0 cm2/g, but the optical opacity is reduced relative to the other used opacities by a factor of 10, so that κ = 0.1 cm2/g.

For simulations with varying opacities, we follow the opacity law of Bell & Lin (1994). The opacity depends primarily on temperature, but also slightly on density, as can be seen in Fig. 1. The plot shows the Rosseland mean and Planck mean opacity. If divided by a factor of 10, the opacities in Fig. 1 reflect the optical opacity. The opacity profile shows several bumps, which are caused by transitions in the opacity regime. For example, the decrease in opacity around 100 K is associated with the melting of ice grains, that dominate the opacity at cooler temperatures. This opacity law is applied for κP and κR, although we point out here the Rosseland mean opacity and the Planck opacity are not the same in reality. However, for the first approximations they are quite similar. Again, the optical opacity κ is 0.1 of the other opacities.

thumbnail Fig. 1

Opacity from Bell & Lin (1994) on a logarithmic scale. Different colours indicate different densities. The bumps in the profile originate from transitions in the opacity regime. This plot features κR and κP as we assume in this paper that κR = κP.

Open with DEXTER

A change in the opacity, will ultimately lead to a change in the discs hydrodynamical structure. Higher values of the Rosseland mean opacity lead to a hotter and vertically more extended disc structure. Semenov et al. (2003) also provided opacity tables for the Rosseland mean opacity and the Planck mean opacity, however the differences between the two do not seem to be very significant, so we use here the same value for κP and κR. The influence of different opacity laws on the disc structure will be investigated in more detail in a future paper, where we will also include a more realistic prescription for the optical opacity.

3. Constant opacity discs

In this section we focus on the structure of discs with a constant opacity, κP = κR = 1.0 cm2/g and κ = 0.1 cm2/g. We split the simulations into an inner disc and an outer disc so we do not have the same time step limitation due to the inner edge of the computational grid. We present in this section viscous simulations with both constant viscosity of ν = 1015 cm2 s-1 and with α viscosity with α = 0.001,0.004,0.008, where , as well as simulations of non-viscous discs. All simulations include stellar heating and the discs all feature Σ0 = 1000   g/cm2. Σ(r) is constant for the simulations with constant viscosity, while Σ(r) changes shape in the discs with α viscosity until an equilibrium is reached.

3.1. Inner disc

3.1.1. Constant-viscosity discs

We expect the structure of the inner disc to be dominated by viscosity. The viscous heating is stronger in the inner parts of the disc compared to the stellar heating because Q+ ∝ r-3 and S ∝ r-2 (see Eq. (16)). In Fig. 2 we present the aspect ratio H/r of discs with constant viscosity. H is defined as H = csK, where cs is the isothermal sound speed computed in the mid-plane of the disc. The discs are evolved from the initially flared structure until they reach thermal equilibrium.

thumbnail Fig. 2

Aspect ratio H/r of the stellar irradiated inner disc calculated using the mid-plane values of cs for discs with constant opacity, κ = 1 cm2/g and constant kinematic viscosity. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium.

Open with DEXTER

The non-viscous disc (blue line in Fig. 2) features a puffed-up inner rim due to stellar irradiation, while the disc behind is shadowed by the rim leading to a vertical contraction and a drop of H/r in the outer parts of the simulations. In theory, if the puffed up inner rim is high enough, the whole outer disc can be shielded from stellar irradiation (Dullemond 2002; Dullemond & Dominik 2004). However, we find that beyond 0.15 AU the discs aspect ratio increases with r, even if the disc remains in the shadow of the inner rim. This is because the heat propagates outwards from the inner rim. Therefore the temperature drops with distance as T = T0rβ. As the disc is in hydrostatic equilibrium, the aspect ratio increases if β < 1.0.

For the viscous disc (red line in Fig. 2) only a very small puffed up inner rim is visible, with a marginally increased H/r near the inner edge of the disc. The otherwise constant aspect ratio is typical of a viscously passive disc, which we show below analytically.

For a purely passive disc (without stellar irradiation), consider the viscous heating Q+ and the radiative cooling Q (blackbody cooling) for an annulus with size A = 2πrδr(16)Please be aware that the disc can cool from both sides, hence the factor 4 in the cooling term. The effective optical depth is given by τeff = 0.5   κΣ. By using the H/r profile for discs in hydrostatic equilibrium

and with τeff = 0.5   Σκ we can equate heating and cooling to find the aspect ratio (17)where Σ = Σ0(r/1   AU)s. Here we have assumed both constant opacity (κ = 1 cm2/g) and constant kinematic viscosity. Therefore, for discs with Σ = Σ0(r/1   AU)-0.5 the result is a disc with a constant aspect ratio, which matches quite well with our numerical simulations.

In total the aspect ratio of the viscous disc is higher compared to the non-viscous disc. This implies that the innermost bump of the disc in the non-viscous disc (which is due to stellar irradiation) does not influence the structure of the viscous disc, as the viscous heating dominates.

3.1.2. α-viscosity discs

For α-viscosity discs we also expect the viscosity to dominate over stellar irradiation in the inner disc. The corresponding H/r profiles are shown in Fig. 3. The black (dotted) lines in the plot show the analytical expectation that neglects stellar irradiation. In fact the aspect ratio of the different α discs vary with the value of viscosity. As the viscosity changes with r, the gradient of the surface density changes as well, influencing the aspect ratio.

thumbnail Fig. 3

Aspect ratio H/r of the stellar irradiated inner disc calculated using the mid-plane values for discs with α-viscosity and constant opacities. The black lines indicate the analytical expectations for the aspect ratio (see text).

Open with DEXTER

For a purely passive disc with α-viscosity, the aspect ratio can be computed in the same way as for a constant viscosity disc. We equate viscous heating Q+ and blackbody cooling Q: (18)where . With the relation of the hydrostatic equilibrium, the heating of the disc is given by (19)By equating heating and cooling, the aspect ratio of the disc can be determined: (20)which clearly indicates that for a large enough gradient of the surface density s, the aspect ratio can decrease with r. This is exactly what happens for our α = 0.008 simulation for which we measure s = 0.85. As the inner disc is dominated by viscosity, this calculation gives an estimate for the aspect ratio of the discs that is in good agreement with our simulation (black lines in Fig. 3).

3.2. Outer disc

Theoretical calculations of Chiang & Goldreich (1997) have shown that the outer disc is flared with H/r ∝ r2/7. This can easily be shown by equating the stellar heating with the cooling at the top of the disc. In our simulations for constant opacity we find exactly this behaviour, as can be seen in Fig. 4 where we display the aspect ratio of discs with viscosity and without. For the simulations with viscous heating an absorbing angle of ϑabs = 6.5° has been used, while for the non-viscous simulations ϑabs = 5.0° has been used, both taken from the results of the inner disc simulations. The aspect ratio of the outer viscous disc matches perfectly with the inner disc (Fig. 2), while there is some small difference for the non-viscous disc. In the non-viscous case the shadowed region behind the innermost rim (in the inner disc simulation) seems not to be captured perfectly by the continuing outer disc simulation. The non-viscous disc simulations show a slight mismatch in the H/r profile between the inner disc and the outer disc. This mismatch cannot be seen in the viscous disc simulations, where H/r matches very well. In reality, accretion discs are viscous, so we are confident that our code reproduces the transition between inner and outer disc very well in that case, which we will focus on in Sect. 4.

thumbnail Fig. 4

Aspect ratio H/r of the stellar irradiated outer disc calculated using mid-plane values for discs with constant opacity κ = 1 cm2/g and constant viscosity. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. These simulations differ from the ones presented in Fig. 2 only by the opening angle and radial extent of the disc.

Open with DEXTER

It can be seen clearly that both discs show a flared profile in the outer regions of the disc, which follows the predicted r2/7 profile. In the outer regions, both disc profiles are nearly the same as stellar heating dominates over viscous heating. In the inner regions of the simulation, small differences between the viscous and non-viscous discs arise, which are due to the importance of viscous heating.

In most 2D simulations of locally isothermal discs (in r − φ plane), a constant H/r profile is assumed, which corresponds to a constant opacity disc without stellar irradiation. However, in radiative discs with temperature dependent opacities the aspect ratio will depend on radius. This is discussed below.

4. Disc structure with temperature dependent opacity

In realistic accretion discs, the opacity is not constant. The opacity depends on the temperature and density (Fig. 1). As the temperature inside an accretion disc drops from a few 1000 K in the inner parts to a few 10 K in the outer parts, we expect the opacity to vary by orders of magnitude as it was already stated in Semenov et al. (2003) and as can also be seen in Fig. 1. As the inner parts of the disc are dominated by viscous heating and as accretion discs are viscously driven, we focus here only on simulations including viscous heating, featuring both constant viscosity and α viscosity. All simulations shown here also include stellar irradiation.

4.1. Inner disc structure

4.1.1. Constant viscosity

In this subsection, all simulated discs have a constant viscosity of ν = 1015 cm2 s-1. In Fig. 5 the aspect ratio and the corresponding opacity κR of the inner disc with Σ0 = 1000 g/cm2 is displayed. Compared to the initial profile the aspect ratio increases, due to viscous heating. This effect is also be seen for a constant opacity (Fig. 2). The drop behind the puffed up inner edge was not visible in the simulations with constant opacity (red line in Fig. 2). As can be seen in Fig. 1, the large features in the opacity profile at high temperature are important for the inner disc.

thumbnail Fig. 5

Aspect ratio H/r of the stellar irradiated inner disc (bottom) calculated using the mid-plane values for Σ0 = 1000 g/cm2 discs with non constant opacity as in Fig. 1 and the corresponding opacity values for κR (top). Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium.

Open with DEXTER

At the inner edge of the grid (r = 0.0416 AU) the temperature drops from  ≈10   000 K to about 3000 K in a few grid cells, which is associated with a large drop in opacity (top in Fig. 5). As the opacity drops, the cooling rate increases, as it is  ∝ 1/κR. An increased cooling also means that the temperature drops even further. As the temperature is linked to H/r, a drop in temperature translates into a drop in aspect ratio. The drop of H/r behind the puffed up inner edge is directly related to a drop in opacity, as we show in Appendix A.2.

If the inner edge were be shifted towards smaller radii, the height of the puffed up inner region would not increase indefinitely, as the opacity does not continue to increase, but rather levels off at high temperatures (see Fig. 1).

As soon as the temperature drops to lower values, the opacity increases again (Fig. 1 and top of Fig. 5), which in turn increases H/r. At around r = 0.83 AU we observe another bump in the H/r profile. Here the temperature crosses the  ≈1000 K bump in the opacity profile, which leads to a decrease of H/r for r > 0.83 AU.

We can clearly relate the transition regions of the opacity to transition regions of the H/r profile of the disc. These changes in the disc’s profile are of crucial importance for irradiated discs. The bump at r = 0.83 AU is higher than the bump at the inner edge, indicating that stellar irradiation will be absorbed at a larger height from mid-plane, which means in turn that the disc will only be heated above a larger angle ϑabs. As the inner region of the disc is dominated by viscous heating, the height and location of the bumps in the disc will change with disc mass and the amount of viscous heating. In contrast, stellar irradiation is relatively unimportant in the inner disc and does not change significantly with disc mass.

4.1.2. α-viscosity

The amount of viscous heating has a crucial influence on the disc structure in the interior of the disc. We now focus on discs with different α viscosities, with , where we use cs in the mid-plane of the disc. The results of simulations of the inner disc with different α and Σ0 values are shown in Fig. 6.

thumbnail Fig. 6

Aspect ratio H/r of the inner disc (with viscous and stellar irradiation) calculated using the mid-plane values for discs with α-viscosity and non constant opacity as in Fig. 1. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium.

Open with DEXTER

As expected the aspect ratio of the disc increases with increasing α parameter due to the increased viscous heating. Also, the aspect ratio increases as we increase Σ0, as more mass inside the disc results in a larger heating of the disc. This has the effect that the peaks in the H/r profile are shifted towards larger radii, as the temperature increases and therefore the transitions in the opacity are reached at larger distances from the star.

The simulation of constant viscosity (with Σ0 = 1000 g/cm2 in Fig. 5) most closely resembles the α = 0.004 simulation. However, the value that matches best might be around α ≈ 0.0055.

With increasing disc mass and viscosity an increase of the height of the peaks in the discs is visible. These bumps in the H/r profile are not directly related to the absorption angle ϑabs of the disc, but it is safe to say that a higher bump in the disc results in a larger absorption angle with important effects on the structure of the outer disc, as discussed below.

4.2. Outer disc

4.2.1. Constant viscosity

We have shown in Sect. 3 that for discs with constant opacity the discs are flared with H/r ∝ r2/7 in agreement with Chiang & Goldreich (1997). However, for non constant opacities we expect another feature in the disc’s structure as the temperature crosses 100 K (Fig. 1). The aspect ratio of the outer disc structure and the corresponding opacity κR for two different surface density values are presented in Fig. 7.

thumbnail Fig. 7

Aspect ratio H/r of the stellar irradiated outer disc (bottom) calculated using the mid-plane values for discs with constant viscosity and non constant opacity as in Fig. 1 and the corresponding opacity values for κR (top). Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. The different colours indicate different Σ0-values. The red line represents the outer disc simulation corresponding to the inner disc simulation shown in Fig. 5.

Open with DEXTER

At r = 1.04 AU the aspect ratio of the Σ0 = 1000 g/cm2 outer disc does not exactly match the aspect ratio computed for the inner disc (red line in Fig. 5). In this region the disc is shadowed from direct stellar irradiation and any radial transfer of energy relies on the re-radiated radiative flux. The outer simulation does not include a source of the re-radiated flux at the inner boundary, leading to the small difference in the inner disc.

The disc with Σ0 = 1000 g/cm2 features another bump in the H/r profile at r ≈ 7.3 AU. At this point in the disc the temperature crosses T ≈ 100 K, which is the transition region for melted ice grains to solid ice grains. The increase of H/r for r < 7.3 AU is related to the increase of opacity for 200 K  > T > 100 K, which can clearly be seen in the top of Fig. 7. The minima and maxima of opacity and H/r are not exactly at the same location, but are slightly shifted. However, a clear trend is still visible. Interestingly, for r > 7.3 AU the H/r profile monotonically decreases and the profile follows that of an non-irradiated disc due to the self-shadowing of the disc.

The fact that the disc is not flared in the outer parts of the disc is a combination of the value of the surface density and of the opacity. In the constant opacity scenario, the outer parts of the disc were flared (Fig. 4) for the Σ0 = 1000 g/cm2 disc. However, in reality for lower temperatures (T < 100 K), the opacity drops to very small values (see top in Fig. 7). Please keep in mind that κ = 0.1κR. As the opacity drops, so does the optical depth τ = ρκΔr of each grid cell, which means that less and less stellar irradiation is absorbed in the outer regions of the disc, as the absorption is proportional to τ (Eq. (6)).

To keep the disc flared, the upper layers in the outer parts of the disc have to absorb stellar irradiation. As the opacity drops in this region of the disc, an increase in density would keep the disc flared, as τ = κρΔr. Indeed, this can be seen in Fig. 7, where the blue line indicates a disc with Σ0 = 3000 g/cm2. The outer parts are flared and follow the 2/7th profile. Also interesting to note is that the bumps in the inner region of the disc in the H/r profile are shifted to larger distances from the star compared to the lower mass disc (red in Fig. 7). This is due to the fact that the higher mass disc produces more viscous heating and is therefore warmer. We compare the disc with Σ0 = 3000 g/cm2 in great detail with a disc that is only heated through viscosity in Sect. 4.3.

4.2.2. α-viscosity

The bumps of the inner disc can shield the outer disc from stellar irradiation implying the outer disc structure is influenced by the inner disc structure, which, in turn, is determined by viscosity and disc mass. In Fig. 8 the aspect ratio for discs with different α-viscosity and disc mass are displayed.

thumbnail Fig. 8

Aspect ratio H/r of the outer disc (with α-viscosity and stellar irradiation) calculated using the mid-plane values for discs with non constant opacity as in Fig. 1. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. These simulations differ from the ones presented in Fig. 6 only by the opening angle and radial extent of the disc.

Open with DEXTER

The Σ0 = 1000 g/cm2 discs are not flared for any viscosity. The outer disc is just too thin to absorb stellar irradiation effectively and viscosity is negligible there, as was seen the constant viscosity disc (Fig. 7). However, we do see the innermost bump in the disc rise as the viscosity increases.

For the Σ0 = 3000 g/cm2 disc with α = 0.001 the equilibrium disc structure features a flared disc profile, which is quite similar to that of a disc with constant opacity (Fig. 7). However, for α > 0.004, the discs are not flared any more. As the viscosity increases the inner bumps of the disc structure and block stellar irradiation from going through to the outer parts of the disc. This effect increases with increasing viscosity.

To have a disc flared in the outer regions, the disc must not only have the right amount of mass, but also the viscosity of the disc should not be too large, as the viscously dominated inner disc can block stellar irradiation from the outer disc. Therefore, simulations of the inner disc are always needed to determine the structure of the outer disc correctly. Whether a disc is shadowed or not, is determined through viscosity and disc mass.

4.3. Passive disc

We now focus on the difference between stellar irradiated discs and discs without stellar irradiation. For this purpose we refer to the case of constant viscosity ν and Σ0 = 3000 g/cm2, as it was already shown in Fig. 7.

In Fig. 9 the two aspect ratios of the disc with Σ0 = 3000 g/cm2 are displayed. In the inner parts, the aspect ratios for both discs are nearly identical, as this region is dominated by viscous heating and not by stellar irradiation. Beyond  ≈20 AU the stellar irradiated disc is flared, while the disc without stellar irradiation collapses to small H/r at large radii. The 100 K bump at r ≈ 9 AU is more prominent in the stellar irradiated disc and also shifted a little bit towards larger distances compared to the non-stellar irradiated disc, as the disc receives more heat which shifts the bump.

thumbnail Fig. 9

Aspect ratio H/r of the outer disc calculated using the mid-plane values for discs with constant viscosity and non constant opacity as in Fig. 1. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. The plot features one disc with and one without stellar heating; both discs have a constant viscosity.

Open with DEXTER

The 2D density map of the discs is shown in Fig. 10. The density profile of the disc without stellar heating does not extend to the same height, because we only use 83° ≤ θ ≤ 90° for this disc. We reduced the vertical extent of the grid in this case because of the smaller disc thickness. In the disc with stellar irradiation, on the other hand, the opening angle of the grid is larger, as we need the disc to absorb the stellar irradiation in the top layers.

thumbnail Fig. 10

Density map (in log-scale) of discs with stellar heating (top) and without stellar heating (bottom). Both discs undergo viscous heating and have Σ0 = 3000 g/cm2. The black dotted line indicates H(r) and the red solid line represents the τl = 1.0 line integrated from the top of the disc.

Open with DEXTER

In Fig. 10 the H(r) line is indicated in black. In both cases the H(r) line rises, however, only for the stellar irradiated disc does H/r increase (Fig. 9). On the other hand, the bumps, which are clearly visible in the H/r profiles are not visible in the H(r) profile, because the increase of H(r) is too strong to notice the small bumps.

Of interest is also the red line in Fig. 10, which refers to the τl = 1.0 line integrated from the top (infinity) of the disc. We define in this case (21)where z is the vertical thickness of the disc and the integration is performed on lines of constant radius r. This line is roughly connected to the flux F (Eq. (2)) of the disc and therefore gives a rough estimate of the location of the photosphere where the disc is cooling. One can clearly see the bumps in the τl = 1.0 line for the stellar irradiated disc. These bumps are at the same location as the ones in the H/r profile, which are not visible in this figure. As τ depends on κR and ρ, it is no surprise to see a drop of density just above the τl = 1.0 line around r ≈ 15.6 AU. The bumps of the τl line can not be seen in the non-irradiated disc. This line just levels off for large distances to the star, where the disc must contract from its initial state to maintain its energy loss to the upper and lower layers.

One of the most important indicators of the disc structure is the temperature of the disc, as the temperature is related to the entropy. Due to the stellar irradiation we expect the upper layers of the disc to be heated from the star, while the mid-plane regions far away from the star are heated vertically from the upper layers. The temperature in mid-plane might therefore be slightly smaller compared to the top layers. The temperature is displayed in Fig. 11. Again we over plot the line for τl = 1.0 (red), H(r) (black) and additionally in blue the τ = 1.0 line, which we define as (22)where the integration is performed along lines of constant θ. The τ line represents where in the disc the stellar heating is absorbed effectively. This line does not match the τl line at 10 AU < r < 25 AU, where τ > τl. This means the stellar heating of the disc is effective higher above the mid-plane compared to the cooling, which results in a higher temperature at the upper layers of the disc, as can be seen clearly in Fig. 11.

thumbnail Fig. 11

Temperature map (log-scale) of a disc with stellar and viscous heating. Indicated in red is the τl = 1.0 line, blue the τ = 1.0 line and in black the H(r) line. The disc has Σ0 = 3000 g/cm2.

Open with DEXTER

The upper layers of the disc are heated by stellar irradiation, while the disc’s interior is also heated through viscosity. However, the heating effect of viscosity is diminished in the upper layers of the disc, as the density is much lower compared to the mid-plane of the disc (see Fig. 10). In the region of the bump in optical depth (10.4 AU  < r < 23.4 AU), the effect of stellar heating is reduced as the opacity is reduced, which leads to less absorption of stellar irradiation. In this region the upper layers of the disc are much hotter compared the regions of the disc right below. This temperature inversion in vertical direction of the disc could have interesting implications to planets on inclined orbits in this part of the disc, as the migration speed is dependent on the gradient of temperature.

In the outer regions of the disc, viscosity is negligible and the disc is heated solely by stellar irradiation. The heat is then transferred from the upper regions to mid-plane, which causes the temperature to show only little vertical dependence.

Observations indicate that protoplanetary discs have an outwards flared profile, which is consistent with our stellar irradiated model. The differences in the disc structure between the two models are dramatic, and we therefore recommend using stellar irradiated disc models for the outer parts of the disc. However, as we have seen in Fig. 7 and Fig. 8, stellar irradiated discs can be self shadowed if the disc mass is small or the viscosity large. In both cases the inner bumps of the discs are so large that they block stellar irradiation from propagating to the outer parts of the disc. For low disc mass models, passive discs follow the same H/r profile as stellar irradiated discs, because the discs are optically thin and can not maintain the flared disc profile.

5. Implications to planet migration

To estimate what is the difference between the torque acting on planets in discs with and without stellar irradiation, we apply the torque formula of Paardekooper et al. (2011) to the discs described in Sect. 4.3. These discs feature Σ0 = 3000 g/cm2, a constant viscosity ν and the opacity law in Fig. 1. The formula captures the behaviour of the torque caused by Lindblad resonances and horseshoe drag on low-mass planets embedded in gaseous discs in the presence of viscous and thermal diffusion. This formula gives the best match to the full 3D radiative simulations of Bitsch & Kley (2011). We do not simulate planets embedded in stellar irradiated discs in this work, as it is beyond the scope of the paper at this point.

The formula of Paardekooper et al. (2011) is very complex and features many details, so that we do not cite every step of the formula at this point. The formula is also additionally displayed in the Appendix of Bitsch & Kley (2011). However, we want to point out a few key items of the formula. The total torque acting on an embedded planet is a composition of its Lindblad torque and its corotation torque: (23)The Lindblad torque depends on the gradients of temperature T ∝ rβ and surface density Σ ∝ rs. It is given in Paardekooper et al. (2011) by (24)where q is the mass ratio between planet and star, h = H/r, ΣP the surface density of the disc at the planets 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, S ∝ rξ, with ξ = β − (γ − 1.0)s. The largest contribution of the corotation torque arises from the entropy related horseshoe drag, which is given by (25)The aspect ratio H/r of the disc with and without stellar irradiation changes from flared to non flared in the outer parts of the disc. This change is related to a change in temperature and the local temperature gradient. In Fig. 11 the outer parts of the disc seem to have a very small radial temperature gradient that reduces the effect on the entropy related horseshoe drag. In the disc without stellar irradiation, on the other hand, planets can migrate outwards to quite large distances, depending on the disc mass (Bitsch & Kley 2011). We therefore expect that the zero-torque radius for planets will be at a smaller distance to the star for stellar irradiated discs compared to discs with only viscous heating, as the temperature gradient is larger in the latter case.

In Fig. 12 the torque acting on planets with different masses in a stellar irradiated disc and a non stellar irradiated disc is displayed. The black lines encircle the regions of outward migration. If a planet is inside the black circles it will migrate outwards, but it migrates inwards for the rest of the positions in the diagram. At the left side of the black circles, a planet would face an unstable torque equilibrium, as the planets would migrate away from the line in both directions, while at the right side of the circles the planet are in a stable equilibrium as they migrate towards the line from both directions.

thumbnail Fig. 12

Torque acting on planets with different masses using the formula of Paardekooper et al. (2011). The top plot features a disc with stellar irradiation and viscous heating (constant ν) and the bottom plot a disc with just viscous heating. The colour coding of the torque has been cut off at Γtot = −0.00015, so that everything that is black in the figure features a large negative torque. The black lines encircle the regions of outward migration.

Open with DEXTER

We only focus here on planets with a few earth masses, as larger planets open gaps in discs and migrate inwards in type-II-migration (Lin & Papaloizou 1986a,b; Crida et al. 2006; Crida & Morbidelli 2007). We only display the torque up to rP = 18.2 AU as the torque acting on planets in the stellar irradiated disc is negative for larger distances to the star. On the other hand, the torque is still positive at that point in the disc with only viscous heating. In fact the torque is still positive up to r ≈ 47 AU. One can also clearly distinguish two different regions of positive torque, indicating areas of outward migration. Finally, the region of outward migration starts at a larger distance to the star for the stellar irradiated disc.

For both disc types the inner region (1.56 < r < 4.7 AU) of positive torque is nearly identical, due to the fact that the disc structure is similar in the inner parts of the disc, as the disc there is dominated by viscous heating and not stellar irradiation. However, as soon as the planetary mass increases above 40MEarth the region of outward migration becomes smaller with increasing planetary mass, and the planets migrate inward for all distances for MP > 60MEarth.

Planets can migrate outwards as long as the libration time-scale is comparable to the cooling time-scale (Baruteau & Masset 2008). For large distances to the central star (Bitsch & Kley 2011) the density is lower in the disc, leading to a shorter cooling time-scale that prevents outward migration.

As the planetary masses increase (MP > 35MEarth) they can start to open up partial gaps inside the disc. The positive corotation torque originates from a region very close to the planet (Kley & Crida 2008; Kley et al. 2009), which is then depleted as the gap starts to form, thus preventing outward migration.

The regions of outward migration in the disc relate to the H/r-profile. As can be seen in Fig. 9, whenever H/r decreases in the disc, the planet migrates outwards. When H/r increases the planet migrates inwards. This means, whenever we have a change in the H/r profile, we change the direction of migration. As the bumps and dips in the H/r profile are influenced by the opacity transitions, the opacity determines the migration. As the transitions of opacity are shifted away from the star with higher viscosity and larger surface density, so is the zero-torque radii in the disc.

To clarify the statement above, keep in mind that the torque depends on the surface density, temperature and entropy gradient, Γ = f(s,β,ξ). (26)which indicates if the disc shows a drop off in H/r, the temperature gradient increases, which therefore increases the entropy gradient as ξ = β − (γ − 1.0)s. Therefore the contribution of the entropy related corotation torque is stronger if H/r decreases. In the flaring part of the disc, where H/r ∝ r2/7, the entropy gradient is smaller as b = −2/7, which leads to β = 3/7, which reduces the entropy gradient ξ, resulting in a negative torque.

On the other hand, it seems that for planet masses with MP < 7MEarth the torque is always negative and that planets migrate inwards at all radii. The reduction of the positive torque for low mass planets is due to the fact that the horseshoe region narrows, thus the horseshoe drag becoming less pronounced, resulting in a smaller net-torque acting on the planet.

For these small planetary masses the inward migration speed is about the same for the passive and stellar irradiated disc. It also seems that the minimum mass a planet needs to sustain outward migration is increased in the case of a stellar irradiated disc, compared to the non stellar irradiated case. The reason for this is a slight change in the entropy gradient of the disc, which is steeper in the non stellar irradiated disc. However, we believe that even smaller mass planets can migrate outwards as the torque formula shows several issues that should be kept in mind:

  • The formula by Paardekooperet al. (2011) is dependent on thesmoothing of the planetary potential, which influences thetorque. A change in the smoothing changes the torque acting onthe planet (see Appendix in Bitsch & Kley 2011).For Fig. 12 we used a smoothing of 0.4.

  • The formula shows some differences to full 3D simulations (Bitsch & Kley 2011), as it was derived from 2D simulations.

  • 3D radiation-hydrodynamical simulations have shown that small mass planets with 5MEarth can still migrate outwards (Kley et al. 2009).

To verify these assumptions, high resolution 3D simulations with embedded small mass planets must be done, as the minimal mass for outward migration has not been investigated with enough detail. We will address this issue in a further paper.

However, the torque formula gives an estimate of the strength of migration in the disc. The difference in the range of outward migration for the two disc types does not change the implication of the existence of the zero-torque radius. At zero-torque distance, planetary embryos and protoplanets can gather and collide and then form larger objects. This is of crucial importance for the growth of massive planets, where a core of a few earth masses needs time to accrete gas to form a gas giant planet (Pollack et al. 1996).

6. Summary

We have investigated the influence of opacity and stellar irradiation on the structure of protoplanetary accretion discs. Utilizing our calculated disc structures we have estimated the torque on embedded planets in discs with and without stellar irradiation by using the theoretical torque formula from Paardekooper et al. (2011).

Before investigating the influence of opacity on the disc structure, we compared non-viscous and viscous discs in the inner and outer regions of the disc for a constant opacity. In the non-viscous disc, a puffed-up inner rim shields the first part of the disc from stellar irradiation, while the outer part of the disc is flared with H/r ∝ r2/7, as was shown in theoretical calculations by Chiang & Goldreich (1997). In the viscous case, in the inner part of the disc, viscous heating dominates over stellar irradiation. By equating viscous heating with radiative cooling, we have shown that the aspect ratio H/r is constant until stellar heating dominates. This was confirmed by our simulations.

We follow the opacity law of Bell & Lin (1994) for our non constant opacity discs. As the opacity changes due to transitions in the opacity law, so does the disc structure. We have seen that the bumps and dips in the opacity law reflect the bums and dips in the discs structure for stellar irradiated and non irradiated discs. If the disc is more massive, the disc produces more viscous heating, which results in a higher temperature at larger distances. But as the temperature determines the opacity, the bumps in the disc are moved outwards with increasing temperature.

The outer parts of the disc (r > 7.8 AU) can be flared if the disc mass is high enough. The effect of the disc mass was not visible in the constant opacity simulations, as a constant opacity leads to a larger optical depth in the upper, less dense regions of the disc. The optical depth is crucial in determining how much stellar irradiation is absorbed. If the disc is optically thin, the stellar irradiation just passes through the disc without any heating effects. Therefore the disc absorbs less heat and can not sustain the flared profile. A minimum mass inside the disc is needed to keep the discs flared.

For increasing viscosity, the aspect ratios of the discs increase. This also means that the bumps due to the opacity in the disc become higher. A higher bump in the disc leads to more absorbed stellar irradiation in the region of the bump. If the bump is high enough, it prevents stellar irradiation from reaching the outer disc, resulting in a non-flared disc. Very high viscosity therefore leads to self-shadowed discs (Fig. 8).

In the inner parts of the disc, where viscous heating dominates, the disc structures for discs with and without stellar irradiation are similar. In the outer regions, the stellar irradiated disc can be flared, which is not the case for the non-stellar irradiated disc. In the outer parts of the stellar irradiated disc, where the viscous heating is negligible, the disc shows a small dependence in T(z), as the heat is transported from the upper layers towards mid-plane.

To estimate the torque acting on planets in discs with stellar irradiation, we apply the torque formula by Paardekooper et al. (2011), where we expect the different disc structure in the outer parts to result in different migration scenarios. The discs feature two regions of outward migration. One in the inner part of the disc (1.56 < r < 4.7 AU) that is nearly identical for the discs with and without stellar irradiation, as viscous heating dominates in that region of the disc. The second region of outward migration is completely different for the two discs. It is much smaller for the case of stellar irradiation, as the very outer parts of the disc have nearly a constant temperature, which favours inward migration. In the case of only viscous heating, outward migration continues to much larger distances from the star (r ≈ 47 AU).

However, important here is that the disc features regions of outward migration for different planetary masses. At the corners of these regions (zero-torque radius), planetary embryos can accumulate and grow until the cores are large enough to accrete gas and form gas giant planets. We note that the results from this torque formula are only an estimate and real 3D simulations with embedded planets need to be done, especially for low mass planets, as they all would migrate towards the star according to the torque formula.

As the differences between stellar irradiated discs and passive discs are dramatic in the outer regions of the disc, we recommend using stellar irradiated disc models for high mass discs. Low mass discs are self-shadowed, so that the structure is well captured by passive disc models as well.

Acknowledgments

B.B. and A.M. have been sponsored through the Helmholtz Alliance Planetary Evolution and Life. W.K. has received support through the German Research Foundation (DFG) through grant KL 650/11 within the Collaborative Research Group FOR 759: The formation of Planets: The Critical First Growth Phase, and gratefully acknowledges the very kind hospitality of the Observatoire de Côte d’Azur. I.D.-D. is supported by the Carl Sagan Postdoctoral program. The calculations were performed on systems of the Computer centre of the University of Tübingen (ZDV) and systems operated by the ZDV on behalf of bwGRiD, the grid of the Baden Württemberg state. Finally, we gratefully acknowledge the helpful and constructive comments of an anonymous referee.

References

Appendix A: Stellar heating

A.1. Calculation of stellar heating in the code

To compute the amount of stellar heating that is deposited in a grid cell, we compute the difference between what arrives and what leaves a grid cell: where i + 1 marks the i + 1th grid cell. The total flux emitted by the star is given by . The stellar flux per surface area on a sphere with radius r is thus (A.3)With the front area of a grid cell given by (A.4)we can compute the flux received by a single grid cell (A.5)which then leads to the stellar heating s of a grid cell As the energy equation (Eq. (1)) is written for energy densities, we have to divide s by the volume of the grid cell V = r2ΔrΔϕ(cosθ1 − cosθ2) to get the same units and ultimately to get S = s/V as it is used in the code: (A.8)This result yields a dependence on resolution. However, this is only important for the optically thick regions (ρκΔr > 1); for optically thin regions we make the approximation (1−eρκΔr)/(Δr) = ρκ. In the inner, optically thick regions of the disc the resolution influences the repartition of the stellar heating, but it does not influence the global net heating. Resolution studies indicate that the height of the puffed up inner edge varies by  ≈3% when lowering the resolution by 33%. As the disc is radially optically thick, the stellar flux is absorbed in the first cells; beyond that the simulations match perfectly for all tested resolutions.

A.2. Innermost bumps and relation to opacity

The innermost bump shown in Fig. 5 and Fig. 6 can be mistaken as a puffed up inner rim caused by stellar irradiation, which would be in contradiction to what is stated in Eq. (11). However, the temperature in that region is a few 1000 K. At that temperature the opacity profile shows a kink (see Fig. 1), indicating that the disc structure will change. In Fig. 5 the change of H/r can be related to a change in opacity. As the opacity drops for slightly larger distances than rmin, so too does H/r. As the opacity determines the disc structure, we show in Fig. A.1 simulations of the inner disc with different rmin for the same disc configuration as in Fig. 5.

thumbnail Fig. A.1

Inner bump of the disc featuring different rmin for the same disc configuration as in Fig. 5. The top panel shows κ, while H/r is shown in the bottom panel.

Open with DEXTER

For distances larger than 0.0416 AU, the opacity drops (top in Fig. A.1) and therefore H/r drops as well (bottom in Fig. A.1). For distances closer to the star, the opacity drops (as already indicated by Fig. 1, because of the higher temperature) and H/r decreases again. This clearly shows that the innermost bump is related to a transition in opacity, as the aspect ration of the disc drops for smaller and larger distances from 0.0416 AU.

As the opacity for even larger temperatures decreases continuously (due to the assumed power-law in the opacity for these temperatures) our simulations cover the innermost rim, indicating that no additional rim is inside.

A.3. Matching of inner and outer disc

The models of the inner and outer disc can be attached to each other. To do so, the regions of the simulations have to overlap. In Fig. A.2 we present the overlapping region of the α = 0.008 simulation with Σ0 = 3000 g/cm2, which are actually the simulations presented in Fig. 6 and Fig. 8.

thumbnail Fig. A.2

Matching of the inner disc simulations with the outer disc simulations for an α = 0.008 disc with Σ0 = 3000 g/cm2.

Open with DEXTER

Clearly H/r of the simulations of inner and outer disc match very well. We are therefore confident that the approximation of θabs, as presented in Eq. (12), for the stellar heating is sufficient. If the matching was done in a region of the disc that is dominated by stellar irradiation, the matching might fall into a shadowed region of the disc that would then not be captured by the outer disc leading to a mismatch of the two simulations. However, if, as in our simulations the matching occurs in a viscously dominated region of the disc, the outer discs inner edge will adapt to the outer edge of the inner edge due to the viscous heating.

Appendix B: Energy equation

For simplicity the derivation of solving the energy equation (Eq. (1)) is shown in Cartesian coordinates in the implemented 3D algorithm, although we used 2D simulations in the r-θ plane in this work. The radiation energy equation is then given by (B.1)where the upper indices (n and n + 1) indicate the time step and are the averaged diffusion coefficients, which are given by (B.2)The term (Tn + 1)4 is non linear and makes the scheme difficult to invert, yet it is much easier to solve a linear system implicitly. Assuming that the changes in temperature are small in each time step (Commercon et al. 2011), we can write (B.3)Now, we need to have Tn + 1 as a function of Tn, and to solve Eq. (B.1). We use now the energy density equation, where we omit advection and compressional heating, (B.4)With ϵ = ρcvT we get (B.5)With our approximation for the temperature (Eq. (B.3)) we find (B.6)where We can now plug this all into the radiation energy equation and solve for . We arrive at a matrix equation: where the superscript n + 1 has been omitted on the left hand side. The matrix elements are given by: Written in matrix notation we have (B.9)which is a matrix that is very similar in nature to the one for the one-temperature energy equation and can be solved with the same matrix solver as in Kley et al. (2009). After the matrix iteration, the new temperature of the system can be calculated with Eq. (B.6), followed by the computation of ϵ = ρcvT.

All Figures

thumbnail Fig. 1

Opacity from Bell & Lin (1994) on a logarithmic scale. Different colours indicate different densities. The bumps in the profile originate from transitions in the opacity regime. This plot features κR and κP as we assume in this paper that κR = κP.

Open with DEXTER
In the text
thumbnail Fig. 2

Aspect ratio H/r of the stellar irradiated inner disc calculated using the mid-plane values of cs for discs with constant opacity, κ = 1 cm2/g and constant kinematic viscosity. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium.

Open with DEXTER
In the text
thumbnail Fig. 3

Aspect ratio H/r of the stellar irradiated inner disc calculated using the mid-plane values for discs with α-viscosity and constant opacities. The black lines indicate the analytical expectations for the aspect ratio (see text).

Open with DEXTER
In the text
thumbnail Fig. 4

Aspect ratio H/r of the stellar irradiated outer disc calculated using mid-plane values for discs with constant opacity κ = 1 cm2/g and constant viscosity. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. These simulations differ from the ones presented in Fig. 2 only by the opening angle and radial extent of the disc.

Open with DEXTER
In the text
thumbnail Fig. 5

Aspect ratio H/r of the stellar irradiated inner disc (bottom) calculated using the mid-plane values for Σ0 = 1000 g/cm2 discs with non constant opacity as in Fig. 1 and the corresponding opacity values for κR (top). Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium.

Open with DEXTER
In the text
thumbnail Fig. 6

Aspect ratio H/r of the inner disc (with viscous and stellar irradiation) calculated using the mid-plane values for discs with α-viscosity and non constant opacity as in Fig. 1. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium.

Open with DEXTER
In the text
thumbnail Fig. 7

Aspect ratio H/r of the stellar irradiated outer disc (bottom) calculated using the mid-plane values for discs with constant viscosity and non constant opacity as in Fig. 1 and the corresponding opacity values for κR (top). Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. The different colours indicate different Σ0-values. The red line represents the outer disc simulation corresponding to the inner disc simulation shown in Fig. 5.

Open with DEXTER
In the text
thumbnail Fig. 8

Aspect ratio H/r of the outer disc (with α-viscosity and stellar irradiation) calculated using the mid-plane values for discs with non constant opacity as in Fig. 1. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. These simulations differ from the ones presented in Fig. 6 only by the opening angle and radial extent of the disc.

Open with DEXTER
In the text
thumbnail Fig. 9

Aspect ratio H/r of the outer disc calculated using the mid-plane values for discs with constant viscosity and non constant opacity as in Fig. 1. Starting from a flared disc profile the discs are evolved until they reach thermal equilibrium. The plot features one disc with and one without stellar heating; both discs have a constant viscosity.

Open with DEXTER
In the text
thumbnail Fig. 10

Density map (in log-scale) of discs with stellar heating (top) and without stellar heating (bottom). Both discs undergo viscous heating and have Σ0 = 3000 g/cm2. The black dotted line indicates H(r) and the red solid line represents the τl = 1.0 line integrated from the top of the disc.

Open with DEXTER
In the text
thumbnail Fig. 11

Temperature map (log-scale) of a disc with stellar and viscous heating. Indicated in red is the τl = 1.0 line, blue the τ = 1.0 line and in black the H(r) line. The disc has Σ0 = 3000 g/cm2.

Open with DEXTER
In the text
thumbnail Fig. 12

Torque acting on planets with different masses using the formula of Paardekooper et al. (2011). The top plot features a disc with stellar irradiation and viscous heating (constant ν) and the bottom plot a disc with just viscous heating. The colour coding of the torque has been cut off at Γtot = −0.00015, so that everything that is black in the figure features a large negative torque. The black lines encircle the regions of outward migration.

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

Inner bump of the disc featuring different rmin for the same disc configuration as in Fig. 5. The top panel shows κ, while H/r is shown in the bottom panel.

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

Matching of the inner disc simulations with the outer disc simulations for an α = 0.008 disc with Σ0 = 3000 g/cm2.

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.