A&A 426, 755-765 (2004)
DOI: 10.1051/0004-6361:20035896

Simulations of vertical shear instability in accretion discs

R. Arlt¹ - V. Urpin²

1 - Astrophysikalisches Institut Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany
2 - A.F. Ioffe Institute for Physics and Technology, 194021 St. Petersburg, Russia

Received 18 December 2003 / Accepted 20 June 2004

Abstract
The nonlinear evolution of the vertical shear instability in accretion discs is investigated using three-dimensional hydrodynamic simulations. A vertical dependence of the angular velocity destabilizes the disc and leads to the generation of velocity fluctuations enhancing the angular momentum transport. The instability emerges in the numerical models for large radial perturbation wave numbers. The growth time is a few tens of orbital revolutions.

Key words: accretion, accretion disks - hydrodynamics - instabilities - turbulence

1 Introduction

The standard accretion disk implies sufficiently strong turbulence to enhance angular momentum transport. Differential rotation has ever been suspected of the most promising source of turbulence despite the radial dependence of the angular velocity satisfies the Rayleigh stability criterion. The conditions in astrophysical disks, however, have a number of differences compared to simple shear flows (vertical and radial stratification, baroclinity, vertical shear, etc.). In such conditions, the Rayleigh criterion of stability is of limited applicability in accretion disks. Therefore, some linear instabilities can occur even if the Rayleigh criterion is fulfilled.

The origin of turbulence in disks is often attributed to the well known magneto-rotational instability (Velikhov 1959; Kurzweg 1963) which can arise in magnetic, differentially rotating fluids. The properties of this instability has been analyzed in detail for stellar conditions (see Fricke 1969; Acheson 1978, 1979) where the instability becomes operative if the angular velocity decreases from the pole to the equator. The role of the magneto-rotational instability in the turbulization of astrophysical disks has been recognized by a number of authors (e.g. Safronov 1969; Balbus & Hawley 1991). Three-dimensional numerical simulations (e.g. Balbus & Hawley 1998, and references therein; Brandenburg et al. 1995; Torkelsson et al. 1997; Arlt & Rüdiger 2001) indicate that the magneto-rotational instability can provide an efficient mechanism of the angular momentum transport.

Most likely, however, the magneto-rotational instability is not the only instability that occurs in such complex bodies as accretion disks. Turbulence of hydrodynamic origin is an option which still demands critical examination. The problem of instability for increasing angular momentum with axis distance was re-examined by Richard & Zahn (1999) - the interesting case for astrophysical disks. Longaretti (2002) used phenomenological arguments to suggest that shearing sheet flows should be turbulent and that the lack of turbulence in simulations is likely caused by a lack of resolution. Note, however, that an important paper is neglected in these publications which shows analytical and experimental stability for increasing angular momentum in a Taylor-Couette experiment (Schultz-Grunow 1956). More work on that subject is required and is already in progress, e.g. at CEA, Saclay.

It seems that hydrodynamic stability of accretion disks is not well understood. For example, the vertical stratification despite being usually stable can provide a catalyzing effect which under certain conditions induces a linear non-axisymmetric instability of anticyclonically sheared flow (Molemaker et al. 2001). In a recent paper, Dubrulle et al. (2002) re-examined this instability by making use of a different approach and found that it can occur in accretion disks. Klahr & Bodenheimer (2003) have found that accretion disks can also be subject to an instability caused by a negative radial entropy gradient and resulting in vigorous turbulence enhancing angular momentum transport. A sound wave instability was found by Drury (1985) and Papaloizou & Pringle (1984) which requires a reflecting boundary. Since we do apply reflecting boundaries, we will come back to this instability in the one but last section of our paper.

One more potential source of instability - the one studied here - is provided by a dependence of the angular velocity on the vertical coordinate. This dependence seems to be inevitable in accretion disks (see, e.g., Kippenhahn & Thomas 1982; Urpin 1984; Kley & Lin 1992). It has first been argued by Kippenhahn & Thomas (1982) that a slight baroclinity is necessary to fulfill hydrostatic and thermal equilibrium in accretion disks. The dependence of $\Omega$ on z is relatively weak but it is a driving force for the well known Goldreich-Schubert  instability (Goldreich & Schubert 1967). Recently, hydrodynamic stability of accretion disks under the action of the vertical shear has been considered by Rüdiger et al. (2002). The authors considered linear and non-linear small scale perturbations, using the ZEUS-3D code, and concluded that the vertical shear cannot overcome the stabilizing effect of stratification if perturbations are adiabatic. However, small scale motions in disks are essentially non-adiabatic, and the exchange of heat between perturbations and the surrounding medium substantially reduces the influence of the buoyancy force and decrease the stabilizing effect of stratification. Therefore, the behavior of perturbations in adiabatic and non-adiabatic cases is qualitatively different (Urpin & Brandenburg 1998; Urpin 2003).

In the present Paper, we perform simulations of the vertical shear instability in accretion disks in a non-adiabatic approximation. The ZEUS-3D code was used for the integration of the hydrodynamic equations. We calculate the main parameters of turbulence caused by the vertical shear instability and argue that the vertical shear existing in accretion disks is sufficient to enhance turbulent angular momentum transport.

The paper is organized as follows: in Sect. 2, the properties of the instability are briefly reviewed and discussed, while the numerical setup is described in Sect. 3. The results from simulations with and without perturbations are presented in Sect. 4. Our findings are summarized in Sect. 5.

