A&A 391, 781-787 (2002)
DOI: 10.1051/0004-6361:20020853

Hydrodynamic stability in accretion disks under the combined influence of shear and density stratification

G. Rüdiger1 - R. Arlt1 - D. Shalybkov1,2


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 February 2002 / Accepted 4 June 2002

Abstract
The hydrodynamic stability of accretion disks is considered. The particular question is whether the combined action of a (stable) vertical density stratification and a (stable) radial differential rotation gives rise to a new instability for nonaxisymmetric modes of disturbances. The existence of such an instability is not suggested by the well-known Solberg-Høiland criterion. It is also not suggested by a local analysis for disturbances in general stratifications of entropy and angular momentum which is presented in our Sect. 2. This confirms the results of the Solberg-Høiland criterion also for nonaxisymmetric modes within the frame of ideal hydrodynamics but only in the frame of a short-wave approximation for small m. As a necessary condition for stability we find that only conservative external forces are allowed to influence the stable disk. As magnetic forces are never conservative, linear disk instabilities should only exist in the magnetohydrodynamical regime which indeed contains the magnetorotational instability as a much-promising candidate.

To overcome some of the used approximations in a numerical approach, the equations of the compressible adiabatic hydrodynamics are integrated, imposing initial nonaxisymmetric velocity perturbations with m=1 to m=200. Only solutions with decaying kinetic energy are found. The system always settles in a vertical equilibrium stratification according to pressure balance with the gravitational potential of the central object.

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


1 Introduction: The Solberg-Høiland criterion

According to the Rayleigh criterion, Keplerian disks are hydrodynamically stable configurations, because the angular momentum per unit mass, $j=R^2 {\it\Omega}$, increases with radius, i.e.
 
$\displaystyle {{\rm d}j^2 \over {\rm d}R} > 0$     (1)

with R being the distance from the rotation axis. They are unstable configurations, however, under the influence of a (weak) magnetic field if the disk plasma has a sufficiently high conductivity (Balbus & Hawley 1991; Balbus 1995; Brandenburg et al. 1995; Ziegler & Rüdiger 2000). There remains the question whether also other, nonmagnetic influences exist which in combination with the basic (negative) shear ${\rm d}{\it\Omega}/{\rm d}R$ allow instabilities in the linear regime. We recall three recent discussions. Urpin & Brandenburg (1998) consider the destabilizing action of any vertical dependence of ${\it\Omega}$ on z which may more than compensate the stabilization of the positive radial gradient of j. There is also the suggestion of Klahr & Bodenheimer (2001) that the negative radial entropy stratification in thin Keplerian disks may act in the sense of a destabilization. Note that with the standard-alpha description of accretion disks one finds the positive value
 
$\displaystyle {\partial S \over \partial R} \simeq {1\over 2} {C_{\rm v} \over R},$     (2)

indicating, however, stability. Molemaker et al. (2001) and Yavnesh et al. (2001) recently pointed out that a vertical density stratification in the Taylor-Couette flow may have a similar destabilizing effect as a global vertical magnetic field for an electrically conducting fluid. The sufficient condition that a (nonaxisymmetric) disturbance will be unstable is given to be the same as for an unstratified ideal MHD Taylor-Couette flow in the presence of a vertical magnetic field, i.e.
 
$\displaystyle {{\rm d} {\it\Omega}\over {\rm d}R} <0.$     (3)

If this were also true for nonmagnetic accretion disks - which always have vertical density stratifications - they would be hydrodynamically unstable, so that a turbulence can develop transporting the angular momentum outwards.

