A&A 407, 779-786 (2003)
DOI: 10.1051/0004-6361:20030867

On hydrodynamic shear turbulence in stratified Keplerian disks: Transient growth of small-scale 3D vortex mode perturbations

A. G. Tevzadze¹ - G. D. Chagelishvili¹ - J.-P. Zahn² - R. G. Chanishvili¹ - J. G. Lominadze ¹

1 - Center for Plasma Astrophysics, Abastumani Astrophysical Observatory, Kazbegi 2a, 380060 Tbilisi, Georgia
2 - LUTH, Observatoire de Paris, 92190 Meudon, France

Received 24 March 2003 / Accepted 3 June 2003

Abstract
This is a sequel to Paper I (Chagelishvili et al. 2003), where we presented the so-called bypass concept for the onset of turbulence in shearing flows. According to this concept, which was worked out during the last decade by the hydrodynamic community for spectrally stable flows, vortical perturbations undergo transient growth by extracting energy from the shear (a linear process), thereby reaching an amplitude which is sufficient to allow for non-linear interactions which, by positive feedback, sustain turbulence. In Paper I we described this transient growth for 2D perturbations in a Keplerian disk; we showed that their kinematics was the same as in plane-parallel flow, and thus that they were not modified by the presence of the Coriolis force. In the present paper, we pursue our goal of applying the bypass scenario to astrophysical disks: we investigate the linear dynamics of 3D small-scale vortical perturbations for single spatial harmonics, in stably stratified, differentially rotating disks, again in the framework of a nonmodal analysis. We find that these 3D perturbations also undergo substantial transient growth, and that they reach a peak amplitude that is comparable to that of 2D perturbations, as long as their vertical scale remains of the order of the azimuthal scale. When the vertical wave-number exceeds the azimuthal one, the amplification rate is reduced, but this may be more than compensated to by the huge Reynolds number and the high shear rate characterizing astrophysical Keplerian disks. Whereas in 2D the Coriolis force had no impact on the transient growth, in 3D this force somewhat constricts the characteristics of the perturbation dynamics in disk flows, and the initial transient growth is followed by some reduction in amplitude. These differences are quantitative, rather than of fundamental character. But the 3D case presents two interesting novelties. In plane parallel flow, the perturbations do not decay after their transient amplification, but their energy stays on a plateau before being dissipated through viscous friction. More importantly, especially for the astrophysicist, in disk flow the 3D vortex mode perturbations excite density-spiral waves, whose energy also settles on a plateau before viscous dissipation. These local vortex mode perturbations fit naturally into the bypass concept of hydrodynamic shear turbulence, which was first developed for plane-parallel flows. We submit that these perturbations will also play an important role in the onset and in the maintenance of turbulence in Keplerian disks.

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

1 Introduction

The prevailing explanation for the accretion process that operates in astrophysical disks is that these disks are turbulent, and that the turbulence is responsible for the inflow of matter and the outward transport of angular momentum. The question that then arises is that of the instability that produces this turbulence. This is not a trivial problem, because cylindrical rotation with a Keplerian profile is linearly (spectrally) stable, since angular momentum increases outward and since there is no extremum of vorticity. A solution has been obtained by Balbus & Hawley (1991), who showed that a weak magnetic field allows for a linear instability, called magneto-rotational instability (MRI). Subsequent numerical investigations followed that instability into the nonlinear (turbulent) regime, and they demonstrated that the turbulent transport then proceeds in the right direction. Thus for most astrophysicists, the question is settled: it is the MRI that is responsible for the turbulent transport.

However, this linear instability does not preclude the existence of other instabilities, and it is impossible to predict at the present stage which one is responsible for the actual turbulence in accretion disks. Also, a significant portion of planetary disks is not ionized enough to allow for the MRI. Thus it makes sense to explore alternate routes that may lead to turbulence, and to start with the hydrodynamic case. One is encouraged in this approach by the recent detection of turbulence in a Couette-Taylor experiment (Richard 2001), under conditions similar to those of a Keplerian disk, namely when angular momentum increases and angular velocity decreases outward.

One possible route to hydrodynamic turbulence is through transient growth of small amplitude perturbations, which extract their energy from the shearing flow. In this scenario, the amplitude of some vortex mode perturbations increases algebraically by a linear mechanism until it reaches the level where the nonlinear interactions, by redistributing energy among modes and providing positive feedback, are able to sustain turbulence. This concept, termed "bypass'', has been intensively developed in recent years by the fluid dynamics community; for a bibliography of the subject, we refer the reader to Paper I (Chagelishvili et al. 2003).

In Paper I, we analyzed the dynamics of 2-dimensional perturbations, with velocity components in the r (radial, shearwise) and $\varphi$ (azimuthal, streamwise) directions. We showed that, in the shearing sheet approximation (i.e. in the limit of large wavenumber, when curvature can be neglected), the energy equations are identical in the plane-parallel (Cartesian) case and in the Keplerian case, after suitable renormalization of the pressure. The same is true for the equations describing the temporal evolution of a spatial Fourier harmonics (SFH), meaning that the Coriolis force does not interfer with the linear dynamics. Singling out a vortex mode with initial wave-vector $[k_r(0), k_\varphi]$, we found that its energy grows by a factor $(k_r(0)/ k_\varphi)^2$, before decaying and being dissipated through viscous friction.

