---
abstract: 'We present the Photon-Plasma code, a modern high order charge conserving particle-in-cell code for simulating relativistic plasmas. The code is using a high order implicit field solver and a novel high order charge conserving interpolation scheme for particle-to-cell interpolation and charge deposition. It includes powerful diagnostics tools with on-the-fly particle tracking, synthetic spectra integration, 2D volume slicing, and a new method to correctly account for radiative cooling in the simulations. A robust technique for imposing (time-dependent) particle and field fluxes on the boundaries is also presented. Using a hybrid OpenMP and MPI approach the code scales efficiently from 8 to more than 250.000 cores with almost linear weak scaling on a range of architectures. The code is tested with the classical benchmarks particle heating, cold beam instability, and two-stream instability. We also present particle-in-cell simulations of the Kelvin-Helmholtz instability, and new results on radiative collisionless shocks.'
author:
- 'Troels Haugb[ø]{}lle'
- Jacob Trier Frederiksen
- '[Å]{}ke Nordlund'
bibliography:
- 'references.bib'
title: 'Photon-Plasma: a modern high-order particle-in-cell code'
---

Introduction
============

Particle-In-Cell models have gained widespread use in astrophysics as a means to understand plasma dynamics, particularly in collisionless plasmas, where non-linear instabilities can play a crucial role for launching plasma waves and accelerating particles. The advent of of tera- and now peta-flop computers has made it possible to study the macroscopic properties of both relativistic and non-relativistic plasmas from first principles in large scale 2D and 3D models, and sophisticated methods, such as the extraction of synthetic spectra is bridging the gap between models and observations.

While Particle-In-Cell codes were some of the first codes to be developed for computers[@Harlow:1957; @Harlow:1964], and several classic books have been written on the subject[@birdsall:1985; @hockney:1988], modern numerical methods are in use today in the community, which did not exist then, and the temporal and spatial scales of the problems have grown enormously. Furthermore, in the context of astrophysics the modeling of relativistic plasmas has become of prime importance.

In this paper we present the relativistic particle-in-cell [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} in use at the University of Copenhagen, and the numerical and technical methods implemented in the code. The code was initially created during the ’Thinkshop’ workshop at Stockholm University in 2005, and has since then been under continuous development. It has been used on many architectures from SGI, IBM, and SUN shared memory machines to modern hybrid Linux GPU clusters. Currently our main platforms of choice are Blue-Gene and Linux clusters with 8-16 cores per node and infiniband. We have also developed a special GPU version that achieves a 20x speedup compared to a single 3GHz Nehalem core (to be presented in a forthcoming paper). The code has excellent scaling, with more than 80% efficiency on Blue-Gene/P scaling weakly from 8 to 262.144 cores, and on Blue-Gene/Q from 64 to 524.288 threads. The I/O and diagnostics is fully parallelized and on some installations we reach more than 45 GB/s read and 8 GB/s write I/O performance.

In Section II and III we introduce the underlying equations of motion, and the numerical techniques for solving the equations. We present our formulation of radiative cooling in Section IV, and in Section V the various initial and boundary conditions supported by the code. Section VI presents on-the-fly diagnostics, including the extraction of synthetic spectra. Section VII describes the binary collision modules, while Section VIII contains test problems. In Section IX we discuss aspects of parallelization and scalability, and finally in section X we finish with concluding remarks.

Equations of motion
===================

The [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} is used to find an approximate solution to the relativistic Maxwell-Vlasov system $$\label{eq:vlasov}
{\frac{\partial f^s}{\partial t}} + {\bm{u}}\cdot\frac{\partial f^s}{\partial {\bm{x}}} +
  \frac{q^s}{m^s}({\bm{E}}+ {\bm{u}}\times{\bm{B}})\cdot\frac{\partial f^s}{\partial ({\bm{u}\gamma})} = \mathcal{C}$$ $$\begin{aligned}
\label{eq:gauss}
{\nabla\cdot}{\bm{E}}& = \frac{\rho_c}{\epsilon_0} \\
\label{eq:divb} {\nabla\cdot}{\bm{B}}& = 0 \\
\label{eq:faraday} {\frac{\partial {\bm{B}}}{\partial t}} & = - {\nabla\times}{\bm{E}}\\
\label{eq:ampere} \frac{1}{c^2}{\frac{\partial {\bm{E}}}{\partial t}} & = {\nabla\times}{\bm{B}}- \mu_0{\bm{J}}\,,\end{aligned}$$ where $s$ denotes particle species in the plasma (electrons, protons, …), $\gamma = {(1-(u/c)^2)^{-1/2}}$ is the Lorentz factor, and $\mathcal{C} \equiv \partial{f^s}/\partial{t}\big|_{coll}$ denotes some collision operator.

In a completely collisionless plasma $\mathcal{C}$ is zero, but in the code we also allow for binary collisions between particles. As shown below in the tests, discretization effects in the interpolation of fields and sources between the mesh and the particles and the integration of the equations of motion lead to a non-zero, non-physical, collision term, which should be minimized, especially in the case of collisionless plasmas, but also with respect to any collisional modeling term introduced explicitly. The charge and current densities are given by taking moments of the distribution function over momentum space $$\begin{aligned}
\label{eq:rho}
\rho_c({\bm{x}}) = \int \textrm{d}{\bm{u}}\sum_s q^s f^s({\bm{x}},{\bm{u}}) \\
{\bm{J}}({\bm{x}}) = \int \textrm{d}{\bm{u}}\sum_s {\bm{u}}\, q^s f^s({\bm{x}},{\bm{u}}) \,.\end{aligned}$$ To find an approximate representation for this six-dimensional system in the particle-in-cell method so-called macro particles are introduced to sample phase space. Macro particles can either be thought of as Lagrangian markers that measure the distribution function in momentum space at certain positions, or as ensembles of particles that maintain the same momentum while moving through the volume. If the trajectory of a macro particle is a solution to the Vlasov equation, given the linearity, a set of macro particles will also be a solution to the system. Other continuum fields, which only depend on position, are sampled on a regular mesh. Macro particles are characterized by their positions ${\bm{x}}_p$ and proper velocities ${\bm{p}}_p={\bm{u}\gamma}$. They have a weight factor $w_p$, giving the number density of physical particles inside a macro particle, and a physical shape $S$. The shape is chosen to be a symmetric, positive and separable function, with a volume that is normalized to 1. For example in three dimensions it can be written $S({\bm{x}}-{\bm{x}}_p) = S_\textrm{1D}(x-x_p)S_\textrm{1D}(y-y_p)S_\textrm{1D}(z-z_p)$, and $\int S({\bm{x}}-{\bm{x}}_p) \textrm{d}{\bm{x}}= 1$. The full distribution function of a single macro particle is then $$\label{eq:pseudop}
f_p({\bm{x}},{\bm{p}}) = w_p\, \delta({\bm{p}}-{\bm{p}}_p)\,S({\bm{x}}-{\bm{x}}_p)\,.$$ Inserting the above in [Eq. \[eq:vlasov\]]{} and taking moments [@Lapenta:2006; @hockney:1988; @birdsall:1985] we find the equations of motion for a single macro particle, $$\begin{aligned}
\label{eq:pmotion}
\frac{\textrm{d}{\bm{x}}_p}{\textrm{d}t} &= {\bm{u}}_p &
\frac{\textrm{d}{\bm{u}}_p\gamma_p}{\textrm{d}t} &= \frac{q}{m}\left({\bm{E}}_p + {\bm{u}}_p \times {\bm{B}}_p \right)\,,\end{aligned}$$ where the electromagnetic fields are sampled at the macro particle position through the shape functions $$\begin{aligned}
{\bm{E}}_p &= {\bm{E}}({\bm{x}}_p) = \int \!\! d{\bm{x}}\, {\bm{E}}({\bm{x}}) \, S({\bm{x}}-{\bm{x}}_p) \\
{\bm{B}}_p &= {\bm{B}}({\bm{x}}_p) = \int \!\! d{\bm{x}}\, {\bm{B}}({\bm{x}}) \, S({\bm{x}}-{\bm{x}}_p)\,.\end{aligned}$$ The Vlasov equation ([Eq. \[eq:vlasov\]]{}) is linear, and if a single macro particles obeys [Eqs. \[eq:pmotion\]]{} a collection of macro particles, describing the plasma, will also obey it. Using that the shape functions are symmetric, and assuming that the electromagnetic fields are constant inside each cell volume we find $$\begin{aligned}
\label{eq:pfields}
{\bm{E}}_p &= \!\!\!\!\!\! \sum_{{\bm{x}}_c = \textrm{cell vertices}} \!\!\!\!\!\! {\bm{E}}({\bm{x}}_c) \, W({\bm{x}}_c-{\bm{x}}_p) \\
{\bm{B}}_p &= \!\!\!\!\!\! \sum_{{\bm{x}}_c = \textrm{cell vertices}} \!\!\!\!\!\! {\bm{B}}({\bm{x}}_c) \, W({\bm{x}}_c-{\bm{x}}_p)\,,\end{aligned}$$ where the weight function $W$ is given by integrating the shape function over the cell volume $$\label{eq:weightfunction}
W({\bm{x}}_c-{\bm{x}}_p) = \int_{{\bm{x}}_c-\frac{\Delta{\bm{x}}}{2}}^{{\bm{x}}_c+\frac{\Delta{\bm{x}}}{2}} \!\!\!\!\!\! \textrm{d}{\bm{x}}\, S({\bm{x}}- {\bm{x}}_p)\,.$$

In principle any shape function would be valid. However, in practice most PIC codes employ shape functions belonging to a family of basis functions known as $B$-splines. They have a number of useful properties[@Chaniotis2004]:

-   The particle interpolation function, a $B$-spline of order $\mathcal{O}$, is ${\mathcal{O}}-1$ times differentiable continuous, for ${\mathcal{O}}>1$

-   Particle weight functions also become $B$-splines of order ${\mathcal{O}}+1$

-   Their support is bounded; their width in number of mesh points is equal to their order + 1.

-   The sum on mesh points over the support is $\Sigma^\textrm{support}_{i} B_i^\mathcal{O} \equiv 1$, for $\mathcal{O}>0$

In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} one can select particle shape functions from one of the four lowest order $B$-splines (see also [Fig. \[fig:bsplines\]]{}) giving the weight functions:

NGP

:   ’Nearest grid point’, $$W^0(x) =
      \begin{cases}
       1 & \text{if } \left|\delta\right| \leq \frac{1}{2} \\
       0 & \text{otherwise }
      \end{cases}$$

CIC

:   ’Cloud-In-Cell’, $$W^1(x) =
      \begin{cases}
       1 - \left|\delta\right| & \text{if } \left|\delta\right| < 1 \\
       0 & \text{otherwise }
      \end{cases}$$

TSC

:   ’Triangular Shaped Cloud’, $$W^2(x) =
      \begin{cases}
       \frac{3}{4} - \delta^2 & \text{if } \left|\delta\right| < \frac{1}{2} \\
       \frac{1}{2} \left(\frac{3}{2} - \left|\delta\right|\right)^2 & \text{if } \frac{1}{2} \leq \left|\delta\right| < \frac{3}{2} \\
       0 & \text{otherwise }
      \end{cases}$$

PCS

:   ’Piecewise Cubic Spline’, $$W^3(x) =
      \begin{cases}
       \frac{1}{6}\left(4-6\delta^2+3\left|\delta\right|^3\right) & \text{if } 0 \leq \left|\delta\right| < 1  \\
       \frac{1}{6}\left(2-\left|\delta\right|\right)^3 & \text{if } 1 \leq \left|\delta\right| < 2 \\
       0 & \text{otherwise }
      \end{cases}\,,$$

where $W^i(x)$ is the one dimensional weight function, and $\delta \equiv (x_c-x_p) / \Delta x$.

![image](bsplines_2){width="0.4\linewidth"} ![image](tsc_cloud_interpretation){width="0.4\linewidth"}

Normally, there are many more particles than mesh points in a particle-in-cell model. One important benefit of introducing a higher order particle shape function is a reduction in the aliasing effects associated with the under-sampling when interpolating particle data on the mesh. When employing a higher order field solver, also, a higher order particle shape function such that the effective width of the particle shape function response and the field differencing operator response are ’not too far apart’. If the particle shape function has a very low order — say $\mathcal{O}=0$ — the strong aliasing at high frequency may be visible to a high order — say $\mathcal{O}=6$ — differential operator, which will introduce spurious contributions from the Maxwell source term(s). One should thus take care to keep the order of the particle shape function ’high’, if the field solver difference operators have order ’high’. On balance, using the highest order cubic spline interpolation with 64 mesh points in three dimensions to couple the fields and the particles has turned out to be effectively the cheapest method in most applications. The large support leads to better noise properties, and hence a lower amount of particles can be used to reach the same quality compared to using more particles and a lower order $B$-spline. The cost is also partially offset by the increased number of FLOPS per memory transfer, when compared to a lower order scheme integration cycle resource consumption.

Discretization and time integration
===================================

To solve the equations of motion for the coupled particle-field system ([Eq. \[eq:faraday\]]{}, [Eq. \[eq:ampere\]]{}, and [Eq. \[eq:pmotion\]]{}) we need to choose both a spatial and temporal discretization, and use either explicit or implicit spatial derivatives and time integration techniques. To optimally exploit the resolution, and taking into account the symmetries of the Maxwell equations, we use a Yee lattice[@yee] to stagger the fields. The charge density is located at cell centers. To comply with Gauss law ([Eq. \[eq:gauss\]]{}) the electric fields and the current density, both entering with the same spatial distribution in Ampère’s law ([Eq. \[eq:ampere\]]{}), are staggered upwards to cell faces, while magnetic fields, to be consistent with the curl operator in [Eq. \[eq:ampere\]]{}, are placed at cell edges. With this distribution (see [Fig. \[fig:yee\]]{}) the derivatives in [Eqs. \[eq:gauss\]]{}-\[eq:ampere\] are automatically calculated at the right spatial positions, and no interpolations are needed. Because of the spatial staggering the numerical derivatives commute, and the magnetic field evolution conserves divergence to round-off precision. The electric potential $\phi_E$ is located at cell centers, and the magnetic potential $\phi_B$ is at cell corners. For the time integration we stagger the proper velocity of the macro particles ${\bm{u}\gamma}_p$ and the associated current density ${\bm{J}}$ backwards time. In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} the natural order to do the different updates in a time step is (see [Fig. \[fig:time\]]{})