2 Basic equations and the vertical shear instability

We consider a non-magnetic axisymmetric disk of finite vertical extent. The unperturbed angular velocity can generally depend on both r and z coordinates, $\vec\Omega = (\Omega(r, z),0,0)$. Throughout this paper, z, r and $\varphi$ are cylindrical coordinates; the unit vectors in these directions are $\vec{e}_z$, $\vec{e}_r$, and  $\vec{e}_\varphi$. In the unperturbed state, the disk is assumed to be in hydrodynamic equilibrium,

$$% \frac{\nabla p}{\rho} = \vec{G} = \vec{g} + r\Omega^{2} \vec{e}_r,$$ (1)

where $\vec{g}$ is the gravity of the central object. Solving Eq. (1) together with the thermal balance equation, one can generally obtain r- and z-dependences of $\rho$, p, and $\Omega$. It turns out that $\Omega$ always depends on z in the unperturbed state (see e.g. Kley & Lin 1992). For the sake of simplicity, we assume the disk to be isothermal.

In the present paper, we treat numerically only short wavelength disturbances with a wavelength much shorter than the half-thickness of a disk, H. The amplitude of such disturbances can grow due to the vertical shear instability, and they may evolve in a non-linear regime but still remaining short-wave. Short wavelength disturbances are approximately isothermal in accretion disks because of the high thermal conductivity. The isothermal approximation applies for disturbances with the wave vector $\vec{k}=(k_{z}, k_{r}, k_{\varphi})$ satisfying the condition

$$% k > k_{\rm min} = \sqrt{\gamma_{\rm vs}/ \chi},$$ (2)

where $\gamma_{\rm vs}$ is the growth rate of the vertical shear instability, and $\chi$ is the thermal diffusivity. In accretion disks, we have $\gamma_{\rm vs} \sim \Omega H/r \sim c_{\rm s}/r$ (see Urpin 2003) where $c_{\rm s}$ is the sound speed. Then, we have from Eq. (2)

$$% k H \gg k_{\rm min} H \sim \left( \frac{H}{r} \right)^{1/2} \left( \frac{c_{\rm s} H}{\chi} \right)^{1/2}\cdot$$ (3)

In disk models, $\chi$ is typically large, and disturbances with $k H \gg 1$ can usually be treated as isothermal.

The behavior of isothermal disturbances is governed by the continuity and momentum equations

$$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho {\vec u}) = 0$$ (4)

$$\rho \frac{\partial \vec{u}}{\partial t}+\rho (\vec{u} \cdot \nabla) \vec{u} = -\nabla p + \rho \vec{g}.$$ (5)

where ${\vec u}$ is the velocity of disturbances. Our consideration does not include self-gravity in the disk.

For the purpose of illustration, we show in this section that the model governed by Eqs. (4) and (5) leads to the vertical shear instability in a linear regime. Consider the behavior of axisymmetric short-wavelength disturbances with the space-time dependence $\exp(\gamma t - i \vec{k} \cdot \vec{r})$ where $\vec{k}= (k_{z}, k_{r}, 0)$ is the wave vector and $\vec{r} = (z, r, \varphi)$ is a locus; $\vert \vec{k} \cdot \vec{r} \vert \gg 1$. Small perturbations will be marked by the subscript 1, whilst unperturbed quantities will have no subscript. The unperturbed velocity is represented by rotation with the angular velocity $\Omega$. Then, the linearized Eqs. (4) and (5) read

$$% \gamma \rho_{1} = i \rho (\vec{k} \cdot \vec{u}_{1}),$$ (6)

$$% \gamma \vec{u}_{1} + 2 \vec{\Omega} \times \vec{u}_{1} + \v... ... = -\frac{\nabla p_{1}}{\rho} + \vec{G} \frac{\rho_{1}}{\rho},$$ (7)

where $\vec{u}_{1}$, p₁ and $\rho_{1}$ are perturbations of the hydrodynamic velocity, pressure and density, respectively. In our isothermal model, $p_{1} = c_{\rm s}^{2} \rho_{1}$. Then, replacing p₁ and $\rho_{1}$ in Eq. (7), we obtain an equation containing  $\vec{u}_{1}$ alone

$$% \gamma \vec{u}_{1} + 2 \vec{\Omega} \times \vec{u}_{1} + \v... ... \vec{u}_{1}) \left(i \vec{k} c_{\rm s}^{2} + \vec{G}\right).$$ (8)

Under the assumption that the growth rate of the vertical shear instability is small compared to the sound speed, $k c_{\rm s} \gg \gamma$, Eq. (8) yields the following simple dispersion relation,

$$% \gamma^{2} = - Q^{2},$$ (9)

where

$$% Q^{2} = 4 \Omega^{2} \frac{k_{z}^{2}}{k^{2}} + 2 \Omega r \... ...ial r} - k_{r} \frac{\partial \Omega}{\partial z} \right)\cdot$$ (10)

The first term on the r.h.s. of Eq. (10) is equal to the square of the frequency of so-called inertial waves that can exist even in a rigidly rotating fluid, and the second term represents the effects associated with a departure from rigid rotation.

The instability occurs if

Q² < 0, (11)

or

$$% \frac{k_{z}^{2}}{k^{2}} \cdot \frac{\partial}{r^{3} \parti... ...2}} \cdot 2 \Omega r \frac{\partial \Omega} {\partial z} < 0.$$ (12)

