Master Thesis

Possible Plasma Instabilities andElectron Heatings in the
Downstream Region of the X-type
Neutral Line

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

January 2003

Contents
1 Introduction

1

2 Possible Instabilities : Linear Theory

4

2.1 Local Dispersion Relation . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2 Drift Velocity Dependence . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Mass Ratio Dependence of the Critical Velocity . . . . . . . . . . . . . .

6
7

3 Anomalous Electron Heating10

3.1 Distribution Function of Electrons . . . . . . . . . . . . . . . . . . . . . .
3.2 Heating Level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
12

3.3 The Case of Massive Ion with mi /me = 400 . . . . . . . . . . . . . . . .

14

4 Discussion
4.1 Travel Time and Growth Time . . . . . . . . . . . . . . . . . . . . . . . .

16
16

4.2 Comparison with Observations . . . . . . . . . . . . . . . . . . . . . . . .

17

5 Summary and Conclusion

19

A Derivation of the Kinetic Plasma Dispersion Relation

22

A.1 Perturbed Distribution Function δf . . . . . . . . . . . . . . . . . . . . .
A.1.1 Ion Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . .

22
24

A.1.2 Electron Contribution . . . . . . . . . . . . . . . . . . . . . . . .

25

A.2 Dispersion Relation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

B Algorithm in the Full Particle Code29

B.1 Fundamental Equations . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

B.2 Boundary and Initial Conditions . . . . . . . . . . . . . . . . . . . . . . .

31

B.3 Momentum Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

i

B.4 Charge Assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

B.5 Poisson’s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

B.6 Faraday’s Law and Ampère’s Law . . . . . . . . . . . . . . . . . . . . . .

34ii

Abstract
A linear analysis and 2-1/2 dimensional full particle simulations are performed toinvestigate mechanisms of electron heating due to a strong Hall current, which is caused
by a large velocity diﬀerence between electrons and ions in the outﬂow region insidethe diﬀusion region associated with the magnetic reconnection process. The numericalsolution of the kinetic dispersion relation indicates that the prospective unstable mode
varies according to the electron-ion relative velocity Vd ; the kinetic cross-ﬁeld streaming
instability (KCSI) is dominant when Vd is under the substantial fraction of the electron
thermal velocity ve , while the electron cyclotron drift instability (ECDI) is ﬁrst excited
for the case where Vd exceeds that critical velocity. The critical drift velocity, where
the peak growth rate of the ECDI exceeds that of the KCSI, decreases as the mass
ratio increases, since the peak growth rate is almost independent of the mass ratio for
the ECDI while it decreases as the mass ratio increases for the KCSI. Our simulation
results show that electrons are heated parallel to the ambient magnetic ﬁeld when the
KCSI is excited, because they can be easily trapped along the magnetic ﬁeld line by
the obliquely propagating waves associated with the KCSI. On the other hand, we also
show that electrons are heated in the perpendicular direction due to the trapping by
the perpendicularly propagating waves assciated with the ECDI. In particular, electrons
are quickly heated up to a high temperature when the ECDI is dominant. The electron
heating by the ECDI is, however, ineﬀective when Vd /ve < 1 because the number ofresonant electrons, which act to dump the waves, is large in this case.
According to observational data obtained by spacecraft, both instabilities are likely
to be excited so that electrons can be heated in either parallel or perpendicular direction
depending on the controlling parameters. The identiﬁcation of the unstable mode will
be possible by analyzing the frequencies of excited waves observed by spacecraft.

Chapter 1
Introduction
Plasma heating in space and laboratory experiments is one of the long-standing problems in space physics as well as plasma physics. Magnetic reconnection is, in particular,
believed to play an important role in a fast conversion process of the magnetic energy
to plasma kinetic energy and be responsible for the plasma heating in the plasma sheetof the Earth. However, the detailed heating mechanism in reconnection is still poorly
understood.
In a collisionless reconnection process, the ideal magnetohydrodynamic (MHD) condition breaks down near the X-type neutral line and the diﬀusion region develops a
two-scale structure associated with the electron and ion scales, that is, the electron inertial length λe and the ion inertial length λi . Outside the diﬀusion region both electrons
and ions are frozen-in to the ambient magnetic ﬁeld and move together. (In this paper
we deﬁne the diﬀusion region as the ion diﬀusion region otherwise mentioned.) Thisimplies that the studies in the MHD framework can be applied. Indeed, the studies withhybrid simulations [e.g., Shay et al., 1998] and full particle simulations [e.g., Hesse et al.,1999] suggest that the outﬂow velocity of ions out of the diﬀusion region never exceed
the Alfvén velocity in the inﬂow region, as expected by Sweet [1958] and Parker [1963],
who predicted it with the use of a resistive MHD description. Furthermore, this fact isalso conﬁrmed by the observations with the Geotail spacecraft [Nagai et al., 1998]. On
the other hand, inside the diﬀusion region but outside the electron diﬀusion region, thediﬀerence in motion between electrons and ions becomes important. Since the inertial
length of ions λi is much larger than that of electrons λe , ions are easily unmagnetizedeven when electrons still remain frozen-in to the magnetic ﬁeld. This diﬀerence allows
the electron and ion motion to decouple. In the inﬂow region, electrons can continue toaccelerate with the E × B drift motion toward the neutral line until they decouple from1the magnetic ﬁeld in the electron diﬀusion region. Ions, however, have large gyroradii
compared to the scale size of the diﬀusion region so that they don’t execute the E × B
drift motion. This relative motion produces currents, that is, the Hall current system,in the vicinity of the X-type neutral line [Sonnerup, 1979]. Sonnerup [1979] ﬁrst predicted the presence of the Hall current system associated with the magnetic reconnection
theoretically in which the four current loops resulted in a quadrupole structure in theout-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 magnetotail of the earth, using magnetic ﬁeld
and plasma measurements onboard the Geotail spacecraft.
In the outﬂow region just outside the electron diﬀusion region, where the Hall currents
are inward into the neutral line associated with the outward electron ﬂow, electrons arefrozen-in to the magnetic ﬁeld while ions are unmagnetized as mentioned above. Thisimplies that a ﬂuid treatment is still available only to the electrons. Applying the analysis
due to Sweet [1958] and Parker [1963] to the electron diﬀusion region, we can estimate theupper limit of the outﬂowing electron velocity as “the electron Alfvén velocity” deﬁned√
as VAe ≡ BL / µ0 me nps [e.g., Shay et al., 2001; Hoshino et al., 2001a], where BL is
the magnetic ﬁeld in the inﬂow region, nps is the plasma density in the plasma sheet(assuming the quasi neutrality, nps = ni,ps  ne,ps , where suﬃx i and e denote ions andelectrons respectively), me is the mass of an electron, and µ0 is the magnetic permeability
in vacuum. Thus inside the diﬀusion region electrons can ﬂow with very high velocity
beyond the ordinary Alfvén velocity with which both electrons and ions are expected to
ﬂow outside the diﬀusion region. In fact, Shay et al. [2001] conﬁrmed by conducting
hybrid simulation that electrons could be accelerated up to approximately the electron
Alfvén velocity just outside the electron diﬀusion region independent of the electron mass.
Since the electron motion is basically described by the E × B drift, it is expected to exista large relative velocity between electrons and ions across the magnetic ﬁeld. In sucha situation, Hoshino et al. [2001b] suggested that the Buneman-like instability might2

be excited and electrons would be heated by the unstable waves. However they havenot conducted detailed quantitative analyses and it is necessary to identify what kindof instability mode will be actually excited, investigating the kinetic plasma dispersion
relation, so as to know the mechanism of the electron heating.
In this paper, we consider locally the outﬂow region just outside the electron diﬀusion
region and study the possible instabilities excited by a cross-ﬁeld relative motion withsuper Alfvénic velocity between electrons and ions, solving the kinetic plasma dispersion
relation. We next study the anomalous heating of electrons, using 2-1/2 dimensional
electromagnetic full particle simulation with high spatial and time resolution. We showthat the expected instability and the resultant mechanism of the electron heating vary
according to the value of the cross-ﬁeld relative velocity between electrons and ions.

3

Chapter 2
Possible Instabilities : Linear Theory
In order to identify the instabilities excited by a cross-ﬁeld relative motion between
electrons and ions just outside the electron diﬀusion region, we conduct a linear analysis.
We assume the frequency of the unstable waves which are expected to be higher than
the ion cyclotron frequency so that ions are unmagnetized. The lower-frequency modes
than the ion cyclotron frequency will not have enough time to grow to the appreciablelevels inside the diﬀusion region. Thus the candidates of the unstable modes excited bya cross-ﬁeld relative drift between electrons and ions are the kinetic cross-ﬁeld streaming
instability (KCSI) [e.g., Lemons and Gary, 1977; Wu et al., 1983; Winske et al., 1985],
the lower hybrid drift instability (LHDI) [e.g., Krall and Liewer, 1971; Davidson et al.,
1977; Huba et al., 1978], the electron cyclotron drift instability (ECDI) [e.g., Wong, 1970;
Gary and Sanderson, 1970; Forslund et al., 1970], and the ion acoustic instability [e.g.,Gary, 1970; Coroniti et al, 1977]. However, because we consider the local region so thatthe eﬀects of gradients in the magnetic ﬁeld magnitude and density are negligible, the
LHDI essentially doesn’t arise from our linear analysis. The ion acoustic instability is
unlikely to occur in the plasma sheet, where the ion temperature is higher than theelectron temperature [e.g., Baumjohann et al., 1989] so that the ion acoustic waves are
dumped. Thus, in the following study, we concentrate our attention on the two modes,
the KCSI and the ECDI.

2.1

Local Dispersion Relation

To derive a general dispersion relation, in which not only electromagnetic eﬀects areincluded but also the arbitrary direction of wave propagation and the electron-cyclotron
resonance are considered, we follow the procedure by Tsai et al. [1983]. We assume wave
4

frequencies as ωuh  |ω|  ωlh >> ωci in which ions are unmagnetized while electrons2
2 1/2
are still magnetized, where ωuh = (ωpe+ ωce
) is the upper hybrid frequency (ωpe andωce are the plasma frequency and the cyclotron frequency of electrons respectively.),ωlh = (ωce ωci )1/2 is the lower hybrid frequency, ωci is the ion cyclotron frequency, andω = ωr + iγ is a complex frequency. We also assume that the electron motion is basicallydescribed by the E × B drift as V d = −E0/B0 êx while ions are not drifting, where
E0 and B0 are the ambient 0th-order electric and magnetic ﬁelds directed to the y andthe −z directions respectively, and êx is the unit vector directed to the x direction.
Strictly speaking, ions are expected to conduct the meandering motion and drift to they direction in the vicinity of the plasma sheet. However the velocity is at most about
0.2 VA according to the results of a hybrid simulation [Shay et al., 1998], where VA is
the Alfvén velocity in the inﬂow region, and therefore much lower than the electron drift
velocity which is expected to be the substantial fraction of the electron Alfvén velocity
VAe [Hesse et al., 1999; Shay et al., 2001]. Thus we neglect the ion drift motion in thefollowing analysis. We set the wave number vector k in the x − z plane (k = k⊥ êx + k êz )so as to consider the oblique propagation to the magnetic ﬁeld. In addition, we assume2k⊥> k2 >> (∂ ln B0/∂x)2 , (∂ ln n/∂x)2, which justiﬁes the use of the homogeneity
approximation in B0 and n, where n is the local plasma density assuming n = ni  ne .