![Spatial staggering in the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}.[]{data-label="fig:yee"}](staggering2){width="0.8\linewidth"}

![Time staggering and integration order in the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}. ${\bm{u}\gamma}_p$ and ${\bm{J}}$ are staggered $\Delta t/2$ backwards in time, while everything else is time centered.[]{data-label="fig:time"}](timestep){width="\linewidth"}

1.  [Update the proper velocity ${\bm{u}\gamma}_p$ of the macro particles using either the Boris or Vay particle pusher]{}

2.  [Calculate the charge density on the mesh $\rho^i({\bm{x}}_p^i)$, update the macro particle positions to ${\bm{x}}_p^{i+1}$, and recalculate the charge density $\rho^{i+1}({\bm{x}}_p^{i+1})$]{}

3.  [If using a charge conserving method, find the current density from the time derivative of the charge density $\hat\partial_t \rho = [\rho^{i+1}({\bm{x}}_p^{i+1}) - \rho^i({\bm{x}}_p^i)]/\Delta t$, else calculate it directly from the particle flux ${\bm{u}}_p$ interpolated on the mesh.]{}

4.  [Given the time staggered current density, update the magnetic field using an implicit method]{}

5.  [Combine the old and new magnetic field values to make a time centered update of the electric field]{}

In between step 3 and 4 first the particle-, charge- and current boundary conditions are applied, together with an exchange of particles between MPI threads. Then the particle array, local to each thread, is sorted sequentially according to the cell position. The sorting assures optimal cache locality when computing the charge density, and step 1 to 3 can be executed in one sweep, and using small cache friendly buffer arrays for accumulating $\rho_c$ and $\hat\partial_t\rho_c$ or ${\bm{J}}$. The boundary conditions for the magnetic and electric fields are applied while solving their time evolution in step 4 and 5. The different steps are described in more detail below.

The time update shown above corresponds to a second order Leap Frog method. Compared to for example Runge-Kutta methods, or other methods where all variables are time centered, the advantage is that it requires zero extra storage for intermediate steps, and the particle update is symplectic, giving stable particle orbits. In general, symplectic integrators use a pair of conjugate variables. In this case they are the position, electromagnetic fields and charge density and the particle proper velocities and current densities. As detailed below, the leap frog update is only possible because we evaluate the ${\bm{u}}\times{\bm{B}}$ term in the Lorentz force using an implicit evaluation of the velocity. There exists higher order symplectic integrators with several substeps of the general positions and momenta composing a full time update. But only for a subset of these integrators the position is always at the mid-distance in time between the current and new momenta; a necessary condition to be able to update the particle momenta.

We have implemented a symmetric fourth order symplectic integrator[@Yoshida:1990; @Forest:1990; @Candy:1991] in the code. Due to its symmetric nature, a single, full timestep is performed by simply taking four leap frog substeps, making it an almost trivial change to implement, with updates needed only in the main driver. The method significantly increases the long term stability of the evolution for e.g. streaming plasma simulations, at the price of increasing the cost of a single time step by a factor of 4. For three dimensional simulations, where increasing the resolution by a factor of 2 costs a factor of 16 in cpu time this can be very worth while, compared to using the second order leap frog method with a $\sqrt{2}$ higher resolution.

Particle motion
---------------

To move the macro particles forward in time we have to solve [Eq. \[eq:pmotion\]]{}. Because of the time staggering it is straight forward to make an $\mathcal{O}(\Delta t^2)$ precise position update $${\bm{x}}_p^{i+1} = {\bm{x}}_p^{i} + \Delta t\, {\bm{u}}_p^{i+1}\,,$$ where the upper index $i$ indicates the iteration number. The position ${\bm{x}}_p^i$ is evaluated at $t^i$, while the proper velocity ${\bm{u}\gamma}_p^i$ is staggered backwards in time and evaluated at $t^{i-1/2}=t^{i}-\Delta t/2$.

The proper velocity update is a bit more delicate. To find a time centered Lorentz force we need a time centered value for the velocity ${\bm{u}}_p$. This gives an implicit equation for ${\bm{u}\gamma}_p^{i+1}$, and traditionally the Boris particle pusher[@boris] has been used. It is formulated so that first the time centered proper velocity is computed as the average ${\bm{u}\gamma}_p(t^i)=({\bm{u}\gamma}_p^i + {\bm{u}\gamma}_p^{i+1})/2$, and from that the normal velocity ${\bm{u}}_p(t^i)$ is derived. It is still an option in the code, but Vay has showed[@vay:2008] that the Boris pusher is not Lorentz invariant, and gives incorrect solutions in simple relativistic test cases. Instead, the Vay particle pusher calculates ${\bm{u}}_p(t^i)=({\bm{u}}_p^i + {\bm{u}}_p^{i+1})/2$ directly, as the average of the two three velocities. Even though the resulting equations for ${\bm{u}\gamma}_p^{i+1}$ involve the root of a fourth order polynomial, there is an analytic solution, and the end result is a Lorentz invariant proper velocity update. Therefore, the Vay particle pusher is now the standard option for updating the proper velocity. It is $\mathcal{O}(\Delta t^2)$ precise.

Solving Maxwell’s Equations
---------------------------

In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} Maxwell’s equations are solved by means of an implicit scheme for evolving the magnetic field, ${\bm{B}}^{i} \rightarrow {\bm{B}}^{i+1}$, followed by an explicit update of the electric field, ${\bm{E}}^{i} \rightarrow {\bm{E}}^{i+1}$. De-centering of the integrator may be employed, such that the implicit magnetic field term, ${\bm{B}}^{i+1}$, is weighted higher than the explicit term, ${\bm{B}}^i$, by setting a parameter $\alpha$, at run initialization.

With this de-centering of the implicit averaging of ${\bm{B}}$ taken at two consecutive time steps, the discretized forms of Faraday’s and Ampère’s Laws read $$\label{eq:bnplus1}
\frac{{\bf B}^{i+1}-{\bf B}^i}{\Delta
t} = - \hat\nabla^+ \times \left(  \alpha  {\bf E}^{i+1} + \beta  {\bf E}^i
\right)\,,$$ and $$\label{eq:enplus1}
\frac{{\bf E}^{i+1}-{\bf E}^i}{\Delta t} = c^2 \left( \hat\nabla^- \times \left( \alpha {\bf B}^{i+1} + \beta {\bf B}^i \right)
                                           - \mu_0 {\bf J}^{i+1} \right)\,,$$ where, $\alpha$ and $\beta=1-\alpha$ with $\alpha\geq1/2$, quantifies the forward de-centering of the implicit magnetic field term and the high frequency wave damping strength and spectral width that accompanies the choice. The hat and sign denotes that it is the discrete version of the differential operator and that it is applied downwards or upwards. $i$ indicates the iteration. The current is downstaggered in time, and the $i+1$ iteration is already calculated from the macro particle distribution, when starting to solve the Maxwell equations.

Isolating ${\bm{E}}^{i+1}$ in [Eq. \[eq:enplus1\]]{} and taking the curl, we can insert it in to [Eq. \[eq:bnplus1\]]{}. Using the vector identity $\nabla\times\nabla\times{\bm{B}}=\nabla\left({\nabla\cdot}{\bm{B}}\right)-{\nabla^2}{\bm{B}}$, which also holds for the staggered discretized operators, and using ${\nabla\cdot}{\bm{B}}=0$, produces an elliptic equation for ${\bm{B}}^{i+1}$ $$\begin{aligned}
 (1-c^2\alpha^2\Delta t^2\hat\nabla^2){\bm{B}}^{i+1} =& {\bm{B}}^i + c^2 \alpha \beta {\Delta t}^2 \hat\nabla^2 {\bm{B}}^i \nonumber \\
          &      - \Delta t \hat\nabla^+ \times {\bm{E}}^i  \nonumber \\ \label{eq:EllipticB}
          &      - c^2 \alpha \mu_0 {\Delta t}^2 \hat\nabla^+ \times {\bm{J}}^{i+1}\,.\end{aligned}$$ The right hand side of [Eq. \[eq:EllipticB\]]{} contains only known terms (at time $t_i$ for ${\bm{E}}^i$ and ${\bm{B}}^i$, and $t_{i} + \Delta t/2$ for ${\bm{J}}^{i+1}$), but the operator on the left hand side is elliptic, complicating a direct solve. We have implemented a simple iterative solver taking ${\bm{B}}^{i+1}={\bm{B}}^{i}$ as a first guess, with a solution found by successive relaxation. Elliptic equations are non-local and our solver requires repeated updates of boundaries and ghost-zones. However, the limitation to the parallel scalability is not serious, in that convergence is normally reached in 1-10 iterations for most simulation setups, with tolerances on the residual error of about $10^{-6}$. Once the relaxed solution has been found, and provided that the initial simulation setup had ${\nabla\cdot}{\bm{B}}=0$, going forward ${\nabla\cdot}{\bm{B}}=0$ is guaranteed, due the constraint (${\nabla\cdot}{\bm{B}}=0$) being implicitely built into the derivation of [Eq. \[eq:EllipticB\]]{}. Having found ${\bm{B}}^{i+1}$, the electric field is simply updated explicitly by [Eq. \[eq:enplus1\]]{}.

The field integrator is unconditionally stable for $1/2^+ < \alpha < 1^-$. However, the scheme tends to damp out high-frequency waves due to the de-centered implicit nature of the scheme, and the solver is only second order accurate for values $\alpha\approx1/2^+$. Besides providing a tunable stabilization of the field integration scheme, this parameter also determines how large a time step can be chosen, and the damping of high frequency waves. Empirically, with our highest order scheme (6$^\textrm{th}$ order fields, PCS particle assignment), for values $\alpha\geq0.525$ the scheme is numerically stable.

Forming spatial derivatives of the field quantities is done by finite differencing on a uniform (but not necessarily isotropic), mesh. A set of operators identical to those used in the [<span style="font-variant:small-caps;">Stagger</span> code]{}  – also developed and maintained in Copenhagen[@1997LNP...489..179N] – are implemented, with a choice between 2$^\textrm{nd}$, 4$^\textrm{th}$ and 6$^\textrm{th}$ order accuracy in space. Staggering of the variables on a Yee lattice [@yee] leads to highly simplified computations for the difference equations. For example, for Faraday’s Law ([Eq. \[eq:bnplus1\]]{}) the $x$-component, namely $\left[{\partial}_t {\bm{B}}\right]_x \equiv {\partial}_y E_z - {\partial}_z E_y$, along the $y-axis$ reduces to ${\partial}_t B_x = {\partial}_y E_z$. From [Fig. \[fig:yee\]]{} we see that this computation yields the desired value exactly where needed, provided that we compute the central differences at the half-staggered mesh point; a single component of the ${\nabla\times}$ operator is illustrated in [Fig. \[fig:diffoper\]]{}.

![Example of the sixth order difference operation ${\nabla\times}$ in 1D (along the y-axis). Due to the Yee mesh staggered layout of variables, the central difference is computed exactly where needed. The differential operator called `’ddyup’` (in the code) produces the derivative w.r.t. the Y-axis, up-shifted one half mesh point on the axis. For the example in the figure this produces the correct derivative of $B_y$ at the desired mesh point location, $y_{j+1/2}$. In the nomenclature adopted here, this operator is denoted $^6\hat{{\partial}}_{y}^+$. Cell centers are marked in blue, and cell edges in red. Further, compare this figure with the respective components in [Fig. \[fig:yee\]]{}.[]{data-label="fig:diffoper"}](6th_order_diff_example){width="\linewidth"}

In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}, for the 6$^\textrm{th}$ order accurate finite difference first derivative with respect to the $y$-axis the expression reads $$\begin{aligned}
 ^6\hat{{\partial}}_y^+ f_{j+1/2} =       & ~ ~ ~ ~   a \left(\frac{f_{j+1}-f_{j}}{\Delta y}\right) \nonumber \\
                                  &         + b \left(\frac{f_{j+2}-f_{j-1}}{\Delta y}\right) \nonumber \\
                                  &         + c \left(\frac{f_{j+3}-f_{j-2}}{\Delta y}\right)~,\end{aligned}$$ with coefficients $a = {25}/{21}$, $b = {-25}/{384}$ and $c = {3}/{640}$, matching those of the [<span style="font-variant:small-caps;">Stagger</span> code]{}. For our example of Faraday’s Law above, the corresponding expression to compute becomes $$\begin{aligned}
 ^6\hat{{\partial}}_y^+ E_z(y_{j+1/2}) = & ~~~~   \frac{25}{21}  \left(\frac{E_z(y_{j+1})-E_z(y_{j})  }{\Delta y}\right) \nonumber \\
                                 &      + \frac{-25}{384}\left(\frac{E_z(y_{j+2})-E_z(y_{j-1})}{\Delta y}\right) \nonumber \\
                                 &      + \frac{3}{640}  \left(\frac{E_z(y_{j+3})-E_z(y_{j-2})}{\Delta y}\right)\end{aligned}$$ for that component along the $y$-direction.

This choice of coefficients (based on a Taylor expansion) for the higher order differential operators has enabled us to consistently port simulation data from the [<span style="font-variant:small-caps;">Stagger</span> code]{} to the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}, for coupling of MHD simulations and PIC simulations.

