Features

The code is very flexible and can easily be extended to include new setups and algorithms. Some of the implemented features are presented below:

Particle Push

Particle push is the procedure that updates both the position and velocity of each particle based on the electric and magnetic fields at the location of the particle. Position and velocity are leap-frogged over each other, that is they are defined at half a timestep offset from each other. This makes the scheme second order accurate in time.

Boris Push

Usually the code uses the original Boris push in the relativistic version as defined in Boris, 1970. (The paper is annoyingly hard to find. It was titled Relativistic plasma simulation - optimization of a hybrid code and was published in the Proceedings of the Fourth Conference on Numerical Simulation Plasmas.)

Vay Modification

For simulations where particles reach more than mildly relativistic velocities the code also implements the modifications to the Boris push that were proposed by Vay, 2008.

Staggered Grid

Unlike many Eulerian codes that have all components of the fields that are stored in the grid located at the center of each grid cell, it is common in particle-in-cell codes to have different field components at different locations. This is especially useful for the commonly used electromagnetic plasma model, where the update of the field quantities in time requires the use of the curl of the dual field. The most widely used arrangement, which is also used in ACRONYM, is the so called Yee lattice, as proposed by Yee, 1966. In principle the locations of the electric and magnetic fields can be interchanged without losing the advantageous properties, however we have decided on the layout shown below, as it makes the implementation of perfect electric conductors easier.

Electric Field

The electric field components are located at grid edges, displaced to the center of the edge that is aligned with the direction of the field component. That is the x component E_x of the electric field is displaced by half a grid spacing dx along the edge e_x.
For a perfectly electrically conducting wall at x=0, this arrangement locates the y and z component exactly at the surface of the conductor, where they are trivially forced to zero.

Magnetic Field

The magnetic field components are located on the dual grid, that is on the cell faces. The x component B_x of the magnetic field is located at the center of the face defined by e_y and e_z.
Therefore the electric field components E_y and E_z are located on the cell edges around the magnetic field component B_x, which is optimal for the calculation of curl(E)_x that is needed to update B_x.

Electric Current

The components of the electric current (and mass flow, if computed for analysis purposes) are located at the mid points of cell edges.
In the electromagnetic plasma model the update of a component E_i of the electric field needs the curl curl(B)_i of the magnetic field and the corresponding component j_i of the current. Therefore it is beneficial to collocate the current density with the electric field.

Charge Density

The charge density is located in the corner of the grid cell. The center of the grid cell would seem more intuitive, as it would then represent all charges in the cell, but the location in the corner makes the formulation of the Poisson equation (for the electrostatic plasma model) and of the discrete continuity equation (when testing the electromagnetic plasma model) much easier.

Shape Functions

Particle-in-cell codes are, as the name already implies, a mix of an Eulerian gridded description of the fields and a Langrangian description of freely moving particles. Connecting the two parts requires interpolation between the staggered grid described above and the arbitrary position of a particle. Many different interpolation schemes have been devised over the years. Schemes that decompose into products of one-dimensional shape functions can easily be implemented in ACRONYM and actually a large number of them is supported. The most commonly used schemes of 1st, 2nd, and 3rd order are shown in a little bit more detail below. For additional details we refer to the paper Different Choices of the Form Factor in Particle-in-Cell Simulations by Kilian et al., 2013.

Cloud in Cell

While assigning particles simply to the nearest grid point (NGP) is possible, this trivial zeroth order interpolation introduces too much numerical noise to be useful. Assigning a finite size of one grid cell to particles strongly reduces this noise. The overlap integral between the charge cloud representing the particle and the grid cell (area weighting) can be calculated analytically beforehand and the resulting computation is equivalent to a linear interpolation. This CIC method was originally introduced by Birdsall and Fuss, 1969 and made semi-Langrangian plasma simulations feasible.

Triangular Shaped Cloud

While CIC is a big improvement over NGP, it is not perfect. A widely discussed (see e.g. chapter 5 in the standard reference Computer simulations using particles by Hockney and Eastwood, 1988) improvement is the triangular shaped cloud (TSC) that replaces the piecewise defined cloud of CIC with a continuous, piecewise linear function. This increases the effective size of the particle and the number of grid cells involved in the interpolation somewhat, but the reduced noise level and the improved isotropy usually justifies the cost.

Piecewise Quadratic Shape

The next higher order interpolation function is the piecewise quadratic shape (PQS), which is not only continuous, but also differentiable. This improves roll-off of the Fourier transform by 3dB per octave (a function that has a discontinuity in the nth derivative has a Fourier transform that rolls of proportional to 1/k^{n+1}). If the additional cost due to the growing number of grid cells that are involved in the interpolation is justified depends on the tolerable noise level and the dimensionality of the problem.

Current Deposition

