PHYSICS OF PLASMAS 13, 092107 共2006兲

Magneto-flow instability in symmetric field profiles
T. Tatsunoa兲 and W. Dorland
Center for Scientific Computation and Mathematical Modeling, The University of Maryland,College Park, Maryland 20742-3289

共Received 26 June 2006; accepted 28 July 2006; published online 12 September 2006兲
Since Rayleigh’s early work on shear-flow driven instabilities in fluids, it has been known thatsheared flows are usually unstable only in the presence of an inflection point in the velocity profile.
However, in magnetohydrodynamics, there are important instabilities for which no inflection pointis required. In tokamak experiments, strongly sheared flows are associated with transport barriers.
Instabilities that may limit the height and extent of transport barriers are of central importance. Here,
we present linear and nonlinear simulations of an ideal magnetohydrodynamic instability that isdriven by sheared flows without inflection points—instead, the instability mechanism requires
reversed magnetic shear. Several symmetric field profiles are studied. In general, the instability leads
to current profile modifications that push the local minimum value of the safety factor 共qmin兲
upward. The possibility of causing disruption in a relatively slow time scale is pointed out when qmincrosses a rational 共especially integral兲 value. The time scale of the instability is governed by thetransit time of the shear flow, which is typically smaller than that of the Alfvén velocity.
Characteristics of this instability are compared with recent experimental observations.
© 2006 American Institute of Physics. 关DOI: 10.1063/1.2338819兴
I. INTRODUCTION

Shear flow instability in plasmas is a profound problem
that has rich consequences and is not understood well. Shearflow appears in almost all laboratory or astrophysical situations, from tokamaks to active galactic nuclei. The mostwell-known instability incurred by shear flow is the KelvinHelmholtz 共KH兲 instability. The KH instability is associatedan inflection point in the velocity profile, and occurs even
without a magnetic field.1 An investigation of this instability
in plasma, including the effect of the ambient magnetic field,
drift and/or gyromotion, is an area of active research.2
The KH instability is not the only instability driven by
sheared flow in magnetized plasma. It is known that the combination of stable shear flow and an inhomogeneous magnetic field can lead to another ideal magnetohydrodynamic
共MHD兲 instability.3–5 In the case of unmagnetized plane
Couette flow, the isosurface of the perturbed vorticity is
stretched away by the background shear flow since perturbed
vorticity is conserved. However, Stern3 first realized that justas the inflection point in the KH instability breaks thisconstraint,1 an inhomogeneous magnetic field can have asimilar effect, and lead to instability. Kent4 has extended thisidea to the case of a symmetric flow profile, and obtained
several instability conditions in terms of the equilibrium values at the boundary of the domain. Chen and Morrison5 revisited Kent’s work and found a profile as simple as linear
flow and a parabolic field. We will refer frequently to thissimplest case in the present paper.
Previous authors have treated the problem with analytic
methods such as perturbation theory and Nyquist analyses.
They only have obtained the marginal stability condition,and the growth rate has never been given as a function of thea兲Electronic mail: tatsuno@umd.edu

1070-664X/2006/13共9兲/092107/14/$23.00

wavenumber. In this paper, we focus on the single-fluid
MHD equations, take the linear shear flow profile which is
KH stable, and investigate the effect of various profiles ofthe magnetic field in slab geometry. Our approach is primarily numerical. Note that the linear flow profile is itself stablewithout a magnetic field, and any magnetic field profile isideally stable in slab geometry without background flow, so
this instability is truly a combined effect of shear flow andmagnetic field.
The reader might think that magnetorotational instability
共MRI兲 is one of the other cases where the stable shear flow
could be destabilized by the inclusion of a weak magnetic
field.6 However, unlike this instability, MRI has an intrinsic
drive 共centrifugal force兲 other than shear flow itself and thephysics of its nonaxisymmetric perturbation would be closer
to that reported in Ref. 7, which investigates the effect ofshear flow on the plasma with an interchange drive: Shearflow stabilized interchange could reappear when one adds a
magnetic field. Similarly, the so-called “joint instability” isinvestigated in the context of solar physics8 as an extensionof the aforementioned works, where it was found that the
combination of stable shear flow and an inhomogeneous
magnetic field could cause MHD instability. However, theeffect of the curvature makes the problem complicated since
their investigation focuses on the spherical geometry.
Experimentally, an intermediate time-scale instability
was recently observed in the reversed shear discharges of
JT-60U,9 where the growth rate was found to be small compared to an ideal instability but too large for a resistive one.
In more recent experiments,10,11 a considerable amount offlow is often observed, which has a strong shear around theinternal transport barrier. Here, we consider a field profilethat may be relevant to these observations. The inferred flow
profiles might contain inflection points. Our aim is to explore

13, 092107-1© 2006 American Institute of Physics

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-2

Phys. Plasmas 13, 092107 共2006兲

T. Tatsuno and W. Dorland

the magneto-flow instability, which could limit performanceeven in the absence of an inflection point. More detailed
comparisons will be required to determine which instabilities
are in fact dominant. The growth rate of this instability isgoverned by the shear flow. It is larger than that of the typical resistive instability and smaller than that of the typicalmagnetic-field driven ideal instability. Therefore, we believethat it is important to investigate this instability carefully,
albeit in the simplest slab geometry.
We present our formulation and the general properties ofthe linear eigenvalue equation in Sec. II. Since the instability
is global and depends strongly on the details of the field
profile 共we focus on a linear flow profile兲 it is difficult tomake a quantitative discussion by analytic calculation. We
cannot discuss the instability in the infinite domain, since theinfinite difference of the shear flow at both edges prevents
any exponential instability to grow due to its stretching
共shearing兲 effect.12,13 On the other hand, a study of this instability in a periodic domain would inevitably introduce inflection points, which makes it difficult to distinguish thephenomena from KH instability.14 In practice, the overall
dispersion relation has never been obtained by analytic calculation, even in the simplest slab geometry. In light of these
difficulties, we decided to pursue mainly numerical analyses,
which we report here.
Several numerical solutions of the linear eigenvalue
equation are obtained by invoking the shooting code, and areshown in Sec. III. These analyses reveal the general conditions that maximize the growth rate: 共a兲 the phase velocity of
the Alfvén wave in the direction of the background mass
velocity around the edge region should have comparable, but
smaller, values than the local mass flow, 共b兲 the central valueof k · B should be finite and have different sign from bothedges, and thus, 共c兲 k · B should change sign twice in thedomain, where k denotes the wavenumber in the homogeneous direction. These basic features lead one to suspect thatthe qmin surface of an advanced tokamak may be susceptible
to this instability. We show a concrete example of a profilethat may be applicable to this situation.
Results from nonlinear simulations are presented in Sec.
IV. In order to simulate the incompressible motion of theplasma, we developed a pseudospectral code that solves twodimensional reduced MHD equations, in a vorticitystreamfunction and flux function formulation. The nonperiodicity of the problem leads to the choice of Chebyshev
polynomials for the basis functions. The Chebyshev spectral
scheme is nowadays a common tool to solve fluid dynamical
turbulence problems, including boundaries. However, it is
not as popular in plasma physics as it is in fluid mechanics.
The development of the combined reduced MHD-Chebyshev
formulation presented here may have wider applications, andis therefore presented in some detail. In the specific casestudied here, we find that nonlinear evolution of the above
mentioned tokamak example may lead to disruption with atime scale slower than Alfvén time and faster than the resistive one, which might account for the disruption recently
observed in JT-60U.9
Finally, the summary is given in Sec. V.