For thin disks, the radial dependence of $\Omega$ is approximately given by the Keplerian law, $\Omega \propto r^{-3/2}$, and hence, the first term on the r.h.s. of the inequality (12) is positive. The sign of the second term depends on the direction of the wave vector, and this term may cause a destabilizing effect. In thin disks with $z/r \ll 1$, the vertical dependence of $\Omega$ can be calculated analytically expanding all quantities around the midplane into a parabola (see e.g. Urpin 1984; Kley & Lin 1992). Then we have

$$% \partial \Omega / \partial z \approx q \Omega z/r^{2}$$ (13)

where $q=q(r) \sim 0.1$ is the parameter in the series expansion of  $\Omega(z,r)$ around the equator. The term associated with the vertical shear dominates the r.h.s. of the inequality (12) if

$$% k_{r} > \mid (k_{z}/2q) (r/z) \mid.$$ (14)

The instability arises only for perturbations with wavelengths much shorter in the radial direction than vertically.

Note that, for any dependence of $\Omega$ on z, there exists a region in the plane  (k_z, k_r) where the condition of instability, Q² < 0, is satisfied. In the general case, the instability arises if the components of a wave vector satisfy the inequality

$$% \left\vert \frac{k_{r}}{k_{z}} \right\vert > \frac{\vert\pa... ...rt}{\vert 2 \Omega r^{4} \partial \Omega / \partial z \vert},$$ (15)

and k_r and k_z are of the same/opposite sign if the quantity $[\partial (r^{4} \Omega^{2}) / \partial r]/ [\partial \Omega/ \partial z]$ is positive/negative.

Table 1: Overview of numerical configurations. All values are in the arbitrary units of the code.

This simple consideration shows that accretion disks are always unstable as soon as disturbances have a very short radial length-scale. The growth rate is given by

$$% \gamma \sim \vert Q\vert.$$ (16)

Since the condition (12) can be satisfied only for perturbations with $\vert k_{r}/ k_{z} \vert \gg 1$, the growth rate depends on the ratio  k_z/k_r rather than on the wavelength of perturbations. The maximum growth rate is reached at

$$% \frac{k_{z}}{k_{r}} \approx \frac{r^{4} \Omega ( \partial \... ...a / \partial z)}{\partial (r^{4} \Omega^{2})/ \partial r}\cdot$$ (17)

In a Keplerian disk, this ratio reads

$$% \frac{k_{z}}{k_{r}} \approx q \cdot \frac{z}{r},$$ (18)

and the corresponding maximum growth rate is given by

$$% \gamma_{\rm vs} \approx \Omega \left\vert \frac{q z}{r} \right\vert\cdot$$ (19)

The growth time of the instability is thus rather short and can even be as small as the order of magnitude of the time scale of vertical shear.

3 Numerical model

The setup for our three-dimensional simulations applies a global, cylindrical computational domain with full azimuthal range, dimensionless radii from 4 to 6 and from 3 to 7, and a vertical extent of 2 density scale heights on average, which translates to -1 to +1 in our models and in dimensionless units. Two main configurations provided the results of the following investigation; their details are given in Table 1.

The ZEUS-3D code (developed by Stone & Norman 1992a,b; Stone et al. 1992) was used for the integration of the compressible hydrodynamics. The configuration is very close to the one given in Arlt & Rüdiger (2001); we also compute in an isothermal gas. Isothermality does not provide a vertical shear unless meridional motions are present in the disk. In our simulations this will be the relaxation motions around the equilibrium solution and the imposed initial perturbations.

The induction equation is naturally not integrated here. The source of gravitation is that of a point mass M=10⁵ at r=0 and z=0. Self-gravity in the disk is not considered. The total mass of the disk is about 10^-3M or less. The Toomre criterion for gravitational stability is readily fulfilled with a Toomre parameter of $\sim$10² in our cases.

The original ZEUS code provides artificial viscosities for improved shock evolution; we have put the von-Neumann-Richtmyer viscosity to zero here since only subsonic flows are expected. The advection interpolation is the second order van-Leer scheme. The additional linear smoothing by the ZEUS variable  $q_{\rm lin}$ is set to 0.2 for improved stability of the unperturbed equilibrium configuration. The perturbed model runs apply the same value of course.

The conditions for the vertical boundaries are stress-free; no matter can exit the computational domain in vertical direction and the vertical derivatives of u_r and $u_\varphi$ vanish.

The radial boundaries of the two models discussed in detail (S and L) are also closed for the flow, and the radial derivative of u_z vanishes. This condition differs from the setup in Arlt & Rüdiger (2001) where boundaries allowed for accretion. Our last model labelled with "O'' uses a simple version of open radial boundaries, where the values of the last grid cell before the boundary are just copied to the ghost zones. This attempt is rather crude, and we will thus concentrate on the closed models S and L in the following.

For radial boundary conditions on the azimuthal flow, we have to consider the intrinsic rotation due to the central mass. A Keplerian radial dependence of rotation is maintained into the r-boundary zones, based on the innermost (outermost) zone of the integration domain. If the azimuthal velocity at the innermost point of the computational domain is denoted by  $u_\varphi^{(1)}$, we use

$$u_\varphi^{(0)}\phantom{^{-}} = u_\varphi^{(1)} \sqrt{r^{(1)}/r^{(0)}}$$ (20)

$$u_\varphi^{(-1)} = u_\varphi^{(1)} \sqrt{r^{(1)}/r^{(-1)}}$$ (21)

