Review Paper
Stabilizing and Destabilizing Effect of Shear Flow
— Beyond Kelvin-Helmholtz Instability —TATSUNO Tomoya*
Graduate School of Frontier Sciences, The University of Tokyo, Tokyo 113-0033, Japan(Received 27 January 2003 / Accepted 11 August 2003)

Abstract
Recent progress in linear spectral studies of incompressible plasma (fluid) with flows is reviewed.
Mathematically, non-Hermiticity of the generating operator brings about incompleteness of theeigenvalue analysis and unremovable coupling between modes. This may lead to secular growth(algebraic in time) even if all eigenvalues exhibit stable oscillations. Physically, the stretching effect ofshear flow is considered a stabilizing mechanism for various plasma instabilities; however, shear flow
may equally effectively destabilize a class of plasma modes. The condition for stabilization/
destabilization is discussed as well as its physical mechanisms.
Keywords:
linear stability, shear flow, stretching effect, Non-Hermiticity (non-self-adjointness), nilpotent, KelvinHelmholtz instability, Alfvén wave

1. Introduction
Shear flows in plasmas are attracting much interest
because of their suppression effect on turbulence inmagnetically confined fusion devices. The current naïve
suggestion regarding this phenomena holds that thestretching of modes in a shear flow brings about lengthscale reduction leading to the suppression offluctuations. This argument for stability, however,ignores the fact that the available free energy associatedwith a shear flow may be a potent source for thedestabilization of some other class of fluctuations. The
Kelvin-Helmholtz (KH) instability [1-3], for instance, isa well-known example of an instability that feeds on theambient flow-energy. On the other hand, various
discrepancies between theory and experiment regarding
the stability limit of neutral fluids have been reported[4]. This is because rigorous treatment of the shear flow
effects encounters a fatal difficulty arising from the nonHermitian (non-self-adjoint) properties of the problem.
The standard normal mode approach breaks down, andthe theory may fail to provide correct predictions ofevolution even if fluctuation remains in the linear
regime.
Since the stability of an equilibrium has to be
discussed for arbitrary initial conditions, it is importantto find a general solution of the evolution equation
governing linear dynamics. Fortunately, the linearized
equation for a static (flow-less) equilibrium is describedby a Hermitian (self-adjoint) operator when we
determine the evolution equation in terms of thedisplacement vector [5-7]. Thus, as we normally do inconventional eigenvalue analysis, we may simply

author’s e-mail: tatsuno@glue.umd.edu
This article is based on the invited talk at the 19th Annual Meeting of JSPF (Nov. 2002, Inuyama).
* Present address: Institute for Research in Electronics and Applied Physics, University of Maryland, College Park, Maryland
20742-3511, USA

1169

J. Plasma Fusion Res. Vol.79, No.11 (2003) 1169 - 1187

Journal of Plasma and Fusion Research Vol.79, No.11

replace ∂t by –iω and see if the equilibrium may yield anexponentially unstable eigenvalue and corresponding
eigenfunction. We may simply conclude that theequilibrium is unstable if we find at least one eigenvalue
having a positive imaginary part of ω and if thecorresponding eigenfunction has finite energy (nondivergent). The completeness of this spectral method is
guaranteed only for Hermitian operator due to von
Neumann’s theorem [8,9]. However, the question on thestability problem is not answered yet when we do notfind any unstable eigenvalue.
As we already mentioned, the problem is related tothe completeness of the solutions. While it is wellknown that complete solutions to nonlinear equations
are difficult to achieve, the difficulty of solving linear
ones is less recognized. Even if we suppose that all
eigenvalues are obtained, it’s still very difficult toconfirm if the obtained eigenfunctions sufficiently
describe complete dynamics. Initially, the number of
eigenvalues must be infinity for infinite-dimensional
operators (eg., integral or differential operators), andworse, for the non-Hermitian operator, theeigenfunctions are not orthogonal to each other. In thecase of finite-dimensional operators (eg., matrix
operators), the Jordan matrix is the known example inwhich eigenvectors are not sufficient to describe
complete dynamics. It is known that secular growth(algebraic in time) may lead to instability even if all
eigenvalues show stability [10]. Recently such an
algebraic instability has been pointed out to occur in theinfinite-dimensional system due to the frequency
overlapping between continuous spectra [11], and we
may even observe an instability with fractional power intime [12]. We simply do not know what may happen in
the infinite-dimensional non-Hermitian operators!
There are also many unresolved problems related tothe physical understanding of the linear dynamics ofshear flow plasmas. KH instability is one of the most
famous instabilities in parallel shear flow systems, inwhich the background flow is stable unless its curvaturechanges sign in the domain. However, even the
background shear flow having a linear profile (free from
KH instability) can make the equilibrium more unstable
when combined with an Alfvén wave [13]. While thestabilizing mechanisms of shear flow have beenstressed, relatively little attention has been paid to its
destabilizing effect. A shear flow may equallyeffectively destabilize a class of plasma modes even if
the flow profile itself is stable.
In this paper, we describe both the mathematical

November 2003

and physical aspects of the recent developments inlinear stability theories of shear flow plasmas. First, we
introduce the linearized reduced MHD equation
describing the incompressible motions of plasmas inSec. 2. We briefly review the energy principle andproperties of spectra in static (flow-less) plasmas.
In Sec. 3, we consider the surface-wave model(piece-wise linear shear flow profile) of twodimensional vortex dynamics in neutral fluids. We firstrevise the Rayleigh’s analysis of KH instability [14]. He
did not directly manipulate the evolution equation ofvorticity (Rayleigh equation), but his analysis gives us
deep physical insight into the problem. Next, we briefly
review the eigenvalue analysis of the Rayleigh equationand discuss how each term affects the system bycomparison to Rayleigh’s analysis. We also show theexistence of degenerate frequencies (nilpotent) for a
KH-stable wave number. Since the generator of theRayleigh equation is non-Hermitian, the complete
spectra may contain the infinite-dimension analogue of
the Jordan block of a finite-dimension operator. Finally,
we include the effect of gravity. We show thatnumerical integration shows localized secular behavior
in the vorticity due to coupling of two continuous
spectra [11]. This example shows that degeneratefrequency spectra may give rise to the algebraic
instability in non-Hermitian systems even if all
eigenvalues (point or continuous spectra) are real.
In Sec. 4, we show the stabilizing effect of shear
flow due to its stretching effect by invoking Kelvin’s
method of shearing modes [15]. While this method has
been previously called a ‘nonmodal’ approach in a wide
range of literatures [16], it is shown as a particular caseof a generalized ‘modal’ one [12,17,18]. We discuss thetransient and secular behaviors of interchange
fluctuations in Couette flow in an infinite interval. The
combined effect of shear flow mixing and Alfvén wave
propagation overcomes the instability driving force aftersufficiently long periods, and damps all fluctuations of
the magnetic flux with an inverse power of time. On theother hand, electrostatic perturbations can bedestabilized for a sufficiently strong interchange drive.
The time asymptotic behavior in each case is algebraic(non-exponential).
In Sec. 5, we discuss the destabilizing effect ofshear flow [13]. In contrast to its well-known
stabilization of the low-frequency plasma motions, ashear flow may equally effectively destabilize a class ofplasma modes. We focus on the latter quality of the flowby reviewing an incompressible ideal plasma having
1170

Review Paper

Stabilizing and Destabilizing Effect of Shear Flow

Couette flow in a finite interval. In the presence of theflow shear, the growth rate of the perturbation increases
due to the coupling of the Alfvén wave with a RayleighTaylor type instability drive. Marginally stable modes in
the flow-less equilibrium achieve their maximum
growth rate when the maximum flow velocity becomescomparable to the Alfvén velocity. At larger shear flow
velocities, however, the stabilizing “stretching” effectbecomes dominant and the instability is quenched.

function, respectively.
In the linearized equation of motion, we neglect thevariation of the density in the inertial term, but not in
the gravitational force (Boussinesq approximation [1]).
This derivation is almost parallel to that of Strauss [19]so that we do not describe it in detail here. By taking thecurl of the equation of motion and projecting it on the zaxis, we obtain the vorticity evolution equation;

ρ0 ∂t + Vy ∂y ∆φ – V ″y ∂y φ

2. Reduced Magnetohydrodynamics and its
Basic Properties
2.1 Formulation
In this article, we consider the ideal incompressible
motion of plasmas in the fluid description.We restrict
the geometry to a one-dimensional plane slab andhighlight the fundamental structures from bothmathematical and physical points of view.
In the Cartesian coordinates, we consider theequilibria in which the plasma current density is taken
as J0 = J0 ez , with ez denoting the unit vector in the z
direction, while a strong magnetic field is also applied in
the z direction. We restrict the problem by assuming
V0 = (0, Vy (x), 0).

