© 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 high-energy 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 particle-particle interactions, and the first particle-mesh 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⁶; in 2008, 10⁹; and in 2012, 10¹². 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 TRISTAN-MP (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), Photon-Plasma (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, Apar-T. 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 Apar-T. 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 Maxwell-Jü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: coarse-graining, i.e., a reduction in the number of particles by a factor reaching p ~ 10¹⁰, and the representation of the fields on a grid. We introduce the notion of coarse-graining 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 coarse-graining 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 finite-sized particles each representing up to 10¹⁰ real plasma particles, while the Vlasov-Maxwell system models a plasma macrostate described by a continuous fluid in six-dimensional phase-space. 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 Vlasov-Maxwell 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 Vlasov-Maxwell model, and a real plasma.

2. Physical model and numerical implementation

This section presents the numerical scheme used in Apar-T. Broadly speaking, Apar-T is a parallel electromagnetic relativistic 3D PIC code with a staggered grid, where the fields are integrated via Faraday and Maxwell-Ampè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 time-evolution 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 Apar-T 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₀, 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 quasi-neutral, 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 Apar-T 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 round-off 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 = μ₀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 Apar-T

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 Apar-T. 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 time-centered and time-reversible leap-frog scheme, and is second-order 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₀, 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: Maxwell-Boltzmann distribution with anisotropic temperature, boosted Maxwell-Boltzmann distribution, waterbag distribution, or Maxwell-Jüttner distribution. We have not found in the literature a method for initializing the Maxwell-Jü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¹. Files are written according to the .h5part format, and can be read by the advanced visualization and data-analysis 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 cell-averaged 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 Apar-T. 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 data-distribution model of the code, and as such is not restricted to the current parallelization model. It allows the data-IO from memory to hard drive – which is a major reason for slow-down of simulations using big data-sets – to be significantly reduced. For example, a volume-rendering 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 single-step 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 x-direction,  ∑ _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 Maxwell-Boltzmann 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₀ is strong (electron cyclotron pulsation larger than electron plasma pulsation, ω_ce ≫ ω_pe), the particle trajectories are Larmor gyrations in B₀, 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₀. 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 under-resolved 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.

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₀ and unaffected by the magnetic field. On the other hand, the electromagnetic wave with k along B₀ 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₀. 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₀. The transverse electromagnetic wave still exists with particle oscillations along B₀, unaffected by the magnetic field. Another branch appears, which is a deformation of the transverse electromagnetic wave for particle oscillations not along B₀, and starts at ω_k = 0 = ω_R with oscillations in the plane perpendicular to B₀. 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₀, and two peaks at ω = ω_L and ω_R for the component of the momentum perpendicular to B₀. 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.

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 under-resolved 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₀/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.

Simulations with under-resolved 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₀, with n an integer and X₀ the grid size) are allowed to couple when 2X₀ > λ (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 under-resolved 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 counter-streaming instability

Another standard test is to study the linear phase of the counter-streaming 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 counter-streaming electron-positron beams, with velocity and associated Lorentz factor Γ₀. There is no background magnetic field, and the particles are loaded according to a drifting Maxwell-Boltzmann 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 two-stream 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 two-fluid 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 .

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₀, 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₀,k_x,k_z). We then average the power spectrum over all the planes y = cst to obtain the power spectrum PS(t₀,k_x,k_z) =  ∑ _y0 | FT_y = y0(t₀,k_x,k_z) | ². We then repeat this procedure for several t₀. The discrete mode spectrum is sampled with (k_zd_e,k_xd_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_xd_e ≤ 5, with , which is close to the cold-fluid 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.

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.

On the other hand, the growth rates measured on the fastest modes differ from by a more systematic factor, 14% for β₀ = 0.95, 7 ± 1% for β₀ = 0.995, and 6% for β₀ = 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_xd_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 β₀ increases.

We note that e and b_z are zero in the linear two-fluid 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 counter-streaming 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 two-dimensional autocorrelation function, in the x − y plane, of the z-averaged 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²(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_ec²). 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 Maxwell-Jü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 non-trivial 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 Vlasov-Maxwell system around the drifting Maxwell-Jü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.

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 electron-positron 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 Vlasov-Maxwell 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 Maxwell-Boltzmann distributions with a mean bulk velocity U₀. 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₀ or, if U₀ is close to c, to boost every particle with the Lorentz boost corresponding to U₀. We will see, however, that this method is no longer correct when both U₀ 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².

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₀ and follows the distribution f(x,p). In the comoving frame or plasma rest frame ℛ₀, the plasma mean velocity is zero and follows the distribution f₀(x,p). We follow a group of particles. Seen from ℛ, they are in a volume d³x around x and have momentum p with a scatter d³p; seen from ℛ₀ these quantities change respectively to d³x₀, x₀, p₀, and d³p₀. The number of particles in our group is (6)so that to find the link between f and f₀ we have to find a relation between d³x₀d³p₀ and d³xd³p.

We start with the momentum. In the rest frame, our group of particles have momenta spanning a range d³p₀. Seen in the boosted frame, their momenta transform according to the Lorentz transformation, and span a new range d³p. These two volumes are thus linked by the Jacobian of the Lorentz transformation, and it can be shown that (7)where γ and γ₀ are the Lorentz factors associated to p and p₀.

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₀ = p₀/γ₀. We denote by a prime all quantities seen from this frame. In ℛ′, the particles occupy a volume d³x′. Since only one direction is contracted, and since ℛ′ moves relative to ℛ₀ with Lorentz factor , we have the relation d³x₀ = d³x′/γ₀. Similarly, ℛ′ moves relative to ℛ with Lorentz factor , and we have d³x = d³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₀(x₀,p₀), and that we boost each particle with a velocity U₀. 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 γ₀/γ does not cancel.

Consequently, boosting each particle from the rest frame leads to the expected distribution only if γ₀/γ is independent of the particle. We can write this ratio as (11)with p_0,y the y component of p₀ and , so that this is the case only if p_0,y ≪ 1 (or if the boost is non-relativistic, U₀ ≪ c). If p_0,y ≪ 1, then γ₀/γ ~ 1/Γ₀ 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, γ₀/γ is not a constant factor and Eq. (10) does not have the expected dependence on momentum p.

4.3. Loading a drifting Maxwell-Jü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 Maxwell-Jüttner distribution (Jüttner 1911; Cubero et al. 2007; Chacón-Acosta 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 Maxwell-Jüttner distribution is given by (12)with n₀ the uniform particle number density and g₀ the momentum distribution, both in the restframe, μ = mc²/T, p = γv/c, and K₂ 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β₀ (and associated Lorentz factor Γ₀), 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₀(p₀) × n₀/n. Now using n = Γ₀n₀, we arrive at (13)This distribution is normalized to unity: .

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 y-component 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² and Γ₀β₀ = 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 Maxwell-Jü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 Vlasov-Maxwell plasmas

Particle-in-cell 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 Fermi-process, 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 coarse-graining, and the discretization of the equations with the presence of a grid. Each of these steps raises questions:

• 1.

  Coarse-graining: 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 two-point 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 Vlasov-Maxwell 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 Vlasov-Maxwell 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 Vlasov-Maxwell descriptions are examined in Sect. 5.4.

5.1. The plasma parameter Λ: a coarse-graining 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 coarse-graining. Unlike fluid equations³, Eqs. (2a) − (2e) are not invariant under coarse-graining (because of the definition of the current, Eq. (2e)). The prototype of p-dependent quantities is the plasma parameter⁴, 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 coarse-graining, and since n = p × n_sp, one has that (14)with Λ = Λ_p = 1 the real plasma parameter.

In a real plasma Λ ranges from 10⁴ to 10²⁰ (for example Λ ~ 10⁶ in solar coronal loops; 10¹² in the magnetotail, magnetopause, or in typical Crab flares; 10¹⁷ 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³ to 10¹⁹. 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 coarse-graining 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 mean-free 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.

Particle-in-cell plasmas are helped by the superparticle finite sizes, which imply that the two-point 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₀, where X₀ 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 p-dependent). 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₀ < λ_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 coarse-graining 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 coarse-graining, 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¹⁰ 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.

To illustrate the dependance of the collision and fluctuation induced thermalization time, we present simulations that initially have two thermal ion-electron plasmas. The first is cold, with a temperature T_1,e(0) = T_1,i(0) = 1.6    ×    10^-3m_ec² for its electrons and ions, while the second is hot, with T_2,e(0) = T_2,i(0) = 1.8    ×    10^-2m_ec². 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)Particle-in-cell 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³ 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₁ and T₂, thermalization occurs according to (18)with n = n₁ = n₂ 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₀. Clearly, (T₁(t) + T₂(t))/2 is constant and equal to T_∞. It follows that if mass m₁ and m₂ are equal, the thermalization time is also constant and can be written (19)with Λ_∞ = n    [ϵ₀T_∞/(ne²)] ^3/2 the plasma parameter based on the temperature T_∞. It also follows that the temperatures vary exponentially as T₁ = T_∞ − 0.5 [T₂(0) − T₁(0)] exp {  − 2t/t_th } , and similarly for T₂.

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₀, but also with the hot ions on a timescale m_i/m_et₀ = 25t₀. The cold ions are heated by interactions with the hot ions on a timescale (m_i/m_e)^1/2t₀ = 5t₀, and by interactions with the hot electrons on a timescale m_i/m_et₀ = 25t₀. 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 electron-electron or ion-electron 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¹⁰ or more, while it is on the order of Λ_p ≤ 10⁴ in PIC simulations.

5.3. Field fluctuation level dependence: coarse-graining 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 coarse-graining 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₀/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₀ and having a velocity v₀, 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 finite-sized 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₀); S(ka) tends to 1 as k^-1 ≫ a.

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:  ⟨ ϵ₀E²/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₀ 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 = λ_Dk and , we arrive at (23)A primitive of the integral is u − arctan(u). We use k_max = 1/X₀ 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π² ~ 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 non-trivial 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₀ 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 Vlasov-Maxwell 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 Vlasov-Maxwell 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γ_jv_j.

