Issue 
A&A
Volume 558, October 2013



Article Number  A133  
Number of page(s)  24  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201321557  
Published online  18 October 2013 
AparT: code, validation, and physical interpretation of particleincell results
^{1} École Normale Supérieure, CRAL, UMR CNRS 5574, Université de Lyon, Lyon, France
email: mickael.melzani@enslyon.fr
^{2} CSCS Lugano, Switzerland
^{3} NVIDIA Switzerland, Technoparkstr 1, 8005 Zürich, Switzerland
Received: 23 March 2013
Accepted: 25 July 2013
We present the parallel particleincell (PIC) code AparT and, more importantly, address the fundamental question of the relations between the PIC model, the VlasovMaxwell theory, and real plasmas. First, we present four validation tests: spectra from simulations of thermal plasmas, linear growth rates of the relativistic tearing instability and of the filamentation instability, and nonlinear filamentation merging phase. For the filamentation instability we show that the effective growth rates measured on the total energy can differ by more than 50% from the linear cold predictions and from the fastest modes of the simulation. We link these discrepancies to the superparticle number per cell and to the level of field fluctuations. Second, we detail a new method for initial loading of MaxwellJüttner particle distributions with relativistic bulk velocity and relativistic temperature, and explain why the traditional method with individual particle boosting fails. The formulation of the relativistic Harris equilibrium is generalized to arbitrary temperature and mass ratios. Both are required for the tearing instability setup. Third, we turn to the key point of this paper and scrutinize the question of what description of (weakly coupled) physical plasmas is obtained by PIC models. These models rely on two building blocks: coarsegraining, i.e., grouping of the order of p ~ 10^{10} real particles into a single computer superparticle, and field storage on a grid with its subsequent finite superparticle size. We introduce the notion of coarsegraining dependent quantities, i.e., quantities depending on p. They derive from the PIC plasma parameter Λ^{PIC}, which we show to behave as Λ^{PIC} ∝ 1/p. We explore two important implications. One is that PIC collision and fluctuationinduced thermalization times are expected to scale with the number of superparticles per grid cell, and thus to be a factor p ~ 10^{10} smaller than in real plasmas, a fact that we confirm with simulations. The other is that the level of electric field fluctuations scales as 1/Λ^{PIC} ∝ p. We provide a corresponding exact expression, taking into account the finite superparticle size. We confirm both expectations with simulations. Fourth, we compare the VlasovMaxwell theory, often used for code benchmarking, to the PIC model. The former describes a phasespace fluid with Λ = + ∞ and no correlations, while the PIC plasma features a small Λ and a high level of correlations when compared to a real plasma. These differences have to be kept in mind when interpreting and validating PIC results against the VlasovMaxwell theory and when modeling real physical plasmas.
Key words: plasmas / methods: numerical / relativistic processes / instabilities / magnetic reconnection
© ESO, 2013
1. Introduction
Full particle codes are now used by a large number of groups worldwide to study plasmas out of equilibrium in a large variety of environments, including electronic devices, inertial fusion (Dieckmann et al. 2006; Bret et al. 2010), tokamaks, Earth and solar magnetospheres (Hesse et al. 2001; Daughton et al. 2006; Drake et al. 2006; Klimas et al. 2008; Markidis et al. 2012; Baumann & Nordlund 2012), or highenergy astrophysics (Silva et al. 2003; Pétri & Lyubarsky 2007; Zenitani & Hoshino 2008; Cerutti et al. 2012; Jaroschek et al. 2005; Nishikawa et al. 2008; Sironi & Spitkovsky 2011).
Particle simulations actually appeared along with the first computers at universities and in industry around 1950. They first concerned electron beams in vacuum tubes, a device extensively used in the computers themselves. The beams were cold (Hartree 1950) or later hot (Tien & Moshman 1956), and consisted of roughly 300 electron slabs moving in one dimension.
The step to plasma simulations was taken by Buneman (1959). He simulated an electrostatic plasma of 512 ions and electrons in one dimension, and showed that particle codes could be used to study the linear, nonlinear, and saturation phases of instabilities. At the time, the relevance of simulations with so few particles per Debye sphere was not clear and in 1962 Dawson (1962) and Eldridge & Feix (1962) made an important contribution by showing that correct thermal behavior was produced.
All these algorithms used particleparticle interactions, and the first particlemesh codes to introduce a grid appeared only later (Burger 1965; Hockney 1966; Yu et al. 1965). A great deal of literature on the drawbacks and benefits of the grid then appeared, and is now mostly concentrated in the two reference books of Birsdall & Langdon (1985) and Hockney & Eastwood (1988). Refinements of the algorithms quickly appeared (higher order grid interpolation, quiet codes, etc.), as well as code optimizations, at a time when programs were written in assembly language and depended heavily on machine architecture. Fully electromagnetic, relativistic, and 3D codes appeared with the studies of laser induced fusion (Buneman 1976), and closely resemble today’s codes.
In 1967, it was possible for Birdsall (1967) to list the papers concerning simulations, and he noticed that they had grown exponentially since 1956. In 1956, “many particles” meant 300; in 1985, 10^{6}; in 2008, 10^{9}; and in 2012, 10^{12}. For more historical details on PIC simulations, one may consult Birdsall (1967, 1999), Hockney & Eastwood (1988, Sect. 9.1), or the introduction of Birsdall & Langdon (1985).
Today, the latest generation of PIC algorithms consists of large versatile codes, featuring high order integration schemes and efficient parallelization. We can quote codes such as TRISTANMP (Spitkovsky 2005), OSIRIS (Fonseca et al. 2002, 2008), VORPAL (Nieter & Cary 2004), WARP (Grote et al. 2005), ALaDyn (Benedetti et al. 2008), iPIC3D (Markidis et al. 2010), PhotonPlasma (Haugboelle et al. 2013), or Zeltron (Cerutti et al. 2013). They employ various simulation methods. The equations solved can differ: electrostatic codes, electromagnetic codes, Darwin approximation (Huang & Ma 2008). The integration scheme can also vary: integration directly on the fields staggered on a grid, either with a charge conserving scheme or via solution of Poisson equation (Cerutti et al. 2012); or integration of the potential and correction of the discrepancies to charge conservation (Daughton et al. 2006). The time integration can be explicit or implicit (Markidis et al. 2012). Special parts of the numerical scheme can also differ; for example, the order of the field integration or the use of Fourier transforms. The interpolation of particle quantities to grid points and reciprocally can be done by a nearest grid point method (NGP), by a linear weighting (cloud in cell (CIC), or the PIC algorithm in the old terminology), or by a smoother shape (spline interpolation, Esirkepov 2001).
This article documents a new PIC code, AparT. Its first version, Tristan, was written by O. Buneman in 1990 (Matsumoto & Omura 1993). It was then made parallel by Messmer (2001), and used in Messmer (2002) or Paesold et al. (2005). Given the large quantity of simulation methods, it is useful to first detail our algorithm. This is done in Sect. 2 and Appendix A.
In Sect. 3, we present a set of test problems and the results obtained with AparT. The first test is the study of the fluctuation spectra of a thermal plasma. The second and third tests explore the linear and nonlinear stages of the filamentation instability. A last test is the computation of the linear growth rates of the relativistic tearing instability, for which we give the general equilibrium relations for the relativistic Harris configuration with arbitrary temperature and mass ratio between the two species (Appendix C). In Sect. 4 we describe a method for loading a MaxwellJüttner momentum distribution with an arbitrary bulk velocity and temperature. The naive method, which initializes the comobile distribution and then boosts particles individually, is shown to be incorrect, mainly because space contraction is absent from the PIC code.
In Sect. 5, motivated by the above tests, we explore the physical implications of two main differences between PIC plasmas and real plasmas: coarsegraining, i.e., a reduction in the number of particles by a factor reaching p ~ 10^{10}, and the representation of the fields on a grid. We introduce the notion of coarsegraining dependent quantities, i.e., physical parameters depending on the number of real particles represented by a single computer superparticle. The prototype of these quantities is the plasma parameter Λ, which scales for a PIC plasma as 1/p. The dependence on coarsegraining has important consequences for key physical quantities such as the thermalization time, the collision time, or the level of field fluctuations.
We point out more generally that a PIC code simulates a microstate constituted by a restricted number of finitesized particles each representing up to 10^{10} real plasma particles, while the VlasovMaxwell system models a plasma macrostate described by a continuous fluid in sixdimensional phasespace. These two descriptions are not equivalent and in particular PIC systems, with their small numbers of particles per Debye sphere, suffer from abnormally high noise levels and include to an unknown degree particle correlations absent from VlasovMaxwell equations.
We conclude in Sect. 6, stressing that validation and interpretation of PIC simulations requires detailed knowledge of the code as presented here as well as awareness of the above differences between the PIC model, the VlasovMaxwell model, and a real plasma.
2. Physical model and numerical implementation
This section presents the numerical scheme used in AparT. Broadly speaking, AparT is a parallel electromagnetic relativistic 3D PIC code with a staggered grid, where the fields are integrated via Faraday and MaxwellAmpère equations, currents are computed by charge conserving volume weighting (CIC), and fields are interpolated with the same CIC volume weighting method.
2.1. The PIC plasma
The code simulates the timeevolution of charged particles under the action of the electromagnetic fields that they generate, and the evolution of these fields.
Plasmas in nature contain millions to tens of billions of particles per Debye sphere, and relevant microphysical phenomena spread over numerous Debye lengths. It is impossible to track these particles one by one. Rather, the numerical particles represent numerous real particles, and are consequently called superparticles.
A superparticle represents either p real ions (having then a rest mass m_{sp} = p × m_{i} and a charge q_{sp} = p × q_{i}), or p real electrons (having then a rest mass m_{sp} = p × m_{e} and a charge q_{sp} = p × q_{e}). The ratio of ion to electron charge is always q_{i}/q_{e} = −1, while that of their rest masses m_{i}/m_{e} can be freely specified. Pair plasmas can thus be simulated. In AparT the number of real particles per superparticles p is the same for all superparticles at all times, but other codes can introduce superparticle splitting (Fujimoto & Machida 2006; Haugboelle et al. 2013; Cerutti et al. 2013).
We denote the physical size of a grid cell by X_{0}, a reference number of superparticles per cell by (including both ion superparticles and electron superparticles), and its associated number density of electrons by . Initially, the plasma is assumed to be quasineutral, in the sense that we load the same number of ion superparticles and electron superparticles in each cell. We have the relation (1)The equations governing the superparticle plasma are the equation of motion with the Lorentz force for each superparticle, and Maxwell equations coupled to the superparticle motions by the current: Here, c is the speed of light, e the electric field, b = cB is c times the magnetic field; q_{sp}, m_{sp}, γ_{sp}, v_{sp}, and x_{sp} are the charge, mass, Lorentz factor, velocity, and position of the superparticle number sp. The fields are stored on a grid, and a consequence is that the superparticles are seen by the grid as having a finite shape S, linked to the interpolation scheme used in the code.
In AparT we do not integrate the two other Maxwell equations because if they hold initially, then fulfilling the equation of conservation of charge and Eqs. (2c) and (2d) at all times insures that they remain correct to roundoff errors. That the current is indeed computed in a charge conserving way is detailed in Appendix A.1. Initially, the fields and the charge density are correctly built by setting a magnetic field satisfying ∇ ∧ b = μ_{0}j, a null electric field, and by placing the superparticles by pairs with one ion superparticle and one electron superparticle on top of each other so that the charge density is zero.
2.2. Numerical implementation of AparT
The code is based on the PIC program Tristan written by Buneman (Matsumoto & Omura 1993), parallelized by Messmer (2001). We largely modified its structure, that now uses Fortran modules and allows for more flexibility, including the possibility of switching between different initial conditions or boundaries by changing entries in a configuration textfile.
The need for a deep understanding of simulation methods to interpret their results, as is the case for the test problems of Sect. 3, motivated a detailed description of AparT. This section presents important points. More details (numerical scheme and parallelization efficiency) can be found in Appendix A.
Briefly, the global integration scheme of Eqs. (2a) − (2e) is a timecentered and timereversible leapfrog scheme, and is secondorder accurate in time and space. The electric and magnetic fields are stored on a staggered Yee lattice, which allows for a simple integration explicit in time of Eqs. (2c) and (2d) (without the current). The current is computed with the volume change of the superparticles in the grid cells, and added in a time explicit way to the integration of the electric field (Eq. (2d)).
2.2.1. Temporal and spatial discretization, normalization
The spatial discretization of the code is , a fraction n_{x} of the electron skin depth , where the electron plasma pulsation is , with − e and m_{e} the electron charge and rest mass. The timestep Δt is a fraction n_{t} of the electron plasma period: , with . We stress that the superscript is used for quantities based on the reference density .
Spatial quantities are normalized by the cell length X_{0}, and normalized quantities are then denoted with a tilde. For example, the electron Debye length (with ) has for normalized counterpart .
2.2.2. Superparticle volume
The use of a grid for PIC algorithms implies that the fields are known at grid nodes, and that information relative to the superparticles (charge and current) need to be interpolated on the grid. This interpolation is equivalent to considering the superparticles as clouds of charge of finite extension (Langdon 1970). The shape of the cloud then determines the interpolation formula.
In our case, a superparticle is assumed to be a cube of volume , and the current it produces is calculated by the change of the volume of the superparticle in the cell containing its center. This interpolation scheme is equivalent to linear weighting (CIC) and is exactly charge conserving. Details of the numerical implementation of the current computation can be found in Appendix A; in Sect. 5 the implications of the superparticle finite sizes are discussed.
2.2.3. Particle initialization in momentum space
Several momentum distributions can be loaded: MaxwellBoltzmann distribution with anisotropic temperature, boosted MaxwellBoltzmann distribution, waterbag distribution, or MaxwellJüttner distribution. We have not found in the literature a method for initializing the MaxwellJüttner distribution when both the bulk velocity and the temperature are relativistic, so we present one in Sect. 4.
2.2.4. Input, output, and data analysis
Data output and data analysis can be very time consuming for large scale simulations. To reduce data storage and writing time, we implemented parallel output in HDF5 format^{1}. Files are written according to the .h5part format, and can be read by the advanced visualization and dataanalysis software VisIt (Childs et al. 2012), which is fully parallel, but also by a reader in Python. The .h5part files can either contain the whole simulation data, and then be used to restart the simulation, or be lighter with only a fraction of the particles. These lighter files then include cellaveraged quantities related to the particles, such as the mean particle velocity or number, temperature, highest energy, etc.
We implemented the VisIt Libsim in situ library (Whitlock et al. 2011) into AparT. In this way, VisIt is able to connect to the simulation while it is running and to access the solver’s data at the current timestep. It can then perform data visualization and data analysis without the need to write data on the hard drive. This feature is fully parallelized, by exploiting the datadistribution model of the code, and as such is not restricted to the current parallelization model. It allows the dataIO from memory to hard drive – which is a major reason for slowdown of simulations using big datasets – to be significantly reduced. For example, a volumerendering of 3D data performed in situ takes less time than dumping 10 GB of data to the disk. VisIt in situ is also well suited to monitor ongoing simulations and to singlestep through the execution and is, in this way, of great help for debugging.
Finally, a set of test problems has been implemented. Based on a Python script, these problems can be automatically run to check code sanity after modifications.
3. Examples and code validation
3.1. Cold plasma modes
A first test is to simulate a thermal plasma at rest and to observe its modes of oscillation. It is an easy test if we focus on the pulsations for modes of zero wavevector, k = 0. To do so, we compute at each timestep t_{j} the sum of the momentum of the particles along a given direction, for example the xdirection, ∑ _{sp}γ_{sp}(t_{j}) v_{sp,x}(t_{j}), where the summation runs over all the electron superparticles in the simulation (of Lorentz factor γ_{sp} and velocity v_{sp}). This sum is equivalent to the volume integral of the momentum, and is thus equal to the k = 0 Fourier mode of the electron momentum, with a spectral resolution of 2π/ [box size in units of electron skin depth] . We then perform a Fourier transform in time to extract the pulsations of oscillation, (3)and similarly along y and z.
For all the simulations of this section, the initial state consists of a homogeneous plasma at rest, with superparticles loaded according to a classical MaxwellBoltzmann distribution. A uniform background magnetic field along z is set up for the magnetized plasma case. Periodic boundaries are used.
When the background magnetic field B_{0} is strong (electron cyclotron pulsation larger than electron plasma pulsation, ω_{ce} ≫ ω_{pe}), the particle trajectories are Larmor gyrations in B_{0}, unperturbed by collective effects such as Langmuir oscillations. This example then probes the accuracy of the particle motion integrator. When the background magnetic field is weaker (ω_{ce} ≲ ω_{pe}), the dynamics is set by Langmuir oscillations possibly modified by B_{0}. These oscillations involve the creation of electric fields by local charge imbalance. The fields set the particles into motion, and the particles then oscillate because of their finite inertia. Several parts of the algorithm are thus involved: electric field production and propagation, as well as particle motion.
3.1.1. No background magnetic field
With no background magnetic field, the only cold modes are the electromagnetic transverse wave of dispersion relation and the electrostatic Langmuir oscillation at the plasma frequency ω_{P} = ω_{pe}(1 + m_{e}/m_{i})^{1/2}. The latter is modified by thermal effects to a wave of dispersion relation . Consequently, we expect F_{a}(ω) to peak at ω_{La}(k = 0) = ω_{Tr}(k = 0) = ω_{P}.
Our simulations span a large range of parameters: n_{x} and n_{t} (spatial and temporal resolution, see Sect. 2.2.1) vary between 5 and 50 and between 300 and 2000, respectively; ρ_{sp} (the number of superparticles per cell) varies between 4 and 32; and the thermal velocity of the electrons between 0.04c and 0.1c (the ions have the same temperature as the electrons). This results in Debye lengths between 0.2 and 2 cells. We checked that the simulation box size, comprised between 10 and 30 Debye lengths, does not influence the results.
We use a mass ratio m_{i}/m_{e} = 49, which results in ω_{P} = 1.010 ω_{pe}. Our simulations last 100 T_{pe}, so that the frequency resolution is Δω = 2π/(100T_{pe}) = 0.01ω_{pe}.
For all these parameters, we find that the position and width of the frequency peak are always the same as in Fig. 1. It coincides with the theoretical plasma pulsation, which is expected because our temporal spectra are for the wavenumber k = 0, and because ω(k = 0) = ω_{P} for the two modes present in this situation.
This is the case even for simulations where the Debye length is not resolved. However, an underresolved Debye length leads to more numerical heating and can trigger instabilities in situations less trivial than a thermal plasma at rest (see Appendix B), so that we have not pushed our investigations too far in this direction.
The main difference between the simulations is that less resolved ones present noisier spectra, and thus more fluctuations. The increase of fluctuation level with decreasing resolution is a universal feature of PIC simulations and is explored in more detail in Sect. 5.
Fig. 1 Power spectra of the total momentum of the electrons  F_{a}  ^{2} (Eq. (3)) with a = x, y, or z in the unmagnetized case. The inset is a zoom around the peak. Here n_{x} = 25, n_{t} = 500, ρ_{sp} = 16, v_{th,e} = 0.04c, T_{i} = T_{e}, box size of 25 × 25 × 25 cells, duration of 100 T_{pe}. 
Theoretical versus experimental pulsations for a magnetized cold plasma.
3.1.2. With a background magnetic field
In a uniform and cold magnetized plasma, the plasma modes depend solely on the ratio of the electron cyclotron pulsation ω_{ce} = eB/m_{e} to the electron plasma pulsation ω_{pe} (for a fixed m_{i}/m_{e}). This ratio sets the relative importance of individual particle motion (ω_{ce}) against collective effects (ω_{pe}). In addition, the background magnetic field favors a direction, thus making the mode spectrum anisotropic (however, the mode pulsations for k = 0 remain independent of the direction of the wavevector).
For wavevectors parallel to the magnetic field, the Langmuir oscillation remains unchanged because the oscillations of the particles are longitudinal, and thus along B_{0} and unaffected by the magnetic field. On the other hand, the electromagnetic wave with k along B_{0} separates into two branches, one starting at and the other at (see, e.g., Fitzpatrick 2011, chap. 4), in both cases with particles oscillating in the transverse plane, i.e., perpendicular to B_{0}. Two other branches appear, but they start at ω_{k = 0} = 0 and will thus not appear in F_{a}(ω).
For wavevectors perpendicular to the magnetic field, the presence of the background magnetic field deforms the Langmuir oscillation ω(k) = ω_{P} into a branch starting from ω_{k = 0} = ω_{L}, with particle oscillations in the plane perpendicular to B_{0}. The transverse electromagnetic wave still exists with particle oscillations along B_{0}, unaffected by the magnetic field. Another branch appears, which is a deformation of the transverse electromagnetic wave for particle oscillations not along B_{0}, and starts at ω_{k = 0} = ω_{R} with oscillations in the plane perpendicular to B_{0}. Another branch appears, starting at ω_{k = 0} = 0.
All in all, we expect to find a peak at ω = ω_{P} for the component of the momentum parallel to B_{0}, and two peaks at ω = ω_{L} and ω_{R} for the component of the momentum perpendicular to B_{0}. We ran the set of simulations described in Table 1, with a ratio ω_{ce}/ω_{pe} ranging from 0.5 to 4, n_{x} from 16 to 128, n_{t} from 1000 to 4000, and ρ_{sp} from 4 to 16, and we did find the required pulsation peaks for F_{a}(ω) (see Fig. 2 for a sample spectrum). The positions and widths of these three peaks are almost constant within our parameter range.
Fig. 2 Power spectra of the total momentum of the electrons  F_{a}  ^{2} (Eq. (3)) with a = x, y, or z in the magnetized case. Parameters are given in Table 1 (case labeled Fig. 2). The purple lines are the pulsations , , and . 
We note that the peak positions and widths are not changed even for cases where the thermal Larmor radius is not resolved. It is expected that the resolution of the Larmor radius by the grid is of no importance to describe particle trajectories in constant fields, because the interpolation of these fields from grid points to superparticle position gives the same result regardless of the grid size if the fields are constant. The relevant constraint is instead that r_{ce} should be resolved along the trajectory, v_{sp}Δt < r_{ce}, with v_{sp} the superparticle velocity. This relation is equivalent to .
However, for simulations with underresolved Larmor radii the electric field energy starts to behave abnormally after some tens of plasma pulsations, and energy conservation curves present an exponential heating (see Sect. 3.1.3) that can lead to dramatic consequences. This parameter range must be avoided.
3.1.3. Energy conservation
Independent of the strength of the background magnetic field, we observe a linear increase of the total energy with time that is due to interactions of particles with the grid (see Appendix B.2). The growth rate is independent of the size of the timestep (from n_{t} = 2000 down to Courant condition Δt ~ X_{0}/c, equivalent to n_{t} = 2πn_{x}). Its dependence on the number of superparticles per cell is quite precisely given by 1/ρ_{sp}. However, its dependence on the spatial resolution n_{x} and thermal spread v_{th,e} is less clear. In particular, it does not depend only on the product . The rate increases with increasing v_{th,e}, and decreases with increasing n_{x}. Some examples are given for reference in Table 2.
Energy conservation for simulations of a thermal plasma with no background magnetic field.
Simulations with underresolved Larmor radii show an exponential (instead of linear) increase of the total energy starting after roughly 40 T_{pe}. This numerical instability is believed to arise because field perturbation at wavelength λ = 2r_{ce} and their aliases ( ± λ + n × 2π/X_{0}, with n an integer and X_{0} the grid size) are allowed to couple when 2X_{0} > λ (see Appendix B.2). An inspection of the energy curves shows that the energy gain is for the kinetic energy. We note that it did not disturbed the spectra of Sect. 3.1.2 because they were computed before the heating reached a significant level. It is interesting to note that this numerical instability is develops slowly, so that particles with underresolved Larmor radii in constant fields can be included in simulations if they spend a small amount of time before being heated or before reaching areas with smaller magnetic fields.
3.2. Linear growth rates of the counterstreaming instability
Another standard test is to study the linear phase of the counterstreaming instability. We use relativistic streaming velocities to validate the behavior of the algorithm for relativistic particle motions. Moreover, since magnetic fields are generated for this range of parameters, this test also probes the integration of b.
The initial setup consists of two unmagnetized and cold counterstreaming electronpositron beams, with velocity and associated Lorentz factor Γ_{0}. There is no background magnetic field, and the particles are loaded according to a drifting MaxwellBoltzmann distribution. This situation is unstable, and the kinetic energy of the beam is converted into particle thermal kinetic energy and electromagnetic field energy, the initial perturbation coming from fluctuations due to the finite superparticle number. The linear instability spectrum is described by a branch comprising an electrostatic longitudinal twostream mode and a transverse electromagnetic filamentation mode, the general case being an oblique mixed mode (Bret et al. 2010).
3.2.1. Theoretical model
We take the growth rates derived analytically by Michno & Schlickeiser (2010) on the basis of a cold twofluid model and denote this result as theoretical. Our thermal velocity v_{th}, identical for both species, is low enough to insure that thermal effects are negligible (Bret et al. 2010, Eq. (28)), but high enough to have a resolved Debye length and to avoid numerical instabilities (Appendix B). Our parameters are chosen such that the transverse filamentation mode always dominates. The fastest growing modes are those at large wavenumbers perpendicular to the beams, i.e., (with d_{e} = c/ω_{pe}), and that grow according to with (4)The k_{ ⊥ }dependance of the growth rate is plotted in Fig. 3. We see that all modes above a few quickly reach the maximum growth rate .
Fig. 3 Wavenumber dependance of the growth rate τ^{theory}(k_{ ⊥ }), represented here as the deviation . We recall that . From Eq. (69) of Michno & Schlickeiser (2010). The squares are the resolution in k_{ ⊥ } for a box of transverse size 9, 4.5, and 3 d_{e}. 
Fig. 4 Top: energy curves for the filamentation instability (case labeled Fig. 4 in Table 3). They are normalized by the total initial energy , which is mostly the kinetic energy of the particles. For example . The curve “energy conservation” is . After 8T_{pe}, the situation is more or less steady. Bottom: autocorrelation scale of the current amplitude. 
Fig. 5 Top: growth of individual Fourier modes for b_{x} (for the run labeled Fig. 4 in Table 3). The modes shown are those for i = 0,1,2, and j = 0,1,...,25, and one mode every 100 modes for the remaining. We note that the graphic has been cut and that the weak modes actually fill a continuum down to an energy of 10^{7}. The sum of all 320 × 90 modes is shown in orange. We recall that mode (i,j) corresponds to (k_{z}d_{e},k_{x}d_{e}) = 2π × 20(i/640,j/180). Bottom: growth map of the Fourier modes, in units of T_{pe}. 
3.2.2. Method of measurement
We measure the growth rates of the magnetic fields b_{x} and b_{y} with two methods. The first is a direct measure on the total energy curve, e.g., (see Fig. 4, top, for an illustration). It gives an effective growth rate that we denote by , equal to 0.48T_{pe} in this case.
The second consists in following the time evolution of the Fourier modes of the fields. At a fixed time t_{0}, we compute the 2D Fourier transform of the fields in a plane (x,z) with a fixed y, that we denote by FT_{y = y0}(t_{0},k_{x},k_{z}). We then average the power spectrum over all the planes y = cst to obtain the power spectrum PS(t_{0},k_{x},k_{z}) = ∑ _{y0}  FT_{y = y0}(t_{0},k_{x},k_{z})  ^{2}. We then repeat this procedure for several t_{0}. The discrete mode spectrum is sampled with (k_{z}d_{e},k_{x}d_{e}) = 2πn_{x}(i/N_{z},j/N_{x}), where N_{z} and N_{x} are the total number of cells in the z and x directions, and i = 0...N_{z}/2, j = 0..N_{x}/2. The spectral resolution in the direction perpendicular to the beam is thus Δk_{ ⊥ }d_{e} = 2πn_{x}/N_{x} = 2π/(box width in d_{e}). The squares in Fig. 3 represent this spectral resolution for the different box sizes that we use.
3.2.3. Results
Figure 5 is an example of the temporal evolution of the modes of b_{x} for the same simulation as in Fig. 4. The sum of all modes grows at the same effective growth rate as the total energy in b_{x} (to within ± 1%), . However, the fastest growing modes are those for k_{z} = 0 and 0 ≤ k_{x}d_{e} ≤ 5, with , which is close to the coldfluid result . It is seen from Fig. 5 that the large difference between the effective growth rate and the growth rate of the fastest modes is due to a significant contribution of all the modes during the whole linear phase. The fastest mode thus never dominates the total energy in the linear phase. We suspect that this is due to the large noise level present in PIC simulations.
Fig. 6 Difference between the growth rate measured on the total energy curve and the theoretical growth rate of the fastest modes . Except when labeled otherwise, v_{th} = 0.1c. Each symbol corresponds to a transverse box size, and each color to fixed (n_{x},n_{t}). 
These results hold for all the test simulations that we conducted, which are summarized in Table 3. The effective growth rates measured on the total energy present various levels of discrepancies with , between 14% and 67%. Figure 6 shows the dependance of these discrepancies. There is a small sensitivity with respect to the timestep (n_{t}) and the box size, and an important influence of the spatial resolution (n_{x}). There is a systematic decrease in the difference when the superparticle number per cell ρ_{sp} is increased (all other parameters are kept constant). Since the fluctuation level in the PIC plasma decreases with increasing ρ_{sp}, this indicates that the high fluctuation level excites all the modes and prevents the fastest ones from dominating the energy.
Fig. 7 Three snapshots of the filamentation instability, at times 2.5, 4, and 5.5T_{pe} (from left to right) for the simulation of Fig. 4. Colors represent the current magnitude per cell (units are number of superparticles per cell times their mean velocity), while arrows are the magnetic field. Lengths are in cell number. 
Theoretical versus experimental values of the filamentation growth rate τ.
Theoretical versus experimental values of the tearing growth rate τ.
On the other hand, the growth rates measured on the fastest modes differ from by a more systematic factor, 14% for β_{0} = 0.95, 7 ± 1% for β_{0} = 0.995, and 6% for β_{0} = 0.999. These systematic differences can be explained by looking at the mode spectrum. In all the simulations, the fastest modes are for k_{z} = 0 and 0 ≤ k_{x}d_{e} ≤ 5 − 15 (see, e.g., Fig. 5, bottom). Given the spectral resolution Δk_{ ⊥ }d_{e} = 2π/(box length in d_{e}), these modes actually cover a portion of k_{ ⊥ }d_{e} where the curves τ^{theory}(k_{ ⊥ }) of Fig. 3 vary significantly and do not yet reach . It explains the sign and order of magnitude of the difference . It also explains the increasingly better agreement when β_{0} increases.
We note that e and b_{z} are zero in the linear twofluid theory, and the fact that they are not zero in our simulation (see Fig. 4) reflects an early nonlinear evolution or the effects of fluctuations and correlations absent from the linear model but present in PIC simulations. These differences are discussed further in Sect. 5.
3.3. Nonlinear evolution of the filamentation instability: filament merging
We now consider the nonlinear phase of the filamentation instability. We study the same counterstreaming configuration as in Sect. 3.2, with the setup labeled by Fig. 4 in Table 3. The energy curves are shown in Fig. 4.
As a diagnostic, we focus on the filament growing and merging processes. The linear phase of the filamentation instability produces current filaments. Since they are threaded by parallel currents, they attract each other and, starting from the end of the linear phase, start to merge to produce larger and larger filaments (Medvedev et al. 2005). This is clearly visible in Fig. 7.
We measure the size of the filaments by computing the twodimensional autocorrelation function, in the x − y plane, of the zaveraged current amplitude (z is the direction of the beams). This autocorrelation function is azimuthally averaged to obtain a radial function corr(r). We normalize corr(0) to 1. The scale of the filaments is then taken to be five times the radius where corr(r) = 0.8 (taking this scale as the radius at which corr(r) first vanishes yields the same results).
The results are shown in Fig. 4 (bottom). We clearly see two regimes: one during the linear growth of the filamentation instability (from t = 2.2 to 3.6T_{pe}) where the filament correlation length is set by the wavelength of the fastest growing mode and remains constant, and one in the nonlinear regime (after t = 3.6T_{pe}) where the filaments merge and thus quickly increase their size. In the second case, the growth is roughly linear with time, which agrees with the PIC simulation results of Dieckmann (2009) for a similar setup. After t = 6.5T_{pe}, the filament growth stops. However, we suspect that the periodic boundaries start influencing the dynamics at this point.
3.4. Linear growth rates of the relativistic tearing instability
We also study the linear phase of the tearing mode for a relativistic Harris sheet in a pair plasma. Contrary to the preceding case, this example provides a test of the algorithm in a situation where thermal effects are essential.
The equilibrium consists of a magnetic field (5)sustained by a population of electrons and positrons of density ∝ 1/cosh^{2}(x/L) flowing with opposite bulk velocities U_{e} = −U_{i} in the ± y directions. We denote the associated Lorentz factor by Γ_{e}, and the temperature of the two species by Θ = T/(m_{e}c^{2}). The exact relations between the different parameters to satisfy the equilibrium are given in Appendix C. In particular, one should be careful to distinguish between quantities in the frame moving with one species (denoted with a prime) and quantities in the simulation frame. For example, contraction of the electron density leads to .
The superparticles are loaded according to a drifting MaxwellJüttner distribution with the method described in Sect. 4. There is no initial perturbation, and the instability grows out of the fluctuations produced by the finite number of superparticles.
As can be seen in Table 4, the bulk velocities and the temperatures of electrons and positrons are relativistic. Loading these distributions in a PIC code is a nontrivial task, and we have developed a special method for this. It is presented in Sect. 4.
The simulation domain is periodic along z and y. Reflecting boundaries for particles and fields are present along the x direction. There are no background particles, only that of the current sheet.
An example of the energy evolution is presented in Fig. 8. After some time, the system becomes unstable and the magnetic field starts reconnecting. As expected, b_{z} dwindles while b_{x} rises, which corresponds to the formation of magnetic islands. We measure the linear growth rate on b_{x} as . For comparison, we use the linear growth rates derived by Pétri & Kirk (2007) by linearizing the VlasovMaxwell system around the drifting MaxwellJüttner distribution C.1 and the magnetic field of Eq. (5).
The results are summarized in Table 4. Discrepancies with Pétri & Kirk (2007) range between 5% and 8%. These growth rates vary by less than 3% when the numerical resolution is doubled (i.e., when ρ_{sp}, n_{x}, and n_{t} are doubled all together). We restrict our analysis to total energy curves because contrary to the case of the filamentation instability, the linear growth of the field energy spans several orders of magnitude and the fastest modes have enough time to dominate the total energy.
Fig. 8 Energy curves for the tearing instability (labeled Fig. 8 in Table 4). Shown are the energy curves, for example ^{∫}dVB_{x}(t)^{2}/(2μ_{0}). They are normalized by the total initial energy (which is mostly the energy in b_{z}). The curve labeled “energy conservation” is . 
We note that when the simulation is launched, an electromagnetic wave is seen to propagate from the sheet in the ± x directions. This wave is a necessary consequence of the fact that at t = 0 we load the superparticles by pairs of electronpositron at the same location, and we set a zero electric field everywhere (see Sect. 2.1). The system then has to relax from this very peculiar state: in less than one plasma period, charge screening is established and an electric field appears. It is this field that partly propagates outside of the sheet. It is then reflected on the ± x boundaries and propagates back to the sheet, causing the oscillations in e_{y} seen in Fig. 8. We have checked that their incidence does not influence the linear growth by using different domain sizes.
It is more puzzling that the sheet contracts slightly just after the beginning of the simulation and the current magnitude at its center rises by ≲ 10%. This may be related to the fact that our algorithm does not solve the VlasovMaxwell system in a strict sense (see also Sects. 5 and 6.3).
4. Loading a relativistic particle distribution
This section deals with the general problem of loading particles with momenta that reproduce a given distribution function.
Very common cases are waterbag and MaxwellBoltzmann distributions with a mean bulk velocity U_{0}. A simple method is then to load the particles in the frame comoving with the plasma, which is fairly easy because in this frame the distributions are isotropic, and then to add to every particle the velocity U_{0} or, if U_{0} is close to c, to boost every particle with the Lorentz boost corresponding to U_{0}. We will see, however, that this method is no longer correct when both U_{0} and the rest frame distribution are relativistic, mainly because boosting particles in a PIC code does not boost space. We present here another method, correct in all cases^{2}.
4.1. Transformation of the distribution function
We start by explaining how the particle distribution changes from one frame to another (see e.g. Mihalas & Mihalas 1984; Pomraning 1973).
We consider a frame ℛ where the plasma has a mean velocity U_{0} and follows the distribution f(x,p). In the comoving frame or plasma rest frame ℛ_{0}, the plasma mean velocity is zero and follows the distribution f_{0}(x,p). We follow a group of particles. Seen from ℛ, they are in a volume d^{3}x around x and have momentum p with a scatter d^{3}p; seen from ℛ_{0} these quantities change respectively to d^{3}x_{0}, x_{0}, p_{0}, and d^{3}p_{0}. The number of particles in our group is (6)so that to find the link between f and f_{0} we have to find a relation between d^{3}x_{0}d^{3}p_{0} and d^{3}xd^{3}p.
We start with the momentum. In the rest frame, our group of particles have momenta spanning a range d^{3}p_{0}. Seen in the boosted frame, their momenta transform according to the Lorentz transformation, and span a new range d^{3}p. These two volumes are thus linked by the Jacobian of the Lorentz transformation, and it can be shown that (7)where γ and γ_{0} are the Lorentz factors associated to p and p_{0}.
We now consider the space volumes. Because of space contraction/dilatation, the group of particles will occupy a different volume in different frames. We consider the frame ℛ′ comoving with the group of particles. This is possible because the particles all move at nearly the same velocity v_{0} = p_{0}/γ_{0}. We denote by a prime all quantities seen from this frame. In ℛ′, the particles occupy a volume d^{3}x′. Since only one direction is contracted, and since ℛ′ moves relative to ℛ_{0} with Lorentz factor , we have the relation d^{3}x_{0} = d^{3}x′/γ_{0}. Similarly, ℛ′ moves relative to ℛ with Lorentz factor , and we have d^{3}x = d^{3}x′/γ. All in all: (8)From this, we deduce that (9)
4.2. Why boosting particles from the rest frame is incorrect for relativistic distributions
We now come back to PIC simulations. We assume that we load particles uniformly in space, with momenta following f_{0}(x_{0},p_{0}), and that we boost each particle with a velocity U_{0}. The momentum volume elements are then transformed according to Eq. (7), but positions are not changed. Equation (8) does not hold, and we obtain a particle distribution (10)The volume contraction/dilatation is not performed, and the factor γ_{0}/γ does not cancel.
Consequently, boosting each particle from the rest frame leads to the expected distribution only if γ_{0}/γ is independent of the particle. We can write this ratio as (11)with p_{0,y} the y component of p_{0} and , so that this is the case only if p_{0,y} ≪ 1 (or if the boost is nonrelativistic, U_{0} ≪ c). If p_{0,y} ≪ 1, then γ_{0}/γ ~ 1/Γ_{0} and when it is inserted back into Eq. (10), we find the usual result of density contraction.
However, when the particle distribution is relativistic in the rest frame of the plasma, γ_{0}/γ is not a constant factor and Eq. (10) does not have the expected dependence on momentum p.
4.3. Loading a drifting MaxwellJüttner distribution with arbitrary temperature and drift speed
We now present a method for loading the superparticle momenta directly in the frame where the distribution has a bulk velocity.
In the literature, there is some agreement around the fact that the particle distribution of a relativistic plasma in thermodynamic equilibrium is given by the MaxwellJüttner distribution (Jüttner 1911; Cubero et al. 2007; ChacónAcosta et al. 2010; Dunkel & Hänggi 2009) and, even if some alternatives are also debated (Treumann et al. 2011), this is the distribution used in PIC simulations (e.g., Pétri & Lyubarsky 2007; Zenitani & Hoshino 2008; Jaroschek & Hoshino 2009).
In the plasma rest frame, the MaxwellJüttner distribution is given by (12)with n_{0} the uniform particle number density and g_{0} the momentum distribution, both in the restframe, μ = mc^{2}/T, p = γv/c, and K_{2} the modified Bessel function of the second kind. If we now denote by f the distribution in the frame where the plasma moves with a bulk velocity cβ_{0} (and associated Lorentz factor Γ_{0}), by n the particle density, and by g the momentum distribution, both in this same frame, we have f(x,p) = ng(p). Equation (9) then gives g(p) = g_{0}(p_{0}) × n_{0}/n. Now using n = Γ_{0}n_{0}, we arrive at (13)This distribution is normalized to unity: .
Fig. 9 MaxwellJüttner distribution for T = mc^{2} and . Upper plot: contours are drawn from the exact expression 2πrg(r,0,p_{y}), with g from Eq. (13) and (r,0,p_{y}) defined in Appendix D. The background color map is the 2D histogram of for 10^{6} particles generated according to our method. Middle plot: normalized histogram of p_{y} for the particles (red dotted), to be compared to the exact expression in Eq. (D.2) (blue line), and for comparison (black dots) the histogram of p_{y} for particles generated the wrong way (initialization in the rest frame and boost of Γ_{0}). Bottom plot: difference between the red points and the blue curve. 
The main difficulty with Eq. (13) is that the variables p_{x}, p_{y}, and p_{z} are coupled and cannot be chosen independently. The solution is to compute the marginal distribution for the variable p_{y}. With this, one can choose the ycomponent of p independently of the others, and then use the distribution g(p_{ ⊥ },p_{y}) knowing p_{y} to choose the component normal to y. Details are presented in Appendix D, where we provide the expressions for the marginal distribution and the conditional probability, as well as a method for generating the velocity by using the cumulative distributions.
Figure 9 shows an example of the distribution generated with this method for T = mc^{2} and Γ_{0}β_{0} = 1.41, and compares it to the theoretical expectation and to the distribution obtained by boosting individually the superparticles from the restframe. We clearly see the accuracy of our algorithm, and the mismatch between the simpler boosting method and the expected result. We note that this mismatch can have significant consequences. For example, when used for the tearing instability of Sect. 3.4 it leads to large adjustments in the initial conditions, and for the most extreme temperatures to a complete disruption of the current sheet.
We also provide in Table 5 the analytical expressions of various quantities averaged over the MaxwellJüttner distribution, for example, the mean Lorentz factor, mean momentum, or the enthalpy. The method used to derive these expressions is presented in Appendix E. Besides being of general interest, these expressions served to further validate our generated distributions by comparing the analytical results of Table 5 with averages performed over particles generated with our method.
5. Real, PIC, and VlasovMaxwell plasmas
Particleincell simulations have brought tremendous new insights into astrophysical plasmas, for example through studies of kinetic instabilities in their nonlinear phase, kinetic turbulence, particle acceleration via the Fermiprocess, or 3D magnetic reconnection and the associated particle acceleration. However, as we will detail in this section, there remain a number of questions with respect to the degree to which PIC models are able to completely mirror real plasmas.
The modeling of a real plasma by a PIC plasma implies two steps (see the right branch of Fig. 12): the grouping of many real particles into a single superparticle, known as coarsegraining, and the discretization of the equations with the presence of a grid. Each of these steps raises questions:

1.
Coarsegraining: Is plasma behavior still expected with so fewsuperparticles per Debye sphere? Is the noise level too large?With this, do nonlinearities appear sooner than in real plasmas?Does the PIC plasma remain collisionless? And what are welosing when we gather the particles into the superparticles?

2.
Discretization and grid: At least for explicit schemes, they bring with them numerical stability problems, reviewed in Appendix B. Moreover, the interpolation of superparticle quantities to grid points implies a finite volume for the superparticles, which in turn implies a vanishing twopoint force at short distances and thus reduces drastically the influence of collisions; it helps the PIC plasma to be collisionless, but is it enough? And what are the consequences of having superparticles whose sizes reach a significant fraction of the Debye length?
We discuss some of these questions in Sects. 5.1 to 5.3.
The distinguishing feature of the VlasovMaxwell description of a plasma is the absence of collisions and of correlations between particles. Given the two preceding points, we can wonder if a PIC plasma can be described by the VlasovMaxwell system, or if it has too few superparticles per cell and thus correlation levels that are too high for this description to be accurate. The differences between PIC and VlasovMaxwell descriptions are examined in Sect. 5.4.
5.1. The plasma parameter Λ: a coarsegraining dependent quantity
As said earlier and expressed in Eq. (1), a real plasma is represented in the computer by grouping many particles into superparticles. This is what is called coarsegraining. Unlike fluid equations^{3}, Eqs. (2a) − (2e) are not invariant under coarsegraining (because of the definition of the current, Eq. (2e)). The prototype of pdependent quantities is the plasma parameter^{4}, which is close to the number of particles per Debye sphere (we recall that p is the number of particles per superparticles, λ_{D} the Debye length and n the real particle number density).
The plasma parameter Λ also expresses the ratio of the particles’ kinetic energy to their electrostatic potential energy of interaction and, as such, varies as 1/p because kinetic energy is proportional to the superparticles’ mass m_{sp} ∝ p while charge interaction energy involves their charge . This can be seen directly by writing for the superparticle plasma, with n_{sp} the number density of superparticles. The Debye length, being derived from fluid theory, is invariant under coarsegraining, and since n = p × n_{sp}, one has that (14)with Λ = Λ_{p = 1} the real plasma parameter.
In a real plasma Λ ranges from 10^{4} to 10^{20} (for example Λ ~ 10^{6} in solar coronal loops; 10^{12} in the magnetotail, magnetopause, or in typical Crab flares; 10^{17} in AGN jets), while in computer experiments where we have to simulate thousands to millions of Debye spheres, Λ reaches hardly a few tens (for a discussion see, e.g., Bykov & Treumann 2011, Sect. 4). The corresponding number of particles per superparticles then reaches p ~ 10^{3} to 10^{19}. The question of the relevance of PIC simulations for describing collisionless plasmas has thus been asked from the beginning, and concerns both terms, collisionless and plasma, that we now discuss.
5.1.1. Plasma behavior
A weakly coupled plasma is characterized by the predominance of collective effects over individual effects. The ratio of these effects is contained in the plasma parameter Λ_{p}, which can be seen as the ratio of collective behavior (the interaction of one particle with the electromagnetic fields collectively generated by all others, which is coarsegraining independent) to binary effects (which are proportional to ).
Since plasma behavior, with Debye screening and local charge neutrality, requires a high plasma parameter, it is wise to ask how large it should be in a PIC plasma. Birsdall & Langdon (1985, chap. 1) and Hockney & Eastwood (1988) have shown that it is not necessary for this ratio to be as high as in real plasmas, and that a Λ_{p} of about a few suffices for correct plasma behavior.
5.1.2. Collisionless behavior
A plasma behaves collisionlessly if the time and length scales of interest are negligible toward the collision time and the meanfree path, respectively. Since the collision time scales as Λ_{p} times the plasma period, it is not clear whether a PIC computer plasma with a plasma parameter on the order of unity will be collisionless.
Particleincell plasmas are helped by the superparticle finite sizes, which imply that the twopoint force decreases to zero for separations smaller than this size. This fact, albeit degrading the accuracy of single particle dynamics, greatly reduces the relative importance of binary collisions so that in order to correctly simulate a collisionless plasma for scales accessible in simulations, one has to insure that (Birsdall & Langdon 1985; Hockney & Eastwood 1988, chap. 1) (15)(Here, is the normalized Debye length. It can be expressed as with .)
This is all the more true if r_{c} < X_{0}, where X_{0} is the grid size and r_{c} is the effective collision radius for Coulomb encounters, expressed by equating the kinetic energy of the meeting particles to their potential energy of interaction: (r_{c} is also pdependent). The Debye length must be resolved for reasons of numerical stability (Sect. 3.1 and Appendix B), so that we arrive at an optimal ordering r_{c} < X_{0} < λ_{D}, which is allowed only if, again, λ_{D}/r_{c} = Λ_{p} > 1.
However, because the PIC model is a description of a plasma of cloud charges, grazing collisions will still be present and will lead to thermalization (Birsdall & Langdon 1985, Chap. 1; Hockney & Eastwood 1988, Chaps. 1 and 9), a point discussed in the next section.
5.2. Thermalization time: a coarsegraining dependent quantity
In a PIC plasma, the behavior of plasma quantities depending on Λ can be guessed by replacing Λ by Λ_{p}. This is the case for the thermalization time of a plasma by grazing Coulomb collisions (Spitzer 1965) or by electric field fluctuations (Birsdall & Langdon 1985, p. 282), which is on the order of t_{th} ~ T_{P} × Λ (with T_{P} the plasma period). This has two important consequences:

We expect t_{th} to depend on resolution and coarsegraining, roughly as (16)

Since Λ_{p} = Λ/p is several orders of magnitude smaller than the real plasma parameter Λ, we expect the thermalization by grazing collisions and fluctuations to be vastly more efficient in PIC codes than in reality.
This can have important consequences in simulations where thermalization plays a key role. For example in real collisionless shocks, the mean free path for collisions l_{mean free path} is far larger than the shock thickness Δ_{shock} and the thermalization processes are collisionless kinetic instabilities. Since the mean free path in a PIC plasma is smaller by a factor of p ~ 10^{10} than in a real plasma, it is not obvious that the ordering still holds. To truly describe a collisionless shock with a PIC algorithm, one has to be careful that the unphysically fast thermalization by collisions or fluctuations remains slower than thermalization by kinetic instabilities.
Fig. 10 Top: electron temperatures for the hot (2) and cold (1) plasmas, from four simulations with different ρ_{sp}. The curves for ion temperatures are similar, except for an overall time dilatation by a factor ~ (m_{i}/m_{e})^{1/2} = 5. In this figure we use T_{∞}(t) = (T_{1}(t) + T_{2}(t))/2. Except for ρ_{sp} = 4 where there is significant numerical heating, T_{∞}(t) is constant in time. Bottom: halfthermalization time for ions and electrons, versus number of superparticles per cell. For ions, we have plotted t_{th}/(m_{i}/m_{e})^{1/2}. The times reported are measured as the initial slope of the temperature curves in a loglin plot, and thus correspond to t_{th}/2. We see the scaling t_{th} ∝ ρ_{sp}. 
To illustrate the dependance of the collision and fluctuation induced thermalization time, we present simulations that initially have two thermal ionelectron plasmas. The first is cold, with a temperature T_{1,e}(0) = T_{1,i}(0) = 1.6 × 10^{3}m_{e}c^{2} for its electrons and ions, while the second is hot, with T_{2,e}(0) = T_{2,i}(0) = 1.8 × 10^{2}m_{e}c^{2}. The mass ratio is m_{i}/m_{e} = 25. The four species interact via collisions and correlations (no sign of plasma kinetic instabilities were found) and tend to reach the same final temperature (17)Particleincell results are shown in Fig. 10 (top) for the electrons, for four simulations with a number of superparticles per cell (including all species) ρ_{sp} = 4, 16, 64, or 128. The other parameters are kept fixed: n_{x} = 25, n_{t} = 500, and box size of 25^{3} cells. It results in , 61, 243, or 485. The temperatures are measured with , where the sum runs over all the superparticles of a given species. In Fig. 10 (top) we clearly see a slower thermalization as ρ_{sp} increases.
To evaluate the thermalization times, we use the result of Spitzer (1965): for two species at temperature T_{1} and T_{2}, thermalization occurs according to (18)with n = n_{1} = n_{2} the particle number density, , and lnΛ^{c} the Coulomb logarithm, or the logarithm of the ratio of the largest to closest distances used in the collision integral. In a PIC code, lnΛ^{c} = lnλ_{D}/X_{0}. Clearly, (T_{1}(t) + T_{2}(t))/2 is constant and equal to T_{∞}. It follows that if mass m_{1} and m_{2} are equal, the thermalization time is also constant and can be written (19)with Λ_{∞} = n [ϵ_{0}T_{∞}/(ne^{2})] ^{3/2} the plasma parameter based on the temperature T_{∞}. It also follows that the temperatures vary exponentially as T_{1} = T_{∞} − 0.5 [T_{2}(0) − T_{1}(0)] exp { − 2t/t_{th} } , and similarly for T_{2}.
Here we do not find the temperature curves to be strictly exponential, mainly because there are four species. The cold electrons interact with the hot electrons on a timescale t_{0}, but also with the hot ions on a timescale m_{i}/m_{e}t_{0} = 25t_{0}. The cold ions are heated by interactions with the hot ions on a timescale (m_{i}/m_{e})^{1/2}t_{0} = 5t_{0}, and by interactions with the hot electrons on a timescale m_{i}/m_{e}t_{0} = 25t_{0}. Since the cold ions are heated more slowly than the cold electrons, a temperature difference between these two components appears and they also start heating or cooling each other. Nevertheless, given the separation of scales we expect a measure of the slope around t = 0 to reflect the electronelectron or ionelectron thermalization times when measured on the electron or ion temperature curves, respectively.
Results are shown in Fig. 10 (bottom). We see that the relation t_{th} ∝ ρ_{sp} is roughly correct for both electrons and ions. We also underline the difference with a real plasma, where t_{th}/T_{pe} ~ Λ reaches 10^{10} or more, while it is on the order of Λ_{p} ≤ 10^{4} in PIC simulations.
5.3. Field fluctuation level dependence: coarsegraining and finite superparticle size
The previous section showed that the behavior of PIC plasma quantities can be guessed by the substitution Λ → Λ_{p}. While it is true for orders of magnitude estimates, this recipe is, however, not exact, and coarsegraining dependent quantities generally follow other relations than their real counterparts with respect to physical parameters (temperature, Debye length, plasma parameter, etc.). The main reason for this is that the finite volume of the superparticles implies a cutoff of the physical processes at smaller scales, an effect that becomes even more important when the superparticle size is close to the Debye length: typically ranges between one (or less) and ten in simulations.
This section illustrates this double dependence (Λ → Λ_{p} and superparticle size) with a detailed study of the level of electric field fluctuations ε in a PIC thermal plasma.
In a real plasma in thermal equilibrium, it is given by (Callen 2006, Sect. 1.1) (20)where the symbol ⟨ · ⟩ denotes an average over space.
Dieckmann et al. (2004) studied the spectrum of thermal fluctuations in a PIC plasma, but without investigating their levels. Hockney (1971, see also Hockney & Eastwood 1988 measured ratios like ε in a series of 2D simulations of thermal plasmas, and found a good agreement with the empirical formula , where is the superparticle geometrical size in number of cells. Its algorithm was two dimensional, electrostatic, and based on the integration of the Poisson equation.
We perform these simulations with our 3D electromagnetic code and measure the level of energy in the electric field. We use thermal velocities from 0.04c to 0.10c, ρ_{sp} from 2 to 500, and n_{x} from 10 to 128. It results in from 0.4 to 12.8, and in from 0.1 to 75 000. The fluctuation levels do not depend on the timestep (which varies from n_{t} = 2000 down to close to the Courant limit Δt ~ X_{0}/c) nor on box size (which is always bigger than n_{x}). We use a pair plasma, but increasing the mass of the ions would only multiply the fluctuation levels by a constant factor.
The results are summarized in Fig. 11: ε is found to be proportional to , and not exactly to 1/Λ_{p}. To explain this, we generalize the computation of Hockney to three dimensions.
For a plasma in thermal equilibrium, the energy in the electric field at location (x,t) can be evaluated by adding the electric field produced at x by charges at location x_{0} and having a velocity v_{0}, E_{x0,v0}(x,t), (21)where ⟨ · ⟩ means an ensemble average (which coincides with a spatial average); E_{x0,v0}(x,t) is a generalization of the Debye electric field for moving particles (Nicholson 1983, Chap. 9); and Eq. (21) can be evaluated for a plasma of finitesized particles as (Hockney 1971; Birsdall & Langdon 1985) (22)where S(ka) is the Fourier transform of the shape of the superparticles and a the characteristic size of the superparticles (in our case a ~ X_{0}); S(ka) tends to 1 as k^{1} ≫ a.
Fig. 11 Top: field energy levels as a function of . The colorbar is . Blue points at low f have an underresolved Debye length that could explain the mismatch with 1/f. Bottom: field energy levels versus . We clearly see the mismatch between 1/Λ_{p} and the results, even if the trend is correct. The large scatter is a hint that Λ is not a relevant parameter to describe field fluctuations. Top and bottom: each point is the result from a simulation. The field energy levels are measured as the energy in the x electric field, (with α from Eq. (A.6)), divided by the kinetic energy of the superparticles, ∑ _{sp}(γ_{sp} − 1). 
Equation (22) cannot be used as such for a real plasma of point particles (S(ka) = 1) because it includes the electric field at arbitrarily small distances from the charge, which has an infinite energy. It leads to Eq. (20) only if a truncation at small distances is performed, for example k < (αλ_{D})^{1} with α any constant: ⟨ ϵ_{0}E^{2}/2 ⟩ in Eq. (20) is then the energy in the electric field for wavelengths larger than αλ_{D}. Alternatively and to avoid a cutting procedure, we note that the electric field in Eq. (20) can be taken as the total field produced by the particles to maintain the screening Debye clouds (the polarization electric field; see Callen 2006, Sect. 1.1).
In the case of a PIC plasma all processes at scales below the grid size a = X_{0} are ignored, so that the upper bound of the integral is k_{max} = a^{1} and there is no small scale divergence. Since S(ka) ~ 1 for k ≪ a^{1}, and given that the integration stops at k = a^{1}, we will assume that S(ka) = 1.
Changing to spherical coordinates, using u = λ_{D}k and , we arrive at (23)A primitive of the integral is u − arctan(u). We use k_{max} = 1/X_{0} or , while the largest wavelength is given by the simulation domain size and verifies u_{min} ≪ u_{max}. Consequently, we obtain (24)in good agreement with the simulations (Fig. 11, top panel), except for the constant factors 18π^{2} ~ 178 (which can be attributed to an approximate choice of u_{max}).
The two limits are interesting. For a very high resolution, , the field energy decreases as , which is nontrivial and different from what is expected in a real plasma where it decreases as . The empirical formula of Hockney (1971), generalized to 3D, would also predict in the high resolution limit. However, our experiments clearly preclude this dependence, and are compatible with Eq. (24) (see Fig. 11). The presence of a finite superparticle volume and of the grid is retained in our calculation only in the upper bound of the integral: physically speaking, physical processes with k^{1} < X_{0} are smoothed out. This explains the difference between Eq. (24) and that of a real plasma.
For low resolutions, , an expansion of the arctangent shows that the field energy behaves as 1/(54π ρ_{sp}), which is finite and independent of . In a real plasma, ε would go to zero as the screening distance vanishes. That this is not the case here indicates that the screening distance does not vanish, because the finite size of the superparticles also plays the role of the screening mechanism.
5.4. Comparing the PIC and VlasovMaxwell models
We now highlight some differences between PIC models and kinetic models based on Liouville or Klimontovich equations. To do so, we recall how the VlasovMaxwell system is derived from these formalisms.
A plasma is constituted of many charged particles in mutual electromagnetic interaction. Under relativistic conditions, a plasma microstate is fully characterized by the positions and velocities of the N particles and by the value of the fields at all space points (the fields must be treated as independent from the particles because of retarded interactions), plus the necessary boundary conditions. Within the frame of classical electrodynamics, the time evolution of a microstate is described by Maxwell equations and by the equations of motion for the particles under the action of the Lorentz force: (25)\arraycolsep1.75ptHere, c is the speed of light, e_{m} the microscopic electric field; b_{m} = cB_{m} is c times the microscopic magnetic field, q_{j}, m_{j}, γ_{j}, v_{j}, and x_{j} the charge, mass, Lorentz factor, velocity, and position of the particle number j. We also define its momentum p_{j} = m_{j}γ_{j}v_{j}.
At a given time t, a microstate is represented by a point in the 6Ndimensional phasespace, and by the fields. One can then consider the collection of microstates having the same macroscopic properties (which can depend on what one is looking for), place them as points in the 6Ndimensional phasespace, and define the Nparticle distribution function as the number density, at a given time, of these microstates in the 6Ndimensional phasespace. Given that the number of microstates in phasespace is chosen as a continuum, f_{N} is a smooth function (Klimontovich 1982; Nicholson 1983). It defines a macrostate, i.e., an ensemble average of a collection of compatible microstates. The dynamic evolution is then obtained by the Liouville equation, which states that the number of microstates is conserved (26)supplemented by Maxwell equations for the microscopic fields e_{m} and b_{m}with for sources the particles of the corresponding microstate.
When there is no background magnetic field, and when the particles are nonrelativistic, the magnetic field remains negligible toward the electric field which can be computed directly from the particle positions at time t with the use of the electrostatic potential V(r) = q^{2}/(4πϵ_{0}r) (the Coulomb model). The Liouville equation can then be transformed to the infinite BBGKY hierarchy. The first BBGKY level involves the oneparticle distribution function f_{1}(t,w) (with w = (x,p)) and the twopoint distribution function f_{2} via the correlation function g_{2}(t,w_{1},w_{2}) = f_{2}(t,w_{1},w_{2}) − f_{1}(t,w_{1})f_{1}(t,w_{2}): Here C is an integral operator, vanishing with g_{2}. The quantity is the mean potential due to the particle distribution f_{1}. Since f_{1}(w) is the probability of finding a particle near w independently of the positions of all others, this potential does not include shortrange correlations, but only longrange collective effects. It is a macroscopic quantity, just as the fields entering into the VlasovMaxwell system (Eq. (28)).
Approximations can then be made to truncate the BBGKY hierarchy. Evaluations of the righthand side of Eq. (27a) can lead, depending on the hypothesis made, to Boltzmann, Landau, LenardBalescu, or more refined kinetic equations. For a fully ionized plasma, the relevant approximation parameter is the number of particles per Debye sphere, or plasma parameter Λ. For large Λ, collisions and correlations between small numbers of particles are negligible toward interactions of particles with the fields collectively generated by all others. Keeping only these collective interactions is equivalent to a truncation of BBGKY hierarchy at its lowest level, i.e., g_{2} = 0, and leads to the VlasovMaxwell system. Generalized to relativistic plasmas^{5} and to two species, it reads (28)\arraycolsep1.75ptwith f_{1,s}(t,x,p) the oneparticle distribution function, for electrons (s = e) or ions (s = i), defined from f_{N}. It is also a smooth function, and the VlasovMaxwell system describes the evolution of a continuous fluid in the sixdimensional phasespace, where information on the individual nature of the particles has been smoothed out. In particular within this description, the plasma parameter Λ is infinite, the fluctuation and collisioninduced thermalization time is infinite, and the level of electric field fluctuations ε is zero. We note that these are analytical properties. Numerical solutions of the VlasovMaxwellsystem will also not strictly recover the collisionless behavior of the plasma. For instance, numerical diffusion arising necessarily from the discretization of Eq. (28) will also lead to a finite thermalization time. We are, however, not aware of a comprehensive study of such effects for algorithms solving the discretized VlasovMaxwell system.
In contrast, the models underlying PIC simulations follow a different path, illustrated in Fig. 12. It consists in following the time evolution of the microstate constituted by N/p superparticles and the fields, Eqs. (2a) − (2e), with p reaching 10^{10} or more. One of the consequences is that the PIC plasma has a plasma parameter Λ_{p} = Λ/p (from Eq. (14)) far smaller than that of the real plasma, so that it includes relatively large correlation and noise levels.
Fig. 12 Different plasma models. Dashed arrows are transitions showing nontrivial effects. Coarsegraining and discretization are discussed in Sects. 5.1 to 5.3. 
5.5. Higherorder effects of coarsegraining
We have highlighted that the coarsegraining step, the description of a real plasma of N particles by a PIC plasma of N/p superparticles, each containing p real particles, involves a reduction of the plasma parameter Λ by a factor p. This is the main effect of coarsegraining. Higherorder effects arise because after coarsegraining, the model ignores the internal dynamics and correlations of the particles contained within a superparticle.
This can be seen by writing explicitly the grouping: we label the particles either by w_{n}, n = 1...N (with w = (x,p)), or by w_{ij} with i = 1...N/p representing the group number, and j = 1...p the particle number within this group. We denote by the position and velocity of the centerofmass of the group number i. The Nparticle distribution function f_{N} of the real plasma can then be written formally as (29)This equation introduces , the analog of f_{N} but for the centerofmass of the particle groups (i.e., of the superparticles); f_{sp,i}(t,w_{i1},...,w_{ip}), the distribution function of the particles contained within a group, which represents the dynamics and correlations between particles of the same group; and g_{corr}(t,w_{1},...,w_{N}), the correlations ignored by writing f_{N} = f_{N/p} × f_{sp,1}...f_{sp,N/p}, i.e., the correlations between particles of different groups. The PIC approximation then consists in setting g_{corr} = 0 and f_{sp,i}(t,w_{i1},...,w_{ip}) = const.
A complete understanding of these approximations would require developing a BBGKY hierarchy from Eq. (29) and making explicit the electric and magnetic field contributions from the particle groups. This is a complex task. We can, however, stress important consequences of the assumption f_{sp,i}(t,w_{i1},...,w_{ip}) = const.:

The superparticles are assumed incompressible. Thecompressibility due to particle motion within a particle group isthus absent from the coarsegrained plasma.

The velocity dispersion of the particles within a group is ignored in the coarsegrained plasma. The kinetic pressure resulting from this dispersion is thus also absent.

The electric fields present in the PIC plasma are computed from the superparticles. This is equivalent to saying that they are computed by taking into account only the monopole distribution of charge created by the internal arrangement of the particles within a group, with higherorder multipole terms neglected. The same holds for the magnetic field.
6. Discussion and conclusion
6.1. Code validation
We have presented our particleincell code AparT (Sect. 2 and Appendix A), and studied several validation tests.
Computation of the spectra from a magnetized or unmagnetized thermal plasma at rest (Sect. 3.1) has proven very accurate, with the plasma pulsation and the right and left cutoff pulsations precisely recovered (Figs. 1 and 2), even in cases where the Debye length and the Larmor radius are not resolved. This proves that the description of individual particle motions in a constant magnetic field and of collective particle dynamics is accurate, and robust with respect to numerical resolution. In particular we showed that Larmor orbits in a constant magnetic field are well described provided that the cyclotron pulsation is well resolved, independently of the grid size. We found, however, a numerical instability with abnormal behavior in the energy curves and high noise levels for underresolved Debye length or underresolved Larmor radius, so that this parameter range should be avoided.
Simulations of the filamentation instability (Sect. 3.2) showed a good agreement with linear cold theories, provided that the growth rates are computed from the temporal evolution of the Fourier modes. Discrepancies with theory then range from 5% to 13%, and can be explained by the wavenumber dependence of the growth rates. On the other hand, the linear growth rates derived from the total energy curves, or equivalently from the sum of all Fourier modes, present larger discrepancies with linear theory, ranging between 12% to 61%. This is explained by the high level of fluctuations in PIC codes that prevent the fastest growing modes to dominate the total energy before the end of the linear phase of the instability.
Simulations of the tearing instability in a relativistic pair plasma gave linear growth rates within 8% of those found by Pétri & Kirk (2007) with an analytical linear VlasovMaxwell solution, a result not varying significantly when changing the numerical resolution (Sect. 3.4, Fig. 8, and Table 4). This example, where the shape of the velocity distribution is a key feature, is thus in agreement with the VlasovMaxwell description. It also validates the new method used to load the relativistic (in both temperature and bulk velocity) MaxwellJüttner distribution that we present in Sect. 4 and Appendix D, as well as the general relations used for the relativistic Harris equilibrium (Appendix C). We note that in this case the total energy curves can be used for evaluation of the linear growth rates because the linear phase spans several orders of magnitude in field intensity, so that the fastest growing modes have enough time to dominate the energy.
All in all, these tests show that our code AparT is a sound basis for future explorations. They also serve as a base to explore important questions regarding the nature of a PIC plasma, that we summarize in the next section.
6.2. PIC and real plasmas
The widespread use of PIC codes for studying plasmas out of equilibrium calls for a deep understanding of the PIC model, and of its relation with a real plasma and with the VlasovMaxwell description. Section 5 attempted to provide some explanations.
We have seen that the PIC model lies on two building blocks. The first stems from the capability of computers to handle only up to ~ 10^{10} particles, while real plasmas contain from 10^{4} to 10^{20} particles per Debye sphere. This means that a coarsegraining step must be used, whereby of the order of p ~ 10^{10} real particles are represented by a single computer superparticle. The second step is field storage on a grid with its subsequent finite superparticle size.
We have introduced the notion of coarsegraining dependent quantities, i.e., physical quantities depending on p. The prototype of such quantities is the plasma parameter Λ, that behaves as Λ_{p} ∝ 1/p. This vast reduction of Λ induces higher noise levels and correlations, but we have again seen that it does not threaten plasma and collisionless behavior as long as Λ_{p} remains above unity. All coarsegraining dependent quantities can be expressed as the product of a fluid quantity (which is coarsegraining independent) and a parameter expressing a number of particles per fluid volume. Examples of such parameters include , n(c/ω_{pe})^{3}, n(c/ω_{pi})^{3}, ... Their behavior in the PIC plasma can be guessed by taking into account the reduction by a factor p of the number of particles, leading for example to the substitution Λ → Λ_{p} = Λ/p in the relevant analytical expressions. We checked this for the collision and fluctuation induced thermalization time (Sect. 5.2), which is indeed proportional to the number of computer superparticles per cell; the lower the number the shorter the thermalization time. Bret et al. (2013) similarly reduce the parameter n(c/ω_{pe})^{3} by a factor p when applying their theory for the magnetic fluctuation level in a drifting plasma to their PIC simulations. However, the substitution Λ → Λ_{p} is strictly valid only for pointsize particles, and the large finite size of the superparticles, which reaches a fraction of a Debye length, suppresses interactions and fluctuations at shorter wavelengths and modifies these scalings. We have detailed how this works for the electric field fluctuation level in a thermal plasma in Sect. 5.3.
We stress that the reduction of the collision and fluctuation induced thermalization time and of other related timescales (e.g., the slowingdown time of fast particles), by 10 or more orders of magnitude, can have important consequences for the relevance of simulations: one has to insure that collisionless kinetic processes remain more efficient than the artificially enhanced collisional and fluctuation induced PIC effects. Similarly, we have seen in Sect 3.2 that the high level of fluctuations alter the linear spectrum of instabilities by preventing the fastest growing modes to dominate the total energy.
A more subtle effect of coarsegraining is due to the loss of the dynamics of the p particles represented by each superparticle. We intuitively expect that it will lead to the overall loss of compressibility due to superparticle incompressibility, of the contribution to kinetic pressure of the particle velocity spreading within a superparticle, and of the multipole contribution to the electric and magnetic fields created by the distribution of particles within a superparticle (see Sect. 5.5). The relevance of these missing effects remains unclear.
6.3. PIC and VlasovMaxwell plasmas
We have highlighted in Sect. 5.4 that a PIC algorithm simulates a plasma of finitesized charges in their selffields, and does not strictly solve the VlasovMaxwell system. Using the Vlasov equation assumes that the plasma is represented in phase space by a continuous fluid. In this limit of an infinite number of particles, the plasma parameter Λ and the thermalization time are infinite, and the collision frequency and the thermal field fluctuation levels are zero. Using a PIC algorithm amounts to dividing the continuous phase space fluid into discrete elements, and to following their orbits. In this sense, one can say that we integrate the characteristics of the Vlasov equation. However, the newly introduced graininess (which is far higher than that of the original plasma) implies the presence of binary collisions and of correlations between superparticles that is not easy to evaluate, in part because they are reduced by the finite size of the superparticles and the subsequent vanishing of the twopoint force at short distances (see Sect. 5.1). The intricate dependence of Eq. (24) is a hint to this complexity. We note that Birsdall & Langdon (1985, Chap. 12) have derived a generalization of the BalescuGuernseyLenard kinetic equation that includes the use of a grid (and thus of finite sized superparticles), and of the discretization in space and time of the equations. The correlations just mentioned are partly present in this equation, but difficult to extract.
These differences between PIC and VlasovMaxwell plasmas are especially enhanced in the linear phase of instabilities. We see two main points. The first is that nonlinear effects absent from the linear theory, and possibly enhanced by the high noise level of the simulation (Sect. 5.3), may have visible consequences (Birsdall & Langdon 1985, Sect. 13.6; Dieckmann et al. 2006). This example is reported by Daughton (2002) in the context of the drift kink instability of a current sheet: the instability is found to grow faster than predicted by the linear Vlasov theory because the early development of another instability quickly produces nonlinear effects. Bret et al. (2010, Fig. 23) also report significant early nonlinear behavior in counterstreaming situations. We have also reported the presence of field components due to nonlinear effects in Fig. 4.
The second point is that the high level of fluctuations delays the dominance of the fastest growing Fourier modes over the sum of other modes. The consequence is that effective linear growth rates measured from total energy curves appear slower than the growth rates of the fastest modes. This is even more important in instabilities where the linear phase is short, and explains the differences between the effective growth rates and the linear cold theory of the counterstreaming instability measured in Sect. 3.2, with discrepancies reaching 60% or more. It may also explain the differences between theory and measured growth rates of Cottrill et al. (2008), Dieckmann et al. (2006), Haugboelle et al. (2013) for the counterstreaming instability. On the other hand, the differences can be small if the linear phase lasts long enough for the fastest mode to dominate the energy, as is the case for the relativistic tearing instability in Sect. 3.4 or for the Weibel instability of Markidis et al. (2010).
6.4. Modeling astrophysical plasmas
In the light of what has been said so far one may wonder what this all implies for the modeling of astrophysical plasmas. We attempt to give some answers in the following. We have shown that the PIC description of an astrophysical plasma bears some risk because the plasma parameter Λ is always underestimated, leading to systematic errors in the evaluation of important parameters of the plasma such as the collision and fluctuation induced thermalization time. We stress that this does not lessen the important role of PIC algorithms for deepening our understanding of plasma physics. In particular they also have their virtues, for instance the consideration of certain correlations and of direct particle encounters (with the restrictions discussed in Sect. 5). They provide an accurate description of collisionless kinetic processes such as instabilities, and of the induced turbulence and eventual associated thermalization relevant to collisionless environments.
The VlasovMaxwell equations perfectly describe a plasma free of collisions and fluctuations. However, their discretization will again introduce different plasma characteristics. A thorough discussion of these effects is still missing. On the other hand, a collisionless description of astrophysical plasmas is not always correct. On larger spatial and temporal scales, the description of flows and the propagation of nonthermal particles must include collisions to a certain degree. In this regime other models, and in particular FokkerPlanck models, have been shown to give a good description of the plasma and have provided significant results. However, FokkerPlanck models have their own drawbacks, notably that they are local and use dragging and diffusion coefficients not selfconsistently derived.
These discrepancies between real, PIC, and VlasovMaxwell plasmas are complex, and it needs to be discussed in further detail under what circumstances which model and which numerical realization comes closest to a real plasma. In the long term, it may be justified to use models including correlations in a more systematic way, for example the Landau or the LenardBalescu equation on the theoretical side, and P^{3}M algorithms on the numerical side (that include shortrange particleparticle interactions Hockney & Eastwood 1988). Both approaches should then be faced with results from wellcontrolled collisionless plasma experiments, which are presently in their infancy (see, e.g., Grosskopf et al. 2013).
We note that after the submission of our article, Swisdak (2013) published a similar method.
In a fully ionized plasma, all coarsegraining dependent quantities can be expressed as the product of a fluid quantity (which is coarsegraining independent) and a parameter expressing a number of particles per fluid volume. Examples of these parameters include , n(c/ω_{pe})^{3}, n(c/ω_{pi})^{3}, ...
The relativistic VlasovMaxwell system is usually derived by using the Klimontovich formalism, not the Liouville formalism (see, e.g., Nicholson 1983; Klimontovich 1982).
Acknowledgments
We would like to thank Jérôme Pétri for information and discussions about linear growth rates. We also thank the anonymous referee for comments that allowed to significantly improve the manuscript. We acknowledge support from the French Stellar Astrophysics Program PNPS. Some of the simulations have been performed at pôle scientifique de modélisation numérique, PSMN, at ENSLyon, whose staff we thank for their steady technical support. Larger simulations were performed using HPC ressources from GENCI (Grand Équipement National de Calcul Intensif) at CINES and CCRT, under the allocation x2012046960.
References
 Baumann, G., & Nordlund, Å. 2012, ApJ, 759, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Benedetti, C., Sgattoni, A., Turchetti, G., & Londrillo, P. 2008, IEEE Trans. Plasma Sci., 36, 1790 [NASA ADS] [CrossRef] [Google Scholar]
 Birdsall, C. K. 1967, NASA Special Publ., 153, 375 [NASA ADS] [Google Scholar]
 Birdsall, C. K. 1999, in Dynamical Systems, Plasmas and Gravitation, eds. P. G. L. Leach, S. E. Bouquet, J.L. Rouet, & E. Fijalkow, Lect. Notes Phys. (Berlin: Springer Verlag), 518, 383 [Google Scholar]
 Birsdall, C., & Langdon, A. 1985, Plasma Physics via Computer Simulation (NewYork: McGrawHill) [Google Scholar]
 Bret, A., Gremillet, L., & Dieckmann, M. E. 2010, Phys. Plasmas, 17, 120501 [NASA ADS] [CrossRef] [Google Scholar]
 Bret, A., Stockem, A., Fiuza, F., et al. 2013, Phys. Plasmas, 20, 042102 [NASA ADS] [CrossRef] [Google Scholar]
 Buneman, O. 1959, Phys. Rev., 115, 503 [Google Scholar]
 Buneman, O. 1976, Comp. Phys. Comm., 12, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Burger, P. 1965, J. Appl. Phys., 36, 1938 [NASA ADS] [CrossRef] [Google Scholar]
 Bykov, A. M., & Treumann, R. A. 2011, A&ARv, 19, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Callen, J. 2006, Fundamentals of Plasma Physics (Madison: University of Wisconsin) [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2012, ApJ, 754, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Cerutti, B., Werner, G. R., Uzdensky, D. A., & Begelman, M. C. 2013, ApJ, 770, 147 [NASA ADS] [CrossRef] [Google Scholar]
 ChacónAcosta, G., Dagdug, L., & MoralesTécotl, H. A. 2010, Phys. Rev. E, 81, 021126 [NASA ADS] [CrossRef] [Google Scholar]
 Childs, H., Brugger, E., Whitlock, B., et al. 2012, in High Performance VisualizationEnabling ExtremeScale Scientific Insight, 357 [Google Scholar]
 Cottrill, L. A., Langdon, A. B., Lasinski, B. F., et al. 2008, Phys. Plasmas, 15, 082108 [NASA ADS] [CrossRef] [Google Scholar]
 Cubero, D., CasadoPascual, J., Dunkel, J., Talkner, P., & Hänggi, P. 2007, Phys. Rev. Lett., 99, 170601 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Daughton, W. 2002, Phys. Plasmas, 9, 3668 [NASA ADS] [CrossRef] [Google Scholar]
 Daughton, W., Scudder, J., & Karimabadi, H. 2006, Phys. Plasmas, 13, 072101 [NASA ADS] [CrossRef] [Google Scholar]
 Dawson, J. 1962, Phys. Fluids, 5, 445 [NASA ADS] [CrossRef] [Google Scholar]
 Dieckmann, M. E. 2009, Plasma Phys. Control. Fusion, 51, 124042 [NASA ADS] [CrossRef] [Google Scholar]
 Dieckmann, M. E., Ynnerman, A., Chapman, S. C., Rowlands, G., & Andersson, N. 2004, Phys. Scr., 69, 456 [NASA ADS] [CrossRef] [Google Scholar]
 Dieckmann, M. E., Frederiksen, J. T., Bret, A., & Shukla, P. K. 2006, Phys. Plasmas, 13, 112110 [NASA ADS] [CrossRef] [Google Scholar]
 DLMF 2013, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.0.6 of 20130506 [Google Scholar]
 Drake, J. F., Swisdak, M., Che, H., & Shay, M. A. 2006, Nature, 443, 553 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Dunkel, J., & Hänggi, P. 2009, Phys. Rep., 471, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Eldridge, O. C., & Feix, M. 1962, Phys. Fluids, 5, 1076 [Google Scholar]
 Esirkepov, T. Z. 2001, Comput. Phys. Comm., 135, 144 [Google Scholar]
 Fitzpatrick, R. 2011, The Phys. Plasmas [Google Scholar]
 Fonseca, R. A., Silva, L. O., Tsung, F. S., et al. 2002, Lect. Notes Comput. Sci., 2331, 342 [CrossRef] [Google Scholar]
 Fonseca, R. A., Martins, S. F., Silva, L. O., et al. 2008, Plasma Phys. Control. Fusion, 50, 124034 [NASA ADS] [CrossRef] [Google Scholar]
 Fujimoto, K., & Machida, S. 2006, J. Comput. Phys., 214, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Greenwood, A. D., Cartwright, K. L., Luginsland, J. W., & Baca, E. A. 2004, provided by the SAO/NASA Astrophysics Data System, J. Comput. Phys., 201, 665, [NASA ADS] [CrossRef] [Google Scholar]
 Grosskopf, M., Drake, R., Kuranz, C., et al. 2013, High Energy Density Physics, 9, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Grote, D. P., Friedman, A., Vay, J.L., & Haber, I. 2005, in Electron Cyclotron Resonance Ion Sources, ed. M. Leitner, Am. Ins. Phys. Conf. Ser., 749, 55 [Google Scholar]
 Hartree, D. R. 1950, Appl. Sci. Res., 1, 379 [CrossRef] [Google Scholar]
 Haugboelle, T., Fredriksen, J. T., & Nordlund, A. 2013, Phys. Plasmas, 20, 06290406290424 [NASA ADS] [CrossRef] [Google Scholar]
 Hededal, C. 2005, Ph.D. Thesis, Niels Bohr Institute [Google Scholar]
 Hesse, M., Kuznetsova, M., & Birn, J. 2001, J. Geophys. Res., 106, 29831 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W. 1966, Phys. Fluids, 9, 1826 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W. 1971, J. Comput. Phys., 8, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1988, Computer simulation using particles (Institute of Physics Publishing) [Google Scholar]
 Hoh, F. C. 1966, Phys. Fluids, 9, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, J., & Ma, Z.W. 2008, Chin. Phys. Lett., 25, 1764 [NASA ADS] [CrossRef] [Google Scholar]
 Jaroschek, C. H., & Hoshino, M. 2009, Phys. Rev. Lett., 103, 075002 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Jaroschek, C. H., Lesch, H., & Treumann, R. A. 2005, ApJ, 618, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Jüttner, F. 1911, Ann. Phys., 339, 856 [Google Scholar]
 Kirk, J. G., & Skjæraasen, O. 2003, ApJ, 591, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Klimas, A., Hesse, M., & Zenitani, S. 2008, Phys. Plasmas, 15, 082102 [NASA ADS] [CrossRef] [Google Scholar]
 Klimontovich, I. L. 1982, Oxford Pergamon Press International Series on Natural Philosophy, 105 [Google Scholar]
 Langdon, A. B. 1970, J. Comput. Phys., 6, 247 [Google Scholar]
 Markidis, S., Lapenta, G., & Rizwanuddin 2010, Math. Comput. Simul., 80, 1509 [CrossRef] [Google Scholar]
 Markidis, S., Henri, P., Lapenta, G., et al. 2012, Nonlinear Processes in Geophysics, 19, 145 [Google Scholar]
 Matsumoto, H., & Omura, Y. 1993, Computer Space Plasma Physics: Simulation Techniques and Software [Google Scholar]
 Medvedev, M. V., Fiore, M., Fonseca, R. A., Silva, L. O., & Mori, W. B. 2005, ApJ, 618, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Messmer, P. 2001, Ph.D. Thesis, Institute of Astronomy, ETH Zürich, Switzerland, No. 14412 [Google Scholar]
 Messmer, P. 2002, A&A, 382, 301 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michno, M. J., & Schlickeiser, R. 2010, ApJ, 714, 868 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Nicholson, D. R. 1983, Introduction to Plasma Theory, eds. John Wiley & Sons [Google Scholar]
 Nieter, C., & Cary, J. R. 2004, J. Comput. Phys., 196, 448 [Google Scholar]
 Nishikawa, K.I., Mizuno, Y., Fishman, G. J., & Hardee, P. 2008, Int. J. Mod. Phys. D, 17, 1761 [Google Scholar]
 Nishikawa, K.I., Niemiec, J., Medvedev, M., et al. 2011, in AIP Conf. Ser. 1366, eds. V. Florinski, J. Heerikhuisen, G. P. Zank, & D. L. Gallagher, 163 [Google Scholar]
 Paesold, G., Blackman, E. G., & Messmer, P. 2005, Plasma Phys. Control. Fusion, 47, 1925 [NASA ADS] [CrossRef] [Google Scholar]
 Pétri, J., & Kirk, J. G. 2007, Plasma Phys. Control. Fusion, 49, 1885 [NASA ADS] [CrossRef] [Google Scholar]
 Pétri, J., & Lyubarsky, Y. 2007, A&A, 473, 683 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pomraning, G. C. 1973, The equations of radiation hydrodynamics (Pergamon Press) [Google Scholar]
 Silva, L. O., Fonseca, R. A., Tonge, J. W., et al. 2003, ApJ, 596, L121 [Google Scholar]
 Sironi, L., & Spitkovsky, A. 2011, ApJ, 726, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Spitkovsky, A. 2005, in Astrophysical Sources of High Energy Particles and Radiation, eds. T. Bulik, B. Rudak, & G. Madejski, Am. Inst. Phys. Conf. Ser., 801, 345 [Google Scholar]
 Spitzer, L. 1965, Physics of fully ionized gases (New York: Interscience Publication) [Google Scholar]
 Swisdak, M. 2013, Phys. Plasmas, 20, 0621100621104 [NASA ADS] [CrossRef] [Google Scholar]
 Tien, P. K., & Moshman, J. 1956, J. Appl. Phys., 27, 1067 [NASA ADS] [CrossRef] [Google Scholar]
 Treumann, R. A., Nakamura, R., & Baumjohann, W. 2011, Ann. Geophys., 29, 1259 [NASA ADS] [CrossRef] [Google Scholar]
 Trier Frederiksen, J., Haugbølle, T., Medvedev, M. V., & Nordlund, Å. 2010, ApJ, 722, L114 [NASA ADS] [CrossRef] [Google Scholar]
 Vay, J.L. 2008, Phys. Plasmas, 15, 056701 [NASA ADS] [CrossRef] [Google Scholar]
 Whitlock, B., Favre, J. M., & Meredith, J. 2011, In Proc. of the 11th Eurographics Symposium on Parallel Graphics and Visualization (EGPGV’11), 101 [Google Scholar]
 Yu, S. P., Kooyers, G. P., & Buneman, O. 1965, J. Appl. Phys., 36, 2550 [NASA ADS] [CrossRef] [Google Scholar]
 Zenitani, S., & Hoshino, M. 2008, ApJ, 677, 530 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Numerical implementation
This Appendix is the direct continuation of Sect. 2.
A.1. Computation of the current
Before going through the normalization of Eqs. (2a) − (2e), we have to understand how the current is computed. As said in Sect. 2.2.2, interpolation of particle quantities to grid nodes is done by attributing to the superparticles a finite volume . We consider a superparticle, and the cell that contains its center. At time t, the superparticle occupies a volume V_{t} of the cell. The charge in the cell is given by Q_{cell}(t) = q_{sp} × (V_{t}/V_{sp}). The charge continuity equation then gives (A.1)The superparticle volume necessarily intersects three faces of the cell that contains its center: one of perpendicular along x, one along y, and one along z. Consequently, the motion of this superparticle will create a current through these three faces, and we can write (A.2)We have to know which part of the volume variation V_{t + dt} − V_{t} is attributed to each part of the current. The displacement of the superparticle between t and t + dt is denoted (Δx,Δy,Δz). The volume variation depends on this displacement, and the part of it proportional to Δx is attributed to j_{x}, and similarly for the y and z components. More specifically, we can write V_{t + dt} − V_{t} = A_{x}Δx + A_{y}Δy + A_{z}Δz. The areas A_{i} can be evaluated with some geometry (see Sect. A.3.4). Then (and similarly for y and z): (A.3)This way of computing the current ensures that the discrete charge conservation equation is fulfilled, and justifies the advection of the divergence of the fields.
A.2. Normalization
The problem is formulated with as many equations as variables (Eqs. (2a) − (2e), variables e, b, j, x_{sp}, and v_{sp} for each superparticle), and it is possible to normalize the equations in a way independent of any physical quantity. We denote normalized quantities with a tilde.
We choose to normalize lengths by X_{0}. Consequently, the normalized stepsize is unity. Times are normalized by T_{0} = X_{0}/c, and velocities are then naturally normalized to X_{0}/T_{0} = c. For the fields E and B, we use e = E and b = cB. These last two quantities are normalized by e_{0} = b_{0} = m_{e}c^{2}/(eX_{0}). With this, Eqs. (2a), (2b), and (2c) transform into \arraycolsep1.75ptwith m_{s} = m_{e} or m_{i} and q_{s} = −e or e (e is positive).
For the current j, the algorithm computes the quantity A_{x}Δx/V_{sp}. Writing for example the x component of Eq. (2d) gives (A.5)Using Eq. (1), we can write (A.6)where sgn is the sign of the charge.
With this, all the equations are completely independent of any physical quantity related to the simulated problem, and depend only on space discretization (n_{x}) and particle coarsegraining (). Time discretization (n_{t}) will play a role in time integration.
As a final comment, we note that is a priori unrelated to the actual number of superparticles per cell during the simulation. One should, however, make these two values not too far apart, because the size and timesteps are a fraction n_{x} and n_{t} of the skin depth and plasma period of a density plasma. If, for example, the superparticle density in the simulation is twice , then n_{x} cells will now represent two skin depths of the density plasma, and the resolution will decrease. This must be kept in mind in simulations where high density contrasts appear.
A.3. Discrete version of the equations
In this section we drop the tilde over normalized quantities. We denote the time at which they are considered by a superscript and their spatial location on the grid by a subscript.
A.3.1. The main loop
The strategy is to use a leapfrog scheme. It has the advantages of being time centered and reversible, and second order in time and space.
Before the loop, b and v are known at time t − dt/2, and e and x are known at time t. This should also be true for the initial conditions, so that initially we integrate backward the velocities and the bfield by − dt/2. Injected particles (if any) should also be correctly staggered.
The structure of the main loop is the following:

1.
Half advance ofb^{t − dt/2} with ∇ ∧ e^{t}; b is now at time t.

2.
Update of v^{t − dt/2} with b^{t} and e^{t}; v is now at time t + dt/2.

3.
Update of x^{t} with v^{t + dt/2}; x is now at time t + dt.

4.
Half advance of b^{t} with ∇ ∧ e^{t}; b is now at time t + dt/2.

5.
Boundary for b.

6.
Full advance of e^{t} with ∇ ∧ b^{t + dt/2}; e is now at time t + dt.

7.
Boundary for e.

8.
Boundary for particles.

9.
Computation of the currents from v^{t + dt} and x^{t + dt/2}.

10.
Filtering of the currents.

11.
Boundary for the currents.

12.
Add currents to e^{t + dt}.
A.3.2. Integration of the fields
The fields are stored on the grid in a staggered way, with e at the center of the grid edges and b at the center of the grid faces (Fig. A.1). This is the socalled Yee lattice. It allows an easy integration of the fields, and is secondorder accurate in time and space (Birsdall & Langdon 1985, Sect. 15), (A.7)and similarly for the other components of b and for e.
To reduce the effects of Čerenkov emission (see Appendix B.1), we have also implemented a fourth order solver (Greenwood et al. 2004).
Fig. A.1 Grid and fields locations. 
A.3.3. Moving the particles
Integration of the equation of motion for the superparticles is done with the algorithm described by Birsdall & Langdon (1985, Sect. 15.4). It is a relativistic generalization of the leapfrog scheme, time centered, time reversible, and second order accurate. We note however that, as pointed out by Vay (2008), it can have shortcomings for ultrarelativistic particles. In short, Eqs. (A.4a) and (A.4b) are discretized as with u = γv and γ the Lorentz factor. Defining u^{ − } = u^{t − dt/2} + qm_{e}/(em)e^{t}dt/2 and u^{ + } = u^{t + dt/2} − qm_{e}/(em)e^{t}dt/2 and substituting into (A.8a) leads to (A.9)This equation is the classical rotation around a b field (Birsdall & Langdon 1985, Sect. 4.4; Hockney & Eastwood 1988, Sect. 4.7.1), and is solved via (A.10)with the rotation vector Ω = qm_{e}/(2emγ^{t})b^{t}. Finally, we use . In Eq. (A.8a), the fields have to be known at time t. It explains the need for the half advances of b in the general scheme.
The interpolation of the fields at particle positions is done via a trilinear interpolation. We denote by (i,j,k) the nodes of the main grid A (Fig. A.1), and we introduce a second grid B whose cell centers are on (i,j,k). Consider a superparticle at position x = i + δx, y = j + δy, z = k + δz. The superparticle is actually a charge cloud of volume equal to a cell, and this volume intersects the cell of the second grid with center (i,j,k) in a volume V_{i,j,k} = (1 − δx)(1 − δy)(1 − δz), the cell of the second grid with center (i + 1,j,k) in a volume V_{i + 1,j,k} = δx(1 − δy)(1 − δz), and so on. For a quantity f defined at grid points (i,j,k), the weight associated to f_{i,j,k} is V_{i,j,k}, the one associated to f_{i + 1,j,k} is V_{i + 1,j,k}, and so on for a total of 8 points.
However, neither e nor b are defined at grid points (i,j,k) (Fig. A.1), and they must be first interpolated at grid points before applying the above procedure. This is done for example with f_{i,j,k} = 0.5(e_{x,i − 1,j,k} + e_{x,i,j,k}) or f_{i,j,k} = 0.25(b_{z,i,j,k} + b_{z,i − 1,j,k} + b_{z,i − 1,j − 1,k} + b_{z,i,j − 1,k}). Details can be found in Matsumoto & Omura (1993); Messmer (2001).
We note that the superparticle shape used for interpolation of fields to particle position and for interpolation of the current to grid nodes is the same. This is required to avoid the existence of a selfforce on the superparticles and to conserve the total momentum (Birsdall & Langdon 1985, Sect. 8.6).
A.3.4. Computation of the current
The current j_{i,j,k} is defined at the same locations as e_{i,j,k}. For current deposition, we again consider the volumes occupied by the superparticle in the grid B cells. As the superparticle moves, these volumes vary. We denote by (i + δx,j + δy,k + δz) the position of the superparticle at t − dt, and we assume that it moves from (Δx,Δy,Δz) between t − dt and t.
Consider, for example, the volume of the superparticle in the cell of center (i,j,k). Its variation is given by dV = (1 − δx − Δx)(1 − δy − Δy)(1 − δz − Δz) − (1 − δx)(1 − δy)(1 − δz). Defining d_{x} = δx + Δx/2, c_{x} = 1 − d_{x}, and similarly for y and z, one finds (A.11)As explained in Appendix A.1, the part of (A.11) proportional to the displacement along x is attributed to j_{x  i,j,k}, and so on.
A similar treatment is done with the cells that intersect the superparticle volume. These are cells centered in (i + ϵ,j + η,k + ξ), with ϵ, η and ξ equal either to 0 or 1 (8 cells). For each of these cells, only the faces intersecting the superparticle volume are concerned, so that in total there are only 12 currents to update.
Currents can be smoothed before being added to e. This has the effect of reducing electromagnetic noise (see Appendix B.2). This is done in the following way for the current of cell (i,j,k): attribute a weight of 1 to cell (i,j,k); of 0.5 to cells (i ± 1,j,k), (i,j ± 1,k), and (i,j,k ± 1); of 0.25 to cells (i ± 1,j ± 1,k), (i ± 1,j,k ± 1), and (i,j ± 1,k ± 1); of 0.125 to cells (i ± 1,j ± 1,k ± 1); and normalize the sum of the weights to 1.
Fig. A.2 Simulation time duration versus number of cores, the spatial domain extending in proportion to the number of cores. The standard deviation corresponds to 4% of the mean value. 
A.4. Boundaries
Periodic and reflective boundaries are available. The latter simulate a perfect conductor at the domain boundary by imposing the correct values for the electric and magnetic fields (b = e = 0 inside the conductor, b_{normal} = e_{tangential} = 0 at the conductor surface), and by reflecting the particles.
A.5. Parallel efficiency
The code parallelization was performed and tested by Messmer (2001, Chap. 4). It uses Fortran 90 and MPI. The simulation domain is decomposed in subdomains of equal length along the z direction, and all the cells and particles of each subdomain are assigned to a processor. To minimize communications between processors, ghost cells for the fields are added to each subdomain. Communication between neighboring processors occurs at each step involving the boundaries: for particles leaving or entering the domain, for the fields, and for the currents.
The domain is currently decomposed along one direction only. This is relevant for simulations of collisionless shocks where the domain is elongated along the flow direction, or for 2D magnetic reconnection simulations where the presence of the overdense current sheet at the domain center would lead to load balancing issues if a 2D domain decomposition were used.
We have tested the efficiency of this implementation with simulations using 16 superparticles per cell and a domain size of 60 × 60 × (16n_{proc}), where the number of cores varied from n_{proc} = 16 to 256. The corresponding (weak) scaling results, shown in Fig. A.2, are satisfactory. Simulation times scatter with a standard deviation of 4% around a constant value. The scatter is probably due to the uncontrolled node geometry.
Appendix B: Numerical effects
We have said in Sect. 5 that passing from a real plasma to a PIC model implies a discretization of the equations. This step comes with numerical issues that have been largely studied by Birsdall & Langdon (1985) and Hockney & Eastwood (1988). We highlight part of their work here.
B.1. Local numerical effects

Stability of the electric part of the superparticle motion integrator used here requires thatΩΔt < 2, with Ω the pulsation of oscillation of the superparticles (usually the plasma pulsation). The magnetic part is unconditionally stable.

Courant condition for the stability of the field integrator in vacuum is .

The dispersion relation of electromagnetic waves in vacuum is modified by the grid. This modification depends on the angle of propagation with respect to the grid, and waves can have a phase velocity smaller than c (Greenwood et al. 2004). If superparticles with velocity close to c are present, they can overtake light waves and emit Čerenkov radiation. This results in the production of nonphysical fields. The situation can be improved with a higher order interpolation scheme for the fields.
B.2. Global numerical effects
By considering the algorithm as a whole, Birsdall and Langdon were able to identify numerical effects not predicted by the consideration of subparts alone.
For example, the discrete space representation of the continuous quantities introduces a periodicity in Fourier space of period k_{0} = 2π/X_{0} (with X_{0} the grid spacing). A physical mode of wavenumber k = 2π/λ will then have, in the numerical plasma, aliases of wavenumbers k + nk_{0}, and − k + nk_{0}, with n an integer. Instabilities can arise if the physical mode couples resonantly with one of the aliases. This coupling cannot occur if k < − k + k_{0}, i.e., if λ > 2X_{0} (we note that it is NyquistShannon criterion to avoid spectral aliasing). Just as in signal processing, the strength of the aliases can be reduced by lowpass filtering the timeseries, and this is what is done by attributing a cloud shape to the superparticles. Aliases are even more reduced when the superparticle shapes have a fast decaying Fourier transform, that is, when they are smoother.
We mention in particular the following effects due to grid aliasing:

A cold beam of velocity v_{beam} becomes unstable if the Doppler shifted frequency of Langmuir oscillations is near the gridcrossing frequency k_{grid} v_{beam}. The beam is then heated. It is not the case if λ_{D}/X_{0} > 0.046.

λ_{D}/X_{0} > 1/π is needed to avoid an artificial numerical heating of a Maxwellian plasma. Otherwise, the plasma is heated up to the point where λ_{D} reaches X_{0}/π.

The rate of passage of the superparticles through the cell faces produces a highfrequency noise; the rougher the superparticle shapes, the more important is the noise.
Similarly to grid effects, a finite timestep implies that harmonics differing from a multiple of 2π/Δt are not differentiated by the algorithm, and there are time aliases as well.

This implies no other instabilities in the case of a nonmagnetizedMaxwellian plasma.

In a magnetized plasma, artificial coupling of cyclotron harmonics can lead to instabilities.
A last point is that the effects of the grid, as well as other errors, act as a random force F(t) on the superparticles. Consequently, the velocity of a superparticle undergoes a random walk, dv/dt ∝ F(t), and the kinetic energy ⟨ v^{2} ⟩ increases linearly with time. Hockney & Eastwood (1988, Sect. 9.2) shows that this is indeed the cause of plasma selfheating in superparticle simulations. This is also what we find in our thermal simulations (see Sect. 3.1.3).
B.3. Qualitative constraints on timestep and sizestep

The stepsize X_{0} of the grid (which is also roughly the superparticle size), and the timestep Δt, must be smaller than the scales of the phenomena studied. This scale can be an instability wavelength or growth rate, the cyclotron radius or pulsation, gradient scales, etc.

The same is true for the mean distance between superparticles: . This is equivalent to having a high enough number of superparticles per volume . The case λ_{relevant} = λ_{D} applies to the description of plasma behavior.

If thermal effects are important, then one should insure that the distribution function g(p) is well represented on scales where these effects are important. It requires a high enough number of superparticles per relevant volume. Birsdall & Langdon (1985, Sect. 15.19) mention that it is sufficient to have a good representation of the relevant projections of g.

The plasma should remain collisionless: the collision time should be greater than relevant timescales (instability growth rates, etc.).
B.4. Limitations for the computation of photon spectra

The highest frequency represented is2πc/X_{0} = 2πn_{x}ω_{pe}, so that high energy radiation is absent from the code and must be computed separately to extract photon spectra. This is done for example by Hededal (2005), Trier Frederiksen et al. (2010), Nishikawa et al. (2011), and Cerutti et al. (2012) from superparticle motions, with the inclusion of radiative energy losses. However, even in these cases, effects such as plasma frequency cutoff, Raizin effect, or transition radiation are not described because they are due to the backreaction of the plasma particles on the electromagnetic waves, waves that are absent from the PIC code and are only computed afterward.
Appendix C: Relativistic Harris current sheet
Harris configuration is one of the rare fully consistent solutions of the VlasovMaxwell system in a nonhomogeneous case. Its generalization to the relativistic case is not difficult, and was partly done by Hoh (1966) for a MaxwellJüttner distribution and nonrelativistic current speeds. The fully relativistic case appeared later for pair plasmas with the same temperature for both species, for example in Kirk & Skjæraasen (2003) or Pétri & Kirk (2007). Here, we propose a formulation for the relativistic case with an arbitrary temperature and mass ratio for ions (singly ionized) and electrons.
We define a frame ℛ_{0}, where the magnetic field is assumed to have the dependence given by Eq. (5), and the particles sustaining this field are assumed to have a MaxwellJüttner distribution function given in ℛ_{0} by (C.1)with s = i for ions or e for electrons, μ_{s} = 1/Θ_{s} = m_{s}c^{2}/T_{s}, p = γv/c, and K_{2} the modified Bessel function of the second kind; U_{s} is the bulk velocity of species s, and Γ_{s} the associated Lorentz factor. We note that f_{s} is normalized with respect to p to , so that is the density of species s in its comobile frame (noted with a prime), and its density in ℛ_{0}.
Inserting Eq. (C.1) into the Vlasov equation expressed in ℛ_{0}, (C.2)leads to the relation for the comobile number density (C.3)with (C.4)with sgn(q_{s}) the sign of the charge, ω_{cs} =  q_{s}  B_{0}/m_{s}, , and .
We note that the absence of electric field in Eq. (C.2) implicitly assumes that the plasma is quasi neutral in ℛ_{0}, which is true only if the overall charge density in this frame vanishes: (C.5)We now use the MaxwellAmpère equation in ℛ_{0}: . Insertion of n′(x) and B(x) leads, not surprisingly, to a pressure balance between the unmagnetized center of the sheet and the magnetically dominated outer domain: (C.6)Manipulating Eq. (C.6) and defining χ = (1 + T_{i}/T_{e})/2 we obtain (C.7)Given a temperature ratio χ, the four variables , Θ_{e}, and Γ_{e}U_{e}/c are constrained by the two Eqs. (C.4) (for s = e) and (C.7). Consequently, one needs to specify two of them. Then, the ion velocity and temperature are easily deduced with Eq. (C.4) for s = i and with χ; L and ω_{ce} expressed in units of d_{e} and ω_{pe}, that are useful for a setup in a simulation, are then deduced with Γ_{e}. A special case is when the temperatures are equal: then Γ_{i}U_{i} = −Γ_{e}U_{e}, , and .
Appendix D: A method for loading a drifting MaxwellJüttner distribution
This Appendix directly follows the notations of Sect. 4, and presents a concrete way to load a drifting MaxwellJüttner distribution in a PIC code.
Starting from the distribution in Eq. (13), we make a first change of variables (p_{x},p_{y},p_{z}) → (x,y,z) = (p_{x}/γ_{y},p_{y},p_{z}/γ_{y}), where , and we then change to cylindrical coordinates (r,θ,y) with the axis along y. We integrate along θ. Finally, a last change of variables leads to the distribution for the random variables (U,Y) ∈ [1, + ∞ [ × ]−∞, + ∞ [ : (D.1)It is then easy to obtain the marginal distribution for Y: (D.2)From this, we deduce the conditional probability distribution of U given the value y of Y: (D.3)with .
Then, for each particle, one has to generate y = p_{y} according to distribution D.2, compute a_{y}, and generate u according to distribution D.3.
For the first step, we use the method of the inversion of the cumulative distribution. This method is based on the fact that if W is a uniform random variable on [0,1] , if is the cumulative distribution of the distribution f and F^{1} its inverse, then F^{1}(W) follows the distribution f. In practice, one has to choose random numbers w_{i} in [0,1] , and the y_{i} = F^{1}(w_{i}) will be distributed according to f.
There is, however, no analytic expression for the cumulative distribution . We compute numerically J^{1}(t) on a grid of points t_{i} = i/N, i = 1...(N − 1). For each index i, we want to find y_{i} such that . We thus compute numerically the integral up to the point where it reaches t_{i}, and then attribute the value of y to y_{i}. We use the following algorithm:

1.
Choose a maximal integration step δy_{max}, and set δy = δy_{max}. Choose a tolerance tol.

2.
Start from a low enough value y_{0} such that j_{Y}(y_{0}) ≪ 1, and set y = y_{0}. Also set i = 1.

3.
Set J = 0, or if possible .

4.
Compute J = J + δy j_{Y}(y).

J is now an estimation of . If  J − t_{i}  < tol, then the desired y_{i} = J^{1}(t_{i}) is y. Set i = i + 1, and go back to step 4. If J < t_{i}, set y = y + δy, and go back to step 4. If J > t_{i}, set δy = δy/2 and go back to step 4.
We run this algorithm once for the needed μ and Γ_{0}β_{0} and store the y_{i} in a file. Then, during the particle initialization, we choose random integers i between 1 and N − 1 and set p_{y} = y_{i}.
Once y is known, we have to pick a u according to Eq. (D.3). Since a_{y} can be anything between μΓ_{0} and + ∞, we cannot generate a file before the program run, and we have to invert F on the flight. It turns out that we can integrate j_{U  y}. After some basic manipulations, we arrive at (D.4)where x = a_{y}u, l(x) = (1 + x)exp( − x), and v and w are two random numbers between 0 and 1. Inversion of the right side of Eq. (D.4) is easily done with a Newton method because of the smoothness of the function l. Starting from x = a_{y} is a good idea, and one must enforce a minimum number of iterations.
In Fig. 9, we show that the method is accurate.
Appendix E: Analytical expressions for quantities averaged over a MaxwellJüttner distribution
This Appendix details a method used to obtain the analytical expressions of Table 5 for the averaged quantities weighted by a MaxwellJüttner distribution. The notations are those in Sect. 4.
E.1. In the comoving frame
We start with averages over the distribution with zero drift velocity (i.e., the distribution in the frame comoving with the plasma). The calculations are easier when the distribution function is expressed in terms of Lorentz factors γ_{0}: (E.1)We introduce two relations^{6}: (E.3)and (DLMF 2013, Eq. (10.29.4)) (E.4)The method is then to derive Eq. (E.3) as many times as needed with respect to μ, with the help of Eq. (E.4). For example the integral of g_{0} is easily deduced from dI/dμ.
E.2. Drifting distribution
We now consider the distribution with a drift velocity . We assume that we want the average of a quantity M, (E.5)with g given by Eq. (13). We use the change of variables p_{0,x} = p_{x}, , p_{0,z} = p_{z}, which amounts to passing back into the comoving frame. Using d^{3}p/γ = d^{3}p_{0}/γ_{0}, we arrive at (E.6)where ⟨ · ⟩ _{0} means that the average is taken with g_{0}, and where γ and p are to be expressed with comoving quantities (subscript 0). With this last formula, one is left with averages in the comoving frame and can use the previous method.
All Tables
Energy conservation for simulations of a thermal plasma with no background magnetic field.
All Figures
Fig. 1 Power spectra of the total momentum of the electrons  F_{a}  ^{2} (Eq. (3)) with a = x, y, or z in the unmagnetized case. The inset is a zoom around the peak. Here n_{x} = 25, n_{t} = 500, ρ_{sp} = 16, v_{th,e} = 0.04c, T_{i} = T_{e}, box size of 25 × 25 × 25 cells, duration of 100 T_{pe}. 

In the text 
Fig. 2 Power spectra of the total momentum of the electrons  F_{a}  ^{2} (Eq. (3)) with a = x, y, or z in the magnetized case. Parameters are given in Table 1 (case labeled Fig. 2). The purple lines are the pulsations , , and . 

In the text 
Fig. 3 Wavenumber dependance of the growth rate τ^{theory}(k_{ ⊥ }), represented here as the deviation . We recall that . From Eq. (69) of Michno & Schlickeiser (2010). The squares are the resolution in k_{ ⊥ } for a box of transverse size 9, 4.5, and 3 d_{e}. 

In the text 
Fig. 4 Top: energy curves for the filamentation instability (case labeled Fig. 4 in Table 3). They are normalized by the total initial energy , which is mostly the kinetic energy of the particles. For example . The curve “energy conservation” is . After 8T_{pe}, the situation is more or less steady. Bottom: autocorrelation scale of the current amplitude. 

In the text 
Fig. 5 Top: growth of individual Fourier modes for b_{x} (for the run labeled Fig. 4 in Table 3). The modes shown are those for i = 0,1,2, and j = 0,1,...,25, and one mode every 100 modes for the remaining. We note that the graphic has been cut and that the weak modes actually fill a continuum down to an energy of 10^{7}. The sum of all 320 × 90 modes is shown in orange. We recall that mode (i,j) corresponds to (k_{z}d_{e},k_{x}d_{e}) = 2π × 20(i/640,j/180). Bottom: growth map of the Fourier modes, in units of T_{pe}. 

In the text 
Fig. 6 Difference between the growth rate measured on the total energy curve and the theoretical growth rate of the fastest modes . Except when labeled otherwise, v_{th} = 0.1c. Each symbol corresponds to a transverse box size, and each color to fixed (n_{x},n_{t}). 

In the text 
Fig. 7 Three snapshots of the filamentation instability, at times 2.5, 4, and 5.5T_{pe} (from left to right) for the simulation of Fig. 4. Colors represent the current magnitude per cell (units are number of superparticles per cell times their mean velocity), while arrows are the magnetic field. Lengths are in cell number. 

In the text 
Fig. 8 Energy curves for the tearing instability (labeled Fig. 8 in Table 4). Shown are the energy curves, for example ^{∫}dVB_{x}(t)^{2}/(2μ_{0}). They are normalized by the total initial energy (which is mostly the energy in b_{z}). The curve labeled “energy conservation” is . 

In the text 
Fig. 9 MaxwellJüttner distribution for T = mc^{2} and . Upper plot: contours are drawn from the exact expression 2πrg(r,0,p_{y}), with g from Eq. (13) and (r,0,p_{y}) defined in Appendix D. The background color map is the 2D histogram of for 10^{6} particles generated according to our method. Middle plot: normalized histogram of p_{y} for the particles (red dotted), to be compared to the exact expression in Eq. (D.2) (blue line), and for comparison (black dots) the histogram of p_{y} for particles generated the wrong way (initialization in the rest frame and boost of Γ_{0}). Bottom plot: difference between the red points and the blue curve. 

In the text 
Fig. 10 Top: electron temperatures for the hot (2) and cold (1) plasmas, from four simulations with different ρ_{sp}. The curves for ion temperatures are similar, except for an overall time dilatation by a factor ~ (m_{i}/m_{e})^{1/2} = 5. In this figure we use T_{∞}(t) = (T_{1}(t) + T_{2}(t))/2. Except for ρ_{sp} = 4 where there is significant numerical heating, T_{∞}(t) is constant in time. Bottom: halfthermalization time for ions and electrons, versus number of superparticles per cell. For ions, we have plotted t_{th}/(m_{i}/m_{e})^{1/2}. The times reported are measured as the initial slope of the temperature curves in a loglin plot, and thus correspond to t_{th}/2. We see the scaling t_{th} ∝ ρ_{sp}. 

In the text 
Fig. 11 Top: field energy levels as a function of . The colorbar is . Blue points at low f have an underresolved Debye length that could explain the mismatch with 1/f. Bottom: field energy levels versus . We clearly see the mismatch between 1/Λ_{p} and the results, even if the trend is correct. The large scatter is a hint that Λ is not a relevant parameter to describe field fluctuations. Top and bottom: each point is the result from a simulation. The field energy levels are measured as the energy in the x electric field, (with α from Eq. (A.6)), divided by the kinetic energy of the superparticles, ∑ _{sp}(γ_{sp} − 1). 

In the text 
Fig. 12 Different plasma models. Dashed arrows are transitions showing nontrivial effects. Coarsegraining and discretization are discussed in Sects. 5.1 to 5.3. 

In the text 
Fig. A.1 Grid and fields locations. 

In the text 
Fig. A.2 Simulation time duration versus number of cores, the spatial domain extending in proportion to the number of cores. The standard deviation corresponds to 4% of the mean value. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.