II. FORMULATION

We wish to consider incompressible plasma motion inslab geometry with equilibrium sheared flows. The governing equations are␳ 共 ⳵ tv + v · ⵱ v 兲 =1
共⵱ Ã B兲 Ã B − ⵱p,␮0

共1兲⳵tB − ⵱ Ã 共v Ã B兲 = 0,
⵱ · v = 0,

and

共2兲

⵱ · B = 0,

共3兲where the density ␳ is assumed constant.
Dividing all fields into equilibrium and perturbed parts,
the former 共latter兲 is expressed by the uppercase 共lowercase兲
letter in this and subsequent sections. We consider a general
one-dimensional slab equilibrium, including both paralleland perpendicular shear flows. The magnetic field and flow
profiles are given by
B = 关0,By共x兲,Bz共x兲兴,

共4兲

V = 关0,Vy共x兲,0兴,

共5兲

in Cartesian coordinates. The equilibrium pressure P共x兲 ischosen such that

冉

冊

d
B2
P+
= 0.
dx
2␮0

共6兲

The profiles of equilibrium magnetic and velocity fields may
be arbitrary functions of x. As we mentioned in Sec. I, weconcentrate on the linear shear flow profile 共plane Couetteflow兲 in this paper. Since the infinite size of the domain in
the x direction introduces the infinite difference of the background flow between both edges of the system,12,13 we set
our domain x 苸 关−L , L兴 and put a no-slip, ideal conducting
wall at x = ± L.
Upon linearizing the equation of motion, taking its curl,
and assuming the ansatz
f共x,y,z,t兲 = f共x兲exp关i共ky y + kzz − ␻t兲兴,共7兲

for all perturbed quantities, we obtain the eigenvalue ODE:
d
dx

冋冉

⍀2 −

冊 冉 冊册

F2 d vx␮0␳ dx ⍀冉

− 共k2y + kz2兲 ⍀2 −

冊

F2 vx
= 0,
␮ 0␳ ⍀
共8兲

where ⍀ = ␻ − k · V and F = k · B. The boundary condition introduced by the no-slip, ideal wall implies vx = 0 at x = ± L.
It is convenient to normalize 共8兲 by the half system size*:
L and characteristic poloidal Alfvén velocity vpA* 冑*= By / ␮0␳, where By is chosen such that By共x兲 is order unity:
x = Lx̂,

␻=

*
V̂y,
Vy = vpA

*
vpA
␻ˆ ,
L

1
k = k̂,
L

t=

L
*
vpA

t̂,

d 1 d
.
=
dx L dx̂

共9兲

共10兲

Here the values carrying ˆ denote the normalized variables.
The normalized ODE is

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-3

d
dx̂

冋共

Phys. Plasmas 13, 092107 共2006兲

Magneto-flow instability in symmetric field profiles⍀̂2 − F̂2兲d

冉 冊册 共v̂xdx̂ ⍀̂

− k̂2 ⍀̂2 − F̂2兲

v̂x

⍀̂

= 0,

共11兲

冑

where F̂ : = k̂ · B̂ and k̂ : = k̂2y + k̂z2. Hereafter, we omit the ˆ forthe normalized quantities, and the values mean the normalized ones unless we explicitly mention otherwise.
We now summarize the general properties of the ODE
共11兲.
If ky = 0, the effect of shear flow is excluded and the
ODE 共11兲 reduces to the one for static equilibria. Thus weknow the system is stable, and we do not consider this casehere.
Since the magnetic field appears only as a square in 共11兲,
the sign of By or Bz does not affect the instability.
If ␻ = ␻e is an eigenvalue for a certain velocity profile
Vy = Ve共x兲, then ␻ = ␻e + kyV0 becomes an eigenvalue for thevelocity profile Vy = V0 + Ve共x兲 with the same eigenfunction,
where V0 is a constant velocity, since these combinations
yield the same factor, ⍀ = ␻e − kyVe 共Doppler shift兲. Thus, we
may assume Vy共0兲 = 0 without loss of generality.
From the sufficient condition of stability given by several authors,3,15,16 it is clear that the system is stable if there
is a reference frame in which the local Alfvén velocity projected onto the mass flow direction is larger than the local
flow velocity everywhere in the domain. That is, in physicalunits,
兩Vy兩 ⱕ

冏 冏

k储vA
共"x兲,
ky

共12兲

where vA = 兩B 兩 / 冑␮0␳ so that k储vA = F. The inverse may not
be true. This means that the velocity Vy has to exceed thevalue of the Alfvén velocity 共directed along B兲 projected
onto the direction of Vy. Note that this condition does notnecessarily mean that Vy has to be larger than the Alfvén
velocity, defined by By. An example is given in Sec. III B.
Another caveat of condition 共12兲 is that it is not the magnitude of the flow itself but the flow shear that matters, as wasseen in the previous paragraph. On the other hand, it is themagnitude of the magnetic field that matters and not themagnetic shear only. Thus it is sufficient for stability if thecondition 共12兲 is satisfied in a particular reference frame.
It is important to realize that the direction of shear flow
does not need to be parallel to the ambient magnetic field.
When a field profile yields instability for a parallel flow, thenthe same profile with perpendicular flow will also lead to
instability if we may rescale the length. Since we derived alinearized ODE for the ambient field and flow 共4兲 and 共5兲, let
us fix the flow profile with V = 关0 , Vy共x兲 , 0兴 and change thedirection of field to check the effect of parallel/perpendicularflow with respect to the field. If one finds an instability for aparallel flow where B = 关0 , B共x兲 , 0兴 with a wavenumber
共ky , kz兲 = 共ke , 0兲, let its eigenvalue and eigenfunction be ␻eand ␸e 共=vx / ⍀兲. Then, the same field profile yields the instability for the wavenumber 共ky , kz兲 = 共ke , ke兲 / 冑2 with an eigenvalue ␻ = ␻e / 冑2 and eigenfunction ␸e in case of a perpendicular flow where B = 关0 , 0 , B共x兲兴. Thus we mainly
concentrate on the parallel flow in this paper, except in Sec.
III B, where we show an explicit example for the application

FIG. 1. Instability threshold and maximum ci for parabolic field profile
By = ax2.

of a tokamak. We also note that the longitudinal magnetic
field 共Bz兲 enters only through F and thus does not affect themode at all as far as kz = 0.

III. LINEAR INSTABILITY
A. Parallel flow

As was obtained by Chen and Morrison,5 we have plotted ci 共=limky→0 ␻i / ky兲 with respect to a in Fig. 1, where theflow profile is fixed at Vy共x兲 = x and a parallel magnetic fieldof the following form is assumed:
By共x兲 = ax2 .

共13兲In agreement with their analytical estimate, we find that thecritical value above which the instability occurs is a
⯝ 0.834. As long as a ⬎ 0.834, the system is unstable. However, for larger values of a 共a ⬍ 30; ci ⲏ 3.28⫻ 10−8兲, ournumerical results indicate that ci converges to zero withci ⯝ 0.027a−4. The dispersion relation is shown in Fig. 2 fora = 1, namely, for the almost largest ci. The growth rate has amaximum at finite ky 共⯝0.8 in this case兲, and the instability
is quenched for a larger ky, which is similar to the behavior
of KH instability. The real part of ␻ may be finite, in general,but turned out to be zero for all ky in the examples shown inthis paper.
In the following sections, we discuss the property of thisinstability by showing the dependence of the growth rate on
the various ambient field profiles.