= µ1 B0 ⋅ ∇∆ψ + ∇ J 0 × e z ⋅ ∇ψ + g ∂y ρ1 ,(3)0

where ∆ = ∂ x2 + ∂ y2 denotes the two dimensional
Laplacian operator. The continuity equation immediately
reads(∂t + Vy ∂y)ρ1 = –ρ0′ ∂y φ.

(4)

The induction equation is decurled to yield(∂t + Vy ∂y)ψ = B0 · ∇φ.

(5)Equations (3)–(5) constitute a closed set of
equations describing the linearized incompressible
motion of the magnetized plasma. Actually, Eqs. (3)–(5)can be obtained by directly replacing g = 2 /R0 in thehigh β reduced MHD equations describing tokamak
plasmas [20], where R0 denotes the major radius of thetoroidal device.
Let us examine the structure of the governing
equations (3)–(5). In the flow-less case (Vy = 0), the left
hand sides of Eqs. (3)–(5) only contain the time
derivative. The first two terms on the right hand side ofEq. (3) and the right hand side of Eq. (5) correspond to
the Alfvén wave. In homogeneous magnetic field, the∇J0 term disappears and we obtain the simple wave
equation in which the phase velocity is given by the
Alfvén velocity vA = k · B0 / µ 0 ρ 0 . The last term of
the right-hand side of Eq. (3) and right-hand side of Eq.
(4) correspond to the effect of the gravity. When ρ0′ > 0,they yield a stable gravity wave propagating in the y
direction, and when ρ 0′ < 0, they yield interchange(Rayleigh-Taylor) instability. All these terms yield a
Hermitian contribution in energy norm (∫ |∇φ |2 dV) as
we show in Sec. 2.3.
On the other hand, when we add finite flow, the left
hand sides of Eqs. (3)–(5) become convective derivatives.
If the flow is shear-less (Vy = const), these contribution
yield simply the frequency shift (Doppler shift);however, once if it is sheared, the problem becomes
non-Hermitian. The Vy ∂ y terms in the convective(1)

This situation is illustrated in Fig. 1.
Under the above conditions, we may assume thatthe perturbed fields come from two-dimensional
incompressible motions and are written asv 1 = ∇ φ × ez , b 1 = ∇ ψ × e z ,T. Tatsuno(2)

where φ and ψ denote the stream function and the flux

FIg. 1 Equilibrium configurations we consider in thispaper.

1171

Journal of Plasma and Fusion Research Vol.79, No.11

derivative are Hermitian in the vorticity norm ∫ |∆φ |2dV), but not in the energy norm. These terms may
appear to add a simple multiplication operator if we
assume a wave number ky; however, the energy norm,
which is suitable for other terms (including Vy″ term),
makes them non-Hermitian [21]. Note that the left handside of Eq. (3) yields the Rayleigh equation which
describes Kelvin-Helmholtz instability when B0 and ρ0′
are absent.

equation governed by a Hermitian generator is obtainedby means of the spectral resolution method. Since theeigenfunctions of Hermitian operators are independentand orthogonal with each other, we can draw information regarding the stability from the eigenvalue analysis.
Thus we may replace ∂t by –iω and the equilibrium isstable if and only if all ω corresponding to theeigenvalues of generator are real.
Moreover, we obtain a conserved quantity from theHermiticity of the generator. Let us symbolically write
Eq. (7) by

2.2 Hermiticity (selfadjointness) and energy
principle for static equilibria
When the plasma equilibrium is static and uniform(Vy = 0, B 0 = const, and ρ 0′ = 0), the energy of theperturbed field in two-dimensions is conserved [see Eq.
(2)];

d
dt

1 ( ρ v 2 ) + 1 b 2 dV = 0 .
1
1
2µ0
2

ρ0 ∂ t2Ψ = F Ψ,

d
dt

(6)

∇φ 2 dV ,ρ0 ∂t Ψ ⏐∂t Ψ + δ W = 0 ,

(10)

where

δW := –〈Ψ |F Ψ 〉,(11)

by taking the scalar products with ∂tΨ on both sides ofEq. (9). The conserved quantity in Eq. (10) is also called
energy, while it is different from the conventional one in
Eq. (6). We also find from Eq. (9) the equality(7)ρ0 d 2Ψ ⏐Ψ = ρ0 ∂t Ψ ⏐∂t Ψ – δ W .
2 dt 2where we still assume the static equilibrium (Vy = 0).
Equation (7) is a second order evolution equation for
vorticity Ψ = –∆φ. When we consider a one-dimensional
equilibrium (equilibrium quantities depend only on x),
we may assume wave numbers (k y and k z) in thehomogeneous directions (y and z directions). The second
and the third terms in Eq. (7) then yield only
multiplications on φ. Since the stream function reads φ =–∆–1Ψ, we find that the generator of Eq. (7) is Hermitian(selfadjoint) under the normΨ ⏐Ψ : = – Ψ ∆–1 Ψ dV =(9)where the operator F is Hermitian and ρ0 is a positive
constant due to Boussinesq approximation. Then, theHermiticity of F immediately yields

When we introduce an inhomogeneity of the ambient
magnetic field, however, the energy of the perturbed
field contained in Eq. (6) is no longer a conserved
quantity. We thus combine Eqs. (3)–(5) and obtainρ0 ∂2t ∆φ= µ1 B0 ⋅ ∇∆ B0 ⋅ ∇φ0
+ ∇ J 0 × e z ⋅ ∇ B0 ⋅ ∇φ – ρ 0′ g ∂2y φ ,November 2003

(12)

Equations (10) and (12) lead to the energy principle[5-7], which states that the equilibrium is linearly stable
if and only ifδW ≥ 0( Ψ)(13)

is satisfied. The sufficiency of the condition (13)
immediately follows from the conservation law (10).
The necessity is drawn from Eq. (12) [23], whichindicates, if we find a normalized function whosebilinear form δ W is negative, the existence of an
eigenfunction of F whose eigenvalue λ = –ω2 is greateror equal to –δW / ρ0 > 0.

(8)

2.3 Spectra of Alfvén waves in static
equilibria
We consider the equilibrium magnetic field

where the bar denotes the complex conjugate.
The mathematical grounds of Hermitian operators
are fairly well developed in the context of quantum
mechanics [22]. Namely, the spectra of a Hermitian
operator consists of points and continua, theireigenvalues are real, and their eigenfunctions span thewhole function space (von Neumann theorem; see Refs.
[8,9]). Therefore, the general solution of an evolution

B = (0, By(x), Bz),

(14)with B z = const. From the homogeneity of theequilibrium quantities in the y and z directions, we may
assume the wave numbers k = (0, k y, kz). Since the1172

Review Paper

T. Tatsuno

Stabilizing and Destabilizing Effect of Shear Floworder to highlight the essential properties of the point
spectra. In this case, Eq. (15) yields

generator is Hermitian as is seen in Sec. 2.2, we may
consider the eigenvalue λ of the generator whenreplacing ∂t by –iω (λ = –ω2). Then, Eq. (7) gives

2

d
dx

2
d φ
– ky 1 + G 2 φ = 0 ,2dxω 2 – k ||ρ′ gφ
ω 2 – ω 2A d – k 2y ω 2 – ω 2A φ + k 2y ρ0 φ = 0 ,
0
dx

(15)

where k || = k · B / |B| (= const) is the wave number
parallel to the ambient magnetic field. For the
eigenfunctions to satisfy the boundary conditions, we
needwhere we have defined F(x) = k · B(x) and ωA(x) = F(x)/ µ0 ρ .
2.3.1 Alfvén continuum
Let us consider the case in which ρ0′ = 0. First we
find a series of singular solutions. Equation (15) should
not contain any singular solution except the Alfvén
resonance (ω2 – ωA2 = 0) [24]. Suppose that ω2 – ωA2(x)has the zero of the order h (∈N) at x = xs, i.e.,ω2 – ωA2(x) = c(x) (x – xs)h,k||2 – G < ω2.

Since> 0, the Alfvén wave acts to stabilize
interchange modes. If Eq. (21) is satisfied, we obtaineven and odd eigenmodes;

(16)

φ=

cos (n π x / 2)
sin (n π x / 2)

(22)

2

ω 2 = k 2|| –

(17)

ky G
2

k y + n2 π 2 / 4

,

(23)

showing that ω 2 decreases monotonically as |k y |increases.
The upper bound of ω 2 is k ||2 , which is theaccumulation point of ωn2 as n → ∞; i.e.,k||2 – G < ω2 < k||2 .

(24)

