Dynamical Instability of Laminar Axisymmetric Flow of Perfect Fluid

with Stratification.
Zhuravlev V.V..1 , Shakura N.I..1,2
1) Sternberg Astronomical Institute, Moscow, Russia
∗

2) Max Planck Institute for Astrophysics, Garching, Germany

accepted to Astronomy Letters
arXiv:0709.1833v1 [astro-ph] 12 Sep 2007

key words: hydrodynamics – instability – non-homoentropic flows.

Abstract
The instability of non-homoentropic axisymmetric flow of perfect fluid with respect to non-axisymmetric in-
finitesimal perturbations was investigated by numerical integration of hydrodynamical differential equations in two-
dimensional approximation. The non-trivial influence of entropy gradient on unstable sound and surface gravity
waves was revealed. In particular, both decrease and growth of entropy against the direction of effective gravita-
tional acceleration gef f give rise to growing surface gravity modes which are stable with the same parameters in the
case of homoentropic flow. At the same time increment of sound modes either grows monotonically while the rate
of entropy decrease against gef f gets higher or vanishes at some values of positive and negative entropy gradient
in the basic flow. The calculations have showed also that growing internal gravity modes appear only in the flow
unstable to axisymmetric perturbations. At last, the analysis of boundary problem with free boundaries uncovered
that’s incorrect to set the entropy distribution according to polytropic law with polytropic index different from
adiabatic value, since in this case perturbations don’t satisfy the free boundary conditions.

∗ e-mail: slava@xray.sai.msu.ru

1
1 Introdution
This paper is natural extension of the two previous, in which the dynamical instability of axisym-
metric flows was considered in approximation of incompressible perfect fluid (Zhuravlev & Shakura,
2007a - hereafter ZS1) and with finite compressibility (Zhuravlev & Shakura, 2007b - hereafter
ZS2). Papaloisou and Pringle (1984, 1985, 1987) were the first who investigated the instability to
non-axisymmetrical perturbations being under discussion in this paper. Their discovery entailed
numerous explorations, which are mentioned partially in the review by Papaloisou & Lin (1995).
The growth of perturbations in axisymmetric flows with free boundaries is an important problem
for astrophysics and is of special interest in the view of question how angular momentum transports
outwards in accretion flows. Though the realistic astrophysical conditions force us to include the
effects of radiative heating and cooling, magnetic fields, relativity and so forth, it’s still necessary to
study the accretion flow as a shear flow in the simple hydrodynamical approximation. In the most of
such investigations the angular velocity profile was set by a power law Ω ∝ r −q with 1.5 < q ≤ 2.0,
which means that gef f attains the greatest value at the boundaries. In ZS1 and ZS2 alone with the
power law we considered another Ω profile which implies that gef f = 0 at the boundaries. It was
revealed that the absence of effective gravity at the edge of the flow essentially affects the growing
perturbation modes. In particular, it causes the disappearance of unstable sound modes so that the
flow with angular velocity profile being close to the Keplerian one becomes stable to infinitesimal
perturbations. Additionally, the special attention was paid to the influence of vorticity gradient
on the instability. In the present paper we are going to study how stratification of the basic flow
affects the growing sound and surface gravity waves and make an attempt to find unstable internal
gravity waves. Following the other authors we consider the barotropic basic flow, i.e. the entropy
distribution is set in such way that isopycnic surfaces coincide with isobaric ones. It allows to regard
angular velocity as function of radial direction. Discussing the instability of stratified axisymmetric
flows we should refer, for instance, to the paper by Frank & Robertson (1988) who considered the
instability of tori with random initial small perturbations and Kojima et al. (1989) who studied
toroidal as well as cylindrical flows. The latter are two-dimensional analogue of real flows in point-
mass gravitational potential. We should note here that in both cases the obtained results proved to
be similar and 3-dimensional growing modes of perturbations in the toroidal flow turned out to have
no essential dependence on vertical direction. The authors explained this by the fact that Reinolds
stress didn’t contain vertical perturbation velocity for barotropic configurations with Ω(r). Further,
Glatzel (1990) considered the instability of cylindrical and plane-parallel flows in approximation of
small extent of a shear layer. In order to exclude the growing sound and surface gravitational modes,
he assumed an incompressible fluid and rigid boundaries so that stratification was entered by the
density profile. The instability he found was explained as the result of over-reflection and coupling
of internal gravitational modes. Later Ghosh & Abramowicz (1991) studied the cylindrical flow con-
sisting of two fluids with different constant densities. The configuration was Rayleigh-Taylor stable.
Alone with modified branch of growing surface gravity mode being consequence of the free bound-
aries (Blaes & Glatzel, 1986) they found another branch of increments appeared owing to density
discontinuity. This instability is caused by growing internal gravitational mode, which, however, is
fully analogous to the growing surface mode, since the boundary between two fluids differs from the
free boundary just by the finite ratio of densities.
We should also mention about the results of Lovelace et al. (1999) and Li et al. (2000) where in the
2-dimensional geometry the instability of thin Keplerian disks with the local maximum of entropy was
investigated. In particular, the local dispersion relation for non-axisymmetric perturbations similar
to Rossby waves dispersion relation was derived. At last, Klahr & Bodenheimer (2003) studied the
instability of Keplerian disks with entropy decreasing outwards.
As it’s well known, the stability of stratified flows is governed by Richardson criterium. Initially it