FIG. 2. Dispersion relation for Vy共x兲 = x and By共x兲 = x2.

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-4

Phys. Plasmas 13, 092107 共2006兲

T. Tatsuno and W. Dorland

FIG. 3. The quantity ci and maxky ␻i for the field profile By共x兲 = 共1 − ␣1兲x2+ ␣ 1.

1. Dependence on central value of field

The instability remains when the magnetic field is modified to
By = 共1 − ␣1兲x2 + ␣1 ,共14兲

which changes the value of By共0兲 = ␣1 with By共±1兲 = 1 fixed.
The maximum growth rate maxky共␻i兲 and the quantity ci areplotted with respect to ␣1 in Fig. 3. Here the upper bound of
␣1 ⯝ 0.5 for instability comes from the fact that local Alfvén
velocity exceeds the local flow velocity everywhere in thedomain, which corresponds to 共12兲. For negative ␣1, we havelarger ci, as we see from Fig. 3. For ␣1 ⱖ −25, we find no
maximum of ci with respect to ␣1. As we decrease ␣1, however, the upper bound of ky for instability decreases while theci corresponding to ␻i / ky in the limit ky → 0 increases. With
the combined effect of them, the maximum growth rate isachieved at about ␣1 ⯝ −0.7 for ky ⯝ 0.6.

FIG. 5. Growth rate for the parabolic field profile and for 共15兲. Notice thatthe scales of ordinates differ between two cases.

兩Vy共x兲 兩 = 兩x兩. Note that the magnetic field at 兩x 兩 ⲏ 0.4 isweaker for the profile 共15兲 than the parabolic one, while
those in 兩x 兩 ⱗ 0.4 are very close to each other. The dispersion
relations for these two cases are shown in Fig. 5. Notice thatthe scale of the ordinate for 共15兲 is two orders of magnitudelarger than that for the parabolic profile, and so as the maximum growth rate. The profile with weaker magnetic shear in
the outer region yields a much larger growth rate than theparabolic case.
The reason becomes clear when we look at the eigenfunctions for ky = 0.1; see Fig. 6. As is seen from the figure,there are resonant and step-like structures around 兩x 兩 ⬃ 0.5
for the parabolic field profile, where 兩Vy 兩 = By. Since the magnetic field is so strong that the local Alfvén velocity exceeds
local flow velocity at 兩x 兩 ⬎ 0.5, the eigenfunction cannot extend to the edge and is rather confined within 兩x 兩 ⱗ 0.5. On
the other hand, the Alfvén velocity corresponding to 共15兲 is

2. Dependence on edge shear

Chen and Morrison5 mention in their paper that the“magnetic field at the boundaries tends to destabilize theshear flow.” This is certainly true for the analytically tractable case where By共x兲 is parabolic. To check the statement
for general profiles of By, we prepared a profile,
By = 1 −tanh2␣2 x

␣22x2

,

共15兲

to compare with the parabolic one. This profile is shown for␣2 = 2 in Fig. 4, where the dotted line shows the profile of

FIG. 4. Profile of the magnetic field and absolute value of the velocity.

FIG. 6. Eigenfunctions 共a兲 for parabolic profile and 共b兲 for 共15兲 withky = 0.1.

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-5

Phys. Plasmas 13, 092107 共2006兲

Magneto-flow instability in symmetric field profiles

slower than the local flow velocity in the whole domain,
therefore, the eigenfunction can be extended in the wholedomain.
Thus we find that a strong magnetic field whose Alfvén
velocity exceeds local flow velocity has a stabilizing effecton the mode. It is interesting to compare this result with thecase with an interchange-type drive.7 In Ref. 7, the shear
flow stretching effect exhibited a stabilizing effect only when
the flow velocity exceeded the Alfvén velocity. In the present
example, this is analogous to the finding that there is a minimum on the magnetic field strength below which the instability is quenched 共see Fig. 1兲. On the other hand, if themagnetic field is too strong compared to the flow velocity,
the interchange-type instability was stabilized, which was
expressed as the destabilizing effect of the shear flow in Ref.
7 from the opposite viewpoint. In the present case, this may
correspond to the stabilizing effect of the edge magnetic
shear in the parabolic field profile. The important point isthat when the effects of field and flow are comparable, themaximum growth rate is achieved.
3. Larger growth rate

To reach nonlinear saturation, the growth must occur
considerably faster than the decay of the background field,
since inhomogeneity of the magnetic field is one of thesource for the instability. For this purpose, we may point out
two effects that increase the growth rate, according to theresults obtained from previous two subsections: 共1兲 to makethe central By value negative, and 共2兲 to make the edge magnetic shear weak. Paying attention to these facts, we prepareda profile
1
1
␣3
By共x兲 = 共x2 + 1兲 −− 关e␤3共x−1兲 + e−␤3共x+1兲兴,2
cosh2共2x兲 ␤3
共16兲
where ␣3 and ␤3 are introduced so that the current density
vanishes at the edge, which is needed for the following resistive simulation. Here we take ␣3 ⯝ 1.27 and ␤3 = 10. Profile 共16兲 is shown in Fig. 7 together with the flow profile forcomparison, and the corresponding maximum growth rate isabout ␻i ⯝ 0.153 for ky ⯝ 1.
B. Application to tokamak qmin surface

When the velocity is everywhere parallel to magnetic
field 共as discussed in the previous section兲, the mass velocity
must be comparable to the Alfvén velocity for the instability
to occur. However, by adding a perpendicular field, we may
reduce the velocity required for instability to be muchsmaller than both of parallel and perpendicular componentsof the field. Here we here show an example of the tokamak
qmin surface in slab geometry.
Let us assume a parabolic safety factor profile:

冉 冊

aBT
,
RB P

共17兲

qmin = 1.7,

共18兲

q = x2 + qmin =with
R/a = 3,

FIG. 7. Profile 共16兲 that makes the growth rate larger, and the corresponding
dispersion relation.

共m,n兲 = 共2,1兲,
Bz = 5.0,

ky = 1,

共19兲
共20兲

where R and a denote the major and minor radius, respectively; thenV y = ␴x

共21兲

would yield an instability with a growth rate max␴ ␻i⯝ 0.016 for ␴ ⯝ 0.23. In this example, the instability is observed for 0.18ⱗ ␴ ⱗ 0.3. Notice that the growth rate ␻i is
one order of magnitude smaller than the flow shear ␴, which
is significantly smaller even than the poloidal Alfvén velocity.
Figure 8 shows the profile and eigenfunction obtained
from the shooting code for ␴ = 0.23. As is clearly seen in Fig.
8, the velocity is everywhere smaller than both the local
poloidal 共y component兲 and toroidal 共z component兲 field,
however, 兩k · V兩 and 兩k · B兩 have a similar magnitude and thesufficient condition 共12兲 is violated. A somewhat peaked
structure is observed in the eigenfunction around x ⯝ ± 0.4,
corresponding to the resonant radius where 兩k · V 兩 ⯝ 兩k · B兩.
The important point is that, as far as the slab model is concerned, the instability can take place if the difference of flow
is comparable to the difference of poloidal Alfvén velocity inthe region of the sheared velocity.
The sufficient condition for stability given by 共12兲 takes
a simple form when we restrict ourselves to a linear shear
flow and symmetric parabolic field profiles, respectively,

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-6

Phys. Plasmas 13, 092107 共2006兲

T. Tatsuno and W. Dorland

IV. NUMERICAL SIMULATION

In order to investigate the nonlinear behavior of the instability, we have developed a two-dimensional pseudospectral reduced MHD code, which solves