The distribution of the eigenvalues is illustrated in Fig.
2, which represents a typical spectral structure of thediscrete part of the shear Alfvén branch [29]. If k||2, islarger than the drive G, there is no instability even if G
> 0. The instability condition is given by(18)

where g1(x) and g2(x) are analytic functions with g1(xs)≠ 0. When we apply the energy norm (8) to the solution(18), which now reads

2

φ ⏐φ = – φ ∆φ dV ,for n : odd,
for n : evenrespectively. The eigenmodes contain n – 1 nodes (zerosin φ). The dispersion relation is

We thus find that the point x = xs is a regular singularpoint of the ODE (15) for any h ∈ N. By means of the
Frobenius expansion [25], we obtain a logarithmic
singularity in the solution since two solutions of theindicial equation have an integral difference for any h[26];

φ(x) = a1g1(x) + a2[g1(x) log |x– xs| + g2(x)],

(21)

k ||2

where c(x) is an analytic function with finite value at x =xs. For investigating the behavior of the solution in thevicinity of the singular point xs, we take the leading
order of the Taylor expansion (16) and substitute it into
Eq. (15), which yields
2
2
d φ
dφ
+ h
+ ky φ = 0 .
dx 2 x – x s dx

(20)

G > k || 1 +(19)

we then find that Eq. (18) actually gives a non-squareintegrable solution corresponding to the continuous
spectrum [27]. We can show that Eq. (15) with ρ0′ = 0
has no other spectrum than the Alfvén continuum [21].
2.3.2 Interchange mode
When ρ 0′ ≠ 0 in Eq. (15), the singularity is nolonger regular if the zero of Eq. (16) is not simple (h >1). In this case, we do not have a general explicit
representation of the singular solution. Moreover, thelast term in Eq. (15) admits the point spectra.
Let us briefly represent the case ρ0′, By = const in

π2
2

4k y

.

(25)

In more general cases, the following facts are
known. Due to Sturm’s oscillation theorem [25], we will
have an unstable eigenvalue ω2 < 0 which satisfies theboundary condition on both sides if the solution for ω2 =Fig. 2 Distribution of eigenvalues in ω2 space.

1173

Journal of Plasma and Fusion Research Vol.79, No.11

November 2003

0 satisfying the boundary condition only on one edge
has any node in the domain [28]. Furthermore, thenumber of these point spectra are infinite, whichaccumulate on the edge of the continuum ω2 = inf ωA2[29]. When the smallest eigenvalue is positive and themode has no resonant surface inside the plasma, theeigenfunction shows a global stable oscillation, which issimilar to the global Alfvén eigenmode [30].

3. Surface Wave Model of RayleighEquationIn this section, we focus on the Rayleigh equation∂t ∆φ + Vy ∂y ∆φ – Vy″∂y φ = 0,

Fig. 3 Rayleigh’s piece-wise linear shear flow profile.

(26)

under the velocity profile shown in Fig. 3. We may
derive Eq. (26) from Eq. (3) by neglecting magnetic
field (B0 = 0) and gravity (g = 0). It is noted that theassumption of a strong B0z in deriving Eq. (3) does notcreate a problem because it directly takes into accountthe curl of the original fluid equation under theincompressibility condition.

Fig. 4 Helmholtz’ principle, which gives the integral
kernel of the operator (∇ ×)–1.

3.1 Rayleigh’s analysis
Rayleigh had first introduced the piece-wise linear
shear flow profile for the analysis of KH instability [14].
His discussion was directly related to the surface
displacement which separates the distinct vorticity
regions. By following his analysis, we will discuss thephysical aspects of this issue in this subsection.
Suppose a stationary shear flow profile as
illustrated in Fig. 3. The profile of the velocity field
consists of two constant regions (x < –a and a < x) and asingle linear shear (–a < x < a) region. Thus, thevorticity field is constant σ = U/a in –a < x < a and zero
elsewhere.
Suppose that we deform the boundary surface on x
= a by an infinitesimally small amountξ1( y,t) = H1ei(ky–ω t),

Fig. 5 Surface displacement at (a,η) produces thevelocity field at (a + ζ, y).

(27)

dq = σ d A ,2π r

where H1 denotes the complex amplitude. Since thevorticity is discontinuous around x = a, deformation of
the surface induces perturbed velocity into the system.
By a theorem due to Helmholtz [14], theinfinitesimally small area element dA with the z-directed
vorticity σ induces the transverse velocity
1(28)

where r denotes the distance from the element (see Fig.
4). This corresponds to the integral kernel of theoperator (∇ ×)–1 1. In the case of surface deformation
shown by Eq. (27), we obtain the velocity at (a + ζ,y)induced by the surface element ξ1 dη at (a,η) as

In three dimensions, the operator (∇ ×)–1 is described by v (r ) = 1

4π

Integrating it with respect to z from –∞ to ∞, we obtain Eq. (28).
1174Ω (r ′ ) ×r – r′
r – r′

3dr ′ in the infinite domain.

Review Paper

T. Tatsuno

Stabilizing and Destabilizing Effect of Shear Flowrelationσξ dηdq = 1
,
2π r
2

2

(29)

ω2 = σ422where r = (y – η) + b . Taking the x component of thevelocity dq yields

dv 1x =

H1 σ (η – y ) eω t – iω tη–ye dη .
r dq = 2π
r2=

H1 σ – iω te2π(η – y ) e i kηdη .
r2

i H1 σ – k ξ i (ky – ω t )
e
e
,
2

(40)

The unstable branch of Eq. (40) isω =i σ2(30)

e–4ka – 2ka – 1

2

(41)

,

when |(2ka – 1)e2ka | < 1. The corresponding eigenstate
has the relative amplitude of the surface displacement

Integrating Eq. (30) with respect to η yields the velocity
field at (a + ζ , y) induced by the whole surface
displacement on x = av 1x =

22ka – 1 –e–4 ka .

H2
2
= (1 – 2 ka ) e2 ka + i 1 – (1 – 2 ka ) e4 ka
H1(42)= ei (π – θ ) ,

(31)

where

(32)

e– iθ : = (2ka –1) e2ka – i

2

1 – (2ka –1) e4ka . (43)where we have used the formula
For the sake of simplicity, suppose ka = 1/2 (θ = π / 2).
In this case, v1x(a,y,t) induced by the surface ξ1 itself
exactly cancels the convective derivative term. Therefore, the surface displacement is locked in the laboratory
frame due to the velocity field induced by itself (ξ1), andis amplified due to that induced by the other (ξ2).
The schematic view of the unstable eigenstate Eq.
(42) is illustrated in Fig. 6. The perturbed velocity
induced by the deformation H1 at the surface x = a isrepresented by the straight arrows. Since we assumed
the eiky dependence on ξ1 and ξ2, Eq. (42) implies thatthe phase ξ2 advances in y direction with respect to ξ1by an amount π – θ. Thus, the velocity field induced by
the perturbed vorticity field on x = a indicates the
direction of the amplification of the surface
displacement on x = –a, and vice versa. The perturbation
grows exponentially due to the mutual interaction of twophase-shifted surface displacements.
The contour plot of the stream function φ = K Ψcorresponding to the unstable eigenmode is shown in
Fig. 7 (k = 0.5 and a = 1). In this figure, we take the x-∞x sin(α x )dx = π e – αβ .
–∞β2 + x2(33)

On the other hand, the evolution equation for the
perturbed surface is governed by∂t ξ1,2 = v1x – Vy ∂y ξ1,2,(34)

ξ2 (y, t) = H2 ei(k y – ω t)

(35)

where
denotes the surface displacement on x = –a. The firstterm on the right hand side in Eq. (34) implies thetransfer of the surface due to perturbed velocity, and thesecond one, the stretching effect of the ambient flow.
The y component of the perturbed velocity is omitted
since it gives a higher order contribution with the
combination of ∂yξ1,2. We thus estimate the perturbed
velocities on x = ±a by putting ζ = 0 and ζ = ±2a in Eq.
(32), which yields

v 1x (a , y , t ) = i σ e i (ky – ω t ) H1 – H2 e–2ka ,2(36)v 1x (– a , y , t ) = i σ e i (ky – ω t ) H1 e–2ka – H2 . (37)2
Substituting Eqs. (36) and (37) into Eq. (34), we obtain–2ka
1 – 2ka + 2 ω
H2 = 0 ,σ H1 – e(38)

e–2ka H1 – 1 – 2ka – 2 ω
σ H2 = 0 .

(39)Fig. 6 Schematic view for the unstable eigenfunction(42) or (62) of Kelvin-Helmholtz instability.

The consistency condition finally yields the dispersion
1175

