Particle-in-cell
   HOME

TheInfoList



OR:

In
plasma physics Plasma ()πλάσμα
, Henry George Liddell, R ...
, the particle-in-cell (PIC) method refers to a technique used to solve a certain class of
partial differential equations In mathematics, a partial differential equation (PDE) is an equation which imposes relations between the various partial derivatives of a multivariable function. The function is often thought of as an "unknown" to be solved for, similarly to ...
. In this method, individual particles (or fluid elements) in a
Lagrangian Lagrangian may refer to: Mathematics * Lagrangian function, used to solve constrained minimization problems in optimization theory; see Lagrange multiplier ** Lagrangian relaxation, the method of approximating a difficult constrained problem with ...
frame are tracked in continuous
phase space In dynamical system theory, a phase space is a space in which all possible states of a system are represented, with each possible state corresponding to one unique point in the phase space. For mechanical systems, the phase space usually ...
, whereas moments of the distribution such as densities and currents are computed simultaneously on Eulerian (stationary)
mesh A mesh is a barrier made of connected strands of metal, fiber, or other flexible or ductile materials. A mesh is similar to a web or a net in that it has many attached or woven strands. Types * A plastic mesh may be extruded, oriented, ex ...
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 Plasma ()πλάσμα
, Henry George Liddell, R ...
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 classical PIC method invented by Buneman, Dawson, Hockney, Birdsall, Morse and others 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 or semi-Lagrangian schemes. Modern geometric PIC algorithms are based on a very different theoretical framework. These algorithms use tools of discrete manifold, interpolating differential forms, and canonical or non-canonical
symplectic integrator In mathematics, a symplectic integrator (SI) is a numerical integration scheme for Hamiltonian systems. Symplectic integrators form the subclass of geometric integrators which, by definition, are canonical transformations. They are widely used in n ...
s to guarantee gauge invariant and conservation of charge, energy-momentum, and more importantly the infinitely dimensional symplectic structure of the particle-field system. These desired features are attributed to the fact that geometric PIC algorithms are built on the more fundamental field-theoretical framework and are directly linked to the perfect form, i.e., the variational principle of physics.


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 force In physics (specifically in electromagnetism) the Lorentz force (or electromagnetic force) is the combination of electric and magnetic force on a point charge due to electromagnetic fields. A particle of charge moving with a velocity in an elect ...
as the equation of motion, solved in the so-called ''pusher'' or ''particle mover'' of the code, and
Maxwell's equations Maxwell's equations, or Maxwell–Heaviside equations, are a set of coupled partial differential equations that, together with the Lorentz force law, form the foundation of classical electromagnetism, classical optics, and electric circuits. ...
determining the
electric Electricity is the set of physical phenomena associated with the presence and motion of matter that has a property of electric charge. Electricity is related to magnetism, both being part of the phenomenon of electromagnetism, as described by ...
and
magnetic Magnetism is the class of physical attributes that are mediated by a magnetic field, which refers to the capacity to induce attractive and repulsive phenomena in other entities. Electric currents and the magnetic moments of elementary particle ...
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 (or ''macroparticle'') 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 acceleration from the
Lorentz force In physics (specifically in electromagnetism) the Lorentz force (or electromagnetic force) is the combination of electric and magnetic force on a point charge due to electromagnetic fields. A particle of charge moving with a velocity in an elect ...
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 (e.g. implicit Euler scheme) 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. In PIC simulation the leapfrog method is used, a second-order explicit method. Also the ''Boris algorithm'' is used which cancel out the magnetic field in the Newton-Lorentz equation. For plasma applications, the leapfrog method takes the following form: :\frac = \mathbf_, :\frac = \frac \left( \mathbf_k + \frac \times \mathbf_ \right), where the subscript k refers to "old" quantities from the previous time step, k+1 to updated quantities from the next time step (i.e. t_ = t_k + \Delta t), and velocities are calculated in-between the usual time steps t_k. The equations of the Boris scheme which are substitute in the above equations are: :\mathbf_ = \mathbf_ + \mathbf_, :\mathbf_ = \mathbf' + q' \mathbf_k, with :\mathbf' = \mathbf + (\mathbf + (\mathbf \times \mathbf)) \times \mathbf, :\mathbf = \mathbf_ + q' \mathbf_k, :\mathbf = q' \mathbf_k, :\mathbf = 2 \mathbf/(1 + h^2) and q' = \Delta t \times (q/2m). Because of its excellent long term accuracy, the Boris algorithm is the de facto standard for advancing a charged particle. It was realized that the excellent long term accuracy of nonrelativistic Boris algorithm is due to the fact it conserves phase space volume, even though it is not symplectic. The global bound on energy error typically associated with symplectic algorithms still holds for the Boris algorithm, making it an effective algorithm for the multi-scale dynamics of plasmas. It has also been shown that one can improve on the relativistic Boris push to make it both volume preserving and have a constant-velocity solution in crossed E and B fields.


The field solver

