# Particle-in-Cell Method

In the *electromagnetic particle-in-cell method* [1, 2],
the electromagnetic fields are solved on a grid, usually using Maxwell’s
equations

given here in natural units (\(\epsilon_0=\mu_0=c=1\)), where \(t\) is time, \(\mathbf{E}\) and \(\mathbf{B}\) are the electric and magnetic field components, and \(\rho\) and \(\mathbf{J}\) are the charge and current densities. The charged particles are advanced in time using the Newton-Lorentz equations of motion

where \(m\), \(q\), \(\mathbf{x}\), \(\mathbf{v}\) and \(\gamma=1/\sqrt{1-v^{2}}\) are respectively the mass, charge, position, velocity and relativistic factor of the particle given in natural units (\(c=1\)). The charge and current densities are interpolated on the grid from the particles’ positions and velocities, while the electric and magnetic field components are interpolated from the grid to the particles’ positions for the velocity update.

## Particle push

A centered finite-difference discretization of the Newton-Lorentz equations of motion is given by

In order to close the system, \(\bar{\mathbf{v}}^{i}\) must be expressed as a function of the other quantities. The two implementations that have become the most popular are presented below.

### Boris relativistic velocity rotation

The solution proposed by Boris [3] is given by

where \(\bar{\gamma}^{i}\) is defined by \(\bar{\gamma}^{i} \equiv (\gamma^{i+1/2}+\gamma^{i-1/2} )/2\).

The system (8, 9) is solved very efficiently following Boris’ method, where the electric field push is decoupled from the magnetic push. Setting \(\mathbf{u}=\gamma\mathbf{v}\), the velocity is updated using the following sequence:

where \(\mathbf{t}=\left(q\Delta t/2m\right)\mathbf{B}^{i}/\bar{\gamma}^{i}\) and where \(\bar{\gamma}^{i}\) can be calculated as \(\bar{\gamma}^{i}=\sqrt{1+(\mathbf{u}^-/c)^2}\).

The Boris implementation is second-order accurate, time-reversible and fast. Its implementation is very widespread and used in the vast majority of PIC codes.

### Vay Lorentz-invariant formulation

It was shown in Vay [4] that the Boris formulation is not Lorentz invariant and can lead to significant errors in the treatment of relativistic dynamics. A Lorentz invariant formulation is obtained by considering the following velocity average

This gives a system that is solvable analytically (see Vay [4] for a detailed derivation), giving the following velocity update:

where

This Lorentz invariant formulation is particularly well suited for the modeling of ultra-relativistic charged particle beams, where the accurate account of the cancellation of the self-generated electric and magnetic fields is essential, as shown in Vay [4].

## Field solve

Various methods are available for solving Maxwell’s equations on a
grid, based on finite-differences, finite-volume, finite-element,
spectral, or other discretization techniques that apply most commonly
on single structured or unstructured meshes and less commonly on multiblock
multiresolution grid structures. In this chapter, we summarize the widespread
second order finite-difference time-domain (FDTD) algorithm, its extension
to non-standard finite-differences as well as the pseudo-spectral
analytical time-domain (PSATD) and pseudo-spectral time-domain (PSTD)
algorithms. Extension to multiresolution (or mesh refinement) PIC
is described in, e.g., Vay *et al.* [5], Vay *et al.* [6].

### Finite-Difference Time-Domain (FDTD)

The most popular algorithm for electromagnetic PIC codes is the Finite-Difference Time-Domain (or FDTD) solver

The differential operator is defined as \(\nabla=D_{x}\mathbf{\hat{x}}+D_{y}\mathbf{\hat{y}}+D_{z}\mathbf{\hat{z}}\) and the finite-difference operators in time and space are defined respectively as

where \(\Delta t\) and \(\Delta x\) are respectively the time step and the grid cell size along \(x\), \(n\) is the time index and \(i\), \(j\) and \(k\) are the spatial indices along \(x\), \(y\) and \(z\) respectively. The difference operators along \(y\) and \(z\) are obtained by circular permutation. The equations in brackets are given for completeness, as they are often not actually solved, thanks to the usage of a so-called charge conserving algorithm, as explained below. As shown in Fig. 28, the quantities are given on a staggered (or “Yee”) grid [7], where the electric field components are located between nodes and the magnetic field components are located in the center of the cell faces. Knowing the current densities at half-integer steps, the electric field components are updated alternately with the magnetic field components at integer and half-integer steps respectively.