At a given time t, a microstate is represented by a point in the 6N-dimensional phase-space, 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 6N-dimensional phase-space, and define the N-particle distribution function as the number density, at a given time, of these microstates in the 6N-dimensional phase-space. Given that the number of microstates in phase-space 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_mwith for sources the particles of the corresponding microstate.

When there is no background magnetic field, and when the particles are non-relativistic, 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²/(4πϵ₀r) (the Coulomb model). The Liouville equation can then be transformed to the infinite BBGKY hierarchy. The first BBGKY level involves the one-particle distribution function f₁(t,w) (with w = (x,p)) and the two-point distribution function f₂ via the correlation function g₂(t,w₁,w₂) = f₂(t,w₁,w₂) − f₁(t,w₁)f₁(t,w₂): Here C is an integral operator, vanishing with g₂. The quantity is the mean potential due to the particle distribution f₁. Since f₁(w) is the probability of finding a particle near w independently of the positions of all others, this potential does not include short-range correlations, but only long-range collective effects. It is a macroscopic quantity, just as the fields entering into the Vlasov-Maxwell system (Eq. (28)).

Approximations can then be made to truncate the BBGKY hierarchy. Evaluations of the right-hand side of Eq. (27a) can lead, depending on the hypothesis made, to Boltzmann, Landau, Lenard-Balescu, 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₂ = 0, and leads to the Vlasov-Maxwell system. Generalized to relativistic plasmas⁵ and to two species, it reads (28)\arraycolsep1.75ptwith f_1,s(t,x,p) the one-particle distribution function, for electrons (s = e) or ions (s = i), defined from f_N. It is also a smooth function, and the Vlasov-Maxwell system describes the evolution of a continuous fluid in the six-dimensional phase-space, 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 collision-induced thermalization time is infinite, and the level of electric field fluctuations ε is zero. We note that these are analytical properties. Numerical solutions of the Vlasov-Maxwell-system 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 Vlasov-Maxwell 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¹⁰ 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.

