© ESO, 2014

1. Introduction

The origin of the angular momentum transport in accretion discs is still not fully understood. Observationally it has been confirmed that the molecular viscosity is by many orders of magnitude too small to explain the effective mass and angular momentum transport in discs (Pringle 1981). This can be inferred for example from time variations in the disc luminosity in close binary systems, or by correlating the infrared-excess caused by discs around young stars with the age of the system. As a consequence it is assumed that discs are driven by some kind of turbulent transport whose cause is still not known. Despite its unknown origin, the efficiency of the turbulence is usually parameterised in terms of the dimensionless parameter, α, as introduced by Shakura & Sunyaev (1973). Observationally, values of a few times 10^-3 as in protostellar discs to 10^-1 for discs in close binary stars are suggested. For sufficiently well ionised discs the magnetorotational instability (MRI) is certainly the most promising candidate to provide the transport (Balbus 2003). While this may be true for the hot discs in close binary systems or in active galactic nuclei, there is the important class of protostellar discs where at least the thermal ionisation levels are too low to provide a sufficient number of charged particles that can support the MRI (Armitage 2011). In such discs turbulence plays an important role in several aspects. Not only does it determine the lifetime of an accretion disc, but it also influences where and how planets can form and evolve in the disc. A variety of sources such as stellar X-rays, cosmic rays or collisions with beta particles from radioactive nuclei have been invoked to provide the required ionization levels, but recent studies indicate the presence of an extended dead zone where, because of the lack of ionization, no magnetically driven instability may operate. Additionally, recent studies on the origins of turbulence in protostellar discs that include non-ideal magnetohydrodynamical (MHD) effects such as ambipolar diffusion or the Hall effect, indicate that the MRI may even be suppressed strongly in these discs, see the review by Turner et al. (2014) and references therein.

As a consequence, alternative mechanisms that provide turbulence are actively discussed. Typical examples for non-magnetised discs are convective instability (Ruden et al. 1988), gravitational instability (Lin & Pringle 1987), or baroclinic instability (Klahr & Bodenheimer 2003), for further references see Nelson et al. (2013). While any of these may operate under special conditions in the disc, e.g. suitable radial entropy gradients or a sufficiently high disc mass, none seems to have general applicability. Searching for alternatives two linear instabilities have been recently discussed in the literature, both acting on the radial temperature structure of the disc. One is a convective overstability that preferably acts for thermal relaxation times close to the orbital period (Klahr & Hubbard 2014; Lyra 2014), and two is a vertical shear instability (VSI) that operates best for cooling times much shorter than the orbital period or for discs that are adiabatically stratified in the vertical direction. The present paper focuses on the latter instability. Here, the instability is caused by a vertical gradient of the angular velocity, Ω, in the disc. Through linear analysis it has been shown that for a sufficiently strong vertical shear there are always modes that can overcome the stabilizing angular momentum gradient (Rayleigh-criterion) and generate instability (Urpin & Brandenburg 1998; Urpin 2003). This instability is related to the Goldreich-Schubert-Fricke instability that can occur in differentially rotating stars (Goldreich & Schubert 1967; Fricke 1968).

Concerning vertical shears effectiveness with respect to angular momentum transport numerical simulations were performed by Arlt & Urpin (2004) and Nelson et al. (2013). The first authors analysed the instability for globally isothermal discs and found that the instability in this case could only be triggered by applying finite initial perturbation because the equilibrium state of the disc (being strictly isothermal) did not contain a shear in Ω. The maximum values of α obtained by Arlt & Urpin (2004) were around 6 × 10^-6, but the turbulence was decaying in the long run. Nelson et al. (2013) extended these simulations and performed high resolution simulations of the VSI for so-called locally isothermal discs that contain a radial temperature gradient, but are vertically isothermal. Under these conditions the equilibrium state has a vertical gradient in the shear and indeed an instability sets in. As shown by Nelson et al. (2013) the instability has two distinct growth phases, it starts from the surface layers of the disc where the shear is strongest and then protrudes towards the midplane. In the final state the vertical motions in the disc are antisymmetric with respect to the disc’s midplane, such that the gas elements cross the midplane, a feature found for the vertical convective motions in discs as well (Kley et al. 1993). For the efficiency of the VSI induced turbulence Nelson et al. (2013) found a weak angular momentum transport with α = 6 × 10^-4. They also showed that in the presence of a small viscosity or thermal relaxation the instability is weaker and can easily be quenched.

It is not clear what influence radiation transport will have on this instability. Without external heat sources one might expect that, because of radiative cooling and the dependence of the instability on temperature, the instability will die out. Here, we evaluate the evolution of the instability for radiative discs and an ideal equation of state. Additionally, we extend the radial domain and include irradiation from the central star. We perform 2D and 3D hydrodynamical simulations including radiative transport.

This paper is organised as follows. In Sect. 2, we present the physical setup of our disc models and in Sect. 3 the numerical approach. The isothermal results are presented in Sect. 4, followed by the radiative cases in Sect. 5. Stellar irradiation is considered in Sect. 6 and in Sect. 7 we conclude.

2. Physical setup

In order to study the VSI of the disc in the presence of radiative transport we construct numerical models solving the hydrodynamical equations for a section of the accretion disc in two and three spatial dimensions.

2.1. Equations