⳵t ⌬␾ + 兵␾,⌬␾其 − 兵␺,⌬␺其 = ␯ ⌬2␾ ,共24兲

⳵ t␺ + 兵 ␾ , ␺ 其 = ␩ ⌬ ␺ ,共25兲

where 兵P , Q其 = 共⳵y P兲共⳵xQ兲 − 共⳵x P兲共⳵yQ兲 is the standard
Poisson bracket, ␯ and ␩ are 共normalized兲 viscosity and resistivity, and ␾ and ␺ are the streamfunction and flux function, related to vector fields asv = ⵱ ␾ Ã e z,

B = ⵱␺ Ã ez ,共26兲respectively. We impose no-slip and perfectly conducting
boundary conditions on x = ± 1, namely,
vx = 0,

FIG. 8. Background field and velocity profile and eigenfunction vx corresponding to ky = 1.

冉V⬘y 2 ⬍ 4B⬙y 共0兲 By共0兲 +

冊

kz
Bz ,
ky

共22兲

which in turn yields the following stability condition in thephysical units:

冉 冊
dVy
dx

2

2

2
⬍ 4vpA

冉

冊

1
dq n
−
,
2
dx m qmin

共23兲2= B2y 共0兲 / ␮0␳ is the local poloidal Alfvén velocity,where vpA
and we assumed a uniform toroidal field.
We must address two cautions on the stability condition
共23兲. First, 共23兲 is formulated as a local condition under theassumption of linear flow and symmetric parabolic field profiles, while this instability is essentially a global one. Since
the smaller magnetic shear in the edge region brought about
larger growth rate as we saw in Sec. III, a more realistic field
profile might make the equilibrium more unstable thanevaluated from this model. Second, it is derived from thesufficient condition for stability. It is clear that whenever q
crosses a rational number 关or 共nqmin − m兲 becomes negative兴,the condition 共23兲 is violated, however, there are many stable
examples that may violate 共23兲, even for the symmetric parabolic field profile 共note the lower limit of B⬙y obtained byChen and Morrison5兲. Besides, one has to restrict the application to a low m/n mode since this instability does not arise
for a too large ky. Nevertheless, this condition may be usefulto get a rough idea on the regime of the instability. If V⬘y is
far smaller than the right-hand side, then the equilibrium
may be stable for this instability, but if they are comparable,
one may have to check the stability numerically for that specific profiles.

vy = const,

Bx = 0,共27兲

jz = 0,and periodicity in the y direction. Note that we cannot use
free-slip 共⌬␾ = 0兲 boundary condition in order to sustain linear profile of the background shear flow. In principle, finite
resistivity introduces the possibility of a tearing mode instability. We have carefully checked that the observed growthrate and nonlinear evolution show no correspondence to tearing mode physics 共see later兲.
For tokamak applications, this code can make singlehelicity calculation of Strauss’ reduced MHD equations,17
2
2
2
2
2
⳵t ⵜ⬜
U + v · ⵱ ⵜ⬜
U = B · ⵱ ⵜ⬜
A + ␯⵱⬜
ⵜ⬜
U,

共28兲

2
⳵tA = B · ⵱U + ␩ ⵜ⬜
A,

共29兲

by changing the definition of the Laplacian, Poisson bracket,
and flux function to2= ⳵2x + ⳵␨2 ,
ⵜ⬜

共30兲

兵P,Q其 = ⳵␨ P ⳵xQ − ⳵x P ⳵␨Q,

共31兲

␺ = A − ␣Bzx,

共32兲␾ = U,

共33兲

where ␨ = y + ␣z is the helical coordinate. The singlehelicity calculation may be justified because there is only asingle set of unstable low wavenumbers in the example
shown later. The inclusion of another dimension should not
significantly alter the results.
18

A. Numerical scheme

In order to impose nonperiodic boundary conditions in
the x direction, we used Chebyshev polynomials as basis
functions. This admits the use of the fast Fourier transformdue to the property
Tn共cos ␪兲 = cos共n␪兲,

共34兲

where Tn is the nth order Chebyshev polynomial. It is notedthat Chebyshev polynomial expansion has an exponential
convergence for any smooth function defined on a finite nonperiodic domain while the Fourier is algebraic in this case.19

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-7

Phys. Plasmas 13, 092107 共2006兲

Magneto-flow instability in symmetric field profiles

We use Orszag’s two-thirds rule20 to remove the aliasing
effect from nonlinear terms. One disadvantage of using
Chebyshev polynomials may be the accumulation of theround-off errors in calculating derivatives,21 however, thisshould not invalidate our calculation in the range of mode
numbers we take in our calculations 共2N / 3 ⱕ 1152, whereN + 1 is the number of grid points in x兲. The periodic direction 共y兲 is expanded by the Fourier series, as usual.
Time integration is conducted by the third-order AdamsBashforth scheme for nonlinear terms and the CrankNicholson scheme for dissipative ones.22 The time step ischosen at each step so that it avoids numerical instability.
See Sec. A 1 for the AB3 formula for an uneven time step.
The reason we used a different scheme depending on linear/
nonlinear terms is basically because of the numerical stability. The dissipative terms bring the highest order restrictions
on time steps in connection with the spatial grid, which is of
O共⌬x2兲, where ⌬x is the smallest grid size. In the Chebyshev
pseudospectral scheme, we use the Chebyshev-Lobatto
共extrema兲 grid, which is concentrated at the edge as⌬x ⬃ O共1 / N2兲. Thus the time step allowed by an explicit
scheme would be extremely small 关of the order O共1 / N4兲兴,which is impractical. On the other hand, using the CrankNicholson scheme for nonlinear terms is not simple since
one needs to solve the nonlinear boundary value problem at
each time step, which would require iteration of the matrix
inversion.
The boundary conditions written in 共27兲 are on velocity
and magnetic fields for which we are not directly solving
equations. Equations 共24兲 and 共25兲 demand two boundary
conditions on vorticity ⌬␾ and flux function ␺ at each time
step, and two on streamfunction ␾ when solving the Poisson
equation. By converting the vy condition with the method of
Coutsias et al.,23 we imposed, in the nth Fourier space,␾n = 0,each ky. On the other hand the integral operation is describedby the band-diagonal matrix due to the three-term recursion
relation:冕

Tn共x兲dx =Tn−1共x兲
Tn+1共x兲
−+ C共n ⱖ 2兲.
2共n + 1兲 2共n − 1兲

Thus, by multiplying the integral matrix prior to inversion,
one may carry out all matrix inversions by O共N兲 operations
for each ky. A more detailed description is given in Sec. A 2.
Boundary conditions are imposed by the Sherman-Morrison
formula.25
As is often referred in spectral simulations, the FourierChebyshev amplitude corresponding to the highest wavenumber component may be regarded as an empirical measureof the numerical error. Using this fact, we may automaticallycontrol the resolution, which greatly enhances the performance in linear instability and/or decaying turbulence simulations. The code checks the ratio between the amplitude ofthe largest wavenumber component and that of the largest
component of the vorticity in the Fourier-Chebyshev space.
When it becomes larger than ulim 共smaller than llim兲, the resolution is enhanced 共reduced兲. We have made these values
adjustable, and the resolution of exact power of 2 is automatically assigned, while other mixed-radix grid points are
given as inputs.
B. Linear instability and its quasilinear effectWe first show the linear instability of the profiles:
Vy共x兲 = x,

共41兲

