Investigation of Kinetic Physics of Magnetic Reconnection under a Finite
Guide Field
by
Xiang Lu
A dissertation submitted to the Graduate Faculty of
Auburn University
in partial ful?llment of the
requirements for the Degree of
Doctor of Philosophy
Auburn, Alabama
December 12, 2011
Keywords: Plasmas, Magnetic reconnection, Particle simulation
Copyright 2011 by Xiang Lu
Approved by
Yu Lin, Chair, Professor of Physics
Stephen Knowlton, Professor of Physics
Edward Thomas, Jr., Professor of Physics
Abstract
Magnetic reconnection, the process by which the magnetic energy is converted
into heat and ?ow energy, is believed to be a very important process in space and lab-
oratory plasmas. Among various reconnection models, the Petschek model suggests a
promising explanation for the fast time scale phenomena caused by the reconnection.
For the antiparallel (zero guide ?eld) cases, the Petschek model contains a small di?u-
sion region and four standing slow shocks that bound the out?ow region. In contrast
to the pure antiparallel case, theoretical and simulation studies show that in general
cases with a ?nite ?eld, the Petschek model is modi?ed with complicated wave struc-
tures. Such models, however, only correspond to the two-dimensional (2D) magnetic
reconnection with a long, extended X-line. Two remaining issues in the Petschek
model are still poorly understood. One is the onset mechanism for the reconnection,
and the other is the wave properties in the out?ow region of general three-dimensional
(3D) reconnection. Tearing mode instabilities are thought to play a fundamental role
to trigger the reconnection. In this thesis, we extend the recently-developed gyroki-
netic electron and fully kinetic ion (GeFi) simulation model to nonuniform plasmas,
and use it to study tearing mode instabilities under a ?nite guide ?eld. By remov-
ing the rapid gyro-motions of electrons in the calculation, the main advantage of the
current GeFi model is that the realistic mass ratio between electrons and ions is al-
lowed to be employed in the simulation code. Through an eigenmode analysis, the
improved GeFi model is benchmarked against the linear theory of the tearing instabil-
ity. Furthermore, our simulation results show that the ion kinetics play an important
role in the linear growth rate of the tearing mode instability as well as in the non-
linear saturation of tearing modes. For the cases with multiple tearing modes, the
ii
interactions between them lead to the coalescence of neighboring magnetic islands,
and form larger saturated magnetic islands, which has been suggested as a necessary
condition to trigger a fast reconnection. In order to explore the properties of low
frequency waves generated by the reconnection under a ?nite guide ?eld, a 3D hy-
brid simulation is carried out for the large-scale structure of the reconnection layer.
For the case with an in?nitely long X-line, quasi-steady rotational discontinuities are
formed behind a leading plasma bugle that propagates away from the reconnection
site. Field-aligned structures are also observed in the transition region between the
steady layer and the leading bulge. For the cases with a limited extent of X-line
with length < 30di, the perturbations caused by the reconnection propagate along
the magnetic ?eld line, and the wavefront of propagation has the properties of kinetic
Alfv?en wave. There are no evident discontinuities that form to bound the out?ow
region. Due to the propagation of the waves, a layer of accelerated plasmas is found
beyond the extent of the X-line. The simulation indicates that the wave structures
in the reconnection are greatly modi?ed from that in the 2D Petschek reconnection
model when the X-line has a ?nite length.
iii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
1 An Overview of Magnetic Reconnection . . . . . . . . . . . . . . . . . . . 1
1.1 Steady State Reconnection Models . . . . . . . . . . . . . . . . . . . 3
1.1.1 Sweet-Parker Model . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Petschek Model . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Tearing Mode Instability . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.1 Resistive Tearing Mode . . . . . . . . . . . . . . . . . . . . . . 7
1.2.2 Collisionless Tearing Modes . . . . . . . . . . . . . . . . . . . 10
1.3 The Structure of Reconnection Layer . . . . . . . . . . . . . . . . . . 13
1.3.1 Summary of MHD Waves and Discontinuities . . . . . . . . . 13
1.3.2 Previous Research on the Structure of Reconnection Layer . . 16
1.4 Numerical Models of Magnetic Reconnection . . . . . . . . . . . . . . 18
1.5 Motivations and Objectives of the Thesis . . . . . . . . . . . . . . . . 21
2 Improved GeFi Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.1 Formulations of the GeFi Model . . . . . . . . . . . . . . . . . . . . . 24
2.1.1 Kinetic Equations of Particles . . . . . . . . . . . . . . . . . . 24
2.1.2 Maxwell Equations of GeFi Model . . . . . . . . . . . . . . . . 27
2.2 Simulation Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.1 The ?f Method . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.2.2 Advance Particles . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.2.3 Operators < .. > and ?.. . . . . . . . . . . . . . . . . . . . . . . 35
iv
2.2.4 Field Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3 Benchmark for Uniform Plasmas . . . . . . . . . . . . . . . . . . . . . 37
2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3 Simulations of Tearing Mode Instabilities under a Finite Guide Field . . 40
3.1 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2 Eigenmode Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Benchmark for a Simple Current sheet . . . . . . . . . . . . . . . . . 45
3.4 Linear Simulation of the Tearing Mode in the Harris current Sheet . . 47
3.5 Nonlinear Simulation of Tearing Modes in the Harris current Sheet . 59
3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4 Generation of Low Frequency Electromagnetic Waves by Magnetic Recon-
nection under a Finite Guide Field . . . . . . . . . . . . . . . . . . . . . 67
4.1 3D Hybrid Simulation Model . . . . . . . . . . . . . . . . . . . . . . . 67
4.2 Properties of Kinetic Alfv?en Waves . . . . . . . . . . . . . . . . . . . 70
4.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.3.1 Case with an In?nite X-line: 2?0 = ? . . . . . . . . . . . . . 71
4.3.2 Case with a Finite X-line: 2?0 = 8di . . . . . . . . . . . . . . . 77
4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
v
List of Figures
1.1 A schematic representation of reconnection is shown. Magnetic ?eld lines
become reconnected in the di?usion region and form a X-point. Two
islands are generated and bounded by the separatrix. Picture is cited
from http://www.aldebaran.cz/astrofyzika/plazma/reconnection.html. . 2
1.2 Sweet-Parker model of reconnection is shown with a steady di?usion region
of length 2L and width 2l. vi is the in?ow speed, and vo the out?ow speed.
Picture is cited from http://www.scholarpedia.org/article/reconnection. 4
1.3 Petschek model of reconnection is shown with a small di?usion region.
The out?ow region is bounded by two pair slow shocks. Picture is cited
from http://www.scholarpedia.org/article/reconnection. . . . . . . . . . 5
1.4 Magnetic ?eld-lines in the vicinity of a magnetic island generated by tear-
ing modes is shown. A magnetic island is bounded by separatrix. Picture
is cited from http://farside.ph.utexas.edu/teaching/plasma/lectures1. . . 8
1.5 A simple slab current sheet is shown with antiparallel magnetic compo-
nents Bx and a strong guide ?eld By. Picture is cited from Drake and Lee
[1977a]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.6 Directions ofkandBnear the center of the current sheet is shown. Picture
is cited from Drake and Lee [1977a]. . . . . . . . . . . . . . . . . . . . . 12
1.7 Phase velocity and group velocity of MHD waves. Picture is cited from
http://solar.physics.montana.edu/magara/Research/Topics/MHDwaves. 14
vi
2.1 Coordinates transformation from particle phase space to guiding-center
phase space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.2 ?b??2?b for Bx0By0 = 1 (black) and Bx0By0 = 0.5 (red). The horizontal line is the
normalized z-axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3 The four points gyro-averaging scheme is shown where b is the local mag-
netic direction and r1 is a point on the gyro-circle. . . . . . . . . . . . . 36
2.4 The top plot : a comparison between the dispersion relations of ?Bz for
the fast magnetosonic/whistler branch obtained from the GeFi simulation
(open dots), and the corresponding analytical linear dispersion relation
based on the fully kinetic mode (solid lines). The bottom plot : a compar-
ison between the dispersion relations of ?By for the shear Alfv?en/kinetic
Alfv?en mode branch obtained from the GeFi simulation and the analytical
theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1 The magnetic ?eld (left, and the red solid line on the right) and particle
density (black solid line on the right) of the Harris current sheet with a
tanh and sech pro?le along Z direction, respectively. The magnetic ?eld
consists of two components, Bx0 and By0. Particle density is indicated by
the black line, where the dotted line means the background particle density. 41
3.2 Contour of Ay and eigenfunction of ?Ay with Lz = 10.0?i, Lx = 10.0?i,
L = 0.5?i, By0/Bx0 = 10, Ti = Te, and C1 = ?0.14vthe. . . . . . . . . . . 45
3.3 Relationship between the linear growth rate and Bx0/By0 with Lz =
10.0?i, Lx = 10.0?i, L = 0.5?i, Ti = Te, and C1 = ?0.14vthe . . . . . . . . 46
3.4 The linear growth rate as a function of kxL (top) and L/?s (bottom)
with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 (top), and kxL = 0.41
(bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
vii
3.5 The contours of perturbations magnetic ?eld and the corresponding eigen-
functions with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 kxL = 0.5. . . 49
3.6 The e?ect of inhomogeneity with cases Bx0/By0 = 0.05 (left column) and
Bx0/By0 = 1 (right column), for Ti/Te = 0.01 and L = 0.5?s. The top
row : eigenfunctions ?Ay (solid lines from GeFi simulation, dashed lines
from eigenmode calculation; the middle row : the contour of ?Ay from
simulations; the bottom row : the contour of ?Ay with ?H?B turned o?. . 50
3.7 The growth rate as a function of Ti/Te for three cases L = ?s,1.5?s,3.0?s
with Bx0/By0 = 0.05, kxL = 0.5. . . . . . . . . . . . . . . . . . . . . . . 51
3.8 The growth rate as a function of kxL for three di?erent guiding ?elds
Bx0/By0 = 5.0, 3.0, 1.0 with Ti/Te = 1.0, L = 0.5?s. . . . . . . . . . . . . 53
3.9 The growth rate as a function of L with Ti/Te = 1.0, Bx0/By0 = 5.0,
kxL = 0.41. ? : cases with ions; ? : cases without ions; solid lines: the
growth rate ratio between cases with ions and cases without ions. . . . . 55
3.10 The growth rate as a function of Ti/Te with L = 0.5?s, Bx0/By0 = 3.0,
kxL = 0.41. ? : with electrostatic e?ect; ? : without electrostatic e?ect;
the solid line : the growth rate ratio between them. . . . . . . . . . . . . 56
3.11 The contours of typical perturbation structure with L = 0.5?s, Bx0/By0 =
3.0, kxL = 0.41, Ti/Te = 2.0. . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.12 The comparison among tearing mode growth rate as a function of Bx0/By0
obtained from GeFi simulation (open circles), the eigenmode calculation
(solid line) and Drake and Lee analytical asymptotic matching Eq.3.23
(dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
viii
3.13 The time evolution of the magnitude of eigenfunction A? (left) and the
contour plot of Ay at time t = 20??1i (right) with L = 0.5?s, Bx0/By0 =
10.0, Ti/Te = 1.0, and kxL = 0.41. . . . . . . . . . . . . . . . . . . . . . 61
3.14 The time evolutions of the magnitude of eigenfunction ne (top left), ni (top
right), A? (bottom left), and the contour plot of Ay at time t = 20??1i
(bottom right) including ions. . . . . . . . . . . . . . . . . . . . . . . . . 62
3.15 Contours of vector potential Ay at di?erent times for the case with four
unstable modes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.16 The time evolution of the magnitude A? for the four modes. . . . . . . . 65
4.1 The contour plots of Bx, By, Bz, Vx, Vy, Vz, T, N, and con?guration of
magnetic ?eld lines. The solid line corresponds to the position at z = 40di,
along which there appears a quasi-steady structure. The solid line at
z = 40di indicates the boundary of the leading bugle. The solid line
z = 40di shows the transit boundary between the plasma bugle and the
steady state layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.2 The spatial pro?les of Bx, By, Bz, N, Vx, Vy, Vz, T?, T?, J?, and B at
z = 40di and hodograms of the tangential magnetic ?eld (Bz-By). . . . . 73
4.3 Top ?gures: the contours of variations of Bx and N along the y direction
in the x-z plane at y = 0; bottom ?gures: the contours of variations of Bx
and N along the y direction in the y-z plane at x = 4di. The solid line
along z = 52di in the top left ?gure indicates the boundary of the leading
bugle. The solid lines in the top ?gures indicate the y-z plane at x = 4di;
the solid lines in the bottom ?gures show the local initial magnetic ?eld
direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
ix
4.4 The contours of variations along the y direction of Bx, By, Bz, N, Vx, Vy,
Vz, Ex, Ey, Ey, J?, E? in the y-z plane at x = 3di when t = 150??1i . . . . 76
4.5 The con?guration of magnetic ?eld lines near the current sheet at t =
50??1i (the top left) and t = 50??1i (the top right), and the contour plots
of Bx (the bottom left) and Vx (the bottom right) in the xz plane of
y = 25di at t = 50?i. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.6 The relationship between time t and the distance L from the center of the
di?usion region (0,0,0) to the position of wave front. . . . . . . . . . . . 79
4.7 The contour plots of Bx and density N in the x-z plane of y = 16di at
t = 80?i. The red dot stands for the position (4, 16, -20). . . . . . . . . 80
4.8 The time evolution of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ez, E?, B,
Ex?b, B?, and the hodogram of Bx ?Bx?b at the position (4, 16, -20) as
indicated by the red dot in Fig.4.7. . . . . . . . . . . . . . . . . . . . . . 81
4.9 The contour plots of Vz at x = 0 (left) and x = 2di (right) at t = 80??1i . 82
4.10 The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane
of y = 0 when t = 150??1i . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.11 The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z plane
of y = 20di when t = 150??1i . . . . . . . . . . . . . . . . . . . . . . . . 85
x
List of Tables
1.1 Structure of reconnection layer under various conditions [Lin and Lee,
1994a]. SS: slow shock, SE: slow expansion wave, RD: rotational discon-
tinuity, CD: contact discontinuity, IS: intermediate shock, TIDS: time-
dependent intermediate shock. . . . . . . . . . . . . . . . . . . . . . . . . 17
xi
Chapter 1
An Overview of Magnetic Reconnection
Magnetic reconnection is a process that changes the topology of magnetic ?eld
lines. It is believed to take place between two plasma regions with antiparallel mag-
netic ?eld components. During magnetic reconnection, sheared magnetic ?eld lines
will be e?ectively annihilated through magnetic di?usion, and then magnetic energy
is transferred into plasma heat and kinetic energy. Magnetic reconnection provides
the free energy for phenomena such as solar ?ares, magnetosphere substorms, the
tokamak sawtooth collapse, and the ion heating in reversed-?eld pinch [Yamada et
al., 2010]. Fig.1.1 is a schematic representation of magnetic reconnection. antipar-
allel magnetic ?eld lines break at the di?usion region and form an X-point at the
center. Two magnetic islands are generated and bounded by the separatrix. The
in?ow plasmas are accelerated and ?ow out from the di?usion region.
In a plasma characterized by the ideal MHD Ohm?s law E+V?B = 0, where
V is the velocity, E and B are the vector electric and magnetic ?elds, magnetic ?eld
lines follow the plasma motion as if they were ?frozen in? the ?ow. This property
conserves the ?eld line topology, and thus reconnection can never happen. Varieties
of mechanisms can lead to a violation the ideal Ohm?s law. Resistivity is the most
commonly cited dissipation that causes reconnection. Models based on resistivity
include the Sweet-Parker model [Parker, 1957; Sweet, 1958], in which the extent of
the dissipation layer in the out?ow direction is much bigger than the width of the
layer, and the Petschek model [Petschek, 1964], in which the dissipation layer is highly
localized in both the in?ow and out?ow direction. Theoretical calculations show the
magnetic reconnection rate, a direct measure of the speed of the mixing process, of
1
Figure 1.1: A schematic representation of reconnection is shown. Magnetic
?eld lines become reconnected in the di?usion region and form a X-point. Two
islands are generated and bounded by the separatrix. Picture is cited from
http://www.aldebaran.cz/astrofyzika/plazma/reconnection.html.
the Sweet-Parker model is too slow to explain the energy release processes that are
thought to be caused by the magnetic reconnection. The Petschek model contains
a set of four MHD discontinuities which in the compressible ?ow may be identi?ed
as slow mode shocks across which the plasma is accelerated into the outgoing ?ow
regions. Instead of dissipating magnetic energy, the di?usion region in the Petschek
model mainly plays a role as the perturbation source of slow modes. The Petschek-
type reconnection is the ?rst fast reconnection model, in which the reconnection
rate is fast enough to explain observed dynamic times of phenomena caused by the
reconnection.
To validate the Petschek model, two questions should be answered. One is what
is the physics in the dissipation region, i.e., the triggering mechanism of reconnection;
the other is what is the structure of out?ow region referred as the reconnection layer.
The tearing mode instability is one candidate to trigger the magnetic reconnection in
2
both the resistive MHD regime [Furth et al., 1963, 1973] and the collisionless kinetic
regime [Drake & Lee, 1977a, 1977b]. On the other hand, theoretical and numerical
studies [Lin and Lee, 1994a, 1994b] showed that the structure of the reconnection
layer is not as simple as predicted by Petschek but varies with di?erent parameters.
This thesis is composed of two parts, which addresses the kinetic physics of
magnetic reconnection on two di?erent spatial and temporal scales, from the di?usion
region to the large-scale reconnection layer. First, we are going to improve the exiting
gyrokinetic electron and fully kinetic ion simulation (GeFi) model [Lin et al., 2005a;
Wang et al., 2010], and use it to study the tearing mode instability. Second, a three-
dimensional (3D) hybrid model [Lin et al., 2005b; Lin et al., 2010] will be used to
investigate the electromagnetic waves generated by the reconnection. In the following,
we ?rst give a brief review of the steady state reconnection models, the tearing mode
instability, and the structure of the reconnection layer, is given. In the end, the results
of this thesis are summarized.
1.1 Steady State Reconnection Models
1.1.1 Sweet-Parker Model
The Sweet-Parker model is sketched in Fig.1.2. Sweet [1957] and Parker [1958]
determined the speed vi with which ?eld lines are carried into a steady di?usion region
of length 2L and width 2l. First of all, for a steady state, the convective in?ow of the
magnetic ?eld is balanced by the outward di?usion, so that
vi = ?l, (1.1)
where ? is the magnetic di?usivity.
3
Figure 1.2: Sweet-Parker model of reconnection is shown with a steady di?usion
region of length 2L and width 2l. vi is the in?ow speed, and vo the out?ow speed.
Picture is cited from http://www.scholarpedia.org/article/reconnection.
Second, the rate 4?Lvi at which mass is entering the sheet must equal the rate
4?lvo at which it is leaving, so that, if the density is uniform,
Lvi = lvo, (1.2)
where vo is the out?ow speed. The width l may now be eliminated from Eqs.(1.1)
and (1.2) to give the in?ow speed as v2i = ?vo/L. If the plasma is accelerated along
the sheet by a Lorentz ?j ? ?B force, the out?ow speed vo is the Alfv?en speed at the
in?ow, namely, vo = vAi = Bi/?4?? and the reconnection rate is
vi = vAiR
mi 1/2
(1.3)
where the magnetic Reynolds number Rmi = LvAi/?, or in dimensionless form Mi =
R?1/2mi , with Mi = vi/vAi being the in?ow Alfv?en Mach number.
In the Sweet-Parker mechanism, we identify the sheet length L with the global
length and Rmi therefore with the global magnetic Reynolds number. Since in practice
4
Rmi ? 1, the reconnection rate is very small: for instance, in the solar corona where
Rmi lies between 106 and 1012, the ?elds reconnect at between 10?3 and 10?6 of the
Alfv?en speed. Note the phenomena caused by the reconnection such as solar ?ares
usually have the Alfv?en time scale ?A = 1/vA where vA is the local Alfv?en speed.
Important extra e?ects are the compressibility, which slows down the out?ow
when ?o > ?i, and a pressure gradient along the di?usion region, which can also slow
the out?ow when the out?ow pressure po is large enough.
1.1.2 Petschek Model
Figure 1.3: Petschek model of reconnection is shown with a small di?usion re-
gion. The out?ow region is bounded by two pair slow shocks. Picture is cited from
http://www.scholarpedia.org/article/reconnection.
The con?guration of the Petschek model is sketched in Fig.1.3. Two waves
(slow mode shocks) stand in the ?ow on either side of the out?ow region, where the
direction of B reverses, marking the boundaries of the plasma out?ow regions. A
small di?usion region still exists, but has much smaller length compared with that
in the Sweet-Parker model. The upstream plasma ?ows into the shock where the
5
magnetic energy is converted into the energy of out?owing plasma jets; the outward
wave propagation velocity in the ?z direction balances the plasma in?ow velocity
in a steady state. The reconnection rate should be equal to the shock propagation
velocity which is independent of the plasma resistivity and can be much larger than
theSweet-Parkermodel. Atthesametime, themagneticdi?usionisstillthedominant
mechanism of energy conversion in the di?usion region near the center z = 0 of the
reconnection site. The match of reconnection rate between the di?usion region and
the shock region leads to a very narrow di?usion length.
In the switch-o? limit of the slow shock, the in?ow speed can be estimated as
vzi = vphi = |Bzi|?4??
i
, (1.4)
where vphi is the phase speed of slow shocks, Bzi is the z component of magnetic
?elds in the in?ow region and is approximated as the normal ?eld of the shocks. The
subscripts ?i? and ?o? stand for the in?ow (upstream) and the out?ow (downstream).
Note that the switch-o? limit corresponds to the strongest slow shock, across which
the tangential magnetic ?elds change to 0. To simplify the calculation, we assume
that the plasma is incompressible, i.e., ? = constant. The MHD mass momentum of
conservation laws give,
vzix = vxo?(x), (1.5)
d
dx(?v
2
xo?(x)) =
BziBi
4? , (1.6)
where Bi is the magnitude of magnetic ?eld in the in?ow region. Insert Eq (1.5) into
Eq (1.6), and do a little algebra,
d
dx(
x2
?(x)) =
1
M2Ai
Bzi
Bi =
Bzi
Bi , (1.7)
6
where MAi the upstream Mach number and MAi = 1 for the switch-o? limit of slow
shocks. After integrating Eq.1.7, we have
x
?(x) =
Bzi
Bi . (1.8)
Since the downstream velocity can be expressed as
vxo = vAi, (1.9)
we have vzi ? BziBi vAi. It has indicated BziBi ? 0.1 by simulations. Therefore the
reconnection rate is vziv
Ai
? 0.1, much faster than in the Sweet-Parker model.
There has been some controversy whether the Petschek model describes the cor-
rect physics. It gives a description of reconnection that is only very weakly dependent
on the properties of the reconnection region and thus on the resistivity. Numerical
simulations [Biskamp, 1986], however, showed that with a spatially constant resistiv-
ity the Petschek-type reconnection does not occur. In order to have the Petschek-type
reconnection, the resistivity should be enhanced in the di?usion region [Lottermoser
et al., 1997]. For collisionless plasmas as in most of the space environments, the de-
termination of the anomalous ? in the di?usion region is a hot topic in reconnection
research. Current driven instabilities such as the modi?ed two-stream instability [Wu
et al., 1983] the low-hybrid drift instability [Daughton, 2003; Ricci et al., 2005] may
play a key role in the formation of anomalous resistivity.
1.2 Tearing Mode Instability
1.2.1 Resistive Tearing Mode
In the resistive (i.e., collisional) plasma, a current sheet tends to di?use at a
slow rate with a time-scale of ?d = l2/?, where 2l is the width of the current sheet
7
and ? = ???1 is the magnetic di?usivity. During the process of magnetic di?usion,
magnetic energy is converted into heat at the same slow rate by the Ohm?s law. As
pointed out in the ?rst section, however, the magnitude of ?d is often far too large to
explain the time scale of dynamical physical processes.
Furth [1963] ?rst found that the di?usion can drive an instability with a time
scale much faster than the di?usion time, referred as the tearing mode instability.
Considering the perturbation as the one in Fig.1.4, the tension force of reconnected
?eld lines tends to pull new loops of ?eld up and down away from the X-point, while
the magnetic pressure gradient tends to push plasma in from the sides towards the
X-point. The ?eld lines at the sides are curved and so possess a restoring magnetic
tension force.
Figure 1.4: Magnetic ?eld-lines in the vicinity of a magnetic island generated by
tearing modes is shown. A magnetic island is bounded by separatrix. Picture is cited
from http://farside.ph.utexas.edu/teaching/plasma/lectures1.
8
Detailed theoretical studies showed [Furth, 1963; Wesson, 1966] that the tearing
mode occurs with a wavelength greater than the width of the sheet (kl < 1). The
smallest allowable wavelength grows in a time ?3/5d ?2/5A , whereas the longest wave-
length has the fastest growth rate ?,
? = 1??
d?A
. (1.10)
The presence of a guiding ?eld does not a?ect the existence of tearing modes.
At any particular location, the maximum growth of tearing modes occurs where the
wave vector k is perpendicular to the equilibrium ?eld B0, that is,
k?B0 = 0. (1.11)
With the addition of a uniform guiding ?eld, this condition can be satis?ed at the
center of the current sheet.
One may have noticed that the typical time scale of tearing modes ?t = ??d?A is
in general much larger than the di?usion time ?d, but much smaller than the Alfv?en
time ?A, or the typical fast reconnection time scale, i.e., ?A << ?t << ?d. So the
tearing mode is still not fast enough to explain the phenomena that are thought to
be caused by the reconnection.
The geospace environment magnetic (GEM) reconnection challenge was to aim
to discover the physics of fast reconnection [Birn et al., 2001]. They found that
nonlinear evolution of a large magnetic island in a current sheet can lead to a fast
reconnection with the rate of the order ??1A . The Hall terms in the generalized or
kinetic Ohm?s law were believed to be important. For high ? plasmas it was clari?ed
that the whistler wave plays a key role in the fast reconnection physics [Pritchett,
2001; Birn et al., 2001]. At lower ? plasmas in the presence of a guide ?eld, it was
9
shown that the kinetic Alfv?en wave drives the fast reconnection [Ricci, 2008; Rogers,
2001].
However, the GEM reconnection challenge avoided the issue on the triggering
of the fast reconnection, or more precisely on how a magnetic island forms initially
with the size comparable to the system?s spatial scale. One may argue that the role
of tearing modes may present a mechanism that form a large magnetic island. Un-
fortunately, the tearing mode often saturates at a small level, and the width of the
magnetic island is not large enough to trigger the fast reconnection [Rutherford, 1973;
Daughton and Karimabadi, 2005]. However, it is also found that a chain of magnetic
islands created by multiple tearing modes can lead to an instability, namely the coa-
lescence [Finn and Kaw, 1977; Pritchett, 1992]. During the coalescence, neighboring
islands attract each other due to the parallel currents inside islands and then form
larger magnetic islands [Karimabadi et al., 2005a; Pritchett, 2005].
1.2.2 Collisionless Tearing Modes
Plasmas in the cosmic and tokamak environments are generally collisionless. So
the resistive MHD approach to the tearing mode may be invalid. Drake and Lee
[1977a] applied a kinetic theory to the tearing mode instability for a current sheet
with a strong guide magnetic ?eld (perpendicular to the antiparallel components of
the ?eld). They found a perturbed current produced by the induced electric ?eld
around the central layer can drive the tearing mode instability, and thereby cause the
magnetic reconnection.
To estimate the growth rate of the collisionless tearing mode, the simple slab
current con?guration as in Fig.1.5 was studied. First, consider a perturbation with
the wave number kx, and the growth rate ?. The perturbed magnetic ?eld can be
10
Figure 1.5: A simple slab current sheet is shown with antiparallel magnetic compo-
nents Bx and a strong guide ?eld By. Picture is cited from Drake and Lee [1977a].
expressed with the vector potential Ay as
Bz = ikxAy, (1.12)
Bx = ?Ay?z . (1.13)
Then the induced electric ?eld can be written as
Ey = ??Ay. (1.14)
It has a component parallel to the local magnetic ?eld E|| = Eycos(?), where sin(?) =
B0x/B0 = k||/kx(see Fig.1.6). The parallel electric ?eld will generate a large parallel
current at z = 0 where
k?B0 = 0. (1.15)
11
Figure 1.6: Directions of k and B near the center of the current sheet is shown.
Picture is cited from Drake and Lee [1977a].
Note that the induced current is restricted inside a narrow region of |z| < ?.
After integrating the z component of the Ampere?s equation, we have
[Bx]??? = 4?c ?jy, (1.16)
which can be rewritten as
[?Ay?z ]??? = ?4?c ?jy. (1.17)
So ?Ay?z is not continuous across the induced current layer. Outside of the current
layer, E|| = 0, and the ideal MHD equations was used. A simple expression is given
as
?? = 1A
y
[?Ay?x ]???. (1.18)
Then combining the above two equations, one has
??Ay = ?4?c ?jy. (1.19)
12
The induced current is
?y = ?n0e
2Ay
me . (1.20)
According to the cherenkov resonances condition, the growth rate should be
? = k|| ?vthe = kxvthe?l
s
(1.21)
where ls is the magnetic ?eld shear length, and vthe is the electron thermal velocity.
A little algebra gives the growth rate as
? = ?
?k
xvthe
lsd2e , (1.22)
where d?1e = c/?pe is the electron skin depth.
1.3 The Structure of Reconnection Layer
1.3.1 Summary of MHD Waves and Discontinuities
Startting from the ideal MHD equations, one can easily deduce the dispersion
relationship of the linear MHD waves, which are summarized below [Landau and
Lifshitz., 1960].
(a) Fast mode:
CF = 12
?
[(C2s +C2A) +
?
(C2s +C2A)?4C2sC2I];
(b) Intermediate mode (Alfv?en wave):
CI = CAcos(?);
(c) Slow mode:
CS = 12
?
[(C2s +C2A)?
?
(C2s +C2A)?4C2sC2I];
where C denotes the wave phase speed, and the subscript F denotes the fast mode, I
the intermediate mode, and S the slow mode. Cs =
??P
? is the sound speed and ? is
the speci?c heat ratio. CA = B?4?? is the Alfv?en speed and ? is the wave propagation
13
angle with respect to the background magnetic ?eld. Fig.1.7 shows the dispersion-
relation of MHD waves. Phase velocity is the propagation speed of a single-mode
wave while group velocity is the envelope speed of the waves, i.e., the propagation
speed of a group of single-mode waves. From Fig.1.7, one can tell that the fast mode
is isotropic, propagating around all angles with respect to the background magnetic
?eld, while the propagations of the intermediate mode and the slow mode concentrate
in the direction parallel to the background magnetic ?eld.
Figure 1.7: Phase velocity and group velocity of MHD waves. Picture is cited from
http://solar.physics.montana.edu/magara/Research/Topics/MHDwaves.
An MHD discontinuity can be considered as a thin transition region between
two uniform stationary plasma regions across which the magnetic ?eld, plasma den-
sity, pressure and ?ow velocity may have a signi?cant jump. The jump conditions
are obtained from the conservation laws of the MHD formulations referred to as the
Rankine-Hugoniot (RH) jump conditions. It is well known that four types of MHD
discontinuities exist: contact discontinuity, tangential discontinuity, rotational dis-
continuity and MHD shocks. The jump conditions of di?erent discontinuities are
summarized below [Hans and Poedts, 2004].
14
(1) Contact discontinuity:
vn = 0, Bn ?= 0, [Bt] = 0, [vt] = 0, [?] = 0, [P] = 0.
A contact discontinuity separates two regions with di?erent plasma densities but
the pressure, magnetic ?eld, and plasma velocity are continuous. The component of
magnetic ?eld normal to the discontinuity front, Bn, is not equal to zero. The contact
discontinuity is a non-propagating structure.
(2) Tangential discontinuity:
vn = 0, Bn = 0, [Bt] ?= 0, [vt] ?= 0, [?] ?= 0, [P + B28? ] = 0.
A tangential discontinuity links two regions with di?erent tangential magnetic
?elds. The component of magnetic ?elds normal to the discontinuity front Bn is
equal to zero. No ?ow normal to the discontinuity is allowed. Across the tangential
discontinuity the total pressure is balanced , while the plasma densities and tangential
plasma ?ow velocities on the two sides of the discontinuity are di?erent.
(3) Rotational discontinuity:
vn = Bn?4??, [vt] = [Bt]?4??, [?] = 0, [P] = 0.
A rotational discontinuity links two regions with di?erent tangential magnetic
?eld and a ?nite normal magnetic ?eld component, and the normal ?ow is allowed.
The rotational discontinuity is associated with the propagating nonlinear Alfv?en mode
structure through which the normal component of the plasma ?ow velocity relative
to the discontinuity is constant and equal to the normal Alfv?en speed, which is vAn =
Bn?
4??. Across a rotational discontinuity, the plasma density, the pressure, the normal
magnetic ?eld and the normal plasma ?ow velocity remain unchanged.
MHD shocks can be further divided into three types: fast shocks, intermediate
shocks, and slow shocks. They are formed by the steepening of MHD fast waves,
intermediate waves, and slow waves, respectively. The jump conditions of MHD
shocks can be described as the following [Priest and Forbes, 2006] :
15
(1) For a fast shock vn1 > CF1 and CI2 < Vn2 < CF2, the tangential magnetic
?eld does not change direction, and [?] > 0,[P] > 0,[Bt] > 0,[vn] < 0. Across the
fast shock, the plasma density, the pressure, the magnitude of the tangential magnetic
?eld increase. The normal component of plasma ?ow velocity decreases.
(2) For an intermediate shock vn1 > CI1 and vn2 < CI2, the tangential magnetic
?eld changes the direction by 1800, and [?] > 0,[P] > 0,[Bt] ?= 0,[vn] < 0. Across
an intermediate shock, the plasm density and the pressure increase. The normal
component of the plasma ?ow velocity decreases. The magnitude of the tangential
magnetic ?elds is changed.
(3) For a slow shock with Cs < vn1 < CI1 and vn2 < Cs2, the tangential magnetic
?eld does not change the direction, and [?] > 0,[P] > 0,[Bt] < 0,[vn] < 0. Across a
slow shock, the plasma density and pressure increase. The normal component of the
plasma ?ow velocity and the magnitude of the tangential magnetic ?eld decrease.
MHD theory, however, does not give the internal structure of discontinuities.
With the particle kinetic e?ects in a collisionless plasma, the discontinuities discussed
above have a ?nite width on the order of several ion radii [Swift and Lee, 1983; Richter
and Sholer, 1989, 1991]. Moreover, MHD theory fails to explain the dissipation mech-
anism across MHD shocks in collisionless plasmas. Numerous numerical experiments
have con?rmed that the dissipation of MHD shocks in collisionless plasmas arise from
the self consistent wave-particle interactions inside of shocks [Sato, 1979; Ugai, 1984;
Lee and Lee, 1991].
1.3.2 Previous Research on the Structure of Reconnection Layer
Based on the ideal MHD formulation, the structure of the reconnection layer
was studied by solving the one-dimensional (1D) Riemann problem [Heyn et al.,
1988; Biernat et al., 1989; Lin and Lee, 1994a]. In this initial value problem, a
nonzero normal component of the magnetic ?eld is superposed on a 1D current sheet,
16
By = 0 By ?= 0
Symmetric SS + SS? RD + SS + SS? + RD?
Asymmetric RD + SE + CD + SS? RD + SE + CD + SS? + RD?
Table 1.1: Structure of reconnection layer under various conditions [Lin and Lee,
1994a]. SS: slow shock, SE: slow expansion wave, RD: rotational discontinuity, CD:
contact discontinuity, IS: intermediate shock, TIDS: time-dependent intermediate
shock.
which breaks the equilibrium of the current sheet. The initial sheet then evolves to
a layer with several MHD discontinuities, shocks. The temporal development of this
layer is viewed as an approximation of the reconnection layer at increasing distances
from the X-line. These studies indicate that ?ve discontinuities and expansion waves
may develop in a reconnection layer, including rotational discontinuities, intermediate
shocks, slow shocks, slow expansion waves, and a contact discontinuity. The rotational
discontinuities and intermediate shocks change the direction of the magnetic ?eld,
the slow shocks and slow expansion waves change the magnitude of magnetic ?eld
and plasma density, and the central contact discontinuity matches di?erent plasma
densities on the two sides. In addition, two fast expansion waves are also present
in the solution of the Riemann problem, but they quickly propagate away from the
reconnection layer. The results of the Riemann problem are summarized in Table.1.1.
The MHD simulations have been carried out for the resistive structure of recon-
nection layers. Using a 1D resistive MHD code, Lin et al. [1992], Lin and Lee [1994a],
Ugai and Shimizu [1995] have simulated the Reimann problem for the evolution of
an initial current system. It is found that the rotational discontinuities are replaced
by time-dependent intermediate shocks. To deal with the multi-dimensional recon-
nection con?guration, two-dimensional (2D) MHD simulations have been carried out
to study the steady-state reconnection. The MHD discontinuities predicted from the
1D simulations are found to develop in the 2D reconnection layer [Shi et al., 1990;
Sholer, 1989; Yan, 1992]. Slow shocks in the Petschek-type symmetric reconnection
17
layer and intermediate shocks in the reconnection layer with asymmetric magnetic
?eld and plasma density have been obtained in the cases with By = 0. In the cases
with By ?= 0, rotational discontinuities are replaced by time-dependent intermediate
shocks that have been found.
MHD models and simulations, however, do not account for the particle kinetic
e?ects in collisionless plasmas. In an e?ort to include the ion kinetic e?ects, both 1D
hybrid simulations and 2D hybrid simulations, in which ions are treated as discrete
particles and electrons are treated as a massless ?uid, of the reconnection layers
have been performed. 1D hybrid simulations [Lin and Lee, 1995] indicate that the
contact discontinuity disappears because of the mixing of ions along the ?eld lines;
the slow shock and slow expansion waves are modi?ed. The kinetic structure of the
reconnection layer have also been investigated by 2D large-scale hybrid simulations
[Lin and Swift, 1996; Lottermoser et al., 1998; Kruassvarban and Omidi, 1995]. These
hybrid simulations have shown that the structures of MHD discontinuities and shocks
are modi?ed in collisionless plasmas..
Note that previous simulations of the reconnection layer are limited to the 1D
or 2D geometry. In these cases, the normal direction to the fronts of low-frequency
electromagnetic waves are mainly along the direction of the normal magnetic ?eld,
whichisalmostperpendiculartothecurrentsheet. Intherealspaceplasmas, however,
magnetic reconnection is frequently observed to have a 3D structure, for which the
length of the X-line has a short or ?nite extent. To further understand the physics
of the reconnection, we will address the 3D e?ects of the reconnection layer in this
thesis.
1.4 Numerical Models of Magnetic Reconnection
Numerical simulations have proven to be a powerful tool to understand the
physics of magnetic reconnection. Simulation codes can be rougthly divided into
18
four categories: (1) MHD codes, (2) fully-kinetic particle codes, (3) hybrid codes that
treat ions as particles and electrons as a massless ?uid, and (4) gyrokinetic codes that
treat charged particles as gyrokinetic particles.
In the MHD codes, MHD equations are solved numerically with ?nite di?erence
or ?nite elements. Since MHD simulations do not contain the individual particle
dynamics, they only describe the large-scale dynamics of plasmas. However, the
MHD model may still give a correct description of physics on the large-scale bulk
dynamics of collisionless plasmas. There are two reasons for that. First, in the
presence of a magnetic ?eld particles undergo gyro-motions around the local magnetic
?eld, and thus the motions across ?eld lines are prevented. Second, along the magnetic
?eld direction numerous short-scale wave-particle interactions tend to impede the
motion of charged particles. MHD simulations of the magnetic reconnection are
mainly focused on the global structure of reconnection. It usually has a much larger
spatial scale than di or ?i, and much larger time scale than 1/?pi or 1/?i, where di is
the ion inertial length, ?i is the ion gyroradius, ?pi is the ion plasma frequency, and
?i is the ion gyrofrequency. To incorporate the local dynamics due to the separation
between the electron and ion masses, the so-called Hall MHD codes and two ?uid
MHD codes have been developed. Nevertheless, kinetic e?ects are again missing in
these models.
Full-particle simulations, in which electrons as well as ions are treated as fully
kinetic particles, have been utilized for decades to investigate the kinetic behavior
of reconnection [Terasawa, 1981; Leboeuf et al., 1982; Hoshino, 1987; Ding and Lee,
1990; Mandt et al., 1994; Cai and Lee, 1997; Hesse and Winske, 1998; Shay and Drake,
1998; Shay et al., 2001, 2007; Pritchett, 2001; Pritchet and Coroniti, 2004; Swisdak
et al., 2005; Daughton et al., 2006; Karimabadi et al., 2007]. In these codes, the
particle-in-cell (PIC) technique is often used in which the simulated ?macroparticles?
represent many plasma particles. Nevertheless due to large-scale separations both in
19
time (from electron plasma oscillation to reconnection time) and in space (from the
Debye length to the system size), various approximations or compromises have to
be used. Most of full-particle simulations have employed unrealistically low ion-to-
electron mass ratio mi/me, and also limited the domain to a few or ten?s of ?i, and
the simulation time much less than the global Alfv?en time scale.
Hybrid simulation is another kinetic approach often used in reconnection studies,
in which ions are treated as particles, but electrons are treated as a massless resistive
?uid. In general, hybrid codes do not resolve the small spatial and time scales asso-
ciated with the electron mass, and have been used to simulate large-scale structures
of reconnection associated with the ion dynamics [Lin and Swift, 1996; Nakamura
and Scholer, 2000; Petschek, 1964; Lin and Lee, 1994a; Lin and Wang, 2006]. Note
that the triggering mechanism of collisionless reconnection, due to the electron kinetic
e?ects, is not included in hybrid models.
An alternative particle model, called gyrokinetic model, has been broadly used in
the fusion research. For the low frequency physics in strong magnetized plasmas such
as a tokamak plasma, the gyromotion can be averaged out. The nonlinear gyrokinetic
Vlasov-Maxwell equations have been derived by Fireman and Chen [1982] and ?rst
applied in gyrokinetic simulations by Lee [1983]. Since the fast gyromotion of charged
particles have been removed, larger spatial and time steps are allowed in gyrokinetic
simulations than in traditional particle simulations. In the presence of a guide ?eld,
particles are strongly magnetized during the magnetic reconnection. Therefore the
gyrokinetic treatment can also be applied to the reconnection research. Tearing mode
instabilities have been simulated using the gyrokinetic model. Wan et al. [2007]
treated ions as gyrokinetic particles and electrons as drift kinetic particles, while
both ions and electrons are gyrokinetic particles in Ricci et al.?s work [2004]. But the
gyrokinetic and drift kinetic treatment of ions are not suitable to the reconnection
physics since the whistler and lower-hybrid regime is beyond the model capability.
20
Recently an innovative fully nonlinear new particle simulation model has been de-
veloped, in which the electron dynamics is handled by the gyrokinetic (GK) equations
and the ions are treated as fully kinetic (FK) particles. In the gyrokinetic electron
and fully kinetic ion (GeFi) simulation, the rapid electron cyclotron motion is re-
moved, while ?nite electron Larmor radii, wave-particle interactions, and o?-diagonal
components of the electron pressure tensor are retained, as well as the fully kinetic
ion physics. This treatment results in a larger time step and allows the realistic mass
ratio mi/me in the simulation of magnetic reconnection. The computation e?ciency
can be signi?cantly improved as compared with that of the full-particle codes. The
model is particularly suitable for plasma dynamics with wave frequencies ? < elec-
tron gyrofrequency ?e, and for problems in which the wave modes range from Alfv?en
waves to lower-hybrid/whistler waves. The GeFi scheme has been benchmarked in
the uniform plasma [Lin et al., 2005] and important new results have been obtained
for current sheet instabilities in the presence of a ?nite guide ?eld [Wang et al., 2008;
Yoon et al., 2008]. The original GeFi model, however, does not include the e?ects
of nonuniformity of the background magnetic ?eld and plasmas. In this thesis, we
are going to improve the existing GeFi codes for a current sheet, and study the
collisionless tearing mode instability with the realistic mass ratio me/mi = 1/1836.
1.5 Motivations and Objectives of the Thesis
As described above, full particle simulations includes the kinetic electron dy-
namics but cannot properly handle the realistic mass ratio and the large spatial and
temporal scales associated with the reconnection layers, while hybrid simulations
cannot model the triggering physics of reconnection in a self-consistent manner. Our
approach:
? We utilize the GeFi simulation model, in which the electron dynamics is handled
by the gyrokinetic equations and the ions are treated as fully kinetic particles.
21
? In the GeFi model, the rapid electron cyclotron motion is removed, while ?nite
electron Larmor radii, wave-particle interactions, and o?-diagonal components
of the electron pressure tensor (necessary for the reconnection physics) are re-
tained.
? This treatment results in a larger time step and allows us to treat the realistic
mass ratio mi/me.
? The model is particularly suitable for plasma dynamics with wave frequencies
< ?e, and wave numbers k? << k?.
Analytical theories on collisionless reconnection are very limited to cases with
either a zero guide ?eld or a very large guide ?eld, and neglect many kinetic e?ect.
On the other hand, because of its low frequency properties, the exiting fully particle
simulations on tearing modes often employ mi/me < 100.
Our approach:
? Although the GeFi model is promising to handle the electron and ion scales
with the realistic mass ratio in one code essential for the reconnection research,
it needs to be improved to include the e?ects of the background inhomogeneity
because the wavelength of the tearing mode is comparable to the current sheet
width. We will modify the GeFi model.
? We derive the eigenmode equations of the tearing modes based on the drift
kinetic approximation. By solving them, we benchmark the GeFi model in the
nonuniform plasmas.
? Using the GeFi model, various kinetic e?ects on the linear growth of tearing
modes are investigated.
? Nonlinear evolutions of both a single tearing mode and multiple tearing modes
are studied by the GeFi simulations.
22
While the full development of the GeFi model for the nonlinear fast reconnection
is still underway, we also address in this thesis the 3D physics of large-scale kinetic
structure of reconnection with the hybrid model. In the Petschek model, the recon-
nection layer is bounded by shocks and discontinuities. Previous simulations show
the ion kinetic e?ect plays an important role on the structure of reconnection layer.
But these simulations are limited to 1D or 2D con?gurations. The large-amplitude
wave/discontinuity structures in the general cases with a ?nite X-line length have not
been studied.
Our approach:
? We extend the previous hybrid simulations to 3D and investigate the 3D physics
of the reconnection layer.
? Using the 3D hybrid code, the ?nite length X-line e?ects on the generation of
low frequency electromagnetic waves are studied.
23
Chapter 2
Improved GeFi Model
Magnetic reconnection involves a wide range of spatial and temporal scales, and
is usually triggered by micro instabilities. The reconnection changes the topology
of magnetic ?elds on a global scale. Aiming at incorporating multi scale physics of
reconnection in one code, the gyrokinetic electron and fully kinetic ion (GeFi) particle
simulation scheme has been developed [Lin et al., 2005; Wang et al., 2008]. While
the original scheme has been benchmarked for uniform plasmas, further improvement
of GeFi is still required for the simulation of magnetic reconnection, in which the
interested wavelengths are usually comparable to the scale length of the current sheet
nonuniformity. In this chapter, a description of the improved GeFi model is given
?rst. Then I will discuss the algorithm of the GeFi numerical scheme. Finally, the
improved model will be validated for uniform plasmas.
2.1 Formulations of the GeFi Model
2.1.1 Kinetic Equations of Particles
In the GeFi model, ions are treated as fully kinetic particles, and electrons are
modeled by gyrokinetic particles. The Vlasov equation is still valid for ions, while for
electrons we transform the phase space variables to guiding-center variables. Di?erent
from the original GeFi scheme in Lin et al. [2005], both ions and electrons are
advanced by the scalar and vector potential ? and A.
1.Ions
24
In the limit of collisionless plasmas, the time evolution of ion distribution function
is totally determined by the Vlasov equation:
?fi
?t +
dxi
dt ?
?fi
?xi +
dpi
dt ?
?fi
?pi = 0, (2.1)
where pi is the ion canonical momentum. For a singe ion, its motion follows
dxi
dt = vi = (pi ?qA/c)/mi, (2.2)
dpi
dt = ?q ??(??xi ?A/c). (2.3)
2. Electrons
The particle phase space variables (x, v, t) can be related to the guiding-center
phase space variables (R, ?,?,?, t) through the following guiding-center transforma-
tion
R = x+ ??1v?e? (2.4)
? = mv?
2
2B (2.5)
? = 12mv2 +q? = 12m|p? ecA|2 +q?, (2.6)
where ? is the particle gyrofrequency, ? is the magnetic moment, ? represents the
Hamiltonian particle energy in the electromagnetic ?eld, ? is the gyro-angle. e? and
e?, respectively, are the parallel and perpendicular units vectors relative to the equi-
librium magnetic ?eld. The relationship between R and x is illustrated in Figure.2.1,
where ?=v?/?=v?e?/? indicates the gyroradius, v?=v?(e1 sin?+e2 cos?). The
units vectors e1 and e2 together with e? form an orthogonal system, in which e? =
e1?e2. The set of guiding-center variables (R, ?,?,?, t) is equivalent to (R, ?,v?,?,
25
Figure 2.1: Coordinates transformation from particle phase space to guiding-center
phase space.
t) , so the Vlasov equation can be expressed as
?F
?t +
dR
dt ?
?F
?R +
d?
dt ?
?F
?? +
dv?
dt ?
?F
?v? +
d?
dt ?
?F
?? = 0. (2.7)
In the case that the magnetic ?eld varies slowly in time, the magnetic moment
? is an adiabatic invariant, thus it can be treated as a constant in the guiding-center
phase space. To get the electron guiding-center distribution function Fge, one could
average Fge around the gyro-motion [Lee, 1983; Frieman and Chen, 1982]. Noting
Fge is independent of the gyro-angle ?, the Vlasov equation of guiding-centers can be
simpli?ed as
?Fge
?t +
dRe
dt ?
?Fge
?Re +
dve?
dt ?
?Fge
?ve? = 0. (2.8)
We can rewrite the above equation with the electron generalized momentum pe? as
?Fge
?t +
dRe
dt ?
?Fge
?Re +
dpe?
dt ?
?Fge
?pe? = 0. (2.9)
26
To determine a single electron motion, one should take the gyro-averaging of the
Lorenz force, and get
dRe
dt = ve?b
? + c
qeB0b0 ?[qe < ??
? > +??B0] (2.10)
dpe?
dt = ?b
? ?[qe < ??? > +??B0]. (2.11)
where b? = b0 + (ve?/?e)b0 ? (b0 ? ?)b0, ?? = ?? ? ve ? ?A/c, where b0 is the
local background magnetic direction, ?? and ?A are the perturbed scalar and vector
potentials, respectively [Frieman and Chen, 1982; Lin et al., 2005]. The operator
< .. > represents the gyro-averaging, which is carried out numerically on a discretized
gyro-orbit in real space.
2.1.2 Maxwell Equations of GeFi Model
Once the distribution function is given, physical quantities such as the charged
particle density and current density are known by integrating the distribution function
in the v space. Then combining with Maxwell?s equations, the plasma behavior is
self-consistently determined.
Integrating the ion distribution function fi over the velocity moments, the num-
ber density ni and current density ji are
ni =
?
fid3pi, (2.12)
ji = qi
?
vifid3pi = (qi/mi)
?
pifid3pi ?q2i niA/mic. (2.13)
To proceed with the electron distribution, let us introduce the gyrokinetic order-
ing for electrons
?
?e ?
?e
L ? k??e ?
?B
B ? O(?), (2.14)
27
k??e ? 1, (2.15)
where ? is the mode frequency of interest, L is the macroscopic system length, and
k? is the component of the wave vector in the parallel direction with respect to the
equilibrium magnetic ?eld B0 [Rutherford and Frieman, 1968; Taylor and Hastie,
1968]. The total magnetic ?eld B consists of the equilibrium magnetic ?eld B0 and
the perturbed magnetic ?eld ?B, B = B0+?B. ? is a small parameter. The interested
perpendicular wavelength ranges from the scale of electron Larmor radius to the global
scale length. In our GeFi model, the parallel direction is de?ned as b0 = B0/B0.
The electron gyro-averaging charge density and current in the guiding-center
phase space are
< ne >=
?
Fed3v, (2.16)
< je? >=
?
pe?Fedp, (2.17)
To deduce Poisson?s equation of the GeFi model, we start from the general form
?2??? = ?4?(qini +qene), (2.18)
where we have assumed ?2? >> ?2?. Note ne is the electron density rather than the
guiding-center density, which can be calculated from the guiding-center distribution
by
ne = qem
e
?
d3v(? ?Fge/?w)[??? < ?? > +1c < v? ??A >]+ < Ne >, (2.19)
where w = v2/2 [Lin et al., 2005]. By taking |?e??| < 1, Poisson?s equation becomes
?? ?[(1 + ??
2
pe
??2e )????] +
4?qe?ne
?B ?B? = ?R? +?H?, (2.20)
28
where ?R? and ?H? are expressed as
?R? = ?4?[qi?ni +qe < ?Ne >] (2.21)
?H? = 4?qe[4??ne?B2c j??A??( ?ne?B2)?(?A? ?B)], (2.22)
and ?B is the background magnetic ?eld average over the spatial nd temporal scales of
wave perturbations. Compared with equation (6) in Lin et al. [2005], the appearance
of ?H? is due to the usually higher order terms associated with the inhomogeneities
in the background magnetic ?eld and density.
To solve Poisson?s equation, ?B? must be known. In the low frequency limit
? << ?e, the electron force keeps balanced in the perpendicular direction. Using this
approximation, ?B? can be obtained from the equation below
?neqe???? = ?SB +??[ 14?(1 +??e) ?B?B?] +?H?B +??P?, (2.23)
where
?SB = ?? < ?Peg > +?[ji ?B/c] (2.24)
?H?B = 14??[(???e)?(?A? ?B)???4?jc ??A]??HB (2.25)
?HB = 14?{?[(B??)B]???B2/2} (2.26)
??P? = ?[??( P0cB??
e
????)] (2.27)
??e = ??e/2 (2.28)
?? = (1 +?
?
e) ?B
4? ?B? +
mc2
4?qe??(?
?
e??). (2.29)
Compared with equation (20) in Lin et al. [2005], the newly added term ?H?B is
again due to the usually higher order contribution from the inhomogeneities in the
background magnetic ?eld and density.
29
In the GeFi model, A0 is determined by the initial condition, and the perturba-
tion ?A is solved at each simulation step. We decompose ?A into three components
?A = A??b + A2e2 + Azez, where e2 = ?b?ez and ?b = ?B/ ?B. Then they obey the
following equations
(?2 ? ?
2
pe
c2 ?
?2pi
c2 +
?b??2?b)A? + (?b??2e2)A2 =
?4?c (?ji?+ < ?je? >), (2.30)
(?2 +e2 ??2e2)A2 + (e2 ??2?b)A? = ?4?c ?j2, (2.31)
?2Az = ?4?c ?jz, (2.32)
where
?j2 = c4?(???B??b)?e2 (2.33)
?jz = c4?(???B??b)?ez. (2.34)
Compared with the original model [Lin et al., 2005], terms such as ?b??2?bA? are added
in the improved model to include the inhomogeneities of the background magnetic
?eld direction.
Di?erent from the original scheme [Lin et al., 2005], the above improved model
includes the e?ects of nonuniformity in the background magnetic ?eld direction b and
density ne. To appreciate the e?ects of the terms associated with the variation of the
background magnetic ?eld direction, in the following we compare the ?rst and the
fourth term in the Eq.2.30, ?2A? and ?b??2?bA?. Considering a the Harris current
sheet with a sheet normal in the z direction and the antiparallel ?eld component in
the x direction, the magnetic ?eld pro?le with a spatial coordinate normalized to the
30
current sheet half-width can be expressed as
Bx = Bx0 tanh2(z) (2.35)
By = By0 (2.36)
Bz = 0. (2.37)
In most cases, we can estimate the scale length of A? as the half current width, and
therefore the the ?rst term can be estimated as ?2A? ? A?. The pro?le of ?b??2?b is
shown in Fig.2.2. The peak values ?b??2?b for the cases with Bx0By0 = 1 and Bx0By0 = 0.5
are located in the center of the current sheet and have values around 1.0 and 0.22
repectively. It is seen that By0 ? Bx0, the background nonuniformity cannot be
neglected.
Figure 2.2: ?b??2?b for Bx0By0 = 1 (black) and Bx0By0 = 0.5 (red). The horizontal line is
the normalized z-axis.
31
2.2 Simulation Algorithms
In the GeFi simulation program, full particle method and ?f method have been
used. Full particle method is a standard approach for most particle in cell codes,
while ?f method is a newly developed tech to reduce the computational variations
[Lee, 1987; Parker and Lee, 1993; Chen and Parker, 2003].
2.2.1 The ?f Method
The ?f method is developed to get rid of the background noise caused by the
?nite number of particles. Since only the perturbed distribution function is evolved
in the ?f scheme, thermal ?uctuations associated with the background particles are
removed. To illustrate the idea of the ?f method, we ?rst linearize the ion and
electron Vlasov equations. For ions, it is
d?fi
dt +
qi
mi(????+vi ?(???A))
?fi0
?t = 0, (2.38)
and for electrons, it is
d?Fge
dt + [ve0??b
? +?ve?b?
0 +
cb0
qeB0 ?(qe < ???
? >)]? ?Fge0
?R
+[?b? ?(??B0) +b?0 ?(qe < ???? >)]? ?Fge0?p
e?
= 0, (2.39)
where?b? = (?ve?/?e0)b0?(b0??)b0?(ve0???e/?2e0)b0?(b0??)b0,?ve? = ?qe?A?/(mec).
The distribution function fe and Fge, the potentials A and ? have been divided into
the equilibrium part with subscript 0 and the perturbed part with the notation ?,
f = f0 +?f (2.40)
? = ?0 +?? (2.41)
32
A = A0 +?A. (2.42)
The linearized Vlasov equation above can be simply written in a more general form
as
d?f
dt = ?
df0
dt . (2.43)
Now de?ne the weight carried by each particle as
w = ?ff . (2.44)
Since dfdt = 0, we have
dw
dt =
d?f
dt
f . (2.45)
Using f = f0 +?f, one can get
1
f =
1
f0(1?
df
f ) =
1
f0(1?w). (2.46)
Note that in the above approximation, only the ?rst order is kept. It is easy to show
that
dw
dt = ?(1?w)
1
f0
df0
dt . (2.47)
Since f0 is usually given in an analytical form, the time evolution of particle weights
can be calculated by Eq.2.47.
2.2.2 Advance Particles
To advance particles in the code, we need to numerically integrate the ordinary
di?erence equations with the following form
dZ
dt = G[Z,t], (2.48)
33
where Z represents the motion quantities, such as R(t),p?(t) and w(t) of particles
in the simulation. There are many methods available for the integration of an ODE,
both explicit and implicit. Explicit schemes are straight forward, but should satisfy
the Courant condition ?max?t < 1 in order to assure the numerical stability, where
?max refers to the highest frequency in the simulation [Hockney and Eastwood, 1981].
Implicit schemes are generally stable, but require large computation time to solve
linear algebra equations.
In the improved GeFi model, two implicit methods have been employed, the
predictor-corrector method and the fourth order Runge-Kutta method [Burden and
Faires, 2005]. In the predictor-corrector method, we need to know not only the value
Z at time n but also the value Z at time n-1, then Z at time n+1 can be calculated
in two steps:
Predictor : Z(n+1) = Z(n?1) + 2G(n)[Z(n)]?t
Corrector : Z(n+1) = Z(n) + 12[G(n+1)[Z(n+1)] +G(n)[Z(n)]]?t. (2.49)
The fourth-order Runge-Kutta method needs the value Z at time n and involves four
steps as
Zn+1 = Zn + 16?t(k1 + 2k2 + 2k3 +k4) (2.50)
k1 = G(tn,Zn) (2.51)
k2 = G(tn + 12?t,Zn + 12?tk1) (2.52)
k3 = G(tn + 12?t,Zn + 12?tk2) (2.53)
k4 = G(tn +?t,Zn +?tk3). (2.54)
Note the predictor-corrector method has the second order accuracy while the Runge-
Kutta method above is of the forth order.
34
2.2.3 Operators < .. > and ?..
Two operators < .. > and ?.. have been seen in the gyrokinetic equations (Eq.2.20,
Eq.2.23, Eqs.2.30-32). < .. > means the average around the particle gyro-motion,
while ?.. stands for the background ?eld averaged over the spatial and temporal scales
ofperturbations. Inthissection, Iamgoingtointroducehowtorealizetheseoperators
numerically.
We use the so-called four points gyro-averaging scheme to realize < .. > numeri-
cally [Lee, 1983; Lee, 1987]. The basic idea is to ?nd four points uniformly distributed
on the gyro-circle of a guiding-center and then take average of ?eld values at these
four points. To ?nd the positions of these four points, one can follow the steps as
indicated by Fig.2.3. We let e1 = b?ez where b is the local magnetic ?eld direction
and ez is the units vector along the z direction. (If e1 = 0, one can just try ex or ey
instead.) The ?rst point r1 is located along the direction e1 at a distance ?e from the
guiding-center. Since the magnetic moment ? is conserved in the gyrokinetic model,
the gyroradius can be easily calculated as
?e =
?
2?me
B /qe. (2.55)
In the opposite direction ofe1, another pointr2 can be located. Now makee3 = b?e1.
Repeating the above procedure, one can ?nd all four points. Then the gyro-averaging
?eld quantities are calculated as
< f >= 14(f(r1) +f(r2) +f(r3) +f(r4). (2.56)
To realize ?.., we should store the data of ?elds before the current time step.
However, how many previous steps should be kept and how to average them is still a
35
Figure 2.3: The four points gyro-averaging scheme is shown where b is the local
magnetic direction and r1 is a point on the gyro-circle.
ongoing issue. For the linear simulation or the cases with ?B << B, we can simply
treat ?.. as the initial value of it.
2.2.4 Field Solver
As is seen, Maxwell?s equations in the GeFi model are complicated partial dif-
ference equations (PDE) and mutually dependent on each other. To solve them,
iterations are done simultaneously at each time step to obtain ??, ?B?, and ?A.
At each iteration, one has to solve the nonlinear Poison?s equation with the
following form
?2f(x,y,z)?g(x,y,z)f(x,y,z) = C(x,y,z) (2.57)
Two methods have been employed to handle the above task. The ?rst one is ?nite
di?erences, which is straight forward, but requires much computational time solving
36
the linear algebra equations, especially for 2D and 3D cases. The second one is related
to our particular boundary conditions. In our simulations, periodic boundaries are
adopted in the x and y directions, while conductor boundaries are used in the z
direction. If g(x,y,z) is only a function of z, i.e., g(x,y,z) = g(z), we are allowed to
take the fourier transform on the x and y variables and have
?2f(kx,ky,z)?g(z)f(kx,ky,z) = C(kx,ky,z). (2.58)
Such 1D PDE is solved quickly by ?nite di?erences. If g(x,y,z) is not as simple as
above, we divide g(x,y,z) into two parts like g(x,y,z) = g0(z)+g1(x,y,z). Note g0(z)
is determined by the initial condition that is independent on the y and z coordinates,
whereas g1(x,y,z) is treated as a perturbation. Rewrite Eq.2.57 as
?2f(x,y,z)?g0(z)f(x,y,z) = C(x,y,z) +g1(x,y,z)f(x,y,z) (2.59)
Eq.(2.59) is solved by iteration with the same method of Eq.(2.58). Since usually
g1 << g0, it converges very fast.
2.3 Benchmark for Uniform Plasmas
In this section, we are going to benchmark the modi?ed GeFi scheme for a one
dimensional uniform system [Lin et al., 2011]. The wave vector k is assumed to be
along x. The background magnetic ?eld is in the x-z plane and is allowed to point in
various directions relative to k. The top plot of Fig.2.4 presents a comparison between
the dispersion relations of ?Bz for the fast magnetosonic/whistler branch obtained
from the particle-in-cell GeFi simulation, shown as open dots, and the corresponding
analytical linear dispersion relation based on the fully kinetic mode, shown as solid
lines. Cases with ?e = ?i = 0.04, mi/me = 1836, and k?/k? = 0.2,0.06, and 0
are plotted. The linear ?uctuations in the simulation are due to random noises, as
37
in usual particle simulations. In the case with k? = 0, the electromagnetic mode
approaches the quasi electrostatic lower hybrid mode, and the frequency ?/?i =
?LH/?i = ??i?e/?i =
?
mi/me = 42.8, where ?LH is the lower-hybrid frequency
and ?i is the ion gyrofrequency of the background plasma.
Figure 2.4: The top plot : a comparison between the dispersion relations of ?Bz
for the fast magnetosonic/whistler branch obtained from the GeFi simulation (open
dots), and the corresponding analytical linear dispersion relation based on the fully
kinetic mode (solid lines). The bottom plot : a comparison between the dispersion
relations of ?By for the shear Alfv?en/kinetic Alfv?en mode branch obtained from the
GeFi simulation and the analytical theory.
The bottom plot of Fig.2.4 shows a comparison between the dispersion relations
of ?By for the shear Alfv?en/kinetic Alfv?en mode branch obtained from the GeFi
simulation and the analytical theory for k?/k? = 0.06. The analytical solution of the
MHD shear Alfv?en mode, ?/?i = k?vA/?i = (k?/k)
?
2/?i(k?i) = 0.42, is also shown
as the dashed line.
38
It is seen that GeFi simulation results are in excellent agreement with the theo-
retical analysis for k? << k?. The benchmark has been performed for cases with ?e
and ?i ranging from O(10?2) to O(1).
2.4 Summary
In summary, our GeFi particle simulation model, in which the electrons are
treated as gyrokinetic particles and ions are treated as fully kinetic particles, has
been improved and modi?ed to allow the the existence of modes with wavelengths
on the same scale of the background nonuniformity. With fast electron gyro-motion
and Langmuir oscillations removed from the dynamics, the GeFi model could be a
powerful candidate to solve physics with realistic mass ratio mi/me is a global system.
The linearized GeFi scheme is benchmarked for uniform plasmas, and the sim-
ulation results using the ?f method show that the model can accurately resolve the
physics ranging from Alfv?en waves to lower-hybrid/whister waves, for k? << k? and
? << ?e.
39
Chapter 3
Simulations of Tearing Mode Instabilities under a Finite Guide Field
In this chapter, we use the improved GeFi model to simulate the tearing mode
instability. The initial conditions are described in the ?rst section. The asymptotic
approximation of the Drake-Lee theory [1977a] and the eigenmode theory are pre-
sented next. Simulation results for a current sheet in a simple slab geometry are
shown in the third section. The current sheet con?guration that we use here has been
studied by various people, and can be used to benchmark the GeFi model for the re-
connection physics. In the last two sections, linear and nonlinear simulation results of
tearing mode instabilities for symmetric the Harris current sheets with various guide
?elds are presented and discussed, repectively.
3.1 Simulation Model
the Harris current sheet is an equilibrium solution of the Maxwell-Vlasov system,
which has been extensively used in the magnetic reconnection research [Harris, 1962].
The equilibrium magnetic ?eld consists of two components, the uniform guide ?eld
By0 and the tanh pro?le of antiparallel component Bx,
B = ?xBx0 tanh(z/L) + ?yBy0, (3.1)
where L is the half-width of the Harris current sheet. The particle density is described
by
n = nHsech2(z/L) +nb0, (3.2)
40
Figure 3.1: The magnetic ?eld (left, and the red solid line on the right) and particle
density (black solid line on the right) of the Harris current sheet with a tanh and sech
pro?le along Z direction, respectively. The magnetic ?eld consists of two components,
Bx0 and By0. Particle density is indicated by the black line, where the dotted line
means the background particle density.
where nb0 is the background density and nH is the peak density in the current sheet.
The peak density is obtained from the total pressure balance
nH(Ti +Te) = 18?B2x0, (3.3)
Here we have assumed that the temperatures for ions and electrons are constant
everywhere. The equilibrium velocity distributions for the background ions and elec-
trons are Maxwellian. The equilibrium velocity distributions for ions and electrons
of the Harris current sheet in the GeFi model are given by [Lin et al., 2005; Wang et
al., 2008;]
fHi = nh0(2?T
i/mi)3/2
e?mi[v2x+(vy?vdi)2+v2z]/2Ti ?e?vdiqiAy(z)/Ti (3.4)
fHe = nHsech
2(z/L)
(2?Te/me)3/2 exp(?
B2xmev2de
B22Te )exp(1?
mevde
Te?e
dBx
dz ?)
41
?exp{? 12T
e
[2?B +me(v? ? vdeBGB )2]}. (3.5)
We focus on the instabilities in the x-z plane, where the tearing instability is included,
while the current driven instabilities such as the kink instability, the lower hybrid
instability are excluded. The normalization follows: the magnetic ?eld is in the units
of B0 =
?
B2x0 +B2y0, the spatial length is in the units of ?e, the time step is in
the units of ??1e , particle density is in the units of n0, where n0 is calculated from
?0 = n0Te/(B20/8?) and ?0 is an input parameter.
To benchmark our GeFi model, we also do the calculations for a di?erent current
sheet geometry that has been used in the previous tearing mode studies [Drake and
Lee, 1977; Katanuma and Kamimura, 1980; Wang et al., 2005;]. In this current sheet,
the current pro?le along a strong guide ?eld By0 is given
Jy0(z) = ?eneu0(z) = C1enee?z
2
L2 . (3.6)
In addition to the guide ?eld By0, the current Jy0 generates an antiparallel ?eld
Bx0(z) = 12eneC1?0L??Erf(zL), (3.7)
where the Erf(z) is the error a function. In the following, I refer to this model as
?the Drake-Lee current sheet?. Di?erent from the Harris current sheet, the Drake-Lee
current sheet is not an equilibrium con?guration, and ions are initialized without the
drift motion. This current sheet is particularly suitable for the laboratory plasmas,
in which there is usually a strong guide ?eld.
42
3.2 Eigenmode Analysis
The energy exchange between tearing modes and charged particles takes place
e?ectively in the region where the condition of the Cherenkov resonance is satis?ed:
? = k?vthe. (3.8)
Since the tearing mode is a low frequency phenomenon, the condition Eq.(3.8) is
mainly satis?ed in the region where k? ? 0. If the gyroradius of charged particle
is much smaller than the current sheet width, the drift kinetic approximation is
valid. My benchmark of the tearing instability against the analytical theory will be
conducted for this simpli?ed situation [Hoshino, 1988; Wan et al., 2005].
The linearized drift kinetic equation is given as
??fj
?t +v?
Bx0
B0
??fj
?x +v?
??fj
?y = [?v?
?Bz
B0
1
n(z)
dn(z)
dz ?
dv?
dt
2(vy ?Vdj)
v2thj ]f0j, (3.9)
where j stands for species of particles and uj means the drift speed of ?uid j [Hoshino,
1988; Wan et al., 2005]. Note that in the above equation the electrostatic term has
been neglected, since its e?ect has proven to be small for the tearing mode. We can
rewrite ?Bz and dv?dt in terms of ?Ay:
?Bz = ??Ay?x , (3.10)
dv?
dt =
qj
mj
??Ay
?t , (3.11)
Writing ?Ay=?Ay(z)ei(kxx??t) and using
?
?x = ikx,
?
?y = 0,
?
?t = ?i?. (3.12)
43
Eq.(3.9) becomes
?fj = 1
? ? v?Bx0kxB0
?( v?kxB
0nj(z)
dnj(z)
dz +
?
mj
(vy ?Vdj)
v2thj )foj?Ay. (3.13)
?Ay is induced by the perturbed current ?Jz through Ampere?s law,
??Ay = ??0?Jz = ??0 ?
j=i,e
qj
?
dvyvy?fj (3.14)
Combining the above two equations, we have a second order ordinary di?erential
equation
[ ddz?(z) ddz ?k2x?(z)? ?
j=i,e
Vj(z)]?Ay = 0, (3.15)
where
?(z) = 1 + ?
j=1,e
?2PjV 2dj
c2?2j , (3.16)
Vj(z) = ?2?
2
Pj
c2 (?
2[1 +?z(??)] +? Vdj
vthj(1 +
kxBx0Vdj
?B0 )), (3.17)
? = ?B0k
xBx0vthj
, (3.18)
?? = ?(1? kxBx0Vdj?B
0
), (3.19)
and Z(??) is the plasma dispersion a function. Eq.(3.15) is solved for eigenfunctions
by following the usual shooting method, with the boundary condition ?Ay = 0 at the
z-boundaries on both ends. A factor of 2 di?erence is noticed between the expression
V5j in Hoshino?s equation(20) and the corresponding expression Eq.(3.15) for Vj in our
eigenmode analysis. Following the same procedure, Wan et al. [2005] has conducted
the eigenmode analysis for the Drake-Lee current sheet. One should keep in mind
that the drift kinetic approach is valid only when k??i << 1.
44
3.3 Benchmark for a Simple Current sheet
Figure 3.2: Contour of Ay and eigenfunction of ?Ay with Lz = 10.0?i, Lx = 10.0?i,
L = 0.5?i, By0/Bx0 = 10, Ti = Te, and C1 = ?0.14vthe.
Fig.3.2 shows the contour plot of Ay and the corresponding eigenfunction of ?Ay
obtained from the nonlinear simulation with Lz = 10.0?i, Lx = 10.0?i, L = 0.5?i,
By0/Bx0 = 10, Ti = Te, and C1 = ?0.14vthe. A magnetic island is generated at the
center with an almost symmetric structure, representing a signature of the tearing
mode. ?Ay has a double-peak near the center, which is consistent with the tearing
mode unstable condition ?? > 0. Recall the de?nition of ??
?? = (??Ay(?)?z ? ??Ay(??)?z )/?Ay(0) (3.20)
where ? is the width of the singular layer.
The linear growth rate as a function of Bx0/By0 is shown in Fig.3.3 with the
same parameter above except for a varying Bx0/By0. ?Fixed ions? means ions are
only advanced by the initial ?elds, meaning the motion of ions is kept unperturbed.
45
Figure 3.3: Relationship between the linear growth rate and Bx0/By0 with Lz =
10.0?i, Lx = 10.0?i, L = 0.5?i, Ti = Te, and C1 = ?0.14vthe
46
Again the linear growth rate agrees the eigenmode theory well. Two other conclusions
can be reached immediately: the growth rate is increasing with Bx0/By0, and the ions
contribute little to the growing of tearing modes. According to the linear theory, the
linear growth rate is given by
? = ?
?k
xvthe
lsd2e , (3.21)
where ls is the magnetic ?eld shear length, vthe is the electron thermal velocity. Since
ls = By0?Bx0/?z ? By0/Bx0, we have ? ? Bx0/By0. Since the current is carried totally
by electrons for the Drake-Lee current, ions contribute little to the perturbed current
based on ?Jy = ?j=i,e(nj0?vjy + vj0y?nj0), where the ratio of the ?rst term between
ions and electrons is ni0?viy/ne0?vey = me/mi, and the second term of ions is 0.
3.4 Linear Simulation of the Tearing Mode in the Harris current Sheet
the Harris current sheet is the ?rst equilibrium model of one-dimensional current
sheet, which has been widely used in the investigation of reconnection. Note that
since the Drake-Lee current sheet is not in equilibrium, it can only be applied to
cases with By0/Bx0 >> 1. In this section, I show the linear simulation results of
tearing modes in the Harris current sheet.
In order to benchmark the GeFi simulation scheme against the existing eigen-
mode analysis, GeFi simulations with ? = 0, k??i ? 0, ky = 0 and ?B? = 0 are
carried out. These parameters are consistent with the assumptions of eigenmode
analysis, where ions and electrons are both treated as drift kinetic particles requiring
k??i << 1. Since the electrostatic, compressible modes are neglected, we have ? = 0,
?B? = 0. Fig.3.4 shows the linear growth rate as a function of kxL (top) and L/?s
(bottom) with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 (top), and kxL = 0.41
(bottom), where ?s is de?ned as ?s =
?
mi/me?e which is the ion gyroradius when
Ti/Te = 1. Again the simulation results agree with the eigenmode analysis well. In
47
the top ?gure, we can see the growth rate is maximized around kxL = 0.4 and tends
to be 0 for kxL > 1, both of which represent the characteristics of tearing modes.
The bottom ?gure depicts that the growth rate is decreasing with the current width
with a relationship ? ? L?2.4.
Figure 3.4: The linear growth rate as a function of kxL (top) and L/?s (bottom) with
Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 (top), and kxL = 0.41 (bottom).
The contours of the magnetic ?eld perturbations Bx, By, Bz from the eigenmode
theory and the GeFi simulation are shown in the left Fig.3.5 with Bx0/By0 = 0.075,
Ti/Te = 0.01, L/?s = 0.75, kxL = 0.5, while the corresponding eigenfunctions are
shown on the right. Little di?erence exists between the simulations and the eigenmode
analysis. The Bx, By perturbations are symmetric, while Bz shifts in the opposite
directions on each side of the current sheet.
As discussed in the chapter2, the nonuniform terms associated with the magnetic
?eld direction b and the electron density ne have been added into the improved GeFi
48
Figure 3.5: The contours of perturbations magnetic ?eld and the corresponding eigen-
functions with Bx0/By0 = 0.075, Ti/Te = 0.01, L/?s = 0.75 kxL = 0.5.
49
Figure 3.6: The e?ect of inhomogeneity with cases Bx0/By0 = 0.05 (left column)
and Bx0/By0 = 1 (right column), for Ti/Te = 0.01 and L = 0.5?s. The top row :
eigenfunctions ?Ay (solid lines from GeFi simulation, dashed lines from eigenmode
calculation; the middle row : the contour of ?Ay from simulations; the bottom row :
the contour of ?Ay with ?H?B turned o?.
model. As is seen in Fig.2.2, such an e?ect becomes more profound as By0 gets
smaller. The top row of Fig.3.6 shows the absolute values of eigenfunctions of ?Ay
obtained from simulations with Bx0/By0 = 0.05 and Bx0/By0 = 1, for Ti/Te = 0.01
and L = 0.5?s. The solid line in each plot shows the simulation results based on the
50
improved GeFi scheme, while the dashed line corresponding to the run in which ?H?B
is turned o?. The corresponding contours of the eigenfunctions of ?Ay obtained from
the GeFi simulation are shown in the middle row of Fig.3.6, while the runs with ?H?B
turned o? are plotted in the bottom row. Clearly, while the di?erence between the
results with and without the ?H?B is negligible for Bx0/By0 = 0.05, the corresponding
eigenmode structures of ?Ay are, however, signi?cantly di?erent in the case with
Bx0/By0 = 1.0. With ?H?B terms included, the tearing mode shows ?ner structure in
the current sheet. The growth rates of tearing modes for cases Bx0/By0 = 1.0 with
and without the ?H?B term are, respectively, 0.098 and 0.128. Further simulation
indicates that the ?H?B term in general cannot be ignored for linear tearing modes
when By0 ? or < Bx0.
Figure 3.7: The growth rate as a function of Ti/Te for three cases L = ?s,1.5?s,3.0?s
with Bx0/By0 = 0.05, kxL = 0.5.
51
The above benchmark work has been limited to the extremely small tempera-
ture ratio between ions and electrons Ti/Te = 0.01 in order to satisfy the drift kinetic
ion approximation requirement k??i << 1. Since k? ? L?1 (L ? ?S for our appli-
cations) and ?i = ?s
?
Ti/Te, the parameters k??i are evaluated as ? 0.1 and ? 1,
respectively, for cases with Ti/Te = 0.01 and for cases with the more realistic tem-
perature ratio Ti/Te = 1.0. It is expected that as the ratio Ti/Te increases, the drift
kinetic ion approximation is not valid any more, and therefore the agreement between
the eigenmode theory and the GeFi simulation cannot hold. The growth rate as a
function of Ti/Te is plotted in Fig.3.7, where three cases with L = ?s,1.5?s,3.0?s
for Bx0/By0 = 0.05 and kxL = 0.5 are presented. The di?erence between the eigen-
mode theory and the GeFi simulation becomes obvious after the value Ti/Te ex-
ceeds 0.5,1.2,3.0, corresponding to k??i ? 0.7,0.73,0.67, repectively, for cases with
L = ?s,1.5?s,3.0?s. Our simulation results here are qualitatively consistent with the
previous study by Ricci et al. [2005], in which it is pointed out both the ?nite electron
Larmor radius e?ect and the ion Larmor radius e?ect are potentially signi?cant in
the calculation of the growth rate of tearing modes.
To explore the tearing mode properties, the realistic temperature ratio Ti/Te ? 1
should be adopted in the simulations. Fig.3.8 displays the linear growth rate as a
function of kxL with Ti/Te = 1.0, L = 0.5?s for three di?erent guide ?elds Bx0/By0 =
5.0, 3.0, 1.0. The linear growth rate is rapidly reduced as the guide ?eld deceases and
kxL corresponding to the peak growth rate shifts to the short wavelength a bit. It
seems not to make sense that the linear growth rate is bigger for larger guide ?elds.
Nevertheless, we ?nd such di?erence comes from the way of our normalization. Based
on the Drake-Lee kinetic theory, Karimabadi et al. [2005] deduced an asymptotic
formulation of the growth rate of tearing modes in the Harris current sheet as
?/?i = 2?2mim
e
(1 + TiT
e
)1?k
2
xL
2
?? Bx0B
y0
(?eL)3. (3.22)
52
In the above equation the magnetic ?eld is normalized to the antiparallel magnetic
?eld Bx0, di?erent from our normalization to the asymptotic total magnetic ?eld B0.
Transforming it to our units, we have
?/?i = 2?2mim
e
(1 + TiT
e
)1?k
2
xL
2
?? (Bx0B
y0
+ By0B
x0
)(?eL)3. (3.23)
? ? (Bx0By0 + By0Bx0), which is increasing with By0Bx0 when By0Bx0 > 1.
Figure 3.8: The growth rate as a function of kxL for three di?erent guiding ?elds
Bx0/By0 = 5.0, 3.0, 1.0 with Ti/Te = 1.0, L = 0.5?s.
53
The ion kinetic e?ect in the tearing mode is examined in Fig.3.9 with Ti/Te = 1.0,
Bx0/By0 = 5.0, kxL = 0.41 for di?erent current widthes. ? represents cases of ?with
ions?, while ? represents cases of ?without ions? that has the same meaning of
??xed ions? we discussed above, and the line stands for the ratio of growth rate
between them. Two conclusions can be reached immediately. One is that the growth
rate decreases rapidly as the current sheet gets wider. The other is that ions play a
signi?cant role in the tearing mode, raising the growth rate almost by 50%. Daughton
et al. [2005] has demonstrated ions can directly resonant with tearing modes and
contribute to around 30% in?uence in the growth rate of tearing modes for thin
current sheets. The compressibility of ions can also contribute to the tearing mode.
Our simulations show that this term vdi?ni contributes to 10% of the total perturbed
current ?j. Finite ion Larmor radius e?ect is another potentially important factor
for cases with k??i ? 1. One may have noticed that in the previous section our
simulations have shown that the ion e?ect is negligible in the Drake-Lee current
sheet. Such inconsistency is due to that in the Drake-Lee current sheet the current
is totally carried by electrons, while in the Harris current sheet the ratio between the
ion current and the electron current is determined by ji/je = Ti/Te.
Because of the low frequency property of tearing modes, the electrostatic e?ect
was not considered in the kinetic theory by Drake and Lee [1977a]. Nevertheless,
Hoshino [1988] pointed out that the electrostatic ?eld in?uences the induced current,
and modi?es the linear growth rate. His analysis shows that in the presence of a guide
?eld, the electrostatic ?eld generated by the compressibility of electrons will increase
the growth rate, but the ion e?ect is ignored due to its large Larmor radius. However,
previous studies and our simulations have shown that the ?nite ion Larmor radius
e?ect and the compressibility of ions are potentially important in the calculation of
the growth rate, so it is necessary to investigate the electrostatic e?ect using the GeFi
simulation. Fig.3.10 displays the growth rate as a function of Ti/Te with L = 0.5?s,
54
Figure 3.9: The growth rate as a function of L with Ti/Te = 1.0, Bx0/By0 = 5.0,
kxL = 0.41. ? : cases with ions; ? : cases without ions; solid lines: the growth rate
ratio between cases with ions and cases without ions.
Bx0/By0 = 3.0, kxL = 0.41. ? represents cases of ?with electrostatic e?ect?, ?
represents cases of ?without electrostatic e?ect?, and the line stands for the growth
rate ratio between them. We see that the growth rate is increasing with Ti/Te and
the in?uence of electrostatic e?ect depends on the value of Ti/Te. Electrostatic e?ect
reduces the growth rate when Ti/Te < 4, but raises the growth rate when Ti/Te >
4. The ion Larmor radius becomes larger as Ti/Te increases, and then Hoshino?s
assumption is satis?ed.
In Fig.3.11, the contours of typical perturbation structures are plotted with L =
0.5?s, Bx0/By0 = 3.0, kxL = 0.41, Ti/Te = 2.0, where Bx, By, Bz, ni,, ne, ?, A?,
Uix, Uiy, Uiz, Uex, Uey, Uez, and je? are repectively three components of the magnetic
?elds, the ion density, the electron density, the electrostatic potential, the parallel
component of the vector potential, three components of the ion ?uid velocity, three
components of the electron ?uid velocity, and the parallel component of the electron
55
Figure 3.10: The growth rate as a function of Ti/Te with L = 0.5?s, Bx0/By0 = 3.0,
kxL = 0.41. ? : with electrostatic e?ect; ? : without electrostatic e?ect; the solid
line : the growth rate ratio between them.
current. An important quantity in the study of the tearing instability is the singular
layer width?, which is de?ned as the region where electrons can free stream along the
direction of the perturbed electric ?eld and Landau resonate with the wave. From the
contour of Uey, we can see that the singular layer?extends to the current sheet width
L with the scale of ?i, and therefore the Drake and Lee kinetic theory assumption
? << L is clearly violated for our parameters. Considering the current parameters
?e = 0.01 and ?i = 0.02, we have de/?e = 1/??e = 10, and di/?i = 1/??e = 7.
Since the characteristic thickness of the current sheet is comparable to de and ?i, two
factors can contribute to broaden the singular layer: one is the ?nite electron skin
depth e?ect, and the other is the ?nite ion Larmor radius e?ect.
The linear stage of collisionless tearing modes contains important information of
the out-of-plane perturbation By. In the limit of a weak guide ?eld, the perturbation
has a quadrupole structure, which has been identi?ed as a key signature of the Hall
56
Figure 3.11: The contours of typical perturbation structure with L = 0.5?s,
Bx0/By0 = 3.0, kxL = 0.41, Ti/Te = 2.0.
57
e?ect in fully nonlinear fast reconnection [Sonnerup, 1974; Terasawa, 1983; Karam-
abadi et al., 2004;]. Simulations by Daughton et al. [2005] show that the addition of a
guide ?eld complicates the out-of-plane perturbation and compresses the spatial width
from the ion inertial scale down to the electron kinetic scale. Our simulation shows
that the perturbation of By in Fig.3.11 is not of a quadrupole structure but is located
in the center of the current sheet with a certain asymmetry. Since di >> ?i ? L for
our parameters, the Hall e?ect is excluded form the formation of the perturbation of
By. Examinations on the contours of Uix and Uex indicate that the ion and electron
motions follow di?erent spatial scales and thus cause the perturbation of By. Since
electrons are strongly magnetized, the electron motion is mainly guided along the
magnetic ?eld. Due to the di?erent directions of magnetic ?eld below and above the
center of current sheet, Uex, the projection of the electron ?uid velocity on the x axis,
shows a quadrupole structure and has a peak value around the O-type neutral point.
Since the ion Larmor radius is comparable to the current sheet width, ions are able
to move across the ?eld lines and be accelerated directly by magnetic tension forces.
The structure of Uix is symmetric with respect to the X-type neutral point with a
?/2 phase di?erence compared with Uex. The separation between Uex and Uix forms
an in-plane current, and consequently generates the out-of-plane perturbation By.
Nonlinear simulations by Ricci et al. [2005] show the asymmetry of Uex results in the
compression of electron ?owing along the parallel direction, and leads to a density
asymmetry across the dissipation region. Our conclusions here are consistent with
the discussion by Drake and Shay [edited by Birn and Priest, 2007].
The comparison among the growth rate of tearing modes as a function of Bx0/By0
obtained from the GeFi simulation (open circles), the eigenmode calculation (solid
line) and the Drake and Lee analytical asymptotic matching [1977a] Eq.3.23 (dashed
line) is plotted in Fig.3.12. Results based on the asymptotic matching Eq.3.23 devi-
ate signi?cantly from the eigenmode calculation and the GeFi simulation, although
58
the general trend of ? versus By0 is consistent with them. This deviation can be
understood, since in the cases shown in Fig.3.12 the thickness of current sheet L is
comparable to de at the sheet center and thus the conditions assumed in the asymp-
totic matching analysis break down. In addition, the ion e?ects are ignored in the
matching theory.
Figure 3.12: The comparison among tearing mode growth rate as a function of
Bx0/By0 obtained from GeFi simulation (open circles), the eigenmode calculation
(solid line) and Drake and Lee analytical asymptotic matching Eq.3.23 (dashed line).
3.5 Nonlinear Simulation of Tearing Modes in the Harris current Sheet
In this section, we perform nonlinear simulations to study the saturation of tear-
ing modes. It has been suggested that when the magnetic island reaches the spatial
order of the singular layer width, the electron orbits are altered. This results in the
saturation of tearing modes. Two cases, one without ions and the other with ions, are
shown ?rst. In these two cases, only parallel component of the vector potential (i.e.,
A? is dropped due to weak magnetic compression) is kept in the simulation model,
59
which is consistent with the previous theoretical analysis. In the end, the coalescence
of multiple tearing modes is studied.
Fig.3.13 shows the time evolution of the magnitude of eigenfunction A? (left)
and the contour plot of Ay at time t = 20.0 (right). The parameters here are L =
0.5?s, Bx0/By0 = 10.0, Ti/Te = 1.0, and kxL = 0.41. This implies that only the
corresponding fastest growing mode is kept in the system. The tearing mode is in
the linear growing stage before t = 9??1i with the growth rate ? = 0.83, which is
almost the same as ? = 0.812 obtained from the linear simulation. After t = 9??1i ,
the tearing mode saturates and oscillates with a frequency ? = 2?/8 = 0.785?i. Wan
et al. [2005] deduced a relationship between the growth rate ? and the oscillation
frequency ? as ?/? = 1.22. Clearly our result ?/? = 0.94 is di?erent from their
prediction, for their analysis still relies on de << L. At t = 20??1i , the width of the
saturation island is w = 0.44?s, smaller than 2L = 1.0. Considering ? = 0.005, we
have de/?e = 14.14 and de/?i = 0.33, so the saturation width of the magnetic island
is on the order of de. Based on the Drake and Lee theory, the singular layer thickness
in the Harris current sheet can be approximated as [Karamabadi et al., 2005]
? = 1?? 1?(kxL)
2
kxL
?2e
L
B0
Bx0(1 +Ti/Te). (3.24)
For our parameters, the above equation gives ? = 1.1?e, which is much smaller
than the saturation island width w obtained from the simulation. Clearly the ?nite
electron skin length e?ect that is neglected in the Drake and Lee theory can broaden
the singular layer and yield a larger saturation amplitude.
The results of simulation with the same parameters above but including the ion
motion are depicted in Fig.3.14. It shows the time evolutions of the magnitude of
eigenfunctions ne (top left), ni (top right), A? (bottom left), and the contour plot
of Ay at time t = 20??1i (bottom right). Two main di?erences can be observed
60
Figure 3.13: The time evolution of the magnitude of eigenfunction A? (left) and
the contour plot of Ay at time t = 20??1i (right) with L = 0.5?s, Bx0/By0 = 10.0,
Ti/Te = 1.0, and kxL = 0.41.
by comparing with the case without ions: one is the faster oscillation frequency
? = 1.1?i; the other is the wider saturation island w = 1.05?i ? 2L. By including the
ion in?uence that was neglected in the work of Drake and Lee, it has been suggested
that much larger saturation amplitudes comparable to the ion gyroradius could be
obtained [Galeev, 1988; Kuznetstova and Zelenyi, 1990]. The strong correlations
between the time evolution of ?A? and ni in Fig.3.14 proves that ions play a signi?cant
role in the saturation mechanism.
The nonlinear evolution of multiple tearing modes is shown in Fig.3.15. Four
unstable modes are kept in the domain, with the wave number kxL = 0.11 (mode1),
kxL = 0.22 (mode2), kxL = 0.33 (mode3), and kxL = 0.44 (mode4) respectively. Note
that the wavelength of mode1 is the same as the system length in the x direction.
Since mode4 corresponds to the wavelength with roughly the largest growth rate, four
magnetic islands are formed at the initial stage. After t = 8.72??1i , the islands start
61
Figure 3.14: The time evolutions of the magnitude of eigenfunction ne (top left), ni
(top right), A? (bottom left), and the contour plot of Ay at time t = 20??1i (bottom
right) including ions.
62
to coalesce due to the attraction between currents, and two wider islands become
dominated at t = 13.08??1i . Two much smaller plasmoids become evident at t =
15.26??1i . The formation of such secondary island may be due to the thin current
near the x-point [Daughton et al., 2005; Karamababidi et al., 2005]. The coalescence
of two big islands continues to compress the smaller plasmoid and triggers a forced
reconnection. At t = 17.44??1i , a very narrow reconnection layer can be seen between
two large islands. At t = 19.62??1i , the system relaxes back to the state that two
islands are dominated.
Fig.3.16 shows the time evolution of the magnitude of eigenfunction A? for the
four modes. Initially only the magnitude of mode4 corresponding to the fastest grow-
ing tearing mode is increasing. Around t = 9??1i , the coalescence of the islands
transfer the energy from mode4 to mode2. At t = 14??1i , all modes except for mode2
start to increase rapidly. This is due to the generation of small islands near the X-lines
as shown in Fig.3.16 at t = 15.26??1i . The magnitude of all modes are deceasing dra-
matically after t = 17??1i , because the energy is released by the forced reconnection
as discussed above.
3.6 Summary
In this chapter, the GeFi simulation model with the realistic electron/ion mass
ratio is employed to study tearing mode instabilities under a ?nite guide ?eld. The
results are summarized below:
1. The GeFi simulation model is benchmarked with tearing mode instabilities
in a simple current sheet. Our results agree with the eigenmode analysis and the
previous studies.
2. With the drift kinetic approximation, the eigenmode theory for the Harris cur-
rent sheet is deduced. For cases with Ti << Te when the drift kinetic approximation
is valid, the GeFi simulation results agree with the eigenmode analysis.
63
Figure 3.15: Contours of vector potential Ay at di?erent times for the case with four
unstable modes.
64
Figure 3.16: The time evolution of the magnitude A? for the four modes.
65
3. The inhomogeneity terms a?ect the structure of tearing modes, and cannot
be neglected when By0 ? Bx0.
4. For the temperature ratio Ti/Te ? 1, the simulation studies reveal the re-
lationship between the growth rate of tearing modes and kxL, Bx0/By0, Ti/Te, and
L/?s.
5. The electrostatic e?ect is negligible on the calculation of the growth rate of
tearing modes.
6. In the general cases in which the ion gyroradius is comparable to the current
sheet width, the ion e?ects are found to play a signi?cant role on the calculation of
the growth rate of tearing modes. In such cases, the Drake and Lee theory [1977a]
and the eignemode analysis based on the drift kinetic assumption are not invalid.
7. The nonlinear GeFi simulation shows that the ion e?ect can broaden the
width of saturated magnetic islands from the electron dynamic spatial scale to the
ion gyroradius.
8. The case with multiple tearing modes shows that the coalescence of magnetic
islands can form a larger saturation island than the case of a single tearing mode.
During the process of coalescence, small magnetic islands are formed around the X
point. The small-scale magnetic reconnection can occur when the small magnetic
islands are compressed by the large magnetic islands. This process results in fast
energy release.
66
Chapter 4
Generation of Low Frequency Electromagnetic Waves by Magnetic Reconnection
under a Finite Guide Field
In the collisionless plasma, magnetic reconnection is usually triggered by mi-
cro instabilities. The change of topology of the magnetic ?eld lines and and the
perturbations in the di?usion region will lead to the generation of low frequency elec-
tromagnetic waves. Previous 2D simulations of steady-state reconnection have found
these waves to be present in the form of MHD discontinuities, as also observed in the
magnetosphere as well as the interplanetary space [Sonnerup et al., 1981; Phan et al.,
2006]. In this Chapter, we investigate the 3D structure of the reconnection layer. The
parameters chosen will be suitable for the Earth?s magnetosphere. It is found that
kinetic Alfv?en waves are present in the 3D reconnection. First the simulation model
is introduced. Then a brief review of the kinetic Alfv?en wave is given. Our results
show the kinetic e?ects play an important role in the 3D reconnection geometry. The
simulation results will be presented in the third section.
4.1 3D Hybrid Simulation Model
In the hybrid code, ions are treated as discrete particles and electrons are treated
as a massless ?uid. The code was developed by Swift[1995], which has been used in
a general curvilinear coordinate system to model the Earth?s magnetosphere. In this
study we use a Cartesian coordinate system to model the generation of low frequency
electromagnetic waves.
67
In the hybrid code, quasi charge neutrality is assumed. The equation for the ion
motion is given by
mi
qi
dvi
dt = E+vi ?B??(Vi ?Ve) (4.1)
where vi is the ion particle velocity, E is the electric ?eld in units of ion acceleration,
B is the magnetic ?eld in units of the ion gyrofrequency, ? is the collision frequency
which is used to model the resistivity at the X-line, and vi and Ve are the bulk of
?ow velocities of ions and electrons, repectively. The electron momentum equation is
written in the form
E = ?Ve ?B??(Ve ?Vi) +?Pe, (4.2)
where Pe is the electron pressure. The electron ?ow speed is evaluated from Ampere?s
law,
Ve = Vi ? ??B?N , (4.3)
where ? = 4?e2/mic2, e is the electron charge, mi is the ion mass, and N is the
particle number density. The magnetic ?eld is advanced in time using Faraday?s law
?B
?t = ???E, (4.4)
where the electric ?eld is calculated from Eq.4.2.
In the simulation the ion particle velocity is updated at half time steps with a
second-order accuracy. The magnetic ?eld and the particle positions are advanced at
the integer time step. The time step is chosen to satisfy the Courant condition with
respect to the whistler mode. In the simulation, the magnetic ?eld is advanced 10
time steps for every time step the ion velocities are advanced.
Here, the magnetosphere plasmas are considered, in particular for the symmetric
reconnection as usually occurring in the magnetostail. The simulation domain is a
cube with x being normal to the current sheet and z being the antiparallel magnetic
68
?eld direction. Di?erent with the previous chapters, we choose the z direction as the
antiparallel direction, which is generally used in the magnetospause research. The
initial pro?le of the z component magnetic ?eld is given by
Bz(x) = Bz0(x)tanh(x/?) (4.5)
where Bz0 is the z-component magnetic ?eld in the lobes and ? is the half-width of
initial current sheet. The pro?le of the y-component magnetic ?eld, i.e., the guide
?eld in the lobes, is assumed to be
By = By0. (4.6)
The initial x-component magnetic ?eld is Bx = 0. The pro?le of ion thermal pressure
is determined from the total pressure balance across the current sheet. A constant,
isotropic temperature with Ti = Te is assumed, where Ti is the ion temperature and
Te the electron temperature. The ion number density is thus given by
N(x) = N0(1 + 1?
0
[1?B2z0/B20tanh2(x/?)?B2y0/B20]) (4.7)
where B0 = (B2z0 +B2y0)12 is the total magnetic ?eld in each lobe.
We study the evolution of a spontaneous reconnection [Scholer, 1989]. The dy-
namics are not driven. A constant resistivity in time is imposed in the system to
trigger the reconnection. The resistivity is modeled through a collisional term in the
ion equation of motion and the electron momentum equation. The collision frequency
is assumed to be
? = ?0exp(?(x2 +z2)/?20)exp(?y2/?20) (4.8)
where ?0 ? 0.1?2?i and ?i is the ion gyrofrequency in the lobes. We de?ne ?0 as
the half-length of X-line.
69
In our simulation, the time is in units of ??1i . For the results shown in this thesis
the magnetic ?eld is expressed in units of B0, the ion number density in units of N0,
and the temperature in units of the lobe temperature T0, where the subscript ?0?
represents the quantities in the lobes. The velocity is normalized to the lobe Alfv?en
speed VA0, and the spatial coordinate is normalized to di the lobe ion inertial length.
The conductor boundary is used in the x direction, , while the periodic bound-
aries are used in both the z and the y directions. To limit the boundary e?ect, the
simulation stops when the waves propagate near the boundaries. The domain size
X ?Y ?Z is 128di ?128di ?256di with the grid size ?x = 0.5di, ?y = 2.0di, and
?z = 1.0di. Almost 100 particles per cell are placed in the lobe region.
4.2 Properties of Kinetic Alfv en Waves
The Alfv?en wave is one of the most important wave modes in magnetized plas-
mas. It exists due to the balance between the magnetic ?eld tension and ion inertia.
These waves were ?rst theoretically predicted by H.Alfv?en and are now called shear
Alfv?en waves because the perturbations of the magnetic ?eld are perpendicular to the
ambient ?eld and wave vector. When the perpendicular wavelengths are comparable
to the ion gyroradius ?i, the properties of shear Alfv?en waves need to be modi?ed.
In a plasma with intermediate beta values (me/mi < ?e << 1)), the modi?ed waves
are usually refereed as kinetic Alfv?en waves.
Starting from the two ?uid-MHD equations, one can deduce the dispersion rela-
tionship of kinetic Alfv?en waves as [Stix, 1992]
?2 = k2?V 2A(1 + (1 + TiT
e
)k2??2i). (4.9)
Obviously the kinetic Alfv?en wave is dispersive, while the shear Alfv?en wave is non-
dispersive.
70
Another important di?erence between the shear Alfv?en wave and the kinetic
Alfv?en wave is that the kinetic Alfv?en wave is compressible and generates an electric
?eld parallel to the ambient magnetic ?eld by
E? = ? Teen
e
??ne. (4.10)
It is suggested that such electric ?eld plays an important role in the energy transport
processes and particle heating/acceleration of space plasmas.
4.3 Simulation Results
Two typical cases with By0 = 0.5Bz0 are presented in this section. One possesses
an in?nitely long X-line, and the other possesses with a ?nite X-line. The cases with
an in?nite X-line corresponds to ky ? 0. Since the propagation of waves is almost
limited in the x-z plane, results similar to the 2D simulations are expected although
some di?erences are found to exist. For the case with a ?nite X-line, the spatial e?ects
along the y direction are found to lead to ky ?= 0, which results in the propagation of
waves very di?erent from those in 2D simulations.
4.3.1 Case with an In nite X-line: 2?0 = ?
Fig.4.1 shows the contour plots of the magnetic ?eld Bx, By, Bz, the ion ?uid
velocity Vx, Vy, Vz, the temperature T, the particle number density N, and the
con?guration of magnetic ?eld lines in the x-z plane located aty = 0 whent = 150??1i .
A plasma bulge and a trailing quasi-steady reconnection layer are formed in the
out?ow region of magnetic reconnection, as shown in the magnetic ?eld and plasma
plots in Fig.4.1. The magnetic ?eld decreases and the plasma temperature increases
from the lobe to the center of the reconnection layer. Since ky ? 0, the ?eld and
71
plasma properties in the x-z planes of other y positions resemble the results shown
here.
Figure 4.1: The contour plots of Bx, By, Bz, Vx, Vy, Vz, T, N, and con?guration of
magnetic ?eld lines. The solid line corresponds to the position at z = 40di, along
which there appears a quasi-steady structure. The solid line at z = 40di indicates the
boundary of the leading bugle. The solid line z = 40di shows the transit boundary
between the plasma bugle and the steady state layer.
72
Figure 4.2: The spatial pro?les of Bx, By, Bz, N, Vx, Vy, Vz, T?, T?, J?, and B at
z = 40di and hodograms of the tangential magnetic ?eld (Bz-By).
73
Fig.4.2 depicts the spatial pro?les of the magnetic ?eld Bx, By, Bz, density
N, the ion ?ow velocity Vx, Vy, Vz, the parallel temperature T?, the perpendicular
temperature T?, the parallel current density J?, and the magnitude of magnetic ?eld
B at z = 40di, and the hodograms of the tangential magnetic ?eld (Bz-By). Note
that z = 40di is located in the steady reconnection layer as indicated in Fig.4.1.
We can see that two discontinuities bound the steady reconnection layer with the
wave front at x = ?8.5di and the width of ? 2.3di. Across the discontinuities, the
plasma is accelerated so that Vz is changed from ? 0 in the upstream to ? 0.75VA0
in the downstream corresponding to the change of the Alfv?en speed ? 0.7VA0. The
hodogram shows that the tangential ?eld rotates by an angle of around? ?/2 through
each discontinuity. These pair of discontinuities are roughly rotational discontinuities.
Notice that the magnetic ?eld strength slightly decreases and the ion density increases
across the rotational discontinuity. The changes in the magnetic ?eld and ion density
at the rotational discontinuity are due to the increase of the parallel temperature
T? > T? because of the streaming ions along the magnetic ?eld lines. The slow shock,
however, is too weak to be identi?ed in the simulation because of T? > T?. All results
shown here are consistent with the previous 2D simulations.
To study the 3D e?ect, we examine the variations along the y direction. The
variation of a certain quantity f is calculated as f(x,y,z) ? ?f(x,z), where ?f(x,z)
is the average along the y direction at ?xed positions (x,z). The top two ?gures in
Fig.4.3 present the contour plots of variations of Bx and N along the y direction in
the x-z plane at y = 0, and the bottom two ?gures present the contour plots of the
corresponding variations in the y-z plane at x = 4di as indicated by the solid line in
the top ?gures. At time t = 150??1i , one can see that the variations are signi?cant
in the transition region between the leading bulge and the quasi-steady structure. In
the y-z plane at x = 4di the structure of the variation of Bx stretches out along the
initial local magnetic direction b. Since the x direction is perpendicular to the initial
74
magnetic ?eld, the variation of Bx represents the shear waves, while the variation of
N represents the compressible waves.
Figure 4.3: Top ?gures: the contours of variations of Bx and N along the y direction
in the x-z plane at y = 0; bottom ?gures: the contours of variations of Bx and N
along the y direction in the y-z plane at x = 4di. The solid line along z = 52di in the
top left ?gure indicates the boundary of the leading bugle. The solid lines in the top
?gures indicate the y-z plane at x = 4di; the solid lines in the bottom ?gures show
the local initial magnetic ?eld direction.
Fig.4.4 shows the contours of variations along the y direction of Bx, By, Bz, N,
Vx, Vy, Vz, Ex, Ey, Ey, J?, E? in the y-z plane at x = 3di when t = 150??1i . In these
75
structures, the shear components Bx and Vx, Ey, and the parallel current density
J? are well correlated. In Fig.4.3, one may have noticed that the spatial scale of the
variation of Bx is ? 3di or ? 10?i, considering to ?i = 0.1. Since kx?i ? 2?10?i?i ? 0.62,
certain kinetic e?ect are expected, which is justi?ed by the strong correlation between
the parallel electric ?eld E? and Bx.
Figure 4.4: The contours of variations along the y direction of Bx, By, Bz, N, Vx, Vy,
Vz, Ex, Ey, Ey, J?, E? in the y-z plane at x = 3di when t = 150??1i .
76
4.3.2 Case with a Finite X-line: 2?0 = 8di
In reality, the length of X-line cannot be in?nite. Shay et al. [2010] have shown
that 3D structures develop in a spontaneous reconnection, corresponding to spatially
isolated sites of reconnection with scale lengths around 10di. In the following, we
show a case in which the length of the X-line is set as 2?0 = 8di to study the ?nite
length X-line e?ect.
The top two ?gures in Fig.4.5 show the con?gurations of magnetic ?eld lines
near the current sheet at t = 50??1i (left) and t = 80??1i (right). At t = 50??1i ,
some magnetic ?eld lines have been reconnected. Since the length of the X-line is
?nite, the process of reconnection is limited in a ?nite region. At t = 80??1i , the
reconnected ?eld lines are dragged towards ?z directions by the accelerated plasmas
and are twisted at the same time. It is due to that the large-amplitude perturbations
generated by the reconnection and their consequent propagations change the local
topology of magnetic ?eld lines. The contour plots of Bx (left) and Vx (right) in
the x-z plane at y = 25di of t = 80??1i are presented in the bottom two plots.
We can see that the shear perturbations of Bx and Vx are evident in the plane far
from the reconnection region, and they are well correlated. The perturbation regions
are associated with the reconnected ?eld lines, while the ?eld lines outside of the
perturbations stay unreconnected. It indicates that the perturbations caused by the
reconnection propagate into the sites far from the di?usion region along the magnetic
?eld lines. We also notice that the perturbations in the x-z plane of y = 25di are
just located at one side of the current sheet. Crossing the current sheet along the red
line indicated in the top left plot, we observe that the ?eld lines are connected on the
side of current sheet and unconnected on the other side. Such asymmetry is simply
due to the ?eld line geometry geometry of this ?component reconnection?, i.e., the
reconnection with a ?nite guide ?eld.
77
Figure 4.5: The con?guration of magnetic ?eld lines near the current sheet at t =
50??1i (the top left) and t = 50??1i (the top right), and the contour plots of Bx (the
bottom left) and Vx (the bottom right) in the xz plane of y = 25di at t = 50?i.
78
To investigate the properties of propagating waves, we now examine the propaga-
tion speed of the wave front. Fig.4.6 shows the relationship between the time t and the
distance L from the center of the di?usion region (0,0,0) to the position of the wave
front. Obviously t and L are linearly correlated with a linear ?t of L = 0.89t?11.99.
Therefore the propagation speed is 0.89VA0 that is close to the lobe Alfv?en speed VA0.
Note that the propagation is mainly located outside the edge of the current sheet,
where the Alfv?en speed is slight smaller than the lobe Alfv?en speed.
Figure 4.6: The relationship between time t and the distance L from the center of
the di?usion region (0,0,0) to the position of wave front.
In Fig.4.7, the contour plots of Bx and N in the x-z plane of y = 16di at t =
80?i are presented. The contours of Bx exhibit clear structure of perturbations,
while the density perturbations are also observed at the same position. Since the
perturbations propagate from the di?usion region that has the characteristic scale
length ? ?i, it makes sense that the spacial scale length of the perturbation Bx is
? 3di, corresponding to kx?i ? 0.62. We expect the waves to have certain kinetic
e?ects.
79
Figure 4.7: The contour plots of Bx and density N in the x-z plane of y = 16di at
t = 80?i. The red dot stands for the position (4, 16, -20).
Fig.4.8 shows the time evolutions of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ez, E?,
B, Ex?b, B?, and the hodogram of Bx ?Bx?b at the position (x,y,z) = (4,16,?20),
as indicated by the red dot in Fig.4.7. Among them, B? is calculated by B ? b,
where b is the initial magnetic ?eld direction, and the subscript x ? b stands for
the direction perpendicular to x and b. We see that the transverse components Bx,
Vx, and Ey, Ez are well correlated. The ratio Ex?b/Bx is ? 0.8VA0 close to the lobe
Alfv?en speed. During the whole time, the magnitude of the magnetic ?eld B stays
almost constant, as well as the parallel component of the magnetic ?eld B?. The
wave appears an Alfv?en mode. Moreover, the spatial scale length of perturbations
along the x direction is comparable to ?i, and E? is correlated with the perturbation
of the local density that develops after the arrival of these perturbations. The slight
increase of density is also due to the ion kinetic e?ect. Recall that due to the ion
kinetic e?ect, the Alfv?en wave becomes compressible and generates an electric ?eld
parallel to the local magnetic ?eld. Since b points into the paper in the hodogram of
Bx?Bx?b, the waves rotate along the counterclockwise direction and have a lefthand
polarization.
80
Figure 4.8: The time evolution of Bx, By, Bz, N, Vx, Vy, Vz, Ex, Ey, Ez, E?, B, Ex?b,
B?, and the hodogram of Bx ?Bx?b at the position (4, 16, -20) as indicated by the
red dot in Fig.4.7.
81
The contour plots of Vz at x = 0 (left) and x = 2di (right) are shown in Fig.4.9.
In the center of the current sheet, the plasma is accelerated by the tension fore
caused the reconnected ?eld lines, and the accelerated plasma is mainly towards ?z
directions. Outside the current sheet, the Alfv?en wave propagates along the ?eld lines,
and the plasma is accelerated through the wave during its propagation. Therefore
the accelerated plasma is mainly along the ?eld lines as shown in Fig.4.8.
Figure 4.9: The contour plots of Vz at x = 0 (left) and x = 2di (right) at t = 80??1i .
Fig.4.10 shows the contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the
x-z plane of y = 0. Note the x-z plane of y = 0di is in the center of X-line. We
can see that the resulting structures here are very di?erent from the 3D run with an
in?nitely long X-line. No clear leading bulge exists during the reconnection, although
an accelerated plasma layer is formed near the center of current sheet where the
parallel heating is observed. Moreover, we do not ?nd the discontinuity fronts on the
edges of the reconnection layer, as in 2D simulations and 3D simulations for cases
with an in?nite X-line. In 3D simulation with the limited X-line, the component Bx
is a shear perturbation to the local magnetic ?eld. In such case with ky ?= 0, such
82
3D perturbations propagate with a group velocity along the ?eld lines. Therefore, no
steady discontinuities or shocks are formed.
Figure 4.10: The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z
plane of y = 0 when t = 150??1i .
Fig.4.11 shows the contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z
plane of y = 20di. The propagations of perturbations generated by the reconnection
83
accelerate the local plasma and form a structure similar to that in the x-z plane
of y = 0di, although no reconnection actually occurs in this x-z plane. Across the
current sheet, the layers are asymmetric due to the reason discussed already. Again,
the parallel heating occurs in the layers. The enhancement of density is found in the
layers because of the acceleration of the plasma.
Cases with various lengths of X-line are also simulated. It is found that the
critical X-line length for the existence of steady-discontinuities is ? 30di.
4.4 Summary
We have used the 3D hybrid simulation model to study the generation of low
frequency waves by the reconnection under a ?nite guide ?eld. The simulation results
are summarized below.
For the case with an in?nite X-line:
1. 3D simulation results are similar to those of 2D cases. Two quasi-steady rota-
tional discontinuities are formed behind a plasma bulge. The plasmas are accelerated
across the rotational discontinuities.
2. The contour plots show that there are perturbation structures in the transition
region between the leading bulge and the steady discontinuities. The transverse
perturbations lead to the presence of short-wavelength ?eld-aligned structures, in
which E? ?= 0 due to the ion kinetic e?ects.
For the cases with a ?nite X-line:
3. It is found that the perturbations caused by the reconnection propagate along
the local magnetic ?eld lines.
4. The wavefront of perturbations propagates with the local Alfv?en speed. The
shear components dominate the wavefront of perturbation which are seen as locations
far from the reconnection region, and the ratio between the corresponding transverse
electric ?eld and magnetic ?eld components is almost equal to the local Alfv?en speed.
84
Figure 4.11: The contour plots of Bx, By, Bz, Vx, Vy, Vz, N, T? and T? in the x-z
plane of y = 20di when t = 150??1i
85
Thesewavesareidenti?edaskineticAlfv?enwaves, inwhichthewavenumberkx?i ? 1.
Our simulation indicates that the ?nite X-line length leads to a structure of the
reconnection layer that is very di?erent from the cases with an in?nite X-line, in
which steady discontinuity fronts are dominant.
5. Parallel electric ?eld E? is generated in the kinetic Alfv?en waves, and the
structure of E? has a strong correlation with the structures of the shear components.
6. The critical length of X-line that leads to the Petschek-type steady shocks
or discontinuities is found to be ? 30di. Our study indicates that a ?nite extent of
X-line alters the wave structure of the reconnection layer.
The generation of kinetic Alfv?en waves in the general cases of reconnection may
provide an important mechanism for the particle heating and acceleration as these
waves propagate along ?eld lines to the ionosphere.
86
Chapter 5
Summary
Magnetic ?eld reconnection provides an e?cient mechanism for the transfer of
magnetic energy to plasma heat and kinetic energy. The kinetic physics of the fast
reconnection in collisionless plasmas, however, is still poorly understood due to the
complex plasma dynamics that is largely caused by the disparate spatial and tem-
poral scales between the electron and ions. The 2D Petschek reconnection model is
the ?rst fast reconnection model, which provides a mechanism of reconnection that
is fast enough to explain the phenomena caused by the reconnection. In order to
fully understand the validity of the Petschek model, two basic questions should be
answered. One is how the fast reconnection is triggered by micro instabilities, and
the other is how the reconnection layer will be modi?ed or altered in the general 3D
geometries in reality.
To tackle the problem due to the vast di?erence between the electron and ion
scales and investigate the physics of reconnection from the micro to macro scales, we
use the newly developed GeFi particle scheme to simulate the reconnection in a 2D
current sheet with a ?nite guide ?eld and a realistic mass ratio. In this thesis, the
GeFi model of Lin et al. [2005] has been improved for the nonuniform current sheets.
Tearing mode instabilities are believed to play a fundamental role to trigger
the fast reconnection. As a ?rst step to apply the GeFi model to the reconnection
research, we have developed a linear eigenmode theory of the tearing instability to
benchmark the improved GeFi scheme. The GeFi results have also been compared
with the asymptotic matching theory of Drake and Lee [1977a]. Furthermore, the
kinetic e?ects that are missing in previous theoretical studies of the tearing mode
87
have been examined through the GeFi simulations. After the code benchmark, we
have then carried out linear and nonlinear GeFi simulations of the tearing instability
in magnetic reconnection.
While the application of the GeFi model to magnetic reconnection has produced
interesting initial results, the complete investigation of 3D nonlinear magnetic re-
connection using the GeFi model is still beyond the capability of the current GeFi
scheme. To explore the properties of low-frequency waves generated by magnetic re-
connection, a hybrid code has been used to simulate the 3D structure of reconnection
layer under a ?nite guide ?eld.
The main results of this thesis are summarized:
The improved GeFi model:
1. Our GeFi particle simulation model, in which the electrons are treated as
gyrokinetic particles and ions are treated as fully kinetic particles, has been improved
and modi?ed to allow the the existence of modes with wavelengths on the same scale
of the background nonuniformity.
2. The linearized GeFi scheme is benchmarked for uniform plasmas, and the
simulation results using the ?f method show that the model can accurately resolve
the physics ranging from Alfv?en waves to lower-hybrid/whister waves, for k? << k?
and ? << ?e.
The tearing mode instability under a ?nite guide ?eld:
1. The GeFi simulation model is benchmarked with tearing mode instabilities
in a simple current sheet and the Harris current sheet. Our results agree with the
eigenmode analysis and the previous studies under simpli?ed conditions of these the-
oretical models. The inhomogeneity terms a?ect the structure of tearing modes, and
cannot be neglected when By ? Bx.
2. In the general cases in which the ion gyroradius is comparable to the current
sheet width, the ion kinetic e?ects are found to play a signi?cant role on the growth
88
rate of tearing modes. In such cases, the Drake and Lee theory [1977a] and the
eignemode analysis based on the drift kinetic assumption are invalid.
3. The nonlinear GeFi simulation shows that the ion e?ect can broaden the
width of saturated magnetic islands from the electron dynamic spatial scale to the
ion gyroradius. The case with multiple tearing modes shows that the coalescence of
magnetic islands can form a larger saturation island than the case of a single tearing
mode.
The generation of low frequency waves by reconnection under a ?nite guide ?eld:
1. For the case with a in?nitely long X-line, 3D simulation results are similar to
those of 2D cases. In addition to the steady rotational discontinuities in the out?ow
region, wave perturbations are also found between between the leading plasma bulge
and the steady discontinuities. The transverse perturbations lead to the presence
of short-wavelength ?eld-aligned structures, in which E? ?= 0 due to the ion kinetic
e?ects.
2. For the case with a ?nite extent of X-line, it is found that the perturbations
caused by the reconnection propagate along the local magnetic ?eld lines. The wave-
front of perturbations propagates with the local Alfv?en speed. The shear components
dominate the wavefront of perturbation which are seen at locations far from the re-
connection region. These waves are identi?ed as kinetic Alfv?en waves, in which the
wave number k??i ? 1. Parallel electric ?eld E? is generated in the kinetic Alfv?en
waves. Our simulation indicates that the ?nite X-line length leads to a structure of
the reconnection layer that is very di?erent from the cases with an in?nite X-line as
as in the Petschek model, in which steady discontinuity fronts are dominant.
89
Bibliography
[1] Biernat, H.K., Heyn, M.F., Semenov, V.S., Farrugia, C.J., The structure of
reconnection layers: Application to the Earths Magnetopause, J.Geophys.Res.,
94, 287, 1989
[2] Birn, J., et al., Geospace Environmental Modeling(GEM) magnetic reconnection
challenge, J.Geophys.Res., 106, 3715, 2001
[3] Birn, J., E. R. Priest, Reconnection of Magnetic Fields: Magnetohydrodynamics
and Collisionless Theory and Observations, Cambridge, NewYork, 2007
[4] Biskamp, D., Magnetic reconnection via current sheets, Phys.Fluids, 29, 1520,
1986
[5] Burden, R.L., Faires, J.D., Numerical Analysis , Brooks/Cole, 2005
[6] Chen, Y., Parker, S.E., A ?f particle method for gyrokinetic simulations with
kinetic electrons and electromagnetic perturbations, J.Comput.Phys, 189, 463,
2003
[7] Drake, J.F. Lee, Y.C., Kinetic theory of tearing mode instabilities, Phys.Fluids,
20, 1341, 1977a
[8] Drake, J.F. Lee, Y.C., Nonlinear evolution of collisionless and semicllisional tear-
ing modes, Phys.Rev.Lett. 39, 453, 1977b
[9] Daughton, W., Electromagnetic properties of the lower-hybrid drift instability
in a thin current sheet, Phys. Plasmas 10, 3103, 2003
[10] Daughton, W., Karimabadi, H., Kinetic theurgy of collisionless tearing at the
magnetopause, J.Geophys.Res., 109, A0317, 2005
[11] Fireman, E.A., Chen., L., Nonlinear gyrokinetic equations for low-frequency elec-
tromagnetic waves in general plasma equilibria, Phys.Fluids, 25, 502, 1982
[12] Furth, H.P., Killeen, J., Rosenbluth, M.N., Finite resistivity instabilities of a
sheet pinch, Phys.Fluid, 6, 459, 1963
[13] Furth, H.P., Rutherford, P.H., Selberg,H., Tearing mode in the cylindrical toka-
mak, Phys.Fluid, 16, 1054, 1973
90
[14] Heyn, M.F., Biernat, H.K., Semenov, V.S., The structure of magnetic reconnec-
tion layer, J.Geophys.Res., 90, 1781, 1981
[15] Hans, J.P., Poedts, S. Principles of Magnetohydrodynamics: With Applications
to Laboratory and Astrophysical Plasmas, Cambridge Press., 2004
[16] Harris, E.G., On a plasma sheath separating regions of one directional magnetic
?eld, Nuovo Cimento 23, 115, 1961
[17] Hockney, R.W., Eastwood, J.W., Computer Simulation Using Particles, Mcgraw-
Hill, New York, 1981
[18] Hoshino, M., The electrostatic e?ect for the collisionless tearing mode,
J.Geophys.Res., 92, 7368, 1988
[19] Landau, L.D., Lifshitz, E.M., Electrodynamics of Continuous Media, Pergamon
Press., 1960
[20] Katanuma, I., Kamimura, T., Simulation studies of the collisionless tearing in-
stabilities Phys.Fluid, 23, 2500, 1980
[21] Karimabadi, H., Daughton, W., Quest, B.K., Physics of saturation collisionless
tearing mode as a a function of guide ?eld, J.Geophys.Res., 110, A010749, 2005
[22] Karimabadi, H., Huba, J., KraussVarban, D., On the generation and structure
of the quadrupole magnetic ?eld in the reconnection process, Geophys.Rev.Lett.
31, L070806, 2004
[23] KrasussVarban, D., Omidi, N., large-scale hybrid simulations of the magnetotail
during reconnection, Geophys.Res.Lett., 22, 3271, 1995
[24] Karimabadi, H., Daughton, W., Quest, K.B., Antiparallel versus component
merging at the magnetopause: Current bifurcation and intermittent reconnec-
tion, J.Geophys.Res., 10, A03213, 2005b
[25] Karimabadi, H., Daughton, W., Quest, K.B., Physics of saturation of collisionless
tearing mode as a a function of guide ?eld, J.Geophys.Res., 110, A03214, 2005c
[26] KrasussVarban, D., Omidi, N., large-scale hybrid simulations of the magnetotail
during reconnection, Geophys.Res.Lett., 22, 3271, 1995
[27] Lee, W.W., Gyrokinectic approach in particle simulation, Phys.Fluids, 26, 556,
1983
[28] Lee, W.W., Gyrokinectic particle simulation model, J.Comput.Phys., 72, 243,
1987
[29] Lin, Y., Lee, L.C., Structure of reconnection layers in the magnetosphere,
Space.Sci.Rev., 65, 59, 1994a
91
[30] Lin, Y., Lee, L.C., Reconnection layer at the ?ank mangetopause in the presence
of shear ?ow, Geophys.Res.Lett., 21, 59, 1994b
[31] Lin, Y., Wang, X.Y., Chen, L., Lin, Z.H., A gyrokinetic electron and fully kinetic
ion plasma simulation model, Plasma Phys.Control.Fusion, 47, 657, 2005
[32] Lin, Y., Wang, X.Y., Three-dimensional global hybrid simulation of dayside
dynamics associated with the quasi-parallel bow shock, J.Geophys.Res., 110,
A12216, 2005
[33] Lin, Y., Wang, X.Y., Three-dimensional global hybrid simulation of dayside
dynamics associated with the quasi-parallel bow shock, J.Geophys.Res., 115,
A04208, 2010
[34] Lin, Y., Lee, L.C., Chaos and ion heating in a slow shock, Geophys.Res.Lett, 18,
1615, 1991
[35] Lin, Y., Lee, L.C., Kennel, C.F., The role of intermediate shocks in magnetic
reconnection, Geophys.Res.Lett, 19, 229, 1992
[36] Lin, Y., Lee, L.C., Simulation study of the Riemann problem associated with
magnetotail reconnection, J.Geophys.Res., 100, 19,277, 1995
[37] Lin, Y., Swift, D.W., A two dimensional hybird simulation of the magnetotaill
reconnection, J.Geophys.Res., 101, 19,859, 1996
[38] , Lin, Y., Wang, X.Y., Chen, L., Lu, X., and Kong, W., An improved gyrokinetic
electron and fully kinetic ion particle simulation scheme: benchmark with a linear
tearing mode, Plasma Phys. Control. Fusion., 53, 054013, 2011
[39] Lottermoser, R.F., Scholar, M., Undriven magnetic reconnection in magnetohy-
drodynamics and Hall magnetohydrodynamics, J.Geophys.Res., 102, A3, 4875,
1997
[40] Lottermoser, R.F., Scholar, M., Matthews, A.P.,Ion kinetic e?ects in magnetic
reconnection, J.Geophys.Res., 103, 4547, 1998
[41] Masaaki, Y., Kulsrud, R., Ji, H.T., Magnetic Reconnection, Rivews of Modern
Physics, 82, 603, 2010.
[42] Mandt, M.E., Denton, R.E., and Drake, J.F., Transition to whistler mediated
magnetic reconnection, Geophys.Res.Lett, 21, 73, 1994
[43] Parker, E.N., Sweet?s mechanism for merging magnetic ?eld in conducting ?uids,
J.Geophys.Res., 62, 509, 1957
[44] Parker, S.E., Lee, W.W., A fully nonlinear characteristic method for gyrokinetic
simulation, Phys.Fluids B, 5, 77, 1993
92
[45] Petschek, H.E., Magnetic ?eld annihilation, in AAS-NASA Symposium on the
Physics of Solar Flares, NASA Spec. Publ., SP-50, 425, 1964
[46] Priest, E., Forbes, T., Magnetic Reconnection: MHD Theory And Applications,
Cambridge Press., 2006
[47] Prtchett, E., Collisionless magnetic reconnection in a three open system,
J.Geophys.Res., 106, 3783, 2001
[48] Prtchett, E., Onset of saturation of guide-?eld magnetic reconnection,
Phys.Plasmas, 12, 062301, 2005
[49] Richter, P., Sholer, M., on the stability of rotational discontinuities,
J.Geophys.Res., 16, 1257, 1989
[50] Richter, P., Sholer, M., on the stability of rotational discontinuities with large
phase angles, Adv.Space.Res., 11, 111, 1991
[51] Ricci, P., Brackbill, U. B., Doughton, W., and Lapenta, G., New role of the
lower-hybrid drift instability in the magnetic reconnection, Phys. Plasmas 12,
5501, 2005
[52] Ricci, P., Brackbill, U. B., Doughton, W., and Lapenta, G., Collisionless mag-
netic reconnection in the presence of a guide ?eld, Phys. Plasmas 11, 4102, 2008
[53] Rogers, B.N., Denton, R.E., Drake, J.F., Shay, M.A., The role of dispersive waves
in collisionless magnetic reconnection, Phys.Rev.Lett., 97, 195004, 2001
[54] Rutherford, P.H., Frieman, E.A., Drift instabilities in general magnetic ?eld
con?guration, Phys.Fluids, 11, 569, 1968
[55] Rutherford, P.H.,Nonlinear growth of tearing mode, Phys.Fluids, 16, 1903, 1973
[56] Sweet, P.A., The neutral point theory of solar ?ares, in Electromagnetic Phe-
nomena in Cosmic Physics, 135, 1958.
[57] Swift, D.W., Lee, L.C., Rotational discontinuities and the structure of the mag-
netopause, J.Geophys.Res., 88, 111, 1983
[58] Swift, D.W., Use of hybrid code to model the Earth?s magnetosphere, Geo-
phys.Res.Lett., 22, 311, 1995
[59] Sato, T.,Strong plasma acceleration by slow shocks resulting from magnetic re-
connection, J.Geophys.Res., 84, 7177, 1979
[60] Shi, Y., Lee, L.C., Structure of the reconnection layer at the dayside magne-
topause, Planet.Space.Sci., 38, 437, 1990
[61] Scholer, M., Undriven magnetic reconnection in an isolated current sheet,
J.Geophys.Res., 94, 15099, 1989
93
[62] Sonnerup, B., Magnetic ?led reconnection, in Solar System Plasma Physics, 46,
Elsevier, New York, 1974
[63] Taylor, J.B., Hastie, R.J., Stability of general plasma equilibria-I. Formal theory,
Phys.Fluids, 10, 479, 1968
[64] Terasawa, T., Hall current e?ect on tearing mode stability, Geophys.Res.Lett.,
10, 475, 1983
[65] Ugai, M., Self-consistent development of fast reconnection with anomalous
plasma resistivity, Plasmas.Phys.Control.Fusion, 26, 1549, 1984
[66] Ugai, M., Shumizu, T., Computer studies on the fast reconnection mechanism
in a sheared ?eld geometry, Phys.Plasmas, 1, 296, 1993
[67] Wang, X.Y., Lin, Y., Chen, L., Lin, Z.H., A particle simulation of current sheet
instabilities under ?nite guide ?eld, Phys.Plasmas 15, 072103, 2008
[68] Wu, C.S., Zhou, Y.M., Tsai, S.T., and Gao, S.C., A kinetic cross ?eld streaming
instability, Phys. Fluids 26, 1259,1983
[69] Wesson, J.A., Finite resisitivity instability of a sheet pinch, Nuclear Fusion, 18,
130,1966
[70] Wan, W., Chen, Y., and Parker, S.E., Gyrokinetic ?f simulation of the collision-
less and semi-collisional tearing mode instability, Phys.Plasmas 12, 012311,2005
[71] Yan, M., Lee, L.C., Priest, E.R., Fast magnetic reconnection with small shock
angles, J.Geophys.Res., 97, 8277, 1992
94