### Non-Standard Finite-Difference Time-Domain (NSFDTD)

An implementation of the source-free Maxwell’s wave equations for narrow-band
applications based on non-standard finite-differences (NSFD)
was introduced in Cole [8], Cole [9], and
was adapted for wideband applications in Karkkainen *et al.* [10]. At
the Courant limit for the time step and for a given set of parameters,
the stencil proposed in Karkkainen *et al.* [10] has no numerical dispersion
along the principal axes, provided that the cell size is the same
along each dimension (i.e. cubic cells in 3D). The “Cole-Karkkainen”
(or CK) solver uses the non-standard finite difference formulation
(based on extended stencils) of the Maxwell-Ampere equation and can be
implemented as follows [11]:

Eqs. (19) and (20) are not being solved explicitly but verified via appropriate initial conditions and current deposition procedure. The NSFD differential operator is given by

where

with

Here \(G\) is a sample vector component, while \(\alpha\), \(\beta\) and \(\xi\) are constant scalars satisfying \(\alpha+4\beta+4\xi=1\). As with the FDTD algorithm, the quantities with half-integer are located between the nodes (electric field components) or in the center of the cell faces (magnetic field components). The operators along \(y\) and \(z\), i.e. \(D_{y}\), \(D_{z}\), \(D_{y}^{*}\), \(D_{z}^{*}\), \(S_{y}^{1}\), \(S_{z}^{1}\), \(S_{y}^{2}\), and \(S_{z}^{2}\), are obtained by circular permutation of the indices.

Assuming cubic cells (\(\Delta x=\Delta y=\Delta z\)), the coefficients
given in Karkkainen *et al.* [10] (\(\alpha=7/12\), \(\beta=1/12\) and \(\xi=1/48\))
allow for the Courant condition to be at \(\Delta t=\Delta x\), which
equates to having no numerical dispersion along the principal axes.
The algorithm reduces to the FDTD algorithm with \(\alpha=1\) and \(\beta=\xi=0\).
An extension to non-cubic cells is provided in 3-D by Cowan *et al.* [12] and in 2-D by
Pukhov [13]. An alternative NSFDTD implementation that enables superluminous waves is also
given in Lehe *et al.* [14].

As mentioned above, a key feature of the algorithms based on NSFDTD
is that some implementations [10, 12] enable the time step \(\Delta t=\Delta x\) along one or
more axes and no numerical dispersion along those axes. However, as
shown in Vay *et al.* [11], an instability develops at the Nyquist
wavelength at (or very near) such a timestep. It is also shown in
the same paper that removing the Nyquist component in all the source
terms using a bilinear filter (see description of the filter below)
suppresses this instability.

### Pseudo Spectral Analytical Time Domain (PSATD)

Maxwell’s equations in Fourier space are given by

where \(\tilde{a}\) is the Fourier Transform of the quantity \(a\). As with the real space formulation, provided that the continuity equation \(\partial\tilde{\rho}/\partial t+i\mathbf{k}\cdot\mathbf{\tilde{J}}=0\) is satisfied, then the last two equations will automatically be satisfied at any time if satisfied initially and do not need to be explicitly integrated.

Decomposing the electric field and current between longitudinal and transverse components

gives

with \(\mathbf{\hat{k}}=\mathbf{k}/k\).

If the sources are assumed to be constant over a time interval \(\Delta t\),
the system of equations is solvable analytically and is given by (see Haber *et al.* [15] for the original formulation and Vay *et al.* [16]
for a more detailed derivation):

with \(C=\cos\left(k\Delta t\right)\) and \(S=\sin\left(k\Delta t\right)\).

Combining the transverse and longitudinal components, gives

For fields generated by the source terms without the self-consistent dynamics of the charged particles, this algorithm is free of numerical dispersion and is not subject to a Courant condition. Furthermore, this solution is exact for any time step size subject to the assumption that the current source is constant over that time step.

As shown in Vay *et al.* [16], by expanding the coefficients \(S_{h}\)
and \(C_{h}\) in Taylor series and keeping the leading terms, the PSATD
formulation reduces to the perhaps better known pseudo-spectral time-domain
(PSTD) formulation [17, 18]:

The dispersion relation of the PSTD solver is given by \(\sin(\frac{\omega\Delta t}{2})=\frac{k\Delta t}{2}.\) In contrast to the PSATD solver, the PSTD solver is subject to numerical dispersion for a finite time step and to a Courant condition that is given by \(\Delta t\leq \frac{2}{\pi}\left(\frac{1}{\Delta x^{2}}+\frac{1}{\Delta y^{2}}+\frac{1}{\Delta z^{2}}\right)^{-1/2}\).

The PSATD and PSTD formulations that were just given apply to the field components located at the nodes of the grid. As noted in Ohmura and Okamura [19], they can also be easily recast on a staggered Yee grid by multiplication of the field components by the appropriate phase factors to shift them from the collocated to the staggered locations. The choice between a collocated and a staggered formulation is application-dependent.

Spectral solvers used to be very popular in the years 1970s to early 1990s, before being replaced by finite-difference methods with the advent of parallel supercomputers that favored local methods. However, it was shown recently that standard domain decomposition with Fast Fourier Transforms that are local to each subdomain could be used effectively with PIC spectral methods [16], at the cost of truncation errors in the guard cells that could be neglected. A detailed analysis of the effectiveness of the method with exact evaluation of the magnitude of the effect of the truncation error is given in Vincenti and Vay [20] for stencils of arbitrary order (up-to the infinite “spectral” order).

WarpX also includes a kinetic-fluid hybrid model in which the electric field is calculated using Ohm’s law instead of directly evolving Maxwell’s equations. This approach allows reduced physics simulations to be done with significantly lower spatial and temporal resolution than in the standard, fully kinetic, PIC. Details of this model can be found in the section Kinetic-fluid hybrid model.

## Current deposition

The current densities are deposited on the computational grid from the particle position and velocities, employing splines of various orders [21].

In most applications, it is essential to prevent the accumulation of errors resulting from the violation of the discretized Gauss’ Law. This is accomplished by providing a method for depositing the current from the particles to the grid that preserves the discretized Gauss’ Law, or by providing a mechanism for “divergence cleaning” [1, 22, 23, 24, 25]. For the former, schemes that allow a deposition of the current that is exact when combined with the Yee solver is given in Villasenor and Buneman [26] for linear splines and in Esirkepov [27] for splines of arbitrary order.

The NSFDTD formulations given above and in Vay *et al.* [11], Cowan *et al.* [12], Pukhov [13], Lehe *et al.* [14]
apply to the Maxwell-Faraday
equation, while the discretized Maxwell-Ampere equation uses the FDTD
formulation. Consequently, the charge conserving algorithms developed
for current deposition [26, 27] apply
readily to those NSFDTD-based formulations. More details concerning
those implementations, including the expressions for the numerical
dispersion and Courant condition are given
in Vay *et al.* [11], Cowan *et al.* [12], Pukhov [13], Lehe *et al.* [14].

### Current correction

In the case of the pseudospectral solvers, the current deposition algorithm generally does not satisfy the discretized continuity equation in Fourier space:

In this case, a Boris correction [1] can be applied in \(k\) space in the form

where \(\mathbf{\tilde{E}}_{c}\) is the corrected field. Alternatively, a correction to the current can be applied (with some similarity to the current deposition presented by Morse and Nielson in their potential-based model in Morse and Nielson [28]) using

where \(\mathbf{\tilde{J}}_{c}\) is the corrected current. In this case, the transverse component of the current is left untouched while the longitudinal component is effectively replaced by the one obtained from integration of the continuity equation, ensuring that the corrected current satisfies the continuity equation. The advantage of correcting the current rather than the electric field is that it is more local and thus more compatible with domain decomposition of the fields for parallel computation [16].

### Vay deposition

Alternatively, an exact current deposition can be written for the pseudo-spectral solvers, following the geometrical interpretation of existing methods in real space [26, 27, 28].

The Vay deposition scheme is the generalization of the Esirkepov deposition scheme for the spectral case with arbitrary-order stencils [16]. The current density \(\widehat{\boldsymbol{J}}^{\,n+1/2}\) in Fourier space is computed as \(\widehat{\boldsymbol{J}}^{\,n+1/2} = i \, \widehat{\boldsymbol{D}} / \boldsymbol{k}\) when \(\boldsymbol{k} \neq 0\) and set to zero otherwise. The quantity \(\boldsymbol{D}\) is deposited in real space by averaging the currents over all possible grid paths between the initial position \(\boldsymbol{x}^{\,n}\) and the final position \(\boldsymbol{x}^{\,n+1}\) and is defined as