The basis our studies are the Euler Eqs. (1)–(3) describing the motion of an ideal gas. These are coupled to radiation transport (4) for which we use the two temperature approximation applying flux-limited diffusion. The equations then read $\begin{array}{ccc}& & \frac{\mathit{\partial \rho }}{\mathit{\partial t}}\mathrm{+}\mathrm{\nabla }\mathrm{\left(}\mathit{\rho }{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\mathit{\rho }{u}\mathrm{+}\mathrm{\nabla }\mathrm{\left(}\mathit{\rho }{u}{u}\mathrm{\right)}\mathrm{+}\mathrm{\nabla }\mathit{p}\mathrm{=}\mathit{\rho }{{a}}_{\mathrm{ext}}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\mathit{e}\mathrm{+}\mathrm{\nabla }\mathrm{\left[}\mathrm{\right(}\mathit{e}\mathrm{+}\mathit{p}\mathrm{\left)}{u}\mathrm{\right]}\mathrm{=}\mathit{\rho }{u}{{a}}_{\mathrm{ext}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{P}}\mathit{\rho c}\mathrm{(}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{-}\mathit{E}\mathrm{)}\\ & & \frac{\mathit{\partial }}{\mathit{\partial t}}\mathit{E}\mathrm{+}\mathrm{\nabla }{F}\mathrm{=}{\mathit{\kappa }}_{\mathrm{P}}\mathit{\rho c}\mathrm{(}{\mathit{a}}_{\mathrm{R}}{\mathit{T}}^{\mathrm{4}}\mathrm{-}{\mathit{E}}^{\mathrm{)}}\mathit{.}\end{array}$Here ρ is the density; u the velocity; e the total energy density (kinetic and thermal) of the gas; p denotes the gas pressure; the acceleration due to external forces, such as the gravitational force exerted by the central star is given by a_ext; and E and F are the energy density and the flux of the radiation. The last terms on the righ-hand side of Eqs. (3) and (4) refer to the coupling of gas and radiation, i.e. the heating/cooling terms. Here, c stands for the speed of light, a_R is the radiation constant, and κ_P the Planck mean opacity.

We close the equations with the ideal gas equation of state $\mathit{p}\mathrm{=}\mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{e}}_{\mathrm{th}}\mathit{,}$(5)where e_th = e−1 / 2ρu² is the thermal energy density. The temperature of the gas is then calculated from $\mathit{p}\mathrm{=}\mathit{\rho }\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu }{\mathit{m}}_{\mathrm{H}}}\mathit{,}$(6)where μ is the mean molecular weight, k_B the Boltzmann constant, and m_H the mass of the hydrogen atom. In our simulations with radiation transport we use γ = 1.4 and μ = 2.35. To compare to previous studies we performed additional isothermal simulations where we use γ = 1.001 and additionally reset to the original temperature profile in every step. This procedure corresponds to an isothermal simulation, but allows for an arbitrary temperature profile. It also allows to use the feature of slowly relaxing to a given original temperature such as used for example in Nelson et al. (2013). Note that without resetting the temperature the gas remains adiabatic, and the perturbation will die out for our setup.

The radiation flux in the flux-limited diffusion (FLD) approximation (Levermore & Pomraning 1981) is given by ${F}\mathrm{=}\mathrm{-}\mathit{\lambda }\frac{\mathit{c}}{{\mathit{\kappa }}_{\mathrm{R}}\mathit{\rho }}\mathrm{\nabla }\mathit{E,}$(7)where κ_R is the Rosseland mean opacity and λ is the flux-limiter, for which we use the description of Minerbo (1978). For the Rosseland mean opacity we apply the model of Bell & Lin (1994). For simplicity, in this initial study we use the same value for the Planck mean opacity, see also Bitsch et al. (2013).

In some of our studies we add viscosity and stellar irradiation to the momentum and energy equations. This will be pointed out below in the appropriate sections.

2.2. Disc model

To be able to study the onset of the instability we start with a reference model in equilibrium. For this purpose, we follow Nelson et al. (2013) and use a locally isothermal disc in force equilibrium, where for the midplane density we assume a power law behaviour $\mathit{\rho }\mathrm{(}\mathit{R,Z}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{0}}}\right)}^{\mathit{p}}\mathit{,}$(8)and that the temperature is constant on cylinders $\mathit{T}\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{0}}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{0}}}\right)}^{\mathit{q}}\mathrm{·}$(9)To specify the equilibrium state we have used a cylindrical coordinate system (R,Z,φ). However, our simulations will be performed in spherical polar coordinates (r,θ,φ) because they are better adapted to the geometry of an accretion disc. In Eqs. (8) and (9), ρ₀ and T₀ are suitably chosen constants that determine the total mass content in the disc and its temperature. The exponents p and q give the radial steepness of the profiles, and typically we choose p = −3 / 2 and q = −1. Assuming that in the initial state there are no motions in the meridional plane and the flow is purely toroidal, force balance in the radial and vertical directions then leads to the equilibrium density and angular velocity profiles that we use for the initial setup (Nelson et al. 2013) $\mathit{\rho }\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{0}}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{0}}}\right)}^{\mathit{p}}\hspace{0.17em}\exp \left[\frac{\mathit{GM}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\left(\frac{\mathrm{1}}{\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\mathrm{-}\frac{\mathrm{1}}{\mathit{R}}\right)\right]\mathit{,}$(10)and $\mathrm{\Omega }\mathrm{\left(}\mathit{R,Z}\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{K}}{\left[\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}{\left(\frac{\mathit{H}}{\mathit{R}}\right)}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{q}\mathrm{)}\mathrm{-}\frac{\mathit{qR}}{\sqrt{{\mathit{R}}^{\mathrm{2}}\mathrm{+}{\mathit{Z}}^{\mathrm{2}}}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{·}$(11)Here, ${\mathit{c}}_{\mathrm{s}}\mathrm{=}\sqrt{\mathit{p}\mathit{/}\mathit{\rho }}$ denotes the isothermal sound speed, ${\mathrm{\Omega }}_{\mathrm{K}}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}\mathit{/}{\mathit{R}}^{\mathrm{3}}}$ the Keplerian angular velocity, and H = c_s/ Ω_K is the local pressure scale height of the accretion disc. We note that the Z dependence of Ω in the equilibrium state is the origin of the VSI because the vertical shear provides the opportunity for fluid perturbations with a wavenumber ratio k_R/k_Z above a threshold to tap into a negative gradient in the angular momentum as the perturbed fluid elements move away from the rotation axis (Nelson et al. 2013). The angular velocity given by Eq. (11) is also used to calculate the Reynolds stress tensor, for details see below.

2.3. Stability

Nelson et al. (2013) repeated the original analysis in Goldreich & Schubert (1967) for a locally isothermal and compressive gas for an accretion disc using the local shearing sheet approximation at a reference radius r₀. They derived the same stability criterion as Urpin (2003) and obtained the following growth rate of the instability ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{-}{\mathit{\kappa }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{(}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{k}}_{\mathit{Z}}^{\mathrm{2}}\mathrm{+}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{0}}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{k}}_{\mathit{R}}{\mathit{k}}_{\mathit{Z}}\frac{\mathit{\partial }\mathit{V̅}}{\mathit{\partial z}}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{(}{\mathit{k}}_{\mathit{Z}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{R}}^{\mathrm{2}}\mathrm{)}\mathrm{+}{\mathit{\kappa }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}{\mathit{N}}_{\mathrm{0}}^{\mathrm{2}}}\mathit{,}$(12)where κ₀ is the epicyclic frequency, c₀ the sound speed, and N₀ is the Brunt-Vaisaila frequency at the radius r₀; denotes the mean deviation from the Keplerian azimuthal velocity profile, and k_R and k_Z are the radial and vertical wavenumbers of the perturbations in the local coordinates.

For negligible N₀, small H₀/R₀, and , as seen in their numerical simulations, Nelson et al. (2013) find $\mathit{\sigma }\mathrm{~}\mathit{q}\mathrm{\Omega }\frac{\mathit{H}}{\mathit{R}}\mathit{,}$(13)which implies that the growth rate per local orbit to first order depends on the temperature gradient as given by q and on the absolute temperature, because of H/R. We will compare our numerical results with these estimates.