to find azimuthal-velocity values for the boundary zones inside the annulus of the computational domain. The outer radial boundary is treated accordingly for  $u_\varphi^{(n+1)}$ and  $u_\varphi^{(n+2)}$. We did not fix the absolute value of Keplerian rotation, since radial pressure gradients may alter the effective gravitational potential (usually decrease the revolution speed slightly in accretion disks). The azimuthal boundary conditions are fully periodic.

A rough representation of a stratified disk with a linear radius dependence for the density scale-height was used for the initial conditions. The initial velocity field is a mere Keplerian rotation with  $\partial\Omega/\partial z=0$. We leave this configuration for free development under the influence of the gravitational potential. The meridional relaxation motions are enough to generate a slight vertical shear on which our simulations rely.

The constant sound speed is set to $c_{\rm s}=10$ in the arbitrary units of the code, corresponding to 7% of the Keplerian velocity in the middle of the computational domain, at r=5. When the equilibrium state is reached, the ratio of the half-thickness to the radius is $\approx$0.1. We will show an unperturbed model in the following Section. The damping of the small initial deviations in the starting configuration from the exact equilibrium solution is an indication for the stability of the numerics, but also shows the numerical diffusion resulting from the integration scheme.

The time step is nearly constant at 3 $\times$ 10^-4which is about  $0.0014~P_{\rm orb}$, the rotation period of the disk at r=5. It is limited by the velocity of sound waves, but the large Keplerian velocity also defines a time-step nearly as small as the above.

The initial conditions for the smaller of the two long runs involve an axisymmetric velocity perturbation of the form

  

$$u_z = A \sin (2\pi r/0.04)$$ (22)

$$u_r = A \cos (2\pi r/0.04).$$ (23)

The amplitude A is $0.01 c_{\rm s}$. Additionally, relaxation motions of the initial density stratification, whose scale-height changes linearly with r, lead to dependences on z with a wavelength of 2, i.e. the vertical extent of the computational domain. This is why Table 1 lists $k_z=\pi$ for model S.

For the larger resolution in $(z, r, \varphi)$ of 81 $\times$ 500 $\times$ 161, the initial perturbations were

   

$$u_z = A_1 \sin 2\pi z$$ (24)

$$u_r = A_1 \sin (2\pi r/0.064)$$ (25)

$$u_\varphi = A_2 \sin 2\varphi.$$ (26)

We thus combined the wavenumbers $k_z = 2\pi$ and k_r = 98 giving a ratio of k_r/k_z = 15.6. The amplitude A₁ is again chosen to provide perturbations of 1% of the thermal sound speed, $A_1=0.01 c_{\rm s}$. The azimuthal perturbation is added to allow for a quicker non-axisymmetric evolution of the disk; it does not significantly contribute to the instability. We set $A_2 = 10^{-4} c_{\rm s}$.

4 Results

4.1 Unperturbed models and diagnostics

Before actually describing the results for the instability, we would like to add some remarks on the unperturbed model. In Fig. 1, the time dependence of the average of the square of velocity fluctuations is plotted for the equilibrium model. For the sake of simplicity, we will often call this quantity "kinetic energy'' of the fluctuations, but we have to note that the density stratification is not involved in the measure. The values are scaled to the square of the constant sound speed, $c_{\rm s}^2$, in these plots. The average "energy density'' in the fluctuations is computed for the individual components. It is the average square of the velocity component minus the azimuthal average at the specific point in the (z,r)-plane:

$$% \langle u_z'^2\rangle=\left\langle\left[u_z(z,r,\varphi)- \overline{u_z}(z,r)\right]^2\right\rangle_{z,r,\varphi}$$ (27)

where $\overline{u_z}(z,r) = \langle u_z(z,r,\varphi)\rangle_\varphi$. The subscripts to "$\rangle$'' denote the directions of averaging. The other energies  $\langle u_r'^2\rangle$ and  $\langle u_\varphi'^2\rangle$ are computed likewise. This procedure excludes all axisymmetric parts of the flow from the average, mainly the differential rotation, but also in u_z and u_r.

The simulation shown in Fig. 1 does not involve any initial disturbances to the model S except slight deviations from the exact equilibrium solution. The typical kinetic energy in the fluctuations is extremely small and is decreasing all the time. The values reduce from $\sim$ $10^{-15}c_{\rm s}^2$ to $\sim$ $10^{-17}c_{\rm s}^2$, i.e. in units of thermal energy, during the time of 85 orbital revolutions of the inner edge of the computed annulus.

	Figure 1: Energies of the velocity fluctuations for the unperturbed model S in units of the thermal energy. The lower axis measures orbital periods of the inner edge at r=4, the upper axis at the outer edge at r=6.
Open with DEXTER

	Figure 2: The $\alpha _{\rm SS}$-parameter for the unperturbed model S. The average $\alpha _{\rm SS}$ is -4 $\times$ 10-11. The horizontal axes are as in Fig. 1.
Open with DEXTER

We may also check for the transport of angular momentum through the $r\varphi$-component of the stress tensor