2
was formulated for plane-parallel flows (see, for instance, Howard, 1961). In this case, the sufficient
condition for stability is Ri > 1/4 in each point, where Ri - is the Richardson number. The gen-
eralization of the Richardson criterium for baroclinic rotating flows, when angular velocity depends
also on the vertical direction, was found by Fujimoto (1987) for incompressible fluid and by Hanawa
(1987) for compressible fluid.
In the present work we follow other authors and consider the simple 2-dimensional geometry, i.e.
assume the cylindrical flow and perturbations with zero vertical velocity. In the first section we’ll
write down the basic equation for infinitesimal perturbations, formulate the free boundary conditions
and set the main features of the basic flow - angular velocity profile and entropy distribution. Then
we’ll touch on a numerical method, in particular, discuss the question when the perturbations can
satisfy the boundary conditions and finally will present our results.

2 Equation and boundary conditions.

2.1 Derivation of main equation
An isentropic motion of fluid is governed by the system of equations: (see, for example, Brekhovskikh
& Goncharov, 1982):
∂v ∇p
+ (v∇)v = − − ∇Φ,
∂t ρ
∂ρ
+ ∇(ρv) = 0, (1)
∂t
!
d p
= 0,
dt ργ
where Φ is an external gravitational potential; The last equation represents the conservation of
entropy in the particle of fluid with the perfect gas equation of state - p = K eS/cv ργ , with γ = cp /cv
- adiabatic index, S cv - specific entropy and heat capacity for constant volume.
The linearization procedure gives the equations for small perturbations:
∂δv ∇δp δρ ∇p
+ (v∇)δv + (δv∇)v = − + ,
∂t ρ ρ ρ
∂δρ
+ ∇(ρδv) + ∇(δρv) = 0 (2)
∂t
! !
∂ δp δρ p δp δρ p p
γ
−γ γ
+ (v∇) γ − γ γ
+ (δv∇) γ = 0,
∂t ρ ρ ρ ρ ρ ρ ρ
where the third one implies the relation between perturbations of density and pressure.

In our approach all quantities are functions of r, so for Eulerian perturbations with δvz = 0 in
the cylindrical coordinates the system (2) will take the form:
!
∂δvr ∂δvr 1 ∂δp δρ dp
+Ω − 2Ωδvϕ = − − ,
∂t ∂ϕ ρ ∂r ρ dr
∂δvϕ ∂δvϕ 1 d 1 ∂δp
+Ω + (Ωr 2 )δvr = − , (3)
∂t ∂ϕ r dr rρ ∂ϕ
∂δρ ∂δρ 1 ∂ 1 ∂
+Ω + (rρδvr ) + (ρδvϕ ) = 0,
∂t ∂ϕ r ∂r r ∂ϕ

3
 ! ! !
∂ δρ ∂ δρ ∂ p
δp − γp +Ω δp − γp + ργ δvr =0
∂t ρ ∂ϕ ρ ∂r ργ
We’ll be looking for a solution in the form of normal modes:

δv = v̄(r)e−i(ωt−mϕ) ,
δp = p̄(r)e−i(ωt−mϕ) (4)
δρ = ρ̄(r)e−i(ωt−mϕ)
here m - is the azimuthal wavenumber, ω = ωr + i ωi - complex frequency, which for growing modes
has the positive imaginary part, i.e. increment of mode. Possible ω has to be determined from the
solution of the boundary problem coming to the integration of ordinary differential equation with
the boundary conditions in the boundary points r1 and r2 for functions v̄, p̄ ρ̄, which, in general,
are complex.
The substitution (4) into (3) yields:
!
1 dp̄ ρ̄ dp
(−iω + imΩ)v̄r − 2Ωv̄ϕ = − − ,
ρ dr ρ dr
1 d  2 im
(−iω + imΩ)v̄ϕ + Ωr v̄r = − p̄ (5)
r dr rρ
1 d imρ
(−iω + imΩ)ρ̄ + (rρv̄r ) + v̄ϕ = 0,
r dr r
! !
ρ̄ d p
(−iω + imΩ) p̄ − γp + ργ v̄r =0
ρ dr ργ
From the last expression in (5) we’ll get the relationship between p̄ and ρ̄:
ρ̄ 1 p̄ N2
= 2 − ξ¯r , (6)
ρ a ρ gef f
where a2 = γp/ρ - is the squared sound velocity,
GM 1 dp
gef f = (Ω2 − 3
)r = (7)
r ρ dr
- the effective gravitational acceleration and ξ¯r = iv̄r /(ω − mΩ) gives the radial dependence of radial
component of the Lagrangian displacement : ξr (r, t) = ξ¯r (r) e−i(ωt−mϕ) . At last,
!
2 1 dρ gef f
N = gef f − 2
ρ dr a
- is the squared buoyancy, or Brunt-Vaasala, frequency which is the characteristic internal oscillation
frequency in medium at rest due to the entropy gradient and gravitational force. Clearly, for constant
entropy dp = a2 dρ, N 2 = 0, and the relation between the Eulerian perturbations of density and
pressure comes to the expression one has for homoenropic configuration, when the perturbed flow as
well as the unperturbed one are barotropic: p̄ = a2 ρ̄.
Then, the equation of state p = Ke(γ−1) s ργ , where s - entropy in the units of the universal gas
constant ℜ and squared sound velocity a2 = γp/ρ, gives the following differential relations:
da2 γ − 1 ds
= (γ − 1) gef f + a2 , (8)
dr γ dr