3. Numerical model

To study the VSI in the presence of radiative transport we perform numerical simulations of a section of an accretion disc in two and three spatial dimensions using spherical polar coordinates (r,θ,φ), and a grid which is logarithmic in radial direction, keeping the cells squared. We solve Eqs. (1) to (4) with a grid-based method, where we use the PLUTO code from Mignone et al. (2007) that utilises a second-order Godunov scheme, together with our radiation transport (Kolb et al. 2013) in the FLD approximation, see Eq. (7).

The simulations span a region in radius from r = 2−10 AU, this is the range where the dead zone can be expected (Armitage 2011; Flaig et al. 2012). Here, we use a larger radial domain as Nelson et al. (2013) did because we intend to study the global properties of the instability over a wider range of distances. Additionally, this larger range is useful because we need some additional space (typically ≈1 AU) to damp possible large scale vortices in the meridional plane that show up at the inner radial boundary of the domain (see below). The origin of these vortices is possibly that the instability moves material along cylindrically shaped shells, a motion that is not adapted to the used spherical coordinates, such that the midplane is cut out at the inner boundary. Vortices can also arise if the viscosity changes apruptly, a situation mimicking a boundary. Additionally, in some cases the wavelengths are large, such that the coupling between different modes cannot be captured in a small domain. We also use a wide range because with radiation transport the growth rates are expected to depend on the opacity, which is a function of ρ and T and thus of the radius. In the meridional direction (θ) we go up to ± 5 scale heights above and below the equator in the isothermal case, and we use the same extension for the radiative simulation, where it corresponds to more scale heights. For the 3D simulations we used in the azimuthal direction (φ) a quarter circle, from 0 to π/ 2.

We use reflective boundaries in the radial direction. In the meridional direction we use outflow conditions for the flow out of the domain and reflective conditions otherwise. For the radiation transport solver we set the temperature of the meridional boundary to 10 K, which allows the radiation to escape freely. We use damping of the velocity near the inner radial boundary within 2–3 AU to prevent the creation of strong vortices arising through the interaction with the reflecting boundary, which can destroy the simulation. This is done by adding a small viscosity of ν = 2 × 10^-7 with a linear decrease to zero from 2 AU to 3 AU (similar to the damping used in de Val-Borro et al. 2006).

We assume that the disc orbits a solar mass star and we apply a density of ρ₀ = 10^-10 g/cm³ at 1 AU. Because the surface density decays with r^-0.5, we get a surface density Σ = 80 g/cm² at 5 AU. To study the mass dependence we vary ρ₀ for the radiative models. To seed the instability we add a small perturbation of up to 1% of the sound speed to the equilibrium velocity, see Eq. (11).

Because our radiation transport solver is only implemented in full 3D (Kolb et al. 2013), we use two grid cells in the azimuthal direction for the 2D axisymmetric simulations using radiation transport.

4. Isothermal discs

Fig. 1 Kinetic energy of the motion in the meridional plane at different radii in an inviscid disc. The kinetic energy at the different locations is in each case averaged over a radial interval with length 0.5 AU. We note that the unit of time is given in local periods at the centre of the specified interval. Hence, it is different for each curve, but this allows easy comparison.

Before studying full radiative discs, we first perform isothermal 2D simulations to compare our results and growth rates to those of Nelson et al. (2013). Then we will extend the simulation to full 3D using a quarter of a disc and discuss the dependence on resolution and viscosity.

4.1. Growth rates

To analyse the possible growth and instability of the initial equilibrium state, we analyse the time evolution of the kinetic energy in the meridional plane ${\mathit{e}}_{\mathrm{kin}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}\mathit{\rho }\mathrm{(}{\mathit{u}}_{\mathit{r}}^{\mathrm{2}}\mathrm{+}{{\mathit{u}}_{\mathit{\theta }}^{\mathrm{2}}}^{\mathrm{)}}\mathit{,}$(14)at different radii. The obtained growth of e_kin of a run with q = −1 and p = −3 / 2 is shown at different radii in Fig. 1 for an inviscid disc model with a grid resolution of 2048 × 512. We note that the time is measured in local orbits (2π/ Ω(r_i)) at the corresponding centres of the intervals, r_i. We measure a mean growth rate of 0.38 per orbit for the kinetic energy (light blue line in Fig. 1), which is twice the growth rate (σ) of the velocity. We calculate the growth rate by averaging the kinetic energy at the different r_i over an interval with length 0.5 AU. Our results compare favourably with the growth rates from Nelson et al. (2013) who obtained 0.25 per orbit averaged over 1−2 AU for q = −1. Averaging over this larger range leads to a reduced growth because the rate at 2 AU, measured in orbits at 1 AU, is smaller by a factor of 2^1.5 = 2.8, and so their result is a slight underestimate.

Fig. 2 Velocity in the meridional direction, uθ, in units of local Kepler velocity for an isothermal run without viscosity. The panels refer to snapshots taken at time 100, 210 and 750 (top to bottom), measured in orbital periods at 1 AU. In units of local orbits at (2.5, 3.5, 4.5) AU this refers to (25, 15, 10) (53, 32, 22) (190, 115, 79) orbits, from top to bottom.

A closer look at Fig. 1 reveals two distinct growth phases. An initial strong linear growth phase with a rate of 0.38 per orbit lasting about 20 local orbits, and a slower second phase with a rate of 0.10 per orbit (grey line in Fig. 1). To understand these regimes, we present in Fig. 2 the velocity in the meridional direction, u_θ, in 2D contour plots at different times. The top panel reveals that the first phase corresponds to symmetric (mirror symmetry with respect to the equatorial plane) disturbances that grow from the top and bottom surface layers of the disc. Here, the gas does not cross the midplane of the disc. When the disturbances meet in the disc’s midplane they develop an anti-symmetric phase with lower growth rates where the gas flow crosses the midplane of the disc as shown in the middle panel. The converged phase shown in the lower panel then shows the fully saturated global flow. Figure 2 indicates that in the top panel the whole domain is still in the anti-symmetric growth phase, in the middle panel only the smaller radii show symmetric growth, while in the lower panel the whole domain has reached the final equilibrium, in accordance with Fig. 1.

We point out that the growth rate per local orbit (~σ/ Ω) is independent of radius in good agreement with the relation (13), for constant H/R. We will show later that the growth rate is also independent of resolution.

4.2. Comparison to 3D results and Reynolds stress