5.5. Higher-order effects of coarse-graining

We have highlighted that the coarse-graining 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 coarse-graining. Higher-order effects arise because after coarse-graining, 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 center-of-mass of the group number i. The N-particle 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 center-of-mass 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₁,...,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 coarse-grained plasma.

• The velocity dispersion of the particles within a group is ignored in the coarse-grained 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 higher-order multipole terms neglected. The same holds for the magnetic field.

6. Discussion and conclusion

6.1. Code validation

We have presented our particle-in-cell code Apar-T (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 under-resolved Debye length or under-resolved 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 Vlasov-Maxwell 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 Vlasov-Maxwell description. It also validates the new method used to load the relativistic (in both temperature and bulk velocity) Maxwell-Jü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 Apar-T 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 Vlasov-Maxwell 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¹⁰ particles, while real plasmas contain from 10⁴ to 10²⁰ particles per Debye sphere. This means that a coarse-graining step must be used, whereby of the order of p ~ 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 coarse-graining 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 coarse-graining dependent quantities can be expressed as the product of a fluid quantity (which is coarse-graining independent) and a parameter expressing a number of particles per fluid volume. Examples of such parameters include , n(c/ω_pe)³, n(c/ω_pi)³, ... 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)³ 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 point-size 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 slowing-down 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 coarse-graining 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 Vlasov-Maxwell plasmas

