The Electrostatic Particle-in-Cell method

The Particle-in-Cell method of Plasma Simulation

I. Overview

The basic idea behind the PIC method of plasma simulation is extremely simple and follows closely the intuitive picture of a plasma. The plasma is a collection of charged particles which interact with each other and with external fields The fields obey the Maxwell equations and the particles follow trajectories determined by Newton’s Law with the force given by the Lorentz equation. The difference between a simulated plasma and a real plasma lies in the representation of the charges, the fields and the space-time in which the phenomena occur. A real plasma usually consists of electrons and singly charged positive ions, though the ions may also be multiply charged, or be negatively charged. In a simulation each charged particle is a homogeneous collection of a large number of real-plasma charges; it is always multiply charged and has a greater mass – but the charge to mass ratio remains the same as that of the real charge. The large number of charges in a plasma are replaced by a much smaller number of these particles. Henceforth by the word “particle” we will mean an entity existing only in computer memory instead of the classical point charge which is thought to occupy real space.

The electromagnetic fields in the simulation are not continuous, neither in space nor in time This is a consequence of the discretisation of the spatial dimensions of the system. The physical volume is divided into cells by lines which run parallel to the boundaries. The intersections of these lines define a set of points called mesh points or grid points.

Each mesh point specifies a location to which the fields and charge densities are assigned after solving the discretised field equations and the discretised equation of motion; the cell itself specifies a volume through the boundaries of which the current densities are calculated. That is why this picture is referred to as the Particle-in-Cell Model (PIC). The particles, whose coordinates are continuous, may occupy positions anywhere within the mesh. The forces acting on them are calculated in terms of the fields at the neighboring mesh points. The particles move through the mesh in finite timesteps. During a timestep the fields are kept constant, at the end of which the discretised field equations are solved again to update the field distribution. Since this update is done over the whole plasma at the same time, it means that propagation effects which depend on the dE/dt term in the Maxwell equations above are neglected. The use of mesh points is a standard technique in the numerical solution of differential equations (Forsyth and Wasow, 1960). In particle simulation, however, the mesh acquires an added significance. It sets a lower limit to the spatial resolution of particle-particle and particle-field interactions. In particular, a plasma with a Debye length which is less than the mesh spacing cannot be accurately simulated for many multiples of the electron plasma frequency. The reason for this is that the Debye length is the distance over which charge separation effects can occur in a plasma, and if the mesh spacing is larger than this, these effects will not be reflected in the fields which are calculated and assigned to the mesh points. Another effect which is a consequence of the mesh is that the particles acquire a finite “size”. This is a result of the scheme used to assign the charge density to the mesh points. This is done by first finding out the number of particles which are near a mesh point and then calculating the fraction of the particles’ charge which should be assigned to it. This fraction is usually some function of the distance of a particle from that point. A mesh point thus starts “seeing” a particle at a distance which depends on the nature of the charge assignment scheme (Langdon, 1970).

A direct result of the mesh-mediated particle size is that the coulomb force law between two particles gets modified. As two particles of the same sign approach each other the force on one due to the other at first rises and then falls. When the particles overlap the force is zero. Hence particles can freely pass through each other, leading to the name Cloud-in-Cell (CIC), by which name the PIC model is also called.

This modification in the force is significant only at distances comparable to the Debye length. At large distances the force approaches the Coulomb force. The mesh also introduces a periodicity in the potential which results in an unphysical coupling between the field and particles. This may be just a minor nuisance or it may result in a situation in which the particles keep on acquiring more and more energy from the field, ultimately leading to a meaningless simulation. The storage requirements for a three dimensional electromagnetic particle simulation can just be met by present day desktops computers, but the time needed to follow the particles for, say, many hundreds of plasma periods would be of the order of several weeks in a typical time sharing computer installation.

Many physical plasmas of interest, however, have some properties which can be exploited to reduce the scale of the simulation. If the plasma extends over a large volume, as in outer space, then a simulation of a relatively small region inside it would be representative of the whole plasma. The plasma may have some spatial symmetry, such as being uniform over one dimension, as is found in long discharge tubes, in which case it may be possible to reduce the dimensionality of the simulation. In some cases other assumptions can be made in order to reduce the number of equations which have to be solved, as discussed below.