In addition to the 2D simulation we ran an equivalent 3D case using a quarter of a disc with a resolution of 512 × 128 × 128 grid cells. We will use this to discuss the validity of the 2D results, in particular the estimates on the turbulent efficiency factor α. In Fig. 3 we compare the growth of the meridional kinetic energy for the 3D and the 2D simulation. After a slower start, the 3D simulation shows very similar growth and reaches the same final saturation level.

Fig. 3 Growth of the kinetic energy for the quarter of a disc and the 2D equivalent. The kinetic energy is averaged from 4 AU to 5.5 AU.

To estimate any possible angular momentum transfer caused by the turbulent motions induced by the instability we calculate the corresponding Reynolds stress (Balbus 2003) ${\mathit{T}}_{\mathit{r\phi }}\mathrm{=}\frac{\mathrm{\int }\mathit{\rho \delta }{\mathit{u}}_{\mathit{r}}\mathit{\delta }{\mathit{u}}_{\mathit{\phi }}\mathrm{d}\mathit{V}}{\mathrm{\Delta }\mathit{V}}\mathrm{=}\hspace{0.17em}\mathrm{⟨}\mathit{\rho \delta }{\mathit{u}}_{\mathit{r}}\mathit{\delta }{\mathit{u}}_{\mathit{\phi }}\mathrm{⟩}\mathit{,}$(15)where δu_r and δu_φ are defined as the fluctuations of the velocity field from the mean flow and ΔV is the volume of the integrated domain. To calculate a coordinate dependent stress we integrate only over thin slices with a thickness of one cell in the apropriate direction. While δu_r is just the radial velocity, u_r, at the point of interest because the initial u_r was zero, δu_φ is difficult to calculate, as one has to subtract the mean background rotational velocity. Armitage (2011) defines it as the difference to the Kepler rotation, while strictly speaking it is the deviation from the unperturbed equilibrium state that is not Keplerian in our case, see Eq. (11). In 3D simulations it is mostly calculated by averaging over the azimuthal direction (Flock et al. 2011; Fromang & Nelson 2006), but this instability is nearly axisymmetric (see Fig. 4), so this is not appropriate here and the correct way is to average over time to obtain the steady-state velocity. However, this is computationally inconvenient because this time average is not known a priori. In Fig. 5 we show that the time averaging method leads to the same results as the equilibrium method using the analytic Eq. (11), and we use the latter for our subsequent simulations.

To calculate the dimensionless α-parameter, T_rφ has to be divided by the pressure. To show the radial and vertical dependence of α it is useful to use different normalisations. We divide the Reynolds stress in Eq. (15) by the midplane pressure to illustrate the dependence on the meridional (vertical) coordinate, thus making it independent of the number of scale heights of the domain. The stress as a function of the radius, T_rφ(R), is divided by the vertical averaged pressure, making it again independent of the numbers of scale heights. This procedure corresponds to a density weighted height integration (Balbus 2003).

Fig. 4 Vertical velocity in the midplane of the disc for the 3D model after 4000 orbits. The nearly axisymmetric property of the instability is clearly visible.

In Fig. 5 we present the different methods for calculating the Reynolds stress, T_rφ, for the simulation of a quarter of a disc with a resolution of 512 × 128 × 128 and the same initial conditions as in the 2D case. We can see that indeed the axisymmetric property of the instability leads to incorrect results if one only averages over the azimuthal direction. All further results for the isothermal discs are calculated with the equilibrium method. This allows us to approximate the Reynolds stress even in a transient disc and calculate the stress continuously during the whole runtime of the simulation, strongly reducing the amount of data needed to be written to the hard drive because the Reynolds stress can now be calculated independent of the other time steps. In addition, the computations show that the stresses of the reduced 2D simulations yield stresses comparable to the full 3D case and can be used as a proxy for the full 3D case. In Fig. 4 we show the vertical velocity in the midplane of the disc for the 3D model. As shown, the motions are only very weakly non-axisymmetric.

Fig. 5 Reynolds stress (code units) from 3−10 AU averaged over 41 time steps, each step 100 orbits apart beginning with orbit 1000, calculated with different averaging methods. For “mean time” the steady-state , needed to calculate the Reynolds stress at each step, was calculated through averaging over the 41 time steps. For “mean phi” the steady-state velocity was calculated by averaging over the azimuthal direction at each time step. For “equilibrium” is calculated analytically by using the equilibrium Eq. (11) at each step. For the 2D model we used the equilibrium method as well.

4.3. Resolution

In this section we look at the effect of resolution. We start with a resolution of 256 × 64, where the instability exists, but clearly is not resolved, and go by doubling the resolution in several steps up to a resolution of 2048 × 512, where the computations start to be expensive. In Fig. 6 we show on the left the Reynolds stress divided by the midplane pressure as a function of vertical distance. This is then averaged over the radius from 3−8 AU. On the right we plot the Reynolds stress divided by the pressure, where both pressure and stress have been averaged over the meridional direction.

From this plot it is not clear if the values for α converge to a specific level for higher resolution. Nevertheless, it gives a first impression on the strength of turbulent viscosity caused by this instability being relatively weak with α-values a few times 10^-4, which is slightly smaller than the value of 6 × 10^-4 found by Nelson et al. (2013).

Fig. 6 Radial and vertical distribution of the Reynolds stress. Left: the Reynolds stress divided by the midplane pressure over the vertical direction. Right: reynolds stress divided by the mean pressure over the radius for different resolutions. Both are averaged over 4001 time steps from orbit 1000 to 5000. The model res2048 corresponds to the results shown in Fig. 1.

Fig. 7 Mean wavenumber of the instability over the radius for different numerical resolutions in the saturated phase. Upper panel: inviscid case with ν = 0, lower panel: viscous case with ν = 5 × 10-7 (dimensionless).

In Fig. 7 we show the wavelength of the perturbation as a function of radius for different numerical resolutions, where the wavelength has been estimated by measuring the distance between two sucsessive changes of the sign of the vertically averaged vertical momentum after the instability is saturated (see Fig. 2, third panel, or Fig. 10 along the radius axis, beginning with orbit 1000). This does not, of course, reveal the full spectrum, but at this point we are more interested in the characteristic mean wavelength. We note that the wavelength in the growth phase can be smaller. In all shown resolutions one wavelength is resolved with 15–50 grid cells, while larger radii are better resolved. Despite the variation with radius one notices in Fig. 7 that the wavelength clearly depends on the numerical resolution. One possible cause for this is the lack of physical viscosity. Because the (intrinsic) numerical viscosity of the code decreases with increasing resolution, this may explain the missing convergence, in particular since the growth rates depend of the wavenumbers of the disturbances, see Eq. (12). We repeated the run with an intermediate resolution of 1440 × 360 with reduced precision by using a first order instead of a second order spatial interpolation. This clearly increased the wavelength (by about 40%) indicating that the problem is caused by the numerical viscosity.