Charge conserving current density
---------------------------------

Assuming the charge and current density on the mesh to be volume averaged, just like the fields, and inserting the phase space density of a set of macro particle ([Eq. \[eq:pseudop\]]{}) into [Eq. \[eq:rho\]]{} for the charge density, we find the charge density in a cell as $$\label{eq:rho2}
\rho_c({\bm{x}}_c) = \!\!\!\! \sum_{\textrm{particles}} \!\!\!\! q_p \, w_p \, W({\bm{x}}_p-{\bm{x}}_c)\,,$$ where $q_p$ and $w_p$ are the charge and number density of the macro particle. Another reason for choosing volume averaging, just like with the electromagnetic fields, is that the computation of the fields at the particle positions and the charge density on the mesh has to use the same interpolation technique, or particles can induce self-forces[@birdsall:1985; @hockney:1988]. In principle, one could use a similar definition for the current density $$\label{eq:jj2}
{\bm{J}}_c({\bm{x}}_c) = \!\!\!\! \sum_{\textrm{particles}} \!\!\!\! q_p \, w_p \, {\bm{u}}_p W({\bm{x}}_p-{\bm{x}}_c)\,,$$ but for the discretized version of Gauss law ([Eq. \[eq:gauss\]]{}) to hold true the correspondingly discretized charge conservation law, $$\label{eq:cc1}
\hat\partial_t \rho_c + {\hat{\nabla}\cdot}{{\bm{J}}} = 0\,,$$ has to be satisfied. In general, with the above definitions for the charge and current densities, this will not be the case. While it holds for particles moving inside a cell, particles that move through a cell boundary in a single time step can violate charge conservation. Methods have been developed[@eastwood:1991; @villasenor:1992; @esirkepov:2001; @umeda:2003] for second order field solvers that instead of using [Eq. \[eq:jj2\]]{} to compute ${\bm{J}}$ use different schemes to directly find ${\bm{J}}$ via [Eq. \[eq:cc1\]]{}. In particular Esirkepov showed[@esirkepov:2001] how to make a unique linear decomposition of the change in charge density $\hat\partial_t \rho_c$ for arbitrary shape functions into different spatial directions, and decouple [Eq. \[eq:cc1\]]{} into a set of differential equations, each involving only one component of the current density $$\label{eq:jj3}
\hat\partial_i^- J_i = \left.\hat\partial_t \rho_c\right|_i \equiv \mathcal{D}\rho_i \,.$$ The directional time derivatives of the charge density $\mathcal{D}\rho_i$ are computed on a macro particle basis, as described by Esirkepov[@esirkepov:2001], and for a second order differential operator, $$^2\hat\partial_i^- J_i = \frac{J_i^0 - J_i^-}{\Delta x_i}\,,$$ it is straight forward to compute $J_i$ from a single particle with a simple prefix sum, assuming that the contribution of a particle to the current density has compact support on the mesh and far away is zero. In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} the current density is found using this method when a second order field solver is in use. Note that *all* PIC codes based directly on the published charge conserving methods[@eastwood:1991; @villasenor:1992; @esirkepov:2001; @umeda:2003] work with second order field solvers *only*, since the order of the discretized differential operators for Gauss law and charge conservation has to be the same. The Esirkepov method has a very concise formulation and is well suited to generalize to higher order field solvers, but it is not clear how to couple the methods by Eastwood[@eastwood:1991], Villasenor and Buneman[@villasenor:1992], or Umeda[@umeda:2003] to higher order field solvers.

For higher order field solvers it is more complicated to solve [Eq. \[eq:jj3\]]{}. For example a sixth order differential operator involves sixth mesh points $$\label{eq:jj6}\footnotesize{
^6\hat\partial_i^- J_i\! = \frac{1}{\Delta x_i}\!\left[c(J_i^{++}\!\!\!\!\!- J_i^{---}) +
b(J_i^{+}\!\! - J_i^{--}) + a(J_i^0\!\! - J_i^-)\right]}$$ and the simple prefix sum that can be used at second order turns into a linear set of equations. It is also not clear what the boundary conditions, even for a single particle with compact support, should be. In the simplest case, for a single particle, the matrix will look something like $$\footnotesize{
\left( \begin{array}{cccccccc}
 a &  b     &  c & 0      &    &\ldots & 0 \\
-a &  a     &  b & c      & 0  &\ldots & 0 \\
-b & -a     &  a & b      & c  &\ldots & 0 \\
 0 & \vdots &    & \ddots &    &\vdots & 0 \\
 0 & \ldots & -b & -a     &  a &  b    & c \\
 0 & \ldots & -c & -b     & -a &  a    & b \\
 0 & \ldots &  0 & -c     & -b & -a    & a
\end{array}\right) \left(\begin{array}{c}
J_i^{5-} \\
J_i^{4-} \\
J_i^{3-} \\
\vdots \\
J_i^{3+} \\
J_i^{4+} \\
J_i^{5+}
\end{array}\right) = \left( \begin{array}{c}
 0 \\
 0 \\
\mathcal{D}\rho_i^{3-} \\
\vdots \\
\mathcal{D}\rho_i^{3+} \\
 0 \\
 0 \\
\end{array}\right)} \,,$$ where we for simplicity have absorbed $\Delta x_i$ in the definition of $a$, $b$, and $c$. There are several problems with using the above set of equations: i) it is not clear how the “no current far away” condition is implemented, and because of the few points involved the cutoff will have some impact on the result ii) the equation for a single macro particle with cubic interpolation easily involves a 10$\times$10 matrix, and iii) as formulated above, the problem is highly unstable, because the coefficients are alternating and of very different size ($a=25/21$, $b=-25/384$, $c=3/640$). Furthermore to solve the above system of equations for three directions and every single macro particle is both very costly, and will introduce noticeable round-off error on the mesh, when all contributions from all particles are summed up. In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} we have developed a fast and stable alternative, similar in structure to the one published recently[@londrillo:2010] for the Aladyn code. Instead of solving the current density for each macro particle, taking advantage of the linearity of the problem, we sum up $\mathcal D\rho_i$ directly on the mesh. Then [Eq. \[eq:jj6\]]{} can be solved on the whole domain. There is still the problem of the alternating coefficients, but that can be dealt with, by first solving for the difference $\Delta J_i \equiv J_i^0 - J_i^-$, and then using a prefix sum, just like in the second order method, to find $J_i(x_i)$. In terms of the differences the differential operator becomes $$\label{eq:jj6d}
^6\hat\partial^-_i J_i\! = \!\tilde c (\Delta J_i^{++}\!+ \Delta J_i^{--}) +
                                            \tilde b (\Delta J_i^+ + \Delta J_i^-) + \tilde a\Delta J_i^0\,,$$ where the coefficients are $\tilde a = (a+b+c)/\Delta x_i$, $\tilde b = (b+c)/\Delta x_i$, $\tilde c = c/\Delta x_i$, and the corresponding linear system, which now has to be solved on all of the mesh, is a symmetric penta-diagonal system $$\label{eq:penta}
\footnotesize{
\left( \begin{array}{cccccccc}
\tilde a & \tilde b & \tilde c & 0        &          & \ldots    & 0 \\
\tilde b & \tilde a & \tilde b & \tilde c & 0        & \ldots    & 0 \\
\tilde c & \tilde b & \tilde a & \tilde b & \tilde c & \ldots    & 0 \\
0        & \vdots   &          & \ddots   &          & \vdots    & 0 \\
0        & \ldots   & \tilde c & \tilde b & \tilde a & \tilde b  & \tilde c \\
0        & \ldots   & 0        & \tilde c & \tilde b & \tilde a  & \tilde b \\
0        & \ldots   &          & 0        & \tilde c & \tilde b  & \tilde a \\
\end{array}\right) \left(\begin{array}{c}
\Delta J_i^0 \\
\Delta J_i^1 \\
\Delta J_i^2 \\
\vdots \\
\Delta J_i^{n-2} \\
\Delta J_i^{n-1} \\
\Delta J_i^n
\end{array}\right) = \left( \begin{array}{c}
\mathcal{D}\rho_i^{0} \\
\mathcal{D}\rho_i^{1} \\
\mathcal{D}\rho_i^{2} \\
\vdots \\
\mathcal{D}\rho_i^{n-2} \\
\mathcal{D}\rho_i^{n-1} \\
\mathcal{D}\rho_i^{n}
\end{array}\right) \,.
}$$ Pentadiagonal systems have stable and efficient solvers[@penta:1996], and it costs a negligible amount of cpu time compared to the rest of the code to find the current density on the mesh, ${\bm{J}}$, from the decomposed time derivative in the charge density, $\bm{\mathcal{D}\rho}$. Furthermore, because the evaluation is done directly on the mesh, the accumulated errors in the charge conservation are smaller than if the current density is computed on a macro particle basis. To avoid having to specify boundary conditions for $\bm{\mathcal{D}\rho}$ and introduce new terms into the matrix in [Eq. \[eq:penta\]]{} we use an extra virtual cell layer. Before a timestep, by definition, the charge density will always be zero in the virtual cell layer, and the current density can then be computed correctly with a “zero-current” ansatz for the current density on the lower boundary. Afterwards the normal boundary conditions for the current density can be applied, just like when using [Eq. \[eq:jj2\]]{} for computing the current density. This method also completely decouples the solution, when running the code with multiple MPI threads. If instead of a sixth order field solver a fourth order field solver is used the system [Eq. \[eq:penta\]]{} becomes tridiagonal.

In the original version of the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} [Eq. \[eq:jj2\]]{} is used to compute the current density. The error introduced into Gauss law ([Eq. \[eq:gauss\]]{}) is mostly on the Nyquist scale, and a simple iterative Gauss-Seidel filtering technique is used to correct ${\bm{E}}$. The module also computes the error, and has been used to validate the charge conserving methods. In most applications, the relative error with the old method can be kept on a $10^{-4}$-$10^{-5}$ level by running the filter 5 to 10 times per iteration, but there is no unique way to correct the electric field, and the divergence cleaning introduces tiny electric fluctuations, which couple back to the particles through the Lorentz force. Apart from the higher cost of the elliptic filter when running on many cores, the method is worse at conserving energy and has a larger numerical heating rate for cold plasma beams than a charge conserving method, and without a current filter it can be unstable. When running the code with charge conservation and in single precision at high resolution (e.g. $\sim\!\!1000^3$ cells and 10-100 billion particles), over time the numerical round-off noise will eventually build up errors both in Gauss’ law and in the solenoidal nature of the magnetic field ([Eq. \[eq:divb\]]{}). A single filtering step on the electric and magnetic fields every $\sim$5$^\textrm{th}$ iteration is enough to keep the relative errors at the $10^{-5}$ level. When running the code in double precision we have not seen any need to use divergence cleaning, and if the initial and boundary conditions obey [Eqs. \[eq:gauss\]]{} and \[eq:divb\], then the relative error typically stays below $\sim\!\!10^{-7}$.

Radiative cooling
=================

Very high energy electrons loose momentum by emitting radiation. The emission in itself is a very valuable diagnostic and the extraction of the spectrum is treated below, but the energy loss for the high energy particles is normally not computed in particle-in-cell codes. Building on the work on radiative losses done by Hededal[@Hededal:2005b] in our old PIC code, we have developed a numerically stable method for correctly calculating radiative losses in the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}. The radiated power from a single particle $P_\textrm{rad}$ can be written as[@Hededal:2005b] $$\label{eq:pradv}
P_\textrm{rad} = \frac{\mu_0 q^2}{6\pi c} \left[ \gamma^4 \dot{\bm{u}}^2 +
\frac{\gamma^6}{c^2} (\dot{\bm{u}}\cdot{\bm{u}})^2\right]\,,$$ where dot denotes the time derivative.

Denoting the proper velocity ${\bm{p}}={\bm{u}\gamma}$, and using the identities ${\bm{p}}\cdot\dot{\bm{p}}= \gamma^4 {\bm{u}}\cdot\dot{\bm{u}}$, and $\dot{\bm{p}}= \gamma \dot{\bm{u}}+ c^{-2}\gamma^3 {\bm{u}}({\bm{u}}\cdot\dot{\bm{u}})$, we can rewrite [Eq. \[eq:pradv\]]{} to $$\begin{aligned}
\nonumber
P_\textrm{rad} & = \frac{\mu_0 q^2}{6\pi c} \left[ \gamma^2 \dot{\bm{p}}^2 -
\frac{1}{c^2} (\dot{\bm{p}}\cdot{\bm{p}})^2\right] \\ \label{eq:prad}
& = \frac{\mu_0 q^2}{6\pi c} \left[ \dot{\bm{p}}^2 + \frac{1}{c^2}
                                   (\dot{\bm{p}}\times{\bm{p}})^2\right] \,.\end{aligned}$$ Notice how [Eq. \[eq:prad\]]{} only involves proper velocities, and is numerically stable at both high and low Lorentz factors, while the first version relies on a cancellation between $\gamma^2$ and ${\bm{p}}^2$, and [Eq. \[eq:pradv\]]{} uses three velocities, prone to numerical errors.

The radiative cooling always acts in the opposite direction to the proper velocity vector, and therefore the change in the length of the proper velocity vector is directly related to the change in the kinetic energy and $P_{rad}$: $$\begin{aligned}
\nonumber
\dot{\bm{p}}_\textrm{rad} &= \frac{{\bm{p}}\cdot\dot{\bm{p}}_\textrm{rad}}{{\bm{p}}^2} {\bm{p}}\\ \label{eq:udot}
  &= - \frac{\mu_0 q^2}{6 \pi m c} \frac{\gamma}{{\bm{p}}^2}
   \left[ \dot{\bm{p}}^2 + \frac{1}{c^2} (\dot{\bm{p}}\times{\bm{p}})^2\right] {\bm{p}}\,,\end{aligned}$$ where we have used that ${\bm{p}}\cdot\dot{\bm{p}}= c^2 \gamma\,\dot\gamma$, and $m c^2 (\gamma - 1\dot{)\,} = -P_\textrm{rad}$.