There is indeed some analogy between accretion disks and experimental Taylor-Couette flows. But the question remains why such a hydrodynamic "radial-shear+vertical-stratification'' instability might exist in accretion disks even though the extensive numerical simulations in local shearing-box systems have always shown that Keplerian disks[*] are stable to infinitesimal and finite disturbances (Balbus et al. 1996). Most of the mentioned discussions are covered by the Solberg-Høiland conditions for dynamic stability which, for the geometrical constellation studied here, take the form

 
$\displaystyle {1\over R^3} {\partial j^2 \over \partial R} - {1\over C_p \rho} \nabla p \cdot \nabla S > 0,$     (4)


 
$\displaystyle {\partial p \over \partial z} \left({\partial j^2 \over \partial ... ...l z} - {\partial j^2 \over \partial z} {\partial S \over \partial R}\right) < 0$     (5)

(Tassoul 2000). If $\partial p/\partial z < 0$ (as is usual for accretion disks), the latter relation turns into
 
$\displaystyle {\partial j^2 \over \partial R} {\partial S \over \partial z} > {\partial j^2 \over \partial z} {\partial S \over \partial R},$     (6)

which for Keplerian disks with ${\it\Omega}= {\it\Omega}(R) \propto R^{-1.5}$ simply leads to the usual Schwarzschild criterion,
 
$\displaystyle {\partial S \over \partial z} > 0,$     (7)

for stability. However, if the angular velocity ${\it\Omega}$ depends on the vertical coordinate z then Eq. (6) allows instability in the standard case of positive $\partial j/\partial R$ and $\partial S/ \partial z$ in case they are large enough. Equation (4) becomes
 
$\displaystyle {1\over R^3} {\partial j^2 \over \partial R} + {g\over C_p} {\par... ...over \rho C_p} {\partial p \over \partial R} {\partial S \over \partial R} > 0,$     (8)

which with (7) and for ${\rm d}S/{\rm d}R > 0$ ensures stability. Here we have used $-g_z = g= {\it\Omega}^2 z > 0$. In the opposite case of a negative entropy gradient one at least needs
 
$\displaystyle \left\vert{\partial S\over \partial R}\right\vert > {\rho C_{\rm p} {\it\Omega}^2 \over \vert\partial p/\partial R\vert}$     (9)

in order to violate Eq. (8) in the most optimistic case of vertical isentropy. If, on the other hand, the standard stratifications of accretion disks are j = j(R) and S=S(z), with ${\rm d}j^2/{\rm d}R > 0$ then the Solberg-Høiland criterion for stability reduces to
 
$\displaystyle 1+ {z\over C_p} {{\rm d}S \over {\rm d}z} > 0, \ \ \ \ \ \ \ \ \ {{\rm d}S \over {\rm d}z} > 0,$     (10)

(in the upper hemisphere). Obviously ${\rm d}S/{\rm d}z > 0$ is a sufficient condition for stability of Kepler disks after the Solberg-Høiland criterion also for the combined action of vertical density stratification and differential Kepler rotation. The same result have been derived in earlier papers by Livio & Shaviv (1977), Abramowicz et al. (1984) and Elstner et al. (1989). The local analysis for ring-like disturbances by Elstner et al. demonstrates in all details that in the framework of an ideal hydrodynamics, i.e. for vanishing Shakura-Sunyaev alpha, the traditional Schwarzschild criterion ${\rm d}S/{\rm d}z > 0$ ensures stability for the Kepler disk. If the disk can be considered as isothermal in the vertical direction then it follows
 
$\displaystyle {{\rm d}S \over {\rm d}z} = - {\cal R\over \mu} {{\rm d} \log \rho \over {\rm d}z}$     (11)

which is always positive for the typical density stratification, i.e. the disk should always be stable.

2 Combined stability conditions

In this section a local stability analysis of the equations of ideal hydrodynamics in cylindrical coordinates ($R,\phi,z$) is presented. Self-gravitation phenomena are excluded. Sound waves and nonaxisymmetric disturbances, however, are in particular allowed to exist. The three components of the momentum equation are
   