Fig. 8 Histogram: colour coded is the logarithm of the probability for the occurrence of a wavelength at a radius normalised at each radius by the sum of all wavelengths for the specific radius. The black lines are proportional to the radius to the power of 2.5 and the lines are a factor of 2 apart from each other. The dashed line has linear slope. One can see that the instability jumps successively between different modes for the wavelength with corresponding jumps in frequency at the same radius.

Figure 7 indicates a strong reduction of the wavenumber with radius. To further explore this dependence of the wavelength on the radius, we performed an additional simulation with an extended radial domain from 2 AU to 50 AU. Again, we estimate the wavelength by measuring the distance between two sign changes in the vertical averaged vertical momentum. This time we show all the wavelengths that were detected by this method in Fig. 8, where we show how often a certain wavelength was captured, normalised to the specific radius where it was measured. An interesting behaviour can be observed. While the global radial wavelength does indeed depend linearly on the radius, locally it clearly deviates from this dependence and instead depends on the radius to the power of 2.5. This can also be seen in the simulation with smaller domain, but there it cannot be clearly distinguished from the interaction with the boundary.

This supplies us with an explanation for the resolution dependence of the instability. Since the modes cannot become arbitrarily small (because of the finite grid) or large (because of the limited vertical scale hight) there will be jumps between different modes. The viscosity and the Kelvin-Helmholtz instability, which can be observed in the simulations with high resolution, are the candidates for a physical cause for this cut off at small wavelengths.

Fig. 9 Fourier power spectrum of the temporal evolution of the instability after saturation (see Fig. 10). Analysed is the averaged meridional momentum of the simulation without viscosity and resolution 1024 × 256. Colour coded is the logarithm of amplitude of the frequency.

Because of the radius dependence of the wavenumber a spatial Fourier transform is not applicable. Additionally, as we show below, the wavelength in radial direction is not constant in time and phase jumps also can occur. However, to obtain more insight into the dynamics of the system, we show in Fig. 9 the results of a Fourier analysis in time of the vertical momentum of the simulation with resolution of 1024 × 256 (Fig. 10, along the time axis). To reduce the problems that phase jumps pose for the analysis (see below), we step through the data with a Hanning window over 1000 orbits and then average over those 5% of the resulting spectra that give the highest amplitude. We can see a dominant frequency at 0.022Ω_K at the inner region; this frequency is halved at the outer region beginning at about 5 AU. These jumps in the frequency domain coincide with the jumps in wavenumber. When the wavenumber jumps up, the frequency falls down, indicating an inverse relationship. On each branch the frequency is constant, while the wavenumber varies as ∝r^-2.5. We can understand this relationship starting from Eq. (12) from which one obtains for stable inertial oscillations (see Eq. (36) in Nelson et al. 2013) ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{~}\mathrm{-}{\mathrm{\Omega }}^{\mathrm{2}}\frac{{\mathit{k}}_{\mathit{Z}}^{\mathrm{2}}}{{\mathit{k}}_{\mathit{R}}^{\mathrm{2}}}\mathrm{·}$(16)The vertical scale is given by the local disc’s scale height H ~ r and hence k_Z ~ r^-1. In the quasistationary phase, we observed k_R ~ r^-2.5 (see Fig. 8), leading to an oscillation frequency independent of the radius, which we also observed (see Fig. 9).

Fig. 10 Large scale time development of the instability. Shown is the vertically averaged momentum in the meridional direction for the inviscid isothermal simulation with a resolution of 1024 × 256 (red curve in top panel of Fig. 7).

To obtain further insight into the spatio-temporal behaviour of the flow dynamics in Fig. 10 we show the vertically averaged momentum in the vertical direction as a function of space and time. In this global overview we observe waves that appear to travel slowly from larger to smaller radii. As noticed already in the Fourier analysis in Fig. 9, there is a transition between 4–5 AU with a change in wavelength of the perturbations and occasional phase jumps. Coupled to this is a change in the typical inward speed of the waves. They move more slowly when farther away from the star. As inferred roughly from Fig. 10, the wave speed at r = 6 AU is about 0.5 AU per 250 orbits, while at 4 AU it is about 1 AU. However, there is some dependence of this speed on time and space.

Near the outer boundary we sometimes see a region with standing waves, indicating that the radial domain should not be too small. This region is mostly only a few wavelengths in size (less then 1 AU), but can sometimes also reach a few AU into the domain. Reflections with the outer boundary play a role here as well, as can be seen in Fig. 10 for example at t ≈ 500 or 3600. We note that in contrast to our treatment at the inner radial boundary, we did not apply a damping region at the outer boundary.

To check if the viscosity is important for the wavelength, we add a small viscosity of ν = 5 × 10^-7. As expected, this leads to a wavenumber that is independent of the resolution, as shown in the bottom panel of Fig. 7. The wavelength is of the order of 0.2 AU at a radius of 4 AU after the instability is saturated.

With the wavelength fixed, the Reynolds stress also shows no strong dependence on the resolution as can be seen in Fig. 11 top panel. The inner region is strongly suppressed because we also increased the damping from 2 AU to 3 AU. With that we conclude that a small viscosity is necessary in order to introduce a physical lengthscale for the smallest unstable wavelength. To further explore the role of viscosity we repeat the simulation for different viscosities. This is done with a resolution of 1440 × 360. The growth rate is then calculated by fitting a linear function to the logarithm of the kinetic energy, which was at each point averaged over 100 grid cells. The results in the lower panel of Fig. 11 indicate that for the two lowest viscosities (10^-8 and 10^-7) the stresses are given by the numerical viscosity. For the intermediate case (10^-7) the stresses are larger while for very large values the effect of the increased damping near the inner boundary influences the results.

5. Discs with radiation transport

The isothermal discs discussed above do not capture the full physics, and most importantly the transport of energy is missing. In this section we include radiative transport and the heating/cooling interaction of the gas with the radiation. In the first set of models we start from the isothermal models as described above and switch on the radiation according to Eqs. (3) and (4); in a second series of models (in Sect. 6) we include irradiation from the central star.

For the simulations with radiative transport we use a resolution of 1024 × 256 and the same spatial extent and initial conditions as in the isothermal case. In Fig. 12 we show the midplane temperature averaged from 4–5 AU and the meridional kinetic energy when radiation is included, for two different values of the disc density ρ₀. In both cases the kinetic energy initially has larger amplitudes than in the previous isothermal simulations because now the disc is no longer in hydrostatic equilibrium initially, and small motions in the meridional plane set in (lower panel in Fig. 12). For the same disc density as before, ρ₀ = 10^-10, the disc cools off quickly as soon as the instability begins to be active, at around t = 10. The reason lies in the efficient radiative cooling in this case, in particular near the surface layers where the optical depth is small and the instability most active. Hence, any turbulent heating will be radiated away quickly.