1
1
␣3
− 关e␤3共x−1兲 + e−␤3共x+1兲兴,
By共x兲 = 共x2 + 1兲 −22cosh 共2x兲 ␤3
共42兲

共35兲

兺m Bmn共⌬␾兲mn = 0,

共36兲

at x = ± 1 for ky ⫽ 0 共n ⫽ 0兲, and⳵x␾0 = const

at

x = 1,

共37兲

⳵x共⌬␾兲0 = 0at

x = ± 1,

共38兲for ky = 0 共n = 0兲, where 共⌬␾兲mn denotes the mth Chebyshev
and nth Fourier amplitude of the vorticity and Bmn denotes
the numerical matrix described in Ref. 23. The arbitrarinessof constant for the potential ␾ is determined by ␾00 = 0. The
boundary condition on magnetic field 共and current density兲 isreadily converted to the one on ␺:␺ = const 共in time and along each wall兲.

共39兲

We use a preconditioner for the inversion of the matrices
in solving Poisson equation and Crank-Nicholson time stepping, as introduced by Coutsias et al.24 Since the derivative
matrix for Chebyshev polynomials is upper-triangular, thesimple inversion of the Laplacian needs O共N2兲 operation for

共40兲

which was introduced in Sec. III A 3. As we discussed before, this particular magnetic field profile was chosen to
maximize the growth rate large, since this instability has arelatively small intrinsic growth rate. The parameters ␣3 and␤3 are introduced in order for the current density to vanish at
the edge, and ␣3 ⯝ 1.27 and ␤3 = 10 in the current run, respectively. We take Ly = 2␲ since ky ⯝ 1 gives the largest growthrate ␻i ⯝ 0.153 in this case.
Figure 9 shows the time evolution of the L2 norm of
perturbed streamfunction 储␾1储, enstrophy, and total energy
when we put the initial perturbation of the form
共⌬␾兲1 =

冉冊

10−2
2␲共x + y兲sin,2cosh 共10x兲
Ly

共43兲

with various diffusivity, where the subscript 1 denotes theperturbed quantity. The enstrophy and energy are normalized
by their initial values. The maximum numbers of grid points
used in these simulations are summarized in Table I for a
sufficient convergence. The number of grid points drasticallyincreases as ␯ = ␩ becomes smaller than 3 ⫻ 10−4, which suggests a qualitative difference between case 共a兲 ␯ = ␩ ⱕ 2⫻ 10−4 and case 共b兲 ␯ = ␩ ⱖ 3 ⫻ 10−4. Figure 9 also shows acritical difference between cases 共a兲 and 共b兲: The 储␾1储 plot

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-8

Phys. Plasmas 13, 092107 共2006兲

T. Tatsuno and W. Dorland

FIG. 9. Time evolutions of 储␾1储, enstrophy, and total energy in the simulations with various diffusivity andnumber of grid points for the case with␯ = ␩ = 2 ⫻ 10−4.

almost overlaps and the enstrophy increases drastically aftert ⲏ 80 for 共a兲, and the energy dissipation rate jumps between
共a兲 and 共b兲 after t ⲏ 100. Case 共b兲 is thought to be dissipation
dominated, and the nonlinear evolution of the instability does
not show interesting behavior. However, as diffusivity goessmaller, the nonlinear behavior becomes turbulent, and theflux surface is completely broken in a wide region. In thisand subsequent sections, we focus on the case with ␯ = ␩= 2 ⫻ 10−4 since other simulations in case 共a兲 show qualitatively similar behavior.
The straight line in the evolution of 兩兩␾1 兩 兩 shows the
exponential growth, with growth rate 0.15. We notice that thegrowth rate of the simulation agrees well with that obtained
from the shooting code, and it approaches the ideal value
0.153 共for ␯ = ␩ = 0兲 as the dissipation rate becomes smaller.
Since the nonlinear evolution of the instability leads to a
decaying turbulence simulation, it is necessary to consider
the diffusion rate of the background field. It is governed by
the decay of the eigenfunction of ⳵2x with the smallest wavenumber. Currently we have a diffusion equation for ␺, which
is an odd function, under the fixed boundary condition. Thus
the smallest wavenumber eigenfunction corresponding to the
profile 共42兲 is sin共␲x兲 since Lx = 2, which in this case gives

rise to the decay rate of the magnetic field of ␲2␩. The e-foldtime of the background magnetic field is, thus, estimatedabout t ⬃ 500 for ␩ = 2 ⫻ 10−4, and by that time, the instability will completely lose its drive by the diffusion of the background magnetic field.
Let us see the quasilinear effect of the instability. It is
obtained by substituting linear eigenmode into nonlinear
terms and by calculating its backreaction on the backgroundfield. Here we will work on the original vector form equationin order to see the effect on vy and By:

⳵tv0 = − 共v1 · ⵱v1兲0 + 关共⵱ Ã B1兲 ⫻ B1兴0 ,共44兲⳵tB0 = 关⵱ Ã 共v1 Ã B1兲兴0 ,共45兲

where the subscripts 0 and 1 denote the background 共ky = 0兲and linearly unstable 共ky = 1兲 components of the fields,
respectively.
By substituting
˜ Ã e z,v1 = ⵱␾

B1 = ⵱˜␺ Ã ez ,

共46兲

˜ =␾
˜ r cos共ky y兲 − ␾
˜ i sin共ky y兲,␾共47兲˜␺ = ˜␺ cos共k y兲 − ˜␺ sin共k y兲,ryiy

TABLE I. Maximum numbers of grid points in the simulations.

into 共44兲 and 共45兲 and taking the ky = 0 component of the rhs,
we obtain

case

␯=␩

Grids 共x ⫻ y兲

ulim

共a兲

1.e-4
1.5e-42.e-41728⫻ 4096
1024⫻ 2304
768⫻ 1280

5.e-6
1.e-5
1.e-4

共v1 · ⵱v1y兲0 = −ky ˜ 2 ˜ ˜ 2 ˜
共␾r ⳵x ␾i − ␾i ⳵x ␾r兲,
2

共48兲

共b兲

3.e-4
4.e-4
6.e-4

384⫻ 384
256⫻ 256128⫻ 128

1.e-3
1.e-3
1.e-3

共B1 · ⵱B1y兲0 = −k y ˜ 2˜ ˜ 2˜
共␺r ⳵x ␺i − ␺i ⳵x ␺r兲,
2

共49兲

since

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-9

Phys. Plasmas 13, 092107 共2006兲

Magneto-flow instability in symmetric field profiles关共⵱ Ã B1兲 Ã B1兴0 = − 21 关⵱共B21兲兴0 + 共B1 · ⵱B1兲0 ,共50兲

and the gradient term gives rise to no contribution on theevolution of the y-component field the same as the ⵱p term
共ky = 0兲. By making the same manipulation on the rhs of 共45兲,
we obtain关⵱ Ã 共v1 Ã B1兲兴0y =ky 2 ˜ ˜ ˜ ˜⳵ 共␾r␺i − ␾i␺r兲.
2 x

共51兲

Since the linear eigenfunction is growing with thegrowth rate ␻i, the evolution equations 共44兲 and 共45兲 yield
˜ r ⳵ 2␾˜ ˜ 2 ˜ ˜ 2˜ ˜ 2˜
⳵tv0y = Ce2␻it共␾
x i − ␾i ⳵x ␾r − ␺r ⳵x ␺i + ␺i ⳵x ␺r兲,共52兲
˜ r˜␺i − ␾˜ i˜␺r兲,⳵tB0y = Ce2␻it⳵2x 共␾

