# The Electrostatic Particle In Cell (ES-PIC) Method

Posted on November 8th, 2010
Next Article

This article first appeared on the website of Micropropulsion and Nanotechnology Laboratory (MpNL), my research group at the George Washington University. It’s reprinted here for convenience.

Particle-In-Cell (PIC) is a technique commonly used to simulate motion of charged particles, or plasma. Here I show you the basics of the method, as applicable to studies of phenomena such as solar wind propagation, or analysis of electric thruster plumes. These discharges are characterized by low plasma densities. At low density, plasma behaves more like a collection of discrete particles, than a single continuous fluid. High-density plasmas are simulated using the extension of computational fluid dynamics into electromagnetics, magnetohydrodynamics.

Modeling of plasmas is complicated by the presence of external and self-induced electromagnetic fields, inter-particle interactions, presence of solid objects, and the different characteristic time scales at which ions and electrons propagate. In order to speed up the calculation, simplifying assumptions are generally taken to match the problem at hand. Here, we are assuming that the current generated by the plasma medium is low enough so that the self-induced magnetic field can be ignored. This is a valid assumption for the class of problems we are dealing with here. As such, the set of underlying Maxwell’s equations is reduced and we obtain an Electro-Static PIC code. In addition, we assume that electrons follow the Boltzmann relationship. The simulation then consists of only the heavy particles, ions and neutrals. This simplification has a tremendous impact on the computational speed, since the time integration can be performed on the much larger ion time scale. Finally, we assume that gas densities are sufficiently low so that particle collisions are negligible.

So what is plasma? It is a collection of charged particles (ions and electrons) as well as non-ionized atoms that is globally charge neutral. In other word, assume that you could grab a chunk of plasma. As long as the radius of this volume is larger than a certain characteristic length (called the Debye lenght), the collection will contain approximately equal amount of positive (from ions) and negative (from electrons and negative ions) charge. The ratio of ionized particles to neutral atoms is the ionization fraction. In
electric thrusters, about 90% of the working gas is ionized.

Charged particles interact with each other by attracting particles of opposite charge and repelling those with the same charge. This is the Coulomb force and is given by: . Conceptually, we could simulate plasma by taking a collection of particles representing the real physical ions and electrons and directly computing this force. However, plasma simulations generally require at least 1 million particles in order to reduce numerical errors. Since Coulomb force leads to an n2 problem, computation of a single time step would require at least 1 trillion operations. Quite a daunting task, especially considering that a typical simulation requires several thousand of such steps!

Introducing the Particle-In-Cell (PIC) method. This method also uses computational particles (called macro-particles) to represent the real ions, electrons and neutrals. However, instead of computing the Coulomb force directly, we obtain the force acting on the particles from the electric field. Force acting on the particles is given by the Lorentz force, , which is a size-n problem. Here q is the particle charge, E is the electric field, and v and B are the particle velocity and magnetic field, respectively. In this example we will assume that plasma is not subjected to any magnetic field (self-induced or applied) and thus the second term in the Lorentz force equation can be ignored. As was eluded to previously, it is not computationally feasible to simulate every physical ion, electron or neutral. As such, each computational particle represents a large number (several tens of million in thruster simulations) of real particles. The ratio of real particles per macro-particle is called the specific weight, sp_wt.

## Fundamental Equations

Motion of each macroparticle is governed by Newton’s Second Law,

Here q is the particle charge and m is the particle mass. Electric field is computed from , where $\phi$ is the electric potential. It is given by , the Poisson equation. Here ρ is the charge density and ε0 is the permittivity of free space. In terms of ion and electron number densities, charge density is written as . The subscripts i and e denote ions and electrons, respectively, and Zi is the average ion charge number.

Due to their tiny mass, electrons move at very high velocities (on the order of 5×105 m/s in the types of plasmas considered by MpNL). Resolving electron motion without introducing numerical errors requires the use of extremely small computation time steps (around 10-12 s). With such a small time step, it is not feasible to study any large-scale problems. For instance, it takes couple tens of µs for an ion beam to propagate several thruster lengths from the spacecraft. At the electron time scale, such a simulation would require 10 million time steps. In addition, from the frame of reference of the ions, electrons move instantaneously. Hence it makes sense to simplify the analysis by considering electrons as a fluid. The electron density is then given by the Boltzmann relationship,

## The Algorithm