The most commonly used methods for solving Maxwell's equations (or more generally,
partial differential equation In mathematics, a partial differential equation (PDE) is an equation which imposes relations between the various partial derivatives of a Multivariable calculus, multivariable function. The function is often thought of as an "unknown" to be sol ...
s (PDE)) belong to one of the following three categories: *
Finite difference method In numerical analysis, finite-difference methods (FDM) are a class of numerical techniques for solving differential equations by approximating derivatives with finite differences. Both the spatial domain and time interval (if applicable) are di ...
s (FDM) *
Finite element method The finite element method (FEM) is a popular method for numerically solving differential equations arising in engineering and mathematical modeling. Typical problem areas of interest include the traditional fields of structural analysis, heat ...
s (FEM) *
Spectral method Spectral methods are a class of techniques used in applied mathematics and scientific computing to numerically solve certain differential equations. The idea is to write the solution of the differential equation as a sum of certain " basis functio ...
s With the FDM, the continuous domain is replaced with a discrete grid of points, on which the
electric Electricity is the set of physical phenomena associated with the presence and motion of matter that has a property of electric charge. Electricity is related to magnetism, both being part of the phenomenon of electromagnetism, as described by ...
and
magnetic Magnetism is the class of physical attributes that are mediated by a magnetic field, which refers to the capacity to induce attractive and repulsive phenomena in other entities. Electric currents and the magnetic moments of elementary particle ...
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 In mathematics, a basis function is an element of a particular basis for a function space. Every 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 represen ...
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 A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in th ...
(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 density The number density (symbol: ''n'' or ''ρ''N) is an intensive quantity used to describe the degree of concentration of countable objects (particles, molecules, phonons, cells, galaxies, etc.) in physical space: three-dimensional volumetric number ...
,
current density In electromagnetism, current density is the amount of charge per unit time that flows through a unit area of a chosen cross section. The current density vector is defined as a vector whose magnitude is the electric current per cross-sectional ar ...
, 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 :S(\mathbf-\mathbf), where \mathbf is the coordinate of the particle and \mathbf 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'': :\mathbf(\mathbf) = \sum_\mathbf_i S(\mathbf_i-\mathbf), where the subscript i 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 Maxwell's equations, or Maxwell–Heaviside equations, are a set of coupled partial differential equations that, together with the Lorentz force law, form the foundation of classical electromagnetism, classical optics, and electric circuits. ...
. Above all, the field interpolation scheme should conserve
momentum In Newtonian mechanics, momentum (more specifically linear momentum or translational momentum) is the product of the mass and velocity of an object. It is a vector quantity, possessing a magnitude and a direction. If is an object's mass an ...
. 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


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 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. Th ...
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 methods, or Monte Carlo experiments, are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be determi ...
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 time and length scale phenomena of interest are properly resolved in the problem. In addition, time step and grid size affect the speed and accuracy of the code. For an electrostatic plasma simulation using an explicit time integration scheme (e.g. leapfrog, which is most commonly used), two important conditions regarding the grid size \Delta x and the time step \Delta t should be fulfilled in order to ensure the stability of the solution: :\Delta x < 3.4 \lambda_D, :\Delta t \leq 2 \omega_^, 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 \Delta t \leq 0.1 \omega_^, is typical. Not surprisingly, the natural time scale in the plasma is given by the inverse
plasma frequency Plasma oscillations, also known as Langmuir waves (after Irving Langmuir), are rapid oscillations of the electron density in conducting media such as plasmas or metals in the ultraviolet region. The oscillations can be described as an instability i ...
\omega_^ and length scale by the
Debye length In plasmas and electrolytes, the Debye length \lambda_ (also called Debye radius), is a measure of a charge carrier's net electrostatic effect in a solution and how far its electrostatic effect persists. With each Debye length the charges are in ...
\lambda_D. For an explicit electromagnetic plasma simulation, the time step must also satisfy the CFL condition: :\Delta t < \Delta x / c , where \Delta x \sim \lambda_D, and c is the speed of light.


Applications

Within plasma physics, PIC simulation has been used successfully to study laser-plasma interactions, electron acceleration and ion heating in the auroral
ionosphere The ionosphere () is the ionized part of the upper atmosphere of Earth, from about to above sea level, a region that includes the thermosphere and parts of the mesosphere and exosphere. The ionosphere is ionized by solar radiation. It plays an ...
,
magnetohydrodynamics Magnetohydrodynamics (MHD; also called magneto-fluid dynamics or hydro­magnetics) is the study of the magnetic properties and behaviour of electrically conducting fluids. Examples of such magneto­fluids include plasmas, liquid metals, ...
,
magnetic reconnection Magnetic reconnection is a physical process occurring in highly conducting plasmas in which the magnetic topology is rearranged and magnetic energy is converted to kinetic energy, thermal energy, and particle acceleration. Magnetic reconnectio ...
, as well as ion-temperature-gradient and other microinstabilities in
tokamak A tokamak (; russian: токамáк; otk, 𐱃𐰸𐰢𐰴, Toḳamaḳ) is a device which uses a powerful magnetic field to confine plasma in the shape of a torus. The tokamak is one of several types of magnetic confinement devices being d ...
s, furthermore vacuum discharges, and
dusty plasma A dusty plasma is a plasma containing micrometer (10−6) to nanometer (10−9) sized particles suspended in it. Dust particles are charged and the plasma and particles behave as a plasma. Dust particles may form larger particles resulting in "gra ...
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 is one of the State of matter#Four fundamental states, four fundamental states of matter (the others being liquid, gas, and Plasma (physics), plasma). The molecules in a solid are closely packed together and contain the least amount o ...
and
fluid mechanics Fluid mechanics is the branch of physics concerned with the mechanics of fluids ( liquids, gases, and plasmas) and the forces on them. It has applications in a wide range of disciplines, including mechanical, aerospace, civil, chemical and bio ...
.


Electromagnetic particle-in-cell computational applications


See also

*
Plasma modeling Plasma modeling refers to solving equations of motion that describe the state of a plasma. It is generally coupled with Maxwell's equations for electromagnetic fields or Poisson's equation for electrostatic fields. There are several main types of pl ...
* Multiphase particle-in-cell method


References


Bibliography

* *


External links


Particle-In-Cell and Kinetic Simulation Software Center (PICKSC), UCLA.

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.


open-pic - 3D Hybrid Particle-In-Cell simulation of plasma dynamics
{{Numerical PDE Numerical differential equations Computational fluid dynamics Mathematical modeling Computational electromagnetics