4
 1 dρ gef f γ − 1 ds
= 2 −
ρ dr a γ dr
(8) shows that squared buoyancy frequency is simply
γ −1 ds
N2 = −gef f , (9)
γ dr
So for stratified medium the eq. (6) is modified just by the additional term proportional to the
entropy gradient. Finally, from the system (5) we obtain the equation for infinitesimal perturbations
with respect to p̄/ρ:
2Ωρ N 2
!! " !
d rρ d p̄ 2m d Ωρ
 
− − +
dr D dr ρ ω − mΩ dr D D gef f
 !2  !
2 2 2 2
!
m N 1 1 N d rρ N p̄
ρ 1− +rρ  + −  = 0, (10)
Dr (ω − mΩ)2 a2 D gef f dr D gef f ρ

where D = Dh + N 2 , Dh = κ2 − (ω − mΩ)2 - the quantity signed as ”D” for homoentropic flow

in ZS2. The squared epicyclic frequency is the same as before:
2Ω d  2 
κ2 = Ωr .
r dr
Excluding the difference in notations, the main equation (10) coincides with one used by Lovelace
et al. (1999) for perturbations integrated in the z-direction. For N 2 = 0 the eq. (10) comes to the
equation for the Eulerian perturbation of enthalpy in homoentropic fluid used in ZS2. In the limit
γ → ∞ the equation (10) determines the behavior of small perturbations in stratified incompressible
flow. At the same time the entropy gradient is connected directly with the variable density:
ds 1 dρ
=− , (11)
dr ρ dr
and the squared buoyancy frequency is:
1 dρ
N 2 = gef f , (12)
ρ dr

At last, we write here the expressions for v̄r and v̄ϕ , in terms of p̄/ρ:
N2
" ! ! !#
i d p̄ 2mΩ p̄
v̄r = (ω − mΩ) − − (ω − mΩ) ,
D dr ρ r gef f ρ

1 κ2 N2 N2
" ! !! ! !#
d p̄ p̄ m p̄
v̄ϕ = + − ω − mΩ −
D 2Ω dr ρ gef f ρ r ω − mΩ ρ

2.2 The boundary conditions

We have to impose boundary condition for p̄/ρ. At the free boundary for perturbed as well as for
unperturbed flow h = a2 = p = 0. Here h - is an enthalpy of perfect gas. This implies vanishing of
the Lagrangian perturbation of enthalpy at the boundary of the unperturbed flow ∆h|b = 0. Further,

5
for the fixed fluid particle dh = dp/ρ + T dS = dp/ρ since the motion is isentropic. Consequently,
∆h = ∆p/ρ = 0 at the free boundary. Using the relation between small Eulerian and Lagrangian
quantities (look Tassoul, 1978) we have:
δp 1 dp
+ ξr | =0 (13)
ρ ρ dr
The eq. (13) yields the boundary condition for p̄/ρ. Here it’s necessary to use the relation between
ξ¯r and v̄r have been used to obtain (6) and expression for v̄r which was written at the end of the last
section:
! !
d p̄ gef f 2mΩ p̄
 
gef f − Dh + r1 ,r2 =0 (14)
dr ρ r ω − mΩ ρ
One can check that (14) is exactly the same as the boundary condition for the homoentropic flow
used in ZS2. In particular, the eq. (14) doesn’t contain terms with variable entropy.

3 Angular velocity and entropy profiles

The rotation curve is determined by gravitation and the pressure gradient. In the previous research
(ZS1 and ZS2) we used the two profiles of angular velocity of the basic flow:
!1/2
r K r − r1
 −3 
Ω(r) = Ω0 + sin 2π (15)
r0 r/r0 r2 − r1
r −q
 
Ω(r) = Ω0 (16)
r0
here r1 and r2 - are the boundary points of the flow, r0 - point where the matter rotates with the
Keplerian frequency Ω0 . The coefficient K sets the deviation of Ω(r) from the Keplerian profile. All
frequencies will be in units of Ω0 and all the distances - will be given in units of r0 below. To define
the radial extent of the flow we’ll use the quantity w = (r2 − r1 )/r0 × 100%.
The power law (16) has been traditionally used in the research of axisymmetric flow instabilities.
This profile is characterized by the maximum of gef f at the edge of the flow for q > 3/2. The rotation
according to (15) implies that gef f obeys the sinusoidal law and vanishes at the boundaries. In ZS2 it
was showed that gef f → 0 in r1 and r2 entails the absence of growing sound modes. For consistency
we take here the same rotation profiles although the main attention will be paid to the power law.
To solve the equation (10) one should calculate first a2 (r). Integrating the first expression in (8)
we have: Z 
a2 = e(γ−1) s/γ (γ − 1) gef f e− (γ−1) s/γ dr + C1 , (17)