The Particle In Cell method consists of an initial setup, the main loop, and a final clean up / results output. All computation happens in the loop. The loop consists of the following steps:

1. Compute Charge Density: particle positions are scattered to the grid
2. Compute Electric Potential: performed by solving the Poisson equation
3. Compute Electric Field: from the gradient of potential
4. Move Particles: update velocity and position from Newton’s second law.
5. Generate Particles: sample sources to add new particles
6. Output: optional, save information on the state of simulation
7. Repeat: loop iterates until maximum number of time steps is achieved or until simulation reaches steady state

### Compute Charge Density

Charge density is a scalar variable spatially varying in space. It indicates the number of charge units per unit volume. We compute it by distributing charge of all particles onto the nodes of computational cells, and then dividing by the corresponding node volume.

The 1st order (linear) scatter operation is schematically shown in Figure 2. This figure shows a snippet of the computational mesh, and the simulation particle, shown in gray. Charge of the particle is distributed among the nodes of the cell the particle is in. Hence then name for this method, particle in cell. Please note that the total charge carried by each macro-particle is sp_wt*q. The yellow node receives the largest fraction, due to the particle’s close proximity to this node. The smallest amount is contributed to the green node. The weight factors are given by the four area fractions,

where 1, 2, 3 and 4 correspond to the blue, yellow, green and pink nodes, respectively. hx is the fractional distance of the particle from the cell origin in the x direction. Node volumes equal the cell volumes, ΔxΔy. Exceptions exist along domain boundaries, since here only a half (or quarter for corners) of cells is contributing. More complicated boundary conditions require additional treatment. For instance, periodic boundaries require summation of nodes on the plus (+) and the negative (-) side of the domain. This interpolation technique can also be used to scatter data with non-rectangular cells by mapping the physical coordinates to a regular logical domain.

### Compute Electric Potential

Next, the electric potential is calculated by calling an elliptic equation solver. The Poisson equation can be discretized as

where we used the finite difference with central differencing. This method is modified along mesh boundaries. An elliptic equation describes a boundary value problem and as such boundary conditions must be prescribed along all boundaries. Two types of boundaries exist: Dirichlet and Neumann. The first, Dirichlet, specifies the value of potential. The second boundary specifies the value of the derivative, or electric field. Plasma problems will typically contain a combination of both boundary types. For instance, Figure 3 shows the boundary conditions for the flow of solar wind around a charged plate problem simulated by the sample MATLAB PIC code.

### Compute Electric Field

Electric field is computed by differencing electric potential, or along boundaries . Only the x direction derivative is shown here, and the one sided boundary equation is applicable to the x- face. The remaining derivatives are obtain in an analogous fashion.

### Move Particles

We next integrate particle motion through a time step Δt. The Leapfrog method is commonly used in PIC codes. It is fast and numerically stable. First, we integrate velocity through the time step. Next, position is updated. The name comes from the fact that times at which velocity and positions are known are offset from each other by half a time step. As such, the two quantities leap over each other. This idea is sketched in Figure 4.

Numerically, the particle push is given by

The electric field at the position of each particle is computed by a gather operation, which is inversely analogous to the process used to compute charge density.

After particles are moved to new positions, it is necessary to verify that all particles are still in the computational domain. Two boundary interactions are possible. The particles can either exit the domain, or can collide with solid objects. Computational boundaries are either open (or absorbing, allowing particles to leave), reflective (elastically returning particles into the domain) or periodic
(particles are transported to the opposite side of the domain). The reflective boundary is used to identify planes of symmetry. Details of particle-surface interaction, which can result in erosion of native material, is not considered here.

### Generate Particles

New particles are generated by sampling sources. Sources include spacecraft thrusters or solar wind. This step is problem specific.
In some cases, all particles will be loaded initially, and the simulation will only computes their final distribution state. Typically we are interested in sampling the Maxwellian distribution. As given in Birdsall, this distribution can be approximated as:

where M is some number and Ri is an i-th random number in range [0:1). The value chosen for M controls the accuracy of this method. Birdsall used 12 to prevent entries larger than 6 times the thermal velocity. In the sample MATLAB code, we use M=3.

Finally, since the Leap-Frog method requires velocity and position to be offset from each other, it is necessary to update the velocity of sampled particles by -Δt. This step is frequently omitted. The size of the numerical error introduced by this omission will depend on several factors: the strength of the electric field near the source and size of the time step. This step is also omitted in the sample code, for clarity as it would require duplication of the velocity update code.