$\displaystyle {\partial u_R \over \partial t} + u_R {\partial u_R \over \partia... ...} - { u_{\phi}^2 \over R} = -{1 \over \rho}{\partial p \over \partial R} + g_R,$     (12)
$\displaystyle {\partial u_\phi \over \partial t} + u_R {\partial u_\phi \over \... ...i} u_R \over R} = -{1 \over \rho R}{\partial p \over \partial \phi} + g_{\phi},$     (13)
$\displaystyle {\partial u_z \over \partial t} + u_R {\partial u_z \over \partia... ...ial u_z \over \partial z} = -{1 \over \rho}{\partial p \over \partial z} + g_z.$     (14)

The mass conservation reads
 
$\displaystyle {\partial \rho \over \partial t} + {\partial \over \partial R}(\r... ... \over \partial \phi}(\rho u_\phi) + {\partial \over \partial z} (\rho u_z) = 0$     (15)

and the energy equation is
 
$\displaystyle {\partial S \over \partial t} + u_R {\partial S \over \partial R}... ...er R} {\partial S \over \partial \phi} + u_z {\partial S \over \partial z} = 0.$     (16)

As usual $\rho$ is the density, p is the pressure, S is the specific entropy, $\mbox{\boldmath$g$ }$ is the vector of any acceleration. The unperturbed state is assumed to be uR,0=uz,0=0, $u_{\phi,0}=R {\it\Omega}(R,z)$, $\rho_0=\rho_0(R,z)$, p0=p0(R,z), S0=S0(R,z), gR=gR(R,z), $g_\phi=0$, gz=gz(R,z). The only non-zero equations of the system (12)-(16) are Eqs. (12) and (14) which take the form
 
$\displaystyle g_R={1\over \rho_0} {\partial p_0\over \partial R}-R{\it\Omega}^2, \qquad g_z={1\over \rho_0} {\partial p_0\over \partial z}\cdot$     (17)

Finally, the ideal gas equation of state is used $p={\cal R}/\mu \rho T$, ${\rm d}S=C_{\rm v} {\rm d}\log (p \rho^{-\gamma})$ with $C_{\rm v}$ as the specific heat at constant volume, and $\gamma=C_p/C_{\rm v}=5/3$.

2.1 Linearized equations

The perturbations in the force are neglected. The perturbations uR, $u_{\phi}, u_z, \rho_1, p_1$ to the basic state are assumed small and the linearized system (12)-(16) is
     
$\displaystyle {\partial u_R \over \partial t} + {\it\Omega}{\partial u_R \over ... ..._1 \over \partial R} -{\rho_1 \over \rho_0^2}{\partial p_0 \over \partial R}=0,$     (18)
$\displaystyle {\partial u_\phi \over \partial t} + {1 \over R} {\partial \over ... ...}\over \partial z} u_z +{1 \over \rho_0 R}{\partial p_1 \over \partial \phi}=0,$     (19)
$\displaystyle {\partial u_z \over \partial t} + {\it\Omega}{\partial u_z \over ... ...1 \over \partial z} - {\rho_1 \over \rho_0^2}{\partial p_0 \over \partial z}=0,$     (20)
$\displaystyle {1 \over \rho_0}{\partial \rho_1 \over \partial t} + {\partial u_... ...og\rho_0 \over \partial R} u_R + {\partial \log\rho_0 \over \partial z} u_z =0,$     (21)
$\displaystyle {\partial S_1 \over \partial t} + {\it\Omega}{\partial S_1 \over ... ...i} + u_R{\partial S_0 \over \partial R} + u_z{\partial S_0 \over \partial z}=0.$     (22)

In the frame of our local approximation the equations are considered only in small volume near some reference point ( $R_0,\phi_0,z_0$). The coefficients in the equations are thus constant taken at ( $R_0,\phi_0,z_0$). As usual both the short-wave approximation and the small-m approximation, i.e.
 
$\displaystyle \vert k_RR\vert \gg m, \ \ \ \ \ \vert k_zz\vert \gg m$     (23)