We calculate the local dispersion relation in the electron frame so that the uniform
electric ﬁeld E0 is canceled out and ions possess a bulk velocity Vd êx . Thus the unperturbed distribution functions of electrons and ions are described as

 2
3

vx + vy2 + vz2
1 2
exp −
fe0 = n
,
πve2
ve2


3

(vx − Vd )2 + vy2 + vz2
1 2
exp −
,
fi0 = n
πvi2
vi2

(2.1)
(2.2)

respectively, where vj = (2Tj /mj )1/2 is the thermal velocity, Tj and mj are the temperature and the mass of a particle of species j (=e, i), respectively. Then there exists
the following relation between the dielectric tensor ε(k, ω) and the ﬂuctuating electric ﬁeld δ Ê (see Appendix A for detailed procedure of the derivation), where δE =
δ Ê exp[i(k⊥ x + k z − ωt)],
D(k, ω) · δ Ê = 0,D(k, ω) ≡ ε(k, ω) −

5

ck
ω

2 
I−

kk
k2

(2.3)


,

where c is the velocity of light, and I is the unit tensor. We can derive the elementsof the dielectric tensor ε(k, ω) by combining the linearized Vlasov equation and theMaxwell equations (see Appendix A). The dispersion relation which gives the complex
eigenfrequency is deﬁned by det[D] = 0, and can be manipulated into the following form,
2
2k⊥
k2
2
2
D12 D13 D23 + D12
D22 − k⊥2 D23
D22
k2
D11 +
= 0.
2
k
2
D22 D33 + k⊥2 D23(2.4)

The elements Dij are deﬁned as


∞
2
2

2ωpi
2ωpe
D11 = 1 + 2 2 [1 + ξi Z(ξi )] + 2 2 1 + ξ0 e−µIn Z(ξn ) ,
k vi
k ve
n=−∞
∞  2
2

n
c2 k 2 2ωpe−µ
I + (In − In ) ξ0 Z(ξn ),D22 = 1 − 2 + 2 µe
2 n
ω
ω
2µ
n=−∞



∞ 
2
2
2 2

ω
c2 k 2 2ωpe
k
k⊥
ce
1−n
In Z(ξn ) ,D33 = 1 − 2 + 2 2 2 1 + ξ0 e−µ2ωk ve k
ω k⊥
n=−∞
∞
2
ωpe
k⊥ −µ D12 = −D21 = ie(In − In )Z(ξn ),kveωce k
n=−∞



∞ 
2
2

2ωpe
k
k2
k2
ω
ce
1 + ξ0 e−µ1−n
In Z(ξn ) ,D13 = 2 D31 = − 2 22k⊥k ve k
ω k⊥
n=−∞

∞ 
2
ωpe
ωce k 2k2k⊥ −µ 1−n
(In − In )Z(ξn ),e
D23 = − 2 D32 = i2k⊥kveωce k
ω
k
⊥
n=−∞

(2.5)

2 2
2
where ξn = (ω − nωce )/k ve , ξi = (ω − k⊥ Vd )/kvi , µ = k⊥ve /2ωce, In = In (µ) is the2
modiﬁed Bessel function of the nth order, In = dIn /dµ, θ = tan−1 (k⊥ /k ), k = k⊥+ k2,√∞and Z(ξ) = ( π)−1 −∞ dx exp[−x2]/(x − ξ) is the plasma dispersion function [Fried andConte, 1961]. In Equation (2.4), the ﬁrst term represents the electrostatic eﬀects while
the electromagnetic contributions are lumped into the second term.

2.2

Drift Velocity Dependence

Though many satellite observations in the plasma sheet associated with the magnetic
reconnection have been reported in the last decade [e.g., Baumjohann et al., 1991; Sergeev
et al., 1995; Shinohara et al., 1998; Nagai et al., 1998, 2001], there is still ambiguity about
the values of plasma parameters in the vicinity of the diﬀusion region and it is diﬃcult
6

to precisely predict the electron Alfvén velocity in the inﬂow region (the upper limit ofthe drift velocity in the outﬂow region just outside the electron diﬀusion region) or theelectron thermal velocity in the diﬀusion region. Thus it is necessary to examine how themost unstable mode varies according to the relative drift velocity between electrons andions. We ﬁrst solve the dispersion relation (2.4) numerically for the various values of thedrift velocity. For calculations we select the values of plasma parameters, ωce /ωpe = 0.5,
Ti /Te = 8.0, mi /me = 1836, and βi = 1.0. Figure 1 shows the numerical solutions of thedispersion relation for the various values of the drift velocity. In Figure 1a, the growth
rates γ maximized over wavenumber for a ﬁxed value of the propagation angle θ areplotted as a function of θ, for the ECDI (solid curves) and the KCSI (dashed curves),
and Figure 1b represents the real frequencies of unstable waves corresponding to Figure
1a. As can be seen in Figure 1a, the peak growth rates of both modes become larger as
the drift velocity increases but the enhancement of the growth rates is more rapid for the
ECDI than the KCSI. Thus, though the peak growth rate of the KCSI is larger than thatof the ECDI at Vd /ve = 0.4, the relation of the magnitudes reverses around Vd /ve = 0.6.
Figure 1b indicates that the frequency range of the ECDI and the KCSI is very much
diﬀerent; the KCSI is excited in a frequency a few times larger than the lower hybrid
frequency (ωlh /ωpe = 0.012 in the current parameters), while the ECDI is in a frequency
a little above the electron cyclotron frequency (ωce /ωpe = 0.5). These results imply thatthe high frequency mode, the ECDI, dominates the low frequency mode, the KCSI, in
the regime where the electron-ion relative velocity exceeds the substantial fraction of theelectron thermal velocity.
We also examined the drift velocity dependence in the case of mi/me = 100 for
application to the numerical simulations shown in the next chapter. In Figure 2, thesolutions of Equation (2.4) is shown for the case of ωce /ωpe = 0.5, Ti /Te = 8.0, mi /me =100, and βi = 1.0. Similar to the case of Figure 1, the peak growth rate of the ECDI
becomes larger rapidly as the drift velocity increases, while that of the KCSI shows rather
ﬂat variation with a little enhancement. In this case the critical drift velocity, where thepeak growth rate of the ECDI exceeds that of the KCSI, is about 1.0 ve , which is larger
than 0.6 ve in the case of mi /me = 1836.

2.3

Mass Ratio Dependence of the Critical Velocity

In the previous section, we found that the value of the critical drift velocity wasdependent on the mass ratio between an ion and electron, mi /me . Since in the pre7

vious section we considered only the cases of mi/me = 1836 (the real mass ratio) andmi /me = 100 (for the purpose of application to the simulational study conducted in thenext chapter), in this section we derive the general relation between the critical drift
velocity and the mass ratio. We assume unstable waves associated with the ECDI to beelectrostatic [e.g., Wong, 1970; Gary and Sanderson, 1970]. Thus the dielectric tensor
ε(k, ω) can be degenerated into a scalar as


∞
2
2

2ωpi
2ωpe
εL (k, ω) = D11 = 1 + 2 2 [1 + ξi Z(ξi )] + 2 2 1 + ξ0 e−µIn Z(ξn ) .
k vi
k ve
n=−∞

(2.6)

The unstable waves in association with the ECDI have a little higher frequencies thanthe upper harmonics of the electron cyclotron frequency (viz. ωr > nωce ) and propagate
almost perpendicularly to the ambient magnetic ﬁeld (viz. k  0). These characteristicsof the unstable waves and the condition (ωce /ωpe )2 << 1 allow us to set ξn >> 1 and√µ >> 1. Making use of the relation e−µIn (µ)  1/ 2πµ for the large µ approximation,
we can derive the growth rate of the ECDI as
γECDI 1
ImZ  (ξi )
2 1 [ReZ(ξi ) + ξi ReZ  (ξi )] + √ 1

2
nωce
Ti
πkve Te (ωr −nωce )2

kvi

,

(2.7)

where Z  (ξi ) = dZ(ξi )/dξi = −2[1 + ξi Z(ξi )]. The most unstable wave with the ECDI
has frequency around ωr  k(Vd − 0.7vi) [Forslund et al., 1972], which leads to ξi  −0.7.
Since it is the fundamental mode (n = 1) that has the peak growth rate according to theresults obtained in the previous section, we assume kVd  ωce . Then Equation (2.7) canbe manipulated into more convenient form, and the peak growth rate of the ECDI over
ω − k plane is expressed as the following form,  
 2  3/2  1/2 −1
M
3
γECDI
Vd
Vd
Te
3
Te
mi
 √
.
1+ √
ωce
4 π veTi
2 π ve
Ti
me

(2.8)

On the other hand, the peak growth rates associated with the KCSI are substantial
fractions of the lower hybrid frequency as shown in Figure 1a and 2a, and they are
assumed to be
M
γKCSI
=ξ
ωce



me
mi

1/2
where 0 < ξ < 1.

(2.9)

We deﬁne the value of the parameter ξ as a function of the mass ratio mi /me by thelinear interpolation between mi /me = 100 and 1836, where the numerical solutions of8

the growth rates are obtained in the previous section for the several values of the drift
velocity. We set ξ = 0.4 at mi /me = 100 and ξ = 0.8 at mi /me = 1836. By equating theMMpeak growth rate of the ECDI with that of the KCSI, that is, γECDI= γKCSI, we canlead the critical drift velocity as



 
 
 2
√  1/2  1/2
Ti  12 π Ti
Vd
1
me
.

−
(2.10)
−
ve cr
Te
4ξ
4ξ
3
Te
mi

If the drift velocity exceeds the critical value given in Equation (2.10), the ECDI will be
ﬁrst destabilized, otherwise the KCSI will be alternatively destabilized. In Figure 3, weshow the critical drift velocity normalized by the electron thermal velocity as a functionof the mass ratio for the case of Ti /Te = 8.0. The ECDI will be ﬁrst destabilized in theregion represented by red while the KCSI has the larger growth rate in the region shown inblue. The critical drift velocity decreases as the mass ratio increases until approximately
mi /me = 500. This is because the peak growth rate of the ECDI is almost independentof the mass ratio [Forslund et al., 1972] while that of the KCSI is the substantial fractionof the lower hybrid frequency whose value decreases with increase of the ion mass. Whenthe mass ratio is larger than approximately 500, the critical velocity is almost constantbecause of increase in ξ.