$$% W_{r\varphi} = \frac{{\left\langle{\rho u_r [u_\varphi-\ove... ...ight\rangle_{z,r,\varphi}}}{\langle\rho\rangle_{z,r,\varphi}},$$ (28)

which is typically normalized to the Shakura-Sunyaev parameter  $\alpha _{\rm SS}$ by

$$% W_{r\varphi}=\alpha_{\rm SS} c_{\rm s}^2.$$ (29)

Figure 2 shows the temporal development of  $\alpha _{\rm SS}$. According to that, the angular momentum transport is vanishing with an average of  $\alpha_{\rm SS} \sim -10^{-11}$. The simulations illustrate the fact that the chosen equilibrium model is sufficiently stable during the computational time of >80 revolutions.

The unperturbed run for model L is convincing as well. The plot of the fluctuation energies in Fig. 3 shows an increase in the radial and vertical components initially. The deviations of the true density distribution from the initial one linear in r are larger than for model S and imply more rearrangement flows. After 30 orbital revolutions, the fluctuations settle and decrease by 8 orders of magnitudes in energy within 220 orbits of the inner edge. Again we look at  $\alpha _{\rm SS}$ which is plotted in Fig. 4. The transient maximum in the strength of velocity fluctuations does not feature in the $\alpha _{\rm SS}$-plot, and the average is again $\alpha _{\rm SS}$ $\sim$ 10^-11.

The zero models do build up gradients $\partial \Omega/\partial z$ as well, but the absence of perturbations with proper wave lengths does not lead to a vertical shear instability.

	Figure 3: Energies of the velocity fluctuations for the unperturbed run of model L. The lower axis measures orbital periods of the inner edge at r=3, the upper axis at the outer edge at r=7.
Open with DEXTER

	Figure 4: The $\alpha _{\rm SS}$-parameter for the unperturbed model L. The average $\alpha _{\rm SS}$ is -2 $\times$ 10-11. The horizontal axes are as in Fig. 3.
Open with DEXTER

  
4.2 Perturbed models

Now we go ahead with imposing the initial perturbations given in (22) and (23) to model S of Figs. 1 and 2. The time dependence of the fluctuation energies and  $\alpha _{\rm SS}$ is shown in Figs. 5 and 6 for the case when initial perturbations are imposed to the S model. Although the initial amplitude of the velocity perturbations is relatively small ( $0.01 c_{\rm s}$), the vertical shear instability leads to the generation of significant velocity fluctuations. During a long time of growth, the velocity and density patterns do not actually show turbulent features. We rather observe a growing large-scale flow. Remember that the unperturbed model showed a decay of velocity fluctuations by the intrinsic "numerical viscosity''.

The $\alpha _{\rm SS}$-parameter in Fig. 6 was again calculated by averaging over the computational box. This parameter is determined by the correlation between radial and azimuthal velocity fluctuations and is very small initially because the fluctuations are nearly axisymmetric in the beginning and cannot transfer the angular momentum. However, $\alpha _{\rm SS}$ exceeds a value of 10^-8 after $t \sim 30$ outer orbits. It is first positive, but still very small implying slight outward angular-momentum transport. After 150 orbital periods at the inner boundary or 80 revolutions at the outer boundary, the sign of  $\alpha _{\rm SS}$ is negative in model S and saturates at  $\alpha_{\rm SS}\sim -0.0005$.

	Figure 5: Energies in velocity fluctuations for model S.
Open with DEXTER

	Figure 6: The evolution of the $\alpha _{\rm SS}$-parameter for model S.
Open with DEXTER

	Figure 7: The growth rate for the square of the velocity fluctuations for model S. The graph is smoothed and does not show the short-term oscillations indicated by Fig. 5. The growth rate is given in units of  $\Omega (r_{\rm max})$.
Open with DEXTER

In the standard disk model with positive  $\alpha _{\rm SS}$, the gas moves to the central object in the surface layers of the disk and in the opposite direction near the midplane (Urpin 1984; Kley & Lin 1992; Lee & Ramirez-Ruiz 2002). Because of the closed radial boundaries, we cannot see the details of accretion in our simulations. An improved model will likely show matter flowing inward near the central plane and outward in the surface layers.

The energy of the fluctuations reaches values of the order of $1\%$ of the thermal energy. Initially, the energy of fluctuations grows fairly fast with a growth rate of $\sim$ $0.03~\Omega(r_{\rm max})$ measured in units of the angular velocity at the outer edge of the computational domain. This corresponds to nearly 0.2 per orbital period. After $\sim$100 outer-disk revolutions, however, the growth decreases by an order of magnitude and indicates saturation of the fluctuation strength. The growth rate evolution is shown in Fig. 7 in terms of $\Omega$ at the outer edge of the computational domain. Note that the growth rate of velocity fluctuations, $\gamma$, introduced in Sect. 2 is smaller by a factor of 1/2 since the energy of turbulence $\propto$u²; hence, $\gamma \sim 0.015~\Omega(r_{\rm max})$.

Initially, the energy of azimuthal fluctuations grows most rapidly, followed by the radial ones. The vertical fluctuations show slowest growth up to about 50 outer orbits and exhibit strong variations in their growth rate. The evolution is quasi-steady with average growth rates returning to nearly zero after 105 outer orbits. In the quasi-steady state, the kinetic energy contained in radial fluctuations is maximal, azimuthal fluctuations are slightly smaller, and vertical fluctuations are by a factor of $\sim$5 less efficient.

	Figure 8: Energies in velocity fluctuations for model L.
Open with DEXTER

	Figure 9: The evolution of the $\alpha _{\rm SS}$-parameter for model L.
Open with DEXTER

The dynamics of the model L shown in Figs. 8-10 is qualitatively the same. The model applies the initial conditions (24) to (26). Note that the instability requires the longest time to reach saturation in the region near the outer boundary where $\Omega$ is minimal. As a result, all turbulent characteristics integrated (or averaged) over the computational domain such as kinetic energy, $\alpha _{\rm SS}$, etc., can become quasi-steady only when turbulence reaches saturation in the outermost region. Therefore, the dimensionless timescale of saturation is approximately the same for both the S and L models if measured in units of the outer rotational period (upper abscissa in the plots). In these units, the saturation timescale is about 100 Keplerian periods. However, measured in units of the inner period, the values of this timescale is by the factors  $(6/4)^{3/2} \approx 1.8$ and $(7/3)^{3/2} \approx 3.5$ longer for the models S and L, respectively.

Unlike in model S, the saturation kinetic energy in the model L reaches less than $\sim$$0.1 \%$ of the thermal energy. The reason is probably the longer physical time over which we integrated here. Numerical viscosity had thus more time to damp the fluctuations before they reach their maximum amplitude. Another issue may be the different sizes of computational domains; since the energy of fluctuations is certainly not distributed homogeneously, runs S and L are likely to exhibit different saturation energies on average.

Since model L is more extended in the radial direction than model S, the ratio of kinetic energy of turbulence in the quasi-steady state and the thermal energy depends weakly on r. The energy is distributed non-equally between the velocity components (see Fig. 8) with radial and azimuthal fluctuations containing the main fraction of kinetic energy and vertical fluctuations being less energetic, again by a factor of $\sim$5. The value of  $\alpha _{\rm SS}$ grows to small positive numbers with their maximum at 6 $\times$ 10^-6, and decays but remains positive throughout the rest of the simulation.

	Figure 10: The growth rate for the square of the velocity fluctuations for model L in units of  $\Omega (r_{\rm max})$.
Open with DEXTER

Contrasting with model S, a transient energy peak is observed in both the unperturbed simulation (Fig. 3) and the perturbed run (Fig. 8) of model L. The peak is attributed to reorganization motions caused by the deviations of the initial conditions from the correct hydrostatic solution. The deviations are naturally larger in model L with its larger radial extent, but they do not come along with signatures in the $\alpha _{\rm SS}$ plots.

In Figs. 11, 13 and 15, we plot the spectral characteristics of the flow generated by the vertical shear instability for the model S. For a comparison, the Kolmogorov spectrum with a power of -5/3 is shown by a dotted line in all figures. The Fourier amplitude for any given coordinate is shown as an average over the two other coordinates. Let the direction of decomposition be r, then we obtain $N_z\times N_\varphi$ individual radial spectra. The absolute values of these spectra are averaged; the square of these average amplitudes is shown in the graphs. The spectra shown in Figs. 11, 13 and 15 are taken from a phase of the simulation which is about 20 inner orbits before the saturation phase. Spectra of model L taken from a similar phase of evolution are plotted in Figs. 12, 14 and 16.

The spectra decomposed in the radial and vertical directions reach a quasi-steady regime on a shorter timescale than other turbulent characteristics like the kinetic energy or  $\alpha _{\rm SS}$. These spectra do not change their shape substantially after $t \sim 60$ orbital periods at the inner radius. They do change their absolute values though.

	Figure 11: Spectra of the velocity components of model S with initial perturbations, decomposed in radial direction after 164 inner periods.
Open with DEXTER

	Figure 12: Spectra of the velocity components of model L with initial perturbations, decomposed in radial direction after 223 inner periods.
Open with DEXTER

	Figure 13: Spectra of the velocity components of model S with initial perturbations, decomposed in vertical direction after 164 inner periods.
Open with DEXTER

	Figure 14: Spectra of the velocity components of model L with initial perturbations, decomposed in vertical direction after 223 inner periods.
Open with DEXTER

	Figure 15: Spectra of the velocity components of model S with initial perturbations, decomposed in azimuthal direction after 164 inner periods.
Open with DEXTER

	Figure 16: Spectra of the velocity components of model L with initial perturbations, decomposed in azimuthal direction after 223 inner periods.
Open with DEXTER

The spectra of the vertical and azimuthal velocities decomposed in the radial direction (Figs. 11 and 12) are particularly simple: they are monotonous with a slope close to 5/3 and do not exhibit any noticeable maximum. This implies that hydrodynamic non-linearity rather than other factors governs the radial dependence of the z- and $\varphi$-components of velocity fluctuations. This also holds for the spectra in the quasi-steady regime of the saturated phase of the simulation. Only the spectra of the radial velocity have steeper slopes than the Kolmogorov spectrum. The large gradient is likely caused by the strong differential rotation in the radial direction that is an additional powerful mechanism of distortion of the radial motions in accretion discs.

Spectra decomposed in the vertical direction show a qualitatively similar behavior (Figs. 13 and 14). In this case, spectra of the radial and azimuthal velocity are steeper than a Kolmogorov one and have exponents of $\sim$-3in model S and $\sim$-4 in model L.

Obviously, the dependence of the velocity on $\varphi$ experiences the most dramatic changes from initial disturbances being nearly independent of $\varphi$ to a rather complicated dependence with a number of peaks shown in Figs. 15 and 16. These peaks coincide with even mode numbers m. The $\varphi$-amplitudes of all three velocity components are qualitatively similar and exhibit a maximum at the azimuthal mode number $m \sim 6$-8 that corresponds to a wavelength comparable to the disk radius. The spectra decomposed in the $\varphi$-direction are qualitatively different from the Kolmogorov spectrum: the decrease of power towards larger wave numbers is the steepest among the spectra obtained. For example, fluctuations with the azimuthal wavelength $\sim$r are approximately 10⁴-10⁵ times stronger than fluctuations with the wavelength $\sim$H. The corresponding exponent is $\sim$-5.

This rapid decrease of the spectra towards high k might be explained by a stabilizing influence of the differential rotation on the fluctuations with large azimuthal wavenumbers. The stabilization prevents a non-linear generation of such modes. Therefore, turbulence in the $\varphi$-direction is mainly represented by relatively large structures with the length-scale comparable to ror longer. The vertical shear instability apparently generates anisotropic turbulence in accretion discs.

	Figure 17: Dependence of radial spectra on height above and below the equatorial plane, taken from model L after  $349\tau _{\rm orb}$. The top graph is based on ur, the bottom graph on uz.
Open with DEXTER

Figure 17 shows the vertical distribution of radial spectra. The contour lines represent the power of spectra taken at each vertical distance from the equator z. The upper panel is based on the radial velocity component u_r, the lower, on the vertical component u_z. Both graphs show that high-k_r modes have more power at large equator distances. In other words, small scales develop best in the disk halo. In the equatorial plane, low-k_r modes dominate in u_r; because of the nearly symmetric flow, the spectrum of the equatorial u_z is close to zero there.

4.3 Comparison with sound waves

Note that in non-magnetic discs, there can generally exist one more instability apart from that considered in the present paper (see Papaloizou & Pringle 1984; Drury 1985). The instability arises by over-reflection of sound waves at a corotation radius, combined with a reflection at (at least) one boundary.

A number of points discriminate between the two instabilities. For instance, the growth rate is approximately an order of magnitude larger for the vertical shear instability ($\approx$ $0.02 \Omega$) than for the reflection one. The growth rate in our simulations agrees very well with the value predicted by the analytic theory of the vertical shear instability; cf. Eq. (19). The reflection instability is adiabatic whereas the instability considered in our paper is substantially non-adiabatic. A numerical evaluation of the same model in the adiabatic limit done by Rüdiger et al. (2002) failed to indicate any sign of instability.

	Figure 18: Energies in velocity fluctuations for model O.
Open with DEXTER

Since the reflection instability is caused by reflection of sound waves at the boundaries, the growth rate of the reflection instability depends critically on the separation between the boundaries which is crucial for the formation of an interference pattern of maxima and minima. By contrast, the vertical shear instability is not sensitive to the positions of the boundaries as shown by the comparison of models S and L.

Outflow boundaries in radial direction shed further light on the nature of the instability, since the sound wave instability should be substantially suppressed in the disk with these. The energies of fluctuations of Model O are shown in Fig. 18. The onset of the instability is still visible, but the model loses too much energy through the radial boundaries before saturation can be observed.

We also made a reverse test of whether or not the vertical-shear instability is actually excited in models S, L, and O. The vertical structure of the disk is switched off by cancelling the vertical component of the gravitational force of the central mass. The equilibrium model then shows no density stratification and no vertical shear. If the same initial perturbations as for Model L (and the same numerical parameters) are employed, no instability emerges as shown in Fig. 19. The level of kinetic energy in the radial and azimuthal fluctuations is 10^-8- $10^{-7} c_{\rm s}^2$. The energy in vertical fluctuations is damped exponentially to zero. Since the weak initial radial fluctuations show little change over the time of simulation. This is regarded as the contribution of sound waves to our simulations, which is small compared to the growth of fluctuations in the models of Sect. 4.2 which have vertical shear.

	Figure 19: Energies in velocity fluctuations for a perturbed model with the same parameters as L, but without vertical gravity (cylindrical potential).
Open with DEXTER

The resulting hydrodynamic pictures associated to the two instabilities are very much different as well. The reflection instability leads to global, large-scale 2D patterns in the $r-\varphi$ plane with small scales developing only at the boundaries, and only in the non-linear regime. The simulations by Kaisig (1989) show no evidence for the development of turbulence. The characteristic length-scales in the radial and azimuthal direction are comparable and $\sim$H. This differs completely from the picture produced by the vertical shear instability generating predominantly anisotropic fluctuations with the largest length-scale ($\sim$r) in the azimuthal direction and the shortest length-scale ($\ll$H) in the radial one, in agreement with the predictions of Sect. 2.

All this allows us to conclude that the turbulence simulated in this paper is caused by the vertical shear instability.

5 Summary

We have presented the results of simulations of the vertical shear instability in accretion disks. The instability can arise for any dependence of the angular velocity on the vertical coordinate and, hence, accretion disks are subject to this instability since $\Omega$ depends on z. As it was first shown by Goldreich & Schubert (1967), this instability is associated with the heat transport and occurs only for non-adiabatic disturbances. Therefore, we concentrated on the behavior of short wavelength disturbances which are certainly non-adiabatic because of the high thermal conductivity of disks. The important point of our numerical simulations is a high radial resolution because only disturbances with the wavelength much shorter in the radial direction than vertically can be unstable (Urpin 2003). Our simulations follow both the linear and non-linear evolution of short wavelength disturbances and indicate the presence of instability. The main findings of the present study are summarized as follows:

1)

The vertical shear instability grows on a relatively short time scale. The growth rate $\gamma$ reaches $\sim$ $0.015 \Omega(r_{\rm max})$. This value is in a good agreement with the growth rate predicted by the linear theory, $\gamma \sim \Omega (qH/r) \sim 0.01 \Omega$, and is independent of the separation of the radial boundaries. After about 100 revolutions of the outer boundary of the computational domain, the instability operates in a quasi-steady regime, and the growth rate decreases to nearly zero or - in the very long run of model L - to a tiny negative value which is due to the final viscous decay of the fluctuations.

2)

