The

**Particle-in-Cell** method refers to a technique used to solve a certain class of partial differential equations. In this method, individual particles (or fluid elements) in a

LagrangianIn fluid dynamics and finite-deformation plasticity the Lagrangian specification of the flow field is a way of looking at fluid motion where the observer follows an individual fluid parcel as it moves through space and time. Plotting the position of an individual parcel through time gives the...

frame are tracked in continuous

phase spaceIn mathematics and physics, a phase space, introduced by Willard Gibbs in 1901, is a space in which all possible states of a system are represented, with each possible state of the system corresponding to one unique point in the phase space...

, whereas moments of the distribution such as densities and currents are computed simultaneously on Eulerian (stationary) mesh points.

PIC methods were already in use as early as 1955

,

even before the first

FortranFortran is a general-purpose, procedural, imperative programming language that is especially suited to numeric computation and scientific computing...

compilers were available. The method gained popularity for plasma simulation in the late 1950s and early 1960s by

BunemanOscar Buneman made advances in science, engineering, and mathematics. Buneman was a pioneer of computational plasma physics and plasma simulation....

,

DawsonJohn M. Dawson was an American computational physicist. He received the Aneesur Rahman prize in computational physics. The Rahman prize is the highest honor given by the American Physical Society for work in computational physics. He is credited as the inventor of the field of Plasma acceleration...

, Hockney, Birdsall, Morse and others. In

plasma physicsIn physics and chemistry, plasma is a state of matter similar to gas in which a certain portion of the particles are ionized. Heating a gas may ionize its molecules or atoms , thus turning it into a plasma, which contains charged particles: positive ions and negative electrons or ions...

applications, the method amounts to following the trajectories of charged particles in self-consistent electromagnetic (or electrostatic) fields computed on a fixed mesh.

## Technical aspects

For many types of problems, the PIC method is relatively intuitive and straightforward to implement. This probably accounts for much of its success, particularly for plasma simulation, for which the method typically includes the following procedures:

- Integration of the equations of motion.
- Interpolation of charge and current source terms to the field mesh.
- Computation of the fields on mesh points.
- Interpolation of the fields from the mesh to the particle locations.

Models which include interactions of particles only through the average fields are called

**PM** (particle-mesh). Those which include direct binary interactions are

**PP** (particle-particle). Models with both types of interactions are called

**PP-PM** or

**P**^{3}M.

Since the early days, it has been recognized that the PIC method is susceptible to error from so-called

*discrete particle noise*.

This error is statistical in nature, and today it remains less-well understood than for traditional fixed-grid methods, such as

EulerianNumerical partial differential equations is the branch of numerical analysis that studies the numerical solution of partial differential equations .Numerical techniques for solving PDEs include the following:...

or

semi-LagrangianThe Semi-Lagrangian scheme is a numerical method that is widely used in Numerical Weather Prediction models for the integration of the equations governing atmospheric motion...

schemes.

## Basics of the PIC plasma simulation technique

Inside the plasma research community, systems of different species (electrons, ions, neutrals, molecules, dust particles, etc.) are investigated. The set of equations associated with PIC codes are therefore the

Lorentz forceIn physics, the Lorentz force is the force on a point charge due to electromagnetic fields. It is given by the following equation in terms of the electric and magnetic fields:...

as the equation of motion, solved in the so-called

*pusher* or

*particle mover* of the code, and

Maxwell's equationsMaxwell's equations are a set of partial differential equations that, together with the Lorentz force law, form the foundation of classical electrodynamics, classical optics, and electric circuits. These fields in turn underlie modern electrical and communications technologies.Maxwell's equations...

determining the

electricIn physics, an electric field surrounds electrically charged particles and time-varying magnetic fields. The electric field depicts the force exerted on other electrically charged objects by the electrically charged particle the field is surrounding...

and