In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} the Boris or the Vay pusher is used to advance the particles. To integrate the effect of radiative cooling, for simplicity and to keep the scheme explicit, we assume that the cooling in a single timestep only changes the energy with a minor amount. The particle pushers advances the four velocities from time step $t - \Delta t/2$ to $t + \Delta t/2$, while accelerations are computed time centered at $t$. Below we denote the three times with $-$, $+$, and $0$ respectively.

To get a time centered cooling rate ${\bm{p}}_0$, $\dot{\bm{p}}_0$, and $\gamma_0$ are needed. The proper acceleration $\dot{\bm{p}}_0$ is already naturally time centered, but with the Vay Pusher, for consistency one should use the time centered three velocity to compute $ww_0$ and $\gamma_0$, which however is numerically imprecise when using single precision. To get a more numerically precise, albeit ever so slightly inconsistent, measure for $\dot{\bm{p}}_0$ and $\gamma_0$ we use the time averaged proper velocities $$\begin{aligned}
{\bm{p}}_0 & = \frac{{\bm{p}}_+ + {\bm{p}}_-}{2} &
\gamma_0 & = \sqrt{1 + {\bm{p}}_0^2 c^{-2}} \,.\end{aligned}$$ These values can then be plugged into [Eqs. \[eq:prad\]]{} and \[eq:udot\] to find the change in momentum due to radiative cooling as $$\begin{aligned}
\nonumber
{\bm{p}}_+ & = {\bm{p}}_- +  \dot{\bm{p}}_0 \Delta t \\ \nonumber
& = {\bm{p}}_- + (\dot{\bm{p}}_\textrm{EM} + \dot{\bm{p}}_\textrm{rad})\Delta t \\
       & = {\bm{p}}_+^\textrm{EM} + \dot{\bm{p}}_\textrm{rad} \Delta t\,.\end{aligned}$$ In principle one could find the converged solution to the above non-linear equation system, under the assumption that cooling is a small correction, by iterating the radiative cooling computation a couple of times, while updating ${\bm{p}}_+$ and $\dot{\bm{p}}_0$. But in practice we have found that the direct explicit calculation of the cooling rate done by making no iterations is acceptable.

Initial and boundary conditions
===============================

Setting up consistent boundary and initial conditions for a particle-in-cell code can be non-trivial, due to the mixed particle and cell nature and the staggering in space and time. In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} the initial and boundary conditions are supported through the loading of different boundary and initial condition modules, and over time a number of modules have been developed.

Particle Injection
------------------

Macro particles in a particle-in-cell code sample phase space and are by construction placed using a random generator. Any call to a random generator is related to the position in the mesh, where the pseudo-random number is needed. To make the random number generation scalable, but independent of the parallelization technique, we have implemented two types of random generators: i) We have developed a multi stream variant of the Mersenne Twister random generator[@mersenne] to generate high quality random numbers, and setup one stream per $xy$-slice of cells. This generator is used in e.g. shock simulations where the particle average density and velocity is a function of a single coordinate. ii) In the case of more complicated setups, e.g. when using snapshots from MHD simulations as described below, a simpler random generator is used, where the state is contained in a single 32-bit integer, but each mesh point has a its own state. The Mersenne Twister generator is used to generate the initial seed in each cell for the simple random generator. To sample the velocity phase-space we have implemented cumulative 2D and 3D relativistic Maxwell distributions[@Dunkel:2009] using an inverted lookup table[@birdsall:1985]. It is important to use the correct dimensionality when initializing the velocities, or the corresponding 2D or 3D temperature will be incorrect. The particle positions are correspondingly injected either uniformly in an $xy$-slice or in a single cell. Notice that using an injection method in a full $xy$-slice allow for larger density fluctuations inside the cells, while if using a cell-by-cell injection method there will only be fluctuations on the sub-cell level. Depending on the physical model on hand one or the other method may be more desirable. We typically use a standard technique to inject particles in pairs of different charge type to avoid having free charges initially. But the code is flexible, and has a built in module to correct Gauss law. Therefore it is also possible to inject particles completely at random and correct the electric field to include the non trivial effect of electrostatic fluctuations from the initial non-neutral charge distribution.

Boundary conditions
-------------------

We handle periodic boundaries trivially by padding the domain with ghostzones, a technique also used at the edge of MPI domains, and copying fields from the top of the domain to the bottom and vice versa while updating boundaries. The particles, on the other hand, are simply allowed to stream freely and the position is in all but a few cases calculated modulo the domain size. To maintain uniform numerical precision even for very large domains the position is decomposed internally as an integer cell number, and a floating point number giving the fractional position inside the cell.

Reflective boundaries are implemented using “virtual particles” on the other side of the reflecting boundary, taking in to account the symmetries of the Maxwell equations and the staggering of the mesh. By convention the boundary is placed at the center of the cell (i.e. the charge density is *on* the boundary), below for simplicity taken to be an upper boundary. When the charge and current densities on the mesh are calculated for each particle inside the boundary a corresponding virtual “ghost particle” outside the boundary has to be accounted for. Particles close to the reflecting boundary will contribute to the charge and current densities outside the boundary, while the virtual particles will contribute a corresponding charge and current density inside the boundary. This can most efficiently be calculated by disregarding the boundary at first after a particle update, and calculate the charge and current density as if the particles were streaming freely. The contribution to the charge and current densities inside the boundary from the virtual particles can then be calculated by taking into account the symmetry. It is easily seen that this corresponds, up to a sign, to the contribution of the normal particles outside the boundary $$\begin{aligned}
   \rho_c(x_\textrm{b} - \delta x) &=    \rho_c(x_\textrm{b} - \delta x) +    \rho_c(x_\textrm{b} + \delta x) \\
{\bm{J}}_\parallel(x_\textrm{b} - \delta x) &= {\bm{J}}_\parallel(x_\textrm{b} - \delta x) + {\bm{J}}_\parallel(x_\textrm{b} + \delta x) \\
 J_\perp(x_\textrm{b} - \delta x) &=  J_\perp(x_\textrm{b} - \delta x) -  J_\perp(x_\textrm{b} + \delta x)\,,\end{aligned}$$ where $x_\textrm{b}$ is the position of the boundary, $\delta x$ is the distance; i.e. integer $\Delta x$ for centered quantities, and half integer for staggered, $\parallel$ indicates components parallel and $\perp$ perpendicular to the boundary. The perpendicular component of the current density is staggered and anti-symmetric across the boundary. Only after calculating the current density on the mesh are all particles outside the boundary reflected according to $$\begin{aligned}
x &\to 2\, x_\textrm{b} - x & v\gamma_\parallel &\to v\gamma_\parallel & v\gamma_\perp &\to - v\gamma_\perp
\,.\end{aligned}$$ When the charge and current densities are correctly calculated inside and on the boundary we can start considering the values outside. The staggered fields, i.e. the perpendicular current density and electric field, the parallel magnetic field, and the magnetic potential, may be shown to be antisymmetric across the boundary $$\begin{aligned}
   J_\perp(x_\textrm{b} + \delta x) &=    - J_\perp(x_\textrm{b} - \delta x) \\
   E_\perp(x_\textrm{b} + \delta x) &=    - E_\perp(x_\textrm{b} - \delta x) \\
{\bm{B}}_\parallel(x_\textrm{b} + \delta x) &= - {\bm{B}}_\parallel(x_\textrm{b} - \delta x) \\
   \phi_B(x_\textrm{b} + \delta x) &=    - \phi_B(x_\textrm{b} - \delta x) \,.\end{aligned}$$ Conversely, the centered fields (in the direction of the boundary), charge density, parallel current densities, electric fields, the electric potential, and the perpendicular component of the magnetic field are symmetric across the boundary $$\begin{aligned}
   \rho_c(x_\textrm{b} + \delta x) &=    \rho_c(x_\textrm{b} - \delta x) \\
{\bm{J}}_\parallel(x_\textrm{b} + \delta x) &= {\bm{J}}_\parallel(x_\textrm{b} - \delta x) \\
   \phi_E(x_\textrm{b} + \delta x) &=    \phi_E(x_\textrm{b} - \delta x) \\
{\bm{E}}_\parallel(x_\textrm{b} + \delta x) &= {\bm{E}}_\parallel(x_\textrm{b} - \delta x) \\
   B_\perp(x_\textrm{b} + \delta x) &=    B_\perp(x_\textrm{b} - \delta x) \,,\end{aligned}$$ and we only need to determine the values of the centered fields at the boundary. The charge and current densities are derived from the particle distribution and are therefore already given at all points, including the boundary. $B_\perp(x_\textrm{b})$ is the only unknown component of the magnetic field, and it can be computed from the solenoidal constraint $\hat\nabla^+\cdot {\bm{B}}=0$ calculated at the boundary. Given that the magnetic field is the first to be evolved forward in time it is then possible to self consistently calculate the parallel electric field on the boundary ${\bm{E}}_\parallel(x_\textrm{b})$, directly from the evolution equation.

Outflow boundaries are less constrained than reflecting boundaries, given the extrapolating nature, and various types of damping layers and extrapolations have been considered[@PML2D; @PML3D; @umeda:2001]. We use a damping layer of mesh points — typically 20 — to damp all perpendicular components of the electromagnetic fields, effectively absorbing reflected waves, and combine it with an extrapolating boundary condition (i.e. symmetric first derivative). We allow particles to still generate current and charge densities on the mesh inside the box until they are well outside the boundary. As long as the disturbances in the outgoing flow are small this works well and is stable for extended run times.

Sliding window and injection of particles
-----------------------------------------

In highly relativistic flows it can be advantageous to simulate the plasma in a frame where the region of interest is relativistic. This constrains the time such a region can be followed, because the computational domain has to be continuously expanding or has to have an enormous aspect ratio in the flow direction. This is the case for example for relativistic collisionless shocks. An alternative to this is to use a sliding window as the computational domain centered on the region of interest and moving with the same velocity maintaining it at the center of the box[@movingframe]. We have implemented this technique in the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}, together with a moveable open boundary and particle injector. If the window is moving with a velocity $v$ relative to the lab frame then in each iteration we check if $v (t_\textrm{i} - t_{\textrm{old}}) > \Delta x$ and in that case we move the box a full mesh point, if necessary. I.e. if $v = 0.5 c$ and $\Delta t = 0.5 \Delta x / c$ the code will roughly move the window one point in every 4 iterations. The move is implemented by removing one cell at one end of the box, translating everything one point, and injecting an extra layer of inflow particles at the other end. This technique was used to successfully capture the long term evolution of an 3D ion-electron collisionless shock[@haugboelle:2011].

Embedded particle-in-cell models {#sec:embed}
--------------------------------

MHD models have successfully been applied for many years to study the large scale structure of plasmas from laboratory length scales to the largest scales in the universe. MHD can reach over such enormous scales, because a statistical description of the plasma is employed, where the microscopic state is captured using statistical quantities like the temperature, viscosity and resistivity. On the other hand, the kinetic description of a plasma used in a particle-in-cell code is ideally suited to investigate non-thermal processes, such as particle acceleration and particle-wave interactions, and can be used to model and understand a much broader range of plasma instabilities than MHD. The drawback is that exactly because the plasma is described in kinetic terms, explicit particle-in-cell codes, like the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}, have to respect microscopic constraints and resolve the Debye length, the plasma skin depth, and the light crossing time of a single cell. Recently, we have developed a technique to couple the two approaches using the results from MHD simulations to supply initial and boundary conditions for our PIC code. This has enabled us to for the first time make a realistic particle-in-cell description of active coronal regions[@baumann:2012a; @baumann:2012b], and investigate the mechanism that accelerates particles in the solar corona. The MHD code is used to evolve the global plasma over several solar hours, and a snapshot just before e.g. a major reconnection event is used to study a small time sequence, of the order of tens of solar seconds, using the PIC code.

Given an MHD snapshot, typically a smaller cutout from a larger simulation of the region of interest, we first interpolate to the resolution that is to be used in the PIC code. The interpolated magnetic fields are corrected with a divergence cleaner that uses the same numerical derivative operator that are used in the PIC code, to assure that the initial magnetic fields are solenoidal to roundoff precision. In a given cell the density in the MHD snapshot is used to set the weight of individual particles. The particles are placed randomly inside the cell, but in pairs, so that initially there are no free charges. The velocity of the particles have three contributions, $$v\gamma = v\gamma_{\textrm{bulk}} + v\gamma_{\textrm{thermal}} + v\gamma_{\textrm{current}}\,.$$ The bulk momentum is taken from the MHD snapshot. The thermal velocity is sampled from a Maxwell distribution using the MHD temperature, and finally the current speed is found from the ideal MHD current $\mu_0{\bm{J}}= \mu_0 \sum_i q^i n^i v^i_\textrm{current} = \nabla \times {\bm{B}}$. The average momentum has to correspond to the bulk momentum in the MHD snapshot. Taking into account the mass ratio, then for example in a neutral two component proton-electron plasma, with $n = \rho_{\textrm{MHD}} / (m_e + m_p)$, the weighting is $$\begin{aligned}
v^e_\textrm{current} &= - \frac{m_p}{\mu_0 |q| \rho } \nabla \times {\bm{B}}&
v^p_\textrm{current} &=   \frac{m_e}{\mu_0 |q| \rho } \nabla \times {\bm{B}}\,,\end{aligned}$$ and in general the correct way is to use harmonic weighting. Finally, an initial condition has to be specified for the electric field ${\bm{E}}$. One possibility is ${\bm{E}}=0$. It satisfies Gauss law – the plasma is neutral initially – but is inconsistent with the EMF from the MHD equations. If using this initial condition, in the beginning of the run a powerful small scale electromagnetic wave is launched throughout the box, when the $\partial_t {\bm{E}}$ term in Ampère’s law adjusts the electric field on a plasma oscillation time scale. Another possibility is to set ${\bm{E}}= - {\bm{u}}\times {\bm{B}}$, in accordance with the ideal MHD equations. Then there is no guarantee that Gauss law is satisfied. In the code the second choice is used, but small scale features in the electric field are corrected by running the build-in Gauss law divergence cleaner for a few iterations. The remaining difference is adjusted by changing slightly the ion- and electron-density, making the plasma charged. Typically this only leads to small scale changes. The resulting initial electric field then both satisfies Gauss law, and is almost in accordance with the MHD EMF. Apart from using the MHD EMF, we have also options to add the Hall and Battery effect terms.

