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.
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.
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.
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.
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.
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.
Last updated: 1 May, 2004
Copyright (c), 2001-2004
Kartik Patel