magneticA magnetic field is a mathematical description of the magnetic influence of electric currents and magnetic materials. The magnetic field at any given point is specified by both a direction and a magnitude ; as such it is a vector field.Technically, a magnetic field is a pseudo vector;...

fields, calculated in the

*(field) solver*.

#### Super-particles

The real systems studied are often extremely large in terms of the number of particles they contain. In order to make simulations efficient or at all possible, so-called

*super-particles* are used. A super-particle is a computational particle that represents many real particles; it may be millions of electrons or ions in the case of a plasma simulation, or, for instance, a vortex element in a fluid simulation. It is allowed to rescale the number of particles, because the

Lorentz forceIn physics, the Lorentz force is the force on a point charge due to electromagnetic fields. It is given by the following equation in terms of the electric and magnetic fields:...

depends only on the charge to mass ratio, so a super-particle will follow the same trajectory as a real particle would.

The number of real particles corresponding to a super-particle must be chosen such that sufficient statistics can be collected on the particle motion. If there is a significant difference between the density of different species in the system (between ions and neutrals, for instance), separate real to super-particle ratios can be used for them.

#### The particle mover

Even with super-particles, the number of simulated particles is usually very large (> 10

^{5}), and often the particle mover is the most time consuming part of PIC, since it has to be done for each particle separately. Thus, the pusher is required to be of high accuracy and speed and much effort is spent on optimizing the different schemes.

The schemes used for the particle mover can be split into two categories, implicit and explicit solvers. While implicit solvers calculate the particle velocity from the already updated fields, explicit solvers use only the old force from the previous time step, and are therefore simpler and faster, but require a smaller time step. Two frequently used schemes are the leapfrog method

, which is an implicit solver, and the

*Boris scheme*
, which is explicit.

For plasma applications, the leapfrog method takes the following form:

where the subscript

refers to "old" quantities from the previous time step,

to updated quantities from the next time step (i.e.

), and velocities are calculated in-between the usual time steps

.

In comparison, the equations of the Boris scheme are:

with

and

.

#### The field solver

The most commonly used methods for solving Maxwell's equations (or more generally,

partial differential equationIn mathematics, partial differential equations are a type of differential equation, i.e., a relation involving an unknown function of several independent variables and their partial derivatives with respect to those variables...

s (PDE)) belong to one of the following three categories:

- Finite difference method
In mathematics, finite-difference methods are numerical methods for approximating the solutions to differential equations using finite difference equations to approximate derivatives.- Derivation from Taylor's polynomial :...

s (FDM)
- Finite element method
The finite element method is a numerical technique for finding approximate solutions of partial differential equations as well as integral equations...

s (FEM)
- Spectral method
Spectral methods are a class of techniques used in applied mathematics and scientific computing to numerically solve certain Dynamical Systems, often involving the use of the Fast Fourier Transform. Where applicable, spectral methods have excellent error properties, with the so called "exponential...

s

With the FDM, the continuous domain is replaced with a discrete grid of points, on which the

electricIn physics, an electric field surrounds electrically charged particles and time-varying magnetic fields. The electric field depicts the force exerted on other electrically charged objects by the electrically charged particle the field is surrounding...

and

magneticA magnetic field is a mathematical description of the magnetic influence of electric currents and magnetic materials. The magnetic field at any given point is specified by both a direction and a magnitude ; as such it is a vector field.Technically, a magnetic field is a pseudo vector;...

fields are calculated. Derivatives are then approximated with differences between neighboring grid-point values and thus PDEs are turned into algebraic equations.

Using FEM, the continuous domain is divided into a discrete mesh of elements. The PDEs are treated as an

eigenvalue problemThe eigenvectors of a square matrix are the non-zero vectors that, after being multiplied by the matrix, remain parallel to the original vector. For each eigenvector, the corresponding eigenvalue is the factor by which the eigenvector is scaled when multiplied by the matrix...

and initially a trial solution is calculated using