To evolve the model, boundary conditions on all six boundaries are needed. For the plasma they are constructed exactly like the initial conditions. They can easily be made time dependent, by loading several MHD snapshots and interpolating in time. In every time step, when applying the boundary conditions, first all particles in the boundary zones are removed – also particles that have crossed the boundary from the interior of the box – and are then replaced with a fresh plasma, according to the MHD snapshot. The boundary is typically 3 zones broad; enough to allow for the sixth order differential operators on the interior of the mesh, and enough to make a well defined charge and current density with the cubic spline interpolation. This plasma is retained when evolving, and particles from the boundary zones are allowed to cross into the interior of the computational domain. By maintaining a correct thermal distribution in the boundary zones, and simply letting the dynamics decide which particles stream into the box, the inflow maintains a perfect Maxwell distribution, and the resulting plasma is practically identical to what would have been obtained with the open boundary method of Birdsall et al[@birdsall:1985], but is much simpler to implement, and correctly accounts for bulk velocities and currents in and out of the box. If there is a differential between the charge or electric current in the boundary. Inside the domain the $\partial_t {\bm{E}}$ term adjusts the plasma almost instantaneously, and the balance is maintained. This boundary condition is very similar to a perfect thermal bath, but with in- and out-going bulk velocities and currents.

We do not keep the fields fixed at the boundary condition, but instead let them evolve freely, only subject to reasonable symmetry conditions at the boundary, which keep them consistent with the Maxwell equations. To respect the symmetry of the equations we let

-   [$E_\perp$ symmetric, ${\bm{E}}_\parallel$ antisymmetric]{}

-   [$\partial_{\perp} B_\perp$ symmetric, ${\bm{B}}_\parallel$ symmetric]{}

-   [$\rho_c$ and ${\bm{J}}$ specified according to MHD snapshot]{}

-   [$\phi_E$ symmetric, $\phi_B$ antisymmetric]{},

where $\parallel$ are the components parallel with the boundary and $\perp$ the component perpendicular to the boundary. For the relatively short times that we have evolved imbedded PIC simulations[@baumann:2012a; @baumann:2012b] these boundaries are stable.

A severe limitation for coupling PIC and MHD codes is that in many situations the Debye length, plasma frequency and other microphysical length and time scales are many orders of magnitude smaller than the scales of interest. For example, in the solar corona the Debye length is measured in millimeters, while interesting macroscopic scales are measured in megameters. If we were to simulate the true system using an explicit PIC code we would need roughly $(10^8)^3$ cells, which is computationally unfeasible in the foreseeable future. To circumvent this problem we have developed a novel method, in which we rescale the physical units while maintaining the hierarchy of time, length and velocity scales. The [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} is flexible and can be employed with a range of different unit systems. Furthermore all natural constants are maintained in the code. The rescaling technique is discussed in detail in Baumann et al[@baumann:2012a].

Diagnostics
===========

Particle tracking and field slicing
-----------------------------------

Particle-in-Cell simulations are routinely run with billions of particles, with each particle taking up $\sim$50 bytes, and for $10^3$-$10^6$ iterations. To store the full data set from every single iteration would take up petabytes of storage, and is impracticable. Instead, a standard practice when running PIC simulations is to only store every n$^\textrm{th}$ particle in a snapshot, and dump snapshots with a reduced frequency, decreasing the data volume dramatically. But to understand the underlying physics of for example particle acceleration in detail a more fine-grained approach is warranted. To that end we have implemented dumping of field slices and tracking of individual particles. Any particle in the code can be tagged for particle tracking, according to a number of criteria. For example based on its energy, at random, or according to the specific ID of the particle, which is reproducible between runs. The tagged particles are harvested by each MPI thread individually, and together with the position and momentum the local values of the current, density, electric and magnetic field are recorded, by interpolation from the mesh to the particle position. Everything is arranged in a single array that is sent to the master thread. The master thread then dumps the particle records to a single file, appended to in each iteration. On x86 clusters we can sustain tracing a million particles without significant performance degradation, while on clusters with weaker CPUs, such as BlueGene/P, we are limited to $\sim10^5$ particles. To put the particle tracks into context, data on of the field evolution is also needed. To save time-resolved field data we have made a field slicing module, where a large selection of fields (e.g. the electric field ${\bm{E}}$, ${\bm{E}}\cdot{\bm{B}}$, ${\bm{J}}\times{\bm{B}}$ etc) may be stored as 2D slices. The extent of a slice, and the number of field layers in the perpendicular direction to the slice, used for averaging, is user selectable. These two techniques have been used in concert, to understand the mechanism behind particle acceleration in reconnection events in the solar corona [@baumann:2012a; @baumann:2012b], and in collisionless shocks [@haugboelle:2011]. Both diagnostics a interactively steerable: parameters can be changed, and diagnostics can be turned on and off while a simulation is running.

Synthetic Spectra
-----------------

The radiation signature from an large number of of accelerated charges is not easily computed analytically from first principle, for plasmas with complex fields topologies, rich phase space structure, and temporal evolution. Application examples include relativistic outflows and collisionless shocks; more specifically, for example, gamma-ray bursts, where magneto-bremsstrahlung is very likely to constitute a major part of the observational signal.

However, since particle-in-cell codes automatically provide all variables needed for producing a radiation spectrum, namely ${\bm{r}}$ (position), ${\bm{\beta}}$ (${\bm{v}}/c$, velocity), and $\dot{{\bm{\beta}}}$ ($\dot{{\bm{v}}}/c$, acceleration), a radiation spectral synthesizer has been integrated into the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}. We need only designate observer position(s) and match the frequency range to the plasma conditions to complete the setup for the computing the radiation integral ([Eq. \[eq:radpower\]]{} below). During run-time the synthesizer computes the radiation signature for an ensemble of charged particles (in most cases electrons) in the simulation volume,; the formula for the spectrum is given by $$\begin{aligned}
\nonumber
    \frac{d^2W}{d\Omega d\omega} & = \frac{\mu_0cq_e^2}{16\pi^3} \times \\ \label{eq:radpower}
    & \!\!\!\! \left| \int^{\infty}_{-\infty} \frac{\textbf{n} \times \big(
    (\textbf{n}-{\bm{\beta}}) \times
    \dot{{\bm{\beta}}}\big)}{(1-\textbf{n}\cdot{\bm{\beta}})^2}e^{i\omega(t'-\textbf{n}\cdot{\bm{r}}_0(t')/c)}dt'\right|^2\end{aligned}$$ with $t'$ the retarded time and $\textbf{n}$ the direction of the observer. A first comprehensive and thorough study of the spectral synthesis method is given by Hededal[@Hededal:2005b], which also covers a range of test examples. While in that study, the spectral synthesis was done as post-production, in the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} all parts of the integration are done at run-time, with very little overhead, even for large numbers of particle traces[@Trier:2010; @Medvedev:2011].

The discretization of [Eq. \[eq:radpower\]]{} is done in four parts:

1.  *Frequency range* is specified as an interval and is discretized into $N_{\omega}$ bins, typically of order $10^3$, either with linear or (more often in practice) with logarithmic binning.

2.  *Observer positions*, often more than one ($N_{obs}>1$), are specified at run initiation time (input), typically with directionality perpendicular to a sphere centered on the simulation volume, or any important direction.

3.  *Time subsampling* may be chosen; this partitions the integration for every simulation time step into a number of subcycled integration intervals: Subcycling is employed on particles selected for synthetic tracing since, for highly relativistic situations, the retarded electric field can be extremely compressed in spikes (for example in the case of synchrotron motion with $\gamma(v)\ll1$). In such situations the subcycling provides a much cheaper alternative than to restrict the Courant condition for the entire simulation.

4.  *Radiative regions* are defined in a uniform meshing (independently of MPI and simulation mesh geometries), which gives the advantage of offering the possibility to sample very local volumes of the plasma. This may be of interest in simulations with — at the same time — subvolumes of very high and low anisotropy, such as is the case in fully resolved shock simulations, and relativistic streams.

A seamless integration into the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} has made the synthesis module computationally efficient, and due to the embarrassingly parallel nature of the spectra collection procedure plasma simulations have been run with millions of particles used for sampling the synthetic spectra. Particles are chosen for spectral integration before or during a run by tagging for synthetic sampling. The detailed sampling is important in highly relativistic cases with bulk flows – for example when investigating the radiation signature from relativistic collisionless shocks, and streaming instabilities, more generally [@Hededal:2005b; @Trier:2010; @Medvedev:2011].

Binary collision operators
==========================

The classic particle-in-cell framework does not take into account physical collision processes, and all particle interactions are mediated through the electromagnetic fields on the mesh. Low energy electromagnetic waves, and large impact parameter electrostatic scatterings between charged particles can be resolved directly on the grid through particle-wave interaction, but the photon energy is limited by the grid resolution and binary electro static scattering is not correctly represented. To allow for binary interactions, high energy photons, and in general interactions of neutral and charged particles and gas drag forces, we have to model them explicitly in cases where they are of importance, such as in high energy density plasmas, and partially ionized mediums. The [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} supports the inclusion of particle-particle interactions, decay of particles, and allow for neutral particles, in particular photons, in the model. The first implementation of binary interactions—in particular Compton scattering—together with tests of the method was given by Haugb[ø]{}lle[@Haugboelle:2005] and Hededal[@Hededal:2005b]. This first implementation motivated the name for the code, the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}. The Compton scattering module was used to model the interaction of a gamma-ray burst with a circumstellar medium[@Trier:2008.2; @Trier:2008.3; @Trier:2008.1]. Coulomb collisions have later been incorporated into the framework, to study particle acceleration in solar active regions [@Baumann:2012c].

Compton Scattering and splitting of particles
---------------------------------------------

The classic Monte Carlo approach to scattering is based on a cut-off probability: first a probability for the process is computed and then it is compared with a random number. If the random number is lower than the threshold the scattering for the full macro particle pair is carried through, and otherwise nothing happens. This probabilistic approach is straight forward both numerically and conceptually, but it can be noisy, in particular when interaction effects are strong but have low probability.

In the code the natural domain to consider is a single cell, partly because that is by definition the volume of a single macro particle, partly because some interactions (e.g. electro static interactions) are mediated by the grid at larger scales. In a PIC simulation typical numbers are $10-10^3$ particles per species per cell, and a probabilistic approach would result in an unacceptable level of noise. Consider a beam incident on a thermal population: The first generation of scattered particles may be computed relatively precise, but the spectra of later generations will require an excessive amount of macro particles, if they all represent an equal amount of physical particles, given the exponentially lower number density of later generations. Another well known consequence is that the precision scales inversely proportional to *the square root* of the number of particles. This is a problematic limitation, when the higher order generations are important ingredients of the physics, and the scattering process is not just a means to thermalizing or equilibrating the phase space distribution.

![Sketch of the scattering mechanism in the code of an incident macro particle (shown as blue/dark gray) on a target macro particle (shown as red/light gray) resulting in the creation of a scattered macro particle pair.[]{data-label="fig:detailedcollisions"}](scattering.eps){width="48.00000%"}

To circumvent these problems, the Compton scattering module is instead based on an explicit splitting approach for particles using the calculation of cross sections and allowing for individual weights for each macro particle. To implement the scattering of two macro particles we transform to the rest frame of the target particle and compute the probability $P(n)$ that a single incident particle during a single timestep is scattered on the $n$ target particles. If the incident macro particle has weight $m$, then $k=P(n) m$ particles will interact and two new macro particles representing the fraction of scattered particles are created (see [Fig. \[fig:detailedcollisions\]]{}). In the case of Compton scattering, the target will always be the charged particle, while the incident is a photon, and the scattering amplitude is calculated with the Klein-Nishina formula, which covers the full energy range from Thompson to high energy Compton scattering. To make the process computationally more efficient prior probabilities can be applied. If instead of selecting all pairs in a given computational cell, one only selects pairs with a prior probability $Q$, then the weight of the scattered macro particles has to be changed to $k/Q$. By cleverly selecting the prior probability such that for example the Thompson regime is avoided, the computational load can be greatly decreased.

Each scattering leads to the creation of a new macro particle pair, and if untamed, the number of macro particles will grow exponentially. To keep the number of particles under control, we use particle merging in cells where the number of particles is above a certain threshold. The algorithm is described in detail by Haugb[ø]{}lle[@Haugboelle:2005] and Hededal[@Hededal:2005b].

The Compton scattering module of the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} has allowed us to approach exciting new topics in high energy plasma astrophysics, where plasmas are excited and populations modified by photons, and the back-reaction of the plasma on the (kinetic) photons produces interesting and detailed descriptions of for example the production of inverse Compton components[@Trier:2008.2], and photon beam induced plasma filamentation[@Trier:2008.1].

Coulomb scattering
------------------