The present paper is a sequel to Paper I, extending its results to 3-dimensional perturbations, in a medium that is stably stratified in the z-direction; this medium can be either an astrophysical disk or the fluid in a laboratory experiment. When the linear scale of the perturbation is much smaller in the radial and azimuthal directions than in the direction perpendicular to the disk (the vertical direction), i.e. $l_r,~ l_{\phi} <\!\!< l_z,~H,~$ where H is the half thickness of the disk (or the pressure scale-height), the problem reduces to the 2D case treated in Paper I; it was examined by Lominadze et al. (1988) and by Fridman (1989). On the other hand, Ioannou & Kakouris (2001) considered the global perturbations where $l_r,~l_{\phi} >\!\!> H.$ Here we shall focus on perturbations that are influenced by the stratification, and that may come closer to isotropy, namely $l_r,~l_{\phi},~l_z <\!\!< H$.

A few words about the mathematical method used in our study, namely the non-modal analysis, that allowed for the significant advances achieved in this field during the last decade. After changing the variables from the laboratory to a co-moving frame, this non-modal analysis consists of describing the temporal evolution of the perturbations expanded in spatial Fourier harmonics (SFH), without performing any spectral expansion in time. This approach considerably simplifies the mathematical formalism of the problem and, completed by numerical calculations, it reveals important properties that are overlooked in the classical modal analysis.

In addition to the linear mechanism of non-exponential energy exchange between the shear flow and its perturbations, which is the basis of the bypass concept, a novel phenomenon was discovered through this non-modal approach, namely the linear mechanism of wave excitation by vortices (see Chagelishvili et al. 1997, 2000). Indeed, we shall see below that the 3D vortex mode perturbations, after transient growth, excite density-spiral waves, which then participate also in the dynamical processes in Keplerian disks.

The paper is organized as follows. The physical approximations and the mathematical formalism are introduced in Sect. 2. In Sect. 3 we present the numerical analysis of the linear dynamics of 3D vortex mode perturbations, including the excitation of density-spiral waves. We summarize and discuss the results in Sect. 4.

2 Mathematical formalism

Consider a differentially rotating thin inviscid hydrodynamic disk around a central object located at r=z=0. For simplicity the equilibrium disk is assumed to be non self-gravitating, and to be locally isothermal, with constant adiabatic sound speed. The relations between the basic physical characteristics for this disk, such as angular velocity of rotation $(\Omega)$, sound speed  $(c_{\rm s})$, pressure and density scale height (H), local mean values of pressure (P₀) and density $(\rho_0)$, are well known (cf. Ryu & Goodman 1992):

$$\Omega \propto r^{-3/2}, ~~ c_{\rm s}^2 = \gamma {P_0 \over \... ...s} \over \sqrt{\gamma} \Omega} = {c_{\rm s}^2 \over \gamma g},$$ (1)

$${P_0(z) \over P_0(0)} = {\rho_0(z) \over \rho_0(0)} = \exp \left(-{\vert z\vert \over H}\right),$$ (2)

where $\gamma$ is adiabatic index that meets the condition $\gamma \geq 1$. In Eqs. (1), (2) and further on, the vertical gravity $\Omega^2z$ is replaced by the parameter g, which we shall assume constant apart from a change of sign at the mid-plane of the disk. The expressions above are valid also for a laboratory experiment, H being the pressure scale-height.

As in Paper I, we shall write the dynamical equations in the local co-moving Cartesian frame at the particular point ${\vec r} = (r_0,0,0)$:

$$x \equiv r - r_0, ~~~~~ y \equiv r_0 (\phi - \Omega_0 t), ~~~~~ z = z ,$$ (3)

where $\Omega_0=\Omega(r_0)$ is the rotation rate of the frame. The perturbations we shall consider have characteristic length-scales much smaller then the distance to the central object:

$${l_x \over r_0},~ {l_y \over r_0} , ~{l_z \over r_0} \ll 1 .$$ (4)

Therefore we may neglect the curvature of the flow streamlines, but we retain the effect of the background velocity shear in the radial direction (shearing sheet approximation). To simplify the equations, we normalize by $\rho_0(0)$ the pressure and density perturbations p' and $\rho'$, and we filter out the vertical dependence of the velocity perturbation field  ${\vec u^\prime}$ by factorizing it through $\exp ({\vert z\vert/ H})$:

$$P = P_0(z) + {\rho_0(0)}p^\prime , ~~~ \rho = \rho_0(z) + {\rho_0(0)} \rho^\prime .$$ (5a)

$${\vec V} = {\vec V}_0 + \exp \left({\vert z\vert \over H}\right) {\vec u}^\prime .$$ (5b)

The equations governing the linear dynamics of the perturbation are thus