basis functionIn mathematics, a basis function is an element of a particular basis for a function space. Every continuous function in the function space can be represented as a linear combination of basis functions, just as every vector in a vector space can be represented as a linear combination of basis...

s that are localized in each element. The final solution is then obtained by optimization until the required accuracy is reached.

Also spectral methods, such as the

fast Fourier transformA fast Fourier transform is an efficient algorithm to compute the discrete Fourier transform and its inverse. "The FFT has been called the most important numerical algorithm of our lifetime ." There are many distinct FFT algorithms involving a wide range of mathematics, from simple...

(FFT), transform the PDEs into an eigenvalue problem, but this time the basis functions are high order and defined globally over the whole domain. The domain itself is not discretized in this case, it remains continuous. Again, a trial solution is found by inserting the basis functions into the eigenvalue equation and then optimized to determine the best values of the initial trial parameters.

#### Particle and field weighting

The name "particle-in-cell" originates in the way that plasma macro-quantities (

number densityIn physics, astronomy, and chemistry, number density is an intensive quantity used to describe the degree of concentration of countable objects in the three-dimensional physical space...

,

current densityCurrent density is a measure of the density of flow of a conserved charge. Usually the charge is the electric charge, in which case the associated current density is the electric current per unit area of cross section, but the term current density can also be applied to other conserved...

, etc.) are assigned to simulation particles (i.e., the

*particle weighting*). Particles can be situated anywhere on the continuous domain, but macro-quantities are calculated only on the mesh points, just as the fields are. To obtain the macro-quantities, one assumes that the particles have a given "shape" determined by the shape function

where

is the coordinate of the particle and

the observation point. Perhaps the easiest and most used choice for the shape function is the so-called

*cloud-in-cell* (CIC) scheme, which is a first order (linear) weighting scheme. Whatever the scheme is, the shape function has to satisfy the following conditions:

space isotropy, charge conservation, and increasing accuracy (convergence) for higher-order terms.

The fields obtained from the field solver are determined only on the grid points and can't be used directly in the particle mover to calculate the force acting on particles, but have to be interpolated via the

*field weighting*:

where the subscript

labels the grid point. To ensure that the forces acting on particles are self-consistently obtained, the way of calculating macro-quantities from particle positions on the grid points and interpolating fields from grid points to particle positions has to be consistent, too, since they both appear in

Maxwell's equationsMaxwell's equations are a set of partial differential equations that, together with the Lorentz force law, form the foundation of classical electrodynamics, classical optics, and electric circuits. These fields in turn underlie modern electrical and communications technologies.Maxwell's equations...

. Above all, the field interpolation scheme should conserve

momentumIn classical mechanics, linear momentum or translational momentum is the product of the mass and velocity of an object...

. This can be achieved by choosing the same weighting scheme for particles and fields and by ensuring the appropriate space symmetry (i.e. no self-force and fulfilling the

action-reaction lawNewton's laws of motion are three physical laws that form the basis for classical mechanics. They describe the relationship between the forces acting on a body and its motion due to those forces...

) of the field solver at the same time.

#### Collisions

As the field solver is required to be free of self-forces, inside a cell the field generated by a particle must decrease with decreasing distance from the particle, and hence inter-particle forces inside the cells are underestimated. This can be balanced with the aid of

Coulomb collisionA Coulomb collision is a binary elastic collision between two charged particles interacting through their own Electric Field. As with any inverse-square law, the resulting trajectories of the colliding particles is a hyperbolic Keplerian orbit...

s between charged particles. Simulating the interaction for every pair of a big system would be computationally too expensive, so several

Monte Carlo methodMonte Carlo methods are a class of computational algorithms that rely on repeated random sampling to compute their results. Monte Carlo methods are often used in computer simulations of physical and mathematical systems...

s have been developed instead. A widely used method is the

*binary collision model*
, in which particles are grouped according to their cell, then these particles are paired randomly, and finally the pairs are collided.