Coulomb scattering of charged particles in an astrophysical context is for example important in order to understand the solar chromosphere and the lower parts of the corona, where the mean free path is similar to the dynamical length scales in the system. In the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} we have integrated a collision model, where all macro particle pairs in a cell are considered for scattering. We calculate the scattering process in the elastic Rutherford regime. The collision process is implemented using a physical description for each macro particle pair, starting by calculating the time of closest approach. Only if that time is less than $\Delta t/2$ from the current time, ie. inside the current time interval, is the scattering carried through. In the limit of very small time steps this makes the algorithm independent of the size of the time step. When calculating the impact parameter between two macro particles we have to take in to account that each particle represents a large number of physical particles, therefore the impact parameter is rescaled with the typical distance between each physical particle $n^{-1/3}$, where $n$ is the number density. Because macro particles carry a variable weight, we use the geometric mean of the number density of the particle pairs to calculate the effective number density $n = (n_1 n_2)^{1/2}$.

We consider three different regimes, based on the impact parameter: i) If the impact parameter is larger than the local Debye length we assume that Debye screening between the two particles is so effective that no scattering happens. ii) If the impact parameter is so large that the effective scattering angle is less than $\theta_c$ radians (normally taken to be 0.2 in the code) we use a statistical approach: At small angles the scattering angle is inversely proportional to the impact parameter, and we can replace the many small angle scattering by fewer large angle scatterings, comparing the ratio of the cut-off to the impact parameter $b_c / b$ with a random number. If it is lower the scattering is carried through, but using the cut-off impact parameter $b_c$ for greater computational efficiency. iii) If the impact parameter is smaller than the cut-off parameter $b_c$ then we make a detailed computation of the scattering. In both of the two last cases the scattering angle is calculated in the center-of-mass frame as an inelastic Rutherford scattering, which conservers the energy of each macro particle.

The explicitly physical implementation at the macro particle level of Coulomb scattering is conceptually completely different from the Compton scattering module, in particular because the main consequence of Coulomb scattering is the thermalization and isotropization of the plasma.

Test Problems
=============

Testing the accuracy and precision of a particle-in-cell code is particularly difficult, because of the non linear nature of plasma dynamics, Monte-Carlo particle sampling, and the few examples of realistic test problems with analytic counterparts. To facilitate cross comparison with other codes, below we apply the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} to a set of classic test problems, and in some cases compare different order splines and differential operators to highlight the impact of using high order methods. We also give an example of of the more non-standard feature of radiative cooling. Further tests have been published in the case of Compton scattering[@Haugboelle:2005; @Hededal:2005b] and synthetic spectra[@Hededal:2005b].

Numerical heating and collision artifacts
-----------------------------------------

Numerical heating is a well studied feature of all PIC codes[@birdsall:1985]. It happens due to the interaction of particles with the mesh; the so-called grid-collisions. If particle populations with different energy distributions exist in the plasma, the interaction through the mesh will tend to equilibrate the kinetic energy of each particle species. The equilibration of temperatures happens in a laboratory plasma too, albeit normally at a much slower rate, and the difference between a laboratory and the computational plasma is the much smaller number of macro-particles used to represent the plasma in the latter case. It is also worth pointing out that the numerical mesh heating is an equilibration of kinetic energies, and as such much more severe for ion-electron plasmas. It is important to keep this energy equilibration in mind when simulating plasmas with greatly varying temperatures, or when analyzing heating rates due to real physics. If the numerical heating rate is close to the physical rate in question, the results cannot be trusted. It is interesting to note that the heating rate is practically invariant with respect to the numerical technique used, and instead mainly depends on the number of time steps taken.

The test is done in two dimensions, with three velocity components. There is no bulk velocity, and the kinetic energy corresponds to a thermal velocity of both ions and electrons of $v_{th,e} = v_{th,i} = 0.1\,c$ per component, or $E_\textrm{kin}^\textrm{s} = 0.015\,m_\textrm{s}\, c^2$. The size of the box is $12.8^2$ electron skin depths with 10 cells per skin depth for a $128^2$ resolution. We perform two tests with 5 and 50 particles per cell. The mass ratio is $m_i / m_e = 16$. We used TSC or cubic spline interpolation, 2$^\textrm{nd}$, 4$^\textrm{th}$, or 6$^\textrm{th}$ order field solver, and in the case of the 6$^\textrm{th}$ order field solver with cubic interpolation we also use the charge conserving (CC) method for the current density with second or fourth order time stepping.

![image](electron_temperature_5ppc){width="0.23\linewidth"} ![image](ion_temperature_5ppc){width="0.23\linewidth"} ![image](electron_temperature_50ppc){width="0.23\linewidth"} ![image](ion_temperature_50ppc){width="0.23\linewidth"}

A relativistic cold beam {#sec:coldbeam}
------------------------

A particle-in-cell code is not Galilean invariant. When a cold plasma beam travels through the box at constant velocity the electrostatic fluctuations inside the Debye sphere or, if it is less than the mesh spacing, inside a single cell, will have resonant modes with the mesh spacing. This leads to an effective drag and redistribution of kinetic energy from the stream direction to the parallel direction, and general warm up of the beam. When the temperature reaches a critical level the instability is quenched. For relativistic beams there is the additional complication that electromagnetic waves are represented on the mesh. The solver has an effective dispersion relation, and short wavelength waves travel below the speed of light. On the other hand, particles are Lagrangian, and if relativistic they can effectively travel above the speed of short wavelength waves, giving rise to numerical Cherenkov radiation.

A classic method to limit the impact the of the cold beam instability is applying filters to either the current density or the electric field. This may to some extent filter out the effects, but will also filter out some of the physics. Alternatively, higher order field solvers, interpolation techniques, and time stepping can mitigate the effects.

![image](energy_coldbeam){width="0.46\linewidth"} ![image](stream_vel_coldbeam){width="0.46\linewidth"}\
![image](par_temp_coldbeam){width="0.46\linewidth"} ![image](perp_temp_coldbeam){width="0.46\linewidth"}

To test the numerical methods used in the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}, we have made a cold beam test with 9 different versions of the code. We do not apply any filters to the current density, to show the actual performance of the different code versions. Apart from the 8 methods used for the numerical heating test there is also a version where instead of the implicit field solver a simple (but 6$^\textrm{th}$ order) FDTD explicit solver is used. Notice how the explicit solver has numerical heating, while implicit solvers have numerical cooling[@birdsall:1985], and how the stability of the beam is greatly enhanced by the implicit solver, compared to using an explicit solver. But only the combination of the implicit solver with a charge conserving current deposition gives stable beams, with the smallest heating (see figure \[fig:coldbeam\]). It is also interesting to notice that the heating is non-isotropic. It is therefore not the same to initialize a plasma with a low but stable temperature, as to use a very low temperature and let the cold beam stability warm up the beam.

The test is done in 2D2V with a $\Gamma=10$ streaming pair plasma through a 256 $\times$ 256 cell domain with 10 particles per species per cell and a relativistic skin depth $\delta = [(m c^2 \Gamma)/(4\pi n q^2)]^{1/2}$ of ten cells. The initial temperature is measured in the rest frame of the plasma and has a root-mean-square per velocity component of $0.025c$, so that $\sum \textrm{rms}(v_\textrm{th}) = 0.05c$.

By rerunning at higher resolutions we have investigated at what resolution other methods give comparable results to the charge conserving method with fourth order time stepping, by comparing stream velocities, and parallel and perpendicular temperatures at $\omega_p t=1000$ (see table \[tab:coldbeam\]). It is only at higher resolutions that the non-charge conserving methods do not suffer from the catastrophic instability seen in figure \[fig:coldbeam\], and cost-wise the sixth order charge conserving method is marginally the cheapest. Had the test been in 3D, where the cost goes like resolution to the fourth power, the fourth order time integration method would have been the cheapest for this particular beam test. At high enough resolution the cold beam instability is quenched. For the fourth order charge conserving method this happens at a resolution of approximately 3072$^2$, where instead of a large heating rate, and then a new stable temperature, we observe a gradual heating over time. A resolution of 3072$^2$ corresponds to resolving the Debye length, $\lambda_D = v_\textrm{th} / (\Gamma c) \delta_e$, with 1.9 cells. Taking into account that we do not apply any damping, this is in good agreement with the common wisdom that the Debye length has to be resolved by roughly one cell.

[l c c c]{} Method & Res & Cost & $\mu$s/part\
CCo 6$^\textrm{th}\!$ order field; 4$^\textrm{th}\!$ order time & 256$^2$ & 1 & 6.60\
CCo 6$^\textrm{th}\!$ order field & 380$^2$ & 0.9 & 1.72\
Cubic interpolation; 6$^\textrm{th}\!$ order field & 512$^2$ & 1.3 & 1.05\
TSC interpolation; 2$^\textrm{nd}\!$ order field & 950$^2$ & 4.9 & 0.63\

Relativistic two-stream instability
-----------------------------------

Relativistically counter-streaming plasmas have previously been established to be subject to a general instability class; the oblique (or mixed-mode) two-stream-filamentation instability (MMI), which mixes the two-stream (TSI) and filamentation (FI) instabilities. A thorough and exhaustive analysis of the MMI was given by Bret et al[@Bret:2004; @Bret:2007], and has been investigated numerically by several groups (see e.g. Tzoufras et al[@Tzoufras:2006] and Dieckmann et al[@Dieckmann:2006]). Due to its mixed nature, the MMI contains both an electrostatic and an electromagnetic wave component[@Bret:2004].

Potentially, both electrostatic and electromagnetic turbulence (wavemode coupling leading to cascades/inverse cascades in k-space) is possible in such systems. This potential for producing very broadband plasma turbulence (in both E and B fields) is highly relevant to inertial confinement fusion experiments. Other examples where electromagnetic wave turbulence[@Frederiksen:2008] and the general MMI are important are astrophysical jets and shocks from gamma-ray bursts and active galactic nuclei, where ambient plasma streams through a shock interface moving at relativistic speeds.

To test physical scenarios responsible for observational signatures from these astrophysical sources, and from plasma experiments, we must construct plausible shocked outflow conditions[@bib:Medvedev1999; @Frederiksen:2004; @haugboelle:2011] and then subsequently synthesize radiation signatures[@Hededal:2005b; @Trier:2010; @Medvedev:2011] to test the assumed physical conditions against observational evidence. Studying the MMI, in both its linear and non-linear evolution is therefore well motivated.

Growth rates of the (general) MMI, the special case of FI and the TSI, respectively, are calculated as in[@Bret:2004]: $$\begin{aligned}
    \gamma_{MMI} & = & 2^{-4/3}~ \sqrt{3}~ \alpha^{1/3}~ \Gamma(v_b)^{-1/3} \label{eq:mmigrowth}\\
    \gamma_{FI}  & = & \beta~ \alpha^{1/2}~ \Gamma(v_b)^{-1/2}~             \label{eq:figrowth}\\
    \gamma_{TSI} & = & 2^{-4/3}~ \sqrt{3}~ \alpha^{1/3}~ \Gamma(v_b)^{-1}   \label{eq:tsigrowth}\end{aligned}$$ where $\beta \equiv v_b/c$ is the beam velocity, $\alpha \equiv n_b/n_p$ is the beam-to-background density ratio, and $\Gamma(v_b)$ is the beam bulk flow Lorentz factor. For thin, high-Lorentz factor beams, the MMI is the fastest growing mode, dominating over both the FI and the TSI. In the test we have chosen the MMI is dominant, with subdominant FI and TSI components.

We perform six runs with identical initial conditions and physical scaling, using combinations of finite difference operators and particle shape functions as given in Table\[tab:orders\] below. This way, we test the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} for differences/similarities between interpolation schemes.

To capture the MMI as the fastest growing mode, we initialize a simulation volume with a cold thin neutral beam (electrons + ions) through a warm thick neutral background (electrons + ions), with no fields initially. The beam and background densities are $n_b=0.1$ and $n_p=0.9$, respectively. The beam velocity is chosen to have $\Gamma(v_b)=4$. Temperatures of the beam and background are $T_b=0.01$ and $T_p=0.1$, respectively. With these choices of physical properties, the growth rates of the fastest growing MMI, FI, and TSI modes become $\gamma_{MMI}=0.201$, $\gamma_{FI}=0.153$, $\gamma_{TSI}=0.080$.

The computational domain is $\{L_x,L_y,L_z\}=\{12.8\delta_e,12.8\delta_e,12.8\delta_e\}$, with $\{N_x,N_y,N_z\}=\{128,128,128\}$ cells. We use 20 particles/cell/species, or a total of 80 particles/cell (beam+background). The physical constants are scaled as $c=1$, $q_e=1$, $m_e=1$ and $m_i/m_e\equiv1836$.

[cccccc]{} Run & Fields & Particles & charge-conservation & Time order\
1 & 2$^\textrm{nd}$ & tsc & no & 2\
2 & 2$^\textrm{nd}$ & cubic & no & 2\
3 & 6$^\textrm{th}$ & tsc & no & 2\
4 & 6$^\textrm{th}$ & cubic & no & 2\
5 & 6$^\textrm{th}$ & cubic & yes & 2\
6 & 6$^\textrm{th}$ & cubic & yes & 4\

For our choice of run parameters a MMI mode develops with a propagtion wave vector, ${\bf{k}}_{MMI}$, that is oblique with respect to the streaming direction, at an angle given by $\theta_{MMI}=\angle({\bf{k}}_{MMI},{\bf k}_\textrm{beam})=
\textrm{arctan}\left(\sqrt{v_b/v_p}\right)\approx74^\circ$. The corresponding electric field component is almost parallel to the direction of propagation[@Bret:2004].