9Chapter 3
Anomalous Electron Heating
In the previous chapter we have shown that the ECDI dominates the KCSI as theelectron-ion relative velocity exceeds a critical value which varies according to the mass
ratio between an electron and ion, but seems to be a substantial fraction of the electron
thermal velocity. In this section we study the electron heating mechanism due to theunstable waves associated with these instabilities, making use of a 2-1/2 dimensional
electromagnetic full particle code with the particle-in-cell method (see Appendix B for
detailes of this code algorithm). In this code, we assume a doubly periodic boundary inthe x and z directions. The numbers of grid points are Nx × Nz = 256 × 512, and thegrid separations are ∆x × ∆z = 0.1λe × 0.1λe . The number of particles per grid is 36, so
the total number is approximately 5 × 106 . The initial values of the plasma parametersare ωce /ωpe = 0.5, Ti /Te = 8.0, mi /me = 100, and βi = 1.0, which are the same as thosein the case of Figure 2. The conﬁguration of the uniform magnetic ﬁeld B 0 and thedrift velocity V d is also given in the electron frame as assumed in the previous chapter.
We conducted simulations mainly for two cases of the drift velocity; one is Vd /ve = 0.8
where the KCSI is expected to dominate and the other is Vd /ve = 2.0 where the ECDI
is expected to dominate.
In Figure 4, we show the time development of the electric ﬁeld energy for the cases of
the two diﬀerent drift velocities. The total energy and the x, y and z component energies
are represented by red, green, yellow, and blue curves, respectively. Figure 4 indicatesthat the x component of the electric ﬁeld is overwhelming the other components in bothcases while the growth time, in which the total energy of the electric ﬁeld develops tothe saturation level, is much shorter in the case of (b) than (a).
In Figure 5, we present the typical wave spectra in the period during which theamplitude of the unstable waves is growing linearly for the cases of (a) Vd /ve = 0.8 and10(b) Vd /ve = 2.0. The upper panels in both cases are the spectra in the k⊥ − k space, andthe lower panels are those in the ω − k⊥ space. Because the x component of the electricﬁeld develops dominantly, we use only Ex to evaluate these spectra. Figure 5a indicatesthat the unstable waves, whose frequencies are approximately the lower hybrid frequency(ωlh ), propagate obliquely to the ambient magnetic ﬁeld. This mode corresponds to the
KCSI as expected in the linear analysis. On the other hand, in Figure 5b the characteristic
spectra of the KCSI disappears and the ECDI, which propagates perpendicularly to theambient magnetic ﬁeld and whose frequencies are a little above the electron cyclotron
frequency (ωce ), becomes dominant. Notice that the upper harmonics appear near theupper hybrid frequency (ωuh ), nearly the twice of the electron cyclotron frequency inthe current parameters. This is because the ECDI is excited by the resonance between
the drifting ions and the electron Bernstein waves in the electron frame [Forslund et al.,1970].

3.1

Distribution Function of Electrons

The distribution functions of electrons are shown in Figure 6 at the time when theamplitude of the unstable waves has been saturated (see Figure 4) for the case of thetwo drift velocities. Figure 6a indicates that electrons are heated parallel to the ambient
magnetic ﬁeld by the unstable waves associated with the KCSI. This is because theobliquely propagating waves due to the KCSI make potential troughs along the magnetic
ﬁeld, and electrons, which are able to move freely in the ﬁeld-aligned direction, canbe easily trapped by the potential troughs [McBride et al., 1972]. On the other hand,
Figure 6b shows that electrons are heated perpendicularly to the magnetic ﬁeld dueto the ECDI. In this case, the unstable waves propagate almost perpendicularly to themagnetic ﬁeld and obliquely propagating waves are weak (see Figure 5b) so that electrons
are strong aﬀected by the perpendicularly propagating waves. Since electrons are strongly
magnetized, they are not trapped eﬀectively by the potential troughs but some electrons
can be trapped in a portion of their cyclotron motions and be accelerated by the waves
if the wave amplitudes are enough large. Because the accelerated electrons continue thegyromotion, electrons will be heated isotropically in the x and y directions.
From these results we suggest that the mechanism of an electron heating and theconsequent distribution functions depend on the value of the electron-ion relative velocity.

11

3.2

Heating Level

We next consider how the characteristic temperature of electrons varies according
to the electron-ion relative velocity. We deﬁne the characteristic temperature as Te =2
2
me ve
/2 for the case in which the KCSI is dominant, and Te⊥ = me ve⊥/2 for the case22
in which the ECDI is dominant, where ve and ve⊥ are the root mean squares of theparallel and perpendicular velocities of electrons, respectively. In Figure 7a we showsuch characteristic temperature at the time when the amplitude of the unstable waves
has been saturated, as a function of the electron-ion relative velocity. The solid line
denotes the case where the ECDI is dominant and electrons are perpendicularly heated,
and the dashed line represents the case where the KCSI is dominant and electrons areheated parallel to the magnetic ﬁeld. Figure 7b represents the growth time which is thenecessary time for the amplitude of the electric ﬁeld to have been saturated. We ﬁnd thatthe growth time for the case where the ECDI is dominant is very short, compared withthat for the case of the KCSI, and seems to be nearly independent of the drift velocity.
Figure 7a and 7b indicate that electrons can be rapidly heated by the ECDI up to higher
temperature as the drift velocity increases.
The electron heating by the ECDI is possibly caused by the trapping in the potential
troughs of the electrostatic waves propagating perpendicularly to the ambient magnetic
ﬁeld, as is mentioned below. The linear growth of the wave amplitude is expected to be
slowed when the electron thermal velocity approaches the phase velocity of the electricﬁeld Vφ and the number of resonant electrons increases. Because the ECDI is excitedby the Landau resonance between the drifting ions and the electron Bernstein waves inthe electron frame, we can approximate Vφ  Vd . Then the wave amplitudes eventuallysaturate when the electron thermal velocity reaches the drift velocity. After the trappedelectrons are released from the potential troughs, they continue the cyclotron motions sothat the distribution function of electrons will be almost isotropic in the vx − vy plane.
Therefore the electron perpendicular temperature at the time when the wave amplitudes
are saturated is expected to be Te,sat /Te,0 = (Vd /ve )2 . In other words, electrons can gain
energy from the waves by the trapping when Vd > ve . In Figure 8 we show such electron
temperature as a function of the drift velocity. The square are the simulation results
which are approximated by the dashed line. The solid line represents the above mentioned
theoretical prediction. We ﬁnd that the simulation results are in good agreement withthe theory though it is very simple.
We present here, for reference, the detailed heating mechanism suggested by Chen12

and Birdsall [1973], in which it is described how an electron is trapped by an electrostatic
potential trough and accelerated. An electron trapping by an electrostatic wave propagating perpendicularly to a magnetic ﬁeld is more complex than that in the case there
is no magnetic ﬁeld, because a strongly magnetized electron is likely to make an escapefrom a potential trough. Consider an electron whose unperturbed gyromotion guiding
center is at x = 0, and an electrostatic wave propagating perpendicularly (directing to+x) to a magnetic ﬁeld (directing to +z). This electron motion satisﬁes the following
equationdvxe= − Ex − ωce vy ,
dt
me
dvy
= ωce vx.
dt

(3.1)
(3.2)

From the combination of Equations (3.1) and (3.2), the total force on the electron in thex direction is



1
∂Φ
∂
2 2
−eφ(x, t) + me ωce x = − ,Fx = −
∂x
2
∂x

(3.3)

1
2 2
x.
Φ = −eφ(x, t) + me ωce
2

(3.4)

where

φ(x, t) = φ0 cos(kx − ωt) is a wave potential propagating to the x direction, and φ0 isassumed to be constant. In Equation (3.4) the ﬁrst term of the right hand side is fromthe wave potential, on the other hand, the second term represents the magnetic eﬀect
and by itself simply describes the gyromotion. A schematic view of wave and particle
force potentials in the electron frame are shown in Figure 9 cited from Chen and Birdsall[1973]. The black dashed curve is the combined potential Φ. Away from x = 0, therelative eﬀect of the wave potential becomes weaker due to the increasing v × B 0 force
2
and, eventually, the combined potential troughs vanish at |x| = xv  ekφ0/me ωce, wherethe wave electric force equals the magnetic force.
Since we consider here the case of Vd > ve , where the ECDI is dominant, the relation
between k and re is kre < 1, and from the linear theory we found that the fastest growing
wave frequency was approximately the electron cyclotron frequency. These facts indicatethat an electron sees at most two troughs during a gyromotion. If the magnetic ﬁeld isabsent, the bounce period is τb = 2π(me /ek 2 φ0 )1/2 . There are two conditions for electrons
to be eﬀectively trapped by the wave potential troughs in the presence of the magnetic
13

ﬁeld. The ﬁrst is that the time during which an electron sees a single trough, τ1 ∼ λ/2Vd ,must be longer than τb , where λ is the wave length. The second is that the time during
which an electron can stay in a traveling trough, τ2 ∼ xv /Vd , has to be longer than τb . If
the wave amplitude is large enough as the above conditions are satisﬁed, electrons will
be eﬀectively trapped by wave potential troughs and the model described in Figure 9 isjustiﬁed.
The model in Figure 9 represents the following scenario. A trapped electron at the
point A (x  xv ), which has the total energy (the kinetic energy plus the potential energy)ΦA , is trapped by the trough and carried to the point B (x  xv ) by the traveling wave.
At this point, the potential trough vanishes and the electron is released to continue thegyromotion with the total energy ΦB  me (ωce xv )2/2.
If Te,sat /Te,0 = (Vd /ve )2 is assumed, we can set ωce xv  Vd , which gives
E0,sat = Vd B0 ,(3.5)

where E0,sat is the saturation level of the wave amplitude of the fastest growing mode.
The comparison between (3.5) and the simulation results is shown in Figure 10. The
theoretical prediction in (3.5) well agrees with the simulation results, which justiﬁes themodel of Figure 9.

3.3

The Case of Massive Ion with mi /me = 400

So far in this chapter, we have considered the case of mi /me = 100, where the electron
heating occurred perpendicularly to the magnetic ﬁeld due to the ECDI when Vd /ve > 1
and parallel to the magnetic ﬁeld due to the KCSI when Vd /ve < 1. We have found thatelectrons are rapidly heated by the ECDI until the electron thermal velocity reaches theion bulk velocity in the electron frame. In the case of mi/me = 100, the electron heating isdue to the fastest growing mode. However the critical drift velocity is below the electron
thermal velocity in the case of mi /me = 1836, as can be seen in Figures 1(a) and 3. Thus
it is possible for the ECDI to be excited even when the drift velocity is under the electron
thermal velocity. If electrons are heated by the ECDI in this case, the mechanism candiﬀer from that in the case of mi /me = 100, because the upper limit of the electron
thermal velocity is the drift velocity in the latter case. Therefore we need to examine
how electrons are heated by the ECDI when Vd < ve . Because it is diﬃcult to conduct
the full particle simulations for the case of mi /me = 1836, with enough resolution, due to14

