## 1. Introduction

Particle-based plasma simulations, including particle-in-cell (PIC) codes, are ubiquitous as a computational tool to understand plasma evolution, aiding research into nuclear fusion (Roth *et al.*
Reference Roth, Cowan, Key, Hatchett, Brown, Fountain, Johnson, Pennington, Snavely and Wilks2001), proton radiography (Borghesi *et al.*
Reference Borghesi, Sarri, Cecchetti, Kourakis, Hoarty, Stevenson, James, Brown, Hobbs and Lockyear2010) and hadron therapies (Bulanov & Khoroshkov Reference Bulanov and Khoroshkov2002). PIC codes with the ability to model particle–particle Coulomb collisions have become commonplace (Nanbu Reference Nanbu1997; Sentoku & Kemp Reference Sentoku and Kemp2008; Peano *et al.*
Reference Peano, Marti, Silva and Coppa2009; Schmitz, Lloyd & Evans Reference Schmitz, Lloyd and Evans2012; Bobylev & Potapenko Reference Bobylev and Potapenko2013), and quantum electrodynamics is the latest area of physics that they have begun to include (Di Piazza *et al.*
Reference Di Piazza, Müller, Hatsagortsyan and Keitel2012).

The addition of more physical phenomena means that increasingly complex interplays between different physics drives the overall evolution of plasmas in PIC simulations. In many applications of PIC, the objective is to understand the important processes behind an experimental result by running the simulations with similar conditions. Unravelling the exact interplay which causes particular effects to occur requires analysing the large amounts of particle-level data which are output by PIC simulations.

Here, we present a method to aid the analysis of the collisional interactions present in a PIC simulation. Usually, the collisional terms in the Boltzmann equation can only be evaluated from PIC output by integration over a six-dimensional phase space in each cell, which requires careful manipulation of often noisy distribution functions. These collisional terms are $Q_{s}$ and $\boldsymbol{R}_{s}$ , which describe the exchange of energy density in time by species $s$ and the exchange of momentum density in time by species $s$ , respectively. We provide a method to evaluate these terms using a simple sum over the properties of each particle in each cell. This simple sum can, in many case, be computed far more efficiently than the integration of the six-dimensional phase space.

This paper is organised as follows: § 2 sets out the advantages of this approach; § 3 defines the Boltzmann equation and its moments; § 4 states the main result of the work; § 5 contains the derivation of the results; and § 6 presents a computational benchmark of the results.

## 2. Merits of the approach

Evaluating the collisional terms $Q_{s}$ and $\boldsymbol{R}_{s}$ defined in § 3 directly from distribution functions as output by PIC simulations requires two steps which can introduce errors. The first is in calculating continuous, smooth distribution functions from the properties of a discrete population of particles. This inevitably involves a binning procedure with a trade-off between resolving features in the distribution function and accepting noise in the value of $f_{s}$ as a function of $\boldsymbol{v}$ . The proposed method requires no binning. The second involves performing the integration over the six-dimensional phase space. As a six-dimensional integration is a computationally intensive procedure, lower-order integration methods which admit more numerical error may be used. In the equations specified, only a direct sum over particle properties is required to evaluate $Q_{s}$ and $\boldsymbol{R}_{s}$ .

The extent of the computational advantage of directly evaluating the collisional terms, as opposed to performing the integration over phase space, depends upon the number of particles per cell and the binning procedure. A binning procedure with $B$ bins per momentum dimension for the particle distribution functions will require at least $O(B^{6})$ computer operations in each cell for the evaluation of the collisional terms. In the equations described in § 4, evaluating the collision terms for $N$ particles per cell requires $O(N^{2})$ computer operations. Given a typical $B$ of 100, the direct evaluation method is more computationally efficient whenever $N<B^{3}\approx 10^{6}$ particles per cell.

The computations necessary to calculate $Q_{s}$ and $\boldsymbol{R}_{s}$ could also be performed at run-time, as part of the operations of the base code. In this case, the energy and momentum transferred in each collision can be recorded in-line before and after each collision. However, widely used PIC codes do not already include this feature, and it would be an extraneous drain on computer and memory resources for simulations in which particle kinetics were not specifically of interest. Furthermore, most PIC modellers are end-users who may prefer not to modify the base code in this way. The proposed equations require only the output particle data with no modification of the base code; they can be used in post-processing. They also involve only a single evaluation per pair of particles (as opposed to a before and after evaluation for each particle-pair collision in-line).

## 3. Boltzmann equation

The kinetic equation describing plasma evolution is the Boltzmann equation given by

where $f_{s}$ is the distribution of particles of species $s$ defined by $n_{s}(\boldsymbol{x}_{s},t)=\int \,\text{d}^{3}v\,f_{s}(\boldsymbol{x}_{s},t,\boldsymbol{v}_{s})$ , with $n_{s}$ the density of species $s$ . Here $C(f_{s})=\sum _{s^{\prime }}C(f_{s},f_{s^{\prime }})$ is the collision operator. The labels $s$ and $s^{\prime }$ refer to different species. Temperatures are expressed in units of energy throughout unless stated explicitly. Henceforth, Latin letters will refer to particles belonging to a particular species so that $i$ is a particle of species $s$ and $j$ is a particle of species $s^{\prime }$ . Greek letters refer to vector and matrix indices of rank three, i.e. ${\bf\nu},{\bf\mu},$ etc., range over $1,2,3$ . In addition, bold font denotes a vector. Derivatives are denoted by

