Particleincell
The particleincell (PIC) 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,[1] 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 selfconsistent electromagnetic (or electrostatic) fields computed on a fixed mesh. [2]
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 (particlemesh). Those which include direct binary interactions are PP (particleparticle). Models with both types of interactions are called PPPM or P^{3}M.
Since the early days, it has been recognized that the PIC method is susceptible to error from socalled discrete particle noise. [3] This error is statistical in nature, and today it remains lesswell understood than for traditional fixedgrid methods, such as Eulerian or semiLagrangian 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 noncanonical symplectic integrators to guarantee gauge invariant and conservation of charge, energymomentum, and more importantly the infinitely dimensional symplectic structure of the particlefield system. [4] [5] These desired features are attributed to the fact that geometric PIC algorithms are built on the more fundamental fieldtheoretical 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 as the equation of motion, solved in the socalled pusher or particle mover of the code, and Maxwell's equations determining the electric and magnetic fields, calculated in the (field) solver.
Superparticles
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, socalled superparticles are used. A superparticle (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 Lorentz force depends only on the chargetomass ratio, so a superparticle will follow the same trajectory as a real particle would.
The number of real particles corresponding to a superparticle 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 superparticle ratios can be used for them.
The particle mover
Even with superparticles, 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 (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 secondorder explicit method. [6] Also the Boris algorithm is used which cancel out the magnetic field in the NewtonLorentz equation. [7] [8].
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 inbetween the usual time steps .
The equations of the Boris scheme which are substitute in the above equations are:
with
and .
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 multiscale dynamics of plasmas. It has also been shown [9] that one can improve on the relativistic Boris push to make it both volume preserving and have a constantvelocity 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 equations (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 gridpoint 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 functions 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.
Particle and field weighting
The name "particleincell" originates in the way that plasma macroquantities (number density, current density, etc.) are assigned to simulation particles (i.e., the particle weighting). Particles can be situated anywhere on the continuous domain, but macroquantities are calculated only on the mesh points, just as the fields are. To obtain the macroquantities, 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 socalled cloudincell (CIC) scheme, which is a first order (linear) weighting scheme. Whatever the scheme is, the shape function has to satisfy the following conditions: [10] space isotropy, charge conservation, and increasing accuracy (convergence) for higherorder 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 selfconsistently obtained, the way of calculating macroquantities 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 selfforce and fulfilling the actionreaction law) of the field solver at the same time[10]
Collisions
As the field solver is required to be free of selfforces, inside a cell the field generated by a particle must decrease with decreasing distance from the particle, and hence interparticle forces inside the cells are underestimated. This can be balanced with the aid of Coulomb collisions between charged particles. Simulating the interaction for every pair of a big system would be computationally too expensive, so several Monte Carlo methods have been developed instead. A widely used method is the binary collision model,[11] 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 electronneutral ionization collision, to chemical reactions; each of them requiring separate treatment. Most of the collision models handling chargedneutral collisions use either the direct MonteCarlo scheme, in which all particles carry information about their collision probability, or the nullcollision scheme,[12][13] 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 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 onedimensional 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.[10][14] Not surprisingly, the natural time scale in the plasma is given by the inverse plasma frequency and length scale by the Debye length .
For an explicit electromagnetic plasma simulation, the time step must also satisfy the CFL condition:
where , and is the speed of light.
Applications
Within plasma physics, PIC simulation has been used successfully to study laserplasma interactions, electron acceleration and ion heating in the auroral ionosphere, magnetohydrodynamics, magnetic reconnection, as well as iontemperaturegradient and other microinstabilities in tokamaks, furthermore vacuum discharges, and dusty plasmas.
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. [15] [16]
Electromagnetic particleincell computational applications
Computational application  Web site  License  Availability  Canonical Reference 

SHARP  [17]  Proprietary  doi:10.3847/15384357/aa6d13  
ALaDyn  [18]  GPLv3+  Open Repo:[19]  doi:10.5281/zenodo.49553 
EPOCH  [20]  GPL  Open to academic users but signup required :[21]  doi:10.1088/07413335/57/11/113001 
FBPIC  [22]  3ClauseBSDLBNL  Open Repo:[23]  doi:10.1016/j.cpc.2016.02.007 
LSP  [24]  Proprietary  Available from ATK  doi:10.1016/S01689002(01)000249 
MAGIC  [25]  Proprietary  Available from ATK  doi:10.1016/00104655(95)00010D 
OSIRIS  [26]  Proprietary  Closed (Collaborators with MoU)  doi:10.1007/3540477896_36 
PICCANTE  [27]  GPLv3+  Open Repo:[28]  doi:10.5281/zenodo.48703 
PICLas  [29]  Proprietary  Available from Institute of Space Systems and Institute of Aerodynamics and Gas Dynamics at the University of Stuttgart  doi:10.1016/j.crme.2014.07.005 
PIConGPU  [30]  GPLv3+  Open Repo:[31]  doi:10.1145/2503210.2504564 
SMILEI  [32]  CeCILLB  Open Repo:[33]  doi:10.1016/j.cpc.2017.09.024 
The Virtual Laser Plasma Library  [34]  Proprietary  Unknown  doi:10.1017/S0022377899007515 
VizGrain  [35]  Proprietary  Commercially available from Esgee Technologies Inc.  
VPIC  [36]  3ClauseBSD  Open Repo:[37]  doi:10.1063/1.2840133 
VSim (Vorpal)  [38]  Proprietary  Available from TechX Corporation  doi:10.1016/j.jcp.2003.11.004 
Warp  [39]  3ClauseBSDLBNL  Open Repo:[40]  doi:10.1063/1.860024 
WarpX  [41]  3ClauseBSDLBNL  Open Repo:[42]  doi:10.1016/j.nima.2018.01.035 
ZPIC  [43]  AGPLv3+  Open Repo:[44] 
References

F.H. Harlow (1955). "A Machine Calculation Method for Hydrodynamic Problems". Los Alamos Scientific Laboratory report LAMS1956. Cite journal requires
journal=
(help)  Dawson, J.M. (1983). "Particle simulation of plasmas". Reviews of Modern Physics. 55 (2): 403–447. Bibcode:1983RvMP...55..403D. doi:10.1103/RevModPhys.55.403.
 Hideo Okuda (1972). "Nonphysical noises and instabilities in plasma simulation due to a spatial grid". Journal of Computational Physics. 10 (3): 475–486. Bibcode:1972JCoPh..10..475O. doi:10.1016/00219991(72)900484.
 Qin, H.; Liu, J.; Xiao, J.; et al. (2016). "Canonical symplectic particleincell method for longterm largescale simulations of the VlasovMaxwell system". Nuclear Fusion. 56 (1): 014001. arXiv:1503.08334. Bibcode:2016NucFu..56a4001Q. doi:10.1088/00295515/56/1/014001.
 Xiao, J.; Qin, H.; Liu, J.; et al. (2015). "Explicit highorder noncanonical symplectic particleincell algorithms for VlasovMaxwell systems". Physics of Plasmas. 22 (11): 12504. arXiv:1510.06972. Bibcode:2015PhPl...22k2504X. doi:10.1063/1.4935904.
 Birdsall, Charles K.; A. Bruce Langdon (1985). Plasma Physics via Computer Simulation. McGrawHill. ISBN 0070053715.
 Boris, J.P. (November 1970). "Relativistic plasma simulationoptimization of a hybrid code". Proceedings of the 4th Conference on Numerical Simulation of Plasmas. Naval Res. Lab., Washington, D.C. pp. 3–67.
 Qin, H.; et al. (2013). "Why is Boris algorithm so good?" (PDF). Physics of Plasmas. 20 (5): 084503. Bibcode:2013PhPl...20h4503Q. doi:10.1063/1.4818428.
 Higuera, Adam V.; John R. Cary (2017). "Structurepreserving secondorder integration of relativistic charged particle trajectories in electromagnetic fields". Physics of Plasmas. 24 (5): 052104. Bibcode:2004JCoPh.196..448N. doi:10.1016/j.jcp.2003.11.004.
 Tskhakaya, David (2008). "Chapter 6: The ParticleinCell Method". In Fehske, Holger; Schneider, Ralf; Weiße, Alexander (eds.). Computational ManyParticle Physics. Lecture Notes in Physics 739. 739. Springer, Berlin Heidelberg. doi:10.1007/9783540746867. ISBN 9783540746850.
 Takizuka, Tomonor; Abe, Hirotada (1977). "A binary collision model for plasma simulation with a particle code". Journal of Computational Physics. 25 (3): 205–219. Bibcode:1977JCoPh..25..205T. doi:10.1016/00219991(77)900997.
 Birdsall, C.K. (1991). "Particleincell chargedparticle simulations, plus Monte Carlo collisions with neutral atoms, PICMCC". IEEE Transactions on Plasma Science. 19 (2): 65–85. Bibcode:1991ITPS...19...65B. doi:10.1109/27.106800. ISSN 00933813.
 Vahedi, V.; Surendra, M. (1995). "A Monte Carlo collision model for the particleincell method: applications to argon and oxygen discharges". Computer Physics Communications. 87 (1–2): 179–198. Bibcode:1995CoPhC..87..179V. doi:10.1016/00104655(94)00171W. ISSN 00104655.
 Tskhakaya, D.; Matyash, K.; Schneider, R.; Taccogna, F. (2007). "The ParticleInCell Method". Contributions to Plasma Physics. 47 (8–9): 563–594. Bibcode:2007CoPP...47..563T. doi:10.1002/ctpp.200710072.
 Liu, G.R.; M.B. Liu (2003). Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific. ISBN 9812384561.
 Byrne, F. N.; Ellison, M. A.; Reid, J. H. (1964). "The particleincell computing method for fluid dynamics". Methods Comput. Phys. 3 (3): 319–343. Bibcode:1964SSRv....3..319B. doi:10.1007/BF00230516.
 Shalaby, Mohamad; Broderick, Avery E.; Chang, Philip; Pfrommer, Christoph; Lamberts, Astrid; Puchwein, Ewald (23 May 2017). "SHARP: A Spatially Higherorder, Relativistic ParticleinCell Code". The Astrophysical Journal. 841 (1): 52. arXiv:1702.04732. Bibcode:2017ApJ...841...52S. doi:10.3847/15384357/aa6d13.
 "ALaDyn". ALaDyn. Retrieved 1 December 2017.
 "ALaDyn: A HighAccuracy PIC Code for the MaxwellVlasov Equations". GitHub.com. 18 November 2017. Retrieved 1 December 2017.
 "Codes". Ccpp.ac.uk. Retrieved 1 December 2017.
 "Sign in". GitLab. Retrieved 1 December 2017.
 "FBPIC documentation — FBPIC 0.6.0 documentation". fbpic.github.io. Retrieved 1 December 2017.
 "fbpic: Spectral, quasi3D ParticleInCell code, for CPU and GPU". GitHub.com. 8 November 2017. Retrieved 1 December 2017.
 "Orbital ATK". Mrcwdc.com. Retrieved 1 December 2017.
 "Orbital ATK". Mrcwdc.com. Retrieved 1 December 2017.
 "OSIRIS  PICKSC". Picksc.idre.ucla.edu. Retrieved 1 December 2017.
 "Piccante". Aladyn.github.io. Retrieved 1 December 2017.
 "piccante: a spicy massively parallel fullyrelativistic electromagnetic 3D particleincell code". GitHub.com. 14 November 2017. Retrieved 1 December 2017.
 "PICLas".
 "PIConGPU  ParticleinCell Simulations for the Exascale Era  HelmholtzZentrum DresdenRossendorf, HZDR". picongpu.hzdr.de. Retrieved 1 December 2017.
 "ComputationalRadiationPhysics / PIConGPU — GitHub". GitHub.com. 28 November 2017. Retrieved 1 December 2017.
 "Smilei — A ParticleInCell code for plasma simulation". Maisondelasimulation.fr. Retrieved 1 December 2017.
 "SmileiPIC / Smilei — GitHub". GitHub.com. 29 October 2019. Retrieved 29 October 2019.
 Dreher, Matthias. "Relativistic Laser Plasma". 2.mpq.mpg.de. Retrieved 1 December 2017.
 "VizGrain". esgeetech.com. Retrieved 1 December 2017.
 "VPIC". github.com. Retrieved 1 July 2019.
 "LANL / VPIC — GitHub". github.com. Retrieved 29 October 2019.
 "TechX  VSim". Txcorp.com. Retrieved 1 December 2017.
 "Warp". warp.lbl.gov. Retrieved 1 December 2017.
 "berkeleylab / Warp — Bitbucket". bitbucket.org. Retrieved 1 December 2017.
 "WarpX Documentation". ecpwarpx.github.io. Retrieved 29 October 2019.
 "ECPWarpX / WarpX — GitHub". GitHub.org. Retrieved 29 October 2019.
 "Educational ParticleInCell code suite". picksc.idre.ucla.edu. Retrieved 29 October 2019.
 "ricardofonseca / ZPIC — GitHub". GitHub.org. Retrieved 29 October 2019.
Bibliography
 Birdsall, Charles K.; A. Bruce Langdon (1985). Plasma Physics via Computer Simulation. McGrawHill. ISBN 0070053715.
 Hockney, Roger W.; James W. Eastwood (1988). Computer Simulation Using Particles. CRC Press. ISBN 0852743920.
External links
 ParticleInCell and Kinetic Simulation Software Center (PICKSC), UCLA.
 Open source 3D ParticleInCell code for spacecraft plasma interactions (mandatory user registration required).
 Simple ParticleInCell code in MATLAB
 Plasma Theory and Simulation Group (Berkeley) Contains links to freely available software.
 Introduction to PIC codes (Univ. of Texas)
 openpic  3D Hybrid ParticleInCell simulation of plasma dynamics