If we are interested in times and distances which are comparable to the electron plasma frequency and the Debye length, that is short distances and short times, then the current density can be assumed to be zero over the length and time of significance. If, in addition, the magnetic field is constant then we are left with a system called an electrostatic plasma. Representative problems in such plasmas are beam-induced instabilities, damping, diffusion, oscillations, etc.

If, however, we wish to simulate a plasma over long times compared to the electron plasma frequency and long distances compared to the Debye length, then the charge density vanishes because the plasma is neutral on the average. This system is called the magnetohydrodynamic plasma. Representative problems in such plasmas are interaction of the solar wind with the earth’s magnetic field, long wavelength oscillations, etc.

II. The Electrostatic PIC Algorithm

In this section we describe the basic steps which will yield a self-consistent particle simulation model for a two dimensional electrostatic plasma. A plasma of this kind is obtained after making the following assumptions: 1) There is no variation in the plasma and particle quantities (charge density, velocity, position, etc.) and in the electric field along the z-direction. 2) The magnetic field is zero in the x-y plane in which the particles move, and may be constant perpendicular to it. For the present we will assume that it is zero. 3) We restrict our interest to distances of the order of the Debye length. Since charge separation occurs over these distances the charge density is non-zero. 4) We restrict our interest to times of the order of the electron plasma frequency. Since average current densities over these times are zero, there is no self magnetic field. 5) There are no collisions. With these assumptions the Maxwell equations are reduced to the familiar Poisson’s equation which completely specifies the electric field distribution at a given time.

A. Discretisation of space

The first step is to discretise the space in which the plasma exists. As mentioned above, this is done by drawing lines parallel to the boundaries. The intersection of these lines defines the points which represent the discretised space. Assuming cartesian coordinates, the two-dimensional x-y plane becomes an array of points laid out like the elements of a matrix. Each point is specified by a pair of integers [i, j], with j running along the x-axis and i along the y-axis. This defines the mesh; the spacing between the lines is known as the mesh spacing.

The mesh together with the boundaries is known as the “computational box”. In case any other coordinate system is preferred or required the geometry of the computational box has to be changed accordingly.

B. Discretisation of Poisson’s equation

The next step is to discretise Poisson’s equation in the chosen coordinate system.  The differential form of Poisson’s equation in two dimensions, is the well known 5-point formula which expresses the potential of any grid point as the mean of the potentials at the four neighbouring points.

C. Discretisation of the equations of motion

The third step is to discretise Newton’s equation of motion. This is done in three parts. First, an expression is obtained for the acceleration in terms of the velocity; second, the velocity is obtained in terms of the particle position; and third, the electric field is expressed in terms of the potential at the mesh points This is the field at the mesh point at a particular time t, but since the particles are not constrained to reside on mesh points, it is necessary to interpolate the fields from the mesh points to the particle positions.

The interpolation of the electric field from the mesh points to the particle positions is discussed below in Section II.E after dealing with the issue of charge assignment in the next section.

In order to preserve the time-symmetry of the equations of motion, the discretisation results is a set of difference equations which is time-centred, but at the same time, the velocities and positions differ by half a timestep. Because of this the method is called the leap-frog algorithm.

D. Charge assignment

The fourth step is to obtain an expression for the charge density at a mesh point.

This depends on the number of particles which are in its neighborhood and their positions relative to the mesh point. The procedure is to give each mesh point a “domain” of influence. If a particle lies within the domain of a mesh point, part or all of the charge carried by the particle is assigned to that point.

As mentioned earlier, the charge assignment scheme has several effects on the physics of the simulation. Firstly, the particles acquire a finite size; secondly, the force law at short distances gets modified; thirdly, the scheme introduces discontinuities in the plasma which results in “noise”. All these effects depend crucially on the exact nature of the scheme. We discuss below some of the charge assignment schemes in particle simulation (Eastwood and Hockney, 1974; Okuda, 1972; Birdsall, Langdon and Okuda, 1970; Langdon and Birdsall, 1970).

The simplest scheme is called the Nearest Grid Point scheme (NGP). As the name implies, in this scheme all the particles which are within the domain of the mesh point are located and their complete charge allotted to that mesh point. The total charge assigned to the mesh point is divided by the volume of the mesh cell to arrive at the value of p for that mesh point.