Journal of Plasma and Fusion Research Vol.79, No.11

supports the mode for escaping from the stretchingeffect.
When wave number k is small, the phase velocity
v ph of the surface wave is so large that it has a highcapability for the transfer of the vorticity, while thevelocity of the ambient flow is fixed (U). As k becomeslarge, the phase velocity vph becomes small, and the
maximum growth rate is achieved around | vph | ~– |Vy |.
Namely, the mode becomes the most unstable when thesurface wave tends to stop due to the effect ofbackground flow. If the wave number k is too large, thephase velocity v ph becomes so small that it cannot
compete with the background flow. In other words, the xcomponent of the velocity created from the ambient
flow Vy ∂y ξ1,2 is too large in Eq. (34), and v1x of itself
cannot afford to lock the surface. Thus, the instability
will be suppressed.

Fig. 7 Contour plot of the real part of stream function.
The figure is plotted for k = 0.5 and a = 1.

3.2 Eigenvalue analysis
In this section, we define a generalized RayleighequationFig. 8 Surface wave due to vorticity variation.

axis in the horizontal direction. We observe that Vy′vxvy< 0 is satisfied throughout in the finite flow shear region(–1 < x < 1) [see Eq. (143) in Sec. 5.2]. The damped
mode has an opposite phase difference between surface
waves, which yields the opposite inclination of the flow
pattern. While the ambient shear flow is in the positive
y-direction in x > 0 and toward the negative in x < 0, themode structure is inclined in the opposite direction to it.
This may be attributed to the surface wave which wewill describe in the following section.
In order to see the surface wave explicitly, it is
helpful to consider the case of a single surface ofvorticity separation. The surface displacement willproduce the velocity field v1x with π /2 out-of-phase [seeEq. (36)]. Due to this effect, the surface displacement
can propagate as a wave with respect to the rest frameof the fluid (see Fig. 8). The dispersion relation of thesurface wave appears asω =– σ ,2i ∂t Ψ = LΨ(46)

L : = kVy (x) + kW(x)K(47)

in the infinite domain (x ∈ R), where the Green operator
K is represented by a convolution integral+∞(K f )(x ) =–∞e– k x – ξ f (ξ ) dξ .
2k

(48)

In the following, we denote the Green function byK(x,ξ );K (x , ξ ) = e–k x –ξ

2k

.

(49)

In Eq. (47), we take Vy(x) and W(x) to be independent
arbitrary functions. The case in which W(x) = Vy″(x)recovers the physically relevant equation (26).
Firstly, let us assume W(x) = 0 in Eq. (47) andconsideri∂tΨ = kVy (x)Ψ(50)

with a continuous real function V y (x). The formal
eigenvalue λ and the corresponding eigenfunction ϕ of
the generator of Eq. (50) are given by(44)

whose phase velocity isv ph = – σ .
2kNovember 2003(45)

λ = kVy (µ),ϕ = δ (x – µ),(51)

where µ is an arbitrary real number. The solution of Eq.
(50) with initial condition Ψ (x,0), thus, yields

Thus, the vorticity variation becomes a source of the
wave which corresponds to the Vy″ term in Eq. (26). It
is, therefore, deduced that the surface wave propagatesin the opposite direction to the ambient flow, whichΨ (x, t) = e –itkVy (x)Ψ (x, 0).
1176

(52)

Review Paper

T. Tatsuno

Stabilizing and Destabilizing Effect of Shear Flow

Two eigenvalues of the 2 × 2 matrix L 3 denote theKH mode [1], whose dispersion relation reads

Second, we assume
Vy (x) = 0,(53)

W(x) = –σδ (x–a),(54)

λ 2 = σ (1 – 2ka )2 – e–4ka .
42When |(2ka – 1)e–2ka | < 1, the unstable branch of Eq.
(61) yields

then the operator reads∞

L 1 = – σ δ (x – a )e– k x – ξ dξ .
2
–∞

(55)

2
λ– = i U
e–4ka – (2ka – 1) ,2a
ϕ – = δ (x – a ) + e– iθ δ (x + a ) .

The eigenvalue and the corresponding eigenfunction of
L 1 are, respectively,

λ =– σ ,2

ϕ (x ) = δ (x – a ) .

W(x) = –σ [δ (x – a) – δ (x + a)],

(56)

(57)

the corresponding eigenvalues and eigenfunctions of theresulting operator L 2 are

λ = ± σ 1 – e–4ka ,2ϕ (x ) = δ (x – a ) – 1 ± 1 – e–4ka e2ka δ (x + a ) .

(58)
Then, we may formally write the operator L 2 in thematrix form by means of the basis vectors δ (x – a) andδ (x + a),– 1 – e–2ka
1
e–2ka

.

(59)

These oscillations represent the coupled surface modes
which are excited on both surfaces x = ± a. Note thatthese modes are stable for positive k in the absence ofcontinuous spectra. Since the governing equation is thesame as in the case of diocotron instabilities in nonneutral plasmas [31,32], these oscillations are called‘diocotron oscillations.’
The flow profile Vy (x) discussed in Sec. 3.1 isaccompanied by the W(x) of Eq. (57). Let us consider
the consistent profile given in Fig. 3 [Vy″(x) = W(x)] anddenote the resulting operator by L 3. Then, by focusing
on the subspace spanned by the two surface waves δ (x –a) and δ (x + a), we may obtain the following matrix
form for the operator L 3;0
L 3 = kU
+L2.
0 – kU(62)

Because of the eiky dependence of the vorticity, theeigensolution (62) states that the phase of the surface
wave on x = a advances in y direction to that on x = –aby an angle θ (0 < θ < π) [see Eqs. (41)–(43)].
This expression exactly coincides with the situation
illustrated in Fig. 6. According to the surface wave
perturbation of the form of Eq. (62), we plot thecoefficient of the δ-function by the solid curve on x = ain the right figure. Since the phase of the surface wave
on x = a advances in y direction to that on x = –a by
the angle θ, we may plot the amplitude of the surface
wave on x = –a by the dashed curve. Since the ambientvorticity is positive in –a < x < a and zero elsewhere,
we can directly regard the perturbed amplitude on x = a
as an illustration of the surface displacement which
separates the two different vorticity region. However,the surface displacement on x = –a is inverted from theperturbed amplitude of the vorticity because of theambient vorticity field (positive in upward and zero inlower), which yields the illustration of the surface
displacement by the solid curve.
Let us close this section by considering the effect

Eigenmode (56) corresponds to the surface wave (44).
When we include another surface of different vorticity
regions and assume

L2=σ
2

(61)

(60)

Fig. 9 Square of the growth rate when rigid wall is put
along x = ±b. The values of b are normalized by a.

1177

Journal of Plasma and Fusion Research Vol.79, No.11

of the finite domain size by putting a rigid boundary at x
= ±b (b > a). The boundary condition for finite ky yields
φ = 0 on x = ±b. The dispersion relation then reads

ω 2 = k 2 a2 –σ2λϕ (x ) = kVy (x ) ϕ (x )∞– σ δ (x – a )e– k x – ξ ϕ (ξ ) dξ . (68)2–∞2ka (1 + X ) + (1 – X ) e–4ka – 1 + e–4ka
2

2(1 + X ) – (1 – X ) e–4ka

For this eigenvalue problem, we have the following setsof eigenvalues and the corresponding eigenfunctions:
1. For λ0 = kU; the corresponding eigenfunctions
are arbitrary continuous functions ϕ0(x) which
satisfy

,

(63)
where
–2k ( b – a )1 < X = 1 + e–2k ( b – a ) < ∞ .
1– eNovember 2003

∞

(64)

ek ( a – ξ ) ϕ0 (ξ ) dξ = 0 ,and ϕ0 (x ) = 0 (x ≤ a ) .

a

By assuming ε : = ka << 1 with a and b – a kept finite,
we can put 1/X ~– k (b – a) << 1 and Eq. (63) isapproximated to O(ε 2) by2ω 2 ∼ a k 2 (2a – b ) .
σ2 b

2. For λ1 = kU – σ / 2; the corresponding eigenfunction is

ϕ1(x) = δ (x – a).

(65)

ϕµ (x ) = δ (x – µ ) +e– k ( a – µ ) δ (x – a ) .
2k (a – µ ) – 1(71)

However, these eigenvalues are not complete andwe have another eigenfunction in a wider sense(nilpotent). That is ϕ2(x) = δ(x – µ0) where µ0 := a –1/2k. We can easily see that