2D Cartesian geometry:

3D Cartesian geometry:

Here, \(w_i\) represents the weight of the \(i\)-th macro-particle and \(\Gamma\) represents its shape factor. Note that in 2D Cartesian geometry, \(D_y\) is effectively \(J_y\) and does not require additional operations in Fourier space.

## Field gather

In general, the field is gathered from the mesh onto the macroparticles using splines of the same order as for the current deposition \(\mathbf{S}=\left(S_{x},S_{y},S_{z}\right)\). Three variations are considered:

“momentum conserving”: fields are interpolated from the grid nodes to the macroparticles using \(\mathbf{S}=\left(S_{nx},S_{ny},S_{nz}\right)\) for all field components (if the fields are known at staggered positions, they are first interpolated to the nodes on an auxiliary grid),

“energy conserving (or Galerkin)”: fields are interpolated from the staggered Yee grid to the macroparticles using \(\left(S_{nx-1},S_{ny},S_{nz}\right)\) for \(E_{x}\), \(\left(S_{nx},S_{ny-1},S_{nz}\right)\) for \(E_{y}\), \(\left(S_{nx},S_{ny},S_{nz-1}\right)\) for \(E_{z}\), \(\left(S_{nx},S_{ny-1},S_{nz-1}\right)\) for \(B_{x}\), \(\left(S_{nx-1},S_{ny},S_{nz-1}\right)\) for \(B{}_{y}\) and\(\left(S_{nx-1},S_{ny-1},S_{nz}\right)\) for \(B_{z}\) (if the fields are known at the nodes, they are first interpolated to the staggered positions on an auxiliary grid),

“uniform”: fields are interpolated directly form the Yee grid to the macroparticles using \(\mathbf{S}=\left(S_{nx},S_{ny},S_{nz}\right)\) for all field components (if the fields are known at the nodes, they are first interpolated to the staggered positions on an auxiliary grid).

As shown in Birdsall and Langdon [1], Hockney and Eastwood [2], Lewis [29], the momentum and energy conserving schemes conserve momentum and energy respectively at the limit of infinitesimal time steps and generally offer better conservation of the respective quantities for a finite time step. The uniform scheme does not conserve momentum nor energy in the sense defined for the others but is given for completeness, as it has been shown to offer some interesting properties in the modeling of relativistically drifting plasmas [30].

## Filtering

It is common practice to apply digital filtering to the charge or current density in Particle-In-Cell simulations as a complement or an alternative to using higher order splines [1]. A commonly used filter in PIC simulations is the three points filter

where \(\phi^{f}\) is the filtered quantity. This filter is called a bilinear filter when \(\alpha=0.5\). Assuming \(\phi=e^{jkx}\) and \(\phi^{f}=g\left(\alpha,k\right)e^{jkx}\), the filter gain \(g\) is given as a function of the filtering coefficient \(\alpha\) and the wavenumber \(k\) by

The total attenuation \(G\) for \(n\) successive applications of filters of coefficients \(\alpha_{1}\)…\(\alpha_{n}\) is given by

A sharper cutoff in \(k\) space is provided by using \(\alpha_{n}=n-\sum_{i=1}^{n-1}\alpha_{i}\), so that \(G\approx1+O\left(k^{4}\right)\). Such step is called a “compensation” step [1]. For the bilinear filter (\(\alpha=1/2\)), the compensation factor is \(\alpha_{c}=2-1/2=3/2\). For a succession of \(n\) applications of the bilinear factor, it is \(\alpha_{c}=n/2+1\).

It is sometimes necessary to filter on a relatively wide band of wavelength,
necessitating the application of a large number of passes of the bilinear
filter or on the use of filters acting on many points. The former
can become very intensive computationally while the latter is problematic
for parallel computations using domain decomposition, as the footprint
of the filter may eventually surpass the size of subdomains. A workaround
is to use a combination of filters of limited footprint. A solution
based on the combination of three point filters with various strides
was proposed in Vay *et al.* [11] and operates as follows.

The bilinear filter provides complete suppression of the signal at the grid Nyquist wavelength (twice the grid cell size). Suppression of the signal at integer multiples of the Nyquist wavelength can be obtained by using a stride \(s\) in the filter

for which the gain is given by