are applied (cf. Meinel 1983). In this approximation all perturbations can be expressed by the modal representation ${\rm exp}[{\rm i}(k_RR+m\phi+k_zz-\omega t)]$ and the equations take the final form
     
$\displaystyle -{\rm i}(\omega-m{\it\Omega})u_R - 2{\it\Omega}u_\phi +{\rm i}k_R{p_1 \over \rho_0} -{\rho_1 \over \rho_0^2}{\partial p_0 \over \partial R}=0,$     (24)
$\displaystyle -{\rm i}(\omega-m{\it\Omega})u_\phi + {1 \over R}{\partial \over \partial R}(R^2 {\it\Omega}) u_R+R {\partial {\it\Omega}\over \partial z} u_z=0,$     (25)
$\displaystyle -{\rm i}(\omega-m{\it\Omega})u_z +{\rm i}k_z{p_1 \over \rho_0} -{\rho_1 \over \rho_0^2}{\partial p_0 \over \partial z}=0,$     (26)
$\displaystyle -{\rm i}(\omega-m{\it\Omega}) {\rho_1 \over \rho_0}+ {\partial \l... ...}u_R+ {\partial \log\rho_0 \over \partial z}u_z+ {\rm i}k_Ru_R+{\rm i}k_zu_z=0,$      
      (27)
$\displaystyle -{\rm i}(\omega-m{\it\Omega}) \left( {p_1 \over p_0} - \gamma {\r... ...tial S_0 \over \partial R} u_R+ {\partial S_0 \over \partial z} u_z \right) =0.$     (28)

Generally, it is
 
$\displaystyle {1 \over C_{\rm v} }{\partial S_0 \over \partial R} = {\partial \... ... {\partial \over \partial z} \log \left( {p_0 \over \rho_0^\gamma} \right)\cdot$     (29)

The epicyclic frequency $\kappa$ and the adiabatic sound speed $c_{\rm ac}$ are
 
$\displaystyle \kappa^2 = {1\over R^3} {\partial j^2\over \partial R}, \qquad\qquad\qquad c_{\rm ac}^2 = \gamma {p_0 \over \rho_0}\cdot$     (30)

The determinant of the above homogeneous fifth order system must vanish and the resulting dispersion relation is
 
$\displaystyle (\omega- m{\it\Omega}) \cdot D =0$     (31)

with
 
$\displaystyle D=(\omega-m{\it\Omega})^4-2 E (\omega-m{\it\Omega})^2 + F.$     (32)

Here we have introduced the expressions
 
$\displaystyle E={1\over2} \left[ (k_R^2+k_z^2)c_{\rm ac}^2 +{1 \over \rho_0^2}{... ...ial p_0 \over \partial z} {\partial \rho_0 \over \partial z} + \kappa^2 \right]$     (33)

and

 
                        $\displaystyle F= \left( {k_R \over \rho_0} {\partial p_0 \over \partial z} - {k... ...partial \over \partial R} \log \left( {p_0 \over \rho_0^\gamma} \right) \right.$      
$\displaystyle \left. +{{\rm i} \over \rho_0^2} \left({\partial p_0 \over \parti... ...al z} - {k_z \over \rho_0} {\partial p_0 \over \partial R} \right) \right)\cdot$     (34)

For axisymmetric modes we simply have to replace $\omega-m{\it\Omega}$ by $\omega$. This procedure influences the frequencies rather than the stability. The stability of axisymmetric and nonaxisymmetric perturbations is thus defined by the same stability criterion in the approximation used. The first factor in (31) describes the waves with $\omega=m{\it\Omega}$ and will not be considered below. The stability is defined by the factor D. The roots of D are
 
$\displaystyle (\omega-m{\it\Omega})^2= E \pm \sqrt{E^2-F}.$     (35)

According to (33) the coefficient E is real. The flow is thus always unstable for negative E and complex F. Positive E and real F are, therefore, the necessary conditions for stability. According to (34), F is real if and only if
 
$\displaystyle {1 \over \rho_0^2} \left( {\partial p_0 \over \partial z} {\parti... ..._0 \over \partial z} \right) - R {\partial {\it\Omega}^2 \over \partial z} = 0.$     (36)

Inserting (36) in (17), we find
 
$\displaystyle {\partial g_R \over \partial z} = {\partial g_z \over \partial R}$     (37)

as the immediate consequence. Without forces this relation is always fulfilled. Any conservative force is a particular solution of (37). If - as it is in accretion disks - the gravity balances the pressure and the centrifugal force, then Eq. (36) is automatically fulfilled. Note that after the Poincare theorem for rotating media with potential force and ${\it\Omega}={\it\Omega}(R)$ both the density and the pressure can be written as functions of the generalized potential so that (36) is always fulfilled. Generally, the magnetic field is not conservative and can never fulfill the condition (37). This is the basic explanation for the existence of the magnetorotational instability which is driven by (weak) magnetic fields. That instability is typically obtained in a short-wavelength approximation too, but is largely independent of the assumptions made.

2.2 Sufficient conditions for stability

Now we can suppose that the necessary condition (36) is fulfilled. The flow is stable if both roots in (35) are real and positive. The sufficient conditions for stability are therefore
 
$\displaystyle E>0, \qquad\qquad\qquad0<F<E^2.$     (38)

Inserting (36) into (34) we find F to be positive if
 
$\displaystyle {k_R^2 \over k_z^2} N_z^2 + {k_R \over k_z} {2 \over \gamma \rho_... ...l {\it\Omega}^2 \over \partial z} {\partial \rho_0 \over \partial R} \right)>0,$     (39)

