We have developed a 2.5-dimensional electromagnetic particle simulation code using the particle-in-cell (PIC) method to investigate electromagnetic phenomena that occur in space plasmas. Our code is based on the leap-frog method and the centered difference method for integration and differentiation of the governing equations. We adopted the relativistic Buneman-Boris 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 two-stream 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.
The Particle-In-Cell (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 self-consistent 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.5-dimensional 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 two-stream 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.
To handle various plasma processes we solve Maxwell’s equations and relativistic Lorentz force equation as follows:
where
When particles are initially loaded into the simulation domain using uniform random distribution, electric field can be produced by the non-uniformity in charge density caused by the random distribution. Thus, we first solve Poisson equation,
where
Next, we solve Faraday’s induction equation, Eq. (2), using the lead-frog 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 leap-frog scheme, as Eq.(7), the electric field and the magnetic field are determined at different times: the electric field at (n + 1) Δ
Particles are then moved to new positions by sloving the following equation:
Finally, we solve the Ampère-Maxwell 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
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 second-order shape function in two-dimensions, given as follows:
Using the second-order shape function, one particle contributes to 9 adjacent grids. The higher-order 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 leap-frog scheme, as mentioned above. To carry out the leap-frog scheme, all the variables that have time dependence are divided into the variables defined at the full-integer time step, (n + 1)nΔ
For the differentiation in space we use the centered difference scheme. In this scheme, the variables are defined either at a full-integer grid or at a half-integer grid, which is called a staggered grid system. For example, the charge density, electrostatic potential, and
To test the developed code, we simulated the electron two-stream 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.
First, we simulated the electron two-stream instability with the parameters used by Umeda (2004). In this simulation, the length is normalized by the Debye length, λ
To excite the electron two-stream 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
The simulation results are presented in Figs. 1 and 2. Fig. 1 shows the temporal variations of the electric field energy density,
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
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 λ
To excite the instability, the electron temperature anisotropy in the perpendicular and parallel directions to the background magnetic field is applied as .
Fig. 3. shows the temporal variations of the temperature anisotropy,
Fig. 4. shows the distributions
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 ≤ |
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.