3.3 Nilpotent and resonance
We have studied the eigenmodes consisting of theRayleigh equation in the previous section. The
frequency of the surface wave is negative in the fluid
frame; therefore, the eigenvalues of Eqs. (51) and (56)[or Eq. (61)] may overlap when |(2ka–1)e–2ka| > 1 issatisfied. Frequency overlapping between modes in nonHermitian systems may cause algebraic instabilities
even if all eigenvalues are real. In this section, we show
the effect coming from such resonance or frequency
overlapping.
For simplicity, we consider the single surface ofvortical discontinuity (54) in Sec. 3.3.1, and next
consider the double-surface system corresponding to Eq.
(57) in Sec. 3.3.2. Here, the velocity field Vy(x) is taken
as consistent to W(x).
3.3.1 Single surface system
Let us first consider the following case(λ 1 – L 4) ϕ2 (x ) = U e– k ( a – µ 0 ) ϕ1 (x ) .
2a

(72)

Of course (λ1 – L 3)2ϕ2(x) = 0 also holds. These relations
are quite similar to the nilpotent in the finite dimension
operator. Thus we may write the operator L 4 in thefollowing matrix form including the case µ = µ0:

kU – σ2
L3=0– σ e– k ( a – µ )
2
,
kσµ

(73)

where a Jordan block is obtained for kU – σ /2 = kσµ (µ= µ0).
Let us evaluate the time evolution of the perturbation when we have taken this ϕ 2 (x) as an initialcondition. As we can see from Eq. (72), we will haveϕ1(x) component by applying the generator L 4 to theinitial condition ϕ 2(x). Here we may consider theevolution to be closed in the subspace spanned by ϕ1(x)and ϕ2(x). Then, it is natural to expand Ψ as∞

L 4 = kVy (x ) – σ δ (x – a )e– k x – ξ dξ , (66)2–∞(x < a ),(a ≤ x )(70)

3. For λµ = kUµ/a (µ < a and µ ≠ a – 1/2k); thecorresponding eigenfunctions are

If b/a becomes smaller than 2, the KH mode may bestabilized (ω2 becomes positive in k << 1). Dispersion
relation (63) is numerically plotted in Fig. 9, whichclearly shows that the KH mode is stabilized for all
wave numbers when b/a < 2.

Vy (x ) = σ x
U

(69)

(67)

2

Ψ (x , t )=where σ = U/a. By assuming ϕ ∝ e–iλ t, we may obtainthe following eigenvalue problem.

Σ α (t )ϕ (x ) .
i=1

i

i

(74)

Substituting Eq. (74) into the original equation (46), we
obtain1178

Review Paper

We find in the eigenvalues the coupling between
two diocotron oscillationsi ∂t ( α 1 ϕ1 + α 2 ϕ2 ) = λ 1α 1 – U α 2 ϕ1 + λ 1 ϕ2 .
2α e

λ †1 = σ e–2ka sinh ψ ,2
λ †3 = – σ e–2ka sinh ψ ,2(75)
When we decompose Eq. (75) into ϕ1 and ϕ2, we obtainthe time evolution of α1

α 1 (t ) = i σ α 2 (0)t + α 1 (0) e– i λ 1 t ,
2 e

(84)
(85)

and the corresponding eigenfunctions are(76)

where α1,2(0) denote initial values. We find that forα 2 (0) ≠ 0, we have a secular growth due to theresonance of the diocotron mode with one of thesingular eigenfunctions in the continuous spectrum.
3.3.2 Rayleigh’s system
Let us consider the case shown in Fig. 3. The
generator is written asϕ†1(x) = δ (x – a) + e–ψ δ (x + a),(86)

ϕ†3(x) = e–ψ δ (x – a) + (x + a),(87)

where cosh ψ = (2ka – 1)e 2ka. By introducing new
coefficientsβ1eψ – 12=β3sinh ψ – 1 eψ

L 5 = kVy ( x )

α1
α3 ,

(88)

we may diagonalize the diocotron modes, which yields∞– σ δ (x – a ) – δ (x + a )2

e– k x – ξ dξ ,

–∞

i

(77)
where the velocity field is defined by

Vy (x ) =(x ≤ – a )(– a < x < a ) .
(a ≤ x )

–U
σx
U

i

(78)

i

2ka

Note that we consider the case | (2ka – 1)e | > 1, where
KH mode is exponentially stable.
If we set the basis vectors as ϕ1(x) = δ (x – a), ϕµ(x)= δ (x – µ ), and ϕ3(x) = δ (x + a), we can obtain thefollowing matrix representation for the operator L 5;

L5=σ2

T. Tatsuno

Stabilizing and Destabilizing Effect of Shear Flowd β1 λ † β
= 1 1
dt
– σ
e– k ( a – µ ) eψ + e– k ( a + µ ) α 2 ,sinh ψ

(89)

dα 2= kσµ α 2 ,
dt

(90)

d β3 λ † β
= 3 3
dt
+ σ
e– k ( a – µ ) + e– k ( a + µ ) eψ α 2 .
sinh ψ

(91)

The solution of Eq. (90) is given by

α2(t) = α2(0) e–i kσµt.

2ka – 1 – e– k ( a – µ ) – e–2ka
2kµ0. (79)
0e– k ( a + µ ) – (2ka – 1)
e–2ka

(92)

When we assume
kσµ = λ†1 or λ †3 ,(93)

When we expand the perturbed vortex field3

Ψ (x , t )=

Σ α (t )ϕ (x ) ,
i=1

i

i

(80)

regarding ϕ2 as one of the ϕµ corresponding to a certainµ, and substitute Eq. (80) into the Rayleigh equation, we
obtainidα 1 σ=(2ka – 1) α 1 – e– k ( a – µ ) α 2 – e–2 ka α 3 , (81)
dt
2

i

dα 2
= kσµα 2 ,dtidα 3 σ –2 ka α
– k (a + µ )α 2 – (2ka – 1) α 3 . (83)
=
e
1+ e
dt
2

Fig. 10 Schematic view of the spectra for the KH (a)unstable and (b) stable systems. The crosses andwavy lines denote point and continuous spectra,
respectively. Point spectra are embedded in thecontinua for KH stable case. Open circles denotethe vacancy of continua (corresponding tonilpotent).

(82)

1179

Journal of Plasma and Fusion Research Vol.79, No.11

November 2003

then we might have secularity due to the resonance withthe simple oscillator ϕ 2†. On the other hand, we do nothave it when the eigenvalues λ †3 are complex or pure
imaginary, namely when the system is unstable in a KH
sense (see Fig. 10).

U) and vorticity, respectively.
Because Vy″(x) yields delta functions, the perturbed
vorticity Ψ must include delta functions, representing
the surface wave perturbations; We write

3.4 Coupling with gravity wave andsecularity
In this section, we include the density gradient into
the Rayleigh equation, and show that the coupling to thegravity wave will introduce a localized secular behavior.
The continuity equation (4) couples with the equation of
motion (3).
We consider the following stationary state:

where α1(t) and α2(t) represent the amplitudes of thesurface waves and f (x, t) is the continuous part of thevorticity fluctuation. Then, we obtain

Vy (x ) =Ψ(x,t) = α1(t)δ (x + 1) + α2(t)δ (x – 1) + f (x, t), (97)i∂t f = kxf (x,t) + kVy″(x)φ(x,t) + kvg2χ(x,t),dα 1= kα 1 (t ) + kφ (–1, t ) ,
dt
dα 2
i
= kα 2 (t ) – kφ (1, t ) ,dt–α a /gρ e
ρ 0 (x ) = ρ 0 e–α x /g
ρ 0 eα a /g
0

(99)

i

–U
Ux / a ,
U

(98)

(100)

i∂t χ = kxχ(x, t) + kvg2φ(x, t),

(101)

where

(x < – a )(– a < x < a ) ,(a < x )(94)
φ (x , t )= K Ψ= 1 α 1 (t ) e– k x + 1 + α 2 (t ) e– k x – 1
2kwhere U, a, and ρ0 are constants. Let us normalize thevariable by a for length and a/U for time. By
introducing a new variable χ := ρ1/ρ0′ and normalizing it
by /a, we obtaini∂tΨ = kVyΨ + kVy″φ + kvg2χ

(95)i∂t χ = kVy χ + kφ,

(96)

+∞

+

e– k x – ξ f (ξ , t ) dξ ,

(102)

–∞

in the infinite domain. Or in the matrix form, we obtiani ∂tf (x )
χ(x )
α1
α2

=