where
 
$\displaystyle N_z^2=-{1 \over \gamma \rho_0} {\partial p_0 \over \partial z} {\partial \over \partial z} \log \left( {p_0 \over \rho_0^\gamma} \right)$     (40)

and
 
$\displaystyle N_R^2=-{1 \over \gamma \rho_0} {\partial p_0 \over \partial R} {\partial \over \partial R} \log \left( {p_0 \over \rho_0^\gamma} \right)$     (41)

are components of the Brunt-Väisälä frequency. In our short-wave approximation (see (23)) we can neglect the last term on the left hand side of (39). The expression (39) is thus a simple quadratic expression in kR / kz. Two conditions will ensure its positivity: i) the expression is positive for some value kR / kz and ii) the expression has no real roots i.e. the discriminant is negative. The first condition is fulfilled if
 
$\displaystyle N_R^2+N_z^2+\kappa^2>0.$     (42)

This is the Schwarzschild criterion for stability in the presence of rotation. The second condition leads to
 
$\displaystyle {\partial p_0 \over \partial z} \left( \kappa^2 {\partial \over \... ...tial \over \partial R} \log \left( {p_0 \over \rho_0^\gamma} \right) \right)<0.$     (43)

Equations (42) and (43) are exactly equivalent to the Solberg-Høiland conditions (4) and (5). The Schwarzschild criterion (42) ensures that E>0 and the short-wave approximation ensures that F<E2 and these relations do not yield any additional conditions. The Solberg-Høiland conditions (4) and (5) can also be obtained in the Boussinesq approximation. Goldreich & Schubert (1967) have already formulated the dispersion relation within the framework of the Boussinesq approximation. Nevertheless, our more general consideration allows to find the new necessary condition for stability (36).

2.3 Limiting cases

Let us consider two interesting limiting cases. Without rotation the Schwarzschild criterion (4) takes the form
 
NR2+Nz2>0,     (44)

where NR and Nz are given by (40) and (41). The stratification is thus always stable if p0 and S0 have opposite stratifications - which indeed is the standard case. Nevertheless, Eq. (44) allows the stability also for the case when p0 and S0 have the same stratifications in one direction and opposite stratifications in another. Note that according to (44) the situation with pressure stratification without any density stratification is always unstable. Without density stratification the sufficient stability condition (4) for a rotating flow is
 