where gef f is defined by the expr. (7), and C1 for specific w is determined from the condition that a2
must vanish at r1 and r2 . Let’s notice that for profile(15) the boundary points are defined beforehand,
since r0 = (r2 + r1 )/2. For (16) r1 and r2 must be calculated together with the constant C1 .
For the certain parameters a2 (r) was determined by the Chebyshev polynomials approximation of
the tabular function obtained by the numerical integration of the equation (17). Note that we could
avoid difficulties with a2 calculation using an appropriate angular velocity profile. But in this case
we would miss the possibility to compare the new results with the previous ones found in SZ1 and
ZS2 in the limit of uniform entropy distribution.

6
 Besides, we’d lose a chance to explore the effects induced by stratification purely, i. e. with the
same distribution of effective gravity and vorticity in the basic flow.
In the limit of γ → ∞ one doesn’t have to find explicit dependencies of the flow quantities (such
as p or ρ) on r to solve (10). However, in case of power law rotation profile we need to determine
the boundary points, what can be done using the stationarity condition. Since gef f = ρ1 dpdr
, we have:
Z
p= gef f ρ dr + C2 (18)

To integrate the equation (18), one should know ρ(r). For the incompressible stratified medium we
get ρ(r) by integration of the equation (11) with the certain entropy profile.

Entropy profile
The entropy distribution was defined as follows:
s(r) = s0 − s1 (r/r0 − 1)2 , (19)
where as before, r0 - is the point in the flow with Keplerian rotation, so that gef f = 0, s0 , s1 - are the
constants. This form is convenient because when s1 6= 0 the entropy gradient is directed in the same
way with respect to gef f everywhere in (r1 , r2 ), in the other words N 2 has the same sign going to
the zero in the vicinity of r0 . Moreover, the function (19) is monotonous on the each side of r0 , what
allows to set entropy before calculation of the boundary points necessary for the power law rotation
profile. Note that for the fixed w the result of calculations doesn’t depend on s0 . Indeed, adjusted
for the dep. (19) squared sound speed (17) is:
a2 (r) = F1 (r) + L(s0 , C1 ) F2 (r)
Now, if in the certain boundary point a2 = 0, we find L(s0 , C1 ). It means, that s0 controls only C1 ,
and leaves a2 to be unchanged. Then, since the equation (10) contains only the derivatives of s(r),
we concluded that eigenvalues of the boundary problem are not affected by the value of s0 . So s1 is
the only parameter that characterizes the stratification of the flow.

Opposite to non-axisymmetric an axisymmetric perturbations obey the Heiland’s criterium (Tas-

soul, 1978). Namely, for considered stationary flow the necessary and sufficient condition of it’s
stability to perturbations with axial symmetry is the inequality:
κ2 + N 2 > 0 (20)
Thus the stratification with s1 < 0 contributes to the stabilization of the flow (to axisymmetric
perturbations), and vice versa, with s1 > 0 - to the instability in comparison with the homoentropic
case. For Ω(r) ∝ r −q with q < 2 and s1 > 0 the flow with sufficiently small w will always satisfy to
Heiland’s criterium. However, the increase of w will cause the break of (20) close to the boundaries,
since the entropy gradient goes up while moving away from r0 .
One more criterium is valid for stratified flows. It is the Richardson criterium. We’ve mentioned
the authors who generalized it for baroclinic axisymmetric flows. For barotropic configuration we
consider here it yields that
N2 1
Ri = 2 > (21)
(rdΩ/dr) 4
everywhere in (r1 , r2 ) is the sufficient condition of the stability (see also Glatzel, 1990). Ri - is the
Richardson number. It’s clear that (21) is never fulfilled here, since for any s1 Ri → 0 while r → r0 .

7
4 Numerical calculation
We’ll solve the boundary problem consisting of the equation (10) and boundary conditions (14) at the
interval (r1 , r2 ) for compressible and incompressible fluid. Let’s rewrite (10) in the form convenient
for integration: !′′ !′
N2
!
p̄ 1 D′ gef f p̄
+ − + + 2 −
ρ r D gef f a ρ
N2 m2 N2
" !! !
2m Ω′ Ω gef f D′
+ − − + 2 1− + (22)
ω − mΩ r r a2 D gef f r (ω − mΩ)2
!′
N2 N2
!# !
D 1 D ′ gef f p̄
− − − + 2 =0
a2 gef f gef f r D a ρ
with the finite compressibility and
!′′ !′
N2 N2
! " !!
p̄ 1 D′ p̄ 2m Ω′ Ω D′
+ − + − − + +
ρ r D gef f ρ ω − mΩ r r D gef f
!′
m2 N2 N2 N2
! !# !
1 D′ p̄
1 − − − − =0 (23)
r2 (ω − mΩ)2 gef f gef f r D ρ
in the incompressible limit. In the equations (22) and (23) the derivative on r is marked by a stroke.
Since we’re interested in the growing modes the integration can be implemented along the real axis
according to the Lin’s rule (1945). Each of the equations (22) and (23) have been separated into
the real and imaginary parts, so that the system of four first order equations consisting of the real
quantities have been integrated.
For ωi 6= 0 the equation (23) doesn’t contain singular points at real axis so the boundary problem
have been solved as in ZS1. The four specified linearly independent vectors have been used as initial
conditions in r1 to obtain the corresponding solutions in r2 . The boundary conditions in r1 and r2
give then four algebraic equations and the eigenvalues should be obtained by setting the determinant
of corresponding matrix equal to zero. The finite compressibility involves the sound speed to vanish
at the boundaries. The main equation contains now singularities in r1 and r2 (for power law rotation
it’s the first order poles and for (15) profile it’s the second order poles) so the numerical algorithm
have to be adjusted appropriately.
As in ZS2 the solution should be expanded to the generalized series in the vicinity of r1 and r2 :
∞
p̄
= (r − r1,2 )µ ai (r − r1,2 )i
X
(24)
ρ i=0

