Development and Test of 2.5Dimensional Electromagnetic PIC Simulation Code
 Author: Lee SangYun, Lee Ensang, Kim KhanHyuk, Seon Jongho, Lee DongHun, Ryu KwangSun
 Publish: Journal of Astronomy and Space Sciences Volume 32, Issue1, p45~50, 15 March 2015

ABSTRACT
We have developed a 2.5dimensional electromagnetic particle simulation code using the particleincell (PIC) method to investigate electromagnetic phenomena that occur in space plasmas. Our code is based on the leapfrog method and the centered difference method for integration and differentiation of the governing equations. We adopted the relativistic BunemanBoris method to solve the Lorentz force equation and the Esirkepov method to calculate the current density while maintaining charge conservation. Using the developed code, we performed test simulations for electron twostream instability and electron temperature anisotropy induced instability with the same initial parameters as used in previously reported studies. The test simulation results are almost identical with those of the previous papers.

KEYWORD
particle simulation , PIC method , twostream instability , whistler instability

1. INTRODUCTION
The ParticleInCell (PIC) method is a technique commonly used to simulate the dynamics of charged particles interacting with electromagnetic fields. Since the 1950s, the PIC method has been used for selfconsistent kinetic simulations along with the Vlasov simulation method. The PIC method is useful to study the dynamics of plasmas that takes place on spatial and temporal scales that magnetohydrodynamics (MHD) cannot handle properly.
In this paper, we present the development of a 2.5dimensional relativistic electromagnetic PIC code. The code is developed to study electromagnetic phenomena that occur in space plasmas. To test this newly developed code, we performed simulations for both electrostatic and electromagnetic phenomena and compared the results with those of two widely used codes (Umeda 2004, Gary et al. 2011). For the test of the electrostatic phenomena, we applied the code to electron twostream instabilities. For the electromagnetic phenomena, whistler instabilities driven by electron temperature anisotropy were simulated.
In section 2, we briefly introduce the numerical schemes adopted in the code, and in section 3, we present test simulation results in comparison with those reported previously. Discussions and conclusions are presented in section 4.
2. METHOD
To handle various plasma processes we solve Maxwell’s equations and relativistic Lorentz force equation as follows:
where
E ,B ,J , andρ are electric field, magnetic field, current density, and charge density, respectively. Lorentz equation,q _{i,}m _{i,}v _{i}, andγ _{i} are charge, rest mass, velocity, and Lorentz factor of thei th particle, respectively. Maxwell’s equations describe the variation of electromagnetic fields, which the motion of charged particles through Lorentz force equation. Particles contribute to Maxwell’s equations by the charge and current densities.When particles are initially loaded into the simulation domain using uniform random distribution, electric field can be produced by the nonuniformity in charge density caused by the random distribution. Thus, we first solve Poisson equation,
where
φ is the electrostatic potential used to obtain the initial electric field. To solve the Poisson equation we use the successive overrelaxation (SOR) method (Press et al. 1992)Next, we solve Faraday’s induction equation, Eq. (2), using the leadfrog scheme to advance the magnetic field in time delineated by the following discretized equation.
Using the electric and magnetic fields we solve the relativistic Lorentz force equation, Eq. (5). As a result of using the leapfrog scheme, as Eq.(7), the electric field and the magnetic field are determined at different times: the electric field at (n + 1) Δ
t , and the magnetic field at at (n + 1)Δt . To obtain the values at the same time, (n + 1)Δt , we average the magnetic fields at nΔt and (n + 1)Δt . To solve the relativistic Lorentz (Birdsall & Langdon 2004).Particles are then moved to new positions by sloving the following equation:
Finally, we solve the AmpèreMaxwell equation, Eq. (1), using the current density and the magnetic field to advance the electric field in time. The estimation of the current density is complicated and many methods have been suggested to correctly calculate the current density (Villasenor & Buneman 1992, Esirkepov 2001, Umeda et al. 2003). Generally, wwe can calculate the current density as
= ∑_{i}J n_{i}q_{i}v_{i} , However, this does not satisfy the charge conservation in discrete grid systems. To satisfy the charge conservation, we have to correct the electric field using some method, or the electric field must have a divergencefree part only. Alternatively, various methods have been suggested to satisfy the charge continuity equation as given below:To slove the Poisson equation and calculate the current densiy it is necessary to estimate the charge density onthe grids from particle distributions The contribution of a particle to the charge density on a grid is determined by using shape functions. In this code, we use secondorder shape function in twodimensions, given as follows:
Using the secondorder shape function, one particle contributes to 9 adjacent grids. The higherorder shape function has the advantage of reducing numerical errors, but at the cost large calculation time.
The time integration of all the differential equations is achieved by using the leapfrog scheme, as mentioned above. To carry out the leapfrog scheme, all the variables that have time dependence are divided into the variables defined at the fullinteger time step, (n + 1)nΔ
t . For example, the magnetic field, current density, and velocity are defined at the fullinteger time step while the electric field and position at the halfinteger time step.For the differentiation in space we use the centered difference scheme. In this scheme, the variables are defined either at a fullinteger grid or at a halfinteger grid, which is called a staggered grid system. For example, the charge density, electrostatic potential, and
E_{z} are defined at the fullinteger grid both in the x and ydirections. On the other hand,E_{x} andB_{y} are defined at the halfinteger grid in the xdirection and the fullinteger grid in the ydirection,E_{y} andB_{x} are defined at the fullinteger grid in the xdirection and the halfinteger grid in the ydirection, andB_{z} is defined at the halfinteger grid both in the x and ydirections.3. TEST OF THE DEVELOPED CODE
To test the developed code, we simulated the electron twostream instability and the whistler instability using the same parameters as described in previously reported papers and compared the results with those of the previous papers.
3.1. Electron TwoStream Instability
First, we simulated the electron twostream instability with the parameters used by Umeda (2004). In this simulation, the length is normalized by the Debye length, λ
_{De} , and time by the electron plasma frequency,w_{pe} , satisfying λ_{De} =v_{the} /w_{pe} , where is the electron thermal speed. The system size isL_{X} ×L_{Y} = 256λ_{De} × 256λ_{De} , and the number of grids isN_{X} ×N_{Y} = 256 × 256, which makes the grid size Δx = Δy = 1.0λ_{De} . The number of particles in a cell is set to be 128 for the electrons and protons, respectively. The ratio of the electron plasma frequency and the electron gyrofrequency is set to bew_{pe} / Ω_{pe} , and the electron thermal velocityv_{the} = 0.05c . The background magnetic field is oriented in the Xdirection , and the magnitude of the magnetic field,B_{0} , is determined by the electron gyrofrequency. Boundary conditions are set to be periodic in both the X and Ydirections.To excite the electron twostream instability, electron beams drifting oppositely to each other were initially injected uniformly throughout the simulation domain. The densities of the two beams were the same. The bulk speeds of the beams were ±2
v_{the} , which makes the relative speed between the beams 4v_{the} . To make the simulation conditions identical with Umeda (2004), we assumed that protons are immobile. This assumption is acceptable because the mass of protons is much larger than that of electrons and the electron twostream instability develops in the time scale of electrons.The simulation results are presented in Figs. 1 and 2. Fig. 1 shows the temporal variations of the electric field energy density,
2>, magnetic field energy density, 2>, kinetic energy density, 2/2>, and total energy density averaged over the whole simulation domain. Because the twostream 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 twostream instabilities. After the electric field energy reached the maximum, it gradually decreased while the kinetic energy increased. 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
w_{pe}t = 60,w_{pe}t = 260, andw_{pe}t = 1020. Because the two electron beams drift oppositely along the background magnetic field, the electrostatic waves propagate parallel toB_{0} producing sinusoidal oscillations of the potential (Fig. 2a), which are represented as electron holes in the phase space. After the instability reached saturation atw_{pe}t ~ 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.3.2. Electron Temperaturedriven Whistler Instability
Next, we followed the initial configurations of Run 3 in Gary et al. (2011), who investigated the electron temperaturedriven whistler instability using their PIC code. The simulation system is normalized by λ
_{e} in space andw_{pe} in time, satisfying λ_{e} =c /w_{pe} , where λ_{e} is the electron inertial length. The system size isL_{X} ×L_{Y} = 12.85λ_{e} × 12.85λ_{e} , the number of grids isN_{X} ×N_{Y} = 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 gyrofrequency is set to bew_{pe} / Ω_{ce} = 4.0, where Ω_{ce} =eB /m_{e} is the electron gyrofrequency. Ions are assumed to be an immobile background as in the previous case. Periodic boundary conditions are used for both the x and ydirections. The background magnetic field is applied along the xdirection, .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 theB_{z} component was dominant, which corresponds to the growth of theE_{x} andE_{y} 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.Fig. 4. shows the distributions
B_{z} (x ,y ) in the configuration space andB_{z} (k _{},k _{⊥}) in the kspace atω_{pe}t = 740 when the instability is close to saturation. It is seen that the whistler waves propagate obliquely. The spatial distribution of theB_{z} 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.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}/v_{the} ≤ 6 in the direction parallel to the background magnetic field. The generation of suprathermal electrons in the parallel direction reduced the temperature anisotropyT _{⊥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 nonMaxwellian 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).4. DISCUSSION AND CONCLUSION
In this paper, we briefly introduced the 2.5dimensional 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 twostream 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 higherorder shape functions or by implementing other numerical schemes.

[]

[]

[]

[]

[]

[]

[]

[]

[]

[]

[Fig. 1.] 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

[Fig. 2.] Distribution of electrostatic potential at (a) wpet= 60, (b) wpet = 260, and (c) wpet= 1020.

[Fig. 3.] 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, zcomponents and the magnitude, respectively.

[Fig. 4.] (a) Spatial distribution of Bx, and (b) the corresponding spectrum in kspace at wpet = 720.

[Fig 5.] 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).