$\displaystyle \kappa^2> {1 \over \rho_0^2 c_{\rm ac}^2} \left( {\partial p_0 \o... ...er \rho_0^2 c_{\rm ac}^2} \left( {\partial p_0 \over \partial z} \right)^2\cdot$     (45)

This is the classical Rayleigh criterion (1) for stability generalized to the compressible case. Obviously, the right-hand side is positive and the criterion (45) is more restrictive than the standard one.

3 Numerical simulations

The main shortcoming in our presentation is the short-wave assumption (23) under which all the above considerations are only valid. Without this assumption, Fourier modes are no longer solutions of the linearized equations and we have to switch to much more complicated mathematics involving integral equations (see Rüdiger & Kitchatinov 2000). We prefer to present numerical simulations for the stability of (isothermal) density-stratified Keplerian disks under the influence of finite disturbances. Along this way also finite-amplitude disturbances can be considered within the framework of a nonideal hydrodynamics. The above analytical study has shown that stratified disks are stable; under the assumption of a small scale in perturbations compared with the disk dimensions. We present a few numerical integrations of the hydrodynamic equations for the unrestricted case. The computations are based on the work of Arlt & Rüdiger (2001) which required a stable, stratified Keplerian disk as a prerequisite. We adapt these simulations with adiabatic evolution and relatively strong nonaxisymmetric perturbations. The setup for our three-dimensional simulations applies a global, cylindrical computational domain with full azimuthal range, dimensionless radii from 4 to 6, and a vertical extent of 2 density scale heights on average, which is -1 to +1 in dimensionless units. The ZEUS-3D code (developed by Stone & Norman 1992a,b; Stone et al. 1992) was used for the integration of the hydrodynamic equations. The modifications to the code are very close to those given in Arlt & Rüdiger (2001); the induction equation is naturally not integrated here.
 
$\displaystyle \frac{\partial\rho}{\partial t}+{\rm div}(\rho {\mbox{\boldmath$u$ }})=0,$     (46)


 
$\displaystyle \rho \frac{\partial{\mbox{\boldmath$u$ }}}{\partial t}+ \rho({\mb... ...$u$ }}\nabla) {\mbox{\boldmath$u$ }} = - \nabla p + \rho~ \mbox{\boldmath$g$ },$     (47)


 
$\displaystyle \frac{\partial e}{\partial t}+{\rm div}(e{\mbox{\boldmath$u$ }})=-p~{\rm div}{\mbox{\boldmath$u$ }},$     (48)

where ${\mbox{\boldmath$u$ }}$, $\rho$, e, and p have the usual meanings of velocity, density, energy density per unit volume, and pressure. The source of gravitation is that of a point mass M=105 at R=0. The simulations are performed in a non-rotating reference frame; they do not include self-gravity in the disk. The polytrophic exponent is $\gamma=5/3$. The original ZEUS code provides artificial viscosities for improved shock evolution; we have put these terms to zero here. The advection interpolation is the second order van-Leer scheme. The conditions for the vertical boundaries are stress-free; no matter can exit the computational domain in vertical direction and the vertical derivatives of uR and $u_\phi$ vanish. The radial boundaries are also closed for the flow, and the radial derivative of uz vanishes. For the azimuthal flow, we have to adopt a modified boundary condition though. A Keplerian rotation profile is maintained into the boundary zones, based on the last zone of the integration domain. If the first azimuthal velocity of the computational domain is denoted by $u_\phi^{(1)}$, we use
$\displaystyle u_\phi^{(0)}\phantom{^{-}}$ = $\displaystyle u_\phi^{(1)} \sqrt{R^{(1)}/R^{(0)}}$ (49)
$\displaystyle u_\phi^{(-1)}$ = $\displaystyle u_\phi^{(1)} \sqrt{R^{(1)}/R^{(-1)}.}$ (50)