Charge deposition is relatively straight forward once a shape function is selected. For an electrostatic particle-in-cell code it is sufficient to calculate the electric field and continue. An electromagnetic code on the other hand needs the electric current as a source term in the Maxwell solver. Naive deposition (interpolating q*v at the position of the particle using the shape function directly) leads to a current distribution that does not satisfy the discrete continuity equation. This in turn leads to ghost charges, violation of Gauss law and phantom forces. While it is possible to correct the computed current distribution to remove these effects (and the code does support divergence cleaning following the scheme by Marder, 1987) it is also possible to construct current deposition schemes that directly conserve continuity of electric charges.
The code supports two different such schemes (details below). Additionally it can compute and deposit quantities such as mass density, mass flow, number density, and moments of the distribution function, both as a total over all species and separated per species. These quantities are not needed for the temporal evolution, but can be valuable diagnostics when investigating the physics of a problem.

Esirkepov

The default current deposition scheme follows the method proposed by Esirekepov, 2001. It decomposes the particle motion into axis parallel motions that are weighted in such a way that the correct total current is produced and charge conservation is satisfied. This method works with all shape functions and shows good performance.

Umeda

Alternatively the code also supports the current deposition algorithm by Umeda, 2003. This algorithm decomposes the particle motion into straight segments with the constraint that the particle can only cross cell boundaries at the middle of cell faces. Only CIC and TSC shape functions are supported. It should be possible to extend this to other shape functions, but to our knowledge this is beyond the published state of the art. Main advantage of the deposition method is a different distribution of numerical noise in k-omega space that might be advantageous to some simulations.

Field Solver

The ACRONYM code supports a large number of different field solvers, that have different tradeoff between computational cost, numerical noise, and accuracy. They can be roughly divided into three classes: Electrostatic solver, where the electric field is determined by the charge density through Poissons equation; Electromagnetic solver, where both electric and magnetic fields are updated self consistently from the current density; and the radiation free plasma model, which neglects the displacement current and removes electromagnetic waves by construction. Warm collisionless plasma (in excess of 1 MK) is usually best handled by the electromagnetic solvers. However if applicable, the other solvers can speed up the computations by one or two orders of magnitude, as the restrictive CFL condition on the time step due to electromagnetic waves is removed.

Electrostatic

In the electrostatic plasma model the (longitudinal) electric field is determined by the net charge density of all particles. The Poisson equation that has to be solved is elliptic, i.e. exhibits instantaneous propagation and requires a global field solver. While surprising at first this can be thought of as the limit of long timesteps (much longer than the propagation time of electromagnetic waves through the simulation domain) or the limit of infinite speed of light. Transverse electric fields are not present in this model and the magnetic field is static. This plasma model is only applicable in the absence of large electric current and (equivalently) slow (certainly non-relativistic) particle motion. Depending on the boundary conditions (mixed boundary conditions are not supported in this plasma model) one of two different solvers is chosen.

Spectral Solver

If the domain is periodic, the Poisson equation of the electrostatic plasma model is solved using a spectral solver. The code deposits the charge density of all particles to the grid and computes the Fourier transform of it. In Fourier space the Poisson equation can easily be solved for each different k, allowing for the direct construction of of the Fourier transform of (one component of) the electric field. A subsequent inverse Fourier transform produces the electric field in real space. The global Fourier transform over all spatial dimensions is computed using a modified version of the Sandia FFT that allows to use FFTW 3 as the underlying local 1d transformation. Two features that are described in the literature but not used in all codes are the direct computation of the electric field instead of the potential, which requires two additional inverse Fourier transformations but provides slightly more accurate fields, and the use of the so-called sinc softening.

Solver using Greens Functions

If open / absorbing boundaries are desired in the electrostatic model, a different solver is used. It also uses the charge density, but then convolves that with the Greens function of free space to determine the electric field. The convolution is performed in Fourier space, where it is a simple point by point multiplication. The Fourier transformed of the charge density is determined directly from the charge density in each time step. The Fourier transformed of the Greens function is determined once, at startup by evaluation of an analytic expression for the Greens function on a double sized (to avoid wraparound in Fourier space) grid followed by a Fourier transformation. After the convolution / multiplication in Fourier space an inverse Fourier transform is necessary for each field component to determine the electric field in real space.

Electromagnetic

The most widely used plasma model in ACRONYM is the electromagnetic plasma model where both electric and magnetic fields are updated selfconsistently. In the initial setup the electric field is usually set to zero and the magnetic field is set to some (divergence-free!) profile determined by the specific setup. In each timestep the electric field is updated based on the curl of the magnetic field and the electric current, before the magnetic field is leap-frogged over it using the curl of the new electric field. (Actually the magnetic field update is split into two half-steps to have the electric and magnetic fields available at the same instant in time to update the particles.) While the details of the calculation of the current density that enters as a source term is shifted to the current deposition, the numerical approximation of the curl is retained directly in the electromagnetic field solver. Different finite differences approximations are implemented in the code:

Yee