the restriction of the computer resources, we execute the simulations with mi /me = 400
in this section.
The time evolution of the wave spectra of Ex in the k⊥ − k space is shown in Figure
11, for the case of Vd /ve = 0.95. The ECDI, which has the peak growth rate at (k⊥ , k ) =(±6.2, 0.0) in this condition, is ﬁrst excited and rapidly grows during tωpe  400, as is
expected from the linear theory. However, it is soon dumped and instead a low wave
number (viz. a low frequency) mode, the KCSI, becomes dominant after tωpe  480. In
Figure 12, we present the electron distribution functions at the corresponding time to
the wave spectra in Figure 11. It is interesting that electrons are not heated eﬀectively,
although the ECDI is dominant until tωpe  400. This is because a lot of electrons canbe resonant with the electrostatic waves in the case of Vd /ve < 1, as a result, the unstablewaves are attenuated before they have grown large enough to eﬀectively trap electrons.
Indeed a low saturation level is also deduced from Figure 10 in the case of Vd /ve < 1.
Since the KCSI is dominant after tωpe  480, electrons are eventually heated parallel tothe ambient magnetic ﬁeld.

15

Chapter 4
Discussion
4.1

Travel Time and Growth Time

So far we have assumed that electrons are strongly magnetized while ions are unmagnetized. However this assumption is valid only inside the diﬀusion region, which spreads
from the X-type neutral line with the order of the ion inertial length λi in the z direction
and 10λi in the x direction [e.g., Shay et al., 1998]. Since ions outside the diﬀusion region are also expected to be magnetized, the relative velocity between electrons and ionsnearly disappears, stabilizing both the ECDI and the KCSI. It is therefore importantto check whether there is enough time for the unstable waves to be saturated and heatelectrons up to the level described in Figure 7a while electrons are passing through thediﬀusion region. In Figure 13 we represent the growth time (which is the simulation
result) of the most unstable wave associated with the ECDI (solid line) and the KCSI
(dashed line), which are the same as those in Figure 7b, and the travel time of electrons
drifting across the diﬀusion region (thick solid curve), as a function of the drift velocity. We estimated the travel time as roughly 10λi /Vd , using the same parameters withthose in the simulation. The travel time may be underestimated because the relativeelectron-ion velocity Vd becomes lower during the passage. Figure 13 indicates that thetravel time is greater than the growth time for all points, so that there is enough time
for the two unstable modes to be saturated and heat the electrons up to the level shownin Figure 7a.

16

4.2

Comparison with Observations

The instabilities associated with a cross-ﬁeld relative motion between electrons andions vary according to the relative magnitude of the relative electron-ion velocity to theinitial thermal velocity of electrons. Here we compare the magnitude of the electron
drift velocity, assuming that ions have no bulk velocity, with that of the electron thermalvelocity inferred from observational results and examine which mode is expected to be
excited in the reconnection region formed in the Earth magnetotail. Since the upper
limit of the outﬂowing electron drift velocity just outside the electron diﬀusion region
is expected to be the electron Alfvén velocity VAe in the inﬂow region [e.g., Shay et al.,
2001; Hoshino et al., 2001] as already mentioned, we assume the relative electron-ion
velocity Vd is equivalent to VAe . Thus we compare the electron thermal velocity ve withthe electron Alfvén velocity VAe so as to determine which unstable mode and the resultant
electron heating is operative.
In order to calculate the electron Alfvén velocity we assume BL = (Bx2 + By2 )1/2 = 15
nT and nps ≈ 0.3 cm−3 after Baumjohann et al. [1989] and obtain VAe /c ≈ 0.085. On
the other hand, to estimate the electron thermal velocity in the plasma sheet at the timeof the reconnection, we use, for example, the observational data presented by Shinohara
et al. [1998] (see Table 1, therein), which were obtained by Geotail observation. Though
their data were probably taken outside the diﬀusion region because they observed fast
tailward ion ﬂows with southward Bz at the same time, we assume that the plasma
parameter necessary for calculating of the electron thermal velocity, that is, the electron
temperature is almost common around the diﬀusion region. From their data we obtainseveral values ranging from ve /c = 0.066 (Te = 1.1 keV) to ve /c = 0.119 (Te = 3.6 keV)(see Figure 14). Since the linear theory represents that the critical drift velocity, where
the growth rate of the ECDI exceeds that of the KCSI, is approximately 0.6ve (see Figure1a), the electron Alfvén velocity seems to exceed the critical value for any case so that
the ECDI is expected to be excited. However the electron heating associated with the
ECDI is weak for the case of Vd /ve < 1, as mentioned in Section 3.3. Thus the electron
heating by the KCSI will dominate in the case of Vd /ve < 1 even though the growth rateof the KCSI is lower than that of the ECDI. Therefore we can expect that both modes
are possible to be excited and heat electrons in parallel or perpendicular to the ambient
magnetic ﬁeld.
As for the identiﬁcation of the excited mode through an analysis of the wave dataobtained by spacecraft observations, the frequency of excited wave is quite helpful. The
17

real frequencies of the KCSI with the maximum growth rate for the various drift velocities can be read from Figure 1 to range from 0.04 ωpe to 0.11 ωpe . The wave numbers(k⊥ ) corresponding to these cases can be evaluated to range from 1.01 λ−1
to 0.93 λ−1ee .
Therefore the Doppler-shifted frequencies (ωr − k⊥ v) are in the range from −0.011 ωpe to−0.010 ωpe (see Figure 15), where we assume that the observation will be made in the ion
frame and take v = Vd . In the same way, the Doppler-shifted frequencies of the waves
excited by the ECDI take the values from −0.062 ωpe to −0.039 ωpe , where the values of−1
ωr and k⊥ range from 0.545 ωpe to 0.678 ωpe , and 12.14 λ−1e to 5.74 λe , respectively. Thus

for the observed waves which have the frequencies near the the lower hybrid frequency(ωlh = 0.012 ωpe in Figure 1), we can assume that they are excite by the KCSI. On theother hand, for the waves which have lower frequencies around 0.05 ωpe , it can be assumed
that they are excited by the ECDI.

18

Chapter 5
Summary and Conclusion
In this paper, we examined a possible mechanism of an electron heating due to astrong Hall current, which is caused by a large velocity diﬀerence between electrons andions inside the diﬀusion region. We ﬁrst solved the linear kinetic dispersion relation inorder to identify possible instabilities. We found that the unstable mode varied according
to the electron-ion relative velocity. In the case of the real mass ratio, mi/me = 1836,
the fastest growing mode is the high frequency mode associated with the ECDI when
Vd /ve  0.6, while it is the low frequency mode associated with the KCSI when Vd /ve 
0.6. However the critical value of the drift velocity, where the sort of the most unstablemode changes, increases as the mass ratio decreases. This is because the growth rateof the ECDI is almost independent of the mass ratio, diﬀerent from that of the KCSI,
which decreases as the mass ratio increases since the value is the substantial fraction of
the lower hybrid frequency. This is the character worth attention when we conduct a
numerical simulation, in which it is usually diﬃcult to use the real mass ratio because of
the restriction of the computer resources.
We next conducted the 2-1/2 dimensional electromagnetic full particle simulations to
examine the electron heating due to the unstable waves. We found that electrons areheated perpendicularly to the ambient magnetic ﬁeld due to the ECDI, while the heatingis in the parallel direction by the KCSI. In particular, we showed that electrons could bevery quickly heated up to ve  Vd by the trapping in the electrostatic potential troughs,
when the ECDI was dominant. We also showed that the electron heating associated with
the ECDI could be, however, ineﬀective when Vd /ve < 1 even though the ECDI is ﬁrst
excited. This is because the number of resonant electrons is so large when Vd /ve < 1that the wave amplitude is saturated at low level and soon dumped before it grows large
enough to trap electrons eﬀectively. It is important that the most dominant instability
19and the resultant mechanism of an electron heating vary around Vd /ve = 1.0, because
the electron thermal velocity may be comparable with the electron-ion relative velocity
inside the diﬀusion region according to the satellite observations. The identiﬁcation of
the dominant unstable mode will be possible by an analysis of the wave data obtainedby spacecraft observations, since the frequency ranges of the KCSI and the ECDI are
diﬀerent; the former is nearly the lower hybrid frequency while the latter is around|ω  | = 0.05ωpe , higher than the former. The Buneman-type instability is not necessarily
excited because this instability is a magnetohydrodynamic limit of the ECDI and seems
to be dominant when Vd /ve >> 1.
Because we have concentrated our attention only on a local region inside the diﬀusion
region with high spatial and time resolution, we cannot discuss the dynamics around thediﬀusion region, including the feedback eﬀects of the electron heating to the diﬀusion
region. In order to understand the self-consistent mechanism of the magnetic reconnection including Hall eﬀects near the X-type neutral line, we need to conduct a simulationthat covers more extensive region with high spatial and time resolution especially in thediﬀusion region.

20

Acknowledgments
The author would like to express grateful appreciation to Prof. S. Machida for his
illuminating guidance and discussion. The author also wishes to thank Prof. T. Araki,
Prof. T. Iyemori, Mr. T. Kamei, Dr. M. Takeda, Dr. A. Saito, and Dr. M. Nosé for theirhelpful comments and discussions. The author is grateful to Drs A. Yajima, N. Terada,
Y. Miyashita, and all colleagues of Department of Geophysics, Kyoto University for theirvaluable instruction and technical support. The author ﬁnally would like to thank Prof.
S. Fujita of Meteorological College for his encouragement. Without his guidance to thespace physics, this work would not have been taken up and, of course, accomplished.

21

Appendix A
Derivation of the Kinetic Plasma
Dispersion Relation
In this Appendix we derive a kinetic plasma dispersion relation, including obliquepropagation of unstable waves and the electron cyclotron resonance, when the unperturbed distribution functions of electrons and ions are deﬁned as Equations (2.1) and(2.2), respectively. We assume wave frequencies are much higher than the ion cyclotron
frequency so that ions are unmagnetized while electrons still remain frozen-in to theambient magnetic ﬁeld.

A.1

Perturbed Distribution Function δf

We follow the method described in Nicholson [1983] in order to derive the generalform of the perturbed distribution function. The fundamental equation is the Vlasov
equation,qs∂fs+ v · ∇x fs +(E + v × B) · ∇v fs = 0,
∂t
ms

(A.1)