$$\left({\partial \over \partial t} + 2Ax {\partial \over \part... ...mega_0 u_y^\prime + {\partial p^\prime \over \partial x} = 0 ,$$ (6a)

$$\left({\partial \over \partial t} + 2Ax {\partial \over \part... ...0 + A) u_x^\prime + {\partial p^\prime \over \partial y} = 0 ,$$ (6b)

$$\left({\partial \over \partial t} + 2Ax {\partial \over \part... ... {\partial p^\prime \over \partial z} + g {\rho^\prime } = 0 ,$$ (6c)

$$\left({\partial \over \partial t} + 2Ax {\partial \over \partial y}\right) {\rho^\prime} + \nabla {\vec u}^\prime = 0 ,$$ (6d)

$$\left({\partial \over \partial t} + 2Ax {\partial \over \part... ...rm s}^2 \nabla {\vec u}^\prime + (\gamma-1) g u_z^\prime = 0 ,$$ (6e)

where

$$A \equiv {1 \over 2} \left[ r {\partial \Omega \over \partial r} \right]_{r=r_0}$$ (7)

is the standard Oort constant which characterizes the background velocity shear in the local frame (in Paper I we used another definition for A). In the case of the Keplerian law of rotation:

$$A = -{3 \over 4} \Omega_0 < 0.$$ (8)

Next, following the standard method of non-modal analysis (cf. Goldreich & Lynden-Bell 1965; Goldreich & Tremaine 1978; Nakagawa & Sekiya 1992) we introduce the spatial Fourier harmonics (SFH) of the perturbations with time dependent phases:

$$\psi^\prime ({\vec r},t) = \psi(k_x(t),k_y,k_z,t) \exp (- {k... ...\left( {\rm i} k_x(t) x + {\rm i}k_y y + {\rm i}k_z z \right),$$ (9a)

k_x(t) = k_x(0) - 2A k_y t, (9b)

where $\psi^\prime \equiv (u_x^\prime, u_y^\prime, u_z^\prime, p^\prime, \rho^\prime)$, $\psi \equiv (u_x, u_y, u_z, p, \varrho)$ and $k_H \equiv 1/(2H)$. The common factor $\exp (- {k_H} \vert z\vert)$ is employed to re-scale the vertical axis in the stratified medium (see Lerche & Parker 1967).

Equations (6a)-(6e) and (9a,b) lead to the dynamical equations for the SFH:

$${{\rm d} \over {\rm d} t}u_x - 2\Omega_0 u_y + {\rm i}k_x(t) p = 0 ,$$ (10a)

$${{\rm d} \over {\rm d} t}u_y + 2(\Omega_0 + A) u_x + {\rm i}k_y p = 0 ,$$ (10b)

$${{\rm d} \over {\rm d} t}u_z + ({\rm i} k_z - k_{\rm H}) p + g \varrho = 0 ,$$ (10c)

$${{\rm d} \over {\rm d} t}\varrho + {\rm i}\left[k_x(t) u_x + {\rm i} k_y u_y + \left(k_z + {\rm i}k_H\right)u_z\right] = 0 ,$$ (10d)

$${{\rm d} \over {\rm d} t}p + {\rm i}c_{\rm s}^2\left[k_x(t) u... ...(k_z + {\rm i} k_H\right) u_z\right] + (\gamma-1) g ~ u_z = 0.$$ (10e)

The conservation of the potential vorticity is described in the considered case by the following time invariant solution of Eqs. (10a)-(10e):

$${\cal I} = \left[ k_y u_x - k_x(t) u_y \right] - 2{\rm i} (\O... ...rm H}) \over (\gamma-1)g} \left( p - c_{\rm s}^2 \rho \right).$$ (11)

This invariant signals the existence of a vortex/aperiodic mode in the perturbation spectrum.

To avoid complex coefficients (and consequently complex physical variables) in the dynamical equations, we renormalize the perturbed variables in the following way:

$${\tilde u_z} \equiv (1 + {\rm i}\alpha) u_z,$$ (12a)

$$\equiv {(1 + {\rm i}\alpha) k_z c_{\rm s}^2 \over (\gamma - 1)g} (p-c_{\rm s}^2 \varrho)$$                            s

$${(1 + {\rm i}\alpha) \gamma \over 2(\gamma - 1)}{k_z \over k_H} (p-c_{\rm s}^2 \varrho),$$ = (12b)

$${\tilde p} \equiv {\rm i}p,$$ (12c)

where

$$\alpha \equiv {k_H \over k_z} - {(\gamma-1) g \over {k_z c_{\rm s}^2}}= {k_H \over k_z}{2-\gamma \over \gamma} \cdot$$ (12d)

Hence the substitution of Eqs. (12a-d) into (10a-e) yields the system of ordinary differential equations which governs the dynamics of variables $u_x, u_y, {\tilde u_z}, s$ and  ${\tilde p}$:

$${{\rm d} \over {\rm d} t}u_x = 2\Omega_0 u_y - k_x(t) {\tilde p},$$ (13a)

$${{\rm d} \over {\rm d} t} u_y = -2(\Omega_0 + A) u_x - k_y {\tilde p} ,$$ (13b)

$${{\rm d} \over {\rm d} t} {\tilde u_z} = {2k_H \over \gamma} s - k_z (1 + \alpha^2) {\tilde p} ,$$ (13c)

$${{\rm d} \over {\rm d} t} s = - {k_z c_{\rm s}^2} {\tilde u_z} ,$$ (13d)

$${{\rm d} \over {\rm d} t} {\tilde p} = c_{\rm s}^2 \left[ k_x(t) u_x + k_y u_y + k_z {\tilde u_z} \right] .$$ (13e)

All variables $u_x, ~u_y,~{\tilde u_z}, ~{\tilde p}, ~s$ have the same phase, which we will set to zero, thus making these variables real. The original variables $u_z,~\varrho~$ are complex, and p is pure imaginary; for instance $\varrho$ is expressed by:

$${\varrho} = - {{\rm i}\over c_{\rm s}^2}{\tilde p} - {1-{\rm ... ... {2 (\gamma - 1) \over \gamma c_{\rm s}^2} {k_H \over k_z} s ,$$ (14a)

$${{\rm Re}\varrho} = - {1 \over 1 + {\alpha}^2}{2 (\gamma - 1) \over \gamma c_{\rm s}^2} {k_H \over k_z} s ,$$ (14b)

$${{\rm Im}\varrho} = - {1\over c_{\rm s}^2}{\tilde p} + { \alp... ... {2 (\gamma - 1) \over \gamma c_{\rm s}^2} {k_H \over k_z} s .$$ (14c)

In the new variables the conservation of potential vorticity (Eq. (11)) assumes the simple form:

$${\cal I} = k_y u_x - k_x(t) u_y - {2(\Omega_0 + A) \over c_{\rm s}^2} \left( {\tilde p} - s \right).$$ (15)

The spectral energy density of the perturbations is defined as follows:

$$E(t) = E_{\rm k}(t) + E_{\rm p}(t) + E_{\rm t}(t),$$ (16a)

with

$$E_{\rm k}(t) = {1 \over 2}\rho_0(0) \left( u_x^2 + u_y^2 + {{\tilde u_z}^2 \over {1 + \alpha^2}} \right),$$ (16b)

$$E_{\rm p}(t) = { \rho_0(0) \over 2 c_{\rm s}^2}~ {\tilde p^2},$$ (16c)

$$E_{\rm t}(t) = {2k_H \over \gamma k_z} {1 \over {1 + \alpha^2}} ~ { \rho_0 (0) \over 2 c_{\rm s}^2 }~ s^2 ,$$ (16d)

where $E_{\rm k}$, $E_{\rm p}$ and $E_{\rm t}$ correspond respectively to the kinetic, elastic and thermobaric energy spectral densities of the perturbations (cf. Gossard & Hooke 1975).

Note that the spectral density of the energy would be conserved in the shearless limit; its variation is due to the velocity shear of the flow:

$${{\rm d} \over {\rm d} t} E = - 2A {\rho_0(0)} u_x u_y .$$ (17)

Formally, Eqs. (1), (2), (6a)-(6e) and the further mathematical analysis are well-posed for any value of the ratio k_z / k_H. But for simplicity the vertical gravity $\Omega^2z$ has been replaced by the constant parameter g. Thus when applying our results to astrophysical disks, we have to restrict ourselves to perturbations whose vertical scale is much shorter than the scale-height H; hence $k_H / k_z~\ll~1$.

2.1 Perturbation spectrum in the shearless limit

The dispersion equation for our system may be obtained in the shearless limit (A=0) using the full Fourier expansion of the variables, including time ( $\psi \propto \exp({\rm i} \omega t)$):

$$\omega ~ \left\{ \omega^4 - \left( c_{\rm s}^2 k^2 + 4 \Omega_0^2 \right) \omega^2 \right.$$

$$\left. \quad+ c_{\rm s}^2 \left[ N_{\rm B}^2 k_\perp^2 + 4 \Omega_0^2 \left(k_z^2 + k_H^2\right) \right] \right\} = 0$$ (18)

where $k_\perp^2 \equiv k_x^2+k_y^2$, $k^2 \equiv k_\perp^2 + k_z^2 + k_H^2$. (The local epicyclic frequency is here $2 \Omega_0$.) $N_{\rm B}$ is the Brunt-Väisälä frequency:

$$N_{\rm B}^2 \equiv (\gamma-1) {g^2 \over c_{\rm s}^2} >0 .$$ (19)

This dispersion equation describes three different modes of perturbations:
1. a high frequency rotational-acoustic wave mode with

$$\omega_s^2 ={ 1 \over 2} \left( c_{\rm s}^2 k^2 + 4 \Omega_0^... ...^2 k^2 + 4 \Omega_0^2\right)^2 } \right)^{1 \over 2} \right] ,$$ (20)

or, in $k_H / k_z \ll 1$ approximation, $\omega_{\rm s}^2 \simeq { c_s^2 k^2 + 4 \Omega_0^2 };$
2. a low frequency density-spiral wave mode with

$$\omega_{g\Omega}^2 ={ 1 \over 2} \left( c_{\rm s}^2 k^2 + 4 \... ...^2 k^2 + 4 \Omega_0^2\right)^2 } \right)^{1 \over 2} \right] ,$$ (21)

or, in $k_H / k_z \ll 1$ approximation, $\omega_{g\Omega}^2 \simeq N_{\rm B}^2 {k_\perp^2 / k^2} + 4 \Omega_0^2 {k_z^2 / k^2};$
3. a vortex mode with

$$\omega = 0.$$ (22)

The perturbation field of the vortex mode may be derived from Eqs. (10a)-(10e) and (13a)-(13e). Indeed, the stationary solution of this system is not trivial and reads as follows:

$$u_x = -{ k_y \over 2 \Omega_0} {\tilde p}, ~~ u_y = {k_x \ove... ...z = 0, ~~ \varrho = - {k_z + {\rm i} k_H \over g} {\tilde p} ,$$

$$s = {\gamma k_z \over 2k_H} (1 + \alpha^2) {\tilde p} = {N_{\rm B}^2 - c_s^2 k_z^2 \over N_{\rm B}^2} {\tilde p}.$$ (23)

Notwithstanding the formal similarity, this vortex mode is different from the vortex mode occurring in the simple not-stratified compressible hydrodynamic flow: the difference is in the topological properties, as well as in the genesis of the mode. The perturbations of the present vortex mode are polarized in the flow rotation plane; they are 2-dimensional in the shearless limit ( $u_z \rightarrow 0$ when $A \rightarrow 0$). This mode originates from the combined action of both the vertical gravity and the Coriolis force. In the absence of any one of these forces this mode degenerates into the trivial solution of the system - it disappears. In the case where both forces are absent, the density-spiral wave mode ($g\Omega$ mode) degenerates into the vortex mode that is the key ingredient of turbulence in the simple (plane and non-stratified) hydrodynamic shear flow.

In previous works, the role of this vortex mode has been often underestimated (cf. Dwarkadas & Balbus 1996; Ryu & Goodman 1992; Goodman & Balbus 2001), or it has been confused with fictitious displacements that arise in the Lagrangian formalism (see Friedman & Schutz 1978). This mode has been ignored also in the recent paper by Goodman and Balbus (2001), which retains only the wave modes.

2.2 Velocity shear effects

The velocity shear has a profound effect on the perturbations, since it permits the extraction of energy from the mean flow (see Eq. (18)). This results in the transient growth of the vortex mode perturbations - the subject of our study. This phenomenon has the same nature as the amplification of vortical perturbations in parallel flows with constant shear rate (cf. Moffatt 1967; Criminale & Drazin 1990; Gustavsson 1991; Reddy & Henningsson 1993; Chagelishvili et al. 1994). Moreover, as was found and described in Chagelishvili et al. (1997, 2000), in relatively complex flows, where vortex mode coexists with wave modes (which is the case of Keplerian flow); the velocity shear causes their linear coupling: vortices are able to excite waves even in the linear approximation. According to Chagelishvili et al. (1997), the efficiency of the wave excitation mechanism is determined by the value of the ratio of the shear parameter to the frequency of the considered wave mode:

$$R \equiv \vert 2A\vert/\omega(t^*),$$ (24)

where the wave frequency is taken at time t=t^* where k_x(t^*)=0, when the frequency is minimal and the mode coupling is strongest. We may estimate crudely this efficiency factor by taking the wave frequencies derived in the shearless limit (Eqs. (20)-(21)). For the production of density-spiral waves its value is:

$$R_{g\Omega} {2\vert A\vert \over \sqrt{4 \Omega_0^2 {k_z^2 \over k^2} + N_{\rm B}^2 {k_y^2 \over k^2}}}$$ =

$${2\vert A\vert \over \sqrt{4 \Omega_0^2 {k_z^2 \over k^2} + {\gam... ...2 {k_y^2 \over k^2}}} \approx {\vert A\vert \over \Omega_0} = {3 \over 4} \cdot$$ = (25)

For the production of rotational-acoustic waves it is:

$$R_{\rm s} = {2\vert A\vert \over \sqrt{c_{\rm s}^2 (k_y^2+k_z... ..._0^2}} \simeq {3 k_H \over 2 \sqrt{\gamma(k_y^2+k_z^2)}} \cdot$$ (26)

Let us compare them:

$${R_{g\Omega} \over R_{\rm s}} \simeq \left( {\gamma(k_y^2+k_z^2)\over 4 k_H^2}\right)^{1/2}\cdot$$ (27)

In the considered problem $k_H / k_z \ll 1$ and hence $R_{\rm s} \ll R_{g\Omega}$. Thus one should expect mainly the excitation of density-spiral waves. Our calculations justify this expectation.

3 Numerical analysis

Figure 1: Dynamics of the perturbation SFH in the Keplerian disk ( left column) and plane shear flow ( right column). From top to bottom: $u_x(t)/(\rho _0 c_{\rm s})$, $u_y(t)/(\rho _0 c_{\rm s})$, ${\tilde u_z(t)}/(\rho_0 c_{\rm s})$, ${\tilde p(t)}/P_0$, ${{\rm Re}\varrho / \rho _0}$ and ${{\rm Im}\varrho /\rho _0}$ (heavy and thin curves, respectively; same in the plane flow) and normalized energy E(t)/E(0) (the latter on logarithmic scale). Here kx(0)/ky=-150, kz/ky=0.1. In the AD case: kH/kz = 0.1 and $R_{\rm s} = 0.1 \ll R_{g\Omega } = 3/4$. The initial values the of perturbations correspond to the quasi-two dimensional pure vortex mode. Initial value of the spectral energy density in both cases $E(0)/(\rho _0 c_{\rm s}^2) = 10^{-3}$.

Open with DEXTER

	Figure 2: The same as in Fig. 1, but at kz/ky=1. Initial perturbations correspond to the pure vortex mode with equal vertical and azimuthal length-scale.
Open with DEXTER

	Figure 3: The same as in Fig. 1, but at kz/ky=10. Initial perturbations correspond to the pure vortex mode with length-scale much shorter in the vertical then in the azimuthal direction.
Open with DEXTER

Equations (13), (14), together with the appropriate initial values, pose the initial value problem describing the dynamics of a perturbation SFH in the Keplerian disk. We solve this system with initial values corresponding to pure vortex mode perturbations. The numerical calculations are performed using the Matlab ODE solver, an explicit Runge-Kutta implementation with the 4(5) pair of Dormand-Prince. Potential vorticity is conserved with an accuracy of 10^-7.

Since we are primarily interested in the transient growth of vortex perturbations in AD and its comparison with the similar process in plane shear flow, we analyze the dynamics of plane shear flow in the same manner. For this purpose we set g=0, k_H=0, $\alpha=0$, $\Omega_0=0$ and A < 0  in Eqs. (13), (14). The evolution of vortex mode SFH, in disk flow and plane shear flows, are presented in Figs. 1, 2 and 3 for different ratios of the vertical and azimuthal wave-numbers k_z/k_y: 0.1,  1, 10; in all cases k_x(0)/k_y = - 150. All calculations for the plane case are carried out with the same kinematic shear parameter. (When comparing the figures, note that the graphs are not on the same scale.)

The numerical results demonstrate the effect of rotation and vertical stratification on the transient growth in Keplerian disks. We see that quasi-2D vortex mode perturbations (i.e. with $k_z \ll k_y$) are amplified by nearly the same factor as in plane shear flow (Fig. 1). This amplification is slightly less in the $k_z \sim k_y$ case (Fig. 2), whereas the transient growth is strongly reduced for vortex mode perturbations with $k_z \gg k_y$(Fig. 3).

But the calculations reveal also a novel linear effect which accompanies the evolution of vortex mode perturbations in the Keplerian disk, namely the excitation of density-spiral waves ($g\Omega$ modes, or gravito-inertial waves), which may play an important role in the onset and in the maintenance of turbulence in AD. For a physical description of this excitation we refer to Chagelishvili et al. (1997). The phenomenon is clearly visible in Figs. 1-3, in the graphs of u_z, ${\rm Re} \varrho~$ and E(t)/E(0). Starting with a pure the vortex/aperiodic SFH at k_x(0)/k_y<0, the perturbation energy reaches a peak value at time $t=t^*\equiv k_x(0) / 2 A k_y$. After that maximum, the SFH undergoes nearly periodic oscillations, meaning that the vortex mode excites a density-spiral wave corresponding to the same SFH; the result is a mixed vortex-wave SFH.

While the vortex energy steadily decreases after the maximum, as in the 2D case considered in Paper I, the wave energy remains constant, and is responsible for the plateau in the E(t) graphs. (In that respect, the disk case bears some similarity with the plane case, where the energy also settles on a plateau.) The plateau energy is more then one order of magnitude lower then the peak energy. For k_z/k_y = 0.1 (Fig. 1) the normalized peak energy is $2.23 \times 10^4$ (very close to the theoretical value $(k_x(0)/k_y)^2 = 2.25 \times 10^4$ derived in Paper I) and the plateau energy is 10.9. For k_z/k_y = 1 (Fig. 2) the peak energy is slightly lower: $2.00 \times 10^4$, whereas that of the plateau is much higher than in the preceding case: $2.95 \times 10^2$. This indicates that the excitation of the wave mode is much stronger when $k_z \sim k_y$.

Note that in the plane case $u_x ({\vec k}) ~u_y({\vec k})$ is always positive. Since the flux of angular momentum is determined by

$$<u_x^\prime ({\vec r}) ~ u_y^\prime({\vec r})> \ = \int {\rm d}{\vec k}~ u_x({\vec k})~ u_y({\vec k}) ,$$ (28)

we see that all SFH then contribute to transport directed down the gradient of mean velocity. However in the disk case, as we can see in Fig. 2, $u_y({\vec k})$ changes sign at t=t^*. Thus the product $u_x ({\vec k}) ~u_y({\vec k})$ also changes sign, and becomes negative in the domain of wavenumber space where k_x   k_y >0. According to Eq. (28), this means that in the disk case there is some cancellation in the transport of angular momentum. But the net contribution of each SFH is outward transport, as shown by integration in time of Eq. (17), since net positive energy is extracted from the mean flow (see Fig. 2). The same property carries over in the turbulent regime, because the nonlinear terms just redistribute energy among modes, and therefore

$${{\rm d} \over {\rm d} t} \int {\rm d}{\vec k} ~ E = - 2A ~ {\rho_0(0)} \int {\rm d}{\vec k} ~ u_x ~ u_y.$$ (29)

Thus the net angular momentum flux is directed outward when energy is extracted from the mean flow, to be dissipated by viscous friction. In other words, if shear flow turbulence is sustained, it transports angular momentum outward.

As it appears, the action of the Coriolis force reduces the rate of transient amplification and it somewhat changes its character in 3D. Consider the case k_z = k_y (Fig. 2). One sees that, in the plane flow, the SFH energy increases transiently, but tends monotonically and smoothly to some  $E_{\rm max}$ when $t \to \infty$( $k_x(t)/k_y \to \infty$). In the disk case, the SFH energy reaches its peak value $E_{\rm peak}$ at k_x(t)=0, but after that it decreases, and settles on a plateau with energy $E_{\rm plat}$ at a time where $k_x(t)/k_y \approx 10$. To estimate quantitatively the effect of rotation, we have to compare $E_{\rm max}$ with $E_{\rm peak}$or/and $E_{\rm plat}$. These quantities are determined by the values of k_x(0)/k_y and k_z / k_y: substantial transient growth occurs when k_z is less or of the order of k_y. Let us introduce the parameter

$$\beta \equiv \left(k_x(0) \over k_y \right)^2$$ (30)

and calculate the dependence of $E_{\rm max},~ E_{\rm peak},~ E_{\rm plat}$ on $\beta$, for the same E(0) and k_z = k_y. The results of our calculations are plotted in Fig. 4. We see that the dependence of the maximal amplification rate in the plane and AD case is qualitatively the same, and that it scales linearly with $\beta$, as was shown in Paper I for 2D perturbations. To reach the same maximal amplification in the disk and plane flows (i.e. $E_{\rm peak} \simeq E_{\rm max}$), the parameter $\beta$ is less than one order of magnitude higher in the disk than in the plane flow. To reach $E_{\rm plat} \simeq E_{\rm max}$, $\beta$ ought to be about 200 times higher.

	Figure 4: Maximum normalized energy $E_{\rm max}$ in the plane case, peak energy $E_{\rm peak}$ and plateau energy  $E_{\rm plat}$ in the disk case, vs. $\beta \equiv (k_x(0)/ k_y)^2$, for the same E(0) and kz = ky.
Open with DEXTER

	Figure 5: E(t)/E(0) vs. t for the Keplerian disk and plane cases, for kz = ky and for values of $\beta$ at which $E_{\rm max}/E(0)=E_{\rm plat}/E(0) = 5 \times 10^4$ ( $\beta_{\rm plane} = 2.25 \times 10^3~$ and $\beta _{AD}= 4.41 \times 10^6$).
Open with DEXTER

Let us analyze the above results in terms of the Reynolds number, which we define as

$$R_e = {\Omega_0 H^2 \over \nu} ,$$ (31)

$\nu$ being the kinematic viscosity, and H the height of the disk (or pressure scale height). The smallest perturbation scale whose transient growth can overcome viscous dissipation may be determined by the local balance between transient growth and viscous dissipation (for details see Paper I). Thus, from (17) we have

$$-2 A~ u_x u_y = \nu k^2 u^2 .$$ (32)

We have seen that optimal growth is achieved when both

$$\beta = \left({k_x(0) \over k_y} \right)^2 \gg 1 \qquad \hbox{and} \quad k_y \approx k_z .$$ (33)

At $\vert{k_x/ k_y}\vert \gg 1$, the vortex mode compressibility is negligible and the velocity field is mainly parallel to the disk plane $(u_z \ll u_x,u_y)$. Thus, from the continuity equation it follows that

$$\vert k_x u_x\vert \approx \vert k_y u_y\vert.$$ (34)

Applying these equations to the initial value k_x = k_x(0), we get the following relation between R_e and $\beta$:

$$R_e \approx {2 \over 3} ~\beta^{3/2} .$$ (35)

This permits to transpose the results given above, namely the amplification factors vs. $\beta$ given in Figs. 4, 5, in terms of the Reynolds number R_e. For instance, we see that to achieve the same amplification in plane and disk flow (i.e. $E_{\rm max}=E_{\rm plat}$), the Reynolds number must be about 2000 times higher in the disk flow than in the plane flow.

4 Discussion

In the present paper we focused our attention on the transient growth of small scale vortices in stably stratified Keplerian disk flow. For the sake of comparison, the analysis of the disk flow was performed in parallel with that of the non-rotating plane shear flow. Let us recall the main results concerning linear dynamics.

• The stratification of the disk flow (vertical gravity) is necessary for the existence of the vortex/aperiodic mode in 3 dimensions. In the case of rigid rotation (A = 0) the vortex mode is characterized by a planar velocity field parallel to the disk plane. For $A \not=0$ this mode involves vertical velocity perturbations (u_z).

• The vortex mode SFH, with initial wave-numbers satisfying the condition k_x(0)/k_y < 0 (as A<0), is able to transiently amplify by several orders of magnitude through a linear mechanism. Peak energy is reached at $t = t^*(\equiv k_x(0)/2Ak_y)$. This peak value scales linearly with the parameter $\beta =(k_x(0)/k_y)^2$; it is highest for $\vert k_z/k_y\vert \la 1$, where it is comparable to that of the 2D case considered in Paper I, and it diminishes rapidly for $\vert k_z/k_y\vert\gg1$.

• After reaching its peak, the vortex mode SFH gives rise to a density-spiral wave of same SFH. Further evolution of the vortex and the excited density-spiral wave proceeds independently. Whereas the vortex energy decreases steadily, the wave energy is nearly maintained in time - it forms a plateau in the E(t)graphs. The wave excitation process is most efficient for vortices with nearly equal vertical and azimuthal wave-numbers: $k_z \sim k_y$. In this case, the plateau energy is about two orders of magnitude lower than the peak energy.

• The inhibiting action of the Coriolis force on the transient growth mechanism is more quantitative then qualitative - to reach the same plateau energy as in parallel shear flow, the vortex mode in disk flow requires a Reynolds number three orders of magnitude higher. The contrast is much lower when comparing the peak energies.

What conclusions can we draw from of our study?

The transient amplification of vortices, together with the phenomenon of wave excitation by these vortices, appear to be a promising mechanism leading to the onset of turbulence in astrophysical disks. These vortex mode perturbations fit naturally into the recently developed concept of bypass transition to turbulence for plane smooth shear flows. The detailed description of the concept was presented in Paper I. Here we only recall the four basic steps of the bypass scenario:
(i) the linear "drift" of SFH of the vortex mode in the k-space;
(ii) the transient growth of SFH;
(iii) viscous dissipation;
(iv) nonlinear processes that close the feedback loop of the transition by mode mixing and angular redistribution of SFH in the k-space.

In this scenario the nonlinear processes do not contribute to the energy growth, as they would in genuine nonlinear instability, but they serve to close the loop, providing positive feedback.

In the classical bypass concept, the main contributors to turbulence in shear flows are the vortex mode perturbations. But in Keplerian disks, according to the presented study, after transient growth, vortex mode perturbations are converted into density-spiral waves, and these long-lived waves become ingredients of the system. Hence, most likely, wave-wave nonlinear interactions will contribute significantly to step (iv). Moreover, in 3 dimensions much richer possibilities exist for such interactions than in 2 dimensions, and therefore more opportunities are offered for positive feedback.

To conclude, the linear phase of the bypass concept is now well established, and this scenario appears as a very promising route to turbulence in Keplerian disks. It remains necessary to confirm the existence of positive nonlinear feedback, which is required to sustain turbulence, and this can only be achieved through numerical simulations, performed at high Reynolds number. Thanks to the steady progress of computer performance, this goal should be reached in the not-so-distant future.

Acknowledgements

This work is supported by the International Science and Technology Center grant G-553. G.D.C. would like to acknowledge the hospitality of Observatoire de Paris (LUTH). The authors wish to thank the referee for helpful remarks.

References

 
• Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 NASA ADS  
• Chagelishvili, G. D., Rogava, A. D., & Segal, I. N. 1994, Phys. Rev. E, 50, R4283 In the text NASA ADS  
• Chagelishvili, G. D., Tevzadze, A. G., Bodo, G., & Moiseev, S. S. 1997, Phys. Rev. Lett., 79, 3178 In the text NASA ADS  
• Chagelishvili, G. D., Tevzadze, A. G., Bodo, G., Rossy, P., & Gogoberidze, G. T. 2000, in Proc. Eighth European Turbulence Conference Barcelona 2000, ed. C. Dopazo, 737 In the text  
• Chagelishvili, G. D., Chanishvili, R. G., Hristov, T. S., & Lominadze, J. G. 2002, JETP, 94, 434  
• Chagelishvili, G. D., Zahn, J.-P., Tevzadze, A. G., & Lominadze, J. G. 2003, A&A, 402, 401 (Paper I) In the text NASA ADS  
• Criminale, W. O., & Drazin, P. G. 1990, Stud. Appl. Math., 83, 123 In the text NASA ADS  
• Dwarkadas, D. D., & Balbus, S. A. 1996, ApJ, 467, 87 In the text NASA ADS  
• Fridman, A. M. 1989, Sov. Astr. Lett., 15, 487 [Russ: Pis'ma Astron. Zh. 15, 1122 ] In the text NASA ADS  
• Friedman, J. L., & Schutz, B. F. 1978, ApJ, 221, 937 In the text NASA ADS  
• Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125 In the text NASA ADS  
• Goldreich, P., & Tremaine, S. D. 1978, ApJ, 222, 850 In the text NASA ADS  
• Goodman, J., & Balbus, S. A. 2001 [arXiv:astro-ph/0110229v1] In the text  
• Gossard, E., & Hooke, W. 1975, Waves in the Atmosphere (Amsterdam: Elsevier) In the text  
• Gustavsson, L. H. 1991, J. Fluid Mech., 224, 241 In the text NASA ADS  
• Ioannou, P. J., & Kakouris, A. 2001, ApJ, 550, 931 In the text NASA ADS  
• Lerche, I., & Parker, E. N. 1967, ApJ, 149, 559 In the text NASA ADS  
• Lominadze, J. G., Chagelishvili, G. D., & Chanishvili, R. G. 1988, Sov. Astr. Lett., 14, 364 In the text NASA ADS  
• Moffatt, K. 1967, The interaction of turbulence with strong wind shear in Atmospheric turbulence and radio wave propagation, ed. A. M. Yaglom, & V. I. Tatarskii (Moscow: Nauka Press) In the text  
• Nakagawa, Y., & Sekiya, M. 1992, MNRAS, 256, 685 In the text NASA ADS  
• Reddy, S. C., & Hennigson, D. S. 1993, J. Fluid Mech., 252, 209 In the text NASA ADS  
• Richard, D. 2001, Thèse de doctorat, Université Paris 7 In the text  
• Ryu, D., & Goodman, J. 1992, ApJ, 388, 438 In the text NASA ADS