It is necessary to ensure that the charge of a particle is not allotted to more than one mesh point, and that there are no particles whose charge does not get assigned. In order to do this the domains must be non-overlapping, and the domains of all the mesh points must completely cover the computational box. The simplest way is to define the boundaries of a domain to be half a mesh spacing in the x- and y-directions from a mesh point. The boundaries thus run parallel to the mesh. NGP has the advantage that it is simple to compute. The particle size in this scheme is the same as the mesh cell. However, it introduces more noise in the simulation than other schemes. When a particle crosses over from the domain of one mesh point into the domain of its neighbor, the charge assigned to the first mesh point reduces and that assigned to the second increases discontinuously. The potential at the two mesh points also suffers from this discontinuity and this results in a noisy simulation. NGP is not commonly used nowadays, since it has been replaced by other schemes which allow a smoother variation in the mesh quantities as a particle cruises through it.

One of these schemes is called the Area weighting scheme, also known as the First-order weighting scheme. In this scheme the charge carried by the particle is distributed amongst the four nearest mesh points. The charge assigned to a mesh point is thus a fraction of the charge carried by the particle.  A mesh point starts to “see” a particle from the moment it enters one of the four cells surrounding it. In this procedure though the mesh points do not have domains assigned to them as for the NGP scheme, domains can be assigned and the same scheme can be interpreted in a more productive way.

Let us consider the particle to have a size and shape which is independent of the mesh and the particular charge assignment scheme but intrinsic to itself. For example, let us assume that the particle has a rectangular shape  with a uniform charge distribution across this shape, and calculate the area of overlap between it and the domains of the same four mesh points.

We assign the charge to each point in the same way as done above, with the difference that now in the numerator of each equation we substitute the area of overlap in place of the product of the relative coordinates. We will get the same value for the charge assigned to each mesh point as before. This point of view has the advantage that now we can think in terms of the particle as a cloud with a distinct shape and a particular charge distribution across that shape. This is the Cloud-in-Cell (CIC) model and here the particle size is in size for the same charge assignment. The advantage we get out of this new point of view is that the force law between two particles, which is actually a consequence of the charge assignment scheme, can be considered to be a consequence of only the particle shapes once the domain of the mesh points is defined. Conceptually this makes a great difference to the way we think about charge assignment, because now a potentially infinite number of schemes can be implemented on the mesh simply by varying the shape of a particle and the charge distribution within it. This approach can be formally described in terms of the particle shape factor , which defines the charge distribution within a particle and the limits of the particle. We use this shape factor to get the charge assigned to a mesh point by simply multiplying it by the total charge on the particle. In general, the charge distribution will not be uniform over the shape of the particle.

E. Field interpolation

With the knowledge of charge assignment on the mesh points, it is possible to solve Poisson’s equation and get the potential at each mesh point. Knowing that, we can calculate the field at the particle position and hence the field acting on it. The interpolation of the electric field from the mesh points to the particle position is the reverse of charge assignment. Since we know the fields in the x- and y-directions at the four mesh points which surround a particle, we can interpolate twice to evaluate the field at the particle position. However the calculation is made simple by the fact that this procedure of interpolation is equivalent to the method of area-weighting discussed above. We apply the same weighting factors to each component of the field at each mesh point to get the field component at the particle position.

III. The Computer Implementation

The implementation of the procedure described above into a computer code is a relatively straightforward job. How many particles are enough? What should be the time step and the mesh spacing? There are some criteria which should be followed so that the plasma is suitably represented. The number of particles should be much greater than the number of mesh cells. This ensures that each cell always contains several particles on the average during the simulation If the number of particles is too low, the simulation will be noisy. The mesh spacing should not be very much greater than the Debye length, so that the potentials set up by charge separation can be resolved.The timestep should be less than the period of the plasma frequency, so that the particle oscillations are faithfully reproduced. A particle should move less than one mesh spacing in one timestep in order to preserve the self-consistency of the simulation. This requirement will not be satisfied by all particles during the simulation as some particles are accelerated more than others, but we can ensure that it is true at the start. The program takes its input from a file which specifies the dimensions of the computational box and the position of the electrodes inside, if any. The potentials on the boundaries and the electrodes are also specified. The plasma is defined by specifying its initial position, size, shape and charge density. The number of particles in each species, their temperature, and the drift velocities are also defined. The accuracy to which the solution to Poisson’s equaation is to be iterated at each time step is also specified.

