2
/2>, and total energy density averaged over the whole simulation domain. Because the two-stream instability is an electrostatic instability, the electric field energy rapidly increases as the instability grows while the magnetic field energy remains almost constant. The increased electric field energy comes from the kinetic energy of particles, which decreases in coincidence with the increase of the electric field energy. This is consistent with the theory of the two-stream instabilities. After the electric field energy reached the maximum, it gradually decreased while the kinetic energy increased.
Temporal variations of electric field energy density (red), magnetic field energy density (blue), kinetic energy density (green), and total energy density (black) averaged over the whole simulation domain. The values are normalized by
The drastic variations of energies in the growth phase of the instability produced numerical errors, which caused a loss of the total energy. However, during gradual variations of energies, the total energy was well conserved. In comparison with the results of
Umeda (2004)
, the variations of the energies are similar except for that in the case of
Umeda (2004)
the total energy density gradually decreased even after the decrease in the rapid growth phase.
Fig. 2
shows the distribution of the electrostatic potential at
wpet
= 60,
wpet
= 260, and
wpet
= 1020. Because the two electron beams drift oppositely along the background magnetic field, the electrostatic waves propagate parallel to
B0
producing sinusoidal oscillations of the potential (
Fig. 2a
), which are represented as electron holes in the phase space. After the instability reached saturation at
wpet
~ 60, the wave structures or electron holes merged together until there remained only one structure. In comparison with the results of
Umeda (2004)
, the wave structures are very similar in the early development stage in
Fig. 2a
. In
Figs. 2b
and
2c
, however, the merging of the structures was slightly faster in our results than in
Umeda (2004)
, which might be related to the difference in the variations of the energies.
Distribution of electrostatic potential at (a) wpet= 60, (b) wpet = 260, and (c) wpet= 1020.
- 3.2. Electron Temperature-driven Whistler Instability
Next, we followed the initial configurations of Run 3 in
Gary et al. (2011)
, who investigated the electron temperature-driven whistler instability using their PIC code. The simulation system is normalized by λ
e
in space and
wpe
in time, satisfying λ
e
=
c
/
wpe
, where λ
e
is the electron inertial length. The system size is
LX
×
LY
= 12.85λ
e
× 12.85λ
e
, the number of grids is
NX
×
NY
= 64 × 64, and the number of particles is set to be 9600 particles in a cell for each of the electrons and protons. The ratio of electron plasma frequency and electron gyro-frequency is set to be
wpe
/ Ω
ce
= 4.0, where Ω
ce
=
eB
/
me
is the electron gyro-frequency. Ions are assumed to be an immobile background as in the previous case. Periodic boundary conditions are used for both the x- and y-directions. The background magnetic field is applied along the x-direction,
.
To excite the instability, the electron temperature anisotropy in the perpendicular and parallel directions to the background magnetic field is applied as .
T
⊥e
/
T
||e
, and the parallel plasma beta of electrons is given as
β
||e
= 0.01.
Fig. 3
. shows the temporal variations of the temperature anisotropy,
T
⊥e
/
T
||e
, magnetic energy density,
, and electric energy density, δ
E
2
/ (
cB
0
)
2
, of the fluctuations. As the instaility grew, the temperature anisotropy decreased rapidly while the electric and magnetic energy densities increased. The kinetic energy density decreased in coincience (not shown). Note that the magnetic energy density is an order of magnitude larger than the electric energy density. The growth of the
Bz
component was dominant, which corresponds to the growth of the
Ex
and
Ey
components. In comparison with results of
Gary et al. (2011)
, all the variations are almost the same except for that our results evolve slightly faster. It should be noted that time was normalized by the electron cyclotron frequency in
Gary et al. (2011)
, but in our simulation it was normalized by the electron plasma frequency. The ration between the frequencies is ω
pe
/ |Ω
e
| = 4, where ω
pe
and Ω
e
are the electron plasma frequency and the cyclotron frequency, respectively.
Temporal variations of (a) temperature anisotropy, (b) magnetic energy density, and (c) electric energy density. Red, green, blue, and black colors in (b) and (c) represents the x-, y-, z-components and the magnitude, respectively.
Fig. 4
. shows the distributions
Bz
(
x
,
y
) in the configuration space and
Bz
(
k
||
,
k
⊥
) in the k-space at
ωpet
= 740 when the instability is close to saturation. It is seen that the whistler waves propagate obliquely. The spatial distribution of the
Bz
component and the propagation characteristics of the whistler waves in our simulation results are almost identical with those of
Gary et al. (2011)
. The only difference is that the wave number, k, of our results is slightly smaller than that of
Gary et al. (2011)
, which implies that the wavelength of the generated whistler waves is slightly longer in our results.
(a) Spatial distribution of Bx, and (b) the corresponding spectrum in k-space at wpet = 720.
Fig. 5
shows the parallel velocity distribution of electrons at the initial and final stages of the simulation. As a result of the instability, suprathermal electrons were generated in the range of 4 ≤ |
v
||e
|/
vthe
≤ 6 in the direction parallel to the background magnetic field. The generation of suprathermal electrons in the parallel direction reduced the temperature anisotropy
T
⊥e
/
T
||e
. However, it should be noted that the temperature is not well defined for the distribution seen in
Fig. 5
because the distribution is a non-Maxwellian distribution consisting of the core and suprathermal components. According to
Gary et al. (2011)
, the suprathermal electrons were produced by the acceleration by the parallel component of the electric field through Landau resonance. The velocity distribution of electrons is also consistent with that of
Gray et al. (2011)
.
Velocity distribution of electrons in the direction parallel to the background magnetic field. The dashed line shows the initial velocity distribution and the soild line the distribution at wpet = 7400 (Ωe t = 1850).
4. DISCUSSION AND CONCLUSION
In this paper, we briefly introduced the 2.5-dimensional relativistic electromagnetic PIC code that we have developed and tested the performance of the code by comparing the results obtained from the code with those of two widely used codes reported in previous studies.
First, we performed a simulation of the electron two-stream instability to test the availability of handling the plasma dynamics in electrostatic phenomena. Our simulation results captured most of the physical features presented in
Umeda (2004)
except that merging of the electron holes in the relaxation phase took place slightly faster in our simulation. Note that both codes suffered similar numerical errors during the rapid growth of the instability. The loss of the total energy appears to be related with numerical diffusion and our code has slightly larger numerical diffusion than Umeda’s. This also appears in the merging rate of the electron holes, which is faster in our results than those of Umeda. Thus, we should be cautious for unphysical loss or gain of energies when a process occurs rapidly. However, the loss or gain of energies was not significant enough to affect the main physical processes.
Next, we performed a simulation of the whistler instability driven by electron temperature anisotropy to test for handling the plasma dynamics in electromagnetic phenomena. In this case, our simulation successfully produced the whistler instability and the results of our simulation were also almost identical with those of
Gary et al. (2011)
. A slight difference was found in the wavelength of the generated whistler waves. However, this did not substantially affect the main features of the instability.
In conclusion, the PIC code we have developed successfully reported in previous studies for both electrostatic and electromagnetic phenomena. We will further improve the code in future studies to reduce the numerical errors by applying higher-order shape functions or by implementing other numerical schemes.
Acknowledgements
This work was in part supported by the Basic Science Research Program (NRF-2013R1A1A2010711) and the BK21 Plus Program through the National Research Foundation (NRF) funded by the Ministry of Education of Korea, and the grant CATER2012-5060 from the Center for Atmospheric Science and Earthquake Research in Korea.
Birdsall CK
,
Langdon AB
2004
Plasma Physics via Computer Simulation
CRC Press
New York
Esirkepov TZ
2001
Exact Charge Conservation Scheme for Particle-In-Cell Simulation with an Arbitrary Form-factor
Comput. Phys. Comm.
http://dx.doi.org/10.1016/S0010-4655(00)00228-9
135
144 -
153
DOI : 10.1016/S0010-4655(00)00228-9
Gary SP
,
Liu K
,
Winske D
2011
Whistler anisotropy instability at low electron β: Particle-in-cell simulations
Phys. Plasma.
http://dx.doi.org/10.1063/1.3610378
18
082902 -
DOI : 10.1063/1.3610378
Press WH
,
Teukolsky SA
,
Vetterling WT
,
Flannery BP
1992
Numerical Recipes in Fortran 77
Cambridge University Press
New York
Umeda T
,
Omura Y
,
Tominaga T
,
Matsumoto H
2003
A new charge conservation method in electromagnetic particle-in-cell simulations
Comput. Phys. Comm.
http://dx.doi.org/10.1016/S0010-4655(03)00437-5
156
73 -
85
DOI : 10.1016/S0010-4655(03)00437-5
Umeda T
2004
Study on Nonlinear Processes of Electron Beam Instabilities via Computer Simulations, PhD Dissertation
Kyoto University
Villasenor J
,
Buneman O
1992
Rigorous charge conservation for local electromagnetic field solvers
Comput. Phys. Comm.
http://dx.doi.org/10.1016/0010-4655(92)90169-Y
69
306 -
316
DOI : 10.1016/0010-4655(92)90169-Y