Substituting (24) into (22) and multiplying it by (r − r1,2 )2 one gets the recurrent expressions for
ai and the squared equation for µ (see ZS2). The last one defines two independent solutions of (22)
close to the boundaries with some µ1 and µ2 . The regularity condition assigns what solution should
be chosen. After that the numerical integration from left and right boundary gives two solutions in
some point inside (r1 , r2 ). The coupling condition can be rewritten as the determinant of the certain
matrix equal to the zero as above and again gives an eigenvalues. At that µ1 and µ2 turn out to be
independent on the entropy distribution in spite of the fact that the recurrent expressions for ai differ
from that in ZS2. For the power law rotation µ = 0 as before corresponds to the regular solution
and the relation between a0 (value of p̄/ρ at the boundary) and a1 (value of (p̄/ρ)′ at the boundary)
as before coincides with the boundary condition (14) (As it’s shown in ZS2 for the uniform entropy
distribution).

8
 For rotation with quasi-sine deviation from the Keplerian law we choose again the solution with
Re(µ) > 1, what simultaneously implies that it satisfies to the boundary condition when gef f |r1,2 = 0.
It’s appropriate to mention here one popular way when the stratification of basic flow is set by
the relation p ∝ ρΓ , where Γ 6= γ (see, for instance, the paper by Kojima, 1989). This approach
can’t be implemented here since at the boundaries r → r1,2 , ρ → 0, and entropy and its derivatives
are not finite in r1,2 . So the coefficients of (22) have the singularities also because of the terms that
contain entropy gradients. Then it’s easy to demonstrate that, at least for the power law rotation,
the relation between a0 and a1 with µ = 0 won’t be equivalent for (14), so that the boundary problem
doesn’t have the solution. The boundary condition is equivalent to the regularity condition only if
s(r) with it’s derivatives is finite at the boundary.

5 The results
5.1 Instability of flow with the power law rotation profile.
There are growing sound and surface gravity modes in axisymmetric flow with free boundaries (Pa-
paloisou & Pringle (1984, 1985, 1987), Blaes & Glatzel (1986), Goldreigh et al. (1986), Glatzel
(1987a,b), Sekyia & Miyama (1988), Kojima (1989) and others). An entropy distribution is sup-
posed to modify mentioned branches of instability and possibly the growing internal gravity modes
should appear (Glatzel, 1990, 1991). The results of calculations will be presented here mostly at the
plots, where the curves of increments will be drawn. We also won’t limit ourselves with the region
stable to axisymmetric perturbations according to the Heiland’s criterium. In the curves of ωi (s1 )
the transition from one region to another will be marked by a stroke.

Calculations with finite compressibility In fig. 1 ωi (s1 ) and ωr (s1 ) are displayed for the growing
surface gravity mode with m = 10, which form so called principal or incompressible branch of
instability, since ωi > 0 exist also in the limit γ → ∞. In the absence of stratification for q < 2.0
and γ < ∞ an increment of principal branch disappears at some maximum value of w. As we can
see for small values of w an increment monotonously increases while s1 gets higher, i.e. while the
configuration becomes less stable according to the Heiland’s criterium (20). However, for greater w
minimum appears in the increment curve and then the growing incompressible mode disappears at
certain s1 > 0. For even greater w surface gravity mode stops growing in the flow with s = const, but
at the same time increment emerges both for the positive (s1 < 0), and the negative (s1 > 0) entropy
gradient. The last result is consistent with the calculations of Kojima et al. (1989), who made the
conclusion that entropy growth in the opposite direction to effective gravity (here it’s means s1 < 0)
contributes to widening of instability range in the values of w. The similar effect gives the calculation
of ωi (s1 ) for m = 1 and different values of q (look at fig. 2). One can see that while the vorticity
gradient is increasing (i.e. q is decreasing) stabilization sets first in the homoentropic case. Beside
that, in fig.2 dashed lines denote the additional branches of instability. As it’ll be clear below, these
are growing internal gravity modes caused by stratification. However, these modes were found in the
present work only in the region unstable with respect to axisymmetric perturbations.
In a more intricated manner stratification affects the growing sound modes. All the branches of
sound modes increments for specific w are depicted in fig. 3. We’ve used here the results of ZS2,
where the detailed calculation of sonic instability was implemented, exactly, the calculation of ωi (w),
with m = 10 and q = 1.58. We’ve used the value of w corresponding to the modes coupling for one
of the increments. In the homoentropic flow the modes coupling entails the strong increase of ωi in a
narrow range of w. In fig. 3 the modes coupling with s1 = 0 occurs for branch (5). To the right side
from the vertical dashed line the flow is unstable according to the Heiland’s criterium. The increment