The code of course implements the standard algorithm by Yee, 1966 that also drives the grid layout for the field quantities.
In this algorithm the curl of the electric field that is necessary to update B_x is computed from the four components of E_y and E_z that surround B_x in the staggered grid. The algorithm is straight forward and guarantees that the magnetic field does not gain a non-physcial divergence.
This algorithm is 2nd order accurate in dx, but is somewhat anisotropic and exhibits significant numerical dispersion along the coordinate axes.

4th Order

In an attempt to improve the field solver a (naive) solver that is 4th order accurate was implemented. This solver uses additional grid location of E_y and E_z to compute curl(E) that enters the update of B_x. The big advantage is that it requires very limited modification of the code as the same grid layout can be used. The major downside is that the higher order of convergence is not visible at typical cell spacings and that smaller cell spacings where the faster convergence would be beneficial are typically prohibitively expensive.

NSFD (CK, CK5)

The big improvement came with the realization that order of convergence is not as important as the error constant in front of the leading term. Here the non standard finite differences schemes proposed by Cole, 1997 and Kärkkäinen, 2006 excel. They use additional grid locations to compute the curl and optimize their relative weights to obtain numerical approximations of the curl operator that combine 2nd order accuracy with other desirable properties. A paper by Vay et al., 2011 discusses a number of potential parameters of such NSFD schemes. ACRONYM implements (in the nomenclature used there) the schemes CK and CK5, that offer maximally long time steps and maximally isotropic propagation, respectively.

M24

In another modification of the Yee scheme Hadi, 1997 has proposed to use the grid locations used by the Yee scheme, the naive 4th order solvers and the grid location required to close the big outer loop. Again different relative weights are possible. Greenwood et al., 2004 have shown that a carefully chosen set of parameters can change the phase speed of electromagnetic waves suitably to remove grid Chernenkov effects. The scheme is not perfect and can therefore not remove grid Chernenkov radiation under all conditions, but the common case of axis parallel beams profits significantly. If Chernenkov radiation is still an issue the solver can be combined with Friedman filtering.

Radiation Free (Darwin)

While the electrostatic plasma model neglects all effects O(v/c) and the electromagnetic model retains the finite speed of light, it is also possible to construct a plasma model that retains effects up to O(v^2/c^2), but ignores higher order effects. Interestingly enough there is a unique construction that does that which is equivalent to removing the displacement current and thus removing electromagnetic radiation. This results in longitudinal fields identical to the electrostatic plasma model, some transverse electric fields (but no light waves), and some corrections to the background magnetic field. Unfortunately the solver is quite complex and follows some unpublished work by Decyk. The biggest challenge is that the equation for the transverse electric field E_T needs the temporal change of the current density d j/dt as an input. Calculating this quantity from standard current deposition is too noisy. It is possible to rewrite the term, but at the cost of extra, numerically expensive, deposition steps.

Boundary Conditions

The code in general allows for several boundary conditions, that can be different for different ends of the simulation domain. (Of course a periodic boundary at the lower end of the x axis needs to be connected to a periodic boundary at the upper end of the same axis.) The main challenge when creating boundary conditions is the need to have full consistency between the boundary formulation for the electromagnetic fields and the boundary conditions for the particles. Discrepancies lead to the formation of charged layers in the best case and violent instabilities in the worst case. Periodic boundary conditions are the easiest and most widely supported boundary conditions. Reflective boundary conditions are only supported by the electromagnetic plasma model. Absorbing boundary conditions are at least supported by electrostatic and electromagnetic solvers.

Periodic Boundary Conditions

Periodic boundary conditions are straight forward and are actually handled internally in the same way as the interior inter-task boundaries, i.e. particles are sent to the new CPU and field values are copied into ghost zones. The wrap around is handled transparently by the geometry of the MPI Cartcom.

Perfect Electric Conductor

Metal walls with perfect electric conductivity are modeled by setting the tangential electric field to zero. The current densities of particles close to the metal wall are adjusted using the current of an imaginary mirror charge on the other side of the wall. Since these boundary conditions are reflective to electromagnetic waves they are complemented on the particle side with specular reflection of incident particles. Unlike a real metal wall no charged particles are absorbed or reflected.

Absorbing Boundary Conditions

The absorbing boundary conditions are implemented using a perfectly matched layer. Many different PML formulations exist. The current implementation uses an unsplit formulation. To improve absorption of low frequency waves a complex frequency shift approach is used, that requires a convolutional PML, i.e. the storage of auxiliary fields in the boundary region. On the particle side, outgoing particles are simply removed. The injection of new particles is handled by filling the ghost zone with new particles, performing one timestep inside the ghost zone and keeping the particles that enter the simulation domain.

Output

Field quantities are written to self-describing HDF5 files in Pixie format using parallel HDF5 (i.e. one HDF5 file per simulation, for all CPUs and all timesteps). Particle data can be written to HDF5 tables or to binary files using MPI-IO.