We repeated the simulation with a higher density, ρ₀ = 10^-9 at 1 AU, to increase the optical thickness. Now the disc does not cool efficiently enough, and the instability begins to set in between t = 10 and t = 20 orbits, very similar to the isothermal models, but then radiative cooling eventually leads again to a cooling of the disc and the instability dies out. From these results it is clear that the instability does not produce enough heat and cannot survive without an external source of heat, for the typical opacities and densities expected in protoplanetary discs. This potential problem was pointed out already by Nelson et al. (2013).

Fig. 11 Vertically averaged Reynolds stress divided by the vertically averaged pressure in the saturated phase. Upper panel: for a viscosity of 5 × 10-7 at different numerical resolutions. Lower panel: fixed resolution of 1440 × 360 and different viscosities. Both are averaged from orbits 1000 to 3000.

Fig. 12 Discs with radiation transport for two different densities, ρ0. Upper panel: the midplane temperature at 4−5 AU as a function of time. Lower panel: kinetic energy in the meridional flow in the discs.

6. Irradiated disc

Here, we extend our models and include irradiation from the central star as an external heat source. Of course there are also other sources possible, for example, the inner region of the disc where the MRI is still active could be important.

6.1. Method of irradiation

We use a simple model for the external heating and consider vertical irradiation from above and below the disc, where the energy flux, F^irr, depends on radius. This procedure avoids the problem of finding a self-consistent solution for the flaring of the disc, as done for an irradiated and internally heated disc by Bitsch et al. (2013).

To obtain a first approximation for the flux in the meridional direction we assume that the angle of incidence of the flux is approximately R_⊙/r, where R_⊙ is the star’s radius. This applies to an infinitely flat disc as well as to the upper and lower surfaces our computational grid in spherical polar coordinates because all three represent planes that cross the centre of the central star. We obtain for the meridional component of the flux ${\mathit{F}}_{\mathit{\theta }}^{\mathrm{irr}}\mathrm{=}{\mathit{F}}_{\mathit{r}}\hspace{0.17em}\frac{{\mathit{R}}_{\mathrm{\odot }}}{\mathit{R}}\mathit{,}$(17)where F_r = F_⊙(R_⊙/r)² is the radial flux from the star at a distance r. Applying this impinging vertical irradiation to the disc leads to a radial temperature profile exponent in the disc of q = −0.55, in good agreement with the models of Chiang & Goldreich (1997). Our procedure does not allow for self-shadowing effects (Bitsch et al. 2013) but should give a physically realistic estimate of the expected temperatures in the disc.

To simplify the calculations and obtain a first order estimate of the effect, we use for the irradiation opacity the Rosseland opacity of Bell & Lin (1994). Hence, in the simulations we use the same opacity for the irradiation, Rosseland and Planck opacity (Bitsch et al. 2013). Numerically, we perform a ray-tracing method to calculate the energy deposited in each cell of the computational grid (Kolb et al. 2013).

6.2. Growth rate

To measure the growth rates of the instability for discs with radiation transport and irradiation we ran models with zero viscosity and a higher resolution case with ν = 10^-7. To be able to compare the growth rates with the previous isothermal cases, we performed additional isothermal simulations using the temperature profile from the simulations with radiation transport and irradiation.

The growth rates for the instability in combination with radiation transport are difficult to capture because the simulation cannot be started in hydrostatic equilibrium because the equilibrium vertical profile is unknown. We use strong damping for the first ten orbits to remove the disturbance caused by the transition to the new density and temperature profile.

Fig. 13 Growth rates for models with radiation transport and irradiation (irr) and isothermal models (iso) that have the same mean temperature, for comparison. For each radius the kinetic energy was smoothed over a range of 10% of its radius before fitting it to an exponential growth.

Fig. 14 Growth of the kinetic energy in the 2D plane with radiation transport and isothermal for comparison.

The results are shown in Fig. 13. We note that this time the growth rates should depend on radius because the growth depends on H/R which is not constant in the radiative cases. From Fig. 13 it is clear that the growth rates for the isothermal models are now lower than in the cases presented above, first because the temperature is lower and second because the radial profile is flatter as before, and both are important for growth. For the irradiated models the growth is again lower, with 0.1−0.2 per local orbit around half the value for the isothermal case. In Fig. 14 we show the evolution of the kinetic energy for the irradiated and corresponding isothermal model. For the inviscid case the final saturation levels agree very well with each other, while for the viscous disc with ν = 10^-7 the instability is weaker in the inner regions of the disc (see below).

6.3. Quasistationary phase

Fig. 15 Upper panel: temperature profile for an irradiated disc in the saturated phase without viscosity and a resolution of 1024 × 256. The dotted line is a run without hydrodynamics, only solving for the radiation energy. Lower panel: vertically integrated optical depth.

In the top panel of Fig. 15 we show the vertical temperature distribution for the saturated state at different radii in the disc for the model without viscosity at a resolution of 1024 × 256. The other models look very similar. In the bulk part of the disc the profile is quite flat with a slight drop towards the midplane. A lower central temperature is to be expected for irradiated dics, detailed radiative transfer models indicate an even larger temperature drop towards the midplane (Dullemond et al. 2002). In the upper layers the temperature falls off because the disc is optically thin and the energy can freely leave the system. This drop of the temperatures towards the surface despite the irradiation is a result of the identical irradiation and Rosseland opacity. If more radiation is allowed to be absorbed in the disc by increasing the irradiation opacity then one can obtain hotter surface layers. For a ten times larger value we find a hot corona similar to Flock et al. (2013) and a cooler midplane. First results seem to indicate a reduction in the Reynolds stress in this case, probably caused by the lower temperature in the bulk of the disc. At this point we leave the details to subsequent studies. The dotted line in Fig. 15 shows the profile for a simulation where we only solve for the radiation energy and disable the hydrodynamic solver. We can infer from this that the flat profile is a result of the combination of turbulent heating and vertical motion. A test simulation with a passive tracer added in the midplane of the disc in the saturated state showed rapid spreading over the whole vertical extent of the disc.

The vertically integrated optical depth is shown in the lower panel, starting from very small values at the disc surfaces it reaches 30–100 at the different radii. The nearly constant vertical temperature within the disc motivates us to use the equilibrium azimuthal velocity for the corresponding isothermal model of the steady state to calculate the Reynolds stress. In Fig. 16 we can see that it is still a good approximation. We note that this time the comparison is done with a 2D-simulation. Additionally shown is the Reynolds stress calculated with the Kepler velocity instead of the equilibrium velocity.