9
 (a)

0.08

5
0.07

0.06
7
0.05
ωi
0.04 9
8
q=1.8
0.03
γ=5/3
m=10 10
0.02
9
0.01 15
10

0
-600 -400 -200 0 200 400 600
S1
(b)

10
5

7
9.99

9.98 8
10

9.97
9

9.96
q=1.8
γ=5/3
9.95
ωr m=10

9.9

9.88 9
9.86

9.84

9.82 10
9.8
-600 -400 -200 0 200 400 600
S1

Figure 1: Increment (pic. (a)) and pattern speed (pic. (b)) of non-axisymmetric modes of perturbations as a function
of entropy gradient for the power law rotation profile. The numbers are the values of w in percents. The solid lines
denote the growing surface gravity modes, the dashed lines - the sound modes. The strokes in the curves mark the
edge of stability to axisymmetric perturbations. For curves corresponding to w = 9%, 10%, 15%, with s1 > 0 the flow
is unstable to axisymmetric perturbations.

10
 (a)

0.5

w=70%
γ=5/3
0.4
m=1

ωi 0.3
1.75
1.7
2.0
0.2 1.77
1.9

1.9

2.0
0.1 1.8
1.75 1.8

1.6
0
-60 -40 -20 0 20 40 60
S1

(b)

1.75
1.77 w=70%
0.95 2.0 γ=5/3
1.9
1.8 m=1

ωr
0.9 1.6

0.85
1.75

1.7

0.8
-60 -40 -20 0 20 40 60
S1

Figure 2: The same as in fig. 1, but for fixed w and the dashed lines denote now the growing internal gravity modes.
The numbers are the values of q. For curves with q = 1.6, 1.7 the flow in unstable to axisymmetric perturbations. The
curves of pattern speed of internal modes are not depicted because ωr > 1.

11
branches fall into three categories: (1),(8),(2),(7),(3) - the perturbations growth disappears both at
some s1 < 0 and at some s1 > 0; (4),(6) - increment gradually increases with s1 ; (5) - corresponds
to the modes coupling and ωi (s1 ) behaves in the most tangled way. Exactly, at some value of
s1 = scr1 > 0 increment (5) vanishes, but comes back at once and grows then with s1 . The same
features of modes coupling arise in dependence ωi (w) for homoentropic flow presented in ZS2. To
check if the presence of s = scr 1 , where ωi → 0 is the common peculiarity for modes coupling, we’ve
calculated few ωi (s1 ) for different w. In can be seen that all determined ωi (s1 ), have a break point
in increment curve (fig. 4).
At last, branches of sonic instability for m = 1 one can see in fig. 5. The calculations here were
mostly carried out in area unstable to axisymmetric perturbations. It was revealed that both at
some s1 < 0 and at some s1 > 0 increment vanishes. In ZS2 we discussed that the growing sound
modes with m = 1 arise as a result of two mechanisms: resonant amplification by the basic flow in
critical layer (mechanism Landau) and coupling with decaying surface gravity mode, what sets it
apart from the above instability being the consequence of coupling of two sound modes. Clearly, the
increment behavior differs from what we’ve seen in fig. 4. The main is that in this case for none of
the branches increment arises again after vanishing at some s1 > scr 1 .
While calculating the presented curves we tried to find the growing internal gravity modes in area
stable to axisymmetric perturbations. However we failed to do that and the only one instability that
arises due to stratification is denoted by dashed lines in fig. 2. We’ve made an attempt then to find
growing internal gravity modes in the incompressible limit, since this branch of instability must be
insensitive to the compressibility.

Calculations in incompressible limit An incompressible limit was used here first of all as an additional
way to check the results presented so far, since the numerical algorithm differs considerably in that
case (look above). In fig.6 we put increment curves with dependence on γ (pics a & b) and with
dependence on s1 ( pics c & d) calculated in approximation concerned here. For γ = 5/3 the
corresponding values of ωi and ωr can be found in fig. 2. The decrease of sound speed suppresses an
internal gravity modes growth (dashed curves) and its increment vanishes at some γ > 1. On the
contrary, the surface gravity modes are unstable right up to γ → 1. Then, as it must be, while s1 → 0
ωi of surface mode diminishes up to it’s value in homogenerous incompressible fluid (Jaroszynski,
1988, SZ1). Concerning internal mode it’s growth disappears at some s1 6= 0.
The parameters in fig.6 correspond to unstable axisymmetric perturbations. We’ve already noticed
that we missed to find growing internal gravity modes in the flow with free boundaries and γ <
∞ that is stable according to Heiland’s criterium. The numerical analysis confirmed this also in
incompressible limit. However, it turned out that if someone impose rigid boundaries the increments
of modes that grow with free boundaries increase significantly exceeding Ω0 .
Moreover, rigid boundaries involve the instability of configuration stable according to Heiland’s
criterium (20). Note that Glatzel (1990) also have found growing internal modes in a problem
with rigid boundaries. At a glance, such influence of boundaries on the unstable internal modes is
embarrassing and the stratification must play the main role. However one should remember that
we consider the amplification of internal oscillations by the shear flow. The neutral internal modes
certainly exist independently on the kind of the boundaries. But the ability to grow (to decay) for
example owing to interaction with the basic flow is controlled not only by a vorticity gradient at
the critical layer (as it’s for homoentropic flow) but additionally by the shape of perturbations field
itself as it was shown, in particular, by Troitskaya & Fabrikant (1989) (see also the monograph by
Stepanyants & Fabrikant, 1996). The latter is confirmed also by Glatzel’s paper (1990), where he
found that the internal modes grow (decay) outside the regions of coupling what can happen only
if perturbations resonantly interact with the basic flow. But he used the simplified angular velocity