To find growth rates of the MMI we calculate the volume integrated electrostatic energy as a function of time $$E_{E,tot}(\theta_{MMI}) = \int_{V} |{\bm{E}}_\perp| sin(\theta_{MMI}) + E_\parallel cos(\theta_{MMI})~dV$$ of the electric field projected on the propagation direction, ${\bm{E}}({\bm{r}})\cdot{\bf{k}}_{MMI}$. Similarly, we also measure the TSI mode growth rate ([Eq. \[eq:tsigrowth\]]{}) by the same calculation, but for the TSI we have $\theta_{TSI} = \angle({\bf{k}}_{TSI},\hat{\bf{z}}=0)$. Our results are summarized in Table \[tab:growth\].

[ccccccccc]{} Run & 1 & 2 & 3 & 4 & 5 & 6 & theory\
$\gamma_{MMI}$ & 0.185 & 0.190 & 0.203 & 0.206 & 0.205 & 0.205 & **0.201**\
$\gamma_{TSI}$ & 0.060 & 0.079 & 0.065 & 0.082 & 0.088 & 0.088 & **0.080**\

From figure \[fig:mmi\_growth\] we see that the initial noise build-up prior to instability dominance is strongest (as expected) for Run1. This results in lower growth rates for lower order runs, since the noise tends to ’flatten« the total energy history, with a lower $\gamma_{MMI}$ as a result. This error in the measurement decreases with increasing run number (Run1$\rightarrow$Run2$\rightarrow$…), with higher order runs’ growth rates less susceptible to noise distortion. This is similar to what was seen in the cold beam tests: in a lower order integration scheme more energy is lost to artificial heating and Cerenkov radiation, making less energy available for the physical instabilities, and changing the parameters (for example temperatures) of the plasma components, modifying the overall setup.

Because of the lower level of noise, the onset of instability is delayed for runs of increasing order, while the peak energies and peak times coincide roughly for all runs; this effect correlates with the increased growth rates of the higher order runs, relatively. Both of the effects mentioned above are caused by the difference in dynamic response of the various schemes. A higher order scheme is desirable since the noise levels, heating, numerical stopping power and dynamical friction are all reduced considerably. The differences between using either TSC or cubic interpolation, second or sixth order field integration, and (no) charge conservation are all clearly reflected in the growth rates of the MMI and TSI.

Nonetheless, growth rates are seen to converge with different methods to a value which deviates from the theoretical prediction by about 2% for the MMI component, and $\sim$10% for the TSI component, for the highest order schemes, and the overall development is in qualitative agreement for all test cases.

![Growth of the volume-totaled electrostatic energy, for the electric field projected onto the direction of propagation of the fastest growing mode, ${\bm{E}}\cdot {\bf{k}}_{MMI}$. Runs 1-6 are compared in the plot. Thickened line segment Run1 (black) gives the fitting interval. Runs have decreasing initial energy for increasing run number designation (Table \[tab:orders\]). Runs 5 (orange) and 6 (yellow) are completely coinciding.[]{data-label="fig:mmi_growth"}](mmi_growth){width=".48\textwidth"}

Concluding this test, we have verified the growth rate of the relativistic two-stream (or oblique or mixed-mode) instability, and found very good agreement with growth rates also when selecting the TSI branch. The slight excesses in values of $\gamma_{TSI}([Run1,...Run6])$ is likely caused by the fact that the relativistic beam is perfectly grid aligned, thus subjected to the finite grid instability, which introduces additional electric field energy in the beam direction, while for lower orders, this is more than compensated by the overall dissipation to all electromagnetic and particle components.

Radiatively cooled collisionless shocks
---------------------------------------

The radiative cooling currently implemented is inspired by the early work of Hededal; he validated it for a simple test case of a radiatively cooling charged particle in a homogeneous magnetic field[@Hededal:2005b]. For a more non trivial application we here for the first time present a series of simulations of radiatively cooled initially non-magnetized collisionless shocks. If we assume that the acceleration of a particle in a collisionless shock is mostly due to a homogeneous magnetic field $B$, we can estimate the synchrotron cooling time for an initial Lorentz factor $\gamma_0$ as $$t^\textrm{syn}_\textrm{cool} = \frac{6 \pi \epsilon_0 m^3 c^3}{q^4 B^2 \gamma_0}$$ Consider a relativistic collisionless shock in the contact discontinuity frame, with the upstream moving at a Lorentz factor $\Gamma$ and having a number density $n$. The kinetic energy density in the upstream is $$E_{kin} = (\Gamma - 1) n \sum m c^2 = (\Gamma - 1) M c^2 n\,.$$ We assume that the magnetic energy density at the shock interface $E_B = B^2 / 2 \mu_0$ is related through some efficiency factor $\alpha\simeq0.1$ to the upstream kinetic energy density. The relativistic plasma frequency in the upstream medium is $\omega_{pe}^2 = q^2 n / (\epsilon_0 m_e \Gamma)$. Putting it all together we can find the synchrotron cooling time at the shock interface, for a particle with a gamma factor $\gamma_0$, to be $$T_\textrm{cool} = \omega_{pe} t^\textrm{syn}_\textrm{cool} =
  \frac{3 \pi^{1/2} \epsilon_0^{3/2} m_e^{5/2} c^3}{\alpha q^3 n^{1/2} \Gamma^{1/2} (\Gamma - 1) M \gamma_0} \,,$$ For a pair plasma $M=2m_e$, and if we consider upstream particles $\gamma_0=\Gamma$ it reduces to $$T_\textrm{cool} =
  \frac{3 \pi^{1/2} \epsilon_0^{3/2} m_e^{3/2} c^3}{2 \alpha q^3 (\Gamma - 1) \Gamma^{3/2} n^{1/2}}\,,$$ We use this definition to label the different runs. Our definition of the cooling time is similar to the one given in Medvedev and Spitkovsky[@Medvedev:2009], but differs because they considered the downstream skin depth: $T^\textrm{syn}_\textrm{cool,our} = 9/4 \Gamma^{1/2} T^\textrm{syn}_\textrm{cool,M-S}$.

![image](epsilon_B){width="0.48\linewidth"} ![image](rcool_phase_space){width="0.42\linewidth"}

We have made long term 2D2V simulations of the shock, using reflecting boundaries and different cooling times, including a run without cooling for reference, and with electron-ion and pair plasmas. The runs were done with cubic interpolation, sixth order field pusher and a 17-point current density filter. We used the sliding window technique to be able to follow the development of the shock up to ${\omega_{pe}\,t}=5000$, and used more than 5 billion particles to model the shocks. In all cases the upstream Lorentz factor is $\Gamma=10$. While below we present results from runs with 20 cells per skin depth, we have used both 10, 20, and 40 cells per skin depth and between 12 and 24 particles per species per cell and find converged results. Examples of the magnetic field and phase space density for different cooling times can be seen in figure \[fig:cooling\]. The shocks can phenomenologically be categorized into three types: shocks in the strong cooling regime, radiative shocks, and weakly radiative shocks. In the strong cooling regime, there are no traces of a power law tail of accelerated particles, and both the evolution, shock jump conditions, and the micro structure itself are qualitatively different from a normal non-radiating shock. In particular, the upstream is completely unperturbed by the existence of the shock. This is found for cooling times $T_\textrm{cool} \lesssim 200$. In a strongly radiating shock the micro structure is disturbed, and the downstream magnetic islands only survive very close to the shock interface, but the upstream does have resemblance to a normal collisionless shock. The power law tail of accelerated particles is completely gone, and the shock jump conditions are altered. This is found for cooling times $200 \lesssim T_\textrm{cool} \lesssim 1000$. The weakly radiative shocks are similar in structure to a non-radiating shock, but with mildly perturbed shock jump conditions and propagation velocity (see figure \[fig:jump\]). The power law index for the high energy population of accelerated electrons is steeper than for a non-radiating shock, with a dependence on the cooling time. The magnetic energy density decreases significantly approximately 150 $\delta_e$ away from the shock interface, but the time a high energy particle spends that close to the shock transition, where the bulk of the cooling happens, is a stochastic function of its angle to the shock normal, and the number of scatterings on fluctuations in the electromagnetic field. In principle, if we simulate for a long enough time, with a long cooling time and with high statistics, the powerlaw tail of accelerated particles will grow over more than a decade in energy, and a well defined cooling break should emerge with two different slopes clearly visible. In practice, given the tangled nature and stochastic propagation, the break will most probably be smooth, significantly smeared out around the characteristic energy, where electrons start to be efficiently cooled. In these exploratory simulations the emergence of a cooling break does not occur. Instead, in the case of weakly radiative shocks, the high-energy part of the particle distribution (PDF) is a powerlaw with an exponential cut-off, but with a steeper powerlaw index than in the non-radiating case (see fig. \[fig:pdf\]). It is well known that in PIC simulations of non-radiating shocks the upstream region affected by high energy particles produced at the shock interface only grows with time, and it has been an open question what the long term structure looks like[@Keshet:2009]. This is different for radiatively cooled shocks, where at large times the shock settles down to a quasi-steady state, making it possible to draw conclusions about the long term behavior. The extent of the upstream is limited, and the powerlaw part of the PDF does not grow in time. Any collisionless shock is radiatively cooling, given long enough time, and from our simulations it is clear that the impact of cooling for weakly radiating shocks is greater than speculated in e.g. Medvedev & Spitkovsky[@Medvedev:2009], where analytic estimates were used. Given the possibility of collisionless shocks to mediate secondary instabilities such as the Bell instability far upstream, by the generation of streaming cosmic rays[@Niemiec:2008], it would be interesting to understand the impact on collisionless shock for the very large cooling times expected in e.g. GRB afterglows. We speculate that using progressively larger cooling times in sufficiently large simulations could give insight in how to scale the solution to arbitrarily long cooling times.

![Top panel: Density ratio between up and down stream at ${\omega_{pe}\,t}= 1600$, as a function of cooling time. Bottom panel: Effective adiabatic index, derived from the shock velocity. The right most point is for a run without radiative cooling.[]{data-label="fig:jump"}](jump_conditions){width="\linewidth"}

![The particle distribution function sampled in the downstream region for the case of no cooling and $T_\textrm{cool}=25600$, firmly in the weakly radiative regime. The power law index is indicated with the dashed line, and given in the legend.[]{data-label="fig:pdf"}](pdf_nocool5_c25600){width="\linewidth"}

The linear magnetized Kelvin-Helmholtz instability
--------------------------------------------------

Earth’s magnetopause constitutes an important region in space. It is the boundary layer, separating the Earth’s magnetosphere from the solar wind. In this region, the Kelvin-Helmholtz instability (KHI) is driven by velocity shears between the magnetosheath and the magnetospheric plasma at low latitude (close to the magnetic equator). We have selected the problem of the linear magnetized KHI to compare with results using our local MHD [<span style="font-variant:small-caps;">Stagger</span> code]{}. Essentially, we use the setup for the SWIFF Magnetopause Challenge [@SWIFF_D2.1_KH] code comparison, except that the amplitude is slightly different and there are half the number of skin depths in the kinetic case. The MHD and PIC code setups are nearly identical, the only difference being the addition of microphysical parameters, and that in the PIC case we must ensure that the initial condition is a kinetic equilibrium respecting constraints like Gauss law.

We compare our PIC results against those obtained with the [<span style="font-variant:small-caps;">Stagger</span> code]{}, also developed and maintained at the Niels Bohr Institute@1997LNP...489..179N. This finite difference mesh based MHD code is a fully 3D resistive and compressible MHD code. The MHD variables are located on staggered meshes, and the discretization is very similar to the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}, with sixth order spatial derivatives, and fifth order interpolation of variables. The time integration of the MHD equations is performed using an explicit 3rd order low storage Runge-Kutta method[@1980JCoPh..35...48W].

The experiment is a periodic 2D3V setup, and the box size is $L_x=90\pi$, $L_y=30\pi$. The initial velocity field $\textbf{V} = V_y(x) \textbf{e}_y$ contains a periodic double shear layer, to avoid boundary effects, with a velocity amplitude $A_{eq}=1$. The sheared velocity jumps are located at $L_x / 4$ and $3/4 \ L_x$, and the transition width is $a=3$. With these parameters, the initial velocity profiles is defined as: $$V_y(x) = \tanh \left( \frac{x - L_x / 4}{3} \right) - \tanh \left( \frac{x - 3/4 L_x}{3} \right)-1~.$$

Initially, $J_{eq} = {\bm{E}}= 0$, ${\bm{B}}= B_0 sin(\theta) \textbf{e}_y + B_0 cos(\theta)  \textbf{e}_z$, with $B_0 = 1$, $\theta = 0.05$. The density, pressure and Alfvén velocity are all unity: $\rho=P=V_A=1$. For the MHD we use an adiabatic equation of state with $\gamma_\textrm{gas} = 5/3$, and impose a small perturbation in the velocity field to seed the KHI $\delta\textbf{V}=\textbf{e}_z\times\nabla\psi$, where $$\psi = \epsilon f(x) \sum_{m=1}^{N_y/4} cos(2 \pi m y / L_y + \phi_m )/m ~,$$ and $$f(x) = \exp \left[- \left(\frac{x - L_x/4} {3}\right)^2 \right] + \exp \left[ - \left(\frac{x - 3/4L_x} {3}\right)^2 \right]$$ and $\epsilon$ is such that $\max(| \delta V |) \simeq 10^{-3}$. $\phi_m$ are random phases.

In the MHD case we resolve the box with $N_x=1536$, $N_y=512$ cells. The PIC setup was prepared with physical initial conditions essentially identical to the MHD setup, using the technique of section \[sec:embed\] to reach an approximate kinetic equilibrium. The objective was to produce the ion-scale KHI, while resolving the electron skin-depth $\delta_e$ and having well separated inertial scales with $m_i/m_e =64$. We use a setup similar to what was used for a hybrid code in the SWIFF comparison[@SWIFF_D2.1_KH], but limit the number of ion skin depths to $45\pi \times 15\pi$, and the resolution to $N_x=6144$, $N_y=2048$. To have a reasonable plasma frequency, we rescale the speed of light to $c=10$, and use $\delta_e = 6 \Delta x$, and 20 particles per cell per species with a total of roughly 500 million particles in the box. Our choice of the shearing jump amplitude ($V_0=\pm1.0$) and width ($a=3.0$) selects a fastest growing mode (FGM) leading to production of two vortices, which pair up and merge during the early and late non-linear stages of the KHI. Figure \[fig:small\_PIC\_KH\_fig\] shows the vortices just prior to, and well after the vortex merging in the PIC model.

