Doctoral Thesis

Studies on Large-Scale Evolutionof Magnetic Reconnection Using
Full Particle Simulations with
Adaptive Mesh Refinement
Technique

Keizo Fujimoto
Department of Geophysics, Graduate School of
Science, Kyoto University

December 20, 2005

Abstract
Magnetic reconnection is widely believed to play an important role in magnetosphericsubstorms and solar ﬂares as a fast conversion process of the magnetic energy to plasma
kinetic energy. However, it has been diﬃcult to describe its large-scale dynamics byusing kinetic simulations due to limited computer resources. The main unresolved issues
in the reconnection processes include the onset mechanism, the reconnection rate in thesteady state, and the energy transport mechanism. The diﬃculty in approaching such
issues lies in the fact that macro-scale structures can be strongly aﬀected by localized
micro-scale processes, in which the kinetic treatment of plasma is required.
In order to overcome this diﬃculty, we ﬁrst develop a new 2-1/2 dimensional electromagnetic particle code that employs adaptive mesh reﬁnement (AMR) technique anda particle splitting algorithm. The AMR technique subdivides and removes cells dynamically in accordance with a reﬁnement criterion and enhances the spacial resolution
locally. On the other hand, the particle splitting algorithm divides ﬁnite-size superparticles locating in ﬁne cells, controling the number of particles per cell in order to suppress
numerical noise arising when the number of particles per cell gets too small. We perform
some test simulations and compare three runs without the AMR, with the AMR, andwith both the AMR and the particle splitting. These results indicate not only that the
AMR and particle splitting algorithms are successfully applied to the conventional particle codes, but also that they are quite eﬀective to achieve high-resolution simulations on
the evolution of the current sheet associated with magnetic reconnection. Thus the newcode allows us to implement a large-scale kinetic simulation of magnetic reconnection.
One of our goals by conducting large-scale simulations is to investigate a long time
evolution of the diﬀusion region and to see what supports fast reconnection in a fully
kinetic system. Our results show that fast reconnection is quickly achieved through theHall eﬀects but the system does not reach steady state and the reconnection rate decreases signiﬁcantly. The key process responsible for slowing the reconnection processes
is the extension of the electron diﬀusion region in association with the enhancement of
the polarization electric ﬁeld directing toward the neutral sheet in the electron inﬂow
region. The polarization electric ﬁeld is caused by the inertia diﬀerence between ions andelectrons, and enhanced by the meandering motions of the background ions. The inﬂow
electrons are forced by the polarization ﬁeld to move toward the out-of-plane direction
by the E × B drift and enhances the out-of-plane current density in the electron inﬂow
region. Since the current density is imposed uniformly along the upstream edge of theABSTRACT

ii

electron diﬀusion region, the magnetic ﬁeld lines become mostly parallel to the x direction near the X-line. As a result, the electron meandering region so the electron diﬀusion
region extends along the x direction. In order to conﬁrm the role of the polarization electric ﬁeld, we compare the simulation results with mi /me = 1 and 100. It is found that
(1) a steady-state reconnection is achieved in the system where the polarization electricﬁeld does not arise, (2) a large reconnection rate is obtained even in the system without
the Hall eﬀects. It is suggested that the anomalous resistivity due to the Buneman-typeinstability might be required to support a steady-state fast reconnection.
Another goal is to reveal energy transport processes associated with magnetic reconnection, that is, where and how plasma is accelerated and heated. In the present study,
we pay attention to an electron heating mechanism in the plasma sheet-lobe boundary
region, where the truncated (i.e., ﬂat-topped) distribution function of electrons is often
observed in the Earth magnetotail during magnetic reconnection. It is found that theelectron two-stream instability between the background cold electrons and the strong
beam electrons with high perpendicular temperature should be responsible for the formation of the ﬂat-topped electrons. Electrons are trapped by the electrostatic potential
wells resulting from the electron two-stream instability and scattered along the ambient magnetic ﬁeld, so that the truncated distribution function is formed in the parallel
direction. We show that the electron heating through this process occurs in the extensive region of the boundary region and the heating level is consistent with satellite
observations. The strong beam electrons consist of two components originating fromdistinct regions. Some electrons come from the opposite boundary region of the plasma
sheet along the magnetic ﬁeld line, forming the cold component of the electron beam.
The others are originating from the lobe region, passing through the electron diﬀusion
region. Furthermore, we demonstrate that the Electrostatic Solitary Waves (ESW) canbe evolved from the electron two-stream instability in association with magnetic reconnection. Such waves should be observed in the boundary region during reconnection.

CONTENTS

iii

Contents
1 General Introduction

1

1.1 Role of Magnetic Reconnection in Solar System . . . . . . . . . . . . . .
1.2 Fluid Models Explaining Fast Reconnection . . . . . . . . . . . . . . . .

1
6

1.3 Kinetic Aspects Beyond the MHD Framework . . . . . . . . . . . . . . .
1.4 Objective of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . .

912

2 Electromagnetic Full Particle Code With Adaptive Mesh Reﬁnement
Technique2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14
14

2.2 Description of the Code . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Data Structure and Cell Generator . . . . . . . . . . . . . . . . .
2.2.2 Particle Splitting and Coalescence . . . . . . . . . . . . . . . . . .

17
17
19

2.2.3
2.2.4

Particle and Force Weighting within the Hierarchical Cell . . . . .
Integration of the Field Equations . . . . . . . . . . . . . . . . . .

21
22

2.2.5 Integration of the Equations of Motion . . . . . . . . . . . . . . .
2.3 Test Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26
26

2.3.1
2.3.2
2.3.3

Landau Damping of the Langmuir Waves . . . . . . . . . . . . . .
Wave Reﬂections at the Boundaries of the Reﬁned Regions . . . .
Nonlinear Evolution of the Plasma Sheet . . . . . . . . . . . . . .

26
28
30

2.4 Summary and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

3 Time Evolutions of the Electron Diﬀusion Region and the Reconnection
Rate in Large and Fully Kinetic System40
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40
42

3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Time Evolution of the Reconnection Rate . . . . . . . . . . . . .

44
44

3.3.2
3.3.3

Structure Change Around the Electron Diﬀusion Region . . . . .
Impact on the Reconnection Rate . . . . . . . . . . . . . . . . . .

46
51

CONTENTS

3.3.4

iv

Comparison With the Case of mi /me = 1 . . . . . . . . . . . . . .

56

3.4 Summary and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

4 A Generation Mechanism of Electrostatic Waves and Subsequent Electron Heating in the Plasma Sheet-Lobe Boundary Region During Magnetic Reconnection
62
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Simulation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62
64

4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Parallel Electron Heating in the Plasma Sheet-Lobe Boundary Re-

66

gion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Instability Responsible for the Electron Heating . . . . . . . . . .

66
71

4.3.3 Origin of the Electron Beam . . . . . . . . . . . . . . . . . . . . .
4.3.4 Electrostatic Solitary Waves Evolved from the Beam Instability .
4.4 Summary and Disccusion . . . . . . . . . . . . . . . . . . . . . . . . . . .

74
76
79

4.3.2

5 Conclusions and Further Studies

82

CHAPTER 1
General Introduction
Magnetic reconnection is known as a fast conversion process of magnetic energy
into plasma kinetic energy and believed to be a phenomenon of considerable importancein space physics and astrophysics. In this part, we brieﬂy review theoretical and observational studies devoted toward understanding the physical mechanism of magnetic
reconnection. In section 1.1, we describe why it is considered to be so important in thesolar system, especially in the Earth magnetosphere, followed by discussions on magnetic reconnection proper. In sections 1.2 and 1.3, some important models explaining
fast energy release due to magnetic reconnection are described. We then summarize themain problems that remain and motivate this study.

1.1

Role of Magnetic Reconnection in Solar System

We ﬁrst introduce the principal physics of magnetic reconnection and explain why it
is such an important process in the solar system, where the particles are highly ionized
and form a plasma state.
We consider a plasma obeying a simple Ohm’s law of the form
E + V × B = η ef f J .

(1.1)

The left-hand side of (1.1) represents the electric ﬁeld in the rest frame of the plasma
ﬂowing with velocity V in the magnetic ﬁeld B, while the right-hand side is an Ohmic
current term proportional to eﬀective scalar resistivity η ef f . The parameter η ef f isexpressed by using either inter-particle or wave-particle collision frequencies asη ef f =mν ef f
,
ne2

(1.2)

where m and n are the mass and the number density of the plasma, and ν ef f is theeﬀective collision frequency. Solving (1.1) with respect to E and substituting it into

CHAPTER 1. GENERAL INTRODUCTION

2

Faraday’s law, and then using Ampère’s law neglecting the displacement current, one
can derive the induction equation which determines how B developes in timeη ef f 2∂B= ∇ × (V × B) +
∇ B.
∂t
µ0

(1.3)

The ﬁrst term describes convection of the magnetic ﬁeld lines with the plasma ﬂow, while
the second term represents diﬀusion of the magnetic ﬁeld through the plasma such as toreduce the gradients. If the ﬁrst term dominates so that the second term is negligible,
the plasma is “fozen-in” to the magnetic ﬁeld such that a ﬂuid elements which is initiallylocated on a ﬁeld line remain on the same ﬁeld line forever. On the other hand, if thesecond term dominates the ﬁrst term, the equation (1.3) becomes a diﬀusion equationand there is no coupling between the magnetic ﬁeld and the plasma ﬂow.
The relative importance of the two terms on the right-hand side of (1.3) is conveniently compared with the magnetic Reynolds number
Rm =

µ0 UL
|∇ × (V × B)|,≈ ef feffη|η /µ0 ∇2 B|(1.4)

where U is a characteristic ﬂow speed and L is a characteristic length of the system. If Rm
is large, convection dominates and frozen-in ﬂow prevails, while if Rm is small, diﬀusion
dominates and the coupling to the plasma ﬂow is weak. If we consider a global structureof the Earth magnetosphere, L is the order of a few Earth radii (Re) (the scale size of themagnetospheric cavity) and U is estimated as 100 km/s (a typical magnetospheric ﬂow
speed), which leads to Rm  1011 . In the solar ﬂares, on the other hand, L ∼ 104 km (thesize of a ﬂare system) and U ∼ 102 m/s (a typical outﬂow speed) result in Rm ∼ 108 . In
both cases, the values of Rm are so large that on these scale sizes the ﬁeld convection is
overwhelmingly dominant and that the eﬀect of magnetic diﬀusion is entirely negligible.
Thus it seems that the diﬀusive process in the magnetosphere and the solar ﬂares isnegligible in considering their global structures.
In such a limit neglecting the diﬀusion term in (1.3), the behavior of plasma is strongly
constrained in such a way that the particles always remain tied to the same magnetic
ﬁeld line as it convects with the plasma ﬂow, and cross-ﬁeld mixing of plasma elements
is suppressed. Therefore when two initially separate plasma regimes come into contact
with each other, one can expect that these two plasmas do not mix, instead a thinboundary layer should be formed between them separating two plasma populations andmagnetic ﬁelds. Generally the boundary layer constitute a current sheet, since the ﬁelds
on either side of the boundary will have diﬀerent strengths and orientations tangential tothe boundary. Actually, a plasma ﬂow that is frozen-in to the magnetic ﬁeld is generated
on the solar surface and prevails in the solar system. The existence of so-called solar

CHAPTER 1. GENERAL INTRODUCTION

3

wind, which carries magnetized plasma, was speculated by Alfvén (1957) and predicted
theoretically by Parker (1958). Since the Earth magnetic ﬁeld contacts with the solar
wind, there indeed exists a boundary layer between them called magnetopause and it
forms a current sheet (Chapman and Ferraro, 1931). Also in the Earth magnetotail
and the solar ﬂares, current sheets are formed between the anti-parallel magnetic ﬁelds
and divide the two regions. Thus the perfect conductive approximation of plasma ﬂows
in the magnetic ﬁeld results in thin current sheets separating distinct plasma regimes.
The important is that it is not guaranteed that the eﬀect of diﬀusion can be neglected in
discussing the physics of these boundaries, even though it may be negligible in describing
the large-scale behavior within each plasma regime. Therefore we next consider theeﬀects of diﬀusion in these boundary layers.
For simplicity, we assume a boundary on which equal and opposite ﬁelds of strength
B0 contact with each other separated by a magnetic neutral sheet. The magnetic ﬁeld
lines in either side are along x axis so that Bx = ±B0 for z ≷ 0. Then a current ﬂowing
along y axis is generated. Plasmas in either side of the boundary layer ﬂow toward
the neutral sheet with speed V , but the ﬂow is negligible in the boundary layer suchthat (1.3) reduces to a diﬀusion equation. Now we set the half width of the equilibrium
current sheet as λ. The plasma inﬂows are generated by the E × B drifts of the particles
so that V  Ey /B0 . For steady state system, Faraday’s law (∇ × E = 0) requires thatthe electric ﬁeld Ey should be spatially uniform. Thus the electric ﬁeld imposed outside
the current sheet
Ey  V B0
can be equated by that in the current sheet where magnetic diﬀusion is dominant (i.e.,the diﬀusion region)
Ey  η ef f Jy .

(1.5)

Using Ampère’s law (Jy  B0 /µ0 λ), we obtainRm =µ0 V λ 1,
η ef f(1.6)

so that the width of the current sheet is suﬃciently small and the diﬀusion becomes
important.
In a practical sense, the eﬀective topology of the magnetic ﬁeld lines is describedin Figure 1.1 during a steady magnetic reconnection. The magnetic ﬁeld lines in eitherside of the boundary layer convect toward the boundary with the plasma ﬂows, thenthey are reconnected to a ﬁeld line convecting from the opposite side of the boundary

CHAPTER 1. GENERAL INTRODUCTION

4

v

E

TKZ

CV
RCT

5G