共53兲

where C is a constant proportional to the square of the initial
amplitude of the eigenmode. Note that temporal oscillationof the eigenfunction 共if there is any兲 cancels out during thecancellation of ky y.
Figure 10共a兲 is a plot of the rhs of 共52兲 and 共53兲, obtained from the numerical eigenfunction by means of theshooting method. We take 401 evenly spaced grid points andderivatives are calculated by means of the second order centered finite difference formula. Figure 10共b兲 shows the difference between the initial fields and solutions of nonlinear
ideal 共␯ = ␩ = 0兲 simulation at t = 25, where ky ⬎ 1 components
are artificially killed during evolution. The difference from共a兲 thus comes from the stable ky = 1 oscillation that is excited initially and the self-consistent change of the eigenfunction due to the change of the background field. These
two figures almost overlap with each other. Figure 10共c兲
shows the difference between solutions of full nonlinear
simulation and that of the one-dimensional 共1-D兲 linear diffusion equation for ␺ 共with the same ␩ and grid points兲 at
t = 25, which includes the effect of dissipation on linear instability and all of the nonlinear effect including the backreaction from higher harmonics. The peaks are broadened
somewhat due to the effect of dissipation but the qualitative
agreement is good.
Here, By = 0 initially takes place around x ⬃ ± 0.39. We
can see that the quasilinear effect of the instability on theequilibrium affects the velocity and magnetic fields roughly
equally. The energy ratio of the quasilinear magnetic/kinetic
field modification is about 1.22 for the shooting result. 共It
affects more on the magnetic field in the energy norm while
the maximum is larger for the velocity in Fig. 10.兲 The velocity is flattened and magnetic field is weakened around theorigin. As can be seen from the figure, most of the effecttakes place where By ⬍ 0, and the area where By is negativeis slowly widened because of the fact that ⳵tB0y ⬍ 0 at
x ⬃ ± 0.39. This is necessary to push up the central value of
By since we are imposing the constraint that 兰By dx = const
by placing the ideally conducting wall.

FIG. 10. Quaslinear effect of the instability: 共a兲 evaluated from 共52兲 and共53兲 with a shooting result, 共b兲 at t = 25 in the ideal simulation with ky ⬎ 1
artificially killed, and 共c兲 at t = 25 in the dissipative simulation with full
nonlinearity.

C. Nonlinear evolution

We now turn to the nonlinear effect of the instability.
Time traces of the y-averaged velocity and magnetic fields
are shown in Fig. 11. As is seen from the figure, the instability brings the positive By region just outside of By ⬃ 0 tolower values, and at the same time it brings up the central
By共0兲, which was pointed out in the previous section. The
instability does not alter the profile of velocity and field in
the outer region 兩x 兩 ⲏ 0.6, and it flattens and makes themapproach zero only in the central region. Note that the flattened flow profile looks as x3, which is KH stable, albeit withan inflection point. The change of the By field at the edge ismostly attributed to the diffusion of the background magnetic
field. According to an independent 1-D simulation of thelinear diffusion equation for ␺, the edge discrepancy turned
out to be within a few percent until t ⱕ 130.
The contour plots of streamfunction ␾ and flux function
␺ at t = 45 are shown in Fig. 12. Stretching and rolling up of
the flux function is observed around two resonant surfaces,

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-10

T. Tatsuno and W. Dorland

Phys. Plasmas 13, 092107 共2006兲

FIG. 11. Time evolution of the average flow and field.

qualitatively similar to a KH instability. The stage looks
rather coherent and the effect of instability is restrictedwithin 兩x 兩 ⱗ 0.5.
By t ⱗ 50, the instability is strongly nonlinear. The sizeof the vortex keeps on growing and vorticity concentration ispushed outward yielding two high-vorticity lines around x⯝ ± 0.5. About the first maximum of the enstrophy 共t ⯝ 57兲,these two high vorticity lines are broken into vortex chains,
as seen from Fig. 13. This stage from the first peak 共t ⯝ 57兲

FIG. 12. Contour plots of the streamfunction ␾ and flux function ␺ at
t = 45.

FIG. 13. Contour plots of the vorticity at t = 50 and 60, and the flux function
at t = 60.

to the second local minimum 共t ⯝ 84兲 of the enstrophy looks
somewhat turbulent. At this stage, however, the turbulent
structure is localized in two regions separated by a central
magnetic field line. Reconnection takes place easily aroundeach separated turbulent region since the equilibrium magnetic field line reverses along x ⯝ ± 0.4. However, this is inhibited at the center.
After t ⯝ 60, the vortex chains start to merge in eachregion and try to entangle the central horizontal magnetic
field line in a clockwise manner. The entanglement of the
field line once reaches maximum around t ⯝ 75, and the field
line reflects these vortex chains due to its tension force 共seeFig. 14兲. The reflected vortices merge further in each region
and form two large vortices of the same polarity. These twolarge vortices entangle the central field line again, and reconnection starts around top and bottom regions along y ⯝ ␲.
At t ⯝ 100, there are two large vortices with the samepolarities 共centered around y ⬃ 0.6␲ and 1.4␲兲 and a small
one with an opposite polarity 共centered around y ⬃ 0 or 2␲兲;
see Fig. 15. Because of the same polarity, a high-vorticity
line is formed and squeezed between them as their sizes are
getting large, which causes the increase of the resolution in
the y direction 共see Fig. 9兲. At this time, the ␺ contour shows

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-11

Phys. Plasmas 13, 092107 共2006兲

Magneto-flow instability in symmetric field profiles

FIG. 14. Entanglement of the central field line.

two large islands extending to 兩x 兩 ⱗ 0.8, indicating a large
energy transport along the field line. Strong current sheets
are observed along horizontal lines x ⯝ ± 0.8 and along a
vertical line y ⯝ ␲. However, the two islands have the opposite polarity, and magnetic reconnection is mainly observed
along horizontal current sheets. It is noted that the second
reconnection event does not take place for smaller diffusivity
关case 共b兲兴 so that the central field line separates the two turbulent regions all the time. However, the critical diffusivity
needed for the reconnection of the central field line is a
subtle issue depending on the amplitude of the initial perturbation, since the background field profile decays away before
the appearance of a turbulent state when one puts a tiny
perturbation in the initial condition. Recall that the decay
time of the background field profile is about 500 for ␩ = 2⫻ 10−4. The average velocity and magnetic field at t = 100 is
flattened in a fairly wide range of the domain, and is close tozero in 兩x 兩 ⱗ 0.5.
At the final peak of the enstrophy 共t ⯝ 113兲, the highvorticity line is broken at the center of the domain, and thetwo vortices start to merge. The enstrophy monotonicallydecreases while the large merged vortex absorbs the smaller
one with the opposite polarity, and the weaker vortex remains at t = 180. The average By is further pushed upward,