The energy contained in turbulent fluctuations can reach values of the order of 1% of the thermal energy. The energy is distributed non-equally among the velocity components, and turbulence is substantially anisotropic. The radial fluctuations are most energetic and contain about 3-4 times more energy than the vertical fluctuations. The azimuthal velocity component contains similar or slightly less energy than the radial one.

3)

The Fourier spectra of velocity fluctuations decomposed in the azimuthal direction tend to exhibit non-monotonous behavior which can be interpreted as the presence of preferable length-scales in the flow. Highest power is found at wavelengths of $\sim$r for all three velocity components. Spectra decomposed in the radial and vertical directions are monotonous and do not exhibit such maximum. The turbulence generated by the vertical shear instability is anisotropic with the largest length-scale in the azimuthal direction and the shortest length-scale in the radial one.

4)

The flow caused by the vertical shear instability can enhance angular momentum transport to a limited degree. The turbulent  $\alpha _{\rm SS}$-parameter reaches the value of $\sim$-0.0005 at the end of simulation S and peaks at 6 $\times$ 10^-6 in model L. The relation between  $\alpha _{\rm SS}$ and the direction of the flow may not be realistic in a model with closed radial boundaries, and we are not drawing conclusions about radial flows. Further simulations with accreting boundaries may provide more information on the accretion structure of disks being subject to the vertical shear instability.