The outer radial boundary is treated accordingly for  $u_\phi^{(n+1)}$ and  $u_\phi^{(n+2)}$. The azimuthal boundary conditions are periodic.

A rough representation of a stratified disk with radius-independent density scale-height, was used for the initial conditions. We leave this configuration for free development under the influence of the gravitational potential. The initial conditions also involve a non-axisymmetric velocity perturbation of the form

uz = $\displaystyle A \sin m\phi,$ (51)
uR = $\displaystyle A \cos (m+1)\phi.$ (52)

The amplitude A is - in terms of the Keplerian velocity in the middle of the annulus, $U_{\rm K}$ - about $7\times 10^{-4}U_{\rm K}$. We have run two models with m=1 and different resolutions and two additional with m=10 and m=200 in order to test short-wave perturbations. The models are summarized in Table 1. The initial configuration is isothermal with an energy density of

\begin{displaymath}e=\frac{1}{\gamma(\gamma-1)} ~\rho~ c_{\rm ac}^2, \end{displaymath} (53)

where the sound speed $c_{\rm ac}$ is 7% of the Keplerian velocity $U_{\rm K}$, and the polytrope exponent is $\gamma=5/3$.
   
Table 1: Overview of simulation configurations. The quantity m denotes the azimuthal mode of the initial perturbation.
Model Resolution m
I $31\times 61\times 351$ 1
II $71\times 71\times 351$ 1
III $31\times 61\times 701$ 10
IV $31\times 61\times 701$ 200

Figure 1 shows the evolution of the directional kinetic energies for the z- and r-components. A gradual decay of fluid motion is observed. Velocity fluctuations can also be measured in terms of the correlation tensor element $Q_{R\phi}$ which is non-dimensionalized with the average pressure giving the Shakura-Sunyaev parameter $\alpha _{\rm SS}$ (cf. Arlt & Rüdiger 2001). Figure 2 shows the evolution of $\alpha _{\rm SS}$ during 30 rotation periods at the inner disk radius. A clear relaxation of the flow is seen. The initially strong velocity fluctuations settle the correct density stratification for the gravitational potential of the point mass (central star). The fluctuations in the vertical direction reach $\pm14$% of the Keplerian velocity $U_{\rm K}$. Those in the radial direction even reach -27%. They are thus not small at all, and even a nonlinear instability might be noticed if existing.
  \begin{figure} \psfig{figure=MS2380f1.ps,width=8.8cm,height=7cm} \end{figure} Figure 1: Kinetic energies in the vertical and radial components of the flow from the run with $31\times 61\times 351$ grid points (Model I). The solid line is the energy of the r-direction, the dashed line that of the z-direction. Times are given in revolutions of the inner boundary of the computational domain.
Open with DEXTER


  \begin{figure} \psfig{figure=MS2380f2.ps,width=8.8cm,height=7cm} \end{figure} Figure 2: Shakura-Sunyaev parameter $\alpha _{\rm SS}$ from the same model as in Fig. 1 (Model I). Times are again in orbital periods of the inner boundary of the domain.
Open with DEXTER

We can read an approximate decay time $\tau_{\rm D}$ from Fig. 1. On average, it amounts to an e-folding time of about 20 orbital periods or $\tau_{\rm D}=3.18$ in non-dimensional units. As the main relaxation motions are vertical, we use the density scale-height for the length scale and obtain a viscosity of 0.03 caused by the finite resolution of the numerical grid. In terms of Reynolds numbers, we achieve

\begin{displaymath}{\rm Re} = \frac{U_{\rm in}\tau_{\rm D}}{H_\rho} = 1520 \end{displaymath} (54)

