Astronomical Notes

Astronomische Nachrichten

1!%"*#(6$ 0!(62*4 ,"#! 2* 03! *"*-2*!(6 !)"-+02"* ": ,(4*!0"6"0(02"*(2*50(&2-20$/9 /(05+*" !+3 .9 8"6-(*#$Q+VQ( N*( ?80Q+V0M8 $*,G=V!V0*+ !+3 E!V2Q,!V08!- E*3Q-0+J4 >2Q )+0;Q(&0V7 *N E!(7-!+34 $*--QJQ A!(. S5#KS4 E!(71-!+34 )?'
@Q8Q0;Q3 S55P E!7 HL4 !88QGVQ3 S55P I=+ SP
A=<-0&2Q3 *+-0+Q S55P '=J U5
7!$ '"6#5!88(QV0*+ 30&.& O 0+&V!<0-0V0Q& O ,!J+QV08 MQ-3& O ,!J+QV*273(*37+!,08& FET"D O G-!&,!&

C=,Q(08!- 0+;Q&V0J!V0*+ *N V2Q V:*130,Q+&0*+!- ,!J+QV08 (Q8*++Q8V0*+ 0& J0;Q+ 0+ V2Q 8*+VQ9V *N V2Q +*+-0+Q!( Q;*-=V0*+*N V2Q E!J+QV*1@*V!V0*+!- R+&V!<0-0V7 FE@RDB %0V2 ! 8!(QN=- 8*,G!(0&*+ V* ;!(0*=& V2Q*(0Q& =&0+J <*V2 *+Q1 !+3 V:*1
30,Q+&0*+!- !+!-7&0&4 0V 0& N*=+3 V2!V ! +Q: &V!<0-060+J QNNQ8V *N V2Q 8Q+V(0N=J!- N*(8Q *+ VQ!(0+J 0+&V!<0-0V7 ,=&V <Q G(Q&Q+V
0+ V2Q &GQ80M8 JQ*,QV(7 *N V2Q E@RB E!J+QV08 (Q8*++Q8V0*+ ,0J2V G-!7 ! .Q7 (*-Q V* V2Q N*(,!V0*+ *N V2Q +*+!90&7,,QV(08&V(=8V=(Q& *<&Q(;Q3 0+ E@R Q9GQ(0,Q+V&B >2Q (Q&=-V& ,!7 !-&* <Q =&QN=- N*( V2Q Q&V0,!VQ *N V2Q !88(QV0*+ (!VQ 0+ ;!(0*=&!&V(*G27&08!- *</Q8V&B

Astron. Nachr. / AN 329, No. 7, 688 –700 (2008) / DOI 10.1002/asna.200811013

T
N
I
R
P
E
R

Astron. Nachr. / AN 329, No. 7, 688 – 700 (2008) / DOI 10.1002/asna.200811013

Secondary tearing mode in the nonlinear evolution of magnetorotational
instability
T. Tatsuno and W. Dorland
Center for Scientific Computation and Mathematical Modeling, The University of Maryland, College Park 20742, Maryland, USA
Received 2008 May 15, accepted 2008 Jun 28
Published online 2008 Aug 30
Key words

accretion disks – instabilities – magnetic fields – magnetohydrodynamics (MHD) – plasmas

Numerical investigation of the two-dimensional magnetic reconnection is given in the context of the nonlinear evolutionof the Magneto-Rotational Instability (MRI). With a careful comparison to various theories using both one- and twodimensional analysis, it is found that a new stabilizing effect of the centrifugal force on tearing instability must be present
in the specific geometry of the MRI. Magnetic reconnection might play a key role to the formation of the nonaxisymmetricstructures observed in MRI experiments. The results may also be useful for the estimate of the accretion rate in various
astrophysical objects.
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim

1 Introduction
Magnetic reconnection is one of the most fundamental
phenomena in plasmas. It often entails vastly multiscale
physics, brings about the change of global topology, converts magnetic energy into particle energy such as heat andacceleration and plays a very important role in the turbulent dynamo phenomena. While a fairly deep study and understanding have been made for the idealized steady reconnecting system by taking various physical effects into account, the understanding of reconnection in the complex
development of the nonlinear plasma phenomena has notbeen achieved even with a simplest model. Here we try toinvestigate such a problem, namely the dynamic magnetic
reconnection in the nonlinear evolution of the backgroundplasma using the simplest model, two-dimensional incompressible magnetohydrodynamics (MHD). We explore thevalidity of several idealized theories on a dynamical system,
and report that the problem is not so straightforward as one
might expect. This is the first step to the understanding of
the dynamic magnetic reconnection in the real physical system (developing system rather than an idealized stationary
system).
In this paper we focus on the nonlinear developmentof the Magnetorotational instability (MRI) for our reconnection study. The first calculation of the MRI goes backto Velikhov (1959) and Chandrasekhar (1960), but it is after Balbus & Hawley (1991) that it has attracted attention
for explaining accretion flows in various astrophysical objects. It has not been widely discussed so far, but in orderfor both the accreting matter and the angular-momentumtransporting matter to be detached from the main body of Corresponding author: tatsuno@umd.edu

the magnetized plasma, magnetic reconnection has to takeplace in the nonlinear evolution of the MRI (Hawley & Balbus 1991; Machida & Matsumoto 2003). In this sense magnetic reconnection should be somehow related to, or mighteven govern the actual accretion rate. Dynamic magnetic reconnection is analyzed in other contexts such as supernova
shock (Tanuma et al. 2001), Parker instability (Tanuma etal. 2003) and turbulent MHD plasmas (Lazarian & Vishniac 1999), however, there is no dedicated analysis in the
nonlinear evolution of MRI.
Recently there are several experimental efforts being
made to clarify the physics of the MRI (Sisan et al. 2004; Ji
et al. 2006; Stefani et al. 2006; Stefani et al. 2007; Wang etal. 2007).
In the spherical experiment at Maryland using liquidsodium (Sisan et al. 2004), coupled velocity and magnetic
field fluctuations are observed above critical applied magnetic field. These fluctuations mostly show nonaxisymmetric structures while the primary MRI mode is expected axisymmetric. Nevertheless the threshold of the applied magnetic field seems to agree quite well with the theoretical predictions of the stability boundary (Tillotson 2007). Cylindrical experiments are also being carried out using liquid
metal (Ji et al. 2006; Stefani et al. 2006; Stefani et al. 2007)
and plasma (Wang et al. 2007). PROMISE device (Stefani
et al. 2006) makes an experiment with a helical field using liquid metal alloy. For a fixed azimuthal field strength
they have increased the axial field and observed increased
level of axial velocity fluctuations in agreement with numerical predictions. They also observed nonaxisymmetricfluctuation in some parameter regimes (Stefani et al. 2007).
These nonaxisymmetric structure may be attributed to thesecondary instability of some sort, one of such an example is the fluid dynamical wavy mode (Davey, Diprima &c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, WeinheimAstron. Nachr. / AN (2008)

689

Stuart 1968). However, as the critical Reynolds number ofthe wavy instability seems to increase in the presence of themagnetic field (Tabeling 1981; Willis & Barenghi 2002),
one might speculate that other MHD instabilities, such as
tearing mode, the three-dimensional version of the reconnection presented here, might create the nonaxisymmetricfluctuation. Thus, it is very important to understand the dynamics of the magnetic reconnection in this sense as well.
In this paper we investigate various (linear) properties of
the two-dimensional (2D) magnetic reconnection occurringin the context of the laboratory MRI experiment. We firstshow the normalized 2D equations for MRI in Sect. 2. Using our upgraded reduced MHD code (Tatsuno & Dorland
2006), we show when and how reconnection event takes
place in Sect. 3 including the general nonlinear evolutionof the MRI. Section 4 is the main part of the paper. We
investigate the validity of various linear reconnection theory (tearing mode) using one-dimensional (1D) approximation (Furth, Killeen & Rosenbluth 1963; Coppi et al. 1976;
Paris 1982) and also make the 2D linearized simulation forthe comparison with other 2D results (Nishikawa & Sakai
1982; Paris 1987). We observe a beautiful 2D tearing eigenfunction during reconnection, and conclude that reconnection is attributed to the tearing instability of the opposing
radial magnetic field which is created by the nonlinear evolution of the MRI. Note, however, that the creation of themagnetic null (X-point) is not necessarily attributed to thetearing mode. In some parameter regimes, X-point is created passively by the MRI but reconnection does not continue unless the surrounding geometry is made tearing unstable. For sufficiently high magnetic Reynolds number weobserved a creation of the secondary island, which is briefly
summarized in Sect. 5 with a comparison to the tearing instability of the current sheet (Bulanov, Sakai & Syrovatskii
1979). Finally, our conclusion is given in Sect. 6.

(9)

We consider a concentric cylinder with rigid ideal conducting walls at radii r = rin and rout . We assume periodicity
in the z direction.
2.2.1 Flow
Noslip boundary condition is imposed on the velocity field:vθ (rin ) = vin

and vθ (rout ) = vout .

(10)

Corresponding to
vr = vz = 0

at both walls,(11)

we may impose
φ00 = 0, ∂r φ0 |rout = 0∂r (rω0 ) = 0 at both walls

(12)
(13)

for kz = 0 component, and homogeneous Dirichlet condition for φ and the converted ω condition as is used in theslab case for kz = 0 (Tatsuno & Dorland 2006).
2.2.2 Magnetic field
Ideally conducting wall yields
n · B = 0,

n×j =0at both walls,(14)

in time and along walls,(15)

at both walls,(16)∂r (rBθ ) = 0

We employ the 2D incompressible MHD equations in ther-z plane of the cylindrical configuration (r, θ, z). Letv = ∇ × [φ(r, z)∇θ] + rvθ (r, z)∇θ,(1)

B = ∇ × [ψ(r, z)∇θ] + rBθ (r, z)∇θ,

(2)



(3)
(4)
(5)

then the normalized 2D MHD equations for the incompressible plasmas yield 2 
 2

 ω
v
Bθ
j
∂t ω+ φ, +∂z θ = ψ, +∂z+ν∇2∗ (rω), (6)rrrr
www.an-journal.org(8)

2.2 Boundary conditions

ψ = const

2.1 Formulation

1
1∂r φ + ∂z2 φ =: ∇2∗ φ,
ω = −(∇ × v)θ = ∂r
r
r


1
1
j = −(∇ × B)θ = ∂r∂r ψ + ∂z2 ψ = ∇2∗ ψ,r
r[f, g] = ∂r f ∂z g − ∂z f ∂r f,(7)

which correspond to

2 Reduced MHD equations for MRI

and define

1
1[φ, rvθ ] = 2 [ψ, rBθ ] + ν∇2∗ (rvθ ),2rr
1∂t ψ + [φ, ψ] = ηrj,
r

 
Bθ
vθ 
= η∇2∗ (rBθ ).
∂t Bθ + φ,
− ψ,
r
r

∂t vθ +

where n represents a unit vector normal to the surface. Notethat ∂t ψ = 0 at wall would yield ∇2∗ ψ = 0 from induction
equation, together with vr = Br = 0.
In the vector form and for kz = 0 component, ∂r Bz = 0
is derived from jθ = 0 condition.
2.3 Description of the code
We have extended the slab code used in Tatsuno & Dorland(2006) to solve the cylindrical problem. This is a pseudospectral code with Chebyshev-Fourier expansion and solves
reduced MHD Eqs. (6)–(9). Spectral accuracy is very usefulfor detailed diagnostics presented here. Normally the reconnection point does not lie exactly on top of the grid and one
has to interpolate various physical quantities. Such interpolations are easily made with high accuracy by a spectral
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
690

T. Tatsuno & W. Dorland: Secondary tearing mode in MRI

scheme, which makes possible the sophisticated measurement of some complicated quantities. For the time integration it uses 3rd order Adams-Bashforth scheme for nonlinear terms and Crank-Nicholson scheme for linear terms.

Details of the upgrade we have made for the cylindricalversion is described in Appendix A.

0.6
0.5
0.4
γ

The automatic resolution adjuster is implemented in thecode, which uses an empirical error estimate at each time
step and changes resolution globally by tapping zeroes in
the Chebyshev (or Fourier) space or just removing small
amplitudes. Time step is also automatically adjusted by the
Courant condition.

0.7

0.3
0.2
0.1

ideal shooting
4
Rm=10 simulation

0
0

2

4

rin = 0.5,

rout = 1.5,

ν = η = 1.5 × 10−4 ,

14

φ
ψ

0.025
0.02
0.015
0.01
0.005

0

Bz = 0.15,

S = 1500,

120Lz = 1,vθ (rin ) = 1,

(17)which correspond to
Rm = 104 ,10

Fig. 1 Growth rates taken from ideal shooting code and from
dissipative simulation.

perturbation amplitude

In this section we show the general evolution of the MRI andshow when and where magnetic reconnection takes place.
We also discuss if saturation mechanism is attributed to it(Goodman & Xu 1994). Parameters are:

8
kz

0.03

3 General evolution of MRI

6

Pm = 1,

5

10

15
t

20

25

30

Fig. 2 (online colour at: www.an-journal.org) Time evolution of
MRI perturbation amplitude.
ψ(t=22)

(18)

1
0.8

where these dimensionless numbers are defined byVθ,in rout
,
Rm =
η

Bz routS=,ηνPm = .
η

0.6 z

(19)

0.4
0.2

Remember that we use velocity unit for Bz in our normalized equation.
We used the Keplerian flow profile (vθ ∼ r−1/2 ). Keplerian flow is not a stationary solution of the reduced MHD
equations with finite dissipation rates, however, we do notcare about the decay of the background flow since the timescale of the dissipation for a global structure is much slower
than the dynamics we are interested in. In all of our simulations we use a magnetic Prandtl number of unity.

0
0.6

0.8

1
r

1.2

1.4

Fig. 3 Magnetic island is observed as a consequence of the reconnection. Contour scale is not equally spaced in order to highlight the island.

3.2 Saturation mechanism
3.1 Linear growth
In order to prove the validity of our code, we show in Fig. 1
the growth rates taken from our 2D simulation as well as
the results from an independent ideal (ν = η = 0) shooting code. For this simulation we took Lz = 2 to have four
unstable modes in the domain to compare with the shooting
result. All four eigenvalues are sufficiently close but lower
than the ideal ones. We also confirmed that the ideal run using our code gives exact growth rates on top of the line fromthe shooting code.
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, WeinheimTime evolution of the L2 norm ||φ − φb || and ||ψ − ψb ||are shown in Fig. 2, where the subscript b denotes the background field (axial field and azimuthal rotation). Saturation
is first observed in φ about t  18 and that of ψ follows
about t  21. As a consequence of the magnetic reconnection, magnetic island is observed at t = 22 as seen in Fig. 3.
Notice the island around (r, z)  (0.95, 0.75). It is truethat the coincidence of the peaking time of ψ curve and theappearance of the island implies some relationships, since
stabilizing effect of the outflow stops at this stage, however,they are not simultaneous in all cases. As we show later,
www.an-journal.org

Astron. Nachr. / AN (2008)

ψX

691

ψO

2

t=0
t=15
t=18
t=20

1.8
1.6

z=0.75Ω

1.4

thickness

1.2
1
0.8
0.6
0.4
0.2
0.6

width

rX
Fig. 4

0.8

1

1.2

1.4

r

rO

Fig. 5 (online colour at: www.an-journal.org) Radial profile of
the averaged angular velocity.

Schematic view of the geometry of the magnetic island.

there are cases where island appears long before the saturation of ψ or even before the peak of φ.
Why is only the outgoing field line stretched here whileboth incoming and outgoing fields are observed symmetricin the shearing slab simulation (Hawley & Balbus 1991)?
It is surely related to the shape of the Taylor vortex whichis narrower in the outgoing direction than the other. In fact
the radial velocity at (r, z)  (0.95, 0.75) and t = 18 isabout vr = 0.087, which roughly agrees with the extension
speed of the field line. Keeping in mind that the local MRI
growth rate is proportional to the local rotation rate (Balbus& Hawley 1991), the MRI grows faster in the inside region
due to its fast time scale, and the inward jet is killed by the
inner boundary there.
If this is the case, there may be observed two or more
outgoing jets simultaneously when higher wave number ismore unstable. In fact we observed two (or three) fingers
by extending the axial length twice, but as is observed by
Hawley & Balbus (1992) and Liu, Goodman & Ji (2006),
they eventually merge to form a single finger at a later time.
Such phenomena are of great interest from the viewpoint of
the magnetic reconnection since merging of fingers involve
numbers of reconnection events. However, we concentrate
on a formation of a single finger in this paper for the simplicity of the analysis, since the axial motion of the finger
introduces the complex motion of the X-points. In our configuration with a single finger, the X-point moves only along
the radial direction and not along the z-axis, which enables
us the detailed diagnostics with a relatively simple implementation.
Figure 4 shows the schematic view of the geometry of
the magnetic island. We implemented measurements of theflux functions at the X- (ψX ) and O-points (ψO ) and theisland size as depicted. The inspection of the X- and Opoints are made by detecting the sign change of Bz along
the z = 0.75 line at every 0.01 time unit. From the time
evolution of the magnetic island we confirmed that the island does not move in the z direction. It should be emphasized that we obtain rX and ψX with a spectrally interpolated accuracy, which may not have been achieved by thefinite difference scheme. In the present simulation the iswww.an-journal.org

land is first observed at t = 20.45, grows until t  24 andeventually disappears at t = 32.76. Another island appears
at t = 37.24 and follows more or less the same processes.
Since our interest is on the dynamical evolution of the island, we do not ask the existence of the magnetic island in
the final stationary state. Such a steady state itself may or
may not exist, which may be separately an interesting subject of research and may be analyzed by tracing a steadystate solution for various plasma parameters finding a bifurcation diagram (Rincon, Ogilvie & Proctor 2007).
From the fact that the saturation of ψ follows that of φand that the magnetic reconnection starts at the saturated
level of ψ, we conclude that the magnetic reconnection is
not attributed to the saturation of the MRI as Goodman &Xu (1994) has proposed in the shearing sheet geometry. In
our geometry, the created radial magnetic field strength iscomparable to the original axial field while Goodman & Xu
assumes Br  Bz . Instead the saturation may be attributed
to the disappearance of the drive as seen in the flat profile of
the averaged angular velocity (See Fig. 5).
3.3 Profiles of various fields
Right before the emergence of the magnetic island, at t =20, the profiles of various field quantities sliced along r 
0.8 are shown in Fig. 6. There is a sharp counter-oriented
radial magnetic field (Br ) centered at z = 0.75, which introduces the following magnetic reconnection event. The
axial field (Bz ) is almost zero around z  0.75. The azimuthal field (Bθ ) also contains the counter-oriented component centered at the same place, which forms anothersource of the reconnection for a nonaxisymmetric perturbation. The radial velocity (vr ) shows a weak outflow representing a jet, and the profile of axial velocity (vz ) showsa weak counter-flow which may correspond to the inflow of
the reconnection.
It is noted that Kelvin-Helmholtz instability does not occur for this configuration. This is because the magnitude ofvr and vz are much smaller than the corresponding magnetic field components, and the spatial scale of vθ is muchlarger than that of the other fields.
Now that we have finished explaining the environmentof the reconnection event, so in the following sections we
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
692

T. Tatsuno & W. Dorland: Secondary tearing mode in MRI
500vrvθ-0.5vz

0F’’ / F @ r=0.797379

vr, vθ, vz @ r=0.797379

0.3
0.25
0.2
0.15
0.1
0.05

-500
-1000
-1500
-2000
-2500
-3000

0

-3500

-0.05
0

0.2

0.4

0.6

0.8

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

1

z

z

Br, Bθ, Bz @ r=0.797379

0.3

Potential F /F of the outer region ODE (22).

Fig. 7

Br
Bθ
Bz

0.2



where we defined F = kx Bx (Biskamp 2000, [4.21]). By
numerically solving (22), we estimate the jump of the logarithmic derivative

0.1
0
-0.1

∆ =

-0.2
-0.3
0

0.2

0.4

0.6

0.8

1

bz (+) − bz (−),
bz (zsheet )(23)

Fig. 6 (online colour at: www.an-journal.org) Various field profiles sliced along r  0.8 at t = 20.

where zsheet is the field reversal plane.
The growth rate may be obtained from (4.66) of
Biskamp (2000) (see also Paris 1982):

proceed to actually characterize them in as much detail as
possible.

8 ∆
−
π kx

z



τA
τη

1/3
=

λ5/6 Γ [(λ − 1)/4],
Γ [(λ + 5)/4]

(24)

where

4 First reconnection
4.1 Slab tearing modelThe first thing we can come out of is the calculation of ∆and the consequent growth rate in comparison to the basic
tearing mode analysis (Furth et al. 1963; Coppi et al. 1976;
Paris 1982). It actually turns out from the more detailed
analysis that slab approximation is a too much simplification, but we here describe what we get from slab analysis
for the comparison to the later sections.
4.1.1 Theory(20)

which of course satisfies ∇ · B 0 = 0, and neglect all background velocity. Invoking slab geometry we also assume
that curvature effect is negligible. We want to obtain ∆ ,the jump of the logarithmic derivative of the perturbed flux
function in the outer region. The ordinary differential equation for the outer region where inertial effect may be neglected comes from the static vorticity equation∇ × (B · ∇B) = 0.

(21)

Decomposing B = B 0 + b, linearizing and assuming
exp(ikx x) dependence yields


d2 bz
F 
2
−
k
+
(22)
bz = 0,x
dz 2
Fc 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim


(25)
(26)
(27)

and Γ is the Gamma function. In the case of small ∆ wemay take the small λ limit and recover the FKR result (Furthet al. 1963)  4/5
4/5
2Γ(5/4)∆−2/5 −3/5
γ
τA τη
.
(28)
πΓ(3/4)kx
The large ∆ limit may correspond to the limit λ → 1 − 0
yielding

Here we take the equilibrium
B 0 = [Bx (z), By (z), 0],λ = γ̂ 3/2 = γ 3/2 τη1/2 τA ,1
τη = 2 ,kx ηkxτA =  ,F−2/3

γ  τη−1/3 τA

,

(29)

which we call CGPR limit (Coppi et al. 1976).
4.1.2 Numerical estimate of ∆In order to invoke the discussion of the previous subsection,
the boundary condition (BC) has to be determined. Let us
find out what can be imposed at where. Since Br looks like
hyperbolic tangent profile around z  0.75 and Br changes
sign just below z = 0.7, we may spot the light on this point.
The function F  /F is shown in Fig. 7 at t = 20 andr  0.8. Notice that F  /F doesn’t depend on the exact
value of kx . The minimum of the well is at z = 0.75 whichwww.an-journal.org

Astron. Nachr. / AN (2008)

Exponential BC So we should consider lower energy
and then we may impose the exponential BC at anywhere
we want (we denote it zstart ) as far as it is outside of thewell since the solution out of the well is going to be almost
exponentially decaying. The gradient of bz is given so thatit
matches the exponent corresponding to the local value ofkx2 + F  /F . For example, the exponential BC given atzstart = 0.6 would be sufficient for −kx2  −100 (The second smallest minima are around −50). Figure 8 shows thedependence of ∆ as a function of zstart for kx = 25 withthe exponential BC and the eigenfunction corresponding tozstart = 0.6. For this setup, we conclude that ∆  92.

95
90

∆’

85
80
75
70
65
0

0.1

0.2

www.an-journal.org

0.4 0.5
zstart

0.6

0.7

0.8

20

bz

15
10
5
0
0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 0.74 0.76
zFig. 8 Values of ∆ for various zstart and the eigenfunction corresponding to zstart = 0.6.
1.2

kx=2

1
0.8
0.6
0.4
0.2
0
0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
z

100

1

kx=30

10
1
0.1
bz

Periodic BC Another fact we may use for the BC is theperiodicity of the system along z. Since the local minimumof the potential is about 40 around z  0.25, we simply
conjecture that resistive layer is not formed for sufficiently
negative −kx2 . Then, from the symmetry around the current
sheet at zsheet = 0.75, we may just flip the sign of bz there
and continue integration to the right. Periodicity brings thesolution to z = 0 when it crosses 1, we may integrate normally through z = 0.25 and match the solution with itself
at the starting point zstart = 0.6.
The scheme seems well-defined even for a fairly small
kx , and the resulting solution is shown in Fig. 9 for kx = 2
and kx = 30. The logarithmic plot is shown for kx = 30,
which represents a very good exponential behavior between
the two field reversal planes. We also show the ∆ obtained
from this scheme with a comparison to the previous exponential assumption. The scheme doesn’t work with too highkx since the exponential behavior becomes very rapid andwe cannot match it even with a double precision computation. On the other hand, it does work for even smaller values of kx than the exponential BC. So these two schemes
are complementary, and the difference is negligible in the

0.3

25

bz

corresponds to the field reversal plane. There is a very mild
concavity extending 0.1  z  0.4, but this should notgive a major effect and hence we neglect it in the following
global calculation.
Since (22) takes exactly the form of the Schrödinger
equation, we may regard F  /F as a potential and −kx2 as
an energy. If the energy −kx2 is larger than the minimumof the potential well F  /F , there is a chance that we may
have oscillatory solution, while smaller, we may not. Thisoscillatory behavior is needed for positive ∆ since otherwise two exponentially growing solutions from outside of
the well would never match at the center of the well withpositive ∆ .
Since the potential is bounded there is clearly a continuum above its maximum if we allow −kx2 to be positive.
This fact brings a difficulty in the numerical estimate of the∆ for small wave number. For example, the energy corresponding to kx = 0.5 is too large (−kx2 = −0.25) and thesolution is strongly affected by the existence of the continuum.

693

0.01
0.001
1e-04
1e-05
0

Fig. 9dition.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9z1

Eigenfunction obtained with the periodic boundary con-overlapping region as shown in Fig. 10. The critical wave
number is obtained from the exponential BC and is aboutkx  40.8.
4.1.3 Results
Now the question is how we use this result. We should ofcourse have that 2π/kmax < Lr where kmax and Lr denote the critical wave number for positive ∆ and the radial
length of the stretched field line. But probably more strinc 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim694

T. Tatsuno & W. Dorland: Secondary tearing mode in MRI
-4

10000

10

periodic BC
exponential BC

fullψ onlyno vθ
vθ,avevθ,rigid

10-6
-8

10

1000

-10

10
||φ1||

∆’

-12

100

10

10-14
10-16
10-18

10

-20

10

10-22

1
0

5

10

15

20

25

30

35

40

10-24

45

kx

γ

shooting
FKR
CGPR

40

50

60

70

80

90

100

Fig. 12 (online colour at: www.an-journal.org) Time evolution
of ||φ1 || with various background setup. The lines labeled ‘full’used full set of the background fields; ‘ψ-only’ retained only ψ0and set φ0 , vθ,0 , Bθ,0 = 0; ‘no vθ ’ retained ψ0 , φ0 and Bθ,0 withvθ,0 = 0; ‘vθ,ave ’ retained ψ0 , φ0 and Bθ,0 with only vθ,0 averaged along z; and ‘vθ,rigid ’ retained ψ0 , φ0 and Bθ,0 with onlyvθ,0 replaced by a rigid rotation corresponding to the angular velocity Ω = 0.7.

4.2 Linearized cylindrical simulation with a fixedbackground
0

5

10

15

20

25

30

35

40

45

kx

Fig. 11 (online colour at: www.an-journal.org) Dispersion relation of the slab tearing instability for the background profile ofFig. 6.

gent condition would be that the highest growth rate must belarger than vr /Lr since otherwise the mode may be transferred away before it grows. We may estimate Lr by picking
a specific value of the flux ψ and see how long it is elongatedradially. Thus we need to obtain the growth rate.
The growth rate is shown with respect to kx in Fig. 11
with a comparison to FKR and CGPR’s analytic results. Our
scheme fails for kx  1.5 since bz (zsheet ) becomes too
close to zero and it erroneously changes sign. Nevertheless,the growth rate approaches the ∆ → ∞ limit, so that we
may use analytic estimate for this region. In any case thisdoes not cause any problem since the largest growth rate iscaptured by our numerical estimate.
The radial length of the stretched field line Lr may beevaluated by the radius at which the background ψ equals aparticular value along z = 0.75 and z = 0.25 chosen well inthe stretched field area. Taking ψ = 0.025 gives us Lr (ψ =0.025)  0.53 and, as seen in Fig. 6, vr (z = 0.75)  0.062.
Thus
γ

30

t

Fig. 10 (online colour at: www.an-journal.org) ∆ obtained withtwo boundary conditions.
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

20

Lr
 7.78  1,
vr

(30)

which shows that the tearing mode grows fast enough before it is transferred away by the radial flow. It is noted that
we don’t have to worry about vθ since the evolution is axisymmetric as far as this analysis is concerned.
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
Here we fix all field quantities at a given time in the full
nonlinear simulation and evolve the linearized equation tosee if there is any instability. We divide the quantities as
f = f0 + f1 , take the field of the nonlinear simulation atsome instant as f0 , and fix it. Then we solve the linearized
reduced MHD equations (same 2D cylindrical ones):ω1  
ω0  2∂t ω1 + φ0 ,+ φ1 ,+ ∂z (vθ,0 vθ,1 ) =
r
r

 r

j1
j0
2
ψ0 ,
+ ψ1 ,+ ∂z (Bθ,0 Bθ,1 )+ν∇2∗ (rω1 ), (31)rr
r
1∂t vθ,1 + 2 ([φ0 , rvθ,1 ] + [φ1 , rvθ,0 ])r
1= 2 ([ψ0 , rBθ,1 ] + [ψ1 , rBθ,0 ]) + ν∇2∗ (rvθ,1 ), (32)r1
∂t ψ1 + ([φ0 , ψ1 ] + [φ1 , ψ0 ]) = ηrj1 ,
(33)
r
 

Bθ,1
Bθ,0
∂t Bθ,1 + φ0 ,+ φ1 ,rr
vθ,1  vθ,0 
+ ψ1 ,+ η∇2∗ (rBθ,1 ),(34)= ψ0 ,
rr
with a tiny initial perturbation on ω1 . The background field
quantities (denoted with subscript zero) do not satisfy theequilibrium equations exactly, however, we may neglect theeffect of their little imbalance as far as their changes are
slow compared to the reconnection dynamics. As a validation test of the linearized version of the code, we have reproduced the growth rates and eigenfunctions of the primary
MRI and tearing instability with tanh(r)-like field profile.
The time evolution of the perturbed quantities are shownin Fig. 12 with the fixed background corresponding to t =20 of the nonlinear simulation. Here we note that t = 20
corresponds to the time right before the appearance of the
X-point (see Sect. 3). The lines labeled ‘full’ used full set
www.an-journal.org

Astron. Nachr. / AN (2008)695

of the background fields; ‘ψ-only’ retained only ψ0 and set
φ0 , vθ,0 , Bθ,0 = 0; ‘no vθ ’ retained ψ0 , φ0 and Bθ,0 withvθ,0 = 0; ‘vθ,ave ’ retained ψ0 , φ0 and Bθ,0 with only vθ,0averaged along z; and ‘vθ,rigid ’ retained ψ0 , φ0 and Bθ,0with only vθ,0 replaced by a rigid rotation.
4.2.1 Stabilizing effect of rigid rotation
As ψ-only line shows, the flux function ψ0 is tearing unstable with a comparable growth rate to the time scale of thereconnection (γ  0.58). The result with full backgroundshows a very slow growth, while killing only vθ immediately recovers the comparable growth rate to the ψ-onlycase. This fact suggests that vθ,0 has a strong stabilizing influence and the slab analysis we made in the previous section may be a too much simplification.
In order to confirm the stabilizing effect of vθ,0 , we haveturned on only the kz = 0 component of the background vθ
(labeled vθ,ave ). The difference between no-vθ and vθ,ave isonly the average (kz = 0) component of the backgroundvθ . From the significant decrease of the growth rate with theaverage vθ component, we may conclude that the stabilizing
effect is not given by the z-shear of vθ , which in fact has amuch larger length scale than the field reversal.
By the inspection of the radial profile of the angular velocity Ω (see Fig. 5), the radial Ω shear is highly localizedaround the inner wall for the averaged vθ model we usedin the linearized simulation. Thus we may conjecture thatthe stabilization may not be attributed to the radial Ω shear
either. Instead the stabilizing effect may be caused by thecentrifugal force.
In order to check our conjecture, we prepared a rigid vθ
profilevθ (r) = 0.7r,(35)which gives rise to the constant angular velocity Ω = 0.7,
taken from the averaged angular velocity around the reconnection site. The result is shown in Fig. 12 with a label vθ,rigid . From the fact that there is no growth at all forthis case, we may draw the following conclusion. The tearing instability is completely stabilized due to the centrifugal force of the background rotation vθ,rigid , and the slow
growth shown in the full and vθ,ave simulations is a remnantof the magnetorotational instability, which is also indicatedby the structure of the eigenfunction (not shown here).
It is interesting to find that centrifugal force seems tohave a stabilizing effect while the similar conclusion isdrawn for the gravitational field in the normal field case
(Nishikawa & Sakai 1982; Paris 1987). The question is if
the reconnection is completely forced. To answer this, wehave to confirm the stabilization of tearing instability at every instant, which may be carried out by changing the time
where background field is fixed.
Table 1 shows a summary of the linearized runs using
the various choice of the retained background fields at several t. Finite growth rates labeled by ∗ have the eigenfuncwww.an-journal.org

Table 1 Growth rates for various combination of the backgroundfield. In the field column, ‘◦’ denotes corresponding field is takenfrom the nonlinear simulation, while ‘–’ denotes it is zeroed out.
The labels ‘ave’ and ‘rig’ for vθ denote the usage of average vθ andrigid rotation. In the growth rate columns, ‘–’ denotes no instability
observed, and the numbers carrying ‘*’ have MRI eigenfunctions.
Case

ψ

Bθ

φ

vθ

γ(t = 21)γ(t = 20)γ(t = 18)

1

◦

–

–

–

0.47

0.58

0.64

2a
2b
2c

◦
◦
◦

–
–
–

–
–
–

◦
ave
rig

0.31

0.32
0.40
0.40

0.35

3

◦

◦

–

–

0.47

0.56

0.52

4a
4b

◦
◦

◦
◦

–
–

◦
rig

0.23

–
–

0.33∗
–

5

◦

◦

◦

–

0.37

0.50

0.13

6a
6b
6c

◦
◦
◦

◦
◦
◦

◦
◦
◦

◦
ave
rig

0.2

0.06∗
0.08∗
–

0.39∗

0.17

–

tion structure corresponding to MRI. As seen in Table 1 wehave tearing instability for ψ-only case, but not with rigid
rotation until t = 20 (cases 1 and 6c). It is noted that for
t = 16 we did not observe tearing instability even with ψonly case. Comparison of cases 1, 2c and 4b for t = 20 and
t = 18 tells us that the application of only rigid rotation
do not completely stabilize the tearing instability, but thatcombined with Bθ does.
4.2.2 Destabilization mechanism
At t = 21, however, we observed a tearing instability witha growth rate γ  0.2 for the full background field. The average angular velocity around reconnection site decreases
only about 10%, but Bθ is only 50% of the t = 20 background. Since the effect of Bθ for finite vθ is big (compare
cases 2a/2c with 4a/4b), we may speculate that reductionof Bθ is affecting the destabilization. In fact by making thelinearized simulation with Bθ (t = 20) and other fields at
t = 21, the growth rate diminishes to γ  0.03. Thus the
background field at t = 20 is almost marginal for tearing
instability.
Figure 13 shows the eigenfunction of the 2D linearized
simulation using the full background field of the nonlinear simulation at t = 21 sliced along r  0.8 and that of
the 1D shooting calculation explained in the previous section using the 1D version of the same background. For the
1D shooting calculation we have taken kx = 20 for which
the growth rate achieves the maximum. The 2D eigenfunction structure suggests that kx  15 along z = 0.75. It isshown that the eigenfunction structure becomes wider and∆ seems smaller for the 2D calculation. The significant
global change of the eigenfunction suggests that the effectof stabilizing mechanism due to centrifugal force and Bθ isglobal and we may need to reconstruct the theory including
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
696

T. Tatsuno & W. Dorland: Secondary tearing mode in MRI
1.2

0.01

shooting: kx=20
2D: ψ-only
2D: full1

Rm=50,000
Rm=40,000
Rm=30,000
Rm=20,000
Rm=10,000
Rm= 5,000
1.5e-3 * t

0.008

ψX - ψO

ψ1

0.8
0.6
0.4

0.006
0.004
0.002

0.2
0
0

0.2

0.4

0.6

0.8

0

1

16

18

20

22

z

24

26

28

30

t

Fig. 13 (online colour at: www.an-journal.org) Eigenfunction
taken from slab shooting analysis and 2D linearized simulation.
Table 2 Critical magnetic Reynolds number for the appearanceof the magnetic island in the nonlinear simulation and the corresponding tearing stability properties with 1D theory and 2D linearized simulation.
Rm

1D tearing

γ of lin sim

2D NL sim

10 000
5 000
4 000
3 000
2 000
1 000
750

unstable
unstable
unstable
unstable
unstable
unstablestable

0.17
0.11
0.08
0.045
—
—
—

islandislandislandno island
no island
no island
no island

the change of the outer solution. From the inspection of the2D linear simulation result, the ∆ is estimated to be∆  12.0,(36)

where we have used the largest gradient and the local minimum of ψ1 around z  0.75 for its jump.
4.2.3 Critical magnetic Reynolds number
The effect of producing magnetic finger becomes weak as
we decrease the Reynolds number due to the increase of themagnetic diffusivity. Thus the magnetic island ceases to appear for a sufficiently small magnetic Reynolds number. In
such a case, we can compare the critical magnetic Reynolds
number of the 2D nonlinear simulation, below which magnetic reconnection does not take place, with the linearized
tearing instability calculation both for 1D theory and 2D
simulation, and make it a validity test of our explanation.
The rough comparison is shown in Table 2. Here for thecritical point of the 1D theory we used the coincidence ofthe quantity in the lhs of (30) to be unity. This condition
may be too generous for the instability but it qualitatively
shows an agreement with the appearance of the magnetic
island in the nonlinear simulation. On the other hand the 2D
linearized simulations show a much better agreement.
4.3 Reconnection rate
Figure 14 shows a plot of the reconnected flux ψX − ψOversus time for various magnetic Reynolds numbers. Herec 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
Fig. 14 (online colour at: www.an-journal.org) Time evolutionof the reconnected flux ψX − ψO at various Rm .

we see two-step behaviors for Rm  20 000. This is because the background field is not unstable to tearing mode
until t  20. The X-point is formed by the nonlinear effectof the overall evolution, however, the system is not tearing
unstable yet. As soon as the background field enters tearingunstable regime, the reconnection speeds-up, leading to thetwo-step behavior in the evolution of the island thickness
and reconnected flux.
Figure 14 also shows that reconnection rate saturates
around Rm  20, 000 and becomes independent of resistivity. From the linear fit, the reconnection rate is evaluatedas
1 d(ψX − ψO )
 0.019,
2
vAu
dt

(37)

where vAu denotes the upstream Alfvén velocity for which
we used a typical value of the maximum Br along r  0.8(Br,max  0.28).
There are other measures for the reconnection rate. The
quantity vin /vAu is one of them, which is
max

vin

z∈(0.7,0.8) vAu

 0.026(38)

at r  0.82 and t = 21 for the Rm = 10 000 run. The
order-of-magnitude of this measurement is consistent withd(ψX −ψO )/dt. This rate corresponds to slow reconnection
as we expect from the usage of the resistive MHD modelwith a homogeneous resistivity.
Another measure may be given by vin /vout , which ismax

z∈(0.7,0.8)

vz (r  0.82)
 0.22,(39)

max |vr (z = 0.75)|

r∈(0.6,1)

for the same situation. Note that vr (z = 0.75) is asymmetric with respect to the X-point and we may obtain the larger
number when we use the maximum value within the island( 0.38). It is interesting to see vout is so different fromvAu , which normally have the same value in Sweet-Parker
reconnection (Sweet 1958; Parker 1963). This is because of
the finite size of the island. This fact is partly because thewww.an-journal.org

Astron. Nachr. / AN (2008)697

ψ(t=36)

-15

10

0.8

||φ1||||ψ1||
exp(0.305 t)

-16

0.78
0.76

z

0.74
0.72

perturbation amplitude10

We find that the secondary island appears above Rm =20 000 as shown in the contour plot of the flux function(Fig. 15). In these parameters, the aspect ratio of the current sheet exceeds 10 around t  35, which coincides withthe simple theoretical argument of the threshold of tearing
instability of the current sheet (A > 2π) (Biskamp 2000).
Notice the difference of critical aspect ratio with the slab
geometry (critical A seems much larger there, see Loureiro
et al. 2005).

-20

10

-21

10

(a)
30

40

50

60

70

80

90

70

80

90

100

t
10-15

||φ1||||ψ1||
exp(0.273 t)-16

perturbation amplitude

10

10-17
-18

10

10-19
10-20
-21

10

10-22
10-23

(b)
30

40

50

60

100

t

Fig. 16 (online colour at: www.an-journal.org) Growth of theperturbed field in the linearized 2D cylindrical simulation: (a)shows the results from rigid simulation with Ω = 0.7 and (b)corresponds to the full background field at t = 35.

φ1 & ψ1 @ r=0.679652

5 Secondary reconnection

10-19

10-23

0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05
rradial simulation box size is finite, but more crucial reason
is the existence of the magnetic field at the outer region as
well. In order to have a radially elongated island, the island
has to push the outer magnetic field, which is resisted by
the tension force of the field line. By only making a larger
domain simulation doesn’t help for enlarging vout .

-18

10

10-22

0.7

Fig. 15 Snapshot of the flux function at t = 36 for Rm =30 000. Secondary magnetic island (inner one) appears for Rm ≥2 · 104 .

-17

10

1.4e-19
1.2e-19
1e-19
8e-20
6e-20
4e-20
2e-20
0
-2e-20-4e-20-6e-20-8e-20

φ1
ψ1

0

0.2

0.4

0.6

0.8

1

z

5.1 Linear analysis
In order to understand the secondary reconnection event
we have made the linearized cylindrical simulation as in
Sect. 4.2. Figure 16 shows the time evolution obtained fromthe linearized simulations for Rm = 20 000 with rigid rotation and full background field.
The time evolution of the perturbation amplitude showsa small difference in oscillation period but the growth rate isalmost same. The place where the secondary reconnection
takes place is around r  0.7 and z  0.75, so we show
the slice along r  0.7 of the eigenfunction for the fullbackground field case in Fig. 17. The flux function suggests
a large ∆ and the instability is identified a strong tearing
instability. As we’ve shown before, the aspect ratio of thecurrent sheet for this case is A  11.7.
5.2 Bulanov’s theory
Bulanov et al. (1979) points out that the expanding outflow
along the current sheet has a stabilizing influence on the
tearing instability, which is attributed to the much greaterwww.an-journal.org

Fig. 17(online colour at: www.an-journal.org) Slice of theeigenfunction for the secondary reconnection event.

critical aspect ratio in the slab geometry. The radial velocity
profile in this case turns out
vr (r  rX , z = 0.75)  1.5(r − rX ) + vr,X

(40)

at t = 35 and anti-symmetric around the X-point. According to the 2D cylindrical simulation of the linearized equations (same calculation with Sect. 4.2) without backgroundφ, we obtained a growth rate of γ  0.82, which is morethan double of that with the finite φ. Thus Bulanov’s stabilizing effect is in fact effective in this case. It is also
pointed out that Table 1 already shows a slight stabilizing
tendency of φ on the tearing instability (Compare cases 3/4b
and 5/6c). However, the vr shear at this stage is very weak
(dvr /dr  0.13–0.45) compared to the present case.
Bulanov’s analysis seems plausible qualitatively, however, the dvr /dr observed in (40) is almost twice as large as
the growth rate without it. The stability condition dvr /dr 
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
698

γ0 is satisfied but the configuration seems still unstable,
where γ0 is the growth rate without vr . Thus the analysis
given by Bulanov et al. is quantitatively unsatisfactory in
our example.

6 Conclusion
We have numerically investigated the 2D magnetic reconnection during the nonlinear evolution of the MRI in thecylindrical geometry with rigid conducting walls. The outgoing radial jet is created in the inner region, which drags
the initial axial field and forms a reversing radial field line
configuration (magnetic finger). From the numerical simulations using our reduced MHD code, magnetic reconnection event is observed in the finger region (Sect. 3). However, unlike the theory by Goodman & Xu (1994) for the
shearing slab geometry, the saturation is attributed to the radial flattening of the average angular velocity profile whoseradial decrease has fed energy to the MRI (Sect. 3.2).
We have focused on the finger region and applied linear
1D and 2D theories (Furth et al. 1963; Coppi et al. 1976;
Paris 1982; Nishikawa & Sakai 1982; Paris 1987). Slab 1D
theory suggests that tearing instability can arise which maylead to the reconnection (Sect. 4.1), however, we need to invoke 2D analysis in order to explain the observation quantitatively (Sect. 4.2).
By the careful comparison to the conventional theories,
we found that a new stabilizing effect must be present for
explaining our results. The saturated MRI includes a radially flattened average angular velocity profile for a fairly
wide regions. It corresponds to a rigid rotation whose centrifugal force has a stabilizing influence on the tearing instability. On the other hand, Kelvin-Helmholtz instability does
not take place since the radial field strength is much stronger
and is stabilized, although the radial outflow jet contains inflection points.
Therefore, the configuration is stable for a much higher
field elongation even with a passively created X-point, andfinally becomes tearing unstable and it is then when the reconnection speeds up (Sect. 4.3). It is interesting to point out
that the creation of an X-point is not necessarily attributed tothe tearing instability, but the reconnection is. The obtained
reconnection rate has no dependence with respect to theresistivity for sufficiently high magnetic Reynolds number
and corresponds to the slow reconnection rate as one might
expect from the usage of the conventional MHD model.
In some parameter regimes we observe a creation of thesecondary island which results from the tearing instabilityof the current sheet (Sect. 5). The threshold of the aspectratio of the current sheet seems much smaller than the slab
cases (Laureiro et al. 2005), and it is closer to the conventional argument of the tearing instability without any flow
effects taken into account. However, the stabilizing effect of
the flow is indeed observed, which suggests that the work by
Bulanov et al. (1979) is qualitatively correct but not quantitatively. The stabilizing effect of the outflow is rather weak.
c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
T. Tatsuno & W. Dorland: Secondary tearing mode in MRIThere are many things to be done to understand the dynamic properties of magnetic reconnection events. First, weneed a theory about the centrifugal stabilization of the tearing instability. The effect of centrifugal (or gravitational)
force perpendicular to the field reversal plane is analyzed
in the early seminal paper (Furth et al. 1963), but that ofthe parallel one is unknown. According to the inspection of
the eigenmode obtained numerically (Fig. 13), this effect is
global and we need to modify the treatment of the outer region as well. To make such a theory may be rather involved,
however, it should be important with various applications to
the most astrophysical objects which include rotation andreconnection, such as accretion disks and earth’s magnetotails. The secondary stability analysis from an exact stationary state may also be a good work in order to support our
results here. In the real turbulent cases including fully turbulent MRI, we need to take the statistics of the reconnection
events analysed here. From the viewpoint of the numerical simulation, it is very interesting to replace the physicalmodel at the reconnection site with kinetic one and make amultiscale computation with a proper boundary condition,
especially for the astrophysical applications.
For more close comparison with the observations of thelaboratory experiments, it is highly desired to make a full
three-dimensional (3D) analysis. In this case, things become
more complicated. It may be especially interesting to ask
whether tearing instability discussed here or 3D fluid instability (Davey et al. 1968) enters in the first place on thebifurcation diagram. In the non-magnetized fluid there is atransition to 3D mode as analyzed by Davey et al. (1968),
however, such a disturbance seems easily stabilized by a superposition of an axial magnetic field (Tabeling 1981; Willis& Barenghi 2002). On the other hand, the magnetic reconnection may be easier for the 3D case than 2D since the MRI
creates azimuthal magnetic field which contains reversing
structure around the radial field reversal plane (See Fig. 6).
Exploration of such a problem may be extremely helpful
to understand the formation of nonaxisymmetric structures
observed in the experiments (Sisan et al. 2004; Stefani et al.
2007).
Acknowledgements. Authors would like to thank Jim Drake,
Richard Paris and François Rincon for valuable discussions.

References
Anderson, E., Bai, Z., Bischof, C., et al.: 1999, LAPACK Users’
Guide, Soc. Ind. Appl. Math., Philadelphia
Balbus, S.A., Hawley, J.F.: 1991, ApJ 376, 214
Biskamp, D.: 1993, Nonlinear Magnetohydrodynamics, Cambridge Univ., Cambridge
Biskamp, D.: 2000, Magnetic Reconnection in Plasmas, Cambridge Univ., Cambridge
Bulanov, S.V., Sakai, J., Syrovatskii, S.I.: 1979, SvJPP 5, 157
Chandrasekhar, S.: 1960, Proc. Natl. Acad. Sci. USA 46, 253
Coppi, B., Galvao, R., Pellat, R., Rosenbluth, M.N.: 1976, SvJPP
2, 533Coutsias, E.A., Hagstrom, T., Hesthaven, J.S., Torres, T.J.: 1996,
Spec. Issue Houston J. Math., p. 21

www.an-journal.org

Astron. Nachr. / AN (2008)

699

Davey, A., Diprima, R.C., Stuart, J.T.: 1968, JFM 31, 17
Furth, H.P., Killeen, J., Rosenbluth, M.N.: 1963, PhFl 6, 459
Goodman, J., Xu, G.: 1994, ApJ 432, 213
Hawley, J.F., Balbus, S.A.: 1991, ApJ 376, 223
Hawley, J.F., Balbus, S.A.: 1992, ApJ 400, 595
Ji, H., Burin, M., Schartman, E., Goodman, J.: 2006, Nature 444,
343
Lazarian, A., Vishniac, E.T.: 1999, ApJ 517, 700
Liu, W., Goodman, J., Ji, H.: 2006, ApJ 643, 306
Loureiro, N.F., Cowley, S.C., Dorland, W., Haines, M.G.,
Schekochihin, A.A.: 2005, Phys Rev Lett 95, 235003
Machida, M., Matsumoto, R.: 2003, ApJ 585, 429
Nishikawa, K.-I., Sakai, J.: 1982, PhFl 25, 1384
Paris, R.B.: 1982, PlPh 24, 1541
Paris, R.B.: 1987, PhFl 30, 102
Parker, E.N.: 1963, ApJS 8, 177
Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T.:
1988, Numerical Recipes, Cambridge Univ., Cambridge
Rincon, F., Ogilvie, G.I., Proctor, M.R.E.: 2007, Phys Rev Lett 98,
254502
Sisan, D.R., Mujica, N., Tillotson, W.A., Huang, Y.-M., Dorland,W., Hassam, A.B., Antonsen, T.M., Lathrop, D.P.: 2004, Phys
Rev Lett 93, 114502
Stefani, F., Gundrum, T., Gerbeth, G., Rüdiger, G., Schultz, M.,
Szklarski, J., Hollerbach, R.: 2006, Phys Rev Lett 97, 184502
Stefani, F., Gundrum, T., Gerbeth, G., Rüdiger, G., Szklarski, J.,
Hollerbach, R.: 2007, NJPh 9, 295
Sweet, P.A.: 1958, NCimS 8, 188
Tabeling, P.: 1981, JFM 112, 329
Tanuma, S., Yokoyama, T., Kudoh, T., Shibata, K.: 2001, ApJ 551,
312
Tanuma, S., Yokoyama, T., Kudoh, T., Shibata, K.: 2003, ApJ 582,
215
Tatsuno, T., Dorland, W.: 2006, PhPl 13, 092107
Tillotson, W.A.: 2007, PhD Thesis
Velikhov, E.P.: 1959, SJETP 36, 995
Wang, Z., Si, J., Liu, W., Li, H.: 2007, preprint
Willis, A.P., Barenghi, C.F.: 2002, JFM 472, 399

A Numerical implementation
For our simulation we used a code developed by Tatsuno & Dorland (2006), however, we needed some upgrade in order to solve
the cylindrical problem. In this section we describe the specificchanges we have made for the cylindrical version. For notation refer Tatsuno & Dorland (2006) or Coutsias et al. (1996).

A.1

Making the band matrix for the inversion becomes a little bit tricky
because of the Jacobian and we end up with a denser matrix. Poisson equation

is written as«„
1
1∂r φ + ∂z2 φ = ∂r2 (rΦ) + ∂r Φ + ∂z2 (rΦ) = ω,
∂r
r
r

(A1)

(A2)

where Φ = φ/r 2 . By multiplying [2] D−2 on both sides from theleft, ∂r and ∂r2 are banded, and we obtain a band matrix on the left
hand side:h´i` 2[2] −2
D(A3)∂r r + ∂r − kz2 r Φ = [2] D−2 ω.

www.an-journal.org

A.2 Time integration
The equations for ω, vθ , and Bθ share their form∂t f = NL + µ∇2∗ (rf ),(A4)

where f is one of the ω, vθ , or Bθ , NL is a symbolic expression for the nonlinear term, and µ is the diffusion coefficient, ν or
η. The diffusion term is treated implicitly, using Crank-Nicholson
scheme. The time derivative and Laplacian terms are discretized
as
«
„ n+1
+ fn
f
f n+1 − f n= NL + µ∇2∗ r
,
(A5)
∆t
2
where the superscript denotes the time step. By manipulation andmultiplication of [2] D−2 , we obtain a banded matrix on the left
hand side:– n+1»∆t ` 2[2] −22 ´ f
D∂ r r + ∂ r − kz r
r−µ
2
r
–
»
∆t
(A6)
= [2] D−2 f n + ∆tNL + µ ∇2∗ (rf n ) .
2
The boundary conditions on f are converted to the ones on f /r,and implemented in the first two rows of (A6).
The equation for ψ has a different form:∂t ψ = NL + ηr∇2∗ ψ(A7)thus its discretization yields
ψ n+1 − ψ n= NL + ηr∇2∗∆t„

ψ n+1 + ψ n
2

«
.

(A8)

Manipulation and multiplication of [2] D−2 on both sides yields»–´∆t ` 2 2[2] −2D∂r r − ∂r r − 1 − kz2 r 2 Ψ n+1
r2 − η
2
«
„
∆t
rj ,
(A9)
= [2] D−2 ψ n + ∆tNL + η
2where we used Ψ = ψ/r 2 andr∇2∗ ψ = ∂r2 (r 2 Ψ ) − ∂r (rΨ ) − Ψ + r 2 ∂z2 Ψ.

(A10)

A.3 Two-step Sherman-Morrison formula

Preconditioning∇2∗ φ = ω

By inverting it and multiplying r 2 , we obtain φ. We used LAPACK
for the inversion of these banded matrix (Anderson et al. 1999).

Tatsuno & Dorland (2006) succeeded to make a band matrix withsingle extra dense row by splitting the matrix equation into odd
and even parts of the Chebyshev coefficients. However, since they
are coupled in the cylindrical geometry, we have two dense rows
related to the boundary conditions. We solve this problem by applying Sherman-Morrison formula twice (Press et al. 1988).
Let A be the perturbed matrix
A = A0 + u1 ⊗ v 1 +u2 ⊗ v 2 ,
|
{z
}

(A11)

A1

where A0 is the band matrix and u1,2 v 1,2 are the vectors. We
want to solve
Ax = b.

(A12)

c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim
700

T. Tatsuno & W. Dorland: Secondary tearing mode in MRI
where

Then, from Sherman-Morrison formula, the solution of(A1 + u2 ⊗ v 2 )x = b(A13)

is obtained by solving
A1 y = b

A1 z 2 = u2

(A14)

as
x=y−

v2 · yz2.
1 + v2 · z2

(A15)

The solution of
A1 y = b(A16)

is obtained by solving

2

(B2)

(B3)

Here 2 F1 is a hypergeometric function, λG the effect of the gravitational force, σ the effect of the normal field, andΛ = [1 + λG (2σ − λG )]1/2 ,(B4)

c = i(λG − σ),λ + iλG
κ=
.
4Λ

(B5)
(B6)

A0 z 1 = u1

(A17)

We only give here a minimal set of description for our purpose.
For the full description of the variables refer to the original paper(Paris 1987).

v 1 · x1
z1.
1 + v1 · z1

(A18)

B.2

A0 x1 = basy = x1 −1 + 2σλG
0
13
2σλλ
1
−
c/Λ
G
@1 −
Λ3/2 41 −F̄ (1)A5
Λ(1 − 4κ)
3 + 4κ
«„17
1 − c/Λ
, n; + κ;.
F̄ (n) = 2 F1
2
4
2
H=

Let us think about the artificial cases with

The solution of
A1 z 2 = u2

Implication of Paris’ dispersion relation(A19)

1. σ = 0 and λG = 0,
2. λG = 0 and σ = 0.

(A20)

It is true that these assumptions break the equilibrium condition,
but they are still helpful in understanding the effect of each term.
In case 1, we obtain

is obtained by solving
A0 z 2 = u2

A0 z 2 = u1asz 2 = z 2 −v 1 · z 2z 2
1 + v 1 · z 2

Λ = 1,

B Analytic estimate of 2D tearing
In this section we discuss the stabilizing effect of the gravitational
force in the two-dimensional magnetic field configuration in orderto point out the similarity to the centrifugal stabilization. This effect is pointed out by the numerical calculation in Nishikawa &Sakai (1982) and implicitly obtained by the analytic calculation of
Paris (1987). In this section we briefly revisit the result of Paris anddiscuss the similarity to the centrifugal force. For the full detail ofthe analysis, refer to the original paper (Paris 1987).

Dispersion relation

When we add a normal in-plane field such as the case of earth’s
magnetotail, the equilibrium has to be sustained by the gravitational force in order to balance the tension force of the magnetic
field. However, when the normal field is weak, to the leading order, the outer region equation can be described by the same second
order ODE without the normal field (Furth et al. 1963).
The effects of the weak normal field and the gravitational force
enter in the resistive layer equation and the proper analysis leads
to the following dispersion relation:∆ =27/2 πΩ Γ(3/4 + κ)H,
1 − (4κ)2 Γ(1/4 + κ)

c 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim


κ=

(A21)

where z 2 turns out to be z 2 = z 1 .

B.1

c = −iσ,1
λ,4

H = 1,(B7)

which exactly recovers the dispersion relation after Eq.(14) ofParis. The effect of the normal field has completely diminishedby only making λG = 0.
On the other hand, case 2 leaves out a finite effect of gravitysinceqc = iλG ,Λ = 1 − λ2G < 1,
λ + iλG,
κ= p
4 1 − λ2GH = Λ−3/2 > 1.

(B8)

Defining the effective ∆ by∆eﬀ =∆orig< ∆orig ,H(B9)where ∆orig is the value matched with outer solutions, we see thatfor a finite λG , (i) the effective ∆ becomes smaller, and (ii) λbecomes smaller for a given κ (due to the smaller denominator),both of which suggest stabilizing influence. It is noted, however,that the second statement may not exactly apply because the imaginary number iλG is also added in the numerator.

(B1)

www.an-journal.org