Fig. 16 Irradiated run. The Reynolds stress was averaged over 41 time steps, from orbit 1000 to 5000 each step 100 orbits apart, calculated with different averaging methods. For “mean time” the mean uφ was calculated through averaging over 40 time steps. For “Kepler” the velocity was calculated by subtracting the Kepler velocity and for the last one uφ is calculated analytically by using the equilibrium Eq. (11). Spatial averages are taken from 4 AU to 10 AU.

Fig. 17 Comparison of the Reynolds stress divided by the mean pressure, averaged over 2000 orbits, beginning with orbit 2000; “irr” stands for the irradiated disc, “iso” for the isothermal disc with analogous initial conditions.

While the growth rates are weaker than in the isothermal case, the kinetic energy in the meridional plane for stable saturated phase in Fig. 17 reaches the same level with radiation transport. The values for α are again between 0.5 × 10^-4 and 2 × 10^-4, depending on the wavelength and thus viscosity, but independent of radiation transport.

The strength of the instability measured in terms of the value of the viscosity under which it still survives is, of course, different. Here, in the irradiated case, even a low viscosity of 10^-7 suppresses the instability in the inner regions of the disc. This is not only a result of the radiation transport, but also of the flat temperature profile. The details will depend on the source of the heating and the opacity, but nevertheless the stability will be weaker than in the purely isothermal case.

Fig. 18 Radiative diffusion time per Orbit at 4 AU for a lengthscale of 0.1 AU.

6.4. Discussion

As we have shown in the previous sections for an irradiated disc there is the possiblity of generating an effective turbulence through the VSI. As pointed out in Nelson et al. (2013) the instability can only be sustained if the diffusion (local relaxation) time is a fraction of the local orbital period. To investigate how this condition is fulfilled in our simulations we analyse the timescale for radiative diffusion for the equilibrium irradiated disc models. In units of the local orbital period this is given by ${\mathit{t}}_{\mathrm{diff}}\mathrm{=}\mathrm{\Delta }{\mathit{x}}^{\mathrm{2}}\frac{{\mathit{c}}_{\mathrm{P}}{\mathit{\rho }}^{\mathrm{2}}{\mathit{\kappa }}_{\mathrm{R}}}{\mathrm{4}\mathit{\lambda ac}{\mathit{T}}^{\mathrm{3}}}\mathrm{·}\frac{\mathrm{\Omega }}{\mathrm{2}\mathit{\pi }}\mathit{,}$(18)where Δx is the characteristic wavelength of the perturbation. In our case the radial diffusion is relevant (Nelson et al. 2013) and we choose here Δx = 0.05r, which is a typical radial wavelength at r = 3 AU. Using Eq. (18) and the results from the simulation we calculate for the optical thin region a very small cooling time per orbit of t_diff = 10^-10 as expected. For the optically thick region we obtain t_diff = 0.11 for our standard density, which is indeed a small fraction of the orbital period as required for the instability to operate, see Fig. 18. The cooling time in the vertical direction is longer, about a few orbital periods as implied by the vertical optical depth (see Fig. 15), but this will keep the disc nearly isothermal, again as required for instability.

In Fig. 19 we illustrate that the instability still resembles closely the locally isothermal case except that the small scale perturbations are missing, even in the optically thin region, where we have very short cooling times. For comparison, Nelson et al. (2013) found that the instability was completely suppressed with relaxation times of t_relax = 0.1, which is the timescale for the flow to relax to the initial isothermal profile. We take this as an indication that physically, radiative diffusion plus irradiation behaves in a different way from a simple model of temperature relaxation as used in Nelson et al. (2013).

As seen in Fig. 15 an increase in the density leads to higher optical depths and longer diffusion times, and consequently to a weaker instability. While doubling the density in a simulation with resolution 2048 × 512 has no clear influence on the kinetic energy and the cooling times in the optical thin regions, the Reynolds stress was clearly weaker by a factor of around 1.5 in the simulation with doubled density (the model in the middle of Fig. 18). In addition the wavelength of the perturbations is decreased.

Fig. 19 Velocity in the meridional direction, uθ, in units of local Kepler velocity for an irradiated run without viscosity at resolution 1024 × 256 (top) and with resolution 2048 × 512 (bottom). Compare with Fig. 2 which has a spatial resolution of 2048 × 512.

A further increase in the density leads also to a strong decrease in the kinetic energy, with again a smaller wavelength. This raises the question whether the simulation with resolution of 2048 × 512 is sufficiently resolved. These results indicate that in very massive discs with long diffusion times (vertical and radial) the disc will behave more adiabatically, and the instability will be quenched. The minimum solar mass nebula at 5 AU corresponds approximately to our model with 2 ρ₀ and the instability might just be operative.

7. Summary and conclusions

We have studied the vertical shear instability as a source of turbulence in protoplanetary discs. For that purpose we have performed numerical simulations solving the equations of hydrodynamics for a grid section in spherical polar coordinates. To study the global behaviour of the instability we have used a large radial extension of the grid ranging from 2 AU to 10 AU.

In a first set of simulations we show that the instability occurs for locally isothermal discs where the radial temperature gradient is a given function of radius. Our results on the growth rates for the instability are in good agreement with the theoretical estimates by Urpin & Brandenburg (1998) and Urpin (2003), and we find two basic growth regimes for the asymmetric and antisymmetric modes as seen by Nelson et al. (2013). After 20 to 30 local orbits the instability saturates and is dominated by the vertical motions, which cover the whole vertical extent of the disc.

Interestingly, we find that the local radial wavelength of the perturbations scales approximately with λ ∝ r^2.5 in the saturated state with a constant frequency. However, on a global scale several jumps occur where the wavelengths are halved, such that the global scaling follows with . We suspect that the instability has the tendency to generate global modes that show the observed wavelength behaviour according to Eq. (16). Because of the radial stratification of the disc, jumps have to occur at some locations.

The waves approximately keep their shape and travel slowly inwards. The two- and three- dimensional simulations yield essentially the same results concerning the growth rates and saturation levels of the instability because of its axisymmetric property. The motions give rise to a finite level of turbulence and we calculate the associated efficiency, measured in terms of α. We first show that, caused by the two-dimensionality, α can be measured directly from the two-dimensional simulations using the proper equilibrium state of the disc. We find that the angular momentum associated with the turbulence is positive and reaches α-values of a few 10^-4. For the isothermal simulations we find that at higher numerical resolution α becomes smaller, but viscous simulations indicate a saturation at a level of about α = 10^-4 even for very small underlying viscosities that are equivalent to α< 10^-6.