![image](test.eps){width="\textwidth"}

The growth rates for the two extremum cases of transverse (${\bm{B}}_0 \perp {\bm{v}}_0$) and parallel (${\bm{B}}_0 \parallel
{\bm{v}}_0$) configurations were investigated theoretically and further calculated by Miura[@Miura:1982.1; @Miura:1982.2]. We may assume that for a weak parallel component, i.e. $B_{0z} \ll B_{0y}$ ($B_{0z} = 0.05 B_{0y}$), the instability evolves almost as the ideal transverse case to a good approximation. It was also determined that the stability criterion for $V_0=2$ is $M_f\equiv{V_0}/{\sqrt{v_A^2+c_s^2}}<2$. Here, $M_f$ is the fast mode magneto-acoustic Mach number.

For our setup $M_f=\sqrt{2}$, with both $v_A^2=1$,$c_s^2=1$, using $V_0=2$ and $a=3$. From Miura 1982[@Miura:1982.1], we find a growth rate for our specific setup to be ${2a}\gamma_{KH}/{V_0}\approx0.162$, which leads to the final result of $$\gamma_{KH,FGM}\approx0.055\pm0.002.$$

To measure the KHI rate growth — in both the MHD and PIC cases — we calculate the quantity $$\left|\tilde{V_x}(x_i,k_y,t_n)\right|^2 = \left| \frac{1}{N_y} 
    \sum^{N_y-1}_{j=0} V_x(x_i,y_j,t_n) e^{-i2\pi\frac{k_yj}{N_y}} \right|^2~,$$ i.e. the power spectra along the $y$-axis of the $x$-component of velocity, $V_x$, at constant $x_i$. This is then averaged in the $x$-direction $$\mathcal{Q}(t_n) = \frac{1}{x_{hi}-x_{lo}} \sum^{x_{hi}}_{x_{lo}} \left|\tilde{V_x}(x_i,k_y,t_n)\right|^2~,$$ for two different sets of $\{x_{lo},x_{hi}\}$, namely i) centered on one shearing layer, 5$\Delta_x$ wide, and ii) across the entire half-volume in the x-direction. We then make a fit to exponential growth to obtain $\gamma_{KI}(mode)$ for the FGM.

[ l c c c ]{} & PIC &   & MHD\
Right & 0.058 &   & 0.058\
Total & 0.050 &   & 0.056\

We find that growth rates from both the MHD and PIC simulations agree well with theory, to within about 10% (see table \[tab:kh\_growth\_rates\_compare\]). Averaging over all layers in the simulation half-plane yields a lower growth rate with the MHD in closest agreement with theory. This is expected, since averaging over only 5 layers does not capture the entire width of the shearing layer which is $W_{shear}\sim50\Delta_x$. For the PIC growth rate results discrepancies in physics and uncertainties in data fitting, due to mode coupling are higher than for the MHD case; and growth rates therefore differ by as much as 10% in the in the volume-averaged case. The explanation for a slightly lower growth rate in the total volume averaged PIC case may be due to enhanced dissipation and intrinsic noise properties of the PIC code, or the development of secondary instabilities.

Concluding this test, we emphasize that a PIC code has been used to run a fully MHD problem in PIC explicit mode, with almost identical growth rates and vortex structures in the linear phase. For the vortex-merging epoch of the PIC run, results deviate qualitatively and quantitatively from MHD, due to kinetic effects and secondary tearing-like instabilities developing during the merging stage.

Parallelization, scaling and performance
========================================

Modern 3D Particle-in-Cell experiments in astrophysics use billions of particles to simulate the macroscopic structure of plasmas. To run these simulations the code has to be massively parallel. The [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} started with a simple domain decomposition along one axis, using MPI. Later it was changed to support a 3D MPI decomposition, and in the current version we use a hybrid parallelization with OpenMP and MPI to scale the code effectively up to hundreds of thousands of cores (with 262.144 cores being the largest case tried so far). The OpenMP hybrid approach is also necessary when running the GPU version of the code, because normally the ratio of CPU cores to GPUs in a GPU system is larger than one.

MPI
---

The MPI parallelization has been designed to be as transparent as possible, and consists of different modules, all collected in a single file. This makes it easy to compile the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} both with and without the MPI library.

-   [Particles: After moving the particles on a node, and applying the physical boundary conditions, on each thread we check sequentially in the $x$-, $y$-, and $z$-direction which particles have moved to other threads, and interchange particles accordingly. A reverse transversion of the linked list structure containing the send particles makes it efficient to store the received particles in the same slots.]{}

-   [Fields: When the physical boundary conditions are applied for the fields, ghostzones are also exchanged between threads.]{}

-   [I/O: The snapshot and restart file format is binary, and the code uses MPI-IO. For the particles each attribute (ie. the $x$-position) is stored in a single block, and the IO-types are all 64-bit, making it possible to store billions of particles in a single snapshot with GB/s of performance. The code stores everything in a few files, and restarts can be done on an arbitrary number of threads.]{}

-   [Synthetic Spectra: The spectra are sampled on sub-volumes of the data, and each thread only allocates data for the sub-volumes that intersect with the local domain. This makes the mechanism scalable to $\mathcal{O}(10^5)$ sub-volumes, without wasting memory. The lowest ranked thread for a given sub-volume writes the data, making the I/O scalable too.]{}

-   [Particle Tracing: Particle tracing is done locally on each thread, while only the master thread writes the I/O. This is a potential bottleneck, and on systems with relatively weak CPUs, e.g. Blue-Gene/P we are limited to tracing $\sim$1 million particles, while for x86-based clusters with global Lustre filesystems we can effectively trace up to $\sim$10 million particles, without significant slowdowns.]{}

-   [2D-slices: The code can extract (averaged) 2D slices on-the-fly while running. While this is very convenient for movie making, for large runs thread contention becomes a problem when reducing the data on a single thread: When thousands of threads try to communicate with a single thread, the slow down can be significant. We use single-sided MPI to effectively circumvent the thread contention, making slice writing scalable to at least 100.000 MPI threads on network transports that support RDMA.]{}

A major problem with simple domain decomposed particle-in-cell codes is that systematic local over-densities of particles easily occur, for example in collisionless shock or laser wake-field simulations. This leads to severe load imbalance. To work around this the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} has a simple dynamic load-balance feature. If enabled, mesh-slices along the $z$-direction can be exchanged in-between neighboring threads, according to a user defined cost formula. In the current version particles and cells are assigned a certain cost, and this makes it possible to dynamically sample local spikes in the particle density (see Figure \[fig:loadbalance\]). When a slice of cells is moved from a thread to the neighbor thread the corresponding field values, particle data, synthetic spectra sub-volumes, and local random number generators are updated. The same mechanism is used to make an optimal distribution of cells when restarting a run. We have tested the MPI performance in a simple, but realistic setup, with two counter streaming plasma beams on the Blue-Gene/P machine JUGENE running in pure MPI mode. Figure \[fig:mpi\_scaling\] shows weak scaling behavior of the code from 8 to 262.144 cores. This was done a few years ago, with a version of the code without charge conservation, and the scaling of the current code is significantly better.

![Density profile of a 2D collisionless shock. The dashed lines indicate the $z$-boundaries of MPI domains. They are updated dynamically to equalize the load.[]{data-label="fig:loadbalance"}](density){width="\linewidth"}

![Weak scaling in pure MPI mode when running on the BG/P machine JUGENE. The setup is a relativistic two-stream experiment with periodic boundaries and a $\Gamma=10$ streaming motion. There are roughly 80 particles per cell with $16^3$ cells per MPI domain. This experiment is not using the charge conserving current deposition, and the scaling of the current code is significantly better.[]{data-label="fig:mpi_scaling"}](mpi_scaling){width="\linewidth"}

![image](openmp_scaling1){width="0.32\linewidth"} ![image](openmp_scaling5){width="0.32\linewidth"} ![image](openmp_scaling3){width="0.32\linewidth"}

OpenMP
------

The 3D MPI domain decomposition has served well, scaling the code, for most applications, to 10.000 cores on x86 clusters and more than 100.000 cores on Blue-Gene/P. To scale the code to a million threads or more on for example Blue-Gene/Q, one must introduce a new layer of parallelization. We have used OpenMP for a number of reasons: It is relatively easy to get almost perfect scaling for a low number of cores, it can be done incrementally, and it is a natural thing to do for the GPU version of the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{}. Furthermore, using a hybrid parallelization, the size of each MPI domain becomes bigger, and the particle load-balance requirements, due to fluctuations and variations in number density, become smaller. This is important for experiments with some level of density fluctuations, even if the scaling per se is good in pure MPI mode. In the current version the most CPU consuming parts of the code have been OpenMP parallelized:

-   [The Mover / Charge deposition is trivially parallelized by allocating one set of charge density fields per thread, and parallelizing the update of particles on a per cell basis.]{}

-   [The field solver consist of simple differential operators, and we use loop-based parallelization for perfect speedup.]{}

-   [The sorting of particles is more challenging: We first parallelize on the number of species, this is embarrassing parallel. Then on a nested level we first partition the particle data with nested parallel sweeps in 128 sets, and afterwards use quick-sort to sort each set of particles in parallel.]{}

-   [Sending particles between threads: We first parallelize on the number of species, this is embarrassing parallel. Then on a nested level we use loop parallelization to select particles to be send, and multiplex MPI communication with OMP sections.]{}

-   [Several other auxiliary routines have been parallelized and made thread safe: The random generators, synthetic spectra sampling, initial- and boundary conditions.]{}

The rest of the code, mostly diagnostics and I/O, can be incrementally OpenMP parallelized as needed. We find excellent scaling inside x86 nodes, with optimal performance at 4 MPI threads per socket / 8 per node, only loosing about 5% when using 1 MPI thread per socket. The 5% is easily regained in higher MPI efficiency and better load-balance when running on a large number of nodes (top panel in [Fig. \[fig:openmp\]]{}). We also find that we gain roughly 15% in performance by running with hyperthreading (32 compared to 16 threads in bottom panel in [Fig. \[fig:openmp\]]{}). Finally, on Blue-Gene/Q it is absolutely crucial to use all 64 hardware threads on a 16 core node. It improves overall performance by a factor of 2. But to effectively use the massive amount of threads exposed in the system, we need OpenMP, since otherwise the domain size is too small. Currently we use 4 or 16 OpenMP threads per node, depending on the scale of the problem. Even though we loose 16% performance by going from 4 to 16 OpenMP threads, for runs with e.g. particle tracing enabled or other IO done only through the master thread, and for very large runs, where load imbalance can be more than 20%, it is advantageous to use 16 OpenMP threads per MPI thread. In simpler and smaller scale runs we use 4 OpenMP threads per MPI thread.

Concluding Remarks
==================

In this paper we present the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} particle-in-cell code, with new numerical methods, physical extensions, and parallelization techniques. Originally the main motivation for developing a new particle-in-cell code from scratch was to implement a modern, modularized and extendible code, with modern numerical techniques. In the paper we have demonstrated how the extension of the classical PIC framework to higher order spatial and temporal derivatives, together with a novel charge conservation scheme, gives a much higher accuracy and better stability than traditional, second order PIC codes, when using the same number of grid points. Conversely, fewer mesh points are needed to reach a given level of fidelity in realistic three dimensional kinetic astrophysical setups, which leads to large savings in computational costs. The excellent scalability of the code, and the ability to do diagnostics on-the-fly, makes it possible to achieve a high scientific turnaround, where numerical experiments with more than 100 billion particles can be performed in a day, with dedicated use of petascale computational resources.

Among the novel physical extensions of the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} code are the description of binary interactions by particle splitting, the ability to include radiative cooling in a self-consistent manner, the possibility of embedding the kinetic model in an MHD snapshot, and the application of a sliding simulation window. Other extensions are in development or are planned for the future, among them nested grids, proper inclusion of neutral particles and dust, self gravity, smooth transition from the kinetic to the MHD regime, and simpler binary collision drag descriptions, relevant for studying fractionation processes.

The [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} has already been in full production mode for some years, the user base is growing steadily, and we expect that it will be used extensively at the Niels Bohr Institute and elsewhere for many years to come. Hopefully the detailed description of the algorithms documented in this paper will also be of use to others, when developing particle-in-cell techniques and applying them astrophysical plasma physics problems in the future.

We thank Juri Poutanen for funding the workshop ’Thinkshop 2005’ in Stockholm 2005; <http://comp.astro.ku.dk/Twiki/view/NumPlasma/ThinkShop2005>, which initialized the development of the [<span style="font-variant:small-caps;">PhotonPlasma</span> code]{} code. Further, we wish to acknowledge Boris Stern for ample help with Compton scattering Large Particle Monte-Carlo techniques. TH and [Å]{}N are supported by the Centre for Star and Planet Formation, which is financed by the Danish National Science Foundation. [Å]{}N and JTF are supported by SWIFF; the research leading to these results has received funding from the European Commission’s Seventh Framework Programme (FP7/2007-2013) under the grant agreement SWIFF (project n$^\texttt{o}$ 263340, www.swiff.eu). Computer time was provided by the Danish Center for Scientific Computing (DCSC), and through the PRACE project ’Ab Initio Modeling of Solar Active Regions’ running on JUGENE and JUQUEEN at the Jülich Supercomputing Centre; we thank the centre staff for their support, in particular in connection with I/O scaling issues, and for help with execution of full machine jobs.