near the inner edge of the annulus where the Keplerian velocity is $U_{\rm in}=158$. The density-scale height at the inner edge is $H_\rho=0.33$ in the initial configuration. A second run with significantly higher resolution - in vertical direction in particular - is shown in Fig. 3. No change in the decaying nature of the perturbation is found. An increase in azimuthal wave number of the perturbation does not provide instability either. Figure 4 is the result of a simulation with a perturbation of very high wave number, $u_r=A\sin 200\phi$, $u_z=A\cos 201\phi$, and doubled azimuthal resolution. The gas again settles towards an equilibrium state and very low kinetic energies in the radial and vertical components.
  \begin{figure} \psfig{figure=MS2380f3.ps,width=8.8cm,height=7cm} \end{figure} Figure 3: Kinetic energies in the vertical (dashed) and radial (solid) components of the flow. The resolution of the model is $71\times 71\times 351$ grid points (Model II). Times are again in orbital periods of the inner boundary of the domain.
Open with DEXTER


  \begin{figure} \psfig{figure=MS2380f4.ps,width=8.8cm,height=7cm} \end{figure} Figure 4: Kinetic energies in the vertical (dashed) and radial (solid) components for the model with m=200 in the initial perturbation (Model IV). The resolution of the computational domain is $31\times 61\times 701$.
Open with DEXTER

4 Conclusions

We have shown in Sect. 1 that the Solberg-Høiland conditions (4) and (5) answer the question whether the combined action of stable (radial) shear and stable (vertical) density stratification can produce a new instability. The clear answer is No, but uncertainties remain about validity of the Solberg-Høiland criterion for the case of nonaxisymmetric disturbances of the compressible accretion disk material. Indeed, in Sect. 2 we were able to reproduce the Solberg-Høiland criterion in the well-known formulation but only with a short-wave approximation and a small-m approximation applied formulated in (23). Self-gravitation (i.e. density waves) has been excluded from the consideration but sound waves have been included. Our dispersion relation (32) is of fourth order contrary to the dispersion relation of second order resulting from the Boussinesq approximation (Goldreich & Schubert 1967) and contrary to the dispersion relation of fifth order resulting from nonideal hydrodynamics (Abramowicz et al. 1984, 1990). Along this way a new necessary condition for stability results which requires the external force which balances pressure and centrifugal force to be a conservative one, i.e. to possess a potential (cf. (36)). It is interesting to note that in the Boussinesq approximation for nonideal fluids the Solberg-Høiland criterion changes to the Goldreich-Schubert-Fricke criterion (Goldreich & Schubert 1967; Smith & Fricke 1975; Urpin & Brandenburg 1998) which no longer allows the existence of a vertical shear for stability, i.e. $\partial {\it\Omega}/\partial z =0$. And, indeed, even in thin accretion disks there is always such a (weak) vertical shear $\partial {\it\Omega}/\partial z$ the meaning of which remains open. We have thus numerically simulated the stability of such accretion disks using the ZEUS-3D code. The numerical integration of the nonlinear adiabatic equations of the fully compressible hydrodynamics with initially nonaxisymmetric velocity perturbations does not lead to the onset of instability although the disk, of course, does possess a small vertical shear. An overview of the simulated configurations can be found in Table 1 which shows that we indeed have also considered very high values of the azimuthal "quantum'' number m which do not fulfil the short-wave condition (23). Based on our numerical simulations, we thus cannot support findings about a fast, hydrodynamic instability in stratified disk flows. Similar simulations including magnetic fields showed the magneto-rotational instability with growth rates near $0.5 P_{\rm orb}^{-1}$. Our hydrodynamic runs cover 30-50 orbits of the inner disk and show only decay in velocity fluctuations. Angular momentum transport hovers at $\alpha_{\rm SS}\lower.4ex\hbox{$\;\buildrel <\over{\scriptstyle\sim}\;$ }\pm10^{-6}$; the temporal average in Fig. 2 is $5\times 10^{-8}$.

References

 
Copyright ESO 2002