Adding radiative transport leads to a cooling from the disc surfaces and the instability dies out subsequently. We then constructed models where the disc is irradiated from above and below which leads to a nearly constant vertical temperature profile within the disc. This again leads to a turbulent saturated state with a similar transport efficiency to the purely isothermal simulations, or possibly slightly higher (see Fig. 17).

In summary, our simulations indicate that the VSI can indeed generate turbulence in discs albeit at a relatively low level of about a few times 10^-4. This implies that even in (magnetically) dead zones the effective viscosity in discs will never fall below this level. Our results indicate that in fully 3D simulations the transport may be marginally larger, but further simulations will have to be performed to clarify this point.

Acknowledgments

Moritz Stoll received financial support from the Landesgraduiertenförderung of the state of Baden-Württemberg. Wilhelm Kley acknowledges the support of the German Research Foundation (DFG) through grant KL 650/8-2 within the Collaborative Research Group FOR 759: The formation of Planets: The Critical First Growth Phase. Some simulations were performed on the bwGRiD cluster in Tübingen, which is funded by the Ministry for Education and Research of Germany and the Ministry for Science, Research and Arts of the state Baden-Württemberg, and the cluster of the Forschergruppe FOR 759 “The Formation of Planets: The Critical First Growth Phase” funded by the DFG.

References

All Figures

	Fig. 1 Kinetic energy of the motion in the meridional plane at different radii in an inviscid disc. The kinetic energy at the different locations is in each case averaged over a radial interval with length 0.5 AU. We note that the unit of time is given in local periods at the centre of the specified interval. Hence, it is different for each curve, but this allows easy comparison.
In the text

	Fig. 2 Velocity in the meridional direction, uθ, in units of local Kepler velocity for an isothermal run without viscosity. The panels refer to snapshots taken at time 100, 210 and 750 (top to bottom), measured in orbital periods at 1 AU. In units of local orbits at (2.5, 3.5, 4.5) AU this refers to (25, 15, 10) (53, 32, 22) (190, 115, 79) orbits, from top to bottom.
In the text

	Fig. 3 Growth of the kinetic energy for the quarter of a disc and the 2D equivalent. The kinetic energy is averaged from 4 AU to 5.5 AU.
In the text

	Fig. 4 Vertical velocity in the midplane of the disc for the 3D model after 4000 orbits. The nearly axisymmetric property of the instability is clearly visible.
In the text

Fig. 5 Reynolds stress (code units) from 3−10 AU averaged over 41 time steps, each step 100 orbits apart beginning with orbit 1000, calculated with different averaging methods. For “mean time” the steady-state , needed to calculate the Reynolds stress at each step, was calculated through averaging over the 41 time steps. For “mean phi” the steady-state velocity was calculated by averaging over the azimuthal direction at each time step. For “equilibrium” is calculated analytically by using the equilibrium Eq. (11) at each step. For the 2D model we used the equilibrium method as well.

In the text

	Fig. 6 Radial and vertical distribution of the Reynolds stress. Left: the Reynolds stress divided by the midplane pressure over the vertical direction. Right: reynolds stress divided by the mean pressure over the radius for different resolutions. Both are averaged over 4001 time steps from orbit 1000 to 5000. The model res2048 corresponds to the results shown in Fig. 1.
In the text

	Fig. 7 Mean wavenumber of the instability over the radius for different numerical resolutions in the saturated phase. Upper panel: inviscid case with ν = 0, lower panel: viscous case with ν = 5 × 10-7 (dimensionless).
In the text

Fig. 8 Histogram: colour coded is the logarithm of the probability for the occurrence of a wavelength at a radius normalised at each radius by the sum of all wavelengths for the specific radius. The black lines are proportional to the radius to the power of 2.5 and the lines are a factor of 2 apart from each other. The dashed line has linear slope. One can see that the instability jumps successively between different modes for the wavelength with corresponding jumps in frequency at the same radius.

In the text

	Fig. 9 Fourier power spectrum of the temporal evolution of the instability after saturation (see Fig. 10). Analysed is the averaged meridional momentum of the simulation without viscosity and resolution 1024 × 256. Colour coded is the logarithm of amplitude of the frequency.
In the text

	Fig. 10 Large scale time development of the instability. Shown is the vertically averaged momentum in the meridional direction for the inviscid isothermal simulation with a resolution of 1024 × 256 (red curve in top panel of Fig. 7).
In the text

	Fig. 11 Vertically averaged Reynolds stress divided by the vertically averaged pressure in the saturated phase. Upper panel: for a viscosity of 5 × 10-7 at different numerical resolutions. Lower panel: fixed resolution of 1440 × 360 and different viscosities. Both are averaged from orbits 1000 to 3000.
In the text

	Fig. 12 Discs with radiation transport for two different densities, ρ0. Upper panel: the midplane temperature at 4−5 AU as a function of time. Lower panel: kinetic energy in the meridional flow in the discs.
In the text

	Fig. 13 Growth rates for models with radiation transport and irradiation (irr) and isothermal models (iso) that have the same mean temperature, for comparison. For each radius the kinetic energy was smoothed over a range of 10% of its radius before fitting it to an exponential growth.
In the text

	Fig. 14 Growth of the kinetic energy in the 2D plane with radiation transport and isothermal for comparison.
In the text

	Fig. 15 Upper panel: temperature profile for an irradiated disc in the saturated phase without viscosity and a resolution of 1024 × 256. The dotted line is a run without hydrodynamics, only solving for the radiation energy. Lower panel: vertically integrated optical depth.
In the text

Fig. 16 Irradiated run. The Reynolds stress was averaged over 41 time steps, from orbit 1000 to 5000 each step 100 orbits apart, calculated with different averaging methods. For “mean time” the mean uφ was calculated through averaging over 40 time steps. For “Kepler” the velocity was calculated by subtracting the Kepler velocity and for the last one uφ is calculated analytically by using the equilibrium Eq. (11). Spatial averages are taken from 4 AU to 10 AU.

In the text

	Fig. 17 Comparison of the Reynolds stress divided by the mean pressure, averaged over 2000 orbits, beginning with orbit 2000; “irr” stands for the irradiated disc, “iso” for the isothermal disc with analogous initial conditions.
In the text

	Fig. 18 Radiative diffusion time per Orbit at 4 AU for a lengthscale of 0.1 AU.
In the text

	Fig. 19 Velocity in the meridional direction, uθ, in units of local Kepler velocity for an irradiated run without viscosity at resolution 1024 × 256 (top) and with resolution 2048 × 512 (bottom). Compare with Fig. 2 which has a spatial resolution of 2048 × 512.
In the text