and it becomes positive almost everywhere in the domain by
t = 130. The average flow is also flattened in a wider region.
We may apply this result to the tokamak qmin surface.18
In this case, the two-dimensional simulation is regarded as asingle helicity calculation, and By = 0 corresponds to the resonant surface of the corresponding mode. Since the tokamak
has a strong toroidal field, it is convenient to consider anormalization to the velocity. In this case, the B*y used for thenormalization in Sec. II is regarded as a value whose corresponding Alfvén velocity gives a comparable magnitude of
the surrounding mass flow velocity. It may be significantly
smaller than both the poloidal and toroidal fields.
Since the sign of By does not matter on the results of theinstability, we may regard the minimum of By as a qmin surface by flipping the By profile. The change of the sign of By
corresponds to the fact that q crosses a rational value for themode we are interested in. The time evolution of the average
field suggests that we may have some difficulty in making q
cross the rational 共or integral兲 value. Moreover, it may leadto a major 共or minor兲 disruption according to the large scale
rearrangement of vy and By 共Note that we are using no-slip
and ideal-wall boundary conditions; the ideal wall implies
that the spatial integral of average By has to be conserved,
which is, in some sense, suppressing the instability.兲
Of course, the occurrence of disruption should be tested
in a larger simulation with a more realistic geometry 共in thiscase flow profile may contain an inflection point兲. We cananticipate what should be observed from our analysis. The
frequency of the instability is zero here since we assumed
symmetry in the flow profile, but it may be finite if the flow
profile is asymmetric. As is easily understood when one
makes a Galilean transform so that central vy is finite, there
appears a finite frequency arising from the Doppler shift, and
it is of the order of the background flow transit frequency.
Because of the global structure of the instability, the fluctuation would appear in low mode numbers. The flux function
will be stretched and rolled up around two rational surfaces,
as is seen in Fig. 12. On the other hand, the diamagnetic
effect might be able to stabilize it in the presence of the
background pressure gradient as is the case of KH
instability.2
V. CONCLUSION

The linear instability incurred by shear flow without an
inflection point is investigated in slab geometry by numerically solving linear eigenvalue equation and by performing a
two-dimensional 共or single-helicity兲 nonlinear simulation.
We fixed the shear flow profile with a spatially linear
function, and mainly focused on the effect of various symmetric field profiles. From the investigation using the shooting code, the typical field profile that gives a large growthrate has a few features: 共a兲 the phase velocity of the Alfvén
wave in the direction of the background mass velocity in theouter region has comparable, but smaller, values than thelocal mass flow; 共b兲 the central value of F 共=k · B兲 can be
finite and has a different sign from both edges, and thus, 共c兲
F changes sign twice in the domain. By adding a component
in the field nonparallel to the flow, we also showed an ex-

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-12

Phys. Plasmas 13, 092107 共2006兲

T. Tatsuno and W. Dorland

FIG. 15. Contour plots of the streamfunction ␾, vorticity ⌬␾, flux function
␺, and current density ⌬␺ at t = 100.

ample that may be relevant to the qmin surface of an advanced tokamak. Notice that 共a兲 does not imply that mass
velocity must exceed local Alfvén velocity. It can be significantly smaller than both toroidal and poloidal Alfvén velocity of a tokamak for the instability to occur.
In the numerical simulations, the growth rate is properly
observed. The quasilinear effect of the instability is to pushback the central value of the magnetic field in the direction
where two zeros of F disappears, but the instability initiallysaturates before the central F changes sign. A weak turbulent
state follows in the nonlinear stage, and after the sequence ofvortex breakup/merges, background flow and field profiles
are flattened in a fairly wide region. The negative region of F
finally disappears due to their nonlinear interactions. The
field line topology changes drastically, suggesting a huge energy transport across negative/positive regions of the magnetic shear, and the negative region of F finally disappears.
When it is applied to the tokamak qmin surface, the disappearance of the negative-F region means the difficulty ofqmin crossing the integral 共or rational兲 value, and the large
scale rearrangement of v and B may imply the disruption.
To understand the role of this instability in tokamak disruptions, it would be necessary to consider more realistic
geometry. There are effects that could significantly alter thegrowth rate and the threshold of the instability. First may be
the existence of an inflection point, and others would include
magnetic curvature and kinetic effects. One may obtain amore realistic condition of instability than 共23兲 by adding
more physics. However, the nonlinear evolution of the instability shown here would be valid once this instability takes
place. As we have seen the breakup of the flux surface andflow profile in the nonlinear evolution, especially at the central region where the transport barrier is thought to be located, MHD instability may be considered a possible candidate of disruption when flow is included.
In the recent experiment in JT-60U, an intermediate timescale instability is observed in the reversed shear discharges,

which is too fast as a resistive instability, and too slow as an
ideal one.9 The instability reported in this paper is a goodcandidate for explaining these observations, as it evolves on
the corresponding 共intermediate兲 time scale. It is driven by
the shear flow, and its characteristic growth time can even be
one order of magnitude longer than the poloidal flow transit
time. We have concentrated on the case without an inflection
point to focus on the simplest dynamics of magneto-flow
instability in this paper, however, it is also true for the casewith the inflection point that the mass flow with neither oftoroidal nor poloidal Alfvén velocity is needed for the instability. The inclusion of inflection point may make the equilibrium more unstable and the growth rate may increase accordingly, however, it is still true that the growth rate isgoverned by the flow, whose inverse is typically muchslower than any Alfvén time. As the importance of shear flow
effect on resistive instability is pointed out in Ref. 26, that in
ideal MHD also needs to be investigated more thoroughly
since a considerable amount of shear flow is expected in theadvanced tokamak.11 Moreover, it may also be applicable toinstabilities observed in the tokamak edge, if the bootstrap
and diamagnetic currents are strong enough to reverse themagnetic shear locally.
ACKNOWLEDGMENTS

We acknowledge fruitful discussions with Dr. Zensho
Yoshida, Dr. Greg Hammett, Dr. Frank Waelbroeck, and Dr.
Swadesh Mahajan.
This research is in part supported by the U.S. DOE Center for Multiscale Plasma Dynamics. A part of the simulations are conducted at SX-7 in the National Institute for
Fusion Sciences, Japan.
APPENDIX: NUMERICAL ALGORITHM
1. AB3 formula for uneven time step

Let the differential equation be

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-13

⳵y
= f.
⳵t

共A1兲

The Adams-Bashforth formula approximates the right-handside by the polynomial f̃共t兲 from a finite number of sample
points, and integrate the extrapolated polynomial to obtain yat the next step:
y 共n+1兲 = y 共n兲 +冕

t共n+1兲

共A2兲

f̃共t兲dt,t共n兲

where the superscript 共n兲 denotes the nth time step. We may
obtain it for the uneven time step by using Lagrangian polynomial for f̃共t兲. Here we show the formula for the third orderscheme, which uses a quadratic polynomial for f̃共t兲:f̃共t兲 =

共t − t共n−1兲兲共t − t共n−2兲兲 共n兲 共t − t共n兲兲共t − t共n−2兲兲 共n−1兲
ff −⌬1共⌬1 + ⌬2兲
⌬ 1⌬ 2+

共t − t共n兲兲共t − t共n−1兲兲 共n−2兲
f,⌬2共⌬1 + ⌬2兲

共A3兲

⌬0 = t共n+1兲 − t共n兲,⌬2 = t共n−1兲 − t共n−2兲 .

Substituting 共A3兲 into 共A2兲 and performing the integration, we obtainy 共n+1兲 = y 共n兲 + C0 f 共n兲 + C1 f 共n−1兲 + C2 f 共n−2兲 ,冋冉1
共⌬0 + ⌬1兲2 ⌬0 + ⌬1 ⌬2+
3
2⌬1 + ⌬2
⌬1
− ⌬1

冉

⌬1 ⌬2
+
3
2

冊册

共A5兲

冊
共A6兲

,

冉

冊

⌬20 ⌬0 ⌬1 + ⌬2,+
C1 = −
2⌬ 1⌬ 2 3
⌬20

冉

共A7兲

冊

共A8兲

Of course, they converge to the normal AB3 coefficients foran even time step,
23
⌬,
12

