Accurately charge-conserving scheme of current assignment based on the current continuity integral equation for particle-in-cell simulations

Abstract

In this paper we present a high-accuracy charge conservation scheme of current assignment for particle-in-cell (PIC) simulations in two dimensional space. The current continuity integral equation is applied to set up a linear equation set about the assigned currents on edges of the cells through which a particle moves during a time step. The currents on the edges can be very accurately assigned by solving the linear equation set. This scheme can be used in Cartesian and cylindrical geometry, and is not affected by the form-factor of quasi-particle. It has been applied to the PIC simulation of Z-pinch of the current-carrying rarefied deuterium plasma shell without resorting to the electric field correction by solving Poisson equation. The simulated results reveal that the relative and absolute errors in satisfying the current continuity integral equation and Poisson equation are very tiny.

Introduction

The principal algorithms of particle-in-cell (PIC) simulation have been well established [1], [2], [3] and extensively used for several decades in many kinds of plasma researches and applications, such as plasma focus and Z-pinch [4], [5], laser inertial confinement fusion (ICF) [6], [7], laser plasma interaction [7], [8], [9], as well as high-power microwave generation [10], just to name a few. The PIC simulations can provide with insights and information which could not be obtained in magnetohydrodynamic simulations and related experiments.

In the traditional procedure of PIC simulation, the electromagnetic fields are advanced by solving the curl equations (i.e. electromagnetic induction equation and Ampere’s law with displacement current). That possibly results in the obtained electric field not to satisfy adequately the Poisson equation, and even in the spurious modes with secular growth. Of course, this problem can be solved by means of repeatedly solving Poisson equation to find a correction to the obtained electric field. However, that will spend much more computer resource and make the program parallelization difficult due to the global character of divergence equations ($\nabla \cdot \overrightarrow{E}=4\pi \rho$ and $\nabla \cdot \overrightarrow{B}=0$) for their solution. Fortunately, one can readily check that $\nabla \cdot \overrightarrow{E}=4\pi \rho$ will always be true if the initial electric field satisfies the Poisson equation and the current continuity equation

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \overrightarrow{\mathrm{J}}=0,$$

where $\rho$ is charge density and $\overrightarrow{J}$ is current density, is maintained to be true during the simulation. One can also check that $\nabla \cdot \overrightarrow{B}$ will remain to be zero if it is same initially. Therefore, the traditional method for advancing the fields in the PIC simulation will be possibly valid after the divergence equations are imposed only as the initial conditions [11].

Thus how to maintain the current continuity equation (1) (i.e. the difference equation form of charge conservation law) to be satisfied during the PIC simulation becomes an important issue. As a particle is pushed forwards in a PIC simulation, a current in the grid is created. Some schemes of the current assignment to edges of the related cells in every time step had been proposed for the goal to ensure that Eq. (1) could be satisfied as accurately as possible [11], [12], [13], [14], [15], [16], [17], [18], [19], and in turn to allow us to avoid solving Poisson equation. They are called “charge conservation methods”. The most general derivations of charge conservation methods were presented in [11], [12] for Cartesian coordinates and [13], [19] for curvilinear meshes in one, two and three dimensions. Basically, the simulated particle (quasi-particle) trajectory is assumed to be a straight line in one time step and split into a few segments which lie appropriately in a few different cells. The current densities for the segments are assigned to the edges of the related cells with some fixed assignment pattern based on the current continuity difference equation and some assignment weightings. Whereas Umeda and his coauthors developed a new charge conservation scheme where a particle trajectory was assumed as a zigzag line [16]. In Esirkepov’s density decomposition scheme [15], a charge flux is decomposed into movements of four particles along the x and y axes in two dimensions, and of twelve particles along the x, y and z axes in three dimensions. Each segment has its own weighting, and the sum of the charge flux segments is equal to $q\overrightarrow{\upsilon }$ ($q$ the charge of particle, $\overrightarrow{\upsilon }$ the velocity of particle). Esirkepov’s scheme is suitable for application to Cartesian geometry and arbitrary higher-order shape-factors. Its computation speed is faster than that of the Villasenor–Buneman method [11], [16]. The satisfying degree of Poisson equation for Villasenor–Buneman charge conservation method was checked [20] in a laser plasma interaction PIC simulation. The relative error

$$g=(\nabla \cdot \overrightarrow{E}-4\pi \rho )∕4\pi \rho$$

is better than 0.6% with typical boundary condition amendment during the simulation.

In this paper, we present a new high-accuracy charge conservation method based on the application of current continuity integral equation to the nodes, at which the charge density is defined, for a Z-pinch PIC simulation in two dimensional space. According to the assumption of particle-in-cell algorithm, in which the charged particle movement in a cell contributes only to the physical quantities of this cell, and to the variable arrangement of Yee grid [21], a linear equation set with respect to the currents on all edges of the related cells, through which a particle moves in a time step, can be set up for all nodes of the cells. The generated current on any edge of the cells can be accurately assigned due to the particle movement by solving the linear equation set. This scheme is suitable for Cartesian and cylindrical geometry, and does not need any special current assignment weighting.