12
 (a)
0.01
q=1.58
0.009 γ=5/3
m=10
0.008 w=24.2% 6

0.007
7
4
ωi 0.006

0.005 5

0.004

0.003

0.002 1
3
8 2
0.001

0
-100 0 100 200 300 400 500 600 700
S1

(b)

11.5 1
2

11 3

4
10.5
ωr 5
q=1.58
10 γ=5/3
m=10
6 w=24.2%
9.5
7

9
8

8.5
-100 0 100 200 300 400 500 600 700
S1

Figure 3: The same as in the fig. 1, but for fixed w and the increment and pattern speed of sound modes only are
presented. The numerals are the numbers of instability branches. The vertical line marks the limit of stability to
axisymmetric perturbations.

13
 (a)
0.008
q=1.58
0.007 γ=5/3
m=10

0.006

0.005
ωi 38.9

0.004
40.2
0.003

0.002 33.15

0.001
27.7
0
-150 -100 -50 0 50 100 150
S1

(b)
12
q=1.58
11.5 γ=5/3 40.2
m=10
11

27.7
10.5
ωr
10
33.15
9.5

8.5
38.9

8
-150 -100 -50 0 50 100 150
S1

Figure 4: The same as in the fig. 3 , but for different w, signed each curve. The branches presented are the
result of sound modes coupling (for s1 = 0). The strokes in the curves denote the edge of stability to axisymmetric
perturbations. For w = 33.15%, 27.7% the flow is stable to axisymmetric perturbations in the whole range of s1 .

14
 (a)

0.1 q=1.58
5
γ=5/3
m=1
w=150% 4

0.01
3
ωi
2
0.001

0.0001

-5 0 5 10 15 20 25 30
S1

(b)

2.5 q=1.58
γ=5/3
m=1
1 w=150%

2 2

ωr 3

4
1.5

-5 0 5 10 15 20 25 30
S1

Figure 5: The same as in the fig. 3 but for m = 1.

15
 ωi
(a) (c)
0.5 0.5

0.4 0.4

0.3 0.3
w=70%
q=1.9
0.2 m=1 0.2
S 1=50 w=70%
q=1.9
0.1 0.1 m=1

0 0
1 10 100 1000 50 40 30 20 10 0
γ S1
1.1 1.1
ωr (b) (d)
1 1

0.9 0.9

0.8 0.8

0.7 0.7
w=70% w=70%
q=1.9 q=1.9
0.6 m=1 0.6 m=1
S 1 =50

0.5 0.5
1 10 100 1000 50 40 30 20 10 0
γ S1

Figure 6: In pics. (a) and (b) - the dependences ωi (γ) and ωr (γ), in pics. (c) and (d) - the dependences ωi (s1 ) and
ωr (s1 ). The solid lines denote the growing surface gravity modes, the dashed lines denote the growing internal gravity
modes.

16
 0.01

0.009 m=10
eq. 1.58
γ=2.5
0.008

0.007 5

0.006
ωi
0.005

0.004 8.5
7
0.003 8.2

0.002
8
0.001 9

0
-700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700
S1

Figure 7: Increment of the surface gravity mode as a function of entropy gradient for the Keplerian rotation with
quasi-sine deviation. The angular velocity profile corresponds to the profile of Ω(r) in the case of power law rotation
with the same w (look comments in text). In the marked range of s1 the flow is stable to axisymmetric perturbations.

profile corresponding to the constant vorticity in the flow, what excludes the interaction of small
perturbations with the flow.
The perturbations profile which is responsible for its growth along with the vorticity gradient in
stratified fluid is the solution of the boundary problem, so that the amplification of internal modes
is determined also by the boundary conditions.

5.2 Instability of rotation with quasi-sine deviation from the Keplerian profile
The instability of flow with such profile was investigated in the previous papers ZS1 and ZS2. Par-
ticularly, it was revealed that the stabilization always occurs at some least K > 0 (look the formula
(15) ), and for γ < ∞ there’s no sonic instability. The last implies that in non-stratified flow only
surface gravity modes can grow. Accordingly, in the fig.7 we can see how this sort of perturbations is
modified by the stratification. The angular velocity profile is set in order to make equal the maximum
enthalpy in the flow with s1 = 0 to maximum enthalpy in the flow with the same radial size w but
rotating with a power law profile with q = 1.58 (see ZS2). Note that calculations were carried out for
γ = 2.5, since the flow turns out to be stable for ordinary γ = 5/3. Several curves ωi (s1 ) correspond
to the different values of w. The dependence on s1 is trivial: ωi increases while the flow becomes
less stable according to the Heiland’s criterium (20). Similar to the power law rotation the negative
entropy gradient (s1 > 0) widens the range of w when the growing surface gravity modes exist.