After reading the input, the program starts by calculating the position of the particles at t = 0 and their velocities at t-halftimestep. This is done by first giving values of velocity and position to each particle inside the computational box, assigning the charges, calculating the fields, and back-stepping the velocities half a timestep. Depending on the input parameters, the initial values of velocity may be distributed within a particular range in order to approximate a velocity distribution function; this is the usual way of initializing a thermal plasma. However, all particles may be given the same velocity, which would be the case of a cold plasma. The particles may be assigned positions at random, or they may be arranged uniformly through the mesh like atoms in a lattice. The initialization of velocity and position is done for each component, x and y. After initializing the particles, we load each mesh point with the appropriate charge densities using the desired charge assignment scheme and solve Poisson’s equation to obtain the potential at each mesh point. Next we scan the array of particles and obtain the electric field components at the position of each particle by using the same weighting scheme as was used for charge assignment. The components of acceleration at the later time for each particle are then calculated. Using this the program finds the components of position for each particle at the later time. If any particle has hit the boundaries or the electrodes, it is removed from the system. This completes one time cycle. At the end of each cycle, the program again finds the solution of Poisson’s equation with the new distribution of charges and continuously repeats the above procedure until the specified number of time steps are over, or all the charges are removed. It may be noted that the particle positions are calculated at integral multiples of the timestep, whereas velocities are found at integral multiples of half a timestep.

The staggering is due to the time-reversibility criterion discussed above. This method is known as the leap-frog algorithm and is accurate to second order.

IV. References

Birdsall, Langdon, and Okuda, 1970. “Finite sized particle Physics applied to Plasma simulation”, by C. K. Birdsall, A. B. Langdon and H Okuda, pages 241-258 in Methods in Computational Physics, volume 9 (Plasma Physics), edited by B . Alder, S. Fernbach, and N. Rotenberg, Academic Press, New York, 1970

Eastwood, Hockney, 1974. “Shaping the force law in Two-Dimensional particle Mesh models”, by J. W. Eastwood and R. W. Hockney, J. Comp. Phys, 16, 342, 1974.

Forsyth, Wasow, 1960. “Finite Difference Methods for Partial Differential Equations”, by G E Forsyth and W R Wasow, Wiley, New York, 1960.

Langdon, A. B., 1970. “Effects of the spatial grid in simulation Plasmas”, by A. B. Langdon, J. Comp. Phys, 6, 247, 1970.

Langdon, Birdsall, 1970. “Theory of Plasma simulation using finite sized particles”, by A. B. Langdon and C. K. Birdsall, Phys Fluids, 13, 2115, 1970.

Okuda, 1972. “Non-physical noise and instabilities in Plasma simulation due to a spatial grid”, H. Okuda, J. Comp Phys, 10, 475, 1972.

Patel, 1989. “Particle Simulation of a Two-dimensional electrostatic plasma”, by Kartik Patel, Report No.:BARC- 1489, 1989.

Patel, 1992. “Concurrent Particle in cell plasma Simulation on a multitransputer parallel computer”, by Kartik Patel, A. N. Khare and A. N. Jethra, Report No.:BARC/1992/E/015, 1992.

Patel, 1994. “Parallel Particle in Cell Plasma Simulation”by K. Patel, Seminar cum Workshop on Supercomputing for Scientific Visualization, BARC, February15-17, 1994.

Patel, 1995. “Two-dimensional particle simulation of plasma expansion between plane parallel electrodes”, by Kartik Patel and V.K.Mago, J.Appl.Phys, vol 78, 4371 (1995)

Patel, 1997. “Observation of ion waves in two-dimensional particle simulation of field assisted plasma expansion”, Kartik Patel, J. Appl.  Phys., vol 81, 6622 (1997)

Jackson, 1978. “Classical Electrodynamics”, by J. D. Jackson, Wiley-Eastern, New Delhi, 1978.