0
kx
k v 2g1 –k x + 1ek K x kx2
1
k K –1 0 – 2 (2k – 1)
–k K 1 0– 12 e– 2kwhere we have introduced a wave number k in ydirection, and vg2 = ρ0′ ga2/ρ0U2 (> 0) and Ψ = –∆φ are
the phase velocity of the gravity wave (normalized by

0
1
2

1
2

e

–k x – 1

1
2

e– 2k

(2k – 1)f (x )
χ(x )
.
α1
α2

(103)
The equations derived here are exactly equivalent
to those in Ref. [11]. Here, we have rederived themfrom a different set of equations in order to demonstratethat this structure has wide generality.
We rewrite λ ± , µ ± more generally as

Fig. 11 The vorticity f(x,t) at t = 65. The initial conditionis ƒ(x,0) = √⎯2cos(πx/2). The amplitude increases
near the resonant surfaces x = µ±(~– ±0.495).

1180

λ±= ± 12

(2k – 1) – e– 4k ,(104)µ± = ± 1
2k(2k – 1) – e– 4k .

2

(105)

2

Review Paper

Observing f (x) at x = µ+ (or x = µ–), Eq. (103) reads

i ddt

renormalized damping rate ν agrees very well with thenumerical results.
We stress here that such an algebraic instability is
not only a mathematical artifact. As we have seen,
algebraic growth is observed when the spectra of theRayleigh equation shows the KH-stable structure [seeFig. 10(b)]. In the infinite domain, KH-unstable(exponential) mode certainly exists for sufficiently smallk (see b = ∞ case in Fig. 9). However, when a rigid wallis placed close enough (b ≤ 2 in Figure 9), the spectra of
the Rayleigh equation becomes that of Fig. 10(b) for all
wave numbers. Thus, we may deduce that there is asystem in which algebraic growth dominates theinstability.

f ( µ+)
χ ( µ+)α1
α2k µ + k v2g=

T. Tatsuno

Stabilizing and Destabilizing Effect of Shear Flow00–k µ + + 1

0

kµ

0

0

– 12 (2k – 1)

00

– 12 e – 2k

1
+ 2

e

1
2

e

1
2
1
2

f ( µ+)
χ ( µ+).
α1
α2–k µ + – 1

e – 2k(2k – 1)

(106)
The generator of Eq. (106) can be written in a Jordan
canonical formTλ+ 1 0 000 λ+ 10 λ+ 0
0
0 λ–
0
0

T

–1

,

4. Couette Flow in the Infinite Domain(Kelvin’s method)
4.1 Time-dependent eigenmodes
Kelvin’s method can resolve, for some classes ofmean flows, the evolution of the system(107)∂t u + V0 · ∇u = A u

by means of a non-orthogonal transform T . Therefore,
we conjecture that f(µ+, t) ∝ t 2 due to the third orderdegeneracy of the frequency.
However, we observe in Fig. 12 that the asymptotic
behavior of f (µ +, t) is not proportional to t2, but f (µ +, t)∝ t. This is because of the integral terms in Eq. (103),
which contribute to the phase mixing damping of thesurface waves. The phase mixing damping of the surface
wave can be evaluated by renormalizing the perturbative
analysis with the assumption v g2 << 1 [33]. The(108)

into new types of modes by means of which bothtransient and secular asymptotic behaviors are
effectively described. Here A denotes a Hermitiandifferential operator (time-independent) defined in a
Hilbert space V, and u (∈V ) denotes a state vector of theperturbed field. Kelvin’s method consists of the combined application of two methods which have been used
extensively in the analysis of wave equations. Precisely
the “Lagrangian” part of Eq. (108), ∂t + V0 · ∇, is solved
by means of the characteristics method, and the“Hermitian” part A by means of the standard spectral
resolution.
The characteristics method is applied to solve thecharacteristic ODE associated with the Lagrangianderivative moving along the characteristic curve of theambient motion, which is given bydx = V ,0dtx (0) = ξ .

(109)

By inverting the modes, which are denoted in
Lagrangian coordinates by ϕ (k,ξ ), we represent themodes in Eulerian coordinates asϕ̃ (t; k, x) = ϕ (k, ξ (t; x)),Fig. 12 Time evolution of the vorticity fluctuation ƒ(x,t )at x = µ+. The amplitude starts to grow inproportion to t 2 and converges to asymptotic
growth ∝ t. Dashed lines denote the analyticalfitting curve for a given damping rate ν.

(110)

where ξ (t; x) denotes the inverse of x(t; ξ ). The
existence of the inverse mapping x(t) |→ ξ is guaranteedin the case of incompressible mean flows. Due to Eq.
(110), ϕ̃ (t; k, x) satisfies the characteristic equation
1181

Journal of Plasma and Fusion Research Vol.79, No.11∂t ϕ̃ (t; k, x) + V0 · ∇ϕ̃ (t; k, x) = 0.

behavior of electromagnetic mode (k · B0 ≠ 0). In thecase of B0 · ∇ ≠ 0, we obtain from Eq. (3)–(5) (see Ref.
[12])

(111)

The essential condition for the applicability of
Kelvin’s method consists in the constraint for thefunctions ϕ̃ (t; k, x) to form the complete set of
eigenfunctions of the operator A [18]. If such a set of
eigenfunctions exists, we can decompose the perturbed
field u by means ofu=u k (t ) ϕ (t ; k, x ) dk .

∂t + Vy ∂y ∆ ∂t + Vy ∂y ψ2=

ψ ( x, t ) =(113)

where the eigenvalue λk also depends on time.
Plugging Eq. (112) into Eq. (108) and exploiting(111) and (113), we obtainu k ( t ) λ k ( t ) ϕ ( t ; k, x ) d k .

ψ k ( t ) ϕ ( t ; k, x ) d k ,ϕ̃ (t; k, x) = exp [ikx x + iky(y – Vy t) + ikz z]= exp [ik̃ x x(t)x + ikyy + ikz z],

(117)

(118)

(114)
and k̃ x (t) = kx – ky σ t. Since Eq. (118) spans the function
space V, Eq. (117) gives a complete expansion by meansof ϕ̃ [18]. Eigenmodes (118) obeys the governing ODE[see Eq. (115)]

Due to the orthogonality of the modes ϕ̃ (t; k, x), theevolution of ûk is governed by the equationd u (t ) = λ (t ) u (t ) .
k
k
dt k

(116)

where each eigenmode can be expressed by thesinusoidal function∂t u k ( t ) ϕ ( t ; k, x ) d k=

B0 ⋅ ∇
ρ′ g∆ψ – 0 ∂2y ψ ,µ 0 ρ0ρ0where ∆ = ∂ x2 + ∂ y2 denotes the two-dimensional
Laplacian. Since the operator on the right hand side is
Hermitian, we can decompose the flux function ψ bymeans of the shearing eigenmodes [see Eq. (112)](112)

We notice that due to Eq. (110) the eigenvalues of A
become time-dependent. The new eigenvalue problem
for A reads
A ϕ̃ (t; k, x) = λk(t)ϕ̃ (t; k, x),November 2003

(115)

2

dψ
d ψ µ
+ (t )+ 1 – S (t ) ψ = 0 ,2
dt
dt

If ϕ̃ (t; k, x) do not satisfy both conditions given bycharacteristic equation (111) and eigenequation (113),
then Eq. (114) will have additional terms representing
the complicated mode coupling, and thus, theapplicability of Kelvin’s method is compromised.
The evolution of û k(t) will not exhibit a simple
exponential dependence due to the time dependence
present in the eigenvalues λk(t). By analyzing this ODE,
we can classify the transient motion and the time
asymptotic behavior of each mode.

(119)

where

µ (t ) = –2σ k y k x (t )22,k x (t ) + k y
2

S (t ) = –

ky G
2

2

,

k x (t ) + k y

and we have omit the subscript k for ψ̂ . We havenormalized the time t by the poloidal Alfvén time τA =a µ 0 ρ 0 /(k · B0), the wave vector k by the characteristic
length scale a, and τG2 = –ρ0 / ρ0′ g and G = τ A2 / τG2 . We
notice that in the absence of shear flow (σ = 0) the usualinterchange instability equation for static equilibrium
can be obtained from Eq. (119) [see Eq. (20)].
When σ ≠ 0, we have µ(t) ≠ 0 and we can draw ananalogy with the dynamics of a damped oscillator with
time-dependent frictional coefficient µ(t). With time,µ(t) becomes always positive, which means a formal
dissipation, and therefore the oscillation energy of the
Alfvén wave [(dψ̂ /dt)2 + ψ̂ 2]/2 decreases monotonically.
In order to study the time asymptotic behavior, we
assume t >> k x /σ ky , 1/σ . In this limit we obtain the

4.2 Asymptotic and transient behavior
In this section, we show the behavior of Kelvin’s
mode [time-dependent eigenmode Eq. (110)] for
reduced MHD equations (3)–(5). We consider a plasma
slab of infinite domain with Vy(x) = σx and B0 = const.
The behavior is quite different depending on whether
the perturbation is coupled with Alfvén wave or not.
Thus, we separately discuss the behavior of theelectromagnetic (k · B0 ≠ 0) and electrostatic (k · B0 = 0)
modes.
4.2.1 Electromagnetic mode
Let us first see the time-asymptotic and transient1182

Review Paper

4.2.2 Electrostatic mode
When the wave vector is purely perpendicular to
the ambient magnetic field (k · B0 = 0), the flux function
ψ decouples from ρ and φ . From Eqs. (3)–(4), thedynamics of the stream function is governed byfollowing ODE;2d ψ + 2 d ψ + 1 – G / σ2 ψ = 0 ,
(120)
2
2
t dt
dt
t
where G = τA2/τG2 denotes the magnitude of the instability
drive term. Since Eq. (120) is the spherical Bessel
equation, the time asymptotic behavior of the modeyieldsψ ∼ 1 sin t – πν + δ ,t2

2∂t + Vy ∂y ∆φ = –(121)

φ ∼ 1 cos t – πν + δ .
t
2

2

d
2
dt

2

+ψ

2

dψ
= – µ (t )
dt

(122)

2

2

2

k x (t ) + k y φ = k y γ G2 φ ,

(125)

2
d φ + 4 d φ + 2–α φ =0,
(126)
2
2
t dt
dt
t
where α = γG2 /σ 2 denotes the ratio between the inter-change destabilizing effect and the flow shear stabilizing
one (Richardson number). The general solution of Eq.
(126) isφ̂ = C1t m+ + C2t m–

(127)

where

m± =

– 3 ± 1 + 4α
.
2

(128)

The asymptotic time behavior is therefore dominated by
the larger index m+. Thus we can state the condition forthe boundedness of φ̃ as

2

.

(124)

where γG2 = – ρ0′ g/ρ0(= τG–2) denotes the characteristic
growth rate of the interchange instability. Here we have
dropped the subscript k for the sake of simplicity. In
order to investigate the asymptotic time behavior ofeach mode, we again assume t >> kx /ky σ and t >> 1/ σ.
Then Eq. (125) becomes

We note here that there is no threshold value for thestabilization of the interchange instability, since thespherical Bessel equation (120) is common to all modes.
Namely, the combined effect of the Alfvén wave
propagation and shear flow mixing always overcomesthe interchange drive and the oscillations of themagnetic flux asymptotically decay proportional to theinverse power of time.
Since an analytic expression is not available for thetransient behavior of each eigenmode, we discuss it byqualitatively analyzing the ODE (119). In the absence of
the instability drive, we will havedψdt

ρ 0′ g 2∂ φ.
ρ0 y

When we represent φ in terms of the shearing mode
given in Eq. (118), we obtain the ODE

where δ denotes a constant phase depending on theinitial conditions. Therefore the mode oscillates withamplitude ψ̂ decaying with the inverse power of time.
The dependence Eq. (121) corresponds to the similar
asymptotic decay in the stream function φ, which reads
in the leading order asddt

T. Tatsuno

Stabilizing and Destabilizing Effect of Shear Flow

(123)

Since the sign of the denominator in µ (t) is always
positive, its behavior will be determined by that of thenumerator. The numerator is written as 2σ 2ky2 t – 2σ ky kx
and according to its initial value we can single out twoclasses of the transients.
When the product σ ky kx is negative, the frictionalcoefficient µ(t) is always positive from the beginning;
therefore, the shear flow acts as a damping force at any
time and the mode shows simple damped behavior. On
the other hand, if the product σ ky kx is positive, thefrictional coefficient µ (t) is initially negative andchanges its sign at the instant t* = kx / σ ky. Therefore themode experiences an initial amplification lasting until
the time t*, which is even faster than it would be in thepresence of only the interchange drive.

ρ′ g
α ≤ 2 ⇒ – 1 0 ≤ σ2 .
2 ρ0

(129)

The condition for the boundedness of φ̃ is improved
compared with the static case ( ρ 0′ ≥ 0) due to thestretching effect of the shear flow. The electrostatic (k ·B0 = 0) mode can be linearly unstable while the electromagnetic (k · B0 ≠ 0) mode is completely stabilized [see(121) and (122)]. However, the asymptotic behavior is
algebraic in time.

5. Couette Flow in the Finite Domain
5.1 Perturbative analysis of eigenvalues
We consider a finite domain [–a, a] in this section
with a linear flow such as (see also Fig. 13)
1183

Journal of Plasma and Fusion Research Vol.79, No.11

Vy(x) = σ x,

and eigenvalues as(130)

φ = φ0 + φ1 + φ2 + …,

where the shear parameter σ is a real constant. The
conventional KH instability is absent because Vy″ = 0[3,14]. Transforming ∂t → –iω and ∇ → ik, we obtain2

2 k y σ k ||2
d φ
– ky φ –22dxΩ Ω 2 – k2
2

2

d φ0
2
– k y 1 + G 2 φ0 = 0 ,dx 2
ω 20 – k ||dφ
ky σ φ + Ω
dx

2

ky G

φ=0.
2

ω = ω0 + ω1 + ω2 + …, (132)

where |ωn|/|ω0| ~ |φn|/|φ0| ~ O(σ n). From the O(1) terms,
we obtain||–

November 2003

(133)

which has been already solved in Sec. 2.3. Eigenvalues
and eigenfunctions are(131)Ω 2 – k ||

2

ky G

2

ω (0n ) = k 2|| –under the normalization in terms of the Alfvén velocity
vA = B0 / µ 0 ρ 0 and the system size a, where Ω = ω –kyVy(x) is the Doppler shifted local frequency.
We now assume σ << 1 and represent a perturbative
analysis for the ODE (131). We expand eigenfunctions2

k y + n2 π 2 / 4

,

(134)

and

cos (n π x / 2)
sin (n π x / 2)

φ (0n ) =for n : odd,
for n : even(135)

respectively. We study the effect of the shear flow for
the mode n = 1, which is the most unstable when V0 = 0.
The terms of O(σ) yields
2

d φ1
2
– k y 1 + G 2 φ1
2dx
ω 20 – k ||2–

2
d φ 0 2ω 0 k y G
ω 1 – k y σx φ0 = 0 .
+
2
dx
2
2
ω 0 – k ||2 k y k || σ

ω 0 ω 20 – k 2||(136)

Multiplying φ 0(1) on both sides of Eq. (136) andintegrating it over the domain, we find(137)

ω1 = 0.
Fig. 13 Slab geometry with gravity.

The next order equation, thus, describes the effectof the flow shear on the eigenvalue.
After the tedious but straightforward calculation
due to the expansion of φ1 in terms of φ0(n), we finallyobtain the second order dispersion relation (see Ref.
[13])
2

2ω 0 ω 2 = –k || σ 22
1 + 4 P 4 + B + 162 P 4 + A
πG2+ k y σ 2 166 P 4 + A + 3B ,
π

(138)

where
∞

A=

Σ

m= 1

Fig. 14 The growth as a function of the shear parameterσ (G = 2.72 and ky = k|| = 0.5). Numerical results
are obtained from the direct solution to theeigenvalue problem Eq. (131) by the shooting
method.

m2
m 2 – 1/ 4

5

∼ 4.219581 ,B = 1 – 22 ∼ 0.130691 ,
3 π

1184

(139)

(140)

Review Paper

Stabilizing and Destabilizing Effect of Shear Flow

P4 + = k y + π .
4
2

2

(141)

When ω0 describes instability (pure imaginary withIm ω0 > 0), the first term on the right hand side of Eq.
(138) brings about a destabilizing effect, while thesecond causes a stabilizing effect. The destabilizing term
contains both k || and G, while the stabilizing one
contains only k y. It is remarkable that k2|| works toincrease the growth rate, and G, decreases it, which is
opposite to the conventional understandings for static
equilibrium. For k || = 0 (or for a neutral fluid), thisdestabilizing effect does not work. There is a thresholdin the ratio of k2|| and G, where the coefficient of σ 2changes its sign.
Variational calculation yields the maximum growthrate with respect to the local maximum flow velocity
Vmax = σFig. 15 Streamlines of the typical eigenfunction. The
parameters are Bz = 0, G = 2.72, k|| = ky = 0.5, andσ = 0.7, and 1.5.

static equilibrium (V y = 0), the eigenfunctions
determined by Eq. (20) are real functions. Then, we
obtain mirror symmetric streamlines, which leads to theexact cancellation of the vx v y term when integrated.
However, an ambient shear flow (Vy′ ≠ 0) yields thecomplex-valued eigenfunction φ, where the integral ofvx vy remains finite because of the breakdown of thesymmetry. Figure 15 shows the streamlines of theeigenfunction for the cases σ = 0.7 and 1.5. We observethat the product vx vy has negative value in most of thedomain, if σ = 0.7. Since Vy′ = σ is a positive constant,this term (representing the work done by the shear flow)
gives a positive contribution in the right hand side ofEq. (143). It is remarkable that the stream line contour is
inclined in the opposite direction to the ambient shear
flow. The mode corresponding to the complex conjugate
eigenvalue (damping mode) has the opposite structure— the streamlines are distorted in the direction of theflow. As σ increases, the ambient flow begins to distort
the mode to the direction of the flow, and finallystabilizes it.

2

Vmax ∼

k || vA ωA
=
,
ky
ky

(142)

where ωA is the Alfvén frequency. Over this criticalvelocity (or σ ), the growth rate starts to diminish,
implying the onset of the shear flow stretching effect.

5.2 Discussion
Shear flows may stabilize instabilities when they
stretch the fluctuation. However, the condition (142)shows that the stabilization occurs only if the local
velocity exceeds the phase velocity of the Alfvén wave.
When the flow shear is weak, modes may be evendestabilized; the growth rate achieves a maximum value,
as a function of the maximum flow velocity, near the
Alfvén velocity [see Eq. (142)], and it is then suppressed due to the stretching effect overcoming the
Alfvénic phase propagation against the flow. We notethat this situation is very different when we consider amathematical model of infinite-domain linear shear flow
(Sec. 4). For such a model, the maximum velocity is
unbounded, and hence, we see only the stabilization of
instabilities [12].
Let us see how some interchange modes aredestabilized in a weak shear flow. Under theincompressibility condition, the energy balance in theperturbed field reads [see also Eq. (6)]ddt=

1 ρ v 2 + 1 b 2 + g ρ2 dv
1
2µ0
2 0
2 ρ′0
V ′y
– ρ0 V ′y v x v y + µ b x b y dv .
0

T. Tatsuno

6. Summary
We have reviewed recent studies for constructing alinear stability theory of shear flow plasmas from theviewpoint of spectral analysis. Mathematically, the nonHermiticity of the operator due to shear flow does not
only imply complex eigenvalues, but also theincompleteness of spectral resolution. Physically, shear
flow not only has a stabilizing effect due to its
stretching mechanism, but also has a destabilizing effectby causing the mode for stationary standing structure in
the face of an ambient wave.
We first revisited Rayleigh’s analysis of KelvinHelmholtz instability for a piece-wise linear shear flow
profile and clarified the necessity of nonzero Vy″ forinstability. The term containing Vy″ carries the surface
wave which supports the mode’s stationary structure.
Eigenvalue analysis of the Rayleigh equation indicatesthe overlap of frequency spectra (nilpotent). This yields,(143)

where v2 = vx2 + vy2 and b2 = b2x + b2y . In the case of a1185

Journal of Plasma and Fusion Research Vol.79, No.11

November 2003

Schmid, Phys. Plasmas 7, 1788 (2000).
[5] I.B. Bernstein et al., Proc. R. Soc. London, Ser. A
244, 17 (1958).
[6] J.P. Freidberg, Ideal Magnetohydrodynamics(Plenum, New York, 1987).
[7] A.E. Lifschitz, Magnetohydrodynamics andSpectral Theory (Kluwer, Dordrecht, 1989).
[8] K. Yosida, Functional Analysis, (Springer Verlag,
Berlin, 1995).
[9] P.D. Lax, Functional Analysis (Wiley, New York,
2002).
[10] Z. Yoshida, Hisenkei Kagaku Nyumon (Iwanami,
Tokyo, 1998) [in Japanese].
[11] M. Hirota, T. Tatsuno, S. Kondoh and Z. Yoshida,Phys. Plasmas 9, 1177 (2002).
[12] T. Tatsuno, F. Volponi and Z. Yoshida, Phys.
Plasmas 8, 399 (2001).
[13] T. Tatsuno, Z. Yoshida and S.M. Mahajan, Phys.
Plasmas, 10, 2278 (2003).
[14] J.W.S. Rayleigh, Proc. London Math. Soc. 9, 57(1880).
[15] W. Thomson (Lord Kelvin), Philos. Mag., Ser. 5
24, 188 (1887).
[16] W.O. Criminale and P.G. Drazin, Stud. Appl. Math.
83, 123 (1990); G.D. Chagelishvili, A.D. Rogava
and I.N. Segal, Phys. Rev. E 50, 4283 (1994); A.D.
Rogava and S.M. Mahajan, Phys. Rev. E 55, 1185
(1997).
[17] F. Volponi, Z. Yoshida and T. Tatsuno, Phys.
Plasmas 7, 2314 (2000).
[18] F. Volponi and Z. Yoshida, J. Phys. Soc. Japan 71,
1870 (2002).
[19] H.R. Strauss, Phys. Fluids 19, 134 (1976).
[20] H.R. Strauss, Phys. Fluids 20, 1354 (1977).
[21] T. Tatsuno, PhD thesis, Kyoto University (2002).
[22] J. von Neumann, Mathematical Foundation of
Quantum Mechanics (Princeton Univ., Princeton,
1996).
[23] G. Laval, C. Mercier and R.M. Pellat, Nucl. Fusion
5, 156 (1965).
[24] E.A. Coddington and N. Levinson, Theory of
Ordinary Differential Equations (McGraw-Hill,
New York, 1955).
[25] E.L. Ince, Ordinary Differential Equations (Dover,
New York, 1956).
[26] Z. Yoshida, Shudan Gensho no Suri (Iwanami,
Tokyo, 1995) [in Japanese].
[27] J.A. Tataronis, J. Plasma Phys. 13, 87 (1975).
[28] W.A. Newcomb, Ann. Phys. 10, 232 (1960).
[29] J.P. Goedbloed and P.H. Sakanaka, Phys. Fluids

when we include the density gradient (gravity wave), thesecular growth of vorticity in time (Sec. 3).
Kelvin’s modal approach is a particular form ofgeneralized spectral method to obtain a general solutionof the system; however, it is restricted to a particular
shear flow profile. According to this analysis, all modes
show asymptotically algebraic behavior in time (Sec. 4).
When we consider a finite domain by imposing a
conducting wall, we observe a destabilizing effect ofshear flow in the presence of the Alfvén wave (Sec. 5).
Many unresolved questions remain:
1. When we fix the flow profile and focus on thefield line bending effect, it affects thedestabilizing mechanism for finite domain
Couette flow [see Eq. (138)], while it also
affects stabilizing in the infinite domain [seeEqs. (122) and (127)].
2. Rayleigh’s criterion insists Vy″ must change its
sign in the domain. We found here that nonzero
Vy″ is responsible for the stationary structure of
the mode; however, we still do not know why itmust change its sign. A complete physicalunderstanding of this phenomenon has not beenachieved.
3. Is resonance between continua universal overthe different sources? For example, does similar
secularity happen if Alfvén and flow continua
are in resonance?
4. It is obvious that exponential is not the only
dependence of instability for non-Hermitian
operators. How rich in variety do the linear
asymptotic behavior show for them?
When we consider multidimensional equilibrium, more
and more questions come to mind. Even area-occupied
spectra are observed for non-Hermitian operators[34,35]. We still do not know if we may dream of thepossibility of spectral resolution for a non-Hermitian
operator, or if we may prove any necessary andsufficient condition for stability such as an energy
principle for static equilibria. These questions are left tobe answered in the future.

References[1] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Clarendon, Oxford, 1961).
[2] G.K. Batchelor, An Introduction to Fluid Dynamics(Cambridge Univ., Cambridge, 1967).
[3] P.G. Drazin and W.H. Reid, Hydrodynamic
Stability (Cambridge Univ., Cambridge, 1981).
[4] L.N. Trefethen et al., Science 261, 578 (1993); P.J.
1186

Review Paper

Stabilizing and Destabilizing Effect of Shear Flow

17, 908 (1974).
[30] K. Appert, R. Gruber, F. Troyon and J. Vaclavik,
Plasma Phys. 24, 1147 (1982).
[31] R.C. Davidson, Physics of Nonneutral Plasmas(Addison-Wesley, Redwood, 1990).
[32] S. Kondoh, T. Tatsuno and Z. Yoshida, Phys.

T. Tatsuno

Plasmas 8, 2635 (2001).
[33] M. Hirota, T. Tatsuno and Z. Yoshida, J. Plasma
Phys. 69, 397 (2003).
[34] Z. Yoshida and Y. Giga, Math. Z. 204, 235 (1990).
[35] A. Lifschitz, Phys. Fluids 9, 2864 (1997).

1187