*a*,

*b*) $$\begin{eqnarray}\frac{\partial }{\partial v_{{\it\mu}}}\equiv \partial _{{\it\mu}},\quad \frac{\partial }{\partial v_{{\it\mu}}^{\prime }}\equiv \partial _{{\it\mu}}^{\prime },\end{eqnarray}$$

where primed and unprimed relate to species $s$ and $s^{\prime }$ respectively.

We utilise the Landau collision operator, which is given by

with collision kernel

where

can be assumed to be constant including the Coulomb logarithm $\ln {\it\Lambda}$ , and

The relative velocity vector is given by $\boldsymbol{u}=\boldsymbol{v}-\boldsymbol{v}^{\prime }$ .

The Landau collision operator is a small-angle collision-only approximation to the full Boltzmann collision operator, meaning that it is obtained by an expansion in ${\rm\Delta}\boldsymbol{v}$ and therefore ignores collisions with large scattering angles (Alexandre & Villani Reference Alexandre and Villani2004). However, the same approximation is used by most, if not all, PIC simulation collision algorithms, though methods to re-incorporate the collisions lost through this method are available (Turrell, Sherlock & Rose Reference Turrell, Sherlock and Rose2015). For most applications, ignoring large-angle collisions is valid because scattering via small angles is approximately $8\ln {\it\Lambda}$ more likely than via large angles.

The evolution of the plasma is described by the moments of the transport equation with the transport of energy given by the second moment,

The definition of these quantities is as follows, with $\boldsymbol{v}_{rel}=\boldsymbol{v}-\boldsymbol{V}_{s}$ and $\unicode[STIX]{x1D644}$ the identity matrix (also denoted by ${\it\delta}_{{\it\mu}{\it\nu}}$ ), and $\boldsymbol{v},\boldsymbol{x},t$ dependence implicit:

Note that $C(f_{s},f_{s^{\prime }})$ appears in the collisional energy exchange terms $\boldsymbol{R}_{s}$ and $Q_{s}$ .

PIC simulations explicitly keep track of particles’ positions and velocities so that within each grid cell a distribution function taking the form of a sum over delta functions may be defined for each species. Let this be

for $\boldsymbol{v}_{i}$ the velocity of the $i$ th particle of species $s$ , $N_{s}$ the total number of computational particles in the space and time element under consideration, and ${\it\delta}(\boldsymbol{v})={\it\delta}(v_{x}){\it\delta}(v_{y}){\it\delta}(v_{z})$ .

## 4. Results

The energy density exchange rate is given by

and the momentum density exchange rate by

where $\boldsymbol{u}_{ij}=\boldsymbol{v}_{i}-\boldsymbol{v}_{j}$ .

## 5. Derivation

The $i$ and $j$ labels will frequently be implicit in the derivation of the results. For Greek indices, the Einstein summation convention of summation over repeated indices is followed. It is assumed that the distribution functions are as defined in § 3. The following relations are useful,

*a*,

*b*) $$\begin{eqnarray}\partial _{{\it\mu}}U_{{\it\mu}{\it\nu}}=-2u_{{\it\nu}}/u^{3},\quad \partial _{{\it\mu}}^{\prime }U_{{\it\mu}{\it\nu}}=2u_{{\it\nu}}/u^{3},\end{eqnarray}$$

*a*,

*b*) $$\begin{eqnarray}\partial _{{\it\mu}}\left(\frac{2u_{{\it\mu}}}{u^{3}}\right)=0,\quad \partial _{{\it\mu}}^{\prime }\left(\frac{2u_{{\it\mu}}}{u^{3}}\right)=0.\end{eqnarray}$$

The collision operator is

in which the second term, in $\partial _{{\it\mu}}f_{s}$ , may be trivially integrated to give

Note that there is an implicit $i$ label on the matrix $\unicode[STIX]{x1D650}$ . The first term, in $\partial _{{\it\mu}}^{\prime }f_{s^{\prime }}$ , can be re-expressed as

where $\hat{n}_{{\it\mu}}$ is the unit normal to the boundary, ${\it\Gamma}$ , of the velocity-space element $\text{d}^{3}v$ . The boundary term vanishes to give

where the dependence of the relative velocity on $j$ is again implicit. Then

The energy exchange rate density using this form of the collision operator is

which can be simplified to

where $V_{{\it\nu}}$ is the ${\it\nu}$ th component of $\boldsymbol{V}_{s}$ . Let the term in $U_{{\it\mu}{\it\nu}}$ be denoted by $(\ast )$ so that

and so

The full expression is

Evaluating the integral gives the result shown in § 4 for $Q_{s}$ .

The momentum density rate is given by