Acknowledgements

The authors thank Günther Rüdiger for stimulating discussions. One of the authors (V.U.) thanks the Deutsche Forschungsgemeinschaft for the support (grant DFG 436 RUS 113/559) and the Astrophysikalisches Institut Potsdam for hospitality.

References

• Acheson, D. J. 1978, Phil. Trans. Roy. Soc. Lond., A289, 459 [NASA ADS] (In the text)
• Acheson, D. J. 1979, Sol. Phys., 62, 23 [NASA ADS] (In the text)
• Arlt, R., & Rüdiger, G. 2001, A&A, 374, 1035 [EDP Sciences] [NASA ADS] (In the text)
• Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] (In the text)
• Balbus, S. A., & Hawley, J. F. 1998, Rev. Mod. Phys., 70, 1 [NASA ADS] [CrossRef] (In the text)
• Brandenburg, A., Nordlund, A., Stein, R., & Torkelsson, U. 1995, ApJ, 446, 741 [NASA ADS] [CrossRef] (In the text)
• Drury, L. 1985, MNRAS, 217, 821 [NASA ADS] (In the text)
• Dubrulle, B., Marie, L., Normand, C., et al. 2003, A&A, submitted
• Fricke, K. 1969, A&A, 1, 388 [NASA ADS] (In the text)
• Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [NASA ADS] [CrossRef] (In the text)
• Kaisig, M. 1989, A&A, 218, 89 [NASA ADS] (In the text)
• Kippenhahn, R., & Thomas, H.-C. 1982, A&A, 114, 77 [NASA ADS] (In the text)
• Klahr, H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef] (In the text)
• Kley, W., & Lin, D. N. C. 1992, ApJ, 397, 600 [NASA ADS] [CrossRef] (In the text)
• Kurzweg, U. H. 1963, J. Fluid Mech., 17, 52 [NASA ADS] (In the text)
• Lee, W., & Ramirez-Ruiz, E. 2002, ApJ, 577, 893 [NASA ADS] [CrossRef] (In the text)
• Longaretti, P.-Y. 2002, ApJ, 576, 587 [NASA ADS] [CrossRef] (In the text)
• Molemaker, M. J., McWilliams, J. C., & Yavneh, I. 2001, Phys. Rev. Lett., 86, 5270 [NASA ADS] [CrossRef] (In the text)
• Papaloizou, J., & Pringle, J. 1984, MNRAS, 208, 721 [NASA ADS] (In the text)
• Richard, D., & Zahn, J.-P. 1999, A&A, 347, 734 [NASA ADS] (In the text)
• Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Safronov, V. S. 1969, Evolution of the protoplanetary cloud and formation of the Earth (Moscow: Nauka) (In the text)
• Schultz-Grunow, F. 1959, Z. Angew. Math. Mech., 39, 101 [NASA ADS]
• Stone, J., & Norman, M. 1992a, ApJS, 80, 753 [NASA ADS] [CrossRef] (In the text)
• Stone, J., & Norman, M. 1992b, ApJS, 80, 791 [NASA ADS] [CrossRef]
• Stone, J., Mihalas, D., & Norman, M. 1992, ApJS, 80, 819 [NASA ADS] [CrossRef] (In the text)
• Torkelsson, U., Brandenburg, A., Nordlund, A., & Stein, R. 1997, IAU Coll., 163, ASP Conf. Ser., 121, 210 (In the text)
• Urpin, V. 1984, SvA, 28, 50 [NASA ADS] (In the text)
• Urpin, V. 2003, A&A, 404, 397 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399 [NASA ADS] [CrossRef] (In the text)
• Velikhov, E. P. 1959, Sov. Phys. JETP, 9, 995 [NASA ADS] [MathSciNet] (In the text)