Particle-in-cell
Encyclopedia
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 Lagrangian
frame are tracked in continuous 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 Fortran
compilers were available. The method gained popularity for plasma simulation in the late 1950s and early 1960s by Buneman
, Dawson
, Hockney, Birdsall, Morse and others. In plasma physics
applications, the method amounts to following the trajectories of charged particles in self-consistent electromagnetic (or electrostatic) fields computed on a fixed mesh.
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 P3M.
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 Eulerian
or semi-Lagrangian
schemes.
as the equation of motion, solved in the so-called pusher or particle mover of the code, and Maxwell's equations
determining the electric
and magnetic
fields, calculated in the (field) solver.
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 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
.
s (PDE)) belong to one of the following three categories:
With the FDM, the continuous domain is replaced with a discrete grid of points, on which the electric
and magnetic
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 problem
and initially a trial solution is calculated using basis function
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 transform
(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.
, current density
, 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 equations
. Above all, the field interpolation scheme should conserve momentum
. 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 law
) of the field solver at the same time.
s between charged particles. Simulating the interaction for every pair of a big system would be computationally too expensive, so several Monte Carlo method
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.
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 frequency
and length scale by the Debye length
.
, magnetohydrodynamics
, magnetic reconnection
, as well as ion-temperature-gradient and other microinstabilities in tokamak
s, furthermore vacuum discharges
, and dusty plasma
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 solid
and fluid mechanics
.
Lagrangian and Eulerian coordinates
In 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 space
Phase space
In 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 Fortran
Fortran
Fortran 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 Buneman
Oscar Buneman
Oscar Buneman made advances in science, engineering, and mathematics. Buneman was a pioneer of computational plasma physics and plasma simulation....
, Dawson
John M. Dawson
John 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 physics
Plasma (physics)
In 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 P3M.
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 Eulerian
Numerical partial differential equations
Numerical 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-Lagrangian
Semi-Lagrangian scheme
The 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 forceLorentz force
In 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 equations
Maxwell's equations
Maxwell'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 electric
Electric field
In 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 magnetic
Magnetic field
A 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 forceLorentz force
In 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 (> 105), 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 equationPartial differential equation
In 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 methodFinite difference methodIn 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 methodFinite element methodThe finite element method is a numerical technique for finding approximate solutions of partial differential equations as well as integral equations...
s (FEM) - Spectral methodSpectral methodSpectral 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 electric
Electric field
In 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 magnetic
Magnetic field
A 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 problem
Eigenvalue, eigenvector and eigenspace
The 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 function
Basis function
In 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 transform
Fast Fourier transform
A 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 densityNumber density
In 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 density
Current density
Current 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
Maxwell'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 momentum
Momentum
In 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 law
Newton's laws of motion
Newton'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 collisionCoulomb collision
A 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 method
Monte Carlo method
Monte 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 oscillation
Plasma 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 lengthDebye length
In 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 ionosphereIonosphere
The 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...
, magnetohydrodynamics
Magnetohydrodynamics
Magnetohydrodynamics 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 reconnection
Magnetic reconnection
Magnetic 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 tokamak
Tokamak
A 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 discharges
Vacuum arc
A 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 plasma
Dusty plasma
A 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 solid
Solid mechanics
Solid 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 mechanics
Fluid mechanics
Fluid 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
- Open source 3D Particle-In-Cell code for spacecraft plasma interactions (mandatory user registration required).
- Simple Particle-In-Cell code in MATLAB
- Plasma Theory and Simulation Group (Berkeley) Contains links to freely-available software.
- Introduction to PIC codes (Univ. of Texas)
- OpenPIC3D - 3D Hybrid Particle-In-Cell simulation of plasma dynamics