### Output

Output from PIC codes includes the spatial distribution of plasma parameters such as potential, charge density, electron temperature, as well as particle data such as velocities and current densities. In addition, time-dependant output of global aggregate diagnostic data, such as total kinetic and potential energy of the simulation, is helpful in diagnosing code performance.

### Repeat

The loop repeats until some condition is satisfied. Simulations with continuous sources are run until a steady state is achieved. The steady state is characterized by the net number of particles in the simulation domain remaining constant between time steps. In other words, the number of new particles generated by sources is balanced by the number of particles leaving the domain through the boundaries.

1. Birdsall, C.K. and Langdon, A.B.,Plasma Physics via computer simulation, 1985
2. Chen, F.F., Introduction to plasma physics and controlled fuison, 1984
3. Lieberman, M.A., and Lichtenberg, A.J., Principles of plasma discharges and material processing, 2005

That’s it. Next see the simple Particle In Cell code in Matlab post for implementation. Feel free to leave a comment if you have any questions.

### 23 comments to “The Electrostatic Particle In Cell (ES-PIC) Method”

1. Etienne Koen
June 27, 2011 at 6:10 am

Hi, I am a PhD student studying plasma physics and using the PIC method. I am having trouble to calculate the v x B term of the Lorentz equation. Could you please give me information on how to solve the equation using the leapfrog method?

Best Regards,

Etienne.

• Lubos
June 27, 2011 at 6:19 am

Hi Etienne,

Actually the next article I will write for the blog will be exactly on this topic. It should be online by the end of the week. Is that ok? But to get you started, you want to use the method of Boris, or some form of it. It basically consists of half acceleration, full rotation, and then another half acceleration. Do you happen to have a copy of Birdsall’s Plasma Physics via Computer Simulations? (If not you should get it, it’s the bible of PIC). Check on page 61 (Implementation of the vxB rotation).

• Matthias
July 5, 2011 at 1:55 am

Great site here.

Since I have to understand the Boris integrator for my master thesis, I have written a short document about it. Maybe it helps you, Etienne, so you don’t necessarily have to get the book: You can obtain it at http://n.ethz.ch/~toggweim/boris.pdf.
Equation (10) is still in implicit form, but you can solve the 3×3 system easily (for example by hand with the Cramer rule). The method is then fully explicit.

By the way, my master thesis is about an integrator for PIC code that uses a global, adaptive time step to increase efficiency. If you have heard of previous attempts or have good ideas please let me know. The main idea is to evaluate the space charge field less often than external fields (in a multiple time stepping fashion), since the space charge solver is the slowest part in the integration procedure.

Regards,
Matthias

• Lubos
July 6, 2011 at 7:52 am

Thanks for sharing your write up, Matthias. By the way, I am still working on the integrator article, I am hoping to get it done this week. The problem is finding a good example to illustrate the difficulty with standard integration.

• Matthias
July 11, 2011 at 10:05 am

I am really looking forward to see that article!
Interestingly Boris can be derived under a different viewpoint than mentioned in the book, namely as composition of a first order method with its adjoint method, where the first order method is obtained with a splitting approach. (Have almost written this down if you’re interested in more details)
What do you mean with “standard” integration? Have you tried just a single particle in a constant magnetic field that flies a circle? Otherwise, a 2D thing where the magnetic field is a restoring force linear in one coordinate (B(x)=-k*x) should be enough to demonstrate the errors in motion.

• Eduardo Orozco
August 13, 2011 at 2:54 pm

Hi Etienne, I’m PhD student too. I use the PIC method to study the electron beam acceleration by a standing electromagnetic field in static inhomogeneous magnetic field based in cyclotron autoresonance phenomenon. I use the Boris method to solve the motion equation. It work very well. If you have any problem still please let me know.

Eduardo.

2. Matthias
July 23, 2011 at 3:38 pm

The link “MATLAB PIC code” is dead. Can you make the respective content available again?

• July 27, 2011 at 6:13 pm

Will do (and will resume posting articles shortly as well). It’s total crunch time here trying to finish a paper for the upcoming International Electric Propulsion Conference and also submit some proposals…

3. August 10, 2011 at 5:11 am

Great hammer of Thor, that is powerfully helpful!