where suﬃx s denotes the particle species (electron or ion). In order to linearize Equation(A.1) around the unperturbed distribution function, we put
fs (x, v, t) = fs0 (v) + δfs (x, v, t),E(x, t) = δE(x, t),B(x, t) = B 0 + δB(x, t).

22

Notice that we study in the electron frame so that the uniform electric ﬁeld is zero. Thenthe linearized Vlasov equation for species s isdδfsqs(δE + v × δB) · ∇v fs0 ,
=−
dt
ms

(A.2)

where d/dt = ∂t + v · ∇x + qs /ms (v × B 0 ) · ∇v is the total diﬀerential along theunperturbed orbits of particles in the phase space. Therefore in order to derive theunperturbed distribution function both sides of Equation (A.2) are integrated along theunperturbed orbits. We pursue the unperturbed orbit (X(t), V (t)) of a particle whose
motion is controled by Newton’s low of motion in the real space. The governing equationsof a particle motion on this orbit aredX(t)= V (t),
dt

(A.3)

qs
dV (t)
V (t) × B 0[X(t), t].
=
dt
ms

(A.4)

and

From Equations (A.3) and (A.4), we know that the unperturbed orbit can be given by
the following form, tdX(t ) 
dt ,
dt
t
 t
dV (t ) dt ,V (t ) = v −
dt
t


X(t ) = x −(A.5)

where the integral constants are determined as they satisfy the relations X(t) = x andV (t) = v. As a result the orbit (X(t ), V (t)) represents that of a particle which approach
the point (x, v) at time t. Here we deﬁne a particle position in the phase space (x, v  )as x = X(t ) and v  = V (t), which mean that the particle on the orbit (X(t), V (t)) is
located at the position (x , v ) at time t. (Notice that capital letters (X, V ) represent
an orbit curve as a function of time t in the seven-dimensional phase space while small
letters (x, v) denotes the particle position on the orbit at particular time.) Then werewrite Equation (A.1) as
qs
d

 
δf
(x
,
v
,
t
)
=
−
[δE(x, t) + v  × δB(x, t )] · ∇v  fs0 (v  ).
s
dt
ms

(A.6)

Integrating both sides of Equation (A.6) from t = −∞ to t = t, we can obtain
 t
qs

 
dt [δE(x, t) + v  × δB(x, t )] · ∇v  fs0 (v  ).
δfs (x, v, t) = δfs (x , v , t = −∞) −
ms −∞
(A.7)
23

Equation (A.7) is a general form of the perturbed distribution function. In the following,
we calculate the wave form solution of Equation (A.7), assuming
δE(x, t) = δ Ê exp[i(k · x − ωt)],δB(x, t) = δ B̂ exp[i(k · x − ωt)],
δf s (x, v, t), = δ f̂ s (v) exp[i(k · x − ωt)],(A.8)

where δ Ê and δ B̂ are constant vectors, k = k⊥ êx + k êz , and ω = ωr + iγ. Since weare here only interested in unstable waves with positive growth rates, we set the initial
ﬂuctuation δfs (t = −∞) to be zero, which allow us to neglect the ﬁrst term in the righthand side of Equation (A.7). Thus Equation (A.7) is manipulated into the followingform,qs
δ fˆs = −
ms

 0
−∞

dτ [δ Ê + v  × δ B̂] · ∇v  fs0 (v  ) exp[i(k · ξ − ωτ )],(A.9)

where ξ = x − x and τ = t − t, and the position (x, v  ) is put on the orbit (X(τ ), V (τ ))which passes the point (x, v) at τ = 0.

A.1.1

Ion Contribution

Since we assume that ions are unmagnetized, the unperturbed orbits of ions arecalculated from Equations (A.3) and (A.4) with B 0 = 0 as
V (t ) = v,

X(t) = x + vτ.

Since v  is independent of τ , it is easy to execute the integration in Equation (A.9) and
we obtain
δ fˆi =δ Êe· ∇v fi0 ,mi i(ω − k · v)

(A.10)

where we put δ B̂ = 0 because of the assumption of unmagnetization. Making use of theisotropic nature of the unperturbed distribution function in the velocity space, we canmanipulate Equation (A.10) into the following form,δ fˆi = i2e (vx − Vd )δ Êx + vy δ Êy + vz δ Êz ∂fi0
,
2
mi
k·v−ω
∂v⊥

2
= (vx − Vd )2 + vy2.
where v⊥

24

(A.11)

A.1.2

Electron Contribution

Electrons conduct the cyclotron motion around the magnetic ﬁeld line. From Equations (A.3) and (A.4) we derive the unperturbed orbits of electrons as
Vx (τ ) = v⊥ cos(ωce (t − t) + θ) = v⊥ cos(ωce τ + θ),
Vy (τ ) = v⊥ sin(ωce (t − t) + θ) = v⊥ sin(ωce τ + θ),
Vz (τ ) = v ,v⊥[sin(ωce τ + θ) − sin θ],ωcev⊥[cos(ωce τ + θ) − cos θ],Y (τ ) = y −ωceZ(τ ) = z + vτ,X(τ ) = x +(A.12)

where v = vz , and θ (0 ≤ θ ≤ 2π) is the initial phase at τ = 0 in the velocity space. In
this case, since v  is a function of τ the execution of the integration in Equation (A.9) isslightly complicated. We eliminate δ B̂ from Equation (A.9) by using Faraday’s low andsubstitute (A.12) into the resultant equation. After some algebraic calculation we obtain∞
2e  nλ v⊥ Jn δ Êx + iv⊥ Jn δ Êy + vJn δ Êz ∂fe0πˆδ fe = iexp[n(θ − ) + λ cos θ],2
me n=−∞
ω − nωce − k v
∂v⊥
2(A.13)
where λ = k⊥ v⊥ /ωce , Jn = Jn (λ) is the Bessel function of order n, and Jn = dJn (λ)/dλ.
In deriving Equation (A.13) we used the following identical equation,eiλ sin θ=∞Jn (λ) exp[inθ],(A.14)n=−∞and the diﬀerential form of (A.14) about λ and θ,eiλ sin θ

sin θ = −i∞Jn einθ ,

n=−∞

A.2

e

iλ sin θ

cos θ

∞

n

λ
n=−∞

Jn einθ .

Dispersion Relation

In order to complete the linearize equation system which gives the linear dispersion
relation, we here introduce Maxwell’s equations. We assume the wave form solutions
given in (A.8) and replace ∂/∂t with −iω and ∇ with ik. Then Faraday’s law and25

Ampère’s law are expressed likek × δ Ê,ω
ω
ik × δ B̂ = µ0 δ ĵ − i 2 δ Ê,cδ B̂ =(A.15)(A.16)

where δ ĵ is the Fourier component of the ﬂuctuating current density δj. δ ĵ is deﬁned by



3
3
δ ĵ = evδ f̂i d v − vδ fˆe d v .
(A.17)
Substituting Equations (A.11) and (A.13) into (A.17) and executing the integration overthe velocity space, we can derive δj and the dielectric tensor by combining it with
Equations (A.15) and (A.16). On calculation of the integration in (A.17) we use thefollowing relations,
 ∞
1
2
2
xJn2(αx) e−x dx = In (α2 /2) e−α /2 ,
2
0 ∞
α
2
2
x2Jn (αx)Jn (αx) e−x dx = [In (α2 /2) − In (α2 /2)] e−α /2 ,
4

 2

0 ∞
n
α2
α2  2
2
3 
2 −x2
2
In (α /2) − In (α /2) e−α /2 ,x [Jn (αx)] e dx =
+
2
2α
4
4
0
where In (x) is the modiﬁed Bessel function of order n, and In (x) = dIn (x)/dx. Aftersome calculation we obtainδ ĵ = σ · δ Ê,(A.18)

where σ is the conductivity tensor and the elements are expressed by the following form,
∞
2
2
2εωpi
ε0ωpe
e−µ  22n In ξ0 Z(ξn ),σxx = −i 2 2 ω[1 + ξi Z(ξi )] sin θ − i
k vi
ω µ n=−∞
∞  2
2

2εωpe
n
−µ
µe
In + (In − In ) ξ0 Z(ξn ),σyy = −i
2
ω
2µ
n=−∞
σzz = −i∞222ε0ωpi2ε0ωpi
2
−µ
e
ω[1
+
ξ
Z(ξ
)]
cos
θ
−
i
In ξ0 ξn [1 + ξn Z(ξn )],ii
k 2 vi2
ω
n=−∞

∞
2

ε0ωpe
−µ
e
n(In − In )ξ0 Z(ξn ),σxy = −σyx =
ω
n=−∞
2
2
2ε0ωpiε0ωpeσxz = σzx = −i 2 2 ω[1 + ξi Z(ξi )] sin θ cos θ − i
k vi
ω



∞
2 −µ 
enIn ξ0 [1 + ξn Z(ξn )],
µ
n=−∞

∞
2 

ε0 ωpe−µ2µe(In − In )ξ0 [1 + ξn Z(ξn )],σyz = −σzy = −
ω
n=−∞

26

2 2
2
where ξn = (ω − nωce )/k ve , µ = k⊥ve /2ωce, ξi = (ω − k⊥ Vd )/kvi , θ = tan−1 (k⊥ /k ), andZ(ξ) is the plasma dispersion function [Fried and Conte, 1961]. We substitute (A.18)
into Equation (A.16) and use Equation (A.15) to eliminate δ B̂. Then (A.16) is reduced
into



δ Êx
Dxx Dxy Dxz
 Dyx Dyy Dyz   δ Êy  = 0,
(A.19)
Dzx Dzy Dzz
δ Êz
where
∞
2
2 −µ 
2ωpi
ωpee
c2 k 2
2
2
cos
θ
+
[1
+
ξ
Z(ξ
)]
sin
θ
+
n2 In ξ0 Z(ξn ),
i
i
2
2
2
2
ω
k vi
ω µ n=−∞
∞2
n2
c2k 2 2ωpe
I + (In − In ) ξ0 Z(ξn ),Dyy = 1 − 2 + 2 µe−µ
2 n
ω
ω
2µ
n=−∞

Dxx = 1 −

∞
2
2

2ωpi
2ωpe
c2 k 222−µIn ξ0 ξn [1 + ξn Z(ξn )],Dzz = 1 − 2 sin θ + 2 2 [1 + ξi Z(ξi )] cos θ + 2 e
ω
k vi
ω
n=−∞
∞
2

ωpe
−µ
n(In − In )ξ0 Z(ξn ),Dxy = −Dyx = −i 2 e
ω
n=−∞
2
2
2ωpi
ωpe
c2 k 2
Dxz = Dzx = 2 sin θ cos θ + 2 2 [1 + ξi Z(ξi )] sin θ cos θ + 2
ω
k vi
ω



∞
2 −µ 
nIn ξ0 [1 + ξn Z(ξn )],
e
µ
n=−∞

∞
2 