In a real plasma, many other reactions may play a role, ranging from elastic collisions, such as collisions between charged and neutral particles, over inelastic collisions, such as electron-neutral ionization collision, to chemical reactions; each of them requiring separate treatment. Most of the collision models handling charged-neutral collisions use either the

*direct Monte-Carlo* scheme, in which all particles carry information about their collision probability, or the

*null-collision* scheme

, which does not analyze all particles but uses the maximum collision probability for each charged species instead.

#### Accuracy and stability conditions

As in every simulation method, also in PIC, the time step and the grid size must be well chosen, so that the shortest time and length scale phenomena are properly resolved in the problem. In addition, time step and grid size have also an impact on the speed and accuracy of the code.

As a rule of thumb, two important conditions regarding the grid size

and the time step

should be fulfilled in order to ensure the stability of the solution:

which can be derived considering the harmonic oscillations of a one-dimensional unmagnetized plasma. The latter conditions is strictly required but practical considerations related to energy conservation suggest to use a much stricter constraint where the factor 2 is replaced by a number one order of magnitude smaller. The use of

is typical. Not surprisingly, the natural time scale in the plasma is given by the inverse

plasma frequencyPlasma oscillations, also known as "Langmuir waves" , are rapid oscillations of the electron density in conducting media such as plasmas or metals. The oscillations can be described as an instability in the dielectric function of a free electron gas. The frequency only depends weakly on the...

and length scale by the

Debye lengthIn plasma physics, the Debye length , named after the Dutch physicist and physical chemist Peter Debye, is the scale over which mobile charge carriers screen out electric fields in plasmas and other conductors. In other words, the Debye length is the distance over which significant charge...

.

## Applications

Within plasma physics, PIC simulation has been used successfully to study laser-plasma interactions, electron acceleration and ion heating in the auroral

ionosphereThe ionosphere is a part of the upper atmosphere, comprising portions of the mesosphere, thermosphere and exosphere, distinguished because it is ionized by solar radiation. It plays an important part in atmospheric electricity and forms the inner edge of the magnetosphere...

,

magnetohydrodynamicsMagnetohydrodynamics is an academic discipline which studies the dynamics of electrically conducting fluids. Examples of such fluids include plasmas, liquid metals, and salt water or electrolytes...

,

magnetic reconnectionMagnetic reconnection is a physical process in highly conducting plasmas in which the magnetic topology is rearranged and magnetic energy is converted to kinetic energy, thermal energy, and particle acceleration...

, as well as ion-temperature-gradient and other microinstabilities in

tokamakA tokamak is a device using a magnetic field to confine a plasma in the shape of a torus . Achieving a stable plasma equilibrium requires magnetic field lines that move around the torus in a helical shape...

s, furthermore

vacuum dischargesA vacuum arc can arise when the surfaces of metal electrodes in contact with a good vacuum begin to emit electrons either through heating or via an electric field that is sufficient to cause field electron emission...

, and

dusty plasmaA dusty plasma is a plasma containing nanometer or micrometer-sized particles suspended in it. Dust particles may be charged and the plasma and particles behave as a plasma, following electromagnetic laws for particles up to about 10 nm...

s.

Hybrid models may use the PIC method for the kinetic treatment of some species, while other species (that are Maxwellian) are simulated with a fluid model.

PIC simulations have also been applied outside of plasma physics to problems in

solidSolid mechanics is the branch of mechanics, physics, and mathematics that concerns the behavior of solid matter under external actions . It is part of a broader study known as continuum mechanics. One of the most common practical applications of solid mechanics is the Euler-Bernoulli beam equation...

and

fluid mechanicsFluid mechanics is the study of fluids and the forces on them. Fluid mechanics can be divided into fluid statics, the study of fluids at rest; fluid kinematics, the study of fluids in motion; and fluid dynamics, the study of the effect of forces on fluid motion...

.

## External links