For a given stride, the gain is given by the gain of the bilinear
filter shifted in k space, with the pole \(g=0\) shifted from the wavelength
\(\lambda=2/\Delta x\) to \(\lambda=2s/\Delta x\), with additional poles,
as given by \(sk\Delta x=\arccos\left(\frac{\alpha}{\alpha-1}\right)\pmod{2\pi}\).
The resulting filter is pass band between the poles, but since the
poles are spread at different integer values in k space, a wide band
low pass filter can be constructed by combining filters using different
strides. As shown in Vay *et al.* [11], the successive application
of 4-passes + compensation of filters with strides 1, 2 and 4 has
a nearly equivalent fall-off in gain as 80 passes + compensation of
a bilinear filter. Yet, the strided filter solution needs only 15
passes of a three-point filter, compared to 81 passes for an equivalent
n-pass bilinear filter, yielding a gain of 5.4 in number of operations
in favor of the combination of filters with stride. The width of the
filter with stride 4 extends only on 9 points, compared to 81 points
for a single pass equivalent filter, hence giving a gain of 9 in compactness
for the stride filters combination in comparison to the single-pass
filter with large stencil, resulting in more favorable scaling with the number
of computational cores for parallel calculations.

C. K. Birdsall and A. B. Langdon. *Plasma Physics Via Computer Simulation*. Adam-Hilger, 1991. ISBN 0 07 005371 5.

R. W. Hockney and J. W. Eastwood. *Computer simulation using particles*. Routledge, 1988. ISBN 0-85274-392-0.

J. P. Boris. Relativistic Plasma Simulation-Optimization of a Hybrid Code. In *Proc. Fourth Conf. Num. Sim. Plasmas*, 3–67. Naval Res. Lab., Wash., D. C., 1970.

J.-L. Vay. Simulation Of Beams Or Plasmas Crossing At Relativistic Velocity. *Physics of Plasmas*, 15(5):56701, May 2008. doi:10.1063/1.2837054.

J.-L. Vay, D. P. Grote, R. H. Cohen, and A. Friedman. Novel methods in the particle-in-cell accelerator code-framework warp. *Computational Science and Discovery*, 5(1):014019 (20 pp.), 2012.

J.-L. Vay, J.-C. Adam, and A. Heron. Asymmetric Pml For The Absorption Of Waves. Application To Mesh Refinement In Electromagnetic Particle-In-Cell Plasma Simulations. *Computer Physics Communications*, 164(1-3):171–177, Dec 2004. doi:10.1016/J.Cpc.2004.06.026.

K. S. Yee. Numerical Solution Of Initial Boundary Value Problems Involving Maxwells Equations In Isotropic Media. *Ieee Transactions On Antennas And Propagation*, Ap14(3):302–307, 1966.

J. B. Cole. A High-Accuracy Realization Of The Yee Algorithm Using Non-Standard Finite Differences. *Ieee Transactions On Microwave Theory And Techniques*, 45(6):991–996, Jun 1997.

J. B. Cole. High-Accuracy Yee Algorithm Based On Nonstandard Finite Differences: New Developments And Verifications. *Ieee Transactions On Antennas And Propagation*, 50(9):1185–1191, Sep 2002. doi:10.1109/Tap.2002.801268.

M. Karkkainen, E. Gjonaj, T. Lau, and T. Weiland. Low-Dispersionwake Field Calculation Tools. In *Proc. Of International Computational Accelerator Physics Conference*, 35–40. Chamonix, France, 2006.

J.-L. Vay, C. G. R. Geddes, E. Cormier-Michel, and D. P. Grote. Numerical Methods For Instability Mitigation In The Modeling Of Laser Wakefield Accelerators In A Lorentz-Boosted Frame. *Journal of Computational Physics*, 230(15):5908–5929, Jul 2011. doi:10.1016/J.Jcp.2011.04.003.

B. M. Cowan, D. L. Bruhwiler, J. R. Cary, E. Cormier-Michel, and C. G. R. Geddes. Generalized algorithm for control of numerical dispersion in explicit time-domain electromagnetic simulations. *Physical Review Special Topics-Accelerators And Beams*, Apr 2013. doi:10.1103/PhysRevSTAB.16.041303.

A. Pukhov. Three-dimensional electromagnetic relativistic particle-in-cell code VLPL (Virtual Laser Plasma Lab). *Journal of Plasma Physics*, 61(3):425–433, Apr 1999. doi:10.1017/S0022377899007515.