ωpe
−µ
(In − In )ξ0 [1 + ξn Z(ξn )].
Dyz = −Dzy = i 2 2µe
ωn=−∞(A.20)The kinetic dispersion relation is deﬁned by det[D] = 0.
In order to facilitate the discussion about the electrostatic and electromagnetic eﬀects
in the dispersion relation, we next derive the dispersion tensor D in a representation basedon the scalar and vector potentials δ φ̂, δ Â, rather than the electric ﬁeld δ Ê, after Tsaiet al., [1984]. Since δ Ê = −ikδ φ̂ + iωδ Â and k · δ Â = 0 (Coulomb gauge), we can write





k
k
δ Êx
δ φ̂
− kk⊥ 0 − k⊥ω δ Êy  = iω 
(A.21)
10   δ Ây  .
 0
k
δ Êz
δ Âz
1
−k 0
Using the potential representations, we can rewrite Equation (A.19) ask/ωδφ̂
D11 D12 D13
 D21 D22 D23   δ Ây  = 0.
D31 D32 D33
δ Âz
27

(A.22)

From Equations (A.19), (A.21) and (A.22), we can show that22k⊥ kkk⊥Dxx +Dxz + 2 Dzz ,
2
2
k
k
k
D22 = Dyy ,2k2
2k⊥ kk⊥D+Dzz ,D33 = 2 Dxx −
xz
k
k2
k2
k
k⊥
D12 = −D21 = − Dxy + Dyz ,kk22k⊥ − k2kkkDxz − Dzz ,D13 = 2 D31 = Dxx −
k⊥
k
kk⊥
k
2
k
k
Dxy + Dyz .
D23 = − 2 D32 =
k⊥
k⊥

D11 =

(A.23)The representation of (2.5) can be obtained by substituting (A.20) into (A.23). In the δ φ̂,δ Â representation, the dispersion relation reduces to D11 = 0 if electromagnetic eﬀect isignored.

28

Appendix B
Algorithm in the Full Particle CodeIn a usual plasma particle simulation, we employ modeled particle called superparticles, which has a ﬁnite spatial scale about the order of the spatial grid size. This isbecause we need to simulate with far fewer particles than in a real case. For example,
provided that plasma includes electrons with number density 1.0 /cc and thermal velocityc (the light velocity), the Debye length is λD ≈ 105 m. If we set the grid separations to
be λD /10, the number of particles par grid is approximately 1012 for a two-dimensional
simulation. Overall the simulation space with 100 × 100 grids, it is above 1016 that isnecessary number of particles to conduct a simulation. It is totally impossible to treatsuch a vast number of particles. Furthermore, the superparticles with broad shapes are
beneﬁcial to moderate numerical noises. In this appendix, we show the detailed algorithm
in our 2-1/2 dimensional full particle simulation code with the superparticles [Birdsall
and Langdon, 1995].

B.1

Fundamental Equations

The position xsj and the velocity v sj of the j-th particle obey Newton’s law,ms

dvsj
= qs [E(xsj ) + v sj × B(xsj )],dt

dxsj= v sj ,
dt

(B.1)
(B.2)

29

where ms and qs represent the mass and the charge of a particle (a superparticle) of
species s, respectively. The electromagnetic ﬁeld is described by Maxwell’s equations,∂B= −∇ × E,∂t
1 ∂E= ∇ × B − µ0 j,
c2 ∂tρ∇·E = ,
ε0
∇ · B = 0,

(B.3)
(B.4)
(B.5)
(B.6)

where the charge density ρ and the current density j is deﬁned at each grid X i ,ρ(X i ) =

qs ns (X i ),(B.7)

qs ns (X i )v s ,

(B.8)

s

j(X i ) =swhere ns is the number density.
In order to normalize Equations from (B.1) to (B.8), we put the normalized quantities
in the following forms,
ms
ωpe
qs
, x∗s =
xs , t∗ = ωpe t, qs∗ =
me
c
e
vs
ω
eE
, E∗ =,v ∗s = , ω ∗ =cωpeme cωpe


2
2
c
c
eB
1
1
∗
∗
∗
, ρ =
ρ, j =j,B =me ωpee ωpe
ec ωpe

2
c
∗
ns ,
ns =
ωpe
m∗s =

(B.9)

where the asterisk denotes that the quantities are dimensionless. Using (B.9), We canrewrite the fundamental equations as the following forms.
Newton’s law,m∗s

dv ∗sj= qs∗ [E ∗ (x∗sj ) + v ∗sj × B ∗ (x∗sj )],
dt∗

dx∗sj
= v ∗sj ,dt∗(B.10)(B.11)

30

Maxwell’s equations,∂B ∗= −∇∗ × E ∗ ,
∗
∂t
∂E ∗
1
= ∇∗ × B ∗ − ∗ j ∗ ,
∗
∂t
n0
1
∇∗ · E ∗ = ∗ ρ∗ ,n0∗∗
∇ · B = 0,

(B.12)(B.13)(B.14)(B.15)whereρ∗ (x∗ ) =qs∗n∗s (x∗ ),(B.16)

qs∗ n∗s (x∗ )v ∗s .

(B.17)sj ∗ (x∗ ) =s

We ﬁrst advance the velocity v ∗sj and the position x∗sj in time, using Equations (B.10)
and (B.11). After all particle velocities and positions are determined, we can next obtainthe charge density and the current density at each grid point from Equations (B.16) and(B.17), respectively. We ﬁnally advance in time the electric ﬁeld and the magnetic ﬁeld
through Maxwell’s equations (B.12) – (B.15). It is noticeable that we separately derive
the longitudinal and transversal components of the electric ﬁeld from (B.14) and (B.13),
respectively, to satisfy the charge conservation.
In the following section, we omit the asterisks to avoid the complication.

B.2

Boundary and Initial Conditions

The boundary condition in our simulation is set to be doubly periodic. The initialconditions should be given with self-consistency to the particle position xsj and velocity
v sj , the ﬂuctuating ﬁelds (E, B 1 ), and the uniform ﬁeld B 0 . The initial distribution ofparticles are assumed to be Maxwellian obtained by making use of Box-Muller method[Press et al., 1992], by which particle velocities with the thermal velocity vt are describedin the following form,
v1 = vt −2 ln x1 cos 2πx2,v2 = vt −2 ln x1 sin 2πx2,(B.18)
where x1 and x2 are the uniform random numbers. The velocities v1 and v2 are independent of each other.
31

B.3

Momentum Equations

The ﬁnite-diﬀerence equations of (B.10) and (B.11) are


n+1/2
n−1/2
n+1/2
n−1/2
qs
v
− vs
+
v
vs
s
s
=
× B n (xns ) ,E n (xns ) +
∆t
ms
2
− xns
xn+1
s
= v n+1/2,s∆t(B.19)(B.20)

where (n−1/2, n, n+1/2) represent time grids at ((n−1/2)∆t, n∆t, (n+1/2)∆t). E(xs )and B(xs ) are the values at the particle position xs . However E and B are deﬁned onlyat spatial grid points X i . Thus it is necessary to assign the values of electromagnetic
ﬁeld deﬁned at grid points X i to those at particle positions xj . In our code, we give aparticle with a broad shape in the following form,
S(x, y) =(∆x − |x|)(∆y − |y|),∆x∆y|x| ≤ ∆x,|y| ≤ ∆y,(B.21)

where S(x, y) is called the shape factor. Then the assignment to a particle position
xsj is conducted from the nearest four grid points (Xi , Yk ), (Xi+1 , Yk ), (Xi , Yk+1 ), and(Xi+1 , Yk+1 ),E(xj , yj ) = E(Xi , Yk )S(Xi − xj , Yk − yj ) + E(Xi , Yk+1 )S(Xi − xj , Yk+1 − yj )+E(Xi+1 , Yk )S(Xi+1 − xj , Yk − yj ) + E(Xi+1 , Yk+1 )S(Xi+1 − xj , Yk+1 − yj ).
This technique is called the particle-in-cell (PIC) method.

B.4

Charge Assignment

In order to obtain the charge density ρ and the current density j at each grid point
X i , it is necessary to assign particle positions xsj and velocities v sj to grid points. In ourcode, the quantities associated with a particle (xsj , v sj ) are assigned to the nearest fourgrid points by making use of PIC method as is conducted in the previous section. Since
the particle velocities are deﬁned at half-integer time, ρ(X i ) and j(X i ) are expressed asnqs S(xnsj − Xi , ysj− Yk ),
ρn (X i ) =s

j n+1/2 (X i ) =

j


s