admin note: this got marked as spam but the comment was too funny not to let it go through…

4. Etienne
August 12, 2011 at 9:18 am

If I want to simulate different populations of plasmas with different densities, do I just allocate a certain ratio of particles to that population to represent the density? Or are there other ways how I can do this?

• September 6, 2011 at 2:06 pm

Hi Etienne, yes, you can easily include multiple particle species with the PIC method. Most of my codes use a fixed specific weight (number of real particle represented by each macroparticle) for each species. This number should be corresponding to the species density. For good statistics, you want to have at least 10-20 particles per cell. This means that you will need to reduce the specific weight for lower density species. Conversely, you will want to increase it for high density species so that you don’t slow down your simulation with too many particles.

5. Sam
March 3, 2012 at 3:23 am

Do you have an upgrading of the Matlab-program in 3D? I have to determine the potential damping of a charge by Debye length. In my Problem, I have a sphere with radius r=Debyelength.
Thank you very much.

• March 3, 2012 at 7:43 am

Hi Sam,

I have a general 3D plasma simulation code but that one is not publicly available. However, it seems your problem would be solved even easier in 1D, by switching to the spherical coordinate system, no?

6. Elnaz
May 5, 2012 at 1:44 am

Hi,I am phd student in laser plasma interaction and I am looking for to find a 2D-Electromagnetic open source PIC code. Please let me know if you have any recommandation in this regards.
All the best

• May 5, 2012 at 6:32 am

Hi Elnaz, I have never used it but I believe that OOPic can be used for EM simulations, see http://ptsg.eecs.berkeley.edu/ for more details.

7. chaudhury
May 15, 2012 at 3:17 am

Hi Lubos, Do you know any article where it is discussed in details the issue of determining the minimum number of particles per cell and the related accuracy issues. I understand this number will vary depending on the nature of the problem and the plasma conditions, and at times we may need to use hundreds of particles per cell. Can you give some insight on this issue. Thanks in advance.

• May 15, 2012 at 5:50 pm

Hi Chaudhury, I don’t know of any paper off the top of my head but the number I have heard is 20 particles per cell. But as always, this will depend on your setup. The particle count gives you two things. First, it decreases the smallest charge unit that can be scattered to the grid. This is especially important if you have a high variation in density, such as when you are doing plume analysis. A higher number will then reduce the gradients along the beam edge. Secondly, the higher number allows you to better resolve the velocity distribution function. This is especially true if you use fixed particle weight (as is for instance the case in my old 3D code Draco). Since we typically load Maxwellian distribution, the larger number gives you a better hope of capturing the high energy electrons that are responsible for much ionization in discharges. But in a sense these particles-per-cell numbers are somewhat meaningless. We typically study electric propulsion plumes, and it is not uncommon to have few thousand particles per cell in the high density region and only a few in the wake. Mesh stretching and variable particle weights can help but these are sometimes tricky to implement correctly.

8. Jimmy
December 12, 2012 at 5:39 pm

Hi Lubos,
I was wondering if the PIC code can be used for description of arcing between electrodes. The arcing may be caused by lightning.

• December 16, 2012 at 9:16 am

Yes in theory, however, how practical it will be will strongly depend on the arc density. PIC is typically more suitable for rarefied flows. This is because as the density increases, the computational requirements for the PIC method grow drastically. Also, the main advantage of kinetic methods is that they don’t assume the Maxwellian distribution function. Rarefied flows may deviate from the Maxwellian since there may not be enough particle interactions to thermalize the flow. In dense flows, thermalization is typically quite rapid, so you may be better off using a fluid based approach, such as magnetohydrodynamics.

9. maha
February 28, 2013 at 3:46 am

Thanh you for the great blog! I have a question, how can we introduce the pressure in the code?

• April 8, 2013 at 4:21 am

This is handled by adding collisions, either with MCC or DSMC methods. With DSMC, you simply add more particles in one area than other and since there is more of them, they will also collide more frequently. The random scattering will then result in the particles diffusing out until you reach an equilibrium situation.

10. Ales Necas
April 7, 2013 at 2:21 pm

Dear Lubos,

Thank you for the fun PIC code. Quick question. Maybe I misunderstood the problem you are solving, but should not the density be zero inside the rectangular object.

Thank you,

Ales

• April 8, 2013 at 4:01 am

Hi Ales, it should be and is, but it may not look like it because of the coarse mesh.