R. Lehe, A. Lifschitz, C. Thaury, V. Malka, and X. Davoine. Numerical growth of emittance in simulations of laser-wakefield acceleration. *Physical Review Special Topics-Accelerators And Beams*, Feb 2013. doi:10.1103/PhysRevSTAB.16.021301.

I. Haber, R. Lee, H. H. Klein, and J. P. Boris. Advances In Electromagnetic Simulation Techniques. In *Proc. Sixth Conf. Num. Sim. Plasmas*, 46–48. Berkeley, Ca, 1973.

J.-L. Vay, I. Haber, and B. B. Godfrey. A domain decomposition method for pseudo-spectral electromagnetic simulations of plasmas. *Journal of Computational Physics*, 243:260–268, Jun 2013. doi:10.1016/j.jcp.2013.03.010.

J. M. Dawson. Particle Simulation Of Plasmas. *Reviews Of Modern Physics*, 55(2):403–447, 1983. doi:10.1103/RevModPhys.55.403.

Q. H. Liu. The PSTD Algorithm: A Time-Domain Method Requiring Only Two Cells Per Wavelength. *Microwave And Optical Technology Letters*, 15(3):158–165, Jun 1997. doi:10.1002/(Sici)1098-2760(19970620)15:3<158::Aid-Mop11>3.3.Co;2-T.

Y. Ohmura and Y. Okamura. Staggered Grid Pseudo-Spectral Time-Domain Method For Light Scattering Analysis. *Piers Online*, 6(7):632–635, 2010.

H. Vincenti and J.-L. Vay. Detailed analysis of the effects of stencil spatial variations with arbitrary high-order finite-difference Maxwell solver. *Computer Physics Communications*, 200:147–167, Mar 2016. URL: https://apps.webofknowledge.com/full{\_}record.do?product=UA{\&}search{\_}mode=GeneralSearch{\&}qid=1{\&}SID=1CanLFIHrQ5v8O7cxqV{\&}page=1{\&}doc=2, doi:10.1016/j.cpc.2015.11.009.

H. Abe, N. Sakairi, R. Itatani, and H. Okuda. High-Order Spline Interpolations In The Particle Simulation. *Journal of Computational Physics*, 63(2):247–267, Apr 1986.

A. B. Langdon. On Enforcing Gauss Law In Electromagnetic Particle-In-Cell Codes. *Computer Physics Communications*, 70(3):447–450, Jul 1992.

B. Marder. A Method For Incorporating Gauss Law Into Electromagnetic Pic Codes. *Journal of Computational Physics*, 68(1):48–55, Jan 1987.

J.-L. Vay and C. Deutsch. Charge Compensated Ion Beam Propagation In A Reactor Sized Chamber. *Physics of Plasmas*, 5(4):1190–1197, Apr 1998.

C. D. Munz, P. Omnes, R. Schneider, E. Sonnendrucker, and U. Voss. Divergence Correction Techniques For Maxwell Solvers Based On A Hyperbolic Model. *Journal of Computational Physics*, 161(2):484–511, Jul 2000. doi:10.1006/Jcph.2000.6507.

J. Villasenor and O. Buneman. Rigorous Charge Conservation For Local Electromagnetic-Field Solvers. *Computer Physics Communications*, 69(2-3):306–316, 1992.

T. Z. Esirkepov. Exact Charge Conservation Scheme For Particle-In-Cell Simulation With An Arbitrary Form-Factor. *Computer Physics Communications*, 135(2):144–153, Apr 2001.

R. L. Morse and C. W. Nielson. Numerical Simulation Of Weibel Instability In One And 2 Dimensions. *Phys. Fluids*, 14(4):830–&, 1971. doi:10.1063/1.1693518.

H. R. Lewis. Variational algorithms for numerical simulation of collisionless plasma with point particles including electromagnetic interactions. *Journal of Computational Physics*, 10(3):400–419, 1972. URL: http://www.sciencedirect.com/science/article/pii/0021999172900447, doi:http://dx.doi.org/10.1016/0021-9991(72)90044-7.

B. B. Godfrey and J.-L. Vay. Numerical stability of relativistic beam multidimensional \PIC\ simulations employing the Esirkepov algorithm. *Journal of Computational Physics*, 248(0):33–46, 2013. URL: http://www.sciencedirect.com/science/article/pii/S0021999113002556, doi:http://dx.doi.org/10.1016/j.jcp.2013.04.006.