2λ
E
&KHHWUKQP4GIKQP:V[RGPGWVTCNNKPG
Figure 1.1: Magnetic ﬁeld lines (light blue lines) and ﬂow vectors (red arrows) aroundthe neutral line during a steady-state magnetic reconnection. Hatched area represents
the diﬀusion region where magnetic diﬀusion dominates the eﬀect of plasma convection.

in the diﬀusion region. The reconnected ﬁeld line is severely bent, so that the magnetic
tension accelerates the inﬂowing plasmas away from the neutral line along the boundary
layer. Thus the magnetic energy is converted to the plasma kinetic energy in the form ofaccelerated and heated plasma jets ﬂowing away from the neutral line. This descriptionof the ﬁeld and plasma behavior was ﬁrst suggested by Dungey (1953).
Next, let’s consider the interaction between the solar wind and the Earth magnetosphere, in which we concluded in the above discussion that the perfect conductivity
approximation was quite valid. If the magnetic ﬁeld in the solar wind has a southward
component, it is likely that the northward magnetic ﬁeld originated from the magnetosphere reconnects with the solar wind magnetic ﬁeld in the dayside magnetopause.
In Figure 1.2, a convection model of the magnetosphere is shown in association withthe dayside reconnection after Dungey (1961). Magnetic reconnection in the daysidemagnetopause gives rise to open ﬁeld lines connecting from the polar regions to the solar wind, allowing the solar wind plasma to enter the magnetosphere. The solar windﬂow carries these ﬁeld lines from the dayside to the night side of the magnetosphereand stretches them out into the magnetotail lobe. The oppositely-directing magnetic
ﬁelds in the northern and southern tail lobes form a current sheet between them. Subsequently magnetic reconnection occurs in the current sheet and reconnects the open

CHAPTER 1. GENERAL INTRODUCTION

5

Figure 1.2: Magnetospheric model of the Earth. Magnetic reconnection occurs at thedayside magnetopause and the magnetotail current sheet, which facilitates the entry of
the solar wind momentum into the magnetosphere driving a large-scale internal convection. (Adopted from Cowley (1985))

ﬁeld lines of the lobes, generating closed ﬁeld lines in the earthward side. These closed
ﬁeld lines are returned back to the near Earth magnetosphere with the plasma ﬂow. If
the dayside reconnection does not occur, the magnetospheric ﬁeld is conﬁned to a cavity
surrounding the Earth as described by Chapman and Ferraro (1931). In such a case, themagnetospheric and solar wind plasmas are separated by the magnetopause, and plasma
mass and momentum transfer across the boundary is weak. Thus magnetic reconnection in the Earth magnetosphere allows eﬃcient entry of solar wind momentum into themagnetosphere, driving a large-scale internal convection.
The importance of magnetic reconnection lies in the fact that a localized breakdown of
the perfect conductivity can lead to global structure changes via magnetic reconnection,even though in the global scale itself the perfect conductivity approximation is quite
valid. Magnetic reconnection not only allows reconﬁgurations of the magnetic ﬁeld totake place toward states of lower energy, but also facilitates exchanges of plasma mass,
momentum, and energy between plasmas originating from diﬀerent regimes.
The magnetospheric convection model driven by the dayside reconnection has been
established to explain large energy releases in association with magnetic substorms (e.g.,
Dungey, 1961; Axford et al., 1965; Petschek, 1966). This model has been supported
by extensive number of satellite observations and correlated ground observations (e.g.,
McPherron, 1970; Rostoker, 1972; Hones, 1973; Coroniti et al., 1980), in which the im-

CHAPTER 1. GENERAL INTRODUCTION

6

portance of the dayside and nightside reconnections for substorm dynamics was demonstrated. Furthermore, subsequent satellite observations by Geotail spacecraft revealed
that magnetic reconnection in the near-Earth magnetotail occured at X = −20 to −30
Re in the GSM coordinate (Nagai et al., 1998), preceding a few minutes to the substorm onsets identiﬁed by the ground observations (Miyashita et al., 1999). Therefore
it is overwhelming today that magnetic reconnection in the dayside and nightside of the
Earth magnetosphere should play an important role in substorm dynamics (Baker et al.,1999).
Magnetic reconnection is also believed to be of importance in the solar ﬂares as aprocess for eﬀective energy releases of the accumulated magnetic energy into the plasma
kinetic energy (e.g., Giovanelli, 1947; Sweet, 1958; Parker, 1963). The evidence of magnetic reconnection in the solar ﬂares has been demonstrated through many X-ray observations (e.g., Gold, 1961; Sweet, 1969) and emphasized its responsibility for plasmaenergization (Lin and Hudson, 1976).

1.2

Fluid Models Explaining Fast Reconnection

In this and next sections, we discuss detailed processes of magnetic reconnection
and show some important models describing the steady-state conversion processes ofmagnetic energy to plasma kinetic energy. The key point in the following models is
how fast energy conversion as observed in the magnetospheric substorms and the solar
ﬂares can be achieved. In this section, we introduce the well established Sweet-Paker
model (Sweet, 1958; Parker, 1963) and the Petschek model (Petschek, 1964), which wereconstructed in the framework of the magnetohydrodynamics (MHD). We will discuss
more recent models beyond the MHD treatment of plasma in the subsequent section.
We start with the Sweet-Parker model shown in Figure 1.3. This model employs
two-dimensional plane system and assumes that plasma acceleration occurs only in thediﬀusion region. Thus all the inﬂow plasma passes through the diﬀusion region and isaccelerated toward the downstream region. Since the magnetic energy is converted to theplasma thermal energy due to Joule heating (via (1.5)) in the diﬀusion region, the plasmapressure therein must be higher than in the surrounding region. If the diﬀusion region
has a ﬁnite length as shown in Figure 1.3 (shaded region), the pressure gradient shouldarise along the boundary, which accelerates the inﬂowing plasma toward the downstream
region. Since either the ﬂow speed in the inﬂow region and the magnetic ﬁeld strength
in the outﬂow region are expected to be very small, Bernoulli’s law along a stream line

CHAPTER 1. GENERAL INTRODUCTION

B1

7

v1
B2
2λ

v2

2L

Diffusion region
Figure 1.3: Sketch of a reconnection geometry in the Sweet-Parker model.

passing through the diﬀusion region gives
P0 +

B12
1
 P  P0 + mnV22 ,
2µ0
2

(1.7)

where P0 is the plasma pressure in the region far from the diﬀusion region, P the plasmapressure in the diﬀusion region, B1 the magnetic ﬁeld strength in the inﬂow region, V2 theﬂow speed in the outﬂow region, and m and n the mass and number density of plasma,
respectively. Note that we assume the plasma density is spacially uniform. Since weconsider the steady-state reconnection, Faraday’s law (∇ × E = 0) requires that theout-of-plane electric ﬁeld E0 should be spacially uniform, so that
E0 = V1 B1 = V2 B2 ,(1.8)

where V1 is the inﬂow speed, and B2 is the magnetic ﬁeld strength in the outﬂow region.
We further assume that the plasma is incompressible, then the conservation of mass gives
V1 L = V2 λ.

(1.9)

We can obtain the outﬂow velocity V2 from (1.7) asV2  √B1≡ VA ,µ0 nm(1.10)

which indicates that the plasma is accelerated up to the order of the Alfvén velocity
deﬁned in the inﬂow region. We consider the reconnection rate, which means the magnetic ﬂux reconnected per unit time and unit length along the X-line. From Faraday’slaw (∂B/∂t = ∇ × E), it is found that the reconnection rate is simply estimated by theCHAPTER 1. GENERAL INTRODUCTION

8

2L

B1

v1

2λ

B2

v2

Diffusion region
Figure 1.4: Sketch of a reconnection geometry in the Petschek model. The diﬀusion
region has been shrunk to a center of the diagram, instead the plasma is accelerated in
the slow-mode shocks (blue lines) extended from the X-line.

electric ﬁeld imposed along the X-line (Dungey, 1978). Using equations (1.8), (1.9), and(1.10), one can obtain
E∗ 

λ
−1/2
 Rm
,
L

(1.11)

where E ∗ = E0 /VA B1 is the normalized electric ﬁeld, and Rm = µ0 v1 L/η is the magnetic
Reynolds number. We used (1.6) for the second transform. Note that the Sweet-Parker
model assumes the classical resistivity η arising from the Coulomb collisions (Spitzer,
1962).
The Sweet-Parker model realizes a fast plasma outﬂow. However, the reconnection
rate is inversely proportional to the square root of Rm deﬁned by the system size L.
Such a Reynolds number is so large in the Earth magnetosphere and the solar coronas,
as we discussed before, that the reconnection rate becomes very small and insuﬃcient to
explain fast energy releases associated with the magnetospheric substorms and the solar
ﬂares.
In order to solve this problem, Petschek (1964) introduced the ﬁeld conﬁguration
shown in Figure 1.4. In the Petschek model, the diﬀusion region is shrunk to a centerof the reconnection region and most of the inﬂowing plasma does not pass through thediﬀusion region in order to be accelerated. Instead, the plasma can be accelerated in
the slow-mode shocks (blue lines) extended from the X-line. The acceleration occurs
due to the J × B force, diﬀerent from the plasma pressure gradient in the Sweet-Parker
model. An approximate balance of the magnetic shear stress at the separatrices and theCHAPTER 1. GENERAL INTRODUCTION

9momentum change yields
mnV1 V2 

B1 B2
.
µ0

(1.12)

Combination of (1.8) and (3.5) gives (1.10) again. However, since the Petschek model
does not depend on the size of the diﬀusion region in terms of plasma acceleration,
the horizontal scale of the diﬀusion region can be small. This is a key to allow fastreconnection. According to Petschek (1964), the upper limit of the reconnection rate is
∗
Emax


π
[ln(MA2 Rm )]−1 ,8(1.13)

where MA = V1 /VA is the Alfvén Mach number of the inﬂowing plasma. The reconnection
rate in this case is inversely proportional to the logarithm of Rm rather than to its squareroot in the Sweet-Parker system, so it can be drastically increased. For typical Rm values,∗
Emaxlies in the range 0.01 to 0.1, compared with 10−5 to 10−4 for the Sweet-Parker
model. Note that the diﬀusion region in the Petschek model is still signiﬁcant since theactual process of reconnecting the magnetic ﬁeld occurs there.
Though the Petschek model achieves much faster reconnection than the Sweet-Parker
model, it is still diﬃcult to produce suﬃciently high resistivity for realizing large reconnection rate in the collisionless plasma. For example, if the size of the diﬀusion regionis assumed to be the order of the ion inertial length, 106 ∼ 107 times larger resistivity than the conventional Coulomb resistivity is required in order to achieve E ∗  0.1,
which is the case in the magnetospheric substorms. Nevertheless the Petschek model iswell consistent with numerical simulations and satellite observations. For example, theslow-mode shocks are obtained in the MHD simulations (Sato, 1979; Scholer and Roth,1987; Ugai, 1993) and detected by the Geotail observations (Saito et al., 1995). The
fact that the outﬂow velocity from the diﬀusion region reaches the Alfvén velocity isalso conﬁrmed by kinetic simulations (e.g., Shay et al., 1998; Hesse et al., 1999) and the
Geotail observations (Nagai et al., 1998). Thus the remaining question is what is thegeneration mechanism of such high resistivity, much higher than the classical resistivity
due to Coulomb collisions. In the next section, we discuss more recent model to produce
such high ‘eﬀective’ resistivity beyond the MHD framework.

1.3

Kinetic Aspects Beyond the MHD Framework

As discussed in the previous sections, the ideal MHD condition breaks down in thevicinity of the X-line and the magnetic diﬀusion becomes dominant, forming the diﬀusion
region. In the MHD framework, the magnetic diﬀusion has attributed to the resistivity

CHAPTER 1. GENERAL INTRODUCTION

10generated by the classical Coulomb collisions (Spitzer, 1962) and it was found that thevalues were too small to drive fast reconnection in the collisionless plasma. Our next
step is to consider more general features of plasma beyond the MHD framework. In thissection we start with the generalized Ohm’s law in a steady state, which describes a
two-ﬂuid plasma,
E+V ×B =−me
11(V e · ∇)V e − ∇ · P e +J × B,enene(1.14)

where me is the electron mass, V e the electron ﬂow velocity, and P e the electron pressure
tensor. We assume that the classical resistivity term (ηJ) is negligible. If any term in
the right-hand side of (1.14) is negligibly small compared with the convection term in
the left-hand side, the ideal MHD condition holds. Thus, in order to determine when
the froze-in condition can break down, it is necessary to estimate the order of magnitudeof each term in the right-hand side.
The ﬁrst term expresses electron inertial eﬀects and the scale size is the order of theelectron inertial length λe . The scale length of the second term (pressure gradient) is theorder of the electron gyroradius λge . The last term represents the decoupling between
ions and electrons, so-called Hall eﬀects, and has a scale size of the ion inertial length
λi . When the scale length of gradient is much larger than any of these characteristic
lengths, the frozen-in condition holds.
In two-ﬂuid system consisting of ions and electrons, the inﬂow ions can be decoupledfrom the ambient magnetic ﬁeld when they approach the neutral sheet as close as λi .
Thus the region within λi is the diﬀusion region. Electrons, however, still remain magnetized at this distance, since their inertia is much smaller than those of ions. Instead,their frozen-in constraint can break down at the distance comparable with the local gyroradius λge , within which electrons are expected to perform the meandering motions(e.g., Pei et al., 2001). This non-magnetized region for electrons is called the electron
diﬀusion region. Therefore the diﬀusion region in the two-ﬂuid system develops twoscale structure as shown in Figure 1.5, diﬀerent from the MHD system (see Figure 1.1).
The electron diﬀusion region with the scale size of λe is found to be embedded within
the ion-scale diﬀusion region. Thus inside the diﬀusion region but outside the electron
diﬀusion region, ions are decoupled from the ambient magnetic ﬁeld but electrons arestill frozen-in to the ﬁeld and continue to move with the E × B drift motion toward
the current sheet. This relative motion between ions and electrons produces currents
in the vicinity of the X-line, forming the Hall current system (Sonnerup, 1979; Terasawa, 1983). Sonnerup (1979) ﬁrst predicted the presence of the Hall current systemtheoretically in which the four current loops resulted in a quadrupole structure in theCHAPTER 1. GENERAL INTRODUCTION

+QP

11

'NGEVTQP

KZ

E

VT
CTC

5GR

E

V1

δ

V2
l

Diffusion Region

Electron Diffusion Region

Figure 1.5: Two-ﬂuid picture around the X-line. The diﬀusion region develops twoscale structure associated with the ion and electron scales, forming the electron diﬀusion
region in side the ion-scale diﬀusion region. Red arrows represent the ion ﬂow vector,
while blue arrows the electron ﬂow vector.

out-of-plane (dawn-dusk) magnetic ﬁeld variations. In the past decade, not a few studies
with hybrid simulations (e.g., Mandt et al., 1994; Hesse and Winske, 1994, 1998; Lin and
Swift, 1996; Nakabayashi and Machida, 1997; Nakamura et al., 1998; Shay et al., 1998,
2001) and full particle simulations (e.g., Tanaka, 1995; Hoshino et al., 1998, 2001a,b;Hesse et al., 1999) have conﬁrmed the presence of the Hall current system and the associated quadrupole magnetic ﬁeld structure in the reconnection process. In these studies
electrons have been ascertained to play as carriers of the Hall currents. Furthermore,
Nagai et al. (2001) recently found the Hall current system (the accelerated electrons)
and its eﬀect on the magnetic ﬁeld in the Earth magnetotail, using magnetic ﬁeld andplasma measurements onboard the Geotail spacecraft.
In the electron diﬀusion region, on the other hand, the kinetic behavior of electrons
becomes important. The reconnection process of the magnetic ﬁeld actually proceeds inthis region. It has been revealed by using two-dimensional kinetic simulations that thedissipation process in the electron diﬀusion region is supported by the electron inertia
and the gradient of the electron pressure tensor in the generalized Ohm’s law (1.14)(Hesse and Winske, 1998; Kuznetsova et al., 1998; Hesse et al., 1999). This indicatesthat the collision eﬀect in a classical sense is replaced by electron inertial eﬀects in theCHAPTER 1. GENERAL INTRODUCTION

12collisionless reconnection. If their model is really appropriate as the dissipation model,
three-dimensional features such as the anomalous resistivity arising from interactions
between waves propagating along the current sheet and plasma particles would be negligible, and the reconnection process might be basically two-dimensional.
Though the dissipation process of magnetic ﬁeld is supported by the electron dynamics, it has little or no eﬀect on the reconnection rate (Biskamp et al., 1997; Shayand Drake, 1998; Birn et al., 2001). Instead ion dynamics can control the reconnection
process so the electron-scale structure, in which the actual reconnection of the magnetic
ﬁeld takes place. They suggest that the high reconnection rate is produced by the Hall
term in the generalized Ohm’s law (1.14), so that the inclusion of the Hall eﬀects should
be a suﬃcient condition to achieve fast reconnection. Though recent simulation studies
comparing MHD, Hall MHD, hybrid, and full particle simulations conﬁrmed the importance of the Hall eﬀect for fast reconnection, the system did not reach steady state andthe reconnection rate decreased after it reached a peak value (Shay et al., 2001; Hesseet al., 2001; Pritchett, 2001; Kuznetsova et al., 2001). They attributed the reason to
the periodic system in the direction of the initial magnetic ﬁeld, in which compressionaleﬀects in the magnetic islands could aﬀect the reconnection process. One of the solutions in order to suppress the boundary eﬀects is to perform the above comparison in alarge system so that the periodicity is negligible therein. However, such simulations withconventional full particle models are still hard to be carried out within the limitation ofcomputer resources. Furthermore, recent hybrid simulations demonstrated that the Halleﬀects were not needed to achieve fast reconnection (Karimabadi et al., 2004), although
their simulations did not include the electron kinetic eﬀects. Therefore it is necessary toreconsider what is supporting fast reconnection in collisionless plasma system.

1.4

Objective of the Thesis

Main problems on the physical processes associated with magnetic reconnection aresummarized in the following three ponts:
1. What trigger magnetic reconnection?
2. What supports fast reconnection in a steady state? Are the Hall eﬀects suﬃcient
for achieving large reconnection rate?
3. What is the energy transport process, that is, how and where can plasma be
accelerated and heated?

CHAPTER 1. GENERAL INTRODUCTION

13

Each of the above processes is expected to include strongly nonlinear coupling among
electron-, ion-, and MHD-scale processes. Thus numerical simulations in self-consistentsystems have been extensively examined to understand these subjects. The diﬃcultyin considering the problems associated with magnetic reconnection lies in the localityof micro-scale processes embedded in macro-scale processes. In order to make a selfconsistent description around the reconnection region, it is necessary to perform a largescale kinetic simulation including from electron kinetic processes to MHD processes.
However, it is still diﬃcult to carry out such simulations by using conventional kinetic
codes (e.g., particle-in-cell code, Vlasov code, and so on) because of the limitation ofcomputer resources.
In this thesis, we ﬁrst develop 2-1/2 dimensional full particle code (in which both ionsand electrons are treated as ﬁnite-size superparticles) with adaptive mesh reﬁnement(AMR) technique. The AMR technique subdivides and removes cells dynamically inaccordance with a reﬁnement criterion and it is quite eﬀective to achieve high-resolution
simulations of phenomena that locally include micro-scale processes. This new code
allows us to perform large-scale simulations of magnetic reconnection, which gives two
advantages on the modeling of magnetic reconnection. One is that it enables us todescribe the long time evolution of the current sheet without caring the eﬀects due tothe periodic boundary condition, which is related to the second problem listed above.
The other is that it is literally possible to pursue kinetic processes in an extensive region
around the X-line, which leads to the understanding of electron heating mechanism inthe plasma sheet-lobe boundary region (relating to the third problem).
In Chapter 2 we will give a detailed description of the new electromagnetic full particle
code with the AMR technique and emphasize the eﬃciency compared with conventional
particle codes. Chapter 3 is devoted to the long time evolution of the reconnection rate.
We will show that the reconnection rate can drop after it reaches a peak value. The Halleﬀects shown in the previous section seem to be still important until the reconnection
rate saturates. However, the meandering motion of the cold ions can lead to a structurechange of the electron diﬀusion region, causing the damping of the reconnection rate. A
generation mechanism of the heated electrons and the solitary waves, which have been
often observed in the Earth magnetotail, is explained in Chapter 4. Finally, in Chapter
5, we summarize this study and present a vision for future studies. Especially it isemphasized that developing the three-dimentional version of our code is inevitable inorder to understand the triggerring mechanism of magnetic reconnection.

CHAPTER 2Electromagnetic Full Particle CodeWith Adaptive Mesh Reﬁnement
Technique2.1

Introduction

Magnetic reconnection plays an important role in a fast conversion process of themagnetic energy to plasma kinetic and thermal energy. It occurs in the accretion disks ofastrophysical bodies, solar ﬂares, and the Earth magnetosphere, so that it is a ubiquitous
phenomenon in the universe.
In the Earth magnetosphere, magnetic reconnection facilitates the entry of particles
and energy from the solar wind into the magnetosphere at the day-side magnetopause.
It is also believed to aﬀect the dynamics of magnetospheric substorms by changing thecurrent sheet structure and the conﬁguration of magnetic ﬁeld lines in the Earth magnetotail. In a reconnection process, the ideal magnetohydrodynamic (MHD) condition
breaks down in the diﬀusion region arising near the X-line, where dissipation processesof the magnetic ﬁeld are dominant. In a suﬃciently collisional plasma, the resistive
MHD theory (e.g., Sweet, 1958; Parker, 1963; Petschek, 1964) is valid for describing thediﬀusion region by parameterizing collisional eﬀects. However in a region where the classical collision rate is very small such as in the magnetosphere, the Sweet-Parker diﬀusion
region is elongated and the reconnection rate is impractically low (Vasyliunas, 1975). Actually in such a collisionless plasma, the diﬀusion region develops a two-scale structureassociated with the electron and ion scales, that is, the electron inertial length λe and theion inertial length λi , respectively. Within a distance from the X-line of the order of λi ,the ions are easily unmagnetized because they have large gyroradii compared with thescale size λi . We call this region the ion diﬀusion region. On the other hand, even closer
to the X-line the electrons also decouple from the ambient magnetic ﬁeld. This region,

CHAPTER 2. FULL PARTICLE CODE WITH AMR

15which is called the electron diﬀusion region, has scale size of λe . Inside the electron diffusion region, where the electrons are unmagnetized and strongly accelerated away fromthe X-line, the nongyrotropic electron pressure can give rise to the reconnection electricﬁeld (Hesse and Winske, 1998), which implies that kinetic eﬀects of the electrons are
essential to dissipation mechanisms of the magnetic ﬁeld in the electron diﬀusion region.
Outside the electron diﬀusion region but inside the ion diﬀusion region, the ion motion
decouples from that of the still frozen-in electrons. This diﬀerence in motion produces
currents, so-called Hall current system, in the vicinity of the X-line (Sonnerup, 1979).
It is suggested that the inclusion of the Hall term realizes much faster reconnection rate
than that in the case of a resistive MHD (Birn et al., 2001, and references therein). Thusa full particle simulation, in which both electrons and ions are treated as particles, isvery eﬀective to describe the dynamics in the vicinity of the X-line.
Most full particle codes employ particle-in-cell (PIC) technique, in which the current
and charge densities are accumulated on a spatial grid from the particle data and theforces on the particles are interpolated from the adjacent spatial grids. Though thefull particle model using PIC technique is conceptually very simple and one does not
introduce any approximation in the basic laws of mechanics and electromagnetism, therestrictions on the grid spacing are very severe and one must set it to be as small as theelectron Debye length λDe , which is the order of 103 m in the central current sheet ofthe magnetotail. On the other hand, magnetic reconnection changes the current sheet
structure dynamically in MHD scale of the order of 109 m. Hence, for modelling of themagnetosphere including the nonlinear evolution of magnetic reconnection, the spatial
dynamic range is required to be 106 . Furthermore, since the electromagnetic PIC code
requires at least dozens of particles per cell to suppress a numerical noise, the numberof particles for each species is over 1013 even in the two dimensional system. However asimulation with 1012 cells and 1013 particles in the two dimensional system is unrealistic
today due to the limitation of computer resources. Usual compromise is to use low mass
ratio of the ion to the electron, although too low mass ratio allows a spurious wave modeto grow (Daughton, 1999). However, even when the mass ratio is 100 instead of 1834 thatis realistic, the required numbers of cells and particles are 1010 and 1011 , respectively. It
is still diﬃcult to conduct such a simulation with currently available computer resources.
In the Earth magnetotail, the number density of plasma is conﬁrmed to be muchhigher in the central plasma sheet than in the lobe region according to the observations
with satellites (e.g., Baumjohann et al., 1989; Baumjohann and Treumann, 1997), sothat the electron Debye length is the smallest in the central plasma sheet. Therefore,
the spatial resolution required for the numerical stability is much higher in the central

CHAPTER 2. FULL PARTICLE CODE WITH AMR

16

plasma sheet than in the lobe region. Since the diﬀusion region appears also in thecentral plasma sheet, the spatial resolution should be high there also in terms of thephysical perspective. In conventional full particle simulations of magnetic reconnection(e.g., Pritchett, 1994; Tanaka, 1995; Hesse and Winske, 1998; Shay et al., 2001; ?), thespatial resolution has been uniform in the whole simulation area in accordance with that
required in the central plasma sheet. Thus, in the lobe region which occupies most of thesimulation area, the spatial resolution has been unnecessarily high. The adaptive meshreﬁnement (AMR) technique is beneﬁcial to reducing the number of cells and saving
computer resources by introducing ﬁner cells hierarchically in only the regions, suchas the central plasma sheet and the diﬀusion region, where higher spatial resolution isrequired (see § 2.2.1).
A block structured adaptive method was extensively developed by Berger and the
coworkers to solve partial diﬀerential equations on a hierarchy of nested cells coveringhigh-resolution regions (Berger and Oliger, 1984; Berger and Colella, 1989). An introduction of the AMR technique to a particle code has been examined for the application to
the N-body simulations describing cosmological matter such as the formation and structure of galaxies (Villumsen, 1989; Jessop et al., 1994; Suisalu and Saar, 1995; Gelato etal., 1997; Kravtsov et al., 1997; Yahagi and Yoshii, 2001; Knebe et al., 2001; Teyssier,2002). In the N-body codes, the mass density is assigned to the neighboring grids using
the cloud-in-cell (CIC) algorithm, which is equivalent to the PIC algorithm for orthogonal coordinates (Birdsall and Langdon, 1991). In the early N-body codes with the
AMR technique (Villumsen, 1989; Jessop et al., 1994; Suisalu and Saar, 1995; Gelato
et al., 1997), they adopted nested cells as like as Berger and Oliger (1984) to increase
force resolution, so that the reﬁned regions were give in rectangular shape. Though therectangular reﬁnements facilitate the data structure, they are not suitable for models
including the complicated structures that are diﬃcult to cover eﬃciently with rectangular regions. On the other hand, a tree structured cells have been also used in order
to increase the spatial resolution locally and achieve a high dynamical range simulation (Barnes and Hut, 1986; Griebel and Zumbusch, 1999; Dreher and Grauer, 2005).
Kravtsov et al. (1997) developed the adaptive reﬁnement tree (ART) code, in which thetree structured cells are connected with each other on all hierarchical levels. The ART
code subdivides all cells that satisfy a reﬁnement criterion regardless of the shape ofthe reﬁned regions, so that it enables us to achieve an eﬃcient reﬁnement to cover acomplicate structure. Although their code is adaptive only in space,recently developed
codes (Yahagi and Yoshii, 2001; Knebe et al., 2001; Teyssier, 2002) are adaptive not
only in space but also in time. Plasma simulations using an electrostatic PIC code with

CHAPTER 2. FULL PARTICLE CODE WITH AMR

17

the AMR have been recently examined in order to describe the ion beam transport in aheavy ion fusion (Vay et al., 2004).
In this paper we present a newly developed two and a half dimensional electromagnetic full particle code with the AMR technique. In this code we deal with two species ofcharged particles, that is, ions (protons) and electrons. The current and charge densities
are assigned onto the hierarchical cell corners using the PIC algorithm, which process
is similar to the mass density assignment by the CIC algorithm used in the N-bodycodes. The electromagnetic ﬁeld is calculated via Maxwell’s equations, unlike the Poisson equation used to derive the gravitational ﬁeld in the N-body codes. The particles(superparticles) are also subdivided in the reﬁned regions in order to suppress the numerical noise that arises due to the decrease in the number of the particles per cell. In
the following sections we describe the details of our code (§ 2.2), and discuss the resultsof test simulations and application to the current sheet evolution (§ 2.3). Finally wesummarize and conclude this work in § 2.4.

2.2

Description of the Code2.2.1 Data Structure and Cell Generator
In order to increase a spatial resolution eﬃciently, the AMR technique subdivides
only cells that satisfy a reﬁnement criterion and add data sets for ﬁner cells hierarchically
onto the uniform base cells that cover the entire simulation area. If a base cell is reﬁned,
it has four child cells that have half the size of the base cell. These child cells can be also
reﬁned in turn and ﬁner cells are generated, and so on. In our code, a reﬁnement level
L in the hierarchy is deﬁned by using the cell size of the level (∆L ) as L ≡ log2 (ly /∆L ),where ly is the vertical size of the two-dimentional simulation area. We use only cellswith integer level.
The data structure in our approach is completely diﬀerent from that used in theconventional electromagnetic particle codes. Cells are treated as independent units organized in reﬁnement trees rather than elements of arrays, so that one can build a very
ﬂexible cell hierarchy that can be easily modiﬁed. The hierarchical cell structure in ourcode is basically the same as the fully threaded tree structure developed by Khokhlov(1998), and supported by a set of pointers as shown in Figure 2.1. Each cell is needed tohave information of the parent, child, and neighboring cells. Since a cell has four childcells, if any, we group these four cells together in order to save memory for pointers andfacilitate parallel computing. We call this group an oct (after an octet in the case of a
cubic cell). Each cell has a pointer to the child oct (OctCh). Each oct has pointers toCHAPTER 2. FULL PARTICLE CODE WITH AMR

OctLevel L

Oct

Cell

Cell

iNb

OctCh
iPr

Cell

Cell

OctCh
iPr

iNb

iNb

Cell Cell Cell Cell

Level L+1

Oct
inp

Particles

18

Oct
inp

inp

inp

pnb

Figure 2.1: Hierarchical data structure in our code. One dimensional structure is shownfor simplicity. The pointers that each cell or oct has are basically the same as those in
the tree structure developed by Khokhlov (1998), but we additionally give two pointers,
inp directed to a representative superparticle from each oct, and pnb directed from eachparticle to the neighboring particle in the same oct.

the parent cell (iPr ), the parent cells of the neighboring octs (iNb), and a representative
particle belonging to the oct (inp). Note that each particle belongs to the ﬁnest oct inwhich the particle is located. Each particle has a pointer to the neighboring particle inthe same oct (pnb), which leads to a beaded structure of particles as described in Figure
2.1 (Teyssier, 2002). In addition to these pointers, each oct has information of the levelof cells included therein, the integral position at a corner of the oct, and the numberof particles in the oct. The level of a cell is very useful information to eﬃciently sweep
through a level. The position of a oct and the number of particles included in the oct
are helpful to construct the beaded structure of particles. Each cell also has physicaldata at a corner of the cell; charge density, current density, and electromagnetic ﬁelds.
A reﬁnement process subdivides a cell that satisﬁes a reﬁnement criterion and generates four child cells. In electromagnetic particle codes, one of the reﬁnement criteria
should be deﬁned by the local electron Debye length in terms of the numerical stability.
However, additional criteria can be introduced in accordance with physics that is incorporated in the code. An example in the case of nonlinear evolution of the plasma sheet isshown in § 2.3.3. When the parent cell is reﬁned, we check if there exist the surrounding
eight cells having the same level as the parent cell. Even if one of eight cells are absent,no splitting is implemented and no child cells are created. Thus, in our tree structure,
no neighboring cells have level diﬀerence greater than one, which assures that there are

CHAPTER 2. FULL PARTICLE CODE WITH AMR

19no sharp gradients on the boundaries of the reﬁned regions.
The reﬁnement procedure is similar to that in Kravtsov et al. (1997). The process
is started from the base level cells and recursively proceeds to reach the dynamic range
level cells at each time step. First, we put ﬂags to cells that satisfy reﬁnement criteria.
At this point, no cells are subdivided. In order to smooth the marked region, we next put
ﬂags to the eight cells surrounding each cell ﬂagged ﬁrstly. Furthermore, we construct
buﬀer region surrounding the reﬁned region by giving other ﬂags to cells that surroundthe originally ﬂagged cells but have no ﬂags for the reﬁnement. The buﬀer regions areused to set boundary conditions for calculating the electromagnetic ﬁelds in the reﬁned
region as mentioned in § 2.2.4. After that, we search again the whole cells with thetarget level. If a cell has a ﬂag but has no child cells, the cell is subdivided and half-sized
reﬁned or buﬀer cells are generated. If a cell has no ﬂag but has child cells, the childcells are removed from the hierarchical structure. If a cell has both a ﬂag and childcells, or has neither a ﬂag nor child cells, no operation is conducted for the cell. Thus
the reﬁnement operation is performed only to the cells that have any changes in termsof the reﬁnement criteria at each time step, so that the reﬁnement procedure is quickly
completed.
The beaded structures of particles must be rearranged in association with the reﬁnement procedure and the movement of particles. When a new oct is created, the particles
located in the parent cell are removed from the beaded structure in the parent oct andreconnected to that in the child oct. On the other hand, when an oct is removed fromthe tree structure, the beaded structure is reconnected to that in the parent oct. Furthermore, particles moving away across oct boundaries are removed from the beaded
structure in the oct and reconnected to those in new octs in which they are located.

2.2.2 Particle Splitting and Coalescence
One of the main problems in developing the electromagnetic PIC codes using the
AMR technique is decrease in the number of particles per cell in the reﬁned regions(Fujimoto and Machida, 2005a). Especially, in the vicinity of the X-line formed in association with magnetic reconnection, spatially high-resolution simulations are required
because the kinetic eﬀects of electrons are expected to be important. However, the number density of plasma in such a region is low due to inﬂow of the tenuous plasma in thelobe region, so that the number of particles per cell is decreased and the numerical noise
is increased.
In order to solve the problem, we subdivide particles (superparticles) residing in thesubdivided cells and control the number of particles per cell. We use the particle splitting-CHAPTER 2. FULL PARTICLE CODE WITH AMR

20

coalescence algorithm developed by Lapenta (2002). The processes should conserve thefollowing quantities,
1. the moments on each grid, that is, the charge and current densities,
2. the total charge and mass of particles,
3. the total momentum and energy of particles,
4. the distribution function of particles.
In the splitting algorithm in the two-dimensional system, a particle with charge q0 , mass
m0 , position x0 , and velocity v 0 can be replaced by four particles. The new particles
labeled by p (p = 1, 2, 3, 4) have charges qp = q0 /4, masses mp = m0 /4, and velocitiesv p = v 0 . The new particles are located at x1,2 = x0 ± ∆r, x3,4 = x0 , y1,2 = y0 ,y3,4 = y0 ± ∆r (see Figure 2.2a). Here, ∆r should be chosen so as to suppress the local
ﬂuctuation of the spatial particle distribution, which can arise from the particle splitting,as much as possible. In this context, ∆r = ∆L / Np is considered to be a natural choicein the two-dimensional case, where Np is the number of particles in the cell with the size
∆L . This algorithm exactly conserve the above mentioned four quantities, if the newparticles are all located in the same cell as the old particle. If one of the new particles
is expected to be placed in a neighboring cell, the splitting is not implemented.
We also conduct the coalescence of particles. This algorithm is also simple. In ourcode, we choose two particles labeled by 1 and 2 that are located in the same cell and havesame charges q1 = q2 = q (that is, same masses m1 = m2 = m), and proximate velocities
v 1 and v 2 which satisfy |v 2 − v 1 | < αvth , where vth is the thermal velocity deﬁned by
vth ≡ 2T /m (T is temperature), and α is a small value which is set as α = 0.1 inour code unless otherwise mentioned. The new particle labeled 0 has charge q0 = 2q,
mass m0 = 2m, position x0 = (x1 + x2 )/2, and velocity v 0 = (v 1 + v 2 )/2 (see Figure2.2b). However, the coalescence process does not exactly conserve the total energy ofparticles and the distribution function (Lapenta, 2002). The diﬀerence in kinetic energyof particles between before and after the coalescence is found to be under 25α2 percent of2, so that the error that arises during one coalescence is within 0.5% of the thermalmvthenergy of particles with mass m in our code (α = 0.1). Nonetheless, we hardly conduct
the particle coalescence (usually once in a few hundred time steps), and never choose theparticles in the most reﬁned cells to avoid the numerical errors in physically importantregions.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

21

!


0



0

"

(a)

(b)

Figure 2.2: Particle splitting and coalescence method in the two-dimensional system.
The splitting process in (a) subdivides a particle labeled 0 into four particles numbered
1 to 4, while the coalescence process in (b) joins two particles numberd 1 and 2 andgenerates a new particle labeled 0 at the midpoint.

2.2.3 Particle and Force Weighting within the Hierarchical Cell
We use the standard PIC algorithm to accumulate the charge and current densities
onto each cell corner and interpolate the force on a cell corner to a particle. The charge
and current that a particle located in a cell carries are assigned to the four cell corners
in the following forms,
ρnl,m =
qsj SL (xnsj − X l,m ),s
n+1/2
j l,m =

(2.1)

j


s

j

n+1
n
n+1/2 SL (xsj − X l,m ) + SL (xsj − X l,m ),qsj v sj2

(2.2)

where qsj , xnsj , and v nsj are the charge, position, and velocity of j-th particle of speciess (ion and electron) at time n∆t (∆t is a time interval and equivalent for all particles),respectively, and X l,m represents the spatial position (l∆L , m∆L ) of a cell corner. SL (x)is the shape function and deﬁned by using the cell size ∆L on each reﬁnement level L in
the form,⎧⎨ (∆L − |x|)(∆L − |y|)if |x| ≤ ∆L and |y| ≤ ∆L
(2.3)
SL (x, y) =∆2L⎩0otherwise.
Each particle is belonging to the ﬁnest oct including the particle, which means that itis only the oct that has a direct connection to the particle without passing through theother octs. However, we solve Maxwell’s equations and the Poisson equation on every

CHAPTER 2. FULL PARTICLE CODE WITH AMR

22

reﬁnement level independently as described in § 2.2.4, so that the charge and current
densities must be assigned onto every cell regardless of whether or not the oct including
the cell has a pointer to a particle. This assignment process using equations (2.1) and(2.2) is very heavy and time consuming. So we do not adopt this manner for all thecells. Alternatively this process is carried out only for the cells belonging to the octs
that have direct pointers to the particles locating in the cells, and the following operationis conducted for the other cells,
ALl,m =

1
{4AL+12l,2m
16 L+1
L+1
L+1+2 (A2l−1,2m + AL+1
2l,2m−1 + A2l+1,2m + A2l,2m+1 )
L+1
L+1
L+1
+ (AL+1
2l−1,2m−1 + A2l−1,2m+1 + A2l+1,2m−1 + A2l+1,2m+1 )},(2.4)

where ALl,m is a physical quantity (charge density and current density in the currentcase) at the spatial position of (l∆L , m∆L ) on the level L cells. The charge and current
densities calculated by the use of equation (2.4) is found to be identical to those derived
by equations (2.1) and (2.2), if they are known on the ﬁner cells. This procedure is very
fast compared with equations (2.1) and (2.2).
The electromagnetic ﬁelds at the positions of the level L particles are interpolated
from the cell corners on the reﬁnement level L in order to calculate the forces that works
on the particles. The interpolation is also implemented by using the PIC algorithm:
E L (xsj ) =lB (xsj ) =
L

(2.5)

B Ll,m SL (xsj − X l,m ),

(2.6)

m


l

E Ll,m SL (xsj − X l,m ),mwhere E Ll,m and B Ll,m are the electric and magnetic ﬁelds at the position of (l∆L , m∆L )on the reﬁnement level L, respectively.

2.2.4 Integration of the Field Equations
We use the complete set of Maxwell’s equations to derive the electromagnetic ﬁelds.
The time advancements of the electric ﬁeld E and magnetic ﬁeld B are described by
Ampère’s law and Faraday’s law, respectively, that is,1∂E= c2 ∇ × B − j,∂t
ε0∂B= −∇ × E,
∂t

(2.7)
(2.8)

CHAPTER 2. FULL PARTICLE CODE WITH AMR

23

where c is the velocity of light, and ε0 is the dielectric constant in vacuum. The time
diﬀerences of equations (2.7) and (2.8) are,∼ n+1/2
E n+1 − E n1= c2 ∇ × B− j n+1/2 ,
∆t
ε0
n+1
n
n+1/2
∼
B
−B
= −∇ × E
,
∆t

(2.9)
(2.10)

where,
E n + E n+1,
2∼ n+1/2
B n + B n+1
B=.
2∼ n+1/2E=(2.11)(2.12)

Here, the superscripts denote the time levels.
The initial ﬁelds E 0 and B 0 are determined self-consistently with the spatial distributions of the initial charge and current densities. To advance the ﬁeld data, oneeliminates B n+1 in equation (2.9) by the use of equations (2.10), (2.11), and (2.12), andobtains,
c2 ∆t2c2 ∆t2n+1∇ × ∇× E∇ × ∇× E n
= 1−
1+
4
4
∆t n+1/2+ c2 ∆t ∇ × B n −j.
ε0(2.13)The right-hand side consists of the known ﬁelds and currents, while the left-hand side has
the elliptic operator. Since the boundaries of the reﬁned regions can have an arbitrary
shape in our algorithm, methods to solve the equation are restricted and we choose theconjugate gradient method which is one of the relaxation methods (Press et al., 1992).
If E n+1 is obtained from equation (2.13), then B n+1 is calculated from equations (2.10)
and (2.11).
However, the assignments of the charge and current densities by using equations (2.1)
and (2.2), or equivalently equation (2.4), may violate the continuity equation,∂ρ+ ∇ · j = 0.
∂t(2.14)

If this equation is not satisﬁed, the straightforward integration of equations (2.7) and(2.8) will lead to the evolution of nonphysical ﬁelds that do not satisfy ∇ · E = ρ/ε0 .
Thus we need to make a correction to ensure that ∇ · E = ρ/ε0 is maintained (Langdon
and Lasinski, 1976). If E is a corrected electric ﬁeld and E is that computed fromequation (2.7), we assume E = E − ∇δφ. Then δφ is calculated by solving the Poisson

CHAPTER 2. FULL PARTICLE CODE WITH AMR

24equation,
∆2L 2 ρ−∇ δφ = 1 +− ∇ · E.
∇6
ε02(2.15)

The ﬁrst term of the right-hand side is worth noticing. The second term in the parentheses arises from the truncation error of the second-order central diﬀerences of ∇ · E.
∇ · E in the diﬀerence equation has the truncation error within the second order in space(∼ o(∆2L )), while the replacement of ∇ · E = ρ/ε0 has no truncation so that no errors
are included. This diﬀerence in truncation error can yield nonphysical electrostatic ﬁeld
because the error can become the source in the Poisson equation. If the diﬀerence formof ∇ · E in one-dimensional system is given by (Ex,i+1 − Ex,i−1 )/2∆L (the subscripts i + 1
and i − 1 denote the positions of cell corners), this term is approximated in the following
form in the spatial fourth-order accuracy,
Ex,i+1 − Ex,i−1∂Ex 
∆2L ∂ 3 Ex +.
(2.16)
2∆L∂x i
6 ∂x3 i
The second term of the right-hand side of the equation expresses the truncation error
within the spatial second order. We add this term explicitly. Since ∂Ex /∂x = ρ/ε0
is exactly satisﬁed, the right-hand side of this equation becomes the ﬁrst term of theright-hand side of equation (2.15). In two- or three-dimensional system, equation (2.15)
is also satisﬁed within the second-order accuracy. The Poisson equation is solved byusing the conjugate gradient method again.
In our code, the electromagnetic ﬁelds are ﬁrst calculated from equations (2.13) and(2.10) on the base level (LB ) cells, because the boundary conditions at the edge of thesimulation area are assumed to be given. The solutions are interpolated onto the buﬀer
cell corners on the next level cells as the boundary conditions as shown in Figure 2.3(Yahagi and Yoshii, 2001). The interpolations are simply performed in the linear form,L
AL+1
2l,2m = Al,m ,
1 L
L
AL+1
2l+1,2m = (Al,m + Al+1,m ),2
1
L+1
A2l,2m+1 = (ALl,m + ALl,m+1 ),2
1
L+1
A2l+1,2m+1 = (ALl,m + ALl+1,m + ALl,m+1 + ALl+1,m+1 ).
4(2.17)

Then the electromagnetic ﬁelds are solved on this level (LB + 1) cells by the use of equations (2.13) and (2.10). These solutions are also interpolated to the buﬀer cell corners on
the next level cells, and the procedure is recursively carried out until the solutions on the
dynamic range level (LD ) cells are obtained. After then, the ﬁeld data on the LD cells

CHAPTER 2. FULL PARTICLE CODE WITH AMR

25

Level L

Level L+1

Buffer cells

Refined cells

Interpolation fromcoarse cells tofine cells

Projection fromfine cells tocoarse cells

Figure 2.3: Interactions of the ﬁeld data between the levels L and L + 1. The solutions
on the level L cells are interpolated onto the buﬀer cell corners on the level L + 1 cells
as the boundary conditions. The electromagnetic ﬁelds solved on the level L + 1 cells
are in turn projected to the level L cells.

are projected onto the LD − 1 cells, and the process are also recursively proceeded until
the data on the LB cells are replaced (see Figure 2.3). The projections are performedby using equation (2.4) which is found to eliminate the aliasing. Finally, the Poissonsolver described in equation (2.15) must be executed to maintain ∇ · E = ρ/ε0 on eachlevel cells. The electromagnetic ﬁelds obtained by this procedure without recalculating
under the corrected boundary conditions, however, do not give proper solutions on eachreﬁnement level, because they are not continuously diﬀerentiable at the boundaries of
the reﬁned regions (Huang and Greengard, 2000; Vay et al., 2004). Nevertheless, thetest simulations described in § 2.3 indicate that numerical errors in association with thecell reﬁnement are almost negligible, so that we consider that the inﬂuence of the inconsistency in the ﬁeld solutions between reﬁned regions is not signiﬁcant. More accurate
treatment of the elliptic equations will be examined in the future work.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

26

2.2.5 Integration of the Equations of Motion
The initial particle distribution in velocity space is chosen to be Maxwellian orshifted-Maxwellian. The Maxwell’s distribution is obtained using the Box-Muller method(Press et al., 1992), in which particle velocities with a thermal velocity vth are describedin the form,
v1 = vth −2 ln x1 cos 2πx2 ,(2.18)
v2 = vth −2 ln x1 sin 2πx2 ,where x1 and x2 are the uniform random numbers. The velocities v1 and v2 are independent of each other.
The particle equations of motion to be integrated are,
dv sj
qsj
(E(xsj ) + v sj × B(xsj )),
=
dt
msj
dxsj
= v sj ,dt(2.19)(2.20)

where, msj is the mass of j-th particle of species s. E and B must be those on the cells
whose level is the same as the particle level. The time diﬀerences of equations (2.19)
and (2.20) are,
n+1/2

v sj

n−1/2

− v sj
∆t

qsj
=
msj

n+1/2

E

n

(xnsj ) +

v sj

n−1/2

+ v sj2

× B n (xnsj ) ,

n
xn+1
n+1/2
sj − xsj= v sj .
∆tAbove equations represent a standard second-order leapfrog integration scheme. We use
n+1/2
the Buneman-Boris method to calculate v sj , in which the electric and magnetic forces
are separately treated (Birdsall and Langdon, 1991).

2.3

Test Simulations

2.3.1 Landau Damping of the Langmuir Waves
Our code is ﬁrstly checked against the Landau damping of the Langmuir waves that
propagate across the reﬁned regions. The initial setting of the Langmuir waves in thetwo-dimensional system is given as one-dimensional plane waves with an appropriate
dispersion relation, which are assumed to propagate in the x direction (the horizontal
direction in the simulation area) and to be uniform in the y direction. The boundariesof the system are periodic in both x and y directions. The number of the base level

CHAPTER 2. FULL PARTICLE CODE WITH AMR

27cells is Nx × Ny = 64 × 64, unless otherwise mentioned, which means LB = 6. In
this and next subsections, we assume the number of reﬁned layers is only one so thatLD = 7 in order to avoid considerable numerical noise due to small number of particles
per cell in the case that we do not conduct the particle splitting for comparison. In
the run including the AMR, we initially subdivide the base level cells that satisfy thecondition, |Xj −lx /2| < 16∆LB , where Xj is the central x-coordinate of the j-th cell, andlx represents the horizontal size of the simulation area. The reﬁned region is assumed
to be ﬁxed to the base level region and not changed, even though a wave ridge moves
from one place to another as time goes on. Thus we have an interest only in whether
the waves propagate properly across the boundaries from the base region to the reﬁned
region and vice versa, and the hierarchical cells with the particle splitting-coalescencealgorithm are adequately applied to the conventional PIC codes. In the runs including
the AMR and particle splitting algorithms, the particles that are located in or enter thereﬁned region can be split into four particles, and the child particles can be carefully
coalesced if they move outside the region. Other initial parameters are set as mi /me =104 , Ti /Te = 1.0, ωce /ωpe = 1.0, and c/vth,e = 102 , where ms and Ts represent themass and temperature of species s, ωce and ωpe are the cyclotron frequency and plasma
frequency of electrons respectively, and vth,e is the electron thermal velocity deﬁned asvth,e ≡ 2Te /me . Adoption of an unrealistic high mass ratio indicates that we have aninterest only in the electron scale processes, but not those of the ion scale. The size ofthe base level cells is set as ∆LB = 1.0 × 10−2λe , where λe is the electron inertial length.
The time step is chosen ∆tωpe = 5.0 × 10−3 in the whole region. The initial number ofparticles is 64 per base level cell, so that the total number is approximately 2.6 × 105 .
Figure 2.4 shows the time evolutions of the ﬁeld energy of the wave component Ex
normalized by the initial values. We have examined three cases of λ = 8∆LB , 12∆LB ,and 16∆LB , where λ is the wavelength. In the runs of λ = 12∆LB , the number ofcells in the base level is set as Nx × Ny = 96 × 64 that is compatible with the cyclic
boundary conditions in this case. Energy is averaged over the base level cells. The
solid, dashed, and dotted lines denote the runs without the AMR, with the AMR but
not including the particle splitting, and including both the AMR and particle splitting,respectively. The thick lines represent the theoretical proﬁles of the growth rate de2/2n0 k v(∂f0 /∂v)v=ω/k , which is deduced from the linear analysis ofscribed as γ = πωpe
the Landau resonance. Here, n0 is the plasma number density averaged over the entire

simulation area, k is the wavenumber, ω is the angular frequency, and f0 is the background part of the electron distribution function assumed to be Maxwellian. We ﬁndthat each simulation result provides a good ﬁt to the theoretical proﬁle at the stage ofCHAPTER 2. FULL PARTICLE CODE WITH AMR

Linear theory
Without AMRWith AMR
With AMR &Particle splitting

Field Energy
1
0.5

n0Ex2/2c228

λ=16∆LB

0.2
0.1

λ=12∆LB

0.05

λ=8∆LB

0.02
0.01

0

10

20

tωpe

30

40

Figure 2.4: Time evolution of the ﬁeld energy of the Langmuir wave Ex , normalized
by the initial values. We examined three cases of wavelength λ = 8∆LB , 12∆LB , and
16∆LB . These energies are averaged over the base level cells. The thick lines denotethe theoretical proﬁles deduced from the linear analysis of the Landau resonance. The
solid, dashed, and dotted lines represent the runs without the AMR, with the AMR but
not including the particle splitting, and including both the AMR and particle splitting
algorithms, respectively.

the linear evolution. We also ﬁnd that the results of ’With AMR’ and ’With AMR &Particle splitting’ runs are in good agreement with ’Without AMR’ run for each wavelength. We can see a little discrepancy at the nonlinear stage of the case λ = 8∆LB .
This is thought to arise from the diﬀerence in noise distribution due to thermal plasmaand can be aﬀected by the short wave reﬂections at the boundaries of the reﬁned region.
However, the discrepancy is not essential because physical phenomena we try to understand using PIC codes should have much larger amplitudes compared with the noiselevel. The eﬀects of the wave reﬂections at the boundaries are described in detail in thenext subsection.

2.3.2 Wave Reﬂections at the Boundaries of the Reﬁned Regions
One of the issues when we combine the PIC and AMR techniques is wave reﬂections
at the boundaries surrounding the reﬁned regions. Wave dispersion relations in theframework using ﬁnite-sized particles are diﬀerent from those in the system using theactual point particles, and dependent on the cloud sizes, that is, the grid spacings in ourcode (Birdsall and Langdon, 1991). Thus if a wave in a coarse cell region propagates
into a reﬁned region or vice versa, the inconsistency in the wave dispersions between

CHAPTER 2. FULL PARTICLE CODE WITH AMR

29"With AMR & Particle splitting"– "Without AMR"

With AMR & Particle splitting

5

5

(b)

(a)
4

3.

3

3

2

2

1

1

0
-150 -100 -50

0.
0

50

100 150

3. 10-2

40-150 -100 -50

050

100 150-3. 10-2

Figure 2.5: (a) Wave spectral density of Ex in the ω − k space for the ’With AMR &Particle splitting’ run in the case of λ = 32∆LB , in which the Landau resonance between
the Langmuir wave and the background electrons is negligible. (b) Diﬀerence in wave
spectral density between the runs of ’With AMR & Particle splitting’ and ’Without
AMR’. The wave lengths are again λ = 32∆LB .

the hierarchical cell layers can cause wave reﬂections at the boundary. In order to check
whether the reﬂections are signiﬁcant, we test again the Langmuir waves that propagate
across the boundaries. The initial settings of the hierarchical cells and plasma condition
are the same as those in the previous subsection. In order to avoid wave damping due
to the Landau resonance, we examine the relatively long wavelength, λ = 32∆LB , ofwhich the phase velocity is much larger than the electron thermal velocity, so that theresonance between the Langmuir waves and background electrons is very weak.
In Figure 2.5(a), the wave spectral density of Ex in the ω − k space is shown for
the run including the AMR and particle splitting processes. Strong peak that arises at(k, ω)  (20λ−1
e , ωpe ) represents the normal Langmuir wave initially loaded. If the wave
reﬂects at the boundaries of the reﬁned region, the reﬂected wave should have strong
peak at (k, ω)  (−20λ−1
e , ωpe ). However, such a peak is not shown in Figure 2.5(a).
The diﬀerence in wave spectral density between the ’With AMR & Particle splitting’
and ’Without AMR’ runs is described in Figure 2.5(b). This diﬀerence is thought torepresent any wave modulation caused by the presence of the reﬁned cells. However, thescale of colour contour in Figure 2.5(b) is two order smaller than that in Figure 2.5(a),which indicates that the modulations including the wave reﬂections at the boundaries
are very small, if any, compared with the source wave amplitude and can be masked by
the thermal ﬂuctuations.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

30

2.3.3 Nonlinear Evolution of the Plasma Sheet
In this subsection, we show the electron dynamics that evolves nonlinearly from theinitial plasma sheet. We use the Harris-type current sheet (Harris, 1962) as an initialcondition, which is given as Bx (z) = −B0 tanh(z/λ), where λ is the half width of theinitial plasma sheet. We set as λ = 0.5λi0 , where λi0 is the ion inertial length deﬁned
by the initial plasma sheet density, nps . In addition to the equilibrium, we put a small
perturbation of the form,
Bxp (x, z) = −2δ/λ cos(2πx/lx ) sech2 (z/λ) tanh(z/λ)
Bzp (x, z) = 2πδ/lx sin(2πx/lx ) sech2 (z/λ),(2.21)

where δ is given as δ = 0.03B0 λi0 in the current case, so that the initial value of thelreconnected magnetic ﬁeld deﬁned by ψ ≡ 0 x |Bz (x, z = 0)|dx is only 0.12 B0λi0 , wherelx is the system size in the x direction. Furthermore, we add the background plasma suchas nb = nb0 tanh2 (z/λ) in order to describe the lobe plasma in the Earth magnetotail.
Here, nb0 is the asymptotic lobe density and nb0 = 0.044nps is loaded in the current
runs. We do not give any perturbation in the initial current density consistent with
the initial ﬁeld perturbation. However, once the simulation starts, the initial current
proﬁle is quickly modiﬁed to adjust the ﬁeld perturbation. As a result, a thinner current
sheet is formed around the center of the simulation area, so that the tearing instability
selectively develops therein. The system size is lx × lz = 15.4λi0 × 15.4λi0, we assume
the periodic boundary in the x direction and the conducting wall in the z direction, thatis, the normal component of the electric ﬁeld (Ez ) and the tangential components of themagnetic ﬁeld (Bx , By ) has no gradients and the charge and current densities are vanishedat the z boundaries. The initial plasma condition is mi /me = 100, Ti,ps /Te,ps = 8.0,
Ti,lobe /Te,lobe = 1.0, Te,lobe /Te,ps = 1.0, and c/vth,e = 5.0, where Ts,ps and Ts,lobe are thetemperatures of the species s at the central plasma sheet and the lobe region, respectively.
Time step is ∆tωci = 8.0 × 10−4 for all particles and reﬁned regions in order to satisfy
the Courant condition on the level LD cells, where ωci = eB0 /mi is the ion cyclotron
frequency in the asymptotic lobe region. The cell shape is assumed to be square in thecurrent runs.
We examine three runs for comparison: ’Without AMR’, ’With AMR’, and ’With
AMR & Particle splitting’ runs. Table I shows information on these three runs, showing
whether the AMR or particle splitting algorithms are used, the base level (LB ), thedynamic range level (LD ), the number of cells, and the number of particles. In the’Without AMR’ run, the whole simulation area is initially covered only by the ﬁne cellswith level L = 10 and the number of particles is large enough to suppress the numericalnoise even in the ﬁne cells. The numbers of cells and particles are not changed during

CHAPTER 2. FULL PARTICLE CODE WITH AMR

RunWithout AMRWith AMR
With AMR &Particle splittingaTable I: Simulation information.
Total number
ParticleAMR
LB LD
of cells
splittingNoNo
10 10 1.0 × 106
YesNo7
10 8.4 × 104 a
Yes

Yes

7

10

9.1 × 104 a

31

Total numberof particles
2.2 × 107
2.2 × 107
4.8 × 106 a

Averaged value over the run.
the run, so that we assume results of the ’Without AMR’ run are identical with those inconventional PIC simulations. In the ’With AMR’ run, the AMR processes are included
but the particle splitting is not implemented. Thus the number of cells is greatly reduced
according to the reﬁnement condition described later, but the number of particles is thesame as the ’Without AMR’ run. Finally, the ’With AMR & Particle splitting’ run
includes both the AMR and particle splitting processes, so that the numbers of cells
and particles are reduced and changed dynamically in association with the evolution of
the plasma sheet. The reﬁnement condition is deﬁned by three physical values. The√
ﬁrst is the local electron Debye length, λDe = vth,e / 2ωpe , which is required to avoid anumerical heating of local plasma, where ωpe is the electron plasma frequency. We use theinitial value of vth,e for calculation of λDe , because electrons are expected to be heated in
the region where the electron dynamics is important, that is, we use the minimum valueof λDe . The second is the out-of-plane electron ﬂow velocity (Vey ), because the inertial
term in the generalized Ohm’s law gets dominant in the region where the electron ﬂow isstrong, so that the dissipation of the magnetic ﬁeld becomes strong and higher resolution

2 + j 2 ), whichis required. The third is the in-plane electron current density (jexz = jexez
is required because strong currents can excite microinstabilities so that high resolution
is needed (e.g., Fujimoto and Machida, 2003). In each time step, if λDe , Vey , and jexz
calculated at the center of a cell satisfy the condition, ∆L ≥ 2.0λDe or Vey > 2.0VA or
jexz > 0.5enps VA , the cell is subdivided and four child cells are produced, otherwise, the√child cells are removed if any. Here, VA = B0 / µ0 nps mi is the Alfvén velocity, and µ0
is the magnetic permeability in vacuum. The process begins from the base level cells
and proceed up to the level LD − 1 cells. However, for the level LD − 1 cells we do not
impose the criterion on λDe in order to avoid making the level LD cell region patchy.
The particle splitting is performed only for the background particles with level except
for LD in the present runs. The particle coalescence algorithm is also implemented forthe background particles once in a few hundred time steps, but we never choose particles
in the level LD cells to avoid the numerical errors in physically important regions.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

32Time evolutions of the reconnected magnetic ﬂux ψ are shown in Figure 2.6 for thethree runs: the ’Without AMR’ (black solid line), ’With AMR’ (black dashed line), and’With AMR & Particle splitting’ (black dotted line) runs. We found that, although theruns with the AMR require slightly more time before the onsets of fast reconnection
occur, each run has the almost identical gradient so the same reconnection rate after
the onset. Furthermore, the saturation levels for the runs with the AMR are clearlyin agreement with that for the ’Without AMR’ run. Thus we can say that the fastreconnection processes described in the three runs are physically identical, so that theadaptation of the AMR and particle splitting algorithms to the conventional PIC codes is
successful. We think that the discrepancy of the onset time is caused by the diﬀerence inthe noise level around the central current sheet at the initial stage in each run. Actually
the initial cell size around the central current sheet in the runs with the AMR is twiceof that in the ’Without AMR’ run, because we do not impose the criterion on λDe forthe level LD − 1 cells. Thus the number of particles per cell is increased four times in
the former, which indicates that the initial noise level around the central current sheet islower in the runs with the AMR. In order to make sure whether this reasoning is correct
or not, we also examined a simulation which imposes the reﬁnement criterion on λDe forthe level LD − 1 cells. In this case, the criterion on λDe is given as ∆L ≥ 1.4λDe , insteadof ∆L ≥ 2.0λDe , only for the level L = LD − 1 cells to avoid the patchy reﬁned regionswith level LD . The results for two cases are shown in Figure 2.6. One describes therun using only the AMR technique (red dashed line) and the other shows that using the
AMR and particle splitting techniques (red dotted line). It is found that these two runs
are in good agreement with the run without the AMR, which indicates that the initialnoise level around the central plasma sheet is very crucial to determine the time when
the onsets of magnetic reconnection occur. Indeed, it is very diﬃcult to discuss the timescale in which magnetic reconnection grows from noise to a signiﬁcant size, because it isstrongly dependent not only on the initial noise level but also on many other conditions
speciﬁc to the system (e.g., Shay et al., 2004). Thus we do not treat such a problem
in this paper. We only compare the description of fast reconnection, which allows us tostart simulations with coarser cells around the central plasma sheet.
Eﬃciency of our code is measured by the elapsed time to complete each run, which
is presented in Figure 2.7. Our code including both the AMR and particle splitting
algorithms is shown clearly to be very eﬃcient compared with the conventional PIC
code without the AMR. This eﬃciency is achieved by reducing the numbers of cells andparticles greatly (see Table I).
Figure 2.8 shows the evolution of the out-of-plane current density (Jy ) for the run in-CHAPTER 2. FULL PARTICLE CODE WITH AMR

Reconnected magnetic flux

7

3

Without AMRWith AMR
With AMR &Particle splittingWith AMR*
With AMR* &Particle splitting

2* The refinement criterion on λDe
is used at the LD-1 level.

6

ψ/B0λi0

33

5
4

1
0
0

5

10

tω ci

15

20

25

Figure 2.6: Time evolutions of the reconnected magnetic ﬂux deﬁned as ψ ≡lx|Bz (x, z = 0)| dx, where lx is the system size in the x direction. The black solid0line shows the result of the non-AMR run, where the whole simulation area is covered byﬁne cells having enough background particles, while the black dashed line denotes therun using only the AMR technique so that the number of particles is the same as in the
non-AMR run. The black dotted line is the result including both the AMR and particle splitting algorithms, so that the number of background particles is greatly reduced
compared with that in the former two cases. In addition to them, two runs in which thereﬁnement criterion on λDe is used also at the LD − 1 level are plotted in red lines.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

34

103

Elapsed time (min)

5
4
3
2
1
0

Without AMRWith AMR

With AMR &Particle splitting

Figure 2.7: Elapsed time to complete simulations for the cases of ’Without AMR’, ’With
AMR’, and ’With AMR & Particle splitting’ runs, which correspond to the solid, dashed,
and dotted lines in Figure 2.6, respectively. Each run has the same time steps. ’With
AMR & Particle splitting’ run is faster than ’Without AMR’ and ’With AMR’ runs bya factor of 5.0 and 4.1, respectively.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

35cluding the AMR and particle splitting algorithms. A strong thin current sheet is formedin the vicinity of the X-line (tωci = 14.72) and then fast reconnection occurs (see Figure2.6). After the initially loaded plasma in the current sheet is carried away from thecentral region, the background lobe plasma enters the diﬀusion region and quasi-steadyelectron diﬀusion region is produced (tωci = 17.92, 20.48). In the electron diﬀusion region, electrons are unmagnetized from the ambient magnetic ﬁeld, and conduct the meandering motions with acceleration due to the out-of-plane electric ﬁeld, which generate
a double-peaked current sheet therein as clearly shown in the ﬁgure. Similar propertiesare described also in other full particle simulations on magnetic reconnection (e.g., Shayet al., 2001; Pritchett, 2001; Hoshino et al., 2001b). The hierarchical cell distributionsin association with the current sheet evolution are shown in the right-hand column of
Figure 2.8. We ﬁnd that the ﬁnest cells are only distributed around the diﬀusion region,
separatrices extending from the X-line, and plasmoid, where the electron kinetic eﬀects
are expected to be important.
Figure 2.9 displays the proﬁles of (a) the number density of electrons (ne ), (b) the outof-plane current density (Jy ), (c) the inﬂow electron velocity (Vez ) and E×B drift velocity(dashed line), and (d) the out-of-plane electric ﬁeld (Ey ), along the x/λi0 = 7.9 axis at
the time tωci = 17.92 for the run including the AMR and particle splitting algorithms.
The electron diﬀusion region can be deﬁned as an area where the electron bulk velocity
is inconsistent with the local E × B drift velocity. Thus the area from z/λi0 = −0.34 to
0.34 is the electron diﬀusion region in this case. The double-peaked proﬁles of ne and
Jy inside the electron diﬀusion region are also shown in Figure 2.9a and b. The peakseparation is estimated at 0.3λi0 , which is approximately same as the electron inertial
length deﬁned by the electron number density averaged within the electron diﬀusion
region, ne /nps  0.071. Shay et al. (2001) has constructed a model around the electron
diﬀusion region, in which vez at the edge of the electron diﬀusion region and Ey within
the region are given by vez ∼ Bd /µ0 ene l and Ey ∼ −Bd2 /µ0ene l respectively, where Bd
is the x component of the magnetic ﬁeld at the edge of the electron diﬀusion region
and l is the length of the region along the x axis. When we estimate l, we assume thatelectrons continue to be accelerated along x in the electron diﬀusion region, so that theelectron bulk velocity Vex should be maximized at the edge of the region. In the current
case, Bd /B0  0.29 and l/λi0  2.0 are estimated from our simulation results, so thatvez /VA  2.1 and Ey /VA B0  −0.61 are derived, which are approximately consistent
with our results described in Figure 2.9c and d. Thus the electron dynamics described in
the run including the AMR and particle splitting algorithms is in good agreement withthe conventional model.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

3.8

3.8

Jy

t= 0.00

2.9

t= 0.00

2.9

1.9

1.9

1.0
Z

36

1.0
Z

0.0

0.0

-1.0

-1.0

-1.9

-1.9

-2.9

-2.9

-3.8
0.0

3.1

6.1

9.2

12.3

-3.8
0.0

15.4

3.1

6.1

X
3.8

3.8

Jy

t= 14.72

2.9

12.3

15.4

9.2

12.3

15.4

9.2

12.3

15.4

0.0

-1.0

-1.0

-1.9

-1.9

-2.9

-2.9

3.1

6.1

9.2

12.3

-3.8
0.0

15.4

3.1

6.1

X
3.8

X
3.8

Jy

t= 17.92

2.9

2.9

1.9

t= 17.92

1.9

1.0

1.0
Z

0.0

0.0

-1.0

-1.0

-1.9

-1.9

-2.9

-2.9

-3.8
0.0

3.1

6.1

9.2

12.3

-3.8
0.0

15.4

3.1

6.1

X
3.8

X
3.8

Jy

t= 20.48

2.9

2.9

1.9

t= 20.48

1.9

1.0

1.0
Z

0.0

0.0

-1.0

-1.0

-1.9

-1.9
-2.9

-2.9
-3.8
0.0

3.1

6.1

9.2

12.3

15.4

X

-2.00

9.2

t= 14.72

1.0
Z

0.0

-3.8
0.0

Z

15.4

1.9

1.0

Z

12.3

2.9

1.9

Z

9.2
X

-1.78

-1.56

-1.34

-1.12

-0.90

-3.8
0.0

3.1

6.1
X

-0.68

-0.46

-0.24

-0.02

0.20

Figure 2.8: Example of the evolution of the hierarchical cell distribution associated withthe current sheet evolution for the run with the AMR and particle splitting processes.
Magnetic ﬁeld and current density (color-coded) evolutions are shown in the left-handcolumn and the hierarchical cells at the times corresponding to them are described in
the right-hand column. Finest cells are only distributed around the diﬀusion region,
the separatrices extending from the X-line, and the plasmoid, where the electron kinetic
eﬀects are expected to be important.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

37

Ey / VAB0

Vez / VA

Jy / enpsVA ne / nps

t = 17.92

0.2
(a)
0.15
0.1
0.05
00-3
(b)
-0.2
-0.4
-0.6
-0.8
-1
-1.2-3
(c)
2
0
-2
0-3
(d)
-0.2
-0.4
-0.6
-0.8
-3

L=7
L=8
L=9
L = 10

-2

-1

0

1

2

3

-2

-1

0

1

2

3

-2

-1

0

1

2

3

-2

-1

0

1

2

3

z/ i0

Figure 2.9: Proﬁles of (a) the number density of electrons, (b) the out-of-plane current
density, (c) the inﬂow electron velocity and E × B drift velocity (dashed line), and (d)the out-of-plane electric ﬁeld, along the x/λi0 = 7.9 axis at the time tωci = 17.92 for thesame run as in Figure 2.8. The green, blue, purple, and red lines represent the proﬁles
on the level L = 7, 8, 9, and 10 cells, respectively.

CHAPTER 2. FULL PARTICLE CODE WITH AMR

2.4

38

Summary and Discussion

We have described a new electromagnetic full particle code with the AMR and particle splitting algorithms. The AMR technique subdivides and removes cells dynamicallyin accordance with the reﬁnement criteria and enhances the spatial resolution. On theother hand, the particle splitting algorithm divides particles that locate in ﬁner cells,
conserving the ratio of mass to charge, the total momentum, energy, and distributionfunction of particles, and the moments on grids. The particle splitting is indispensable
for adapting the AMR to the electromagnetic particle codes, because the numerical noise
increases on the subdivided cells due to the decrease in the number of particles per cell,
especially around the X-line formed in association with magnetic reconnection (Fujimotoand Machida, 2005a).
We conduct several test simulations and compare three runs without the AMR, with
the AMR, and with the AMR and particle splitting. Test simulations on the propagationof the Langmuir waves indicate that the AMR and particle splitting algorithms are
successfully applied to the conventional PIC codes, and we have found that the wave
reﬂections at the boundaries of the reﬁned regions are not essential. The nonlinear
evolution of the Harris-type current sheet and the electron dynamics around the X-line
using our AMR code are basically the same as those in other reconnection simulations.
However, the numbers of cells and particles are greatly reduced in the run including
the AMR and particle splitting algorithms, so that the time to complete the simulationis considerably shortened. Thus we conclude that we have realized eﬀectively highresolution simulations on the evolution of the current sheet by the use of the AMR
technique and particle splitting algorithm. Our code enables us to conduct large-scaleparticle simulations on magnetic reconnection in a wide range of the order 104 – 105 .
Our next challenge is a massively parallel computation for the AMR code using
several computational nodes. Our code at the current version is parallelized only within
a node of shared memory processors using OpenMP. It is worth mentioning that theaccumulation of the charge and current densities onto each cell corner is parallelized
by preparing temporal lists of charge and current densities for each processor. In otherwords, each processor has private lists of charge and current densities, and assigns particle
to the lists independently of the other processors. After a sweep of the particle loop, thecharge and current densities assigned to each processor are summed up respectively andsubstituted into the shared lists. This procedure is eﬀectively parallelized at the expenseof the computational memory. Diﬃculties in parallelizing our code arise mainly fromthe Poisson solver to calculate δφ by equation (2.15). The relaxation methods as usedin the code are not adequate to a parallel computing. This is conceptually because aCHAPTER 2. FULL PARTICLE CODE WITH AMR

39global distribution of the source ﬁeld contributes to each local potential, and technicallybecause the iteration loop repeated until the residual error becomes small enough must
be carried out sequentially. In order to avoid solving the Poisson equation, we can use
the charge conservation method for example, in which the current density is determined
directly from knowledge of charge motion so as to satisfy the charge continuity equation(e.g., Villasenor and Buneman, 1992). The improved version of our code will be discussed
in a separate article.

CHAPTER 3
Time Evolutions of the Electron
Diﬀusion Region and theReconnection Rate in Large andFully Kinetic System
3.1

Introduction

Magnetic reconnection is widely believed to play an important role in the magnetospheric substorm and solar ﬂares as a fast conversion process of the magnetic energy
to plasma kinetic and thermal energies. In reconnection processes, the ideal magnetohydrodynamic (MHD) condition breaks down in the diﬀusion region arising around amagnetic X-line, where dissipation processes of the magnetic ﬁeld are signiﬁcant (Vasyliunas, 1975). In a suﬃciently collisional plasma, the resistive MHD theory (Sweet,
1958; Parker, 1963; Petschek, 1964) is valid for describing the diﬀusion region by parameterizing eﬀects of the classical Coulomb collisions (Spitzer, 1962). However, in thecollisionless plasma as seen in the Earth magnetosphere, the magnetic diﬀusion due tothe classical resistivity is too small to drive fast reconnection. In such a system, we needto replace the standard Ohm’s law in the resistive MHD by the generalized Ohm’s law
written as
E +V ×B = −me
1
1(V e · ∇)V e −∇ · Pe +J × B,
e
ne
ne

(3.1)

where E is the electric ﬁeld, B the magnetic ﬁeld, J the current density, n the plasma
number density, V e the electron ﬂow velocity, P e the electron pressure tensor, and me
the electron mass. The ﬁrst and second terms in the right-hand side arise from electron
inertial eﬀects, and the scale sizes are evaluated as the electron inertial length λe and theelectron gyroradius λge , respectively. The last term represents the decoupling between

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 41

ions and electrons, so-called Hall eﬀects, and has a scale size of the ion inertial length
λi . When the scale length of gradient is much larger than any of these characteristic
lengths, the frozen-in condition holds.
In two-ﬂuid system consisting of ions and electrons, the inﬂowing ions can be decoupled from the ambient magnetic ﬁeld when they approach the neutral sheet as close as
λi . Thus the region within λi is the diﬀusion region where the frozen-in condition breaks
down. Electrons, however, still remain magnetized at this distance, since their inertia
is much smaller than those of ions. Instead, their frozen-in constraint can break downat the distance comparable with the local gyroradius λge , within which electrons are expected to perform the meandering motions (e.g., Pei et al., 2001). This non-magnetizedregion for electrons is called the electron diﬀusion region. Therefore the diﬀusion region
in the two-ﬂuid system develops two-scale structure, embedding the electron diﬀusion
region within the ion-scale diﬀusion region.
Inside the diﬀusion region but outside the electron diﬀusion region, ions are decoupledfrom the ambient magnetic ﬁeld but electrons are still frozen-in to the ﬁeld and continue
to move with the E × B drift motion toward the current sheet. This relative motion
between ions and electrons produces currents in the vicinity of the magnetic X-line,
forming the Hall current system (Sonnerup, 1979; Terasawa, 1983). On the other hand,
the kinetic behavior of electrons becomes important in the electron diﬀusion region.
It has been revealed by using two-dimensional kinetic simulations that the dissipation
process in the electron diﬀusion region is supported by the electron inertia term andthe gradient term of the electron pressure tensor in the generalized Ohm’s law (3.1)(Cai and Lee, 1997; Hesse and Winske, 1998; Kuznetsova et al., 1998; Hesse et al., 1999;
Swisdak et al., 2005). Though reconnecting of the magnetic ﬁeld actually proceeds inthis region, it has been suggested that the electron dynamics should have little or noeﬀect on the reconnection rate (Biskamp et al., 1997; Birn et al., 2001). Instead, ion
dynamics can control the reconnection processes so the electron-scale structure. They
conclude that a high reconnection rate results from the Hall term in the generalized
Ohm’s law (3.1), so that the inclusion of the Hall eﬀects is a suﬃcient condition toachieve fast reconnection. This indicates that three-dimensional eﬀects arising along thecurrent sheet would not be essential for obtaining the high reconnection rate (Huba andRudakov, 2004). Though recent simulation studies comparing MHD, Hall MHD, hybrid,
and full particle simulations conﬁrmed the importance of the Hall eﬀects, the system did
not reach steady state and the reconnection rate decreased after it reached a peak value(Shay et al., 2001; Hesse et al., 2001; Pritchett, 2001; Kuznetsova et al., 2001). They
attributed the depression to the periodic system in the direction of the initial magnetic

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 42

ﬁeld, in which compressional eﬀects in the magnetic islands could aﬀect the reconnection
processes. One of the solutions in order to suppress the boundary eﬀects is to performthe simulations in a large system so that the periodicity is negligible therein.
Such large-scale simulations including the Hall eﬀects were examined using a twodimensional hybrid code (Shay et al, 1999) and a two-dimensional two-ﬂuid code (Shayet al., 2004). Both results also indicate the importance of the Hall eﬀects for fast reconnection and show a quasi-steady reconnection with the reconnection rate of the orderof 0.1VA0B0 , where VA0 and B0 are, respectively, the Alfvén velocity and the magnetic
ﬁeld deﬁned in the asymptotic lobe region, and the reconnection rate is estimated by
the electric ﬁeld strength along the X-line. However, Karimabadi et al. (2004) recently
compared two diﬀerent hybrid simulations with and without the Hall term under anopen boundary condition, and demonstrated that ion kinetic behavior alone could giverise to fast reconnection even in the absence of the Hall eﬀects. Thus the Hall eﬀects
may not be necessarily required for fast reconnection.
The present work shows results from 2-1/2 dimensional full particle simulations ina large system, in which both ions and electrons are treated as superparticles. The
purpose of this study is to see a long evolution of the diﬀusion region in a fully kinetic
system. The main concern is how the ion dynamics aﬀects the electron-scale structureso the reconnection processes. It is still diﬃcult to conduct large-scale simulations using
conventional particle-in-cell (PIC) codes because of limited computer resources. Thus
we employed the adaptive mesh reﬁnement (AMR) technique and the particle splitting
algorithm to the conventional PIC code (Fujimoto and Machida, 2005b). Adoption ofthese techniques enables us to perform eﬀectively high-resolution simulations including
kinetic processes of magnetic reconnection.

3.2

Simulation Model

The present study employs a 2-1/2 dimensional electromagnetic PIC code with
the AMR technique and the particle splitting algorithm, which has been described by
Fujimoto and Machida (2005b). The AMR technique subdivides and removes computing
cells dynamically in accordance with a reﬁnement criterion and quite eﬀective to achieve
high-resolution simulations of phenomena that locally include micro-scale processes. In
our code, the spatial resolution is increased by introducing ﬁner cells hierarchically onto
the uniform base cells that cover the entire simulation area. If a base cell is reﬁned, fourchild cells that have half the size of the base cell are generated. These child cells can bealso reﬁned in turn and ﬁner cells are produced, and so on. One of the main problems in

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 43

developing the electromagnetic PIC codes using the AMR technique is the decrease in
the number of particles per cell in the reﬁned region. In order to solve this problem, we
subdivide particles residing in the subdivided cells and control the number of particles
per cell. We employ the particle splitting algorithm developed by Lapenta (2002), in
which the moments on each grid (the charge and current densities), the total charge,
mass, momentum, and energy of particles, and the distribution function of particles are
retained between before and after the particle splitting. Each cell is needed to haveinformation of the parent, child, and neighboring cells, and the particles residing in
the cell. The communication between cells or between a cell and the particles in thecell is supported by a set of pointers, constructing the fully threaded tree structure(Khokhlov, 1998). A reﬁnement level L in the hierarchical cell structure is deﬁned byusing the cell size of the level (∆L ) as L ≡ log2 (lz /∆L ), where lz is the vertical size ofthe two-dimensional simulation area. We use only cells with integer level. More detailed
description on our code is shown in Fujimoto and Machida (2005b).
The reﬁnement condition in the current study is deﬁned by three physical values. The√
ﬁrst is the local electron Debye length, λDe = vth,e / 2ωpe , which is required to avoid

a numerical heating of local plasma, where vth,e = 2Te /me is the electron thermalvelocity (Te and me are the temperature and mass of electrons, respectively), and ωpeis the electron plasma frequency. We use the initial value of vth,e for calculation of λDe ,because electrons are expected to be heated in the region where the electron dynamics
is important, that is, we use the minimum value of λDe . The second is the out-of-planeelectron ﬂow velocity (Vey ), because the inertial term (the ﬁrst term of the right-handside) in the generalized Ohm’s law (3.1) becomes dominant in the region where theelectron ﬂow is strong, so that the dissipation of the magnetic ﬁeld becomes strong andhigher resolution is required. The third is the in-plane electron current density (Jexz =
2 + J 2 ), which is required because intense currents can excite microinstabilities soJexezthat electron-scale waves are expected to arise. In each time step, if λDe , Vey , and Jexz
calculated at the center of a cell satisfy the condition, ∆L > 2.0λDe or Vey > 2.0VA or
Jexz > 0.5enps VA , the cell is subdivided and four child cells are produced, otherwise,
the child cells are removed if any. Here, VA is the Alfvén velocity deﬁned by the initial
lobe ﬁeld (B0 ) and plasma sheet density (nps ), and e is the electron charge. The particle
splitting is also performed in each time step, but this is done only for the backgroundparticles whose initial distribution is described later. We carry out not only particle
splitting but also particle coalescence in order to reduce the number of particles in
the unreﬁned region. The particle coalescence algorithm is implemented once in a few
hundred time steps for the background particles. In both the particle splitting andCHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 44

coalescence algorithms, we never choose particles in the most reﬁned (ﬁnest) cells to avoid
numerical errors in physically important regions. In the present study, the hierarchicalcell structure consists of four cell layers.
The initial magnetic ﬁeld conﬁguration is given as Bx (z) = −B0 tanh(z/λ), and thecorresponding density is n(z) = nps sech2 (z/λ), where λ is the half width of the initialplasma sheet and set as λ = 0.5λi0 (λi0 is the ion inertial length deﬁned by nps ). In
addition to the equilibrium, a small perturbation is superposed in the form
Bxp (x, z) = 2a0 /λ sech2 ((x − lx /2)/L ) sech2 (z/λ) tanh(z/λ)
Bzp (x, z) = −2a0 /L sech2 ((x − lx /2)/L ) sech2 (z/λ) tanh((x − lx /2)/L ),(3.2)

where, lx is the horizontal length of the simulation area, and a0 and L provide theamplitude and horizontal size of the perturbation, respectively. In the present study,
a0 = 0.15B0λi0 , L = 3.8λi0 , and lx × lz = 122.9λi0 × 30.7λi0 are chosen, unless otherwise
mentioned. Furthermore, the background plasma is loaded as nb (z) = nb0 tanh2 (z/λ)
and nb0 = 0.044nps in order to describe the lobe plasma in the magnetospheric tail.
Although there appears a weak pressure imbalance due to this background proﬁle, it isquickly justiﬁed without any signiﬁcant modiﬁcation of the current sheet structure. Oncethe simulation starts, the initial plasma sheet density and current proﬁles are quickly
modiﬁed to adjust the ﬁeld perturbation. As a result, a thinner current sheet is formednear the center of the simulation area (x = lx /2, z = 0), so that the tearing instability
selectively develops therein. The cell size is ∆LB = 0.12λi0 for the coarsest cells and∆LD = 0.015λi0 for the ﬁnest cells, and the time step is ∆tωci = 8.0 × 10−4 for allparticles and reﬁned regions in order to satisfy the Courant condition on the ﬁnest cells,
where ωci is the ion cyclotron frequency deﬁned by the lobe ﬁeld (B0 ). The initial plasma
condition is mi /me = 100, Ti,ps /Te,ps = 8.0, Ti,lobe /Te,lobe = 1.0, Te,lobe /Te,ps = 1.0, unlessotherwise mentioned, and c/VA = 16.7, where Ts,ps and Ts,lobe are the temperatures ofthe species s at the central plasma sheet and the magnetic lobe, respectively, and c isvelocity of light. We assume the periodic boundary in the x direction and the conducting
wall in the z direction.

3.3

Results

3.3.1 Time Evolution of the Reconnection Ratel /2

The evolution of the reconnected magnetic ﬂux, ψ = 0 x |Bz (x, z = 0)| dx, is plottedin Figure 3.1a. Here, it is assumed that the X-line is formed at the center of the simulation
area (x = lx /2, z = 0). It increases very slowly in the ﬁrst half of the simulation time, but

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 45

5(a)ψZB0 λi0

4
3
2
1
0

0

2.5

5

7.5

10

12.5

15

17.5

5

7.5

10

12.5

15

17.5

0.8

- EyZVAB0

(b)
0.6
0.4
0.2
0
0

2.5

tωci
Figure 3.1: Time evolutions of (a) the reconnected magnetic ﬂux, (b) the out-of-planeelectric ﬁeld at the center of the simulation area (x = lx /2, z = 0).

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 46
develops at drastically high rate after tωci ≈ 12.5, so that fast reconnection is switchedon. The reconnection rate should be given by the instantaneous growth rate of ψ, which
is derived from Faraday’s law as ∂ψ/∂t = −Ey (x = lx /2, z = 0), which expresses theelectric ﬁeld imposed along the X-line. Figure 3.1b shows the time evolution of theout-of-plane electric ﬁeld at the center of the simulation area. The reconnection rate
develops associated with the evolution of the magnetic reconnection until it reaches apeak at tωci ≈ 15, and a fast reconnection with |Ey | ∼ 0.7VA B0 = 0.15VA0B0 is achieved.
However, after then it starts to decrease and fall to nearly the half of the peak value atthe end of the simulation. This time history is similar to those obtained in other particle
simulations having smaller simulation areas, in which the periodicity of the system has
tended to aﬀect the reconnection rate and slow down reconnection processes (Shay etal., 2001; Hesse et al., 2001; Pritchett, 2001; Kuznetsova et al., 2001). In the present
run, however, the simulation area is large enough to describe longer time evolutions ofreconnection without the eﬀects of the periodic boundaries. Actually the plasma density
at the edge of the simulation area (x = 0, z = 0), which is initially set as n = 1.00nps ,is only slightly increased to n = 1.04nps at the end of the simulation. This indicatesthat the compressive eﬀects in the magnetic island as expected in the periodic system
are mostly negligible in the present large system.

3.3.2 Structure Change Around the Electron Diﬀusion Region
The density proﬁles along the z direction at x = lx /2 are shown in Figure 3.2b forions (red solid lines) and electrons (blue solid lines) at two diﬀerent times: tωci = 12.8
(left) and 16.7 (right). At tωci = 12.8, both ion and electron density proﬁles have peaks
at z = 0 that are reminiscent of the initial equilibrium proﬁles. On the other hand, attωci = 16.7, electrons have still a peak at z = 0 while ion peaks appear at the ﬂanks ofthe electron peak and a density hole is formed at z = 0. The double-peaked structurein the ion density proﬁle as seen at tωci = 16.7 indicates that most ions experience
meandering motions near the X-line, staying long around the turning points and passing
fast around z = 0. The red dashed lines in Figure 3.2a denote the density proﬁles of the
background cold ions, which are initially loaded only in the lobe region. Its proﬁle attωci = 12.8 is not yet changed from the initial proﬁle in essence. However, as magnetic
reconnection proceeds, the background ions approach toward the X-line and start the
meandering motions. It is found that the double-peaked structure at tωci = 16.7 is mainly
produced by the background ions meandering around the X-line. Because the hot ions
initially loaded in the current sheet have large velocity in the out-of-plane direction inorder to produce the current supporting the initial magnetic ﬁeld proﬁle, they can easily

JyZenpsVA EzZVAB0n Z nps

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 47
0.35
0.3 (a)
0.25
0.2
0.15
0.1
0.05
0
-3
-2
1 (b)
0.5
0
-0.5
-1
-1.50
-3
-2
-0.5 (c)t=12.8
ni
nbi
ne
-1

0

1

2

3

-1

0

1

2

3

Jiy+Jey
Jiy
Jey

-1
-1.5
-2
-3

-2

-1

0

z/λi0

1

2

0.12
0.1
0.08
0.06
0.04
0.02
0
-3
1
0.5
0
-0.5
-1
-1.50
-3
-0.5

t=16.7

-2

-1

0

1

2

3

-2

-1

0

1

2

3

-2

-1

0

1

2

3

-1
-1.5
3

-2
-3

z/λi0

Figure 3.2: Comparison of structures along the z direction through the center of thesimulation area (x = lx /2, z = 0) at two diﬀerent times: tωci = 12.8 (left) and 16.7
(right). (a) The number densities of ions (red solid lines) and electrons (blue solid lines).
The red dashed lines represent the contribution from the background cold ions. (b) The
electric ﬁeld along the z direction. (c) The out-of-plane current densities (black solidlines), which are calculated by summing up the ion and electron current densities (redand blue lines, respectively).

escape from the diﬀusion region due to the Lorentz force before their meandering motions
become dominant. Thus, as the density of the background cold ions dominates that ofthe plasma sheet ions in the vicinity of the X-line, the double-peaked structure becomesdominant. This structure change occurs around tωci = 15 in the present run.
The diﬀerence in density peak position between ions and electrons are mainly causedby their inertia diﬀerence, and give rise to a strong polarization electric ﬁeld Ez directing
toward the neutral sheet in the electron inﬂow region. Its proﬁles along the z direction are
shown in Figure 3.2b, and the two-dimensional snapshots are also given in Figure 3.3.
We ﬁnd that the meandering motions of the background ions enhance the amplitudeof Ez and broaden its imposed region not only in the z direction but also in the x
direction around the X-line. The time evolution of the electric ﬁeld directing toward theneutral sheet Ez = −z/|z| Ez , averaged over the electron inﬂow region (lx /2 − 0.2λi0 ≤x ≤ lx /2 + 0.2λi0 , 0.25λi0 ≤ |z| ≤ 0.5λi0 ), is plotted in Figure 3.4. Ez develops aftertωci ≈ 13 as magnetic reconnection proceeds. Especially, its growth rate is enhanced
after tωci ≈ 15, indicating that the background ions carrying out the meandering motions
become dominant around the X-line. Figure 3.2c shows the comparison of the out-ofplane current density proﬁles (black lines), which are calculated by summing up the ion
and electron current densities (denoted by red and blue lines, respectively). We ﬁnd that

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 48
the polarization electric ﬁeld Ez forces the inﬂow electrons to execute the E × B drift
toward the out-of-plane direction and enhances the out-of-plane electron current density
in the electron inﬂow region.
In order to see the role of the enhanced electron current in the electron inﬂow region,
we plot in Figure 3.5 the time evolution of each term in the right-hand side of Ampère’slaw, ∂Ey /∂t = −Jy /ε0 + c2 (∂Bx /∂z − ∂Bz /∂x), averaged over the electron inﬂow region(lx /2 − 0.2λi0 ≤ x ≤ lx /2 + 0.2λi0 , 0.25λi0 ≤ |z| ≤ 0.5λi0 ). The term contributed by
the electron current (blue solid line) decreases until tωci ≈ 13.5, which results from thethinning of the electron current sheet. Then this term, however, increases again due tothe enhancement of the out-of-plane electron drift caused by the polarization electricﬁeld. The important point in this ﬁgure is that the contribution from the x-gradientof the magnetic ﬁeld Bz also has a peak at tωci ≈ 13.5. The absolute value of ∂Bz /∂x
in the electron inﬂow region reﬂects a curvature of the magnetic ﬁeld lines, which is
expressed as κ  |Bx /B 2 ∂Bz /∂x| when |Bz /B| << 1, where B = Bx2 + Bz2 . Since Bx(not shown) is almost constant after tωci ≈ 14, κ is mostly proportional to |∂Bz /∂x| atits declining phase. Therefore it is found that the curvature of the magnetic ﬁeld lines
decreases as the out-of-plane electron current in the electron inﬂow region is enhanced.
This means that the magnetic ﬁeld lines in the upstream region tends to be parallel tothe x direction, leading to the extension of the electron meandering region along the x
direction. Note that the left-hand side of Ampère’s law (∂Ey /∂t, not shown) is mostly
zero within the present simulation time in comparison with each term in the right-handside.
Figure 3.6 shows the snapshots of the out-of-plane electron velocity Vey with themagnetic ﬁeld lines (white solid lines) around the X-line. The strong electron velocity
seen near the center of the simulation area indicates that electrons are accelerated by
Ey and perform the meandering motions. Thus this region is mostly interpreted as theelectron diﬀusion region. As the curvature of the magnetic ﬁeld lines decreases in theelectron inﬂow region, the electron meandering region so the electron diﬀusion region
tends to extend along the x direction. This is also evident in Figure 3.7, which describes
time evolutions of (a) the length l, (b) the width δ, and (c) the aspect ratio δ/l of theelectron diﬀusion region. The length l and the width δ of the electron diﬀusion region
are properly determined from the simulation results at a give time. The electron diffusion region should be the region where the dissipation process are mainly supportedby the electron dynamics, that is, the region where electrons are accelerated by the inductive electric ﬁeld Ey and experience the meandering/Speiser motions (Speiser, 1965).
Thus the downstream edge of the electron diﬀusion region is determined by the location

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 49

7.7

Ez

t= 12.8
3.8

Z

0.0

-3.8

-7.7
43.0

50.4

57.7

65.1

72.5

79.9

X
7.7

Ez

t= 14.8
3.8

Z

0.0

-3.8

-7.7
43.0

50.4

57.7

65.1

72.5

79.9

X
7.7

Ez

t= 16.7
3.8

Z

0.0

-3.8

-7.7
43.0

50.4

57.7

65.1

72.5

79.9

X

-1.00

-0.80

-0.60

-0.40

-0.20

0.00

0.20

0.40

0.60

0.80

1.00

Figure 3.3: Snapshots of the electric ﬁeld along the z direction at tωci = 12.8, 14.8, and16.7. Magnetic ﬁeld lines (white solid lines) are superposed on each panel.

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 50

1.75

Ez ZVAB0

1.5
1.25
1
0.75
0.5
0.25
0

0

5

2.5

7.5

10

12.5

15

17.5

tωci
Figure 3.4: Time evolution of the electric ﬁeld directing toward the neutral sheet Ez =−z/|z| Ez , averaged over the electron inﬂow region (lx /2 − 0.2λi0 ≤ x ≤ lx /2 + 0.2λi0 ,0.25λi0 ≤ |z| ≤ 0.5λi0 ).

Each term of Ampere’s law

250-Jy/ε0-Jiy/ε0-Jey/ε0-c2 ∂Bz/∂x-c2 ∂Bx/∂z

200

150

100

50

0
0

2.5

5

7.5

tω ci

10

12.5

15

17.5

Figure 3.5: Time history of each term in the right-hand side of Ampère’s law (∂Ey /∂t =−Jy /ε + c2 (∂Bx /∂z − ∂Bz /∂x)) averaged over the electron inﬂow region. The black solidline denotes the contribution of the current density, which is calculated by summing up
the ion and electron current densities (red and blue solid lines, respectively). The dashedand dashed-dotted lines represent the contribution from the rotation of the magnetic
ﬁeld.

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 51

where the electron outﬂow velocity along the x axis reaches a peak value. Outside theedge, electrons are decelerated and the electron ﬂow energy is partly converted to its
thermal energy. On the other hand, the upstream edge of the electron diﬀusion region
is determined by the turning point of the electron meandering motions. The amplitude of the electron meandering orbit is given by δm = λg (δm ), where λg (z) = vth /ωcis the local Larmor radius of electrons, vth = 2Te /me is the electron thermal velocity, and ωc = −eBx /me is the electron cyclotron frequency. The width δ is deﬁned by+−±+ (−δm), where δmdenotes the location of the upstream edge in the ±z region.
δ = δmThe length l of the electron diﬀusion region increases as time goes on after tωci ≈ 13.5
(Figure 3.7a) when the curvature of the magnetic ﬁeld lines reaches the peak. Since thewidth δ is mostly constant in time (Figure 3.7b), the aspect ratio δ/l decreases aftertωci ≈ 13.5, especially after tωci ≈ 15, so that a long and narrow diﬀusion region isformed (Figure 3.7c).
In Figure 3.8, a schematic ﬁeld structure around the electron diﬀusion region is
described with the polarization electric ﬁeld Ez and the out-of-plane current density
Jy . The electric ﬁeld directing toward the neutral sheet is signiﬁcant in the shadowed
regions. The out-of-plane current density Jy in the shadowed region is enhanced in themanner that the polarization electric ﬁeld Ez forces the inﬂow electrons to perform theE × B drift to the out-of-plane direction and enhances the out-of-plane electron current
density. Since the out-of-plane current density is almost uniform along the upstreamedge of the electron diﬀusion region, the magnetic ﬁeld Bz does not appear near the Xline. Thus the magnetic ﬁeld lines in the electron inﬂow region should be parallel to the
upstream edge, so that the curvature of the magnetic ﬁeld lines is reduced in association
with the enhancement of Jy in the electron inﬂow region. On the other hand, near thedownstream edge of the shadowed region, Jy has a gradient along the x direction, which
gives rise to the magnetic ﬁeld ﬂuctuation δBz (dashed arrows in Figure 3.8) and also
reduces the curvature of the magnetic ﬁeld lines. As a result, the magnetic ﬁeld lines inthe electron inﬂow region tend to be parallel to the x axis, leading to the extension of
the electron meandering region so the electron diﬀusion region along the x direction.

3.3.3 Impact on the Reconnection Rate
We consider here how the structure change of the electron diﬀusion region aﬀects
the reconnection rate, that is, the electric ﬁeld along the X-line. Figure 3.9 shows
the ﬂow velocity of electrons along the z direction at x = lx /2. This proﬁle is usefulto understand the electron behavior around the X-line, when the intense polarizationelectric ﬁeld is imposed. The dashed line expresses the E × B drift velocity in the z

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 52

7.7

Vey

t= 12.8
3.8

Z

0.0

-3.8

-7.7
43.0

50.4

57.7

65.1

72.5

79.9

X
7.7

Vey

t= 14.8
3.8

Z

0.0

-3.8

-7.7
43.0

50.4

57.7

65.1

72.5

79.9

X
7.7

Vey

t= 16.7
3.8

Z

0.0

-3.8

-7.7
43.0

50.4

57.7

65.1

72.5

79.9

X

-2.00

-0.80

0.40

1.60

2.80

4.00

5.20

6.40

7.60

8.80

10.00

Figure 3.6: Snapshots of the out-of-plane electron ﬂow velocity at tωci = 12.8, 14.8, and16.7. Magnetic ﬁeld lines (white solid lines) are superposed on each panel.

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 53

8

(a)
lZλi0

6
4
2
0

0

5

7.5

10

12.5

15

17.5

5

7.5

10

12.5

15

17.5

5

7.5

10

12.5

15

17.5

(b)

0.5

δZλi0

2.5

0.4
0.3
0.2
0.1
0

0

2.5

0.14

(c)

0.12

δZl

0.1
0.08
0.06
0.04
0.02
0

0

2.5

tωci
Figure 3.7: Time evolutions of (a) the length, (b) the width, and (c) the aspect ratio of
the electron diﬀusion region, respectively.

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 54

Electron diffusion region

Ez

Jy

δBz

Figure 3.8: Schematic ﬁeld structure including the polarization electric ﬁeld Ez andout-of-plane current around the electron diﬀusion region. The electric ﬁeld directing
toward the neutral sheet is signiﬁcant in the shaded regions. The out-of-plane currents
in the shaded regions produce the magnetic ﬁeld δBz (dashed arrows) at the edges,
which reduces the curvature of the ﬁeld lines in the electron inﬂow region and extends
the electron diﬀusion region along the x axis.

direction, that is, −Ey /Bx . The E × B drift velocity diverges from the electron ﬂow
velocity at z ≈ ±0.7λi0 , thereby the frozen-in condition for electrons breaks down in theregion between z ≈ ±0.7λi0 . However, the upstream edges of the electron diﬀusion region
where meandering electrons turn back to the neutral sheet is located at z ≈ ±0.15λi0 in
our deﬁnition. Actually, the strong polarization electric ﬁeld imposed inside the ion-scale
diﬀusion region makes the inertial term in the generalized Ohm’s law (3.1) considerablylarge and suppresses the electron inﬂow velocity. In this case, the generalized Ohm’s law(3.1) can be expressed in the form
Ey + Vez Bx  −

∂Vey
me
Vez
,
e
∂z

(3.3)

where Vey  Ez /Bx is a good approximation. Thus we can obtain the electron inﬂow
velocityVez  −Ey
1,
−1 ∂V /∂z B
1 + ωce
ey
x

(3.4)

where ωce is the local electron cyclotron frequency. We call this velocity the “modiﬁed
E × B drift velocity”. This proﬁle is plotted by the dotted line in Figure 3.9. The
modiﬁed E × B drift velocity diverges from the electron inﬂow velocity at z ≈ ±0.15λi0 ,which are consistent with the upstream edges of the electron diﬀusion region.
Assuming that the inﬂow electrons along the z direction are mostly accelerated toward
the x direction by the magnetic tension in the electron diﬀusion region, we can obtain

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 55
4

t=16.7Vez-Ey/Bx-1/(1+1/ωce(Vey/(z) Ey/Bx

Vez /VA

2

0

-2

-4

-3

-2

-1

0

1

2

3

z/λi0
Figure 3.9: Flow velocity of electrons (solid line) along the z direction at x = lx /2 andtωci = 16.7. Dashed line denotes the E × B drift velocity deﬁned by −Ey /Bx , anddotted line is a velocity described by −1/(1 + 1/ωce∂Vey /∂z) Ey /Bx , which expresses the“modiﬁed E × B drift velocity” due to a large Ez .

the approximation−me ne1 Ve1 Ve2 

B1 B2
,
µ0

(3.5)

where nej , Vej , and Bj are, respectively, the electron density, velocity, and magnetic
ﬁeld in the location denoted by j (j = 1 and 2 represent the upstream and downstream+
edges of the electron diﬀusion region, respectively), that is, Ve1 = Vez (x = lx /2, z = δm),+), B2 = Bz (x = l+ , z = 0), and l+ is
Ve2 = Vex (x = l+ , z = 0), B1 = Bx (x = lx /2, z = δmthe location of the right-hand-side edge of the electron diﬀusion region. The out-of-planeelectric ﬁeld Ey is assumed to be spatially uniform around the electron diﬀusion region,
so that


∂Vey
me


Ey  −Ve1 B1  −Ve1 B1 −Ve1
 Ve2 B2 ,
(3.6)
e
∂z 1
where the prime denotes the location at which the electron inﬂow velocity diverges fromthe E × B drift velocity (−Ey /Bx ), that is, Ve1 = Vez (x = lx /2, z = 0.7λi0 ), and
B1 = Bx (x = lx /2, z = 0.7λi0 ) at tωci = 16.7, for example. We further assume Ve1  Ve1 ,which implies for the incompressible ﬂuid that all the electrons entering the region where
the frozen-in condition breaks down reach the electron diﬀusion region deﬁned in theCHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 56
present study. From (3.5) and (3.6) using l|Ve1 |  δ|Ve2 | for continuity, one can derive
Ve2  −

B1
VAe ,
B1

1
δ
Ey  − 1 +
l
ωce,1

(3.7)


∂Vey
∂z

 3/2
1

B1 VAe ,(3.8)√
where VAe = B1 / µ0 ne1 me is called the electron Alfvén velocity. Since B1 /B1  1,
it is found that the electron outﬂow velocity can exceed the electron Alfvén velocity
deﬁned at the upstream edge of the electron diﬀusion region. The theoretical predictionof the out-of-plane electric ﬁeld using (3.8) is plotted as a function of time in Figure 3.10
(dashed line), superposed on the direct measurement from the simulation results (solidline). This theoretical prediction is calculated from the upstream parameters (ne1 , B1 ,and (∂Vey /∂z)1 ) and the size of the electron diﬀusion region (δ and l) determined by
the simulation results at each given time. It increases as time goes on until it reaches
a peak at tωci ≈ 15 and then starts to drop associated with the reduction of the aspectratio (δ/l), that is, with the extension of the electron diﬀusion region. This proﬁle isconsistent with the simulation results (solid line). The damping rate of Ey after tωci ≈ 15
is slightly diﬀerent between the simulation results and the theoretical prediction. We
think that this is caused by the underestimation of δ. In our model, we assume that all
the inﬂow electrons are accelerated toward the x direction only in the electron diﬀusion
region deﬁned in the present study. However the polarization electric ﬁeld Ez directing to
the neutral sheet forces the inﬂow electrons toward the out-of-plane direction, even when
they reside outside the electron diﬀusion region. Thus the acceleration region should be
wider along the z direction than the electron meandering region, though it is diﬃcult tocorrectly determine the eﬀective value of δ.

3.3.4 Comparison With the Case of mi /me = 1
In order to emphasize the role of the polarization electric ﬁeld caused by the ion
meandering motions around the X-line, we compare the results discussed in the previous
sections with those obtained from the simulations under mi /me = 1, in which the decoupling between ions and electrons does not arise, so that the polarization electric ﬁeld iseliminated. At the same time, we check the eﬀects of the conducting wall located at theboundary of the simulation area in the z direction by performing same simulations inlarger system in the z direction. In the runs with mi /me = 1, the temperature ratio in
the initial plasma sheet (Ti,ps /Te,ps ) is also set as unity in order to avoid the electrostaticﬁeld arising from the diﬀerence in Larmor radius between ions and electrons. The system

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 57
0.8

Simulation
Theory- EyZVAB0

0.6
0.4
0.2
0
0

2.5

5

7.5

10

tωci

12.5

15

17.5

Figure 3.10: Comparison of the time evolutions of the out-of-plane electric ﬁeld. The
solid line denotes the direct measurement of the simulation results at the center ofthe simulation area. The dashed line is the theoretical prediction calculated from the
upstream parameters (ne1 , B1 , and (∂Vey /∂z)1 ) and the size of the electron diﬀusion
region (δ and l) determined by the simulation results.

parameters in the case of mi /me = 1 are ∆LB = 0.96λi0 , ∆LD = 0.12λi0 , ∆tωci = 0.01,
and c/VA = 6.0. The other conditions are the same as the case of mi /me = 100.
We examine several system sizes: lx × lz = 122.9λi0 × 30.7λi0 , and 122.9λi0 × 61.4λi0
for the mi /me = 100 runs, and lx × lz = 122.9λi0 × 30.7λi0 , 122.9λi0 × 61.4λi0 , and122.9λi0 × 122.9λi0 for the mi /me = 1 runs.
The time evolution of the out-of-plane electric ﬁeld at the center of the simulation
area (the reconnection rate) is shown in Figure 3.11 for each run. In the case of mi /me =1, we ﬁnd that a steady-state reconnection is achieved in larger system, although thereconnection rate in the smaller system (lx × lz = 122.9λi0 × 30.7λi0 ) tends to drop
due to the eﬀects of the conducting wall boundary. On the other hand, in the case of
mi /me = 100, the time evolution of the reconnection rate is independent of the systemsize, rising until tωci ≈ 15 and then falling down associated with the extension of theelectron diﬀusion region. These results demonstrate the importance of the polarizationelectric ﬁeld Ez in suppressing the reconnection processes. Another important pointis that the asymptotic reconnection rate in the mi /me = 1 runs (|Ey | ∼ 0.45VA B0 =0.09VA0B0 ) is large enough to explain fast dissipative events, even though the Hall eﬀects
are not included. Therefore it is suggested that the Hall eﬀects should not be necessarily
required for fast reconnection. This is well consistent with the results of Karimabadi etal. (2004)

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 58- EyZVAB0

0.6

0.4

0.2
lx*lz=123*30.7, mi/me=100
lx*lz=123*61.4, mi/me=100
lx*lz=123*30.7, mi/me=1
lx*lz=123*61.4, mi/me=1
lx*lz=123*123 , mi/me=1

0
0

5

10

15

tωci

20

25

30

Figure 3.11: Comparison of the time evolutions of the out-of-plane electric ﬁeld at thecenter of the simulation area (the reconnection rate). Red lines denote the runs formi /me = 100, and black lines the runs for mi /me = 1. Two system sizes are examined
under each mass ratio: lx × lz = 122.9λi0 × 30.7λi0 (solid lines), and 122.9λi0 × 61.4λi0(dashed lines). Larger system size with lx × lz = 122.9λi0 × 122.9λi0 is also employed
under mi /me = 1 (shown by black dotted line).

3.4

Summary and Discussion

In this paper, we have investigated long time evolutions of the electron diﬀusion
region embedded in the ion-scale diﬀusion region and the resulting reconnection rate associated with magnetic reconnection in a fully kinetic and large system. We performed
2-1/2 dimensional full particle simulations with the adaptive mesh reﬁnement (AMR)
technique and the particle splitting algorithm (Fujimoto and Machida, 2005b). The
adoption of these techniques to the conventional particle-in-cell (PIC) code enables us
to perform eﬀectively high-resolution simulations including kinetic processes of magnetic
reconnection. Previous reconnection studies using full particle simulations in the periodic
system (Shay et al., 2001; Hesse et al., 2001; Pritchett, 2001; Kuznetsova et al., 2001)
achieved fast reconnection, but the system did not reach steady state and the reconnection rate decreased after it reached a peak value. This is because the compressionaleﬀects in the magnetic islands aﬀect the reconnection processes in small periodic systems.
The full particle code with the AMR and particle splitting employed in the present study
allows us to take a large system, in which the eﬀects of the periodic boundary hardly
aﬀect the reconnection processes within the present simulation time.
We ﬁnd that the reconnection rate increases associated with magnetic reconnection
and reaches a peak value large enough for fast reconnection, but then it decreases astime goes on, even though the periodicity of the system is negligible. The key process

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 59

responsible for slowing magnetic reconnection is the extension of the electron diﬀusion
region associated with the evolution of the strong polarization electric ﬁeld Ez imposedin the electron inﬂow region. The polarization electric ﬁeld directing toward the neutral
sheet is caused by the inertia diﬀerence between ions and electrons, and enhanced by themeandering motions of the background ions. Because the hot ions initially loaded in the
Harris-type current sheet have a large velocity in the out-of-plane direction, they caneasily escape from the diﬀusion region due to the Lorentz force before their meandering
motions become dominant. Thus their impact on the strong Ez is small. The polarizationelectric ﬁeld Ez forces the inﬂow electrons to move toward the out-of-plane direction by
the E × B drift and enhances the out-of-plane current density Jy in the electron inﬂow
region. The role of Jy is to reduce the curvature of the magnetic ﬁeld lines so that they
are mostly parallel to the x direction near the X-line, in such a way that Jy becomes
uniform along the upstream edge of the electron diﬀusion region. As a result, the electron
meandering region so the electron diﬀusion region extends along the x direction.
We performed a Sweet-Parker-like analysis around the electron diﬀusion region, taking into account the breakdown of the electron frozen-in condition in the electron inﬂow
region, and demonstrate that the extension of the electron diﬀusion region can suppress
the reconnection rate. In order to emphasize the role of the polarization electric ﬁeld
and the Hall eﬀects, we compare the simulation runs with mi /me = 1 and mi /me = 100.
It is found that (1) a steady-state reconnection is achieved in the mi /me = 1 cases where
the polarization electric ﬁeld does not appear, diﬀerent from the mi /me = 100 cases, (2)a large reconnection rate with |Ey | ∼ 0.09VA0B0 is obtained in the mi /me = 1 runs, eventhough they do not include the Hall eﬀects, so that the Hall eﬀects are not necessarily
required for fast reconnection.
The polarization electric ﬁeld Ez caused by the ion meandering motions has been
seen in previous studies within large systems using hybrid simulations (Shay et al, 1999;
Arzner and Scholer, 2001). However, their systems achieved a quasi-steady reconnection,
diﬀerent from the present results. This inconsistency is possibly attributed to the absenceof the electron meandering eﬀects in the previous simulations. The damping of thereconnection rate results from the extension of the electron diﬀusion region in whichelectrons carry out the meandering motions. The electron meandering eﬀects should
be incorporated into the term arising from the nongyrotropic electron pressure of thegeneralized Ohm’s law (the second term in the right-hand side of (3.1)) (Cai and Lee,
1997; Hesse and Winske, 1998; Kuznetsova et al., 1998; Hesse et al., 1999; Swisdak etal., 2005). However, both of the previous studies did not include this term.
A long time evolution of the diﬀusion region using full particle simulations has been

CHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 60
15
10

Ion
Electronvx /VA

5
0
-5
-10

-10

-5

0

vy /VA

5

10

15

Figure 3.12: Distributions in the vx -vy velocity space of ions (crosses) and electrons (opentriangles) at the center of the simulation area at tωci = 14.8.

investigated by Pei et al. (2001) under an open system. Their system also reached asteady-state reconnection and the structure was mostly unchanged during the simulation
time. Magnetic reconnection in their system is driven by an external electric ﬁeld,
supplying an external plasma into the system. The reconnection rate in such a driven
system is mainly controlled by the external driving electric ﬁeld (Sato and Hayashi, 1979;
Pei et al., 2001). Thus the size of the diﬀusion region should be determined in the mannerthat the output of the plasma and magnetic ﬁeld ﬂux from the system becomes consistent
with the input from the external region. Thus the reconnection processes in the driven
system are essentially diﬀerent from those in the spontaneous system as employed in thepresent study, in which there is no energy input from the external region.
The question that remains is what can support fast reconnection in actual spaceplasma after the reconnection rate reaches the peak. In fact, the peak value of the reconnection rate obtained in the present runs is ∼ 0.7VA B0 . If we take B0 ∼ 10 nT andnps ∼ 0.5 cm−3 , we can calculate |Ey | ∼ 2.2 mV/m, which is not unrealistic comparedwith observational results (Asano et al., 2004). However, in our model, the reconnection rate is expected to continue decreasing after the meandering motions due to thebackground cold ions have been dominant around the X-line. Thus it is supposed that
other mechanisms that are not included in the model might play an important role inenhancing the reconnection rate. In the two-dimensional system employed in the present
study, plasma instabilities excited along the y direction are suppressed in principle. Actually, since electrons in the electron diﬀusion region can be strongly accelerated toward
the y direction beyond the electron Alfvén velocity deﬁned in the upstream region, aCHAPTER 3. EVOLUTION OF RECONNECTION RATE IN LARGE SYSTEM 61

large velocity diﬀerence between ions and electrons is expected to arise. In Figure 3.12,
distribution in vx − vy space at the center of the simulation area and tωci = 14.8 areshown for ions (crosses) and electrons (open triangles). The velocity diﬀerence between
ions and electrons along the y direction is beyond the thermal spread of the electron
distribution function. In such a condition, the Buneman-type instability can be excitedalong the y direction. Indeed, such an instability and its nonlinear structure have been
demonstrated in the three-dimensional system with a strong guide ﬁeld (Drake et al.,2003), and can give rise to the anomalous resistivity, so that the reconnection rate (|Ey |)is enhanced. In order to investigate the role of the Buneman-type instability in the electron diﬀusion region, it is necessary to perform a large-scale kinetic simulation in thethree-dimensional system.

CHAPTER 4
A Generation Mechanism ofElectrostatic Waves and SubsequentElectron Heating in the Plasma
Sheet-Lobe Boundary Region
During Magnetic Reconnection
4.1

Introduction

Magnetic reconnection is widely believed to play an important role in the magnetospheric substorm and solar ﬂares as a fast conversion process of the magnetic energy
to kinetic and thermal energy of ions and electrons. In fact, electron heating and acceleration arising during a number of phenomena in space are recognized as signaturesof magnetic reconnection. For example, in the Earth magnetotail, bursts of high energy
electrons have been reported in association with topological change of the near-Earthﬁeld line, so that they have been attributed to magnetic reconnection (Terasawa andNishida, 1976; Baker and Stone, 1976). In solar ﬂares, x-ray observations indicate that alarge part of the total energy is released in energetic electrons accelerated via magnetic
reconnection (Lin and Hudson, 1976; Lin et al., 2003). Plasma heating and acceleration
processes in association with magnetic reconnection have been discussed over the last
decades. However, our understanding on the energization mechanisms of electrons still
remains poor.
It is well known that non-Maxwellian electrons selectively heated parallel to the ambient magnetic ﬁeld, having an velocity distribution of the ﬂat-topped form, have been
often observed in the plasma sheet-lobe boundary region located between the magnetotail lobe and the plasma sheet in association with magnetic reconnection (Saito et al.,CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 63

1995; Shinohara et al., 1998; Hoshino et al., 2001a). One of the candidates responsible forthe anomalous heating is the acceleration in the slow-mode shocks (Schwartz et al., 1987;
Saito et al., 1995), which can explain a fast energy conversion from the stored magnetic
energy to the plasma kinetic energy (Petschek, 1964). However, Saito et al. (1995) havedemonstrated that the plasma sheet-lobe boundaries do not necessarily satisfy the slow
mode jump condition even if the electron ﬂat-topped distributions are observed in thedownstream region. This indicates that the conventional slow-mode shocks might notplay a signiﬁcant role in the electron heating in the Earth magnetotail. Shinohara et al.
(1998) also observed the ﬂat-topped electrons associated with an electromagnetic turbulence near the X-type neutral line. They concluded that the electromagnetic turbulence
results from the lower hybrid drift instability (LHDI) and suggested that the electron
heating should be closely related with the LHDI. Hoshino et al. (2001a), on the other
hand, discussed the physical processes producing the ﬂat-topped electrons by comparing the satellite observations and the electromagnetic particle-in-cell (PIC) simulations.
They suggested that the Buneman instability (Buneman, 1958) excited due to the velocity diﬀerence between the stationary ions and the beam electrons accelerated near the
X-type neutral line and/or the electron bump-on-tail instability (Omura et al., 1996)
arising between the weak beam electrons and the stationary hot electrons should be responsible for the electron scattering in the plasma sheet-lobe boundary region. However,they do not specify the nature of the instabilities in detail possibly because the simulation
area with the periodic boundaries is too small to describe the nonlinear feature of such
instabilities. In fact, a small system with the periodic boundaries can aﬀect nonlinear
wave-particle interactions in the boundary region via high energy particle passing acrossthe boundaries.
Wave activities in the plasma sheet-lobe boundary region of the Earth magnetotail have been extensively studied by means of the Geotail satellite observations (e.g.,
Matsumoto et al., 1994; Kojima et al., 1994). Their main ﬁnding is that most of thebroadband noise in the wave spectrum obtained in the boundary region are not continuous noise but are composed of a series of solitary pulses called the electrostatic solitary
waves (ESW). The generation mechanism of the ESW observed in the magnetotail is
explained by the nonlinear mode of the electron two-stream instabilities (e.g., Omura etal., 1996) including bump-on-tail instability or the Buneman instability as reported by
Drake et al. (2003). Therefore the electron beam ﬂowing relatively to the backgroundelectrons or ions is inevitable in order to generate the ESW. Magnetic reconnection canbe one of the strong candidates for the process producing such an electron beam (Kojima
et al., 1994). Actually, the satellite observations on the dayside magnetopause boundary

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 64

reported by Matsumoto et al. (2003) suggest that magnetic reconnection should be related to the ESW activities. Recently, Cattell et al. (2005) have observed large-amplitude
solitary waves and associated electron holes in the plasma sheet-lobe boundary region of
the magnetotail when near-Earth reconnection has occured. Their conclusion is that thedetected ESW result from the Buneman instability, and not from the electron two-stream
instabilities. The observations by Matsumoto et al. (2003) and Cattell et al. (2005) indicate that the electrostatic wave activities may play an important role in the electronenergization in the plasma sheet-lobe boundary region during magnetic reconnection.
In this paper, we study electron heating process in detail expected in the boundary between the magnetotail lobe and the plasma sheet associated with magnetic reconnection
by performing large-scale PIC simulations using the adaptive mesh reﬁnement (AMR)
and particle splitting techniques (Fujimoto and Machida, 2005b). We show that intenseelectrostatic waves are excited in the plasma sheet-lobe boundary region by the electron
two-stream instability between the cold electrons loaded in the lobe and the strong beamelectrons with high temperature. The electrostatic waves propagating along the ambient
magnetic ﬁeld scatter and heat the electrons along the ﬁeld line, and form the ﬂat-topped
electron distribution function. We also demonstrate that the ESW are generated as aresult of the nonlinear evolution of the electron two-stream instability during magnetic
reconnection, even when a strong guide ﬁeld as imposed in Drake et al. (2003) is absent.

4.2

Simulation Model

The present study employs a 2-1/2 dimensional electromagnetic PIC code with the
AMR technique and the particle splitting algorithm, which has been described by Fujimoto and Machida (Fujimoto and Machida, 2005a,b). The AMR technique subdivides
and removes computing cells dynamically in accordance with a reﬁnement criterion andquite eﬀective to achieve high-resolution simulations of phenomena that locally include
micro-scale processes. In our code, the spatial resolution is increased by introducing ﬁner
cells hierarchically onto the uniform base cells that cover the entire simulation area. If
a base cell is reﬁned, four child cells that have half the size of the base cell are generated. These child cells can be also reﬁned in turn and ﬁner cells are produced, and soon. One of the main problem in developing the electromagnetic PIC codes using the
AMR technique is the decrease in the number of particles per cell in the reﬁned region(Fujimoto and Machida, 2005a). In order to solve this problem, we subdivide particles
residing in the subdivided cells and control the number of particles per cell. We employ
the particle splitting algorithm developed by Lapenta (2002), in which the moments on

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 65

each grid (the charge and current densities), the total charge, mass, momentum, andenergy of particles, and the distribution function of particles are retained between before
and after the particle splitting. Each cell is needed to have information of the parent,child, and neighboring cells, and the particles residing in the cell. The communication
between cells or between a cell and the particles in the cell is supported by a set ofpointers, constructing the fully threaded tree structure (Khokhlov, 1998). A reﬁnementlevel L in the hierarchical cell structure is deﬁned by using the cell size of the level (∆L )as L ≡ log2 (lz /∆L ), where lz is the vertical size of the two-dimensional simulation area.
We use only cells with integer level. More detailed description on our code is shown in
Fujimoto and Machida (2005b).
The reﬁnement condition in the current study is deﬁned by three physical values. The√
ﬁrst is the local electron Debye length, λDe = vth,e / 2ωpe , which is required to avoid

a numerical heating of local plasma, where vth,e = 2Te /me is the electron thermalvelocity (Te and me are the temperature and mass of electrons, respectively), and ωpeis the electron plasma frequency. We use the initial value of vth,e for calculation of λDe ,because electrons are expected to be heated in the region where the electron dynamics
is important, that is, we use the minimum value of λDe . The second is the out-ofplane electron ﬂow velocity (Vey ), because the inertial term in the generalized Ohm’s law
becomes dominant in the region where the electron ﬂow is strong, so that the dissipationof the magnetic ﬁeld becomes strong and higher resolution is required. The third is the
2 + j 2 ), which is required because intense
in-plane electron current density (jexz = jexezcurrents can excite microinstabilities so that electron-scale waves are expected to arise.
In each time step, if λDe , Vey , and jexz calculated at the center of a cell satisfy thecondition, ∆L > 2.0λDe or Vey > 2.0VA or jexz > 0.5enps VA , the cell is subdivided andfour child cells are produced, otherwise, the child cells are removed if any. Here, VA is
the Alfvén velocity deﬁned by the initial lobe ﬁeld (B0 ) and plasma sheet density (nps ),and e is the electron charge. The particle splitting is performed only for the backgroundparticles whose distribution is described later. The particle coalescence algorithm is also
implemented for the background particles once in a few hundred time steps. However,in both the particle splitting and coalescence algorithms, we never choose particles in
the most reﬁned (ﬁnest) cells to avoid numerical errors in physically important regions.
In the present study, the hierarchical cell structure consists of four cell layers.
The initial magnetic ﬁeld conﬁguration is given as Bx (z) = −B0 tanh(z/λ), and thecorresponding density is n(z) = nps sech2 (z/λ), where λ is the half width of the initialplasma sheet and set as λ = 0.5λi0 (λi0 is the ion inertial length deﬁned by nps ). In

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 66

addition to the equilibrium, a small perturbation is superposed in the form,
Bxp (x, z) = 2a0 /λ sech2 ((x − lx /2)/L ) sech2 (z/λ) tanh(z/λ)
Bzp (x, z) = −2a0 /L sech2 ((x − lx /2)/L ) sech2 (z/λ) tanh((x − lx /2)/L ),(4.1)

where, lx is the horizontal length of the simulation area, and a0 and L provide theamplitude and horizontal size of the perturbation, respectively. In the present study,
a0 = 0.15B0 λi0 , L = 3.8λi0 , and lx × lz = 122.9λi0 × 30.7λi0 are chosen. Furthermore,
the background plasma is loaded as nb (z) = nb0 tanh2 (z/λ) and nb0 = 0.044nps in orderto describe the lobe plasma in the magnetospheric tail. Although there appears a weakpressure imbalance due to this background proﬁle, it is quickly justiﬁed without any
signiﬁcant modiﬁcation of the current sheet structure. Once the simulation starts, theinitial plasma sheet density and current proﬁles are quickly modiﬁed to adjust the ﬁeld
perturbation. As a result, a thinner current sheet is formed near the center of thesimulation area, so that the tearing instability selectively develops therein. The cell sizeis ∆LB = 0.12λi0 for the coarsest cells and ∆LD = 0.015λi0 for the ﬁnest cells, and thetime step is ∆tωci = 8.0 × 10−4 for all particles and reﬁned regions in order to satisfy the
Courant condition on the ﬁnest cells, where ωci is the ion cyclotron frequency deﬁned
by the lobe ﬁeld (B0 ). The initial plasma condition is mi /me = 100, Ti,ps /Te,ps = 8.0,
Ti,lobe /Te,lobe = 1.0, Te,lobe /Te,ps = 1.0, and c/VA = 16.7, where Ts,ps and Ts,lobe are
the temperatures of the species s at the central plasma sheet and the magnetic lobe,
respectively, and c is velocity of light. We assume the periodic boundary in the x
direction and the conducting wall in the z direction.

4.3

Results

4.3.1 Parallel Electron Heating in the Plasma Sheet-Lobe Boundary Region
In order to investigate the electron heating process in the plasma sheet-lobe boundary
region, we have performed large-scale kinetic simulations so that the non-realistic eﬀects
arising from the periodicity of the system, such as the compression of the magnetic ﬁeld
and plasma in the magnetic island and high energy particles passing across the periodic
boundary are negligible in the extensive region around the X-line. A snapshot of theout-of-plane current density with magnetic ﬁeld lines at tωci = 15.6 is shown in Figure
4.1a. It is clearly shown that a strong current sheet with narrow width is formed aroundthe X-line, supporting fast reconnection resulting from an intense magnetic dissipation
due to the electron inertia eﬀects (e.g., Vasyliunas, 1975; Hesse et al., 1999). Figure

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 67

(a)

7.7

t= 15.6

Jy

3.9

Z

0.0

-3.9

-7.7
30.7

-2.00

(b)

43.0

-1.78

-1.56

55.3

-1.34

-1.12

X

-0.90

67.6

-0.68

79.9

-0.46

-0.24

92.2

-0.02

0.20

7.7

E||
3.9

Z

0.0

-3.9

-7.7
30.7

-0.20

43.0

-0.16

-0.12

55.3

-0.08

-0.04

X

0.00

67.6

0.04

79.9

0.08

0.12

92.2

0.16

0.20

Figure 4.1: (a) Out-of-plane current density (Jy ) and (b) electric ﬁeld parallel to the
local magnetic ﬁeld (E ) at tωci = 15.6. The magnetic ﬁeld lines (white solid lines) aresuperposed in each ﬁgure.

4.1b presents a contour map of the electric ﬁeld parallel to the in-plane magnetic ﬁeld
described by the white solid lines, that is, E = E·b, where b = [(Bx , 0, Bz )/ Bx2 + Bz2 ]T .
It is found that wave-like structures of E are induced almost along the magnetic ﬁeld
lines around the edge of the plasma sheet, that is, in the plasma sheet-lobe boundary
region. It is worth noticing that the wave-active region is inside the magnetic separatrices
where the fast outﬂow of electrons accelerated near the X-line is expected to present.
The magniﬁcation of the right-hand plane of Figure 4.1b is shown in Figure 4.2.
In order to reveal the role of the waves, we have investigated the time evolutions of
the ve − x phase space distribution and E proﬁle along a magnetic ﬁeld line passing
through the wave-active region at tωci = 15.6, where ve = v e · b is the electron velocity
parallel to the in-plane magnetic ﬁeld. Figure 4.3 shows the snapshots of such proﬁles
at (a) tωci = 13.6, (b) 15.6, and (c) 17.2. We choose the ﬁeld line passing through(x/λi0 , z/λi0 ) = (84.0, −1.32) at each time, the trace of which is presented in Figure 4.2
by the white dashed line. The magnetic ﬁeld line at tωci = 13.6 is not yet reconnected. At
this time, the electrons locating along the ﬁeld line are slowly pulled into the diﬀusion
region due to E induced near the X-line. Counter-streaming electrons coming fromCHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 68

z/λi0

4

tωci = 15.6

E||

(a)
(b) (c)2 tω ci =13.6 15.6 17.2

0.2

0
-2
-4
65

70

x/λi0

75

80

-0.2

Figure 4.2: Magniﬁcation of Figure 4.1(b) for the right-hand side of the X-line.
White dashed lines represent the magnetic ﬁeld lines passing through (x/λi0 , z/λi0 ) =(84.0, −1.32) at (a) tωci = 13.6, (b) 15.6, and (c) 17.2.

the left-half plane of Figure 4.1b appear only near the X-line, and the wave activity is
weak in most region. On the other hand, the magnetic ﬁeld line has been reconnected attωci = 15.6 and passes through the downstream region of the fast electron ﬂow generated
near the X-line. At this time, a strong electron beam ﬂowing away from the plasma sheet
appears in the extensive region along the ﬁeld line and the wave activity becomes strong.
It is shown that some of electrons are trapped by the electrostatic potential of the waves
and form the electron holes in the ve − x phase space (see Figure 4.3b). At tωci = 17.2,
the magnetic ﬁeld line reaches further downstream region and the electrons residing
along the line are found to be heated parallel to the ambient magnetic ﬁeld, so thatthe waves are damped due to the Landau resonance with the heated electrons. The
electron heating due to the waves is also evident in Figure 4.4, which shows the electron
distribution function at (x/λi0 , z/λi0 ) = (71.8, −1.32) for each time presented in Figure
4.3. The electron distribution function at tωci = 13.6 consists of one component of acold isotropic electrons with a weak bulk velocity, Ve ≈ −VA , along the magnetic ﬁeld
line. At tωci = 15.6, in addition to this cold component, a strong beam component withVe ≈ 7VA and high perpendicular temperature has arisen. Then electrons are ﬁnallyheated selectively parallel to the ambient magnetic ﬁeld due to the scattering by thewaves and form the ﬂat-topped distribution function.
In order to clarify the heating mechanism, it is important to estimate the heatinglevel resulting from the electron trapping by the electrostatic potential wells along theﬁeld line. In this case, the trapping width vt can be a good benchmark for the heatinglevel, where vt means the absolute value of the velocity in the rest frame of the waves,below which electrons are trapped by the electrostatic potential wells. The trapping

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 69

E||/VAB0 ve||/VA(a) 20 tωci = 13.6
10
0
-10
-20
61.5
0.2
0
-0.2
-0.4
61.5

66.3

71.3

76.3

81.3

84.

66.3

71.3

76.3

81.3

84.

E||/VAB0 ve||/VA(b) 20 tωci = 15.6
10
0
-10
-20
65.1
0.2
0
-0.2
-0.4
65.1

69.6

74.6

79.6

83.9

69.6

74.6

79.6

83.9

E||/VAB0 ve||/VA(c) 20 tωci = 17.2
10
0
-10
-20
67.3
0.2
0
-0.2
-0.4
67.3

71.8

76.8

81.8

83.9

71.8

76.8
x/λi0

81.8

83.9

Figure 4.3: Electron velocity (ve ) and electric ﬁeld (E ) parallel to the in-plane magnetic
ﬁeld are plotted along the ﬁeld lines passing through (x/λi0 , z/λi0 ) = (84.0, −1.32) at(a) tωci = 13.6, (b) 15.6, and (c) 17.2, each of which is represented by the white dashedline in Figure 4.2.

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 70(a)10

t= 13.6
Log N
3.

vz/VA

5

(b)
10 t= 15.6

Log N
3.

5

(c)
10 t= 17.2

0

0

0

-5

-5

-5

-10

-10

-10

-10

-5

0

vx/VA

5

10

0.

-10

-5

5

0

10

0.

Log N
3.

5

-10

vx/VA

-5

0

5

10

0.

vx/VA

Figure 4.4: Electron distribution functions at the position (x/λi0 , z/λi0 ) = (71.8, −1.32)
and the time (a) tωci = 13.6, (b) 15.6, and (c) 17.2.

width vt is related to the amplitude of the wave potential φ0 by,
1
me vt2 ∼ 2eφ0 .
2

(4.2)

The potential amplitude φ0 is represented as,
φ0 =

E0
,
k

(4.3)

where we assume Ex (x) = E0 sin kx and Ex  E in the rest frame of the waves. Thus
vt is estimated as,

eE0
(4.4)
vt ∼ 2me k
In our simulations, the wave amplitude reaches E0 ≈ 0.40VA B0 and the typical wave
length is λw ≈ 1.3λi0 which corresponds to kλi0 ≈ 4.8 (see also Figure 4.7). Thus thetrapping width is calculated using (4.4) as vt ≈ 5.8VA , so that the shoulder width of the
ﬂat-topped electron distribution is estimated as 2vt ≈ 12VA , which is fairly consistent
with our simulation results (see Figure 4.4c). The shoulder energy of the ﬂat-topped
electrons can be also obtained from Esh ∼ me vt2 /2. Taking typical parameters in the
Earth magnetotail, B0 ∼ 20 nT and nps ∼ 0.3 cm−3 , the shoulder energy can be calculated to be Esh ∼ 1 keV, which is consistent with satellite observations (e.g., Shinohara
et al., 1998).
It is also valuable to estimate the bouncing time τb of trapped electrons in order toknow the time scale for the electron heating. We consider a single electron trapped by
the sinusoidal electrostatic potential. This electron obeys the equation of motion,d2 x
me 2 = −eE0 sin kx
dt −eE0 kx,(4.5)

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 71

where we assume a small amplitude of the electron displacement compared to the wave
length. Thus the bouncing frequency of the trapped electron is obtained as

ekE0
ωb =
,
me

(4.6)

so that τb is estimated to beτb ∼ ωb−1 =



me
.
ekE0

(4.7)

−1
Using the same parameters employed in the above discussion, we get τb ∼ 0.072 ωci,which indicates that the time scale for the electron heating is quite shorter than the−1
Alfvén transit time τA ∼ 1.0 ωcithat is thought to govern the structure change of thecurrent sheet associated with magnetic reconnection. Therefore such electron heating

should occur in the extensive region in the plasma sheet-lobe boundary region where theintense wave activity appears. This is conﬁrmed by Figure 4.5, which shows the proﬁlesof the electron temperatures (Te , Te⊥ , and Tey ) at tωci = 17.2 along the magnetic ﬁeldline traced in Figure 4.2. The temperatures are deﬁned by,me(vej − vej )2 f (v e ) d3 v,
(4.8)
Tej =
n
where,
1
vej f (v e ) d3v,
vej =
n

n = f (ve ) d3 v.

(4.9)(4.10)

Here, f (v e ) is the local distribution function of electrons, j (= , ⊥, and y) denotes theeach component, ve = v e · b with b = [(Bx , 0, Bz )/ Bx2 + Bz2 ]T , ve⊥ = v e · (b × êy ), andvey = v e · êy (êy is the unit vector directing toward +y). It is found that the ﬂat-topped
electrons, which have higher temperature in the parallel direction, widely appear along
the ﬁeld line.

4.3.2 Instability Responsible for the Electron Heating
In the previous subsection, we have found that the electron heating parallel tothe ambient magnetic ﬁeld occurs in the plasma sheet-lobe boundary region due to theelectron trapping by the electrostatic potential wells associated with the waves, so thatthe ﬂat-topped electron distribution is formed. In this subsection, we investigate how
the unstable waves are excited in the boundary region.

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 72

Te/miVA2

0.4

Te
Te
Te

0.3
0.2
0.1
0
67.3

71.8

76.8

81.8

83.9

x/λi0
Figure 4.5: Electron temperature proﬁles along the ﬁeld line passing through(x/λi0 , z/λi0 ) = (84.0, −1.32) at tωci = 17.2. Solid, dashed, and dotted lines denoteTe , Te⊥ , and Tey , respectively.

Figure 4.6 shows the time evolution of Ex proﬁle along the x axis at z/λi0 = −1.32,
where the waves are quite active (see Figure 4.2). It is found that the waves are propagating away from the diﬀusion region almost along the magnetic ﬁeld line. The phasevelocity is Vph ∼ 4VA in the early stage and increases as they propagate downstream.
These waves hardly propagate across the magnetic ﬁeld line (i.e., in the z direction in
the boundary region) and the perpendicular component is quite small, indicating thatthey are highly electrostatic.
The wave spectrum of Ex in the ω-k space is shown in Figure 4.7(contours). The
sampling region is 68.0 ≤ x/λi0 ≤ 75.6 at z/λi0 = −1.32 and the time interval is14.2 ≤ tωci ≤ 15.3, which is the region surrounded by the white dashed line in Figure 6
where the number of the trapped electrons are so small that the linear theory is expected
to be available. It is found that the spectrum has a clear peak around kλe0  0.48
and ω  0.11 ωpe0, where ωpe0 is the electron plasma frequency deﬁned by nps , and
λe0 = c/ωpe0 is the electron inertial length. Figure 4.3 and Figure 4.4 indicate that thestrong electron beam along the magnetic ﬁeld line seems to be responsible for the unstablewaves, providing the free energy. In this case, two kinds of instability are prospective.
One is the electron two-stream instability between the strong beam electrons and thealmost stationary electrons, and the other is the Buneman instability excited due to thevelocity diﬀerence between the strong beam electrons and the stationary ions.
The unstable waves excited by the electron two-stream instability have the following
dispersion relation at the linear stage,
1−

2
2
ωpe,f
ωpe,s
−
= 0,
(ω − kVs )2 (ω − kVf )2(4.11)

where ω = ωr + iγ, subscripts s and f denote the slow (or stationary) and fast beamelectrons, respectively, and Vj (j = s, f ) is the beam velocity of component j. In Figure

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 73

16.5

tωci

16

15.5

Ex
0.2
15

14.5
-0.2
68

70

72

74

76
x/λi0

78

80

82

Figure 4.6: Time evolution of Ex proﬁle along the x axis at z/λi0 = −1.32.

0.3

ωr

Ex

0.25

0.3

ω/ωpe0

0.2
0.15
0.1

γ

0.05

0.

0
0

0.2

0.4

0.6

0.8
kλe0

1

1.2

Figure 4.7: Wave spectrum of Ex (contours). The sampling region is 68.0 ≤ x/λi0 ≤ 75.6
at z/λi0 = −1.32, and the time interval is 14.2 ≤ tωci ≤ 15.3, which is the region
surrounded by the white dashed line in Figure 4.6. Bold solid and dashed lines are thetheoretical proﬁles of the wave frequency (ωr ) and growth rate (γ), respectively, deduced
from the linear analysis of the electron two-stream instability.

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 74

4.7, the wave frequency ωr and growth rate γ satisfying the dispersion equation (4.11) aresuperposed by the bold solid and dashed lines, respectively, for the case of Vs = −1.0 VA ,Vf = 7.0 VA , ne,s = 0.031 nps, and ne,f = 0.013 nps, where ne,j (j = s, f ) is the number
density of the electrons involved in component j. These theoretical proﬁles well agreewith the simulation results, that is, the wave spectrum obtained from our simulations
has an intense peak being consistent with a theoretical dispersion relation (bold solidline), where the growth rate from the theory (bold dashed line) reaches a maximum
value. There remains a slight ambiguity in the values of the adopted parameters (i.e.,
Vs , Vf , ne,s , and ne,f ), but the theoretical proﬁles are essentially unchanged even if theparameters are slightly shifted.
On the other hand, the linear dispersion relation of the Buneman instability is represented as,
1−

2
2
ωpe,f
ωpi
−
= 0,
ω2
(ω − kVf )2(4.12)

where ωpi is the ion plasma frequency, and we assume that the ion beam velocity is
negligible compared with the electron beam velocity. The unstable waves excited by2
the Buneman instability are expected to have ωr ∼ (ωpiωpe,f )1/3 = 0.037 ωpe0, and k ∼ωpe,f /Vf = 0.27 λ−1e0 , using ni = 0.044 nps , ne,f = 0.013 nps , and Vf = 7.0 VA , where ni is
the number density of ions (Buneman, 1958). Therefore the dispersion relation of theBuneman instability is inconsistent with the spectral peak obtained from our simulations.
Furthermore, the growth rate is estimated as γ ∼ ωr ∼ 0.037 ωpe0, which is smaller thanthat of the electron two-stream instability. Thus we can conclude that the instability
responsible for the electron heating should be the electron two-stream instability, andnot the Buneman instability. A weak peak that appears at the lower left of the intense
peak probably represents the Buneman instability.

4.3.3 Origin of the Electron BeamIn this subsection, we study how the strong electron beam with high perpendiculartemperature as shown in Figure 4.4b is produced. In order to investigate the origin of theelectron beam, we pick out the electrons forming the beam component with ve > 4.0
at tωci = 15.6 and trace their trajectories backward in time. Figure 4.8 shows the
positions of such electrons at tωci = 15.6 (ﬁlled circles) and 14.0 (blue triangles andred squares), and the magnetic ﬁeld lines at tωci = 14.0 (solid lines). Note that weclassify the electrons chosen at tωci = 15.6 into two categories, that is, the electrons
2
2 < v
+ vey
with low perpendicular energies ( ve⊥th,e0 ) and those with high perpendicular

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 75
4

z/λi0

2

tωci = 15.6
tωci = 14.0, ve 2+vey2 < vth,e0tωci = 14.0, ve 2+vey2 > 2vth,e0

0
-2
-4

50

55

60

x/λi0

65

70

75

Figure 4.8: Positions of electrons at tωci = 14.0 (triangles and squares) traced backwardin time from the positions denoted by the ﬁlled circles at tωci = 15.6. Here, we choose
the electrons forming the beam component with ve /VA > 4.0 at tωci = 15.6 as shownin Figure 4.4b. Of these having low electrons, triangles and squares represent those
2
2
2
2 >
+ vey
perpendicular energy ( ve⊥ + vey < vth,e ) and high perpendicular energy ( ve⊥2vth,e ) at tωci = 15.6, respectively. Magnetic ﬁeld lines at tωci = 14.0 are superposed.
2
2 < 2v
energies ( ve⊥+ vey
2Te,ps /me is the initial electron thermalth,e0 ), where vth,e0 =velocity and vth,e0 = 3.33 VA in the present case. The positions of the former electronsat tωci = 14.0 are indicated by the blue triangles and those of the latter electrons are by
the red squares in Figure 4.8.
We ﬁnd that the origin of the electrons with the low perpendicular energy is clearly
diﬀerent from that with the high perpendicular energy. Most electrons having the highperpendicular energy at tωci = 15.6 have passed through the region quite near the X-line.
Thus we think that those electrons are originating from the lobe region and accelerated
by the inductive electric ﬁeld Ey in the electron diﬀusion region and experience themeandering/Speiser motion (Speiser, 1965). Some of them reﬂect at the magnetic mirror
points depending on their pitch angles after they are ejected from the plasma sheet, andstart the bouncing motion around the central plasma sheet. Such electrons can gain
more kinetic energy by performing the ∇B drift and the curvature drift in the opposite
direction of Ey (Hoshino et al., 2001b) and/or by experiencing the betatron acceleration
resulting from the increase of the magnetic ﬁeld strength in the pileup region of the ﬁeld
lines (Birn et al., 2000). Both of these processes can contribute to the increase of theperpendicular energy of electrons. On the other hand, most electrons forming the cold
beam component with the low perpendicular energy at tωci = 15.6 are located in the
opposite side of the plasma sheet at tωci = 14.0. They are accelerated by E (see Figure4.1b) toward the plasma sheet along the magnetic ﬁeld line and quickly pass through theCHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 76

downstream region in the vicinity of the X-line where strong Ey is imposed. Thus theseelectrons arrive at the opposite side of the plasma sheet without experiencing signiﬁcantacceleration in the perpendicular direction.

4.3.4 Electrostatic Solitary Waves Evolved from the Beam Instability
It is known that the nonlinear evolution of the electron two-stream instabilities canresult in the generation of the ESW (e.g., Matsumoto et al., 1994; Kojima et al., 1994;
Omura et al., 1996), which have been observed in the dayside magnetopause (Matsumoto
et al., 2003) and the near-Earth magnetotail (Cattell et al., 2005) in association withmagnetic reconnection. In this subsection, we check whether the electron two-stream
instability excited in the plasma sheet-lobe boundary region evolves into the ESW. Figure
4.9a shows the scatter plot in the (−ve ) − x phase space along the x axis at z/λi0 = 1.44
and tωci = 17.2. It is found that the electron holes are clearly formed at x/λi0  74.6
and 77.2 associated with the bipolar structures of E presented in Figure 4.9b. Thesebipolar E structures result from the electrostatic potential, because the depressions in
the electron density relative to the ion density appear in Figure 4.9c, in which the ion
and electron densities are plotted by the dashed and solid lines, respectively. Thus we
can understand that the electron holes in the phase space represent the electron trappingin the electrostatic potential wells.
In Figure 4.10, we show the time evolution of Ex proﬁle along the x axis at the same
position in z as Figure 4.9. The bipolar structures seen in Figure 4.9b are found to be apart of the waves excited by the electron two-stream instability in the early stage beforetωci ≈ 15.6. After tωci ≈ 15.6, the bipolar structures are separated from the unstablewaves and start to propagate independently. We consider that the diﬀerence in the phasevelocity attributes to the diﬀerence in the density of the electron beam component. Thisidea comes from the fact that the phase velocity of the waves excited by the electron
two-stream instability increases as the density of the beam component decreases so thatthe plasma frequency in the rest frame of the electron beam is dropped. The density
drop in the beam component is caused by the decrease of the high energy particles, which
are trapped by the electrostatic potential wells and excluded from the beam component.
However, the density in the source region of the electron beam continues to increase as
the reconnected ﬁeld line moves from the X-line toward the pileup region of the ﬁeld
lines. Therefore the wave packets locating near the source region (i.e., near the plasma
sheet along the ﬁeld line) are dragged away from the main part of the unstable waves as
the density of the beam component increases.
The important point in Figure 4.10 is that the bipolar structures are stable for much

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 77-ve||/VA

20
10tωci = 17.2

(a)

0
-10

nZnps

E||/VAB0-20
0.4 (b) 70

72.5

75

77.5

80

82.5

85

72.5

75

77.5

80

82.5

85

0.2
0
-0.2
-0.4
0.12
(c) 70
0.1
0.08
0.06
0.04
0.02
70

Ion
Electron

72.5

75

77.5

80

82.5

85

x/λi0
Figure 4.9: (a) Parallel electron velocity (ve ), (b) parallel electric ﬁeld (E ), and (c) ion
and electron densities (dotted and solid lines, respectively) along the x axis at z/λi0 =1.44 and tωci = 17.2.

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 78
18

tωci

17

Ex

16

0.2
15

-0.2
68

70

72

74

76
x/λi0

78

80

82

Figure 4.10: Time evolution of Ex proﬁle along the x axis at z/λi0 = −1.44.

−1longer time compared with the electron bouncing time (∼ 0.1 ωci). Thus we identify

these bipolar structures with the ESW. The left packet of the ESW located at x/λi0 74.6 in Figure 4.9b, which seems to be more stable than the right one, is propagating
downstream with the velocity VESW  3 VA and the wave length λESW  1.5 λi0, which
are comparable with the wave properties excited from the electron two-stream instability.
The ESW have arisen in association with magnetic reconnection also in earlier simulations using a three-dimensional particle-in-cell code (Drake et al., 2003). However,
the generation mechanism of the ESW in their simulations is diﬀerent from those in our
simulations, in that their ESW are evolved from the Buneman instability. The ESW in
our simulations are also diﬀerent from the satellite observations reported by Cattell et al.
(2005) for the same reason. Nevertheless, the generation mechanism in our study can beone of the strong candidates to explain the ESW associated with magnetic reconnection,
because a large number of observations in the plasma sheet-lobe boundary region havedetected the ESW resulting from the electron two-stream instabilities (e.g., Matsumoto
et al., 1994; Kojima et al., 1994), which, we showed, are responsible for the formationof the ﬂat-topped electrons as often observed during magnetic reconnection (e.g., Saito
et al., 1995; Shinohara et al., 1998). We think that the Buneman instability dominantly
evolves in the regions where electrons are entirely ﬂowing relatively to ions so that the
number of the stationary electrons are negligibly small, such as in the electron diﬀusion
region. Though it is impossible to study such current-aligned instabilities in our two-

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 79

dimensional simulations, it is important to investigate the role of the instabilities in alarge three-dimensional system without the guide ﬁeld.

4.4

Summary and Disccusion

We have demonstrated the electron heating mechanism expected in the plasma sheetlobe boundary region speciﬁcally of the Earth magnetotail. We used the newly developed
2-1/2 dimensional PIC code with the AMR and particle splitting techniques, whichenables us to perform large-scale kinetic simulations including the extensive downstream
region of the X-line. Our main ﬁnding is that the electron two-stream instability between
the background cold electrons and the strong beam electrons with high perpendiculartemperature should be responsible for the formation of the ﬂat-topped electrons, whichhave been often observed in the plasma sheet-lobe boundary region in association withmagnetic reconnection. Both the bump-on-tail instability and the Buneman instability
do not play a signiﬁcant role in the electron heating process in the boundary region.
The heating mechanism of electrons is simple. Electrons are trapped by the electrostatic potential wells resulting from the electron two-stream instability and scatteredalong the ambient magnetic ﬁeld, so that the truncated (i.e., ﬂat-topped) distributionfunction is formed in the parallel direction. The amplitude of the electrostatic turbulence is large enough to produce the ﬂat-topped electrons having a high shoulder energy
up to a few keV. The heating process is quickly completed in comparison with the timescale of the dynamic change of the current sheet, so that, once magnetic reconnection
is triggered, the ﬂat-topped electrons should be observed in the extensive region of theboundary region. The strong beam electrons consist of two components originating fromdistinct regions. Some electrons come from the opposite boundary region of the plasma
sheet along the magnetic ﬁeld line, forming the cold component of the electron beam.
The others are originating from the lobe region, passing through the electron diﬀusion
region. These electrons are strongly accelerated in the perpendicular direction near the
X-line, forming the hot component of the electron beam. We have also revealed that
the ESW, which have been often observed in the plasma sheet-lobe boundary region, areevolved from the electron two-stream instability in association with magnetic reconnection. Thus we conclude that the electron two-stream instability presented in the current
study should be responsible not only for the formation of the ﬂat-topped electron distribution, but also for the generation of the ESW in the plasma sheet-lobe boundary region
when magnetic reconnection has occured.
Another important candidate responsible for the formation of the ﬂat-topped elec-

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 80

trons is the LHDI driven by the cross-ﬁeld current ﬂowing in the boundary region (Hubaet al., 1978; Shinohara et al., 1998). The unstable waves excited by the LHDI propagate
obliquely to the ambient magnetic ﬁeld, so that electrons can be scattered along themagnetic ﬁeld by the waves. Though it is impossible to study the eﬀects of the LHDI
by using our two-dimensional simulations because they omit wave modes propagatingin the direction of the sheet current (i.e., the y direction), it is possible alternatively to
speculate on the instability by resorting to the linear theory. As shown in Figure 4.7,
the growth rate of the electron two-stream instability is estimated as γET S ∼ 0.1 ωpe ,where ωpe is the local plasma frequency of electrons. On the other hand, the growth√rate of the LHDI can be approximated by γLHD ∼ 0.1 ωLH , where ωLH = ωci ωce (ωciand ωce are the local cyclotron frequencies of ions and electrons, respectively, here) isthe local lower hybrid frequency (Davidson et al., 1977). If the relateion ωpe ≈ ωce isassumed, γLHD ∼ me /mi γET S is derived. Therefore the growth rate of the LHDI ismuch smaller than that of the electron two-stream instability, which indicates that the
LHDI could not aﬀect the formation of the ﬂat-topped electrons so much as the electron
two-stream instability.
Slow-mode shocks can be also a candidate for the electron scattering in the parallel
direction (Schwartz et al., 1987). However, we could not obtain a clear structure of theslow-mode shocks in the present simulation within an allowed simulation time. We thinkthat more large-scale simulations lasting for longer time are necessary for the generationof the slow-mode shocks (Arzner and Scholer, 2001). This indicates that it takes much
longer time for the slow-mode shocks to aﬀect the electron heating compared with theelectron two-stream instability. Thus we suggest that the heating in the slow-mode
shocks would not be eﬀective for the formation of the ﬂat-topped electron distribution
in the vicinity of the X-line.
Finally, we make a comment on the net acceleration due to E that the electronsejected from the plasma sheet may experience during passing through the wave active
region along the magnetic ﬁeld line. Such acceleration could be one of the generation
mechanisms for high energy electrons in association with magnetic reconnection. In
Figure 4.11, we show the kinetic energy We that a representative electron can gain fromE while it moves from the center of the plasma sheet (z/λi0 = 0) along the magnetic
ﬁeld line at tωci = 15.6. Here, We is deﬁned by, s=lWe (x(l)) =(−e)E (s) ds,(4.13)

s=0

where s denotes the ﬁeld-aligned coordinate starting at the center of the plasma sheet,and l is the distance from the origin along the coordinate. The electron is strongly

We||/miVA2

CHAPTER 4. WAVE AND ELECTRON HEATING DURING RECONNECTION 81
0.1
0
-0.1
-0.2
-0.3
-0.4
-0.5
-0.6
65.1

69.6

74.6

79.6

83.9

x/λi0
Figure 4.11: Kinetic energy that an electron can gain from the electric ﬁeld while it
moves from the center of the plasma sheet (z/λi0 = 0) along the magnetic ﬁeld line
passing through (x/λi0 , z/λi0 ) = (84.0, −1.32) at tωci = 15.6 as shown in Figure 4.2.

decelerated by E during the fast ejection from the plasma sheet. In the wave active
region, the electron repeats acceleration and deceleration in accordance with the wave
structure, but the net energy gain through the region is absent. Therefore the net
acceleration due to E along the wave active region does not seem to occur in our
simulations. Nevertheless, this problem should be carefully considered. Recent study
by Hoshino (2005) has demonstrated that electrostatic waves having 10 times larger
amplitude than that in our simulations can be excited in association with a driven
reconnection. Such large waves may contribute to the net energy gain of the outﬂowing
electrons.

CHAPTER 5
Conclusions and Further Studies
So far, a large number of simulation studies using PIC codes have been devoted tounderstanding a nonlinear evolution of the current sheet associated with magnetic reconnection. However, it has been diﬃcult to describe a large-scale dynamics of magnetic
reconnection by the kinetic simulations due to limited computer resources. The mainunresolved issues in the reconnection processes include the triggering mechanism, the reconnection rate in the steady-state reconnection, and the energy transport mechanism.
The diﬃculty in approaching such issues lies in the fact that macro-scale structures canbe strongly aﬀected by localized micro-scale processes, in which the kinetic treatment ofplasma is required.
In order to overcome the diﬃculty, we ﬁrst developed a new electromagnetic kinetic
code that employs the AMR technique and the particle splitting algorithm to conventional PIC codes. The AMR technique subdivides and removes cells dynamically inaccordance with a reﬁnement criterion and enhances the spacial resolution locally. On
the other hand, the particle splitting algorithm divides particles locating in ﬁne cells sothat it conserves the mass per charge, the total momentum, energy, and distributionof function of particles, and the moments on grids. The particle splitting is needed inorder to control the number of particles per cell. Too small number of particles percell can lead to a signiﬁcant numerical noise. We performed some test simulations andcompared three runs without the AMR, with the AMR, and with both the AMR and theparticle splitting. These results indicate not only that the AMR and particle splitting
algorithms are successfully applied to the conventional PIC codes, but also that they
are quite eﬀective to achieve high-resolution simulations on the evolution of the current
sheet associated with magnetic reconnection. Thus the new code allows us to implementa large-scale kinetic simulation on magnetic reconnection.
One of our goals by using large-scale simulations is to investigate a long time evolutionof the diﬀusion region and to see what supports fast reconnection in a steady state. Our

CHAPTER 5. CONCLUSIONS AND FURTHER STUDIES

83

results showed that fast reconnection was quickly achieved in a fully kinetic system but
the system did not reach steady state and the reconnection rate went down signiﬁcantly,
even though the eﬀects of the periodic system were almost excluded. We conclude thatthe key process responsible for slowing magnetic reconnection is the extension of theelectron diﬀusion region associated with the evolution of the strong polarization electricﬁeld imposed in the electron inﬂow region. The polarization electric ﬁeld directing toward
the neutral sheet is caused by the inertia diﬀerence between ions and electrons, andenhanced by the meandering motions of the background ions. Because the hot ions
initially loaded in the Harris-type current sheet have a large velocity in the out-of-plane
direction, they can easily escape from the diﬀusion region due to the Lorentz force
before their meandering motions become dominant. Thus their impact on the strong
electrostatic ﬁeld is small. The polarization electric ﬁeld forces the inﬂow electrons tomove toward the out-of-plane direction by the E ×B drift and enhances the out-of-planecurrent density in the electron inﬂow region. The role of the current is to reduce thecurvature of the magnetic ﬁeld lines so that they are mostly parallel to the x direction
near the X-line, in such a way that the current density becomes uniform along the
upstream edge of the electron diﬀusion region. As a result, the electron meandering
region so the electron diﬀusion region extends along the x direction. In order to conﬁrm
the role of the polarization electric ﬁeld, we compare the simulation results with mi /me =1 and 100. It is found that (1) a steady-state reconnection is achieved in the system where
the polarization electric ﬁeld does not arise, (2) a large reconnection rate is obtainedeven in the system without the Hall eﬀects. It is suggested that other mechanismsthat are not included in the present model might play an important role in enhancing
the reconnection rate. One of the candidates is the anomalous resistivity arising from
nonlinear interactions between waves and particles. Indeed the Buneman-type instability
is expected to arise in the y direction, because there exists a large velocity diﬀerence
between ions and electrons beyond the electron thermal velocity in the electron diﬀusion
region. In order to investigate the role of the Buneman-type instability, it is necessary
to perform a large-scale kinetic simulation in the three-dimensional system.
Another goal is to reveal energy transport processes associated with magnetic reconnection, that is, where and how plasma is accelerated and heated. In the present study,
we paid attention to an electron heating mechanism in the plasma sheet-lobe boundary
region, where the truncated (i.e., ﬂat-topped) distribution function of electrons was often observed in the Earth magnetotail during magnetic reconnection. We found that theelectron two-stream instability between the background cold electrons and the strong
beam electrons with high perpendicular temperature should be responsible for the for-

CHAPTER 5. CONCLUSIONS AND FURTHER STUDIES

84

mation of the ﬂat-topped electrons. Electrons are trapped by the electrostatic potential
wells resulting from the electron two-stream instability and scattered along the ambient magnetic ﬁeld, so that the truncated distribution function is formed in the parallel
direction. We showed that the electron heating through this process occured in the extensive region of the boundary region and the heating level was consistent with satellite
observations. The strong beam electrons consist of two components originating fromdistinct regions. Some electrons come from the opposite boundary region of the plasma
sheet along the magnetic ﬁeld line, forming the cold component of the electron beam.
The others are originating from the lobe region, passing through the electron diﬀusion
region. Furthermore, we revealed that the ESW could be evolved from the electron twostream instability in association with magnetic reconnection. Thus such waves should
be observed in the boundary region during reconnection.
The important issue which we could not access through this study is the triggering
mechanism of magnetic reconnection. In the present study, magnetic reconnection is
initialized with a small perturbation in the magnetic ﬁeld, so that the resulting thinnercurrent sheet enhances the growth rate of the tearing instability (Furth et al., 1963; Drakeand Lee, 1977). However, since the current sheet width typically reaches a few Re in the
Earth magnetotail (e.g., Sergeev et al, 1993), it is diﬃcult for the tearing mode alone totrigger fast reconnection within a suﬃciently short time scale. Recent two-dimensional
simulations in the plane orthogonal to the initial magnetic ﬁeld (in which the tearing
instability does not arise) have revealed that the LHDI is excited at the edge of thecurrent sheet and it modiﬁes the initial structure, so that a Kelvin-Helmholtz instability(KHI) develops at much faster rate than predicted by linear theories based on the initial
equilibrium (Shinohara et al., 2001; Lapenta and Brackbill, 2002). Furthermore, recent
three-dimensional simulations have suggested that the modiﬁcation of the initial proﬁlecaused by the LHDI can allow the onset of magnetic reconnection and enhance thegrowth of the tearing instability (Lapenta et al., 2003; Scholer et al., 2003; Ricci et al.,2004). However, the system sizes employed in their simulations are not large enough for
the KHI (which has a MHD scale) to arise because of limited computer resources. Thus
it is an open question how the KHI aﬀects the tearing instability. Three-dimensionalversion of our code would enable us to perform 10 − 50 times larger-scale simulations
than the previous simulations. It is a promising tool to study the relationship between
the tearing instability, the LHDI, and the KHI, and a key to open a new era.

This thesis was written on the basis of the fallowing papers.
1. Fujimoto, K., and S. Machida (2005), Full particle simulation of the plasma sheet
using adaptive mesh reﬁnement (AMR) technique, Adv. Space Res., in press,
doi:10.1016/j.asr.2005.03.096.
2. Fujimoto, K., and S. Machida (2005), Electromagnetic full particle code with adaptive mesh reﬁnement technique: Application to the current sheet evolution, J.
Comput. Phys., in press, doi:10.1016/j.jcp.2005.10.003.
3. Fujimoto, K., and S. Machida (2005), A generation mechanism of electrostatic
waves and subsequent electron heating in the plasma sheet-lobe boundary region
during magnetic reconnection, J. Geophys. Res., submitted.

ACKNOWLEDGMENTS

86

Acknowledgments
The author would like to express grateful appreciation to Prof. S. Machida for usuful
discussions and ﬁnancial support. The author also wishes to thank all colleagues ofDepartment of Geophysics, Kyoto University for their valuable instruction and technicalsupport, and would like to express thanks to the staﬀs of the Academic Center for
Computing and Media Studies of Kyoto University for their technical advices to improve
the performance of the simulation code. Some calculations were performed by FUJITSU
HPC2500 installed at the Academic Center for Computing and Media Studies of Kyoto
University, the Information Technology Center of Nagoya University, and the ResearchInstitute for Sustainable Humanosphere (RISH) of Kyoto University. Finally, the author
acknowledges the Research Fellowships of the Japan Society for the Promotion of Science(JSPS) for Young Scientists. This work has been partly supported by the Grant-in-Aid
for JSPS Fellows.

References
Alfvén H. (1957), On the theory of comet tails, Tellus, 9, 92–96.
Arzner, K., and M. Scholer (2001), Kinetic structure of the post plasmoid plasma sheet
during magnetotail reconnection, J. Geophys. Res., 106, 3827–3844.
Asano, Y., T. Mukai, M. Hoshino, Y. Saito, H. Hayakawa, and T. Nagai (2004), Current
sheet structure around the near-Earth neutral line observed by Geotail, J. Geophys.
Res., 109, A02212, doi:10.1029/2003JA010114.
Axford, W. I., H. E. Petschek, and G. L. Siscoe (1965), Tail of the magnetosphere, J.
Geophys. Res., 70, 1231–1236.
Baker, D. N., and E. C. Stone (1976), Energetic electron anisotropies in the magnetotail:Identiﬁcation of open and closed ﬁeld lines, Geophys. Res. Lett., 3, 557–560.
Baker, D. N., T. I. Pulkkinen, J. Büchner, and A. J. Klimas (1999), Substorms: A global
instability of the magnetosphere-ionosphere system, J. Geophys. Res., 104, 14,601–
14,611.
Barnes, J., and P. Hut (1986), A hierarchical O(N log N) force-calculation algorithm,
Nature, 324, 446–449.
Baumjohann, W., G. Paschmann, and C. A. Cattell (1989), Average plasma properties
in the central plasma sheet, J. Geophys. Res., 94, 6597–6606.
Baumjohann, W., and R. A. Treumann (1997), Basic Space Plasma Physics, Imperial
College Press, London.
Berger, M. J., and J. Oliger, Adaptive mesh reﬁnement for hyperbolic partial diﬀerential
equations, J. Comput. Phys., 53, 484–512.
Berger, M. J., and P. Colella (1989), Local adaptive mesh reﬁnement for shock hydrodynamics, J. Comput. Phys., 82, 64–84.

REFERENCES

88

Birdsall, C. K. and A. B. Langdon (1991), Plasma physics via computer simulation,
Institute of Physics Publishing, London.
Birn, J., M. F. Thomsen, J. E. Borovsky, G. D. Reeves, and M. Hesse (2000), Particle
acceleration in the dynamic magnetotail, Phys. Plasmas, 7, 2149–2156.
Birn, J., J. F. Drake, M. A. Shay, B. N. Rogers, R. E. Denton, M. Hesse, M. Kuznetsova,
Z. W. Ma, A. Bhattacharjee, A. Otto, and P. L. Pritchett (2001), Geospace environmental modeling (GEM) magnetic reconnection challenge, J. Geophys. Res., 106,3715–3719.
Biskamp, D., E. Schwarz, and J. F. Drake (1997), Two-ﬂuid theory of collisionless magnetic reconnection, Phys. Plasmas, 4, 1002–1009.
Buneman, O. (1958), Instability, turbulence, and conductivity in current-carrying
plasma, Phys. Rev. Lett., 1, 8–9.
Cai, H. J., and L. C. Lee (1997), The generalized Ohm’s law in collisionless magnetic
reconnection, Phys. Plasmas, 4, 509–520.
Cattell, C., J. Dombeck, J. Wygant, J. F. Drake, M. Swisdak, M. L. Goldstein, W.
Keith, A. Fazakerley, M. André, E. Lucek, and A. Balogh (2005), Cluster observations of electron holes in association with magnetotail reconnection and comparison
to simulations, J. Geophys. Res., 110, A01211, doi:10.1029/2004JA010519.
Coroniti, F. V., L. A. Frank, D. J. Williams, R. P. Lepping, F. L. Scarf, S. M. Krimigis,
and G. Gloeckler (1980), Variability of plasma sheet dynamics, J. Geophys. Res., 85,
2957–2977.
Cowley, S. W. H. (1985), Magnetic reconnection, in Solar System Magnetic Fields, editedby E. R. Priest, pp. 121-155, Reidel, Dordrecht.
Chapman, S., and V. C. A. Ferraro (1931), A new theory of magnetic storms, Terr.
Magn. Atmos. Electr., 36, 77–97.
Daughton, W. (1999), The unstable eigenmodes of a neutral sheet, Phys. Plasmas, 6,1329–1343.
Davidson, R. C., N. T. Gladd, C. S. Wu, and J. D. Huba (1977), Eﬀects of ﬁnite plasma
beta on the lower-hybrid-drift instability, Phys. Fluids, 20, 301–310.

REFERENCES

89

Drake, J. F., and Y. C. Lee (1977), Kinetic-theory of tearing instabilities, Phys Fluids,
20, 1341–1353.
Drake, J., M. Swisdak, C. Cattell, M. Shay, B. Rogers, and A. Zeiler (2003), Formationof electron holes and particle energization during magnetic reconnection, Science, 299,
873–877.
Dreher, J., and R. Grauer (2005), Racoon:
work for hyperbolic conservation laws,
doi:10.1016/j.parco.2005.04.011.

A parallel mesh-adaptive frame-

Parallel Computing,31,

913–932,

Dungey, J. W. (1953), Conditions for the occurrence of electrical discharges in astrophysical systems, Phil. Mag., Ser. 7, 44, 725–738.
Dungey, J. W. (1961), Interplanetary magnetic ﬁeld and the auroral zones, Phys. Rev.
Lett., 6, 47–48.
Dungey, J. W. (1978), The history of the magnetopause regions, J. Atmos. Terr. Phys.,
40, 231–234.
Fujimoto, K., and S. Machida (2003), An electron heating mechanism in theoutﬂow region from the X-type neutral line, J. Geophys. Res., 108, 1349,
doi:10.1029/2002JA009810.
Fujimoto, K., and S. Machida (2005a), Full particle simulation of the plasma sheet
using adaptive mesh reﬁnement (AMR) technique, Adv. Space Res., in press,
doi:10.1016/j.asr.2005.03.096.
Fujimoto, K., and S. Machida (2005b), Electromagnetic full particle code with adaptive
mesh reﬁnement technique: Application to the current sheet evolution, J. Comput.
Phys., in press, doi:10.1016/j.jcp.2005.10.003.
Furth, H. P., J. Killeen, and M. N. Rosenbluth (1963), Finite-resistivity instabilities ofa sheet pinch, Phys. Fluids, 6, 459–484.
Gelato, S., D. F. Chernoﬀ, and I. Wasserman (1997), An adaptive hierarchical particlemesh code with isolated boundary conditions, Astrophys. J., 480, 115–131.
Giovanelli, R. G. (1947), Magnetic and electric phenomena in the sun’s atomosphereassociated with sunspots, Mon. Notices Roy. Astron. Soc., 107, 338–355.

REFERENCES

90

Gold, T. (1961), Present evidence concerning magnetic ﬁelds and particle ﬂuxes in thesolar system, Space Res., 2, 828–836.
Griebel, M., and G. Zumbusch (1999), Parallel multigrid in an adaptive PDE solver
based on hashing and space-ﬁlling curves, Parallel Computing, 25, 827–843.
Harris, E. G. (1962), On a plasma sheath separating regions of oppositely directed magnetic ﬁeld, Nuovo Cimento, 23, 115–121.
Hesse, M. and D. Winske (1994), Hybrid simulations of collisionless reconnection incurrent sheets, J. Geophys. Res., 99, 11,177–11,192.
Hesse, M. and D. Winske (1998), Electron dissipation in collisionless magnetic reconnection, J. Geophys. Res., 103, 26,479–26,486.
Hesse, M., K. Schindler, J. Birn, and M. Kuznetsova (1999), The diﬀusion region incollisionless magnetic reconnection, Phys. Plasmas, 6, 1781–1795.
Hesse, M, J. Birn, and M. Kuznetsova (2001), Collisionless magnetic reconnection: Electron processes and transport modeling, J. Geophys. Res., 106, 3721–3735.
Hones, E. W., Jr. (1973), Plasma ﬂow in the plasma sheet and its relation to substorms,
Radio Sci., 8, 979–990.
Hoshino, M., T. Mukai, T. Yamamoto, and S. Kokubun (1998), Ion dynamics in magnetic
reconnection: Comparison between numerical simulation and Geotail observations, J.
Geophys. Res., 103, 4509–4530.
Hoshino, M., K. Hiraide, and T. Mukai (2001a), Strong electron heating and nonMaxwellian behavior in magnetic reconnection, Earth Planets Space, 53, 627–634.
Hoshino, M., T. Mukai, T. Terasawa, and I. Shinohara (2001b), Suprathermal electron
acceleration in magnetic reconnection, J. Geophys. Res., 106, 25,979–25,997.
Hoshino, M. (2005), Electron surﬁng acceleration in magnetic reconnection, J. Geophys.
Res., 110, A10215, doi:10.1029/2005JA011229.
Huang, J., and L. Greengard (2000), A fast direct solver for elliptic partial diﬀerential
equations on adaptively reﬁned meshes, SIAM J. Sci. Comput., 21, 1551–1566.
Huba, J. D., N. T. Gladd, and K. Papadopoulos (1978), Lower-hybrid-drift wave turbulence in the distant magnetotail, J. Geophys. Res., 83, 5217–5226.

REFERENCES

91

Huba, J. D., and L. I. Rudakov (2004), Hall magnetic reconnection rate, Phys. Rev.
Lett., 93, 175003.
Jessop, C., M. Duncan, and W. Y. Chau (1994), Multigrid methods for N-body gravitational systems, J. Comput. Phys., 115, 339–351.
Karimabadi, H., D. Krauss-Varban, J. D. Huba, and H. X. Vu (2004), On magnetic
reconnection regimes and associated three-dimensional asymmetries: Hybrid, Hall-lesshybrid, and Hall-MHD simulations, J. Geophys. Res., 109, A09205.
Khokhlov, A. M. (1998), Fully threaded tree algorithms for adaptive reﬁnement ﬂuid
dynamics simulations, J. Comput. Phys., 143, 519–543.
Knebe, A., A. Green, and J. Binney (2001), Multi-level adaptive particle mesh
(MLAPM): A c code for cosmological simulations, Mon. Not. R. Astron. Soc., 325,
845–864.
Kojima, H., H. Matsumoto, T. Miyatake, I. Nagano, A. Fujita, L. A. Frank, T. Mukai,
W. R. Paterson, Y. Saito, S. Machida, and R. R. Anderson (1994), Relation between
electrostatic solitary waves and hot plasma ﬂow in the plasma sheet boundary layer:
GEOTAIL observations, Geophys. Res. Lett., 21, 2921–2924.
Kravtsov, A. V., A. A. Klypin, and A. M. Khokhlov (1997), Adaptive reﬁnement tree:
A new high-resolution N-body code for cosmological simulations, Astrophys. J. Suppl.
Ser., 111, 73–94.
Kuznetsova, M. M., M. Hesse, and D. Winske (1998), Kinetic quasi-viscous and bulk
ﬂow inertia eﬀects in collisionless magnetotail reconnection, J. Geophys. Res., 103,
199–213.
Kuznetsova, M. M., M. Hesse, D. Winske (2001), Collisionless reconnection supported
by nongyrotropic pressure eﬀects in hybrid and particle simulations, J. Geophys. Res.,
106, 3799–3810.
Langdon, A. B., and B. F. Lasinski (1976), Electromagnetic and relativistic plasma
simulation models, Methods Comput. Phys., 16, 327–366.
Lapenta, G. (2002), Particle rezoning for multidimensional kinetic particle-in-cell simulations, J. Comput. Phys., 181, 317–337.
Lapenta, G. and J. U. Brackbill (2002), Nonlinear evolution of the lower hybrid drift
instability: Current sheet thinning and kinking, Phys. Plasmas, 9, 1544–1554.

REFERENCES

92

Lapenta, G., J. U. Brackbill, and W. S. Daughton (2003), The unexpected role of thelower hybrid drift instability in magnetic reconnection in the three dimensions, Phys.
Plasmas, 10, 1577–1587, doi:10.1063/1.1560615.
Lin, R. P., and H. S. Hudson (1976), Nonthermal processes in large solar-ﬂares, Sol.
Phys., 50, 153–178.
Lin, R. P., S. Krucker, G. J. Hurford, D. M. Smith, H. S. Hudson, G. D. Holman, R. A.
Schwartz, B. R. Dennis, G. H. Share, R. J. Murphy, A. G. Emslie, C. Johns-Krull, and
N. Vilmer (2003), RHESSI observations of particle acceleration and energy release in
an intense solar gamma-ray line ﬂare, Astrophys. J., 595, L69–L76.
Lin, Y. and D. W. Swift (1996), A two-dimensional hybrid simulation of the magnetotail
reconnection layer, J. Geophys. Res., 101, 19,859–19,870.
Mandt, M. E., R. E. Denton, and J. F. Drake (1994), Transition to whistler mediatedmagnetic reconnection, Geophys. Res. Lett., 21, 73–76.
Matsumoto, H., H. Kojima, T. Miyatake, Y. Omura, M. Okada, I. Nagano, and M.
Tsutsui (1994), Electrostatic Solitary Waves (ESW) in the magnetotail: BEN wave
forms observed by GEOTAIL, Geophys. Res. Lett., 21, 2915–2918.
Matsumoto, H., X. H. Deng, H. Kojima, and R. R. Anderson (2003), Observation ofelectrostatic solitary waves associated with reconnection on the dayside magnetopause
boundary, Geophys. Res. Lett., 30, 1326, doi:10.1029/2002GL016319.
McPherron, R. L. (1970), Growth phase of magnetospheric substorms, J. Geophys. Res.,75, 5592–5599.
Miyashita, Y., S. Machida, A. Nishida, T. Mukai, Y. Saito, and S. Kokubun (1999),
GEOTAIL observations of total pressure and electric ﬁeld variations in the near andmid-distant tail associated with substorm onsets, Geophys. Res. Lett., 26, 639–642.
Nagai, T., M. Fujimoto, Y. Sato, S. Machida, T. Terasawa, R. Nakamura, T. Yamamoto,T. Mukai, A. Nishida, and S. Kokubun (1998), Structure and dynamics of magnetic
reconnection for substorm onsets with Geotail observations, J. Geophys. Res., 103,
4419–4440.
Nagai, T., I. Shinohara, M. Fujimoto, M. Hoshino, Y. Saito, S. Machida, and T. Mukai(2001), Geotail observations of the Hall current system: Evidence of magnetic reconnection in the magnetotail, J. Geophys. Res., 106, 25,929–25,949.

REFERENCES

93

Nakabayashi, J. and S. Machida (1997), Electromagnetic hybrid-code simulation of magnetic reconnection: Velocity distribution functions of accelarated ions, Geophys. Res.
Lett., 24, 1339–1342.
Nakamura, M. S., M. Fujimoto, and K. Maezawa (1998), Ion dynamics and resaltant velocity space distribution functions in the cource of magnetic reconnection, J. Geophys.
Res., 103, 4531–4546.
Omura, Y., H. Matsumoto, T. Miyake, and H. Kojima (1996), Electron beam instabilities
as generation mechanism of electrostatic solitary waves in the magnetotail, J. Geophys.
Res., 101, 2685–2697.
Parker, E. N. (1958), Dynamics of the interplanetary gas and magnetic ﬁelds, Astrophys.
J., 128, 664–676.
Parker, E. N. (1963), The solar-ﬂare phenomenon and the theory of reconnection andannihilation of magnetic ﬁelds, Astrophys. J., Suppl. Ser., 8, 177-212.
Petschek, H. E. (1964), Magnetic ﬁeld annihilation, in AAS-NASA Symposium on the
Physics of Solar Flares, NASA Spec. Publ. SP-50, 425–439.
Petschek, H. E. (1966), The mechanism for reconnection of geomagnetic and interplanetary ﬁeld lines, in The Solar Wind, edited by R. J. Mackin, and M. Neugebauer, pp.
257–273, Pergamon, Oxford.
Pei, W., R. Horiuchi, and T. Sato (2001), Ion dynamics in steady collisionless driven
reconnection, Phys. Rev. Lett., 87, 235003.
Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery (1992), Numerical
Recipes in FORTRAN, 2nd ed., Cambridge Univ. Press, New York.
Pritchett, P. L. (1994), Eﬀects of electron dynamics on collisionless reconnection in twodimensional magnetotail equilibria, J. Geophys. Res., 99, 5935–5941.
Pritchett, P. L. (2001), Geospace Environment Modeling magnetic reconnection challenge: Simulations with a full particle electromagnetic code, J. Geophys. Res., 106,3783–3798.
Ricci, P., J. U. Brackbill, W. Daughton, and G. Lapenta (2004), Inﬂuence of the lower
hybrid drift instability on the onset of magnetic reconnection, Phys. Plasmas, 11,4489–4500, doi:10.1063/1.1778744.

REFERENCES

94

Rostoker, G. (1972), Polar magnetic substorms, Rev. Geophys. Space Phys., 10, 157–211.
Saito, Y., T. Mukai, T. Terasawa, A. Nishida, S. Machida, M. Hirahara, K. Maezawa, S.
Kokubun, and T. Yamamoto (1995), Slow-mode shocks in the magnetotail, J. Geophys.
Res., 100, 23,567 – 23,581.
Sato, T. (1979), Strong plasma acceleration by slow shocks resulting from magnetic
reconnection, J. Geophys. Res., 84, 7177–7190.
Sato, T., and T. Hayashi (1979), Externally driven magnetic reconnection and a powerfulmagnetic energy converter, Phys. Fluids, 22, 1189–1202.
Sergeev, V. A., D. G. Mitchell, C. T. Russell, and D. J. Williams (1993), Structure of
the tail plasma/current sheet at ∼ 11 RE and its changes in the course of a substorm,
J. Geophys. Res., 98, 17,345–17,365.
Scholer, M., and D. Roth (1987), A simulation study on reconnection and small-scale
plasmoid formation, J. Geophys. Res., 92, 3223–3233.
Scholer, M., I. Sidorenko, C. H. Jaroschek, R. A. Treumann, and A. Zeiler (2003), Onsetof collisionless magnetic reconnection in thin current sheets: Tree-dimensional particle
simulations, Phys. Plasmas, 10, 3521–3527, doi:10.1063/1.1597494.
Schwartz, S. J., M. F. Thomsen, W. C. Feldman, and F. T. Douglas (1987), Electron
dynamics and potential jump across slow mode shocks, J. Geophys. Res., 92, 3165–
3174.
Shay, M. A., and J. F. Drake (1998), The role of electron dissipation on the rate ofcollisionless magnetic reconnection, Geophys. Res. Lett., 25, 3759–3762.
Shay, M. A., J. F. Drake, R. E. Denton, and D. Biskamp (1998), Structure of thedissipation region during collisionless magnetic reconnection, J. Geophys. Res., 103,
9165–9176.
Shay, M. A., J. F. Drake, B. N. Rogers, and R. E. Denton (1999), The scaling of collisionless, magnetic reconnection for large systems, Geophys. Res. Lett., 26, 2163–2166.
Shay, M. A., J. F. Drake, and B. N. Rogers (2001), Alfvénic collisionless magnetic
reconnection and the Hall term, J. Geophys. Res., 106, 3759–3772.
Shay, M. A., J. F. Drake, M. Swisdak, and B. N. Rogers (2004), The scaling of embedded
collisionless reconnection, Phys. Plasmas, 11, 2199–2213.

REFERENCES

95

Shinohara, I., T. Nagai, M. Fujimoto, T. Terasawa, T. Mukai, K. Tsuruda, and T.
Yamamoto (1998), Low-frequency electromagnetic turburence observed near the substorm onset site, J. Geophy. Res., 103, 20,365–20,388.
Shinohara, I., H. Suzuki, M. Fujimoto, and M. Hoshino (2001), Rapid large-scale
magnetic-ﬁeld dissipation in a collisionless current sheet via coupling between KelvinHelmholtz and lower-hybrid drift instabilities, Phys. Rev. Lett., 87, 095001.
Sonnerup, B. U. Ö. (1979), Magnetic ﬁeld reconnection, in Solar System Plasma Physics,
vol. III, edited by L. J. Lanzerotti, C. F. Kennel, and E. N. Parker, pp. 45–108, NorthHolland, Amsterdam.
Speiser, T. W. (1965), Particle trajectries in model current sheets 1: Analytical solutions,J. Geophys. Res., 70, 4219–4226.
Spitzer, L., Jr. (1962), Physics of Fully Ionized Gases, 2nd ed., p. 139, Interscience, NewYork.
Suisalu, I., and E. Saar (1995), An adaptive multigrid solver for high-resolution cosmological simulations, Mon. Not. R. Astron. Soc., 274, 287–299.
Sweet, P. A. (1958), The neutral point theory of solar ﬂares, in Electromagnetic Phenomena in Cosmical Physics, edited by B. Lehnert, pp. 123–134, Cambridge University
Press, London.
Sweet, P. A. (1961), Mechanisms of solar ﬂares, Ann. Rev. Astron. Astrophys., 7, 149–
176.
Swisdak, M., J. F. Drake, M. A. Shay, and J. G. McIlhargey (2005), Transition fromantiparallel to component magnetic reconnection, J. Geophys. Res., 110, A05210,
doi:10.1029/2004JA010748.
Tanaka, M. (1995), Macro-particle simulations of collisionless magnetic reconnection,
Phys. Plasmas, 2, 2920–2930.
Terasawa, T., and A. Nishida (1976), Simultaneous observations of relativistic electron
bursts and neutral-line signatures in the magnetotail, Planet. Space Sci., 24, 855–866.
Terasawa, T. (1983), Hall current eﬀect on tearing mode instability, Geophys. Res. Lett.,
10, 475–478.

REFERENCES

96

Teyssier, R. (2002), Cosmological hydrodynamics with adaptive mesh reﬁnement A newhigh resolution code called RAMSES, Astronomy and Astrophysics, 385, 337–364.
Ugai, M. (1993), Computer studies on the fast reconnection mechanism in a sheared ﬁeld
geometry, Phys Fluids B, 5, 3021–3031.
Vasyliunas, V. M. (1975), Theoretical models of magnetic ﬁeld line merging, 1, Rev.
Geophys., 13, 303–336.
Vay, J. L., P. Colella, A. Friedman, D. P. Grote, P. McCorquodale, and D. B. Seraﬁni (2004), Implementations of mesh reﬁnement schemes for particle-in-cell plasma
simulations, Comput. Phys. Comm., 164, 297–305.
Villasenor, J., and O. Buneman (1992), Rigorous charge conservation for local electromagnetic ﬁeld solvers, Comput. Phys. Comm., 69, 306–316.
Villumsen, J. V. (1989), A new hierarchical particle-mesh code for very large scale cosmological N-body simulations, Astrophys. J. Suppl. Ser., 71, 407–431.
Yahagi, H., and Y. Yoshii (2001), N-body code with adaptive mesh reﬁnement, Astrophys. J., 558, 463–475.