4
C1 → − ⌬,3C2 →f ⬘共x兲 = 兺 bmTm共x兲.

共A10兲

m

Then the integral operation on Chebyshev amplitudes is expressed as−1
al = 兺 Dlm
bm ,

共A11兲

m

where D−1 is a tridiagonal matrix defined from the integral
three-term recursion relation 共40兲:冕

T j dx = 兺 TiD−1
ij .

共A12兲

i

The first row of D−1 is arbitrary and is usually determined by
the boundary condition. Here we define the matrix 关1兴D−1 inwhich the first row of D−1 is replaced by zero. When we
define D by
bl = 兺 Dlmam ,

共A13兲

m

D−2D2 = 关2兴I共⫽D2关2兴D−2兲.

5
⌬,
12

共A9兲

in the limit ⌬0,1,2 → ⌬.

2. Preconditioner and time stepping

Let the Chebyshev expansion of a smooth function f共x兲defined on x 苸 关−1 , 1兴 and its first derivative f ⬘共x兲 be

共A14兲

Since 关1兴D−1 is a tridiagonal matrix with all diagonal components zero, 关2兴D−2 has nonzero components only at diagonal
and the second super-/subdiagonal components. Thus wemay divide 关2兴D−2 into odd and even parts to make two tridiagonal matrices of half size. We also note that all boundary
conditions used here can be divided into odd/even parts without loss of generality.
At each time stepping of the vorticity equation, we firsthave to find streamfunction ␾ from vorticity ␹ 共=⌬␾; notethat the sign is flipped in the definition of ␹兲 by solving the
Poisson equation:
共D2 − k2y,n兲␾n = ␹n

共A15兲

for each Fourier mode n 共in the y direction兲. Instead of inverting the upper-triangular matrix 共D2 − k2y,n兲, we multiply
the second integral operator 关2兴D−2 on both sides and obtain
共 关2兴I − k2y,n 关2兴D−2兲␾n = 关2兴D−2␹n .

⌬0 ⌬1
+
.
C2 =
2
⌬2共⌬1 + ⌬2兲 3

C0 →

l

关2兴

⌬1 = t共n兲 − t共n−1兲,

共A4兲

where

f共x兲 = 兺 alTl共x兲,and 关2兴D−2 = 关2兴D−1共 关1兴D−1兲2, we findwhere we have definedC0 =Phys. Plasmas 13, 092107 共2006兲

Magneto-flow instability in symmetric field profiles

共A16兲

We may replace the first two rows of 共A16兲 by two boundary
conditions, divide it into odd/even parts, and invert the resulting operator on the left-hand side. All matrices appearing
in 共A16兲 are tridiagonal when divided into odd/even parts
共plus one dense row corresponding to the boundary conditions兲, therefore, one can solve it by O共N兲 operations. Besides, the operator 共 关2兴I − k2y,n 关2兴D−2兲 is much better conditioned than 共D2 − ky,n兲: For details about the condition
numbers, see Ref. 24.
The same preconditioner is used for Crank-Nicholson
time integration. The reduced MHD equations 共24兲 and 共25兲basically consist of two parts that are convective 共nonlinear兲terms described by the Poisson bracket and diffusive 共linear兲terms described by the Laplacian. They are symbolically
written as

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp

092107-14

Phys. Plasmas 13, 092107 共2006兲

T. Tatsuno and W. Dorland

⳵y= Ny + Ly,
⳵t

共A17兲

where N and L denote nonlinear and linear terms, respectively; especially, L ⬀ 共D2 − k2y,n兲 for the n th Fourier mode.
As is mentioned in Sec. IV A, we use Adams-Bashforth for
nonlinear, and Crank-Nicholson for linear terms, so the discretization would yield

冉

1−

冊

⌬t
L y 共n+1兲 = y 共n兲 +2

冕t共n+1兲

t共n兲

dt + ⌬t Ly 共n兲 ,
Ny
2

共A18兲

 denotes the extrapolated polynomial described inwhere Ny
the previous section. Instead of inverting the upper-triangularmatrix 共1 − ⌬t L / 2兲 on the left-hand side, we first multiply
关2兴 −2
D on both sides from left and then only the band matrix
remains to be inverted.
J. W. S. Rayleigh, Proc. London Math. Soc. 9, 57 共1880兲.
B. N. Rogers and W. Dorland, Phys. Plasmas 12, 062511 共2005兲.
3
M. E. Stern, Phys. Fluids 6, 636 共1963兲.
4
A. Kent, Phys. Fluids 9, 1286 共1966兲; J. Plasma Phys. 2, 543 共1968兲.
5
X. L. Chen and P. J. Morrison, Phys. Fluids B 3, 863 共1991兲.
6
E. P. Velikhov, Sov. Phys. JETP 36, 995 共1959兲; S. Chandrasekhar, Proc.
1
2

Natl. Acad. Sci. U.S.A. 46, 253 共1960兲; S. A. Balbus and J. F. Hawley,
Astrophys. J. 376, 214 共1991兲.
7
T. Tatsuno, Z. Yoshida, and S. M. Mahajan, Phys. Plasmas 10, 2278
共2003兲.
8
P. A. Gilman and P. A. Fox, Astrophys. J. 484, 439 共1997兲; 510, 1018
共1999兲; 522, 1167 共1999兲.
9
S. Takeji, S. Tokuda, T. Fujita et al., Nucl. Fusion 42, 5 共2002兲.
10
Y. Miura and the JT-60 Team, Phys. Plasmas 10, 1809 共2003兲.
11
Y. Sakamoto, T. Suzuki, S. Ide et al., Nucl. Fusion 44, 876 共2004兲.
12
F. Volponi, Z. Yoshida, and T. Tatsuno, Phys. Plasmas 7, 2314 共2000兲.
13
T. Tatsuno, F. Volponi, and Z. Yoshida, Phys. Plasmas 8, 399 共2001兲.
14
A. Ito, Z. Yoshida, T. Tatsuno et al., Phys. Plasmas 9, 4856 共2002兲.
15
A. Miura and P. L. Pritchett, J. Geophys. Res. 87, 7431 共1982兲.
16
Y. Y. Lau and C. S. Liu, Phys. Fluids 23, 939 共1980兲.
17
H. R. Strauss, Phys. Fluids 19, 134 共1976兲.
18
T. Tatsuno and W. Dorland, J. Plasma Fusion Res. 1, 015 共2006兲.
19
J. P. Boyd, Chebyshev and Fourier Spectral Method 共Dover, New York,
2000兲.
20
S. A. Orszag, J. Atmos. Sci. 28, 1074 共1971兲.
21
K. S. Breuer and R. M. Everson, J. Comput. Phys. 99, 56 共1992兲.
22
J. D. Hoffman, Numerical Methods for Engineers and Scientists, 2nd ed.
共Marcel Dekker, New York, 2001兲.
23
E. A. Coutsias and J. P. Lynov, Physica D 51, 482 共1991兲.
24
E. A. Coutsias, T. Hagstrom, J. S. Hesthaven et al., Spec. Issue Houston J.
Math., 1996, p. 21.
25
W. H. Press, S. A. Teukolsky, W. T. Vetterling et al., Numerical Recipes
共Cambridge University Press, Cambridge, 1988兲.
26
Y. Ishii, M. Azumi, and Y. Kishimoto, Phys. Rev. Lett. 89, 205002 共2002兲.

Downloaded 13 Sep 2006 to 128.8.80.199. Redistribution subject to AIP license or copyright, see http://pop.aip.org/pop/copyright.jsp