The ${\it\nu}$ th component of $\boldsymbol{R}_{s}$ is then

The second term, in $U_{{\it\mu}{\it\nu}}$ , yields a contribution in the ${\it\nu}$ th component of

which leads to the result stated in § 4.

Note that the total energy transfer rate per unit volume due to collisions is given by

where the $s^{\prime }$ label is implicit, and that

## 6. Benchmark

As a benchmark that the proposed equations give the correct values of the exchange rate densities, we evaluate the energy exchange density in a case for which there is a known analytical solution. The scenario is temperature equilibration between two species of ion with Maxwell–Boltzmann distributions,

The particle data used to evaluate $Q_{s}$ according to (4.1) is used to infer the rate of change of temperature of species $s$ from

where the above equation expresses temperature in units of degrees rather than in units of energy. The data used in (4.1) is taken from a particle Monte Carlo code which has been benchmarked for thermal equilibration problems (Turrell *et al.*
Reference Turrell, Sherlock and Rose2015), and which uses Takizuka and Abe’s collision algorithm (Takizuka & Abe Reference Takizuka and Abe1977). Particles populate a single point but have three dimensions in velocity space, a situation sometimes abbreviated as ‘0D3V’, representing a single cell of a PIC simulation. Both species’ begin the simulation in equilibrium with themselves, but not with each other. This is compared with the Landau–Spitzer temperature equilibration rate for two Maxwell–Boltzmann distributions, which is given by (Spitzer Reference Spitzer1967)

in which $s\neq s^{\prime }$ and $T_{s}(t=0)\neq T_{s^{\prime }}(t=0)$ . The derivation of this equation is based on the Landau collision operator (Trubnikov Reference Trubnikov1965) in its equivalent formulation in terms of Rosenbluth potentials (Rosenbluth, MacDonald & Judd Reference Rosenbluth, MacDonald and Judd1957). For two ion species $s$ and $s^{\prime }$ at a similar temperature and density, the temperature equilibration times are in the ratio

so that the intra-species equilibration time is shorter than the inter-species equilibration time regardless of the masses of the two species (Trubnikov Reference Trubnikov1965). Therefore, these species should maintain their own Maxwell–Boltzmann distributions during equilibration with each other.

The results are shown in figure 1 and are expressed as the ratio of the Landau–Spitzer equilibration rate to the rate predicted by (6.2) as a function of the initial temperature ratio between the two ion species. There is excellent agreement between the theory and the computational estimate (which is also dependent on the resolution of the simulation) across a large range of temperature ratios.

Figure 1 provides a baseline scenario which uses $N=N_{s}=N_{s^{\prime }}=3333$ as the number of computational particles of each type of ion, and a computational time step of ${\rm\Delta}t={\it\tau}_{0}/10$ . ${\it\tau}_{0}$ is ${\it\tau}_{ss^{\prime }}(t=0)=1/{\it\nu}_{0}$ . To show the effects of varying $N$ or ${\rm\Delta}t$ on the quality of the match with Landau–Spitzer theory, they are both varied relative to this base case. Figure 2 shows that the match is poorer with fewer computational particles as represented by the specific case of using $N/10$ . This is not a direct consequence of the equations shown in § 4 but a more general consequence of the reduction in accuracy when using fewer computational particles because, for example, the distribution functions of the ions will be more ‘noisy’, and less representative of the Maxwell–Boltzmann distribution assumed by the Landau–Spitzer theory. However, using $10N$ computational particles demonstrates that the equations do converge on the Landau–Spitzer theory as the computational resolution is increased, and that the equations presented are accurate.

Similarly, figure 3 shows the effect of varying the computational time step, ${\rm\Delta}t$ . The results are less sensitive to this parameter as neither the Landau–Spitzer rate nor the equations in § 4 explicitly depend on the computational time step. However, there is an implicit dependence of the simulations: the (Takizuka & Abe Reference Takizuka and Abe1977) scattering routine used in the simulations relies on repeated small-angle collisions between computational particles with a scattering angle ${\it\theta}$ chosen from a distribution whose variance scales as ${\rm\Delta}t$ . At large ${\rm\Delta}t$ , the scattering angles are no longer guaranteed to be small and the theory breaks down. Using ${\rm\Delta}t\gg {\it\tau}/10$ would cause this to happen, and indirectly cause the results of evaluating (6.2) to become inaccurate. Figure 3 demonstrates that for ${\rm\Delta}t={\it\tau}/10$ , varying ${\rm\Delta}t$ by an order of magnitude largely maintains the accuracy of the equations in § 4.

## 7. Conclusion

We have derived equations for plasma collisional energy exchange terms which are less computationally intensive to evaluate than the alternative methods available which use the output of PIC simulations. They involve a direct sum of the properties of particles rather than the integration of distribution functions over a six-dimensional velocity phase space. Results obtained using this method and data taken from particle simulations are in excellent agreement with the known test case problem of thermal equilibration, demonstrating that the derived equations are accurate.

## Acknowledgements

A.E.T. was supported by EPSRC grants EP/L504786/1 and EP/K028464/1 while undertaking this research.