6 Conclusions
Well, we investigated in the two-dimensional approximation the instability of the laminar axisym-
metric stratified flow with free boundaries to infinitesimal non-axisymmetric perturbations. In the
case of the power law rotation profile the essential influence of entropy gradient on the growing sound
and surface gravity modes was revealed. With both negative and positive sign of entropy gradient
relative to gef f direction the range of w and q where growing surface gravity modes exist becomes

17
wider. Specifically, its increment vanishes at greater w and smaller q. For sound modes the incre-
ments fall into three categories. In the first one ωi gradually increases while the flow becomes less
stable according to the Heiland’s criterium (i.e while s1 gets higher). In the second case increment
vanishes both at some positive (s1 < 0) and negative (s1 > 0) entropy gradient. At last, the third
category contains increments being the result of modes coupling. Corresponding ωi (s1 ) has a break
in scr
1 , where ωi → 0. Similar peculiarities in increment behavior has been revealed also in ZS2 for
ωi (w) which also had the points ωi → 0.
Besides the sound and surface gravity modes we tried to find growing internal gravity modes. The
latter was done as well in incompressible limit. It turned out that mentioned sort of perturbations
exists only in the flow that is unstable according to the Heiland’s criterium, i.e. unstable to axisym-
metric perturbations. However, with rigid boundaries the growing internal modes where found also
in the stable to axisymmetric perturbations region (see also Glatzel, 1990). We suppose such crucial
role of the boundaries to be the consequence of specific character of amplification in a stratified
medium. In particular, the transition of energy to perturbations from the basic flow is determined
now not only by a vorticity profile in the critical layer but also by the profile of perturbations, which
in turn depends on the boundaries (Troitskaya & Fabrikant, 1989).
In the case of the Keplerian rotation with quasi-sine deviation we found that the increment of the
surface gravity mode increases gradually while s1 gets higher.
We also discussed in the paper that it’s incorrect to set the stratification in a problem with free
boundaries by a polytropic law with the index different from the adiabatic value, since in this case the
perturbations don’t satisfy the boundary conditions. The last is the consequence of infinite entropy
and it’s derivatives in the boundary points. It reveals that the condition of perturbations regularity
at the boundaries is equivalent to the boundary condition (i.e. ∆p = 0) only if entropy and it’s
derivatives have finite values at the boundaries.
This paper was supported by grant RFFI-NNIO 06-02-16025.

References
[1] Blaes O. M., Glatzel W., MNRAS, v.220, p.253 (1986).
[2] Glatzel W., MNRAS, v.225, p.227 (1987a).
[3] Glatzel W., MNRAS, v.228, p.77 (1987b).
[4] Glatzel W., MNRAS, v.242, p.338 (1990).
[5] Glatzel W., Rev. in Modern Astron., v.4, p.104 (1991).
[6] Goldreigh P., Goodman J., Narayan R., MNRAS, v.221, p.339 (1986).
[7] Ghosh P., Abramowicz M. A., ApJ, v.366, p.221 (1991).
[8] Zhuravlev V., Shakura N., AstL, v.33, n.8, p. 536 (2007) .
[9] Zhuravlev V., Shakura N., AstrL, v. (2007).
[10] Klahr H. H., Bodenheimer P., ApJ, v.582, p.869 (2003).
[11] Kojima Y., MNRAS, v.236, p.589 (1989).
[12] Kojima Y., Miyama S. M., Kubotani H., v.238, p.753 (1989).
[13] Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., ApJ, v.513, p.805 (1999).

18
[14] Li H., Finn J.M., Lovelace R.V.E., Colgate S. A., ApJ, v.533, p.1023 (2000).
[15] Lin C. C., Q. Appl. Math.,, v.3, p.117 (1945).
[16] Tassoul J.-L., Theory of Rotating Stars, Princeton University Press (1978).
[17] Papaloizou J.C.B., Lin D.N.C., Annu. Rev. Astron. Astrophys., v.33, p.505 (1995)
[18] Papaloizou J.C.B., Pringle J.E., MNRAS, v.208, p.721 (1984).
[19] Papaloizou J.C.B., Pringle J.E., MNRAS, v.213, p.799 (1985).
[20] Papaloizou J.C.B., Pringle J.E., MNRAS, v.225, p.267 (1987).
[21] Sekiya M., Miyama S., MNRAS, v.234, p.107 (1988).
[22] Stepanyants Y.A., Fabrikant A.L. Rasprostranenie Voln v Sdvigovykh Potokakh (in Russian),
Moscow, Fizmatlit (1996).
[23] Troitskaya Y.I., Fabrikant A.L., izv. vuzov: Radiofizika (in Russian), v.32, no.10, p.1221 (1989).
[24] Frank J., Robertson J.A., MNRAS, v.232, p.1 (1988)
[25] Fujimoto M.Y., Astron. Astrophys., v.176, p.53 (1987).
[26] Hanawa T., Astron. Astrophys, v.179, p.383 (1987).
[27] Howard L.N., J. Fluid Mech., v.10, p.496 (1961).
[28] Jaroszynski M., Acta Astronomica, v.38, p.289 (1988).

19