We have highlighted in Sect. 5.4 that a PIC algorithm simulates a plasma of finite-sized charges in their self-fields, and does not strictly solve the Vlasov-Maxwell 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 two-point 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 Balescu-Guernsey-Lenard 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 Vlasov-Maxwell 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 counter-streaming 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 counter-streaming 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 counter-streaming 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 Vlasov-Maxwell 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 non-thermal particles must include collisions to a certain degree. In this regime other models, and in particular Fokker-Planck models, have been shown to give a good description of the plasma and have provided significant results. However, Fokker-Planck models have their own drawbacks, notably that they are local and use dragging and diffusion coefficients not self-consistently derived.

These discrepancies between real, PIC, and Vlasov-Maxwell 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 Lenard-Balescu equation on the theoretical side, and P³M algorithms on the numerical side (that include short-range particle-particle interactions Hockney & Eastwood 1988). Both approaches should then be faced with results from well-controlled collisionless plasma experiments, which are presently in their infancy (see, e.g., Grosskopf et al. 2013).

---

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 ENS-Lyon, 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

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₀. Consequently, the normalized step-size is unity. Times are normalized by T₀ = X₀/c, and velocities are then naturally normalized to X₀/T₀ = c. For the fields E and B, we use e = E and b = cB. These last two quantities are normalized by e₀ = b₀ = m_ec²/(eX₀). 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 coarse-graining (). 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 leap-frog 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 b-field 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 so-called Yee lattice. It allows an easy integration of the fields, and is second-order 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).

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 leap-frog 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^tdt/2 and u ⁺  = u^t + dt/2 − qm_e/(em)e^tdt/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 self-force 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.

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 sub-domains of equal length along the z direction, and all the cells and particles of each sub-domain are assigned to a processor. To minimize communications between processors, ghost cells for the fields are added to each sub-domain. 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 over-dense 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 superparti-cle 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 non-physical 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₀ = 2π/X₀ (with X₀ the grid spacing). A physical mode of wavenumber k = 2π/λ will then have, in the numerical plasma, aliases of wavenumbers k + nk₀, and  − k + nk₀, 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₀, i.e., if λ > 2X₀ (we note that it is Nyquist-Shannon criterion to avoid spectral aliasing). Just as in signal processing, the strength of the aliases can be reduced by low-pass filtering the time-series, 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 grid-crossing frequency k_grid   v_beam. The beam is then heated. It is not the case if λ_D/X₀ > 0.046.

• λ_D/X₀ > 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₀/π.

• The rate of passage of the superparticles through the cell faces produces a high-frequency 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 non-magnetizedMaxwellian 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² ⟩  increases linearly with time. Hockney & Eastwood (1988, Sect. 9.2) shows that this is indeed the cause of plasma self-heating 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 step-size X₀ of the grid (which is also roughly the superparticle size), and the time-step Δ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₀ = 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 back-reaction 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 Vlasov-Maxwell system in a non-homogeneous case. Its generalization to the relativistic case is not difficult, and was partly done by Hoh (1966) for a Maxwell-Jüttner distribution and non-relativistic 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 ℛ₀, 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 Maxwell-Jüttner distribution function given in ℛ₀ by (C.1)with s = i for ions or e for electrons, μ_s = 1/Θ_s = m_sc²/T_s, p = γv/c, and K₂ 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 ℛ₀.

Inserting Eq. (C.1) into the Vlasov equation expressed in ℛ₀, (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₀/m_s, , and .

We note that the absence of electric field in Eq. (C.2) implicitly assumes that the plasma is quasi neutral in ℛ₀, which is true only if the overall charge density in this frame vanishes: (C.5)We now use the Maxwell-Ampère equation in ℛ₀: . 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 Γ_eU_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 Γ_iU_i =  −Γ_eU_e, , and .

Appendix D: A method for loading a drifting Maxwell-Jüttner distribution

This Appendix directly follows the notations of Sect. 4, and presents a concrete way to load a drifting Maxwell-Jü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₀ such that j_Y(y₀) ≪ 1, and set y = y₀. 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 Γ₀β₀ 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 μΓ₀ 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_yu, 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 Maxwell-Jüttner distribution

This Appendix details a method used to obtain the analytical expressions of Table 5 for the averaged quantities weighted by a Maxwell-Jü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 γ₀: (E.1)We introduce two relations⁶: (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₀ 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³p/γ = d³p₀/γ₀, we arrive at (E.6)where  ⟨ · ⟩ ₀ means that the average is taken with g₀, 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

All Figures