n+1
n+1
n
n
− Yk )
n+1/2 S(xsj − Xi , ysj − Yk ) + S(xsj − Xi , ysj

qs v sj

2

j

32

.

B.5

Poisson’s Equation

After the spatial distribution of the charge density are derived, we can obtain the longitudinal electric ﬁeld E L from Poisson’s equation,1
ρ(X i ).
n0

∇ · E L |X i =(B.22)

We solve Equation (B.22) by executing Fourier transformation. Since E L can be expressed by a scalar potential φ in the form,
E L (X i ) = −∇φ|X i .

(B.23)1
ρ(X i ).
n0

(B.24)(B.22) is rewritten as∇2φ|X i = −Here we write the position of grid points with the integer l and m,
X i = l ∆xêx + m ∆yêy ,l = 0, 1, · · · , Nx − 1

Nx : the number of grids in the x direction,m = 0, 1, · · · , Ny − 1

Ny : the number of grids in the y direction,and the physical quantities at the grid points as Al,m.
The Fourier transformation in the discrete form is, for example,
ρkx ,ky = ∆x∆yl

ρl,m e−i(kx xl +ky ym ) ,(B.25)mwhere(xl, ym ) = (l∆x, m∆y ),
2πnx 2πny(kx , ky ) =,,Nx ∆x Ny ∆ynx,y = 0, ±1, ±2, · · · , ±

Nx,y
Nx,y
− 1, −.
2
2(B.26)On the other hand, the inverse Fourier transformation in the discrete form is1
ρl,m =(Nx ∆x)(Ny ∆y)

Nx /2−1Ny /2−1nx =−Nx /2 ny =−Ny /233

ρkx ,ky ei(kx xl +ky ym ) .

(B.27)

The ﬁnite-diﬀerence equation of (B.24) is
φl+1,m − 2φl,m + φl−1,m φl,m+1 − 2φl,m + φl,m−1
n0
+
=−
.
2
2
∆x
∆y
ρl,m
And its Fourier transformation isρ2
φk = k ,
Kk
n0
where


2

Kk = kx2 dif 2

kx ∆x
2




+ ky2dif 2

ky ∆y
2

(B.28)
,dif θ ≡sin θ.
θThen the ﬁnite-diﬀerence form of the x component of Equation (B.23) is
Ex,l,m = −φl+1,m − φi−1,m.
2∆x

And its Fourer transformation is
Ex,k = −iκx,k φk ,where κx,k = kx dif (kx ∆x).

(B.29)

From (B.28) and (B.29), we can obtain
Ex,k = −i

κx,k

ρ .
n0 K 2 kk

In the same way, we can get the y component,
κy,kEy,k = −iρ .
n0K 2 kk

(B.30)(B.31)The longitudinal electric ﬁeld E L is derived by the inverse Fourier transforming of (B.30)
and (B.31).

B.6

Faraday’s Law and Ampère’s Law

We ﬁnally advance the ﬂuctuating magnetic ﬁeld and the transversal electric ﬁeld,
using Faraday’s law and Ampère’s law. The ﬁnite-diﬀerence form of Equations (B.12)
and (B.13) is n+1
B n+1 − B n
E T + E nT
= −∇ ×
,
∆t
2
 n+1

B
1 n+1/2− E nT
+ Bn
E n+1
T
=∇×
,
− jT
∆t
2
n0
34

(B.32)(B.33)

where the suﬃx T denotes the transversal component of the quantity. Combining (B.32)
and (B.33), we eliminate B n+1 to obtain E n+1
from the previously given quantities (E nT ,T
n+1/2

B n, j T
),




∆t2
∆t n+1/2∆t2
n+1∇ × ∇× E T = 1 −∇ × ∇× E nT + ∆t∇ × B n −j. (B.34)
1+
4
4
n0 T
Here, making use of the notation ∇ × ∇ × E T = −∇2E T , we rewrite the above equation
as





∆t2 2
∆t2 2
∆t n+1/2
n+1
1−j.
∇ ET = 1 +∇ E nT + ∆t∇ × B n −
4
4
n0 T

(B.35)In Equation (B.35), the ﬁnite-diﬀerence forms of the terms including ∇ operator are∇2 E =

E l+1,m − 2E l,m + E l−1,m E l,m+1 − 2E l,m + E l,m−1
+
,
∆x2
∆y 2
 B

z,l,m+1 −Bz,l,m−1



2∆y



−Bz,l−1,m∇ × B =  − Bz,l+1,m2∆x.
Bx,l,m+1 −Bx,l,m−1
By,l+1,m −By,l−1,m
−
2∆x
2∆y
Thus the Fourier transformation of (B.35) is∆t n+1/2
∆t2 2
∆t2 2
n+1
Kk E
Kk E nk + i∆t κk × B nk −
= 1−
1+
j
,
k
4
4
n0 k
n+1/2

(B.36)
n+1/2

nn) are the Fourier components of (E n+1
where (E n+1 , E nk , B nk , j
T , ET , B , j Tkk2
spectively, and Kk
and κk are deﬁned as the following forms,




kx ∆x
ky ∆y
2
2
2
2
2
+ ky dif,
Kk = kx dif
2
2), re-

kx dif (kx ∆x)κk =  ky dif (ky ∆y)  .
0Then we can calculate B n+1 , the Fourier component of B n+1 , by Fourier transformingkEquation (B.32).
 n+1

E k + E nk
n+1n
= B k − i∆t κk ×
B.
(B.37)k2
n+1) are obtained by inversely Fourier
The advanced electromagnetic ﬁelds (E n+1
T , B
n+1
n+1
transforming (E , B ).
kk35

References
Baumjohann, W., G. Paschmann, and C. A. Cattell, Average plasma properties in thecentral plasma sheet, J. Geophys. Res., 94, 6597-6606, 1989.
Baumjohann, W., G. Paschmann, T. Nagai, and H. Lühr, Superposed epoch analysis of
the substorm plasma sheet, J. Geophys. Res., 96, 11,605-11,608, 1991.
Birdsall, C. K., and A. B. Langdon, Plasma physics via computer simulation, Instituteof Physics, London, 1995.
Coroniti, F. V., and A. Eviatar, Magnetic ﬁeld reconnection in a collisionless plasma,
Astrophys. J. Suppl. Ser., 33, 189-210, 1977.
Davidson, R. C., N. T. Gladd, C. S. Wu, and J. D. Huba, Eﬀects of ﬁnite plasma betaon the lower-hybrid-drift instability, Phys. Fluids, 20, 301-310, 1977.
Forslund, D. W., R. L. Morse, and C. W. Nielson, Electron cyclotron drift instability,
Phys. Rev. Lett., 25, 1266-1270, 1970.
Forslund, D., R. Morse, and C. Nielson, and J. Wu, Electron cyclotron drift instability
and turbulence, Phys. Fluids, 15, 1303-1318, 1972.
Fried, B. D. and S. D. Conte, The plasma dispersion function, Academic Press, New Yorkand London, 1961.
Gary S. P. and J. J. Sanderson, Longitudinal waves in a perpendicular collisionless plasma
shock I. Cold ions, J. Plasma Phys., 4, 739-751, 1970.
Gary S. P., Longitudinal waves in a perpendicular collisionless plasma shock II. Vlasovions, J. Plasma Phys., 4, 753-760, 1970.
Hesse, M. and D. Winske, Hybrid simulations of collisionless reconnection in current
sheets, J. Geophys. Res., 99, 11,177-11,192, 1994.
Hesse, M. and D. Winske, Electron dissipation in collisionless magnetic reconnection, J.
Geophys. Res., 103, 26,479-26,486, 1998.
36

Hesse, M., K. Schindler, J. Birn, and M. Kuznetsova, The diﬀusion region in collisionlessmagnetic reconnection, Phys. Plasmas, 6, 1781-1795, 1999.
Hoshino, M., T. Mukai, T. Yamamoto, and S. Kokubun, Ion dynamics in magnetic reconnection: Comparison between numerical simulation and Geotail observations,J. Geophys. Res., 103, 4509-4530, 1998.
Hoshino, M., K. Hiraide, and T. Mukai, Strong electron heating and non-Maxwellian
behavior in magnetic reconnection, Earth Planets Space, 53, 627-634, 2001a.
Hoshino, M., T. Mukai, T. Terasawa, and I. Shinohara, Suprathermal electron acceleration in magnetic reconnection, J. Geophys. Res., 106, 25,979-25,997, 2001b.
Huba, J. D., N. T. Gladd, and K. Papadopoulos, Lower-hybrid-drift wave turbulence inthe distant magnetotail, J. Geophys. Res., 83, 5217-5226, 1978.
Kall, N. A. and P. C. Liewer, Low frequency instabilities in magnetic pulses, Phys. Rev.
A, 4, 2094-2103, 1971.
Lemons, D. S. and S. P. Gary, Electromagnetic eﬀects on the modiﬁed two-stream instability, J. Geophys. Res., 82, 2337-2342, 1977.
Lin, Y. and D. W. Swift, A two-dimensional hybrid simulation of the magnetotail reconnection layer, J. Geophys. Res., 101, 19,859-19,870, 1996.
Mandt, M. E., R. E. Denton, and J. F. Drake, Transition to whistler mediated magnetic
reconnection, Geophys. Res. Lett., 21, 73-76, 1994.
McBride, J. B., E. Ott, J. P. Boris, and J. H. Orens, Theory and simulation of turbulent
heating by the modiﬁed two-stream instability, Phys. Fluids, 15, 2367-2383, 1972.
Nagai, T., M. Fujimoto, Y. Sato, S. Machida, T. Terasawa, R. Nakamura, T. Yamamoto,T. Mukai, A. Nishida, and S. Kokubun, Structure and dynamics of magnetic
reconnection for substorm onsets with Geotail observations, J. Geophys. Res.,
103, 4419-4440, 1998.
Nagai, T., I. Shinohara, M. Fujimoto, M. Hoshino, Y. Saito, S. Machida, and T. Mukai,
Geotail observations of the Hall current system: Evidence of magnetic reconnection in the magnetotail, J. Geophys. Res., 106, 25,929-25,949, 2001.
Nakabayashi, J. and S. Machida, Electromagnetic hybrid-code simulation of magnetic
reconnection: Velocity distribution functions of accelarated ions, Geophys. Res.
Lett., 24, 1339-1342, 1997.
37

Nakamura, M. S., M. Fujimoto, and K. Maezawa, Ion dynamics and resaltant velocity
space distribution functions in the cource of magnetic reconnection, J. Geophys.
Res., 103, 4531-4546, 1998.
Nicholson, D. R., Introduction to plasma theory, John Wiley & Sons, New York, 1983.
Parker, E. N., The solar-ﬂare phenomenon and the theory of reconnection and annihilation of magnetic ﬁelds, Astrophys. J., Suppl. Ser., 8, 177-212, 1963.
Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical recipes
in FORTRAN, 2nd ed., Cambridge University Press, New York, 1992.
Sergeev, V. A., V. Angelopoulos, D. G. Mitchell, and C. T. Russell, In situ observationsof magnetotail reconnection prior to the onset of a small substorm, J. Geophys.
Res., 101, 19,121-19,133, 1995.
Shay, M. A., J. F. Drake, R. E. Denton, and D. Biskamp, Structure of the dissipation
region during collisionless magnetic reconnection, J. Geophys. Res., 103, 91659176, 1998.
Shay, M. A., J. F. Drake, and B. N. Rogers, Alfvénic collisionless magnetic reconnection
and the Hall term, J. Geophys. Res., 106, 3759-3772, 2001.
Shinohara, I., T. Nagai, M. Fujimoto, T. Terasawa, T. Mukai, K. Tsuruda, and T. Yamamoto, Low-frequency electromagnetic turburence observed near the substorm
onset site, J. Geophy. Res., 103, 20,365-20,388, 1998.
Sonnerup, B. U. Ö., 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, 1979.
Sweet, P. A., The neutral point theory of solar ﬂares, in Electromagnetic Phenomena in
Cosmical Physics, edited by B. Lehnert, pp. 123-134, Cambridge University Press,
London, 1958.
Tanaka, M., Macro-particle simulations of collisionless magnetic reconnection, Phys.
Plasmas, 2, 2920-2930, 1995.
Tsai, S. T., M. Tanaka, J. D. Gaﬀey, Jr., E. H. Da Jornada, C. S. Wu, and L. F.
Ziebell, Eﬀect of electron thermal anisotropy on the kinetic cross-ﬁeld streaming
instability, J. Plasma Phys., 32, 159-178, 1984.

38

Winske, D., M. Tanaka, C. S. Wu, and K. B. Quest, Plasmaheating at collisionless shocks
due to the kinetic cross-ﬁeld streaming instability, J. Geophys. Res., 90, 123-136,
1985.
Wong, H. V., Electrostatic electron-ion streaming instability, Phys. Fluids, 13, 757-760,
1970.
Wu, C. S., Y. M. Zhou, S. C. Guo, D. Winske, and K. Papadopoulos, A kinetic cross-ﬁeld
streaming instability, Phys. Fluids, 26, 1259-1267, 1983.

39

List of Figures
1(a) Peak growth rates γ/ωpe , and (b) the corresponding real frequencies
as a function of propagation angle θ for ECDI (solid curves) and KCSI
(dashed curves), with various values of the drift velocity, Vd /ve . The
plasma parameters are ωce /ωpe = 0.5, Ti /Te = 8.0, mi /me = 1836, andβi = 1.0.

2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43(a) Peak growth rates, and (b) the corresponding real frequencies as afunction of propagation angle θ for ECDI (solid curves) and KCSI (dashedcurves), with various values of the drift velocity. The values of the plasma
parameters are the same as those in the case of Figure 1 except the mass

3

ratio, mi /me = 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Mass ratio dependence of the critical drift velocity (Vd /ve )cr (solid line).

43

The plasma parameters are the same as those in the case of Figure 1 except
the mass ratio mi /me . If the drift velocity exceeds a critical value the
ECDI will be ﬁrst excited (red region), otherwise the KCSI will be ﬁrst
excited (blue region). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
444Time development of the electric ﬁeld energy for the case (a) Vd /ve = 0.8
and (b) Vd /ve = 2.0. Red curves represent the total energy, green, yellow,
and blue curves denote the x, y, and z component energies of the electric
ﬁeld, respectively. The values of the initial parameters are ωce /ωpe = 0.5,

5

Ti /Te = 8.0, mi /me = 100, and βi = 1.0. . . . . . . . . . . . . . . . . . .
Typical wave spectra of Ex while the amplitude of the unstable waves is

45

growing linearly; k⊥ - k space (upper panels) and ω - k⊥ space (lower
panels), for the case of (a) Vd /ve = 0.8 < (Vd /ve )cr and (b) Vd /ve = 2.0 >(Vd /ve )cr , where (Vd /ve )cr is the critical drift velocity normalized by theelectron thermal velocity. The values of the initial parameters are thesame as those in the case of Figure 4. . . . . . . . . . . . . . . . . . . . .

40

46

6

The distribution function of electrons at the time when the amplitude of
the unstable waves has been saturated, for the case of (a) Vd /ve = 0.8 and(b) Vd /ve = 2.0. The values of the initial parameters are the same as thosein the case of Figure 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

747(a) Characteristic temperature of electrons at the time when the amplitudeof the unstable waves has been saturated Te,sat, normalized by the initial
temperature Te,0, is plotted versus the drift velocity Vd normalized by the
initial thermal velocity of electrons ve . The solid line and the squares
represent the heating by ECDI, and the dashed line and the triangles
show the heating by KCSI. The characteristic temperature of electrons
2is Te⊥ = me v⊥/2 for the case of ECDI and Te = me v2 /2 for the caseof KCSI. (b) The growth time is also plotted versus the drift velocity.
The deﬁnitions of both lines and points are the same as (a). The initial

8

parameters are the same as those in the case of Figure 4. . . . . . . . .
Electron heating level at the time when the amplitude of the unstable47

waves has been saturated, for the case where the ECDI is dominant. The
squares are the simulation results which are approximated by the dashedline. The solid line represents the theoretical prediction, Te,sat /Te,0 =(Vd /ve )2 .
9

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

48

Schematic illustration of the electron trapping by an electrostatic wave
propagating perpendicularly to a magnetic ﬁeld, cited from Chen andBirdsall [1973]. This is shown in the electron frame. The electron gyrocenter is assumed to be at x = 0. A trapped electron at the point A,
which has total energy (potential energy plus kinetic energy) ΦA , can be
carried to the point B by propagating wave. At this point the potential
trough vanishes and the electron is released with the total energy ΦB . . .

10Saturation level of wave amplitude. The squares represent the simulation
results. The solid line is the theoretical prediction in which we employ the
relation Te,sat /Te,0 = (Vd /ve )2 and the model described in Figure 9. . . .

11

49

50

Time evolution of the wave spectra of Ex in the k⊥ − k space. The
initial values of the plasma parameters are ωce /ωpe = 0.5, Ti/Te = 8.0,
mi /me = 400, and βi = 1.0. The initial value of the ion drift velocity is
set to be Vd /ve = 0.95. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

51

12

Time evolution of the electron distribution function. The initial values ofthe plasma parameters and the drift velocity are the same as those in the13case of Figure 11. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The growth time of the most unstable wave associated with ECDI (solid

52line and squares) and KCSI (dashed line and triangles), which are thesame as those in Figure 7b, and the travel time of electrons drifting acrossthe diﬀusion region (thick solid curve) are plotted versus the drift velocity
normalized by the initial thermal velocity of electrons.
14. . . . . . . . . .

53

Four electron thermal velocities observed near the X-type neutral line by
Geotail spacecraft associated with substorms. These data set is takenfrom Shinohara et al. [1998]. The electron heating is perpendicular to themagnetic ﬁeld when ve < Vd , while it is parallel when ve > Vd . . . . . . .

15

54



Dopples-shifted angular frequencies (ω = ωr −k⊥ v) derived from the linear
theory for the various values of the drift velocity, assuming that we arein the ion frame. The blue dots represent ω  of the KCSI while the red
dots represent that of the ECDI. The values of plasma parameters areωce /ωpe = 0.5, Ti /Te = 8.0, mi /me = 1836, and βi = 1.0. . . . . . . . . . .

42

55

Growth Rate

Angular Frequency

Figure 1.

Angular Frequency

Growth Rate

Figure 2.

43

8FXG ET

ECDI

KCSI

OKOG
Figure 3.

44

Figure 4.

45

D 8FXGC 8FXG㪋VȁRG

㪎㪅㪌

'Z M

'Z M

㪌

㪉

㫂㺠㩷㪺㪆㱥㫇㪼

VȁRG

㪉㪅㪌
㪇

㪇
㪄㪉㪅㪌

㪄㪉
㪄㪌
㪄㪋

㪇㪅
㪄㪋

㪈

㪄㪉

㪇

㪉㪅
㪉

㪋

VȁRG㪄㪎㪅㪌㪄㪎㪅㪌 㪄㪌 㪄㪉㪅㪌㪉VȁRG

㪇㪅
㪇

㪊㪅
㪉㪅㪌

㪌

㪎㪅㪌

㪈

ȁWJ

㪇

ȁNJ 㪇

ȁEG
ȁNJ

㪄㪇㪅㪌

㪄㪈

㱥㪆㱥㫇㪼

㪇㪅㪌

㪄㪈
㪄㪋

㪇㪅
㪄㪉

㪇

㪇㪅

㪈㪇㪅
㪉

㪋

㪄㪉

㫂㼠㩷㪺㪆㱥㫇㪼

㪄㪎㪅㪌

㪄㪌

㪄㪉㪅㪌

㪇

㪉㪅
㪉㪅㪌

㫂㼠㩷㪺㪆㱥㫇㪼
Figure 5.

46

㪌

㪎㪅㪌

Figure 6.

Figure 7.

47

Figure 8.

48

Ǟ

㪙

㱂㪙
㪄㪼㱢㩷㪂㩷㫄㪼㱥㪺㪼㪉㫏㪉㪆㪉㫄㪼㱥㪺㪼㪉㫏㪉㪆㪉㱂㪘㪘TrappingZGǾ

Wave Propagation
Figure 9.

49

UCV

Figure 10.

50tωpe = 80.

5

5

5

2.5

2.5

2.5

0
-2.5

0
-2.5
-5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

5

-2.5

5

-7.5
-7.5 -5 -2.5 0 2.5kӆcêωpe

7.5tωpe = 320.
7.5

5

5

5

2.5

2.5

2.5

-2.5
-5

kц cêωpe

7.5

0

0
-2.5
-5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

5tωpe = 480.

5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

7.5tωpe = 560.

5

2.5

2.5

2.5

kц cêωpe5kц cêωpe

5

0
-2.5
-5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

5tωpe = 720.

5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

7.5tωpe = 800.

5

5

2.5

2.5

2.5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpekц cêωpe5kц cêωpe

7.5

-5

0
-2.5
-5

5

7.5

7.5

tωpe = 880.

7.5

-2.5

5

-2.5

7.5

0

7.5

0

-5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

7.55tωpe = 640.
7.5

-5

7.5

-2.5

7.5

-2.5

5

0

7.5

0

7.5

-5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

7.55tωpe = 400.

7.5kц cêωpekц cêωpe

0

-5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe

7.5tωpe = 240.

kц cêωpekц cêωpe7.5-5kц cêωpetωpe = 160.

7.5kц cêωpekц cêωpetωpe = 0
7.5

0
-2.5
-5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe.0

Figure 11.

5

7.5

-7.5
-7.5 -5 -2.5 0 2.5kӆ cêωpe2.

tωpe = 80.

0.5

0.5

0.5

0

-0.5

-0.5

0
0.5
vӆêc
tωpe = 240.

-1
-1

1

0

-0.5

-0.5

0
0.5
vӆêc
tωpe = 320.

1

-1
-1

1

1

0.5

0.5

0.5

0

-1
-1

vцêc

1

-0.5

0

-0.5

-0.5

0
0.5
vӆêc
tωpe = 480.

-1
-1

1

-0.5

0
0.5
vӆêc
tωpe = 560.

1

-1
-1

0.5

0.5

0.5

-1
-1

vцêc

1

-0.5

0

-0.5

-0.5

0
0.5
vӆêc
tωpe = 720.

-1
-1

1

-0.5

0
0.5
vӆêc
tωpe = 800.

1

-1
-1

0.5

0.5
vцêc

0.5
vцêc

1

-1
-1

0

-0.5

-0.5

0
vӆêc

0.5

-1
-1

1

-0.5

1

-0.5

1

0
0.5
vӆêc
tωpe = 640.

-0.5

1

-0.5

1

0

1

0

0
0.5
vӆêc
tωpe = 400.

-0.5

1

0

-0.5

0

1

vцêc

vцêc

0

vцêc

vцêc

-1
-1

vцêc

1

-0.5

vцêc

tωpe = 160.

1

vцêc

vцêc

tωpe = 0.
1

0
0.5
vӆêc
tωpe = 880.

0

-0.5

-0.5

0
vӆêc

0

0.5

1

-1
-1

-0.5

0
vӆêc

0.5

30

Figure 12.

1Figure 13.

53

Parallel Heating
XG8F

Perpendicular Heating

Figure 14.

54㪛㫆㫇㫇㫃㪼㫉㪄㫊㪿㫀㪽㫋㪼㪻㩷㱥㪇ȁNJȁ ȁRG

㪄㪇㪅㪇㪉
㪄㪇㪅㪇㪋
㪄㪇㪅㪇㪍

-%5+
'%&+

㪄㪇㪅㪇㪏
㪇㪅㪋

㪇㪅㪌

㪇㪅㪍

㪇㪅㪎

8FXG
Figure 15.

55

㪇㪅㪏

㪇㪅㪐

㪈