2.2.1. Regularisations

At the core of our scheme is an absolutely continuous regularisation of $f^N$. We will consider a spatial spline $\psi _{\eta }$, with scale parameters $\eta =(\eta _1,\ldots,\eta _{d_x})$. The spline must satisfy the following properties:

1. (i) non-negativity: $\psi _{\eta }(x)\geq 0$;

2. (ii) symmetry: $\psi _{\eta }(-x) = \psi _{\eta }(x)$;

3. (iii) unit mean: $\int _{{\mathbb {R}}^{d_x}} \psi _{\eta }(x) \,\mathrm{d}\kern0.07em x = 1$.

We will also consider a velocity spline $\varphi _{\varepsilon }$, with scale parameters $\varepsilon =(\varepsilon _1,\ldots,\varepsilon _{d_v})$, which must satisfy the same properties.

For the sake of concreteness, we choose the shape function

(2.14)\begin{equation} {G}(x) = \begin{cases} 1 - |x|, & \text{if } |x| \leq 1, \\ 0, & \text{otherwise}, \end{cases} \end{equation}

and define the splines

(2.15a,b)\begin{gather} \psi_{\eta}(x) := \frac{1}{\eta^{d_x}} \hat{\psi}_{\eta}(x), \quad \hat{\psi}_{\eta}(x) := \prod_{d=1}^{{d_x}} {G}\left(\frac{x_d}{\eta_d}\right), \end{gather}

(2.16a,b)\begin{gather}\varphi_{\varepsilon}(x) := \frac{1}{\varepsilon^{d_v}} \hat{\varphi}_{\varepsilon}(v), \quad \hat{\varphi}_{\varepsilon}(v) := \prod_{d=1}^{{d_v}} {G}\left(\frac{v_d}{\varepsilon_d}\right). \end{gather}

Here, $\eta ^{d_x}$ stands for the product $\prod _{d=1}^{{d_x}}\eta _d$, and similarly $\varepsilon ^{d_v}$ stands for $\prod _{d=1}^{{d_v}}\varepsilon _d$. We also use the convention $x=(x_1,\ldots,x_{d_x})$ and $v=(v_1,\ldots,v_{d_v})$. This choice of splines is practical, not only because of its simplicity, but also because the compact support of the splines can be exploited for computational efficiency (see § 2.3.1). However, we remark that the properties discussed in §§ 2.4 and 2.5 do not rely on this choice, only on the three properties discussed above.

Given a suitable choice of splines, we can define the regularised solution

(2.17)\begin{equation} \tilde{f}^N(t,x,v) := (\psi_{\eta}\varphi_{\varepsilon} \ast_{x,v} f^N)(t,x,v) = \sum_{p=1}^{N} w_{p} \psi_{\eta}(x - x_{p}(t)) \varphi_{\varepsilon}(v - v_{p}(t)), \end{equation}

where ‘$\ast _{x,v}$’ denotes the double convolution in $x$ and $v$.

2.2.2. Lorentz term

In order to determine the electromagnetic force experienced by each particle, we must define the terms $E(t,x_{p})$ and $B(t,x_{p})$. To do so, information from the particles is first interpolated to a spatial mesh, where Maxwell's equations are solved. That solution is then extrapolated back to the particles.

The first step is to define a regularised electric charge current, specified as an integral of $\tilde {f}^N$:

(2.18)\begin{equation} \tilde J(t,x) := \int_{{\mathbb{R}}^{d_v}} v \tilde{f}^N(t,x,v) \,\mathrm{d} v = \sum_{p=1}^{N} w_{p} v_{p}(t) \psi_{\eta}(x - x_{p}(t)). \end{equation}

This can be evaluated on a spatial mesh $\varOmega _{h} = \{ x_{h} \}_{h}$ (assumed here to be a tensorised grid with mesh size $h_d$ on each dimension), yielding $\tilde J(t,x_{h})$. These mesh values will be used as inputs to solve Maxwell's equations, yielding the electromagnetic fields $\tilde E(t,x_{h})$ and $\tilde B(t,x_{h})$ at the mesh points.

In practice, we only solve Ampère's and Faraday's equations (2.1b); Poisson's and Gauss’ equations (2.1c) are only solved once, in order to determine $\tilde E(0,x_{h})$ and $\tilde B(0,x_{h})$ in a way that is consistent with the datum $\tilde {f}^N(0,x,v)$. For further particulars, see § 2.6.

Once the values of the fields are known on the mesh, they are extrapolated back to the particle positions,

(2.19)\begin{gather} E(t,x_{p}) := \sum_h \tilde E(t,x_{h}) \hat{\psi}_{\eta}(x_{p} - x_{h}) = \sum_h \tilde E(t,x_{h}) \psi_{\eta}(x_{p} - x_{h}) \eta^{d_x}, \end{gather}

(2.20)\begin{gather}B(t,x_{p}) := \sum_h \tilde B(t,x_{h}) \hat{\psi}_{\eta}(x_{p} - x_{h}) = \sum_h \tilde B(t,x_{h}) \psi_{\eta}(x_{p} - x_{h}) \eta^{d_x}. \end{gather}

This defines the Lorenz force acting on the $p\textrm {th}$ particle, $E(t,x_{p}) + v_{p} \times B(t,x_{p})$. If the terms (2.19) and (2.20) are seen as a spline reconstruction over a mesh,

(2.21a,b)\begin{equation} E(t,x_{p}) = \sum_h \tilde E(t,x_{h}) \hat{\psi}_h(x_{p} - x_{h}), \quad B(t,x_{p}) = \sum_h \tilde B(t,x_{h}) \hat{\psi}_h(x_{p} - x_{h}), \end{equation}

it becomes evident that the choice of mesh size $h_d=\eta _d$ is, in fact, the only reasonable one, for any other would not scale appropriately.

2.2.3. Collision term

We now construct a regularised collisional velocity field $\mathcal {U}_{\eta,\varepsilon }[f](x,v)$. By analogy with (2.1), we define, for an arbitrary distribution $f$,

(2.22a)\begin{gather} \mathcal{U}_{\eta,\varepsilon}[f](x,v) = \iint_{\varOmega\times{\mathbb{R}}^{d_v}} \psi_{\eta} (x - x_*) A (v - v_*) b _{\eta,\varepsilon} [f] (x,x_*,v,v_*) f (x_*,v_*) \,\mathrm{d} v_* \,\mathrm{d}\kern0.07em x_* , \end{gather}

(2.22b)\begin{gather}b _{\eta,\varepsilon} [f] (x,x_*,v,v_*) = \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f] (x,v) - \boldsymbol{\nabla}_{v_*} \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f] (x_*,v_*), \end{gather}

(2.22c)\begin{gather}{\mathcal{H}_{\eta,\varepsilon} [f]} = { -\mathcal{S}_{\eta,\varepsilon} [f] ,\quad \mathcal{S}_{\eta,\varepsilon} [f] ={-}\iint_{\varOmega\times{\mathbb{R}}^{d_v}} f \log(\psi_{\eta}\varphi_{\varepsilon} \ast_{x,v} f) \,\mathrm{d} v\,\mathrm{d}\kern0.07em x ,} \end{gather}

where $A$ is the collision kernel (2.1i).

The departure from (2.1) is twofold: first, system (2.22) is well-defined for discrete measures; second, the entropy and collision operators have been altered to include spatial dependencies. Note that the discrete collision operator has been defined in terms of the total entropy $\mathcal {S}$, rather than the function of space $H(x)$ that appears in the continuous Landau operator. The physical interpretation of (2.22) is that first we delocalise in space, to consider collisions in a neighbourhood of each point $x$, and then we localise the interactions in order to associate each particle with a distinct position. This regularised collision somewhat resembles the Enskog collision operator, which has been studied in the context of the Boltzmann equation (Villani Reference Villani2006).

To define the collisional component of the scheme, we will evaluate $\mathcal {U}_{\eta,\varepsilon }$ for $f^N$ at the particle coordinates $(x_{p},v_{p})$. First, we note that the regularised functional (2.22c) can be rewritten as

(2.23)\begin{equation} \mathcal{H}_{\eta,\varepsilon} [f] = \iint_{\varOmega\times{\mathbb{R}}^{d_v}} f \log(\tilde f)\,\mathrm{d} v\,\mathrm{d}\kern0.07em x, \end{equation}

using $\tilde f$ as shorthand for $(\psi _{\eta }\varphi _{\varepsilon } \ast _{x,v} f)$. It is straightforward to compute

(2.24)\begin{equation} \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f] (x,v) = \log(\tilde f) + \left( \frac{f}{\,\tilde f\,} \ast_{x,v} (\psi_{\eta}\varphi_{\varepsilon}) \right) \end{equation}

and

(2.25)\begin{equation} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f] (x,v) = \frac{\boldsymbol{\nabla}_v \tilde f}{\tilde f} + \left( \frac{f}{\,\tilde f\,} \ast_{x,v} (\psi_{\eta}\boldsymbol{\nabla}_v\varphi_{\varepsilon}) \right). \end{equation}

We now structure the evaluation of the collision term in three sequential steps.

1. (Step I) Compute the values of $\tilde {f}^N$ and $\boldsymbol {\nabla }_v\tilde {f}^N$ at each particle,

   (2.26a)\begin{gather} \tilde{f}^N(x_{p},v_{p}) = \sum_{q=1}^{N} w_{q} \psi_{\eta} (x_{p} - x_{q}) \varphi_{\varepsilon} (v_{p} - v_{q}), \end{gather}
   (2.26b)\begin{gather}\boldsymbol{\nabla}_v \tilde{f}^N(x_{p},v_{p}) = \sum_{q=1}^{N} w_{q} \psi_{\eta} (x_{p} - x_{q}) \boldsymbol{\nabla}_v \varphi_{\varepsilon} (v_{p} - v_{q}). \end{gather}
   (Step II) Compute the values of $\boldsymbol {\nabla }_v ({\partial \mathcal {H}_{\eta,\varepsilon }}/{\partial f}) [f^N]$ at each particle,
   (2.26c)\begin{align} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) & = \frac{\boldsymbol{\nabla}_v \tilde{f}^N(x_{p},v_{p})}{\tilde{f}^N(x_{p},v_{p})} + \left( \frac{f^N}{\tilde{f}^N} \ast_{x,v} (\psi_{\eta}\boldsymbol{\nabla}_v\varphi_{\varepsilon}) \right) (x_{p},v_{p}) \nonumber\\ & = \frac{\boldsymbol{\nabla}_v \tilde{f}^N(x_{p},v_{p})}{\tilde{f}^N(x_{p},v_{p})} + \sum_{q=1}^{N} w_{q} \left( \frac{ \psi_{\eta} (x_{p} - x_{q}) \boldsymbol{\nabla}_v \varphi_{\varepsilon} (v_{p} - v_{q}) }{ \tilde{f}^N(x_{q},v_{q}) } \right) . \end{align}
   (Step III) Compute the velocity field,
   (2.26d)\begin{equation} \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p}) =\sum_{q=1}^{N} w_{q} \psi_{\eta} (x_{p} - x_{q}) A_{p,q} b_{p,q}, \end{equation}
   where
   (2.26e)\begin{equation} A_{p,q} = A (v_{p} - v_{q}) , \quad b_{p,q} = \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) - \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{q},v_{q}) . \end{equation}

Section 2.3 discusses the practical implementation of these steps: the use of a cell list and a random batch approach can greatly reduce the computational cost of the method.

Remark 2.2 An alternative regularisation of the functional is

(2.27)\begin{equation} \iint_{\varOmega\times{\mathbb{R}}^{d_v}} (\psi_{\eta}\varphi_{\varepsilon} \ast_{x,v} f) \log(\psi_{\eta}\varphi_{\varepsilon} \ast_{x,v} f) \,\mathrm{d} v\,\mathrm{d}\kern0.07em x, \end{equation}

reminiscent of the one employed in Carrillo et al. (Reference Carrillo, Hu, Wang and Wu2020) for the homogeneous Landau equation. Under this regularisation, the scheme will achieve very similar properties (see §§ 2.4 and 2.5). However, the calculation of the corresponding $\boldsymbol {\nabla }_v ({\delta \mathcal {H}_{\eta,\varepsilon }}/{\delta f}) [f^N] (x_{p},v_{p})$ term will require the evaluation of an integral. The need for numerical quadrature makes this approach less convenient.

2.2.4. Choice of regularisation parameters

The method involves two parameters, $\eta$ and $\varepsilon$. However, it is neither obvious how their value should be chosen, nor how it relates to the number of particles $N$. We present here a convention to make this choice more intuitive.

The computation of the Lorentz terms involves a spatial mesh, see §§ 2.2.2 and 2.6. Without loss of generality, we suppose the spatial domain is of the form $\varOmega = (0,L_1) \times \cdots \times (0,L_{d_x})$ for $L_1,\ldots,L_{d_x} > 0$. We assume that a uniform mesh is employed, and that each dimension is discretised with $N_{x_d}$ cells. We then choose

(2.28)\begin{equation} \eta_d := \frac{L_d}{N_{x_d}} \end{equation}

for each dimension $d$. The size of the spatial spline can therefore be understood as consequence of the spatial resolution of the scheme.

We apply similar thinking for the velocity splines. The method does not require a mesh in velocity; however, we imagine a fictitious mesh, and inspire a similar parameter choice. We stress, however, that this mesh is never used in computation, it is simply a thinking tool. Armed with this caveat, we approximate the velocity domain as $(-L_v,L_v)^{d_v}$, for some $L_v>0$. For each dimension, we choose a number of fictitious cells $N_{v_d}$, and let

(2.29)\begin{equation} \varepsilon_d := \frac{L_v}{N_{v_d}} \end{equation}

for each dimension $d$.

Given now the full mesh (space and velocity combined), we populate it with a number of particles per cell $N_c$. The total number of particles shall therefore be

(2.30)\begin{equation} N := N_{x_1} \times \cdots \times N_{x_{d_x}} \times N_{v_1} \times \cdots \times N_{v_{d_v}} \times N_c, \end{equation}

where $N_c$ is a measure of the granular resolution of the scheme. Note that, while we speak of ‘particles per cell’, the particles are not confined to individual cells and may roam the numerical domain. Therefore, the number of particles within a specific cell in phase space is not fixed a priori.

To think of a notion of convergence of the scheme, we fix $N_c$ (typically $N_c=8$) and let $N_{x_d}$ and $N_{v_d}$ grow (see § 3.1).

2.3.1. Cell list optimisation

Cell lists (Allen & Tildesley Reference Allen and Tildesley1990, Chapter 5.3.2) are data structures commonly used to efficiently find all pairs of particles within a given distance of each other. Their use in particle simulations can reduce the complexity of computing short-range particle interactions from $O(N^2)$ to $O(N \log N)$.

First, each particle is located within an auxiliary mesh (an $O(N)$ operation), and a list of which particle belongs in each mesh cell is stored. The size of the cells is chosen in relation to the size of the support of the splines, so that terms of the form $\psi _{\eta } (x_{p} - x_{q})$ can only be non-zero if both particles lie in the same cell, or in cells which are neighbouring in space. Since the lists of particles in each cell are precomputed, looking up each particle's ‘neighbours’ is trivial. The same logic applies in the velocity dimensions, with the term $\varphi _{\varepsilon } (v_{p} - v_{q})$.

The complexity of the three steps of § 2.2.3 reduces as follows.

1. (Step I) This step involves $N$ particles with a compactly supported interaction; therefore, the evaluation cost is reduced from $O(N^2)$ to $O(N \log N)$.

2. (Step II) This step is as the previous one.

3. (Step III) This step involves an interaction which is local in space only (the Landau operator is, after all, non-local in velocity). Therefore, the cost reduction is lesser. The typical number of particles within one spatial cell is $O(N_{v_1} \times \cdots \times N_{v_{d_v}} \times N_c)$; the evaluation of the collision operator will be of quadratic order on this quantity. Accounting for the spatial cells, which do benefit from the localisation, the overall cost is estimated as

   (2.31)\begin{equation} O( N_{x_1} \times \cdots \times N_{x_{d_x}} \times \log(N_{x_1} \times \cdots \times N_{x_{d_x}}) \times (N_{v_1} \times \cdots \times N_{v_{d_v}} \times N_c)^2 ), \end{equation}
   which will be the dominant complexity of the numerical scheme.

We remark that the cell list optimisation does not alter the method, as it only discards contributions which are exactly zero.

2.3.2. Random batch optimisation

The random batch method (Jin, Li & Liu Reference Jin, Li and Liu2020; Carrillo et al. Reference Carrillo, Jin and Tang2022b) is also a staple of particle simulations. Given the $N$ particles, we assign them randomly to $R$ batches of equal size $N/R$. Steps I and II of § 2.2.3 are computed as before, but the collision operator of step III is replaced by

(2.32)\begin{equation} \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p}) =\frac{R(N-1)}{N-R}\sum_{q\in {\mathcal{B}}_p} w_{q} \psi_{\eta} (x_{p} - x_{q}) A_{p,q} b_{p,q}, \end{equation}

where ${\mathcal {B}}_p$ is the batch which contains the $p\textrm {th}$ particle. This term acts as an unbiased estimator of the original velocity field.

The complexity of the third step of § 2.2.3 is reduced by the batching. Since the batch size is $N/R$, the interaction is quadratic, and it must be computed for each batch, the overall cost is $O(N^2\times R^{-1})$. If we combine this technique with the cell list, the resultant complexity is estimated as

(2.33)\begin{equation} O( N_{x_1} \times \cdots \times N_{x_{d_x}} \times \log(N_{x_1} \times \cdots \times N_{x_{d_x}}) \times (N_{v_1} \times \cdots \times N_{v_{d_v}} \times N_c)^2 \times R^{{-}1} ). \end{equation}

This remains the dominant complexity of the numerical scheme.

The random batch optimisation does alter the numerical method. However, it is easy to prove that the structural properties of the method presented in §§ 2.4 and 2.5 survive the batching; see Carrillo et al. (Reference Carrillo, Jin and Tang2022b) for details, where the random batch method was used to accelerate the particle method for the homogeneous Landau equation (Carrillo et al. Reference Carrillo, Hu, Wang and Wu2020). Furthermore, the use of batches does not significantly affect the accuracy of the method, see § 3.1.

2.4.1. Collision invariants

The conservation properties of the method arise naturally from the definition of the discrete collision operator (2.26). Specifically, we are able to emulate the weak form (2.2) of the Landau operator at the discrete level, which shows that the scheme preserves the usual collision invariants: $1$, $v$ and $\tfrac {1}{2}|v|^2$.

The discrete weak form is easily found. By analogy with

(2.34)\begin{equation} \int_{{\mathbb{R}}^{d_v}} g(v) Q[f,f] \,\mathrm{d} v = \int_{{\mathbb{R}}^{d_v}} g(v) \boldsymbol{\nabla}_v\boldsymbol{\cdot}(f\mathcal{U}[f]) \,\mathrm{d} v ={-}\int_{{\mathbb{R}}^{d_v}} \boldsymbol{\nabla}_v g(v) \boldsymbol{\cdot} (f\mathcal{U}[f]) \,\mathrm{d} v, \end{equation}

we compute, for a test function $g(x,v)$,

(2.35)\begin{align} & -\sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v g(x_{p},v_{p}) \boldsymbol{\cdot} \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p})\nonumber\\ & \quad={-}\sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v g(x_{p},v_{p}) \boldsymbol{\cdot} \sum_{q=1}^{N} w_{q} \psi_{\eta} (x_{p} - x_{q}) A_{p,q} b_{p,q}\nonumber\\ & \quad={-}\sum_{p,q=1}^{N} w_{p} w_{q} \psi_{\eta} (x_{p} - x_{q}) \boldsymbol{\nabla}_v g(x_{p},v_{p}) \boldsymbol{\cdot} A_{p,q} b_{p,q}\nonumber\\ & \quad={-}\frac{1}{2}\sum_{p,q=1}^{N} w_{p} w_{q} \psi_{\eta} (x_{p} - x_{q}) [ \boldsymbol{\nabla}_v g(x_{p},v_{p}) - \boldsymbol{\nabla}_v g(x_{q},v_{q})] \boldsymbol{\cdot} A_{p,q} b_{p,q}, \end{align}

symmetrising on the last line, and exploiting the symmetry $\psi _{\eta } (x_{p} - x_{q})$ and the antisymmetry of $b_{p,q}$. This term vanishes trivially for $g=1$ as well as $g=v$. To see that it also vanishes for $g=\tfrac {1}{2}|v|^2$, one uses the fact that $A_{p,q}$ projects $b_{p,q}$ onto the perpendicular to $v_{p} - v_{q}$. We therefore conclude as follows.

Theorem 2.3 (Discrete collision invariants)

The functions $g=1$, $g=v$ and $g=\tfrac {1}{2}|v|^2$ are collision invariants of the C-PIC collision operator:

(2.36)\begin{equation} -\sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v g(x_{p},v_{p}) \boldsymbol{\cdot} \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p}) = 0. \end{equation}

2.4.2. The H-Theorem

The discrete H-theorem in phase space also holds, and is readily obtained from the discrete weak form. At the continuous level, one would choose the test function $g=\log f+1={\partial \mathcal {H}}/{\partial f}$. Instead, we choose $g=({\partial \mathcal {H}_{\eta,\varepsilon }}/{\partial f}) [f^N]$, to find

(2.37)\begin{align} \mathcal{D}_{\eta,\varepsilon} & := \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) \boldsymbol{\cdot} \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p})\nonumber\\ & = \frac{1}{2}\sum_{p,q=1}^{N} w_{p} w_{q} \psi_{\eta} (x_{p} - x_{q}) [ \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) - \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{q},v_{q})] \boldsymbol{\cdot} A_{p,q} b_{p,q}\nonumber\\ & = \frac{1}{2}\sum_{p,q=1}^{N} w_{p} w_{q} \psi_{\eta} (x_{p} - x_{q}) b_{p,q} \boldsymbol{\cdot} A_{p,q} b_{p,q}, \end{align}

which is non-positive due to the positive semidefinite property of the matrix $A_{p,q}$. We therefore conclude as follows.

Theorem 2.4 (Discrete H-theorem)

The H-theorem holds for the C-PIC entropy and variation:

(2.38)\begin{equation} \mathcal{D}_{\eta,\varepsilon} = \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) \boldsymbol{\cdot} \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p}) \geq 0. \end{equation}

2.5.1. Conservation of charge and energy

We may compute the global evolution in time of a test function $g(x,v)$:

(2.39)\begin{align} & \frac{\mathrm{d} }{\mathrm{d} t} \iint_{\varOmega\times{\mathbb{R}}^{d_v}} g(x,v) f^N(t,x,v) \,\mathrm{d} v\,\mathrm{d}\kern0.07em x = \frac{\mathrm{d} }{\mathrm{d} t} \sum_{p=1}^{N} w_{p} g(x_{p},v_{p})\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} [ \boldsymbol{\nabla}_x g(x_{p},v_{p}) \boldsymbol{\cdot} \dot{x}_{p} + \boldsymbol{\nabla}_v g(x_{p},v_{p}) \boldsymbol{\cdot} \dot{v}_{p} ]\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} [ \boldsymbol{\nabla}_x g(x_{p},v_{p}) \boldsymbol{\cdot} v_{p} + \boldsymbol{\nabla}_v g(x_{p},v_{p}) \boldsymbol{\cdot} ( E(t,x_{p}) + v_{p} \times B(t,x_{p}) - \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p}) )]. \end{align}

Mass and charge ($g=1$) are trivially conserved at the global level.

The evolution of the kinetic energy ($g=\tfrac {1}{2}|v|^2$) is given by

(2.40)\begin{align} & \frac{\mathrm{d} }{\mathrm{d} t} \frac{1}{2} \iint_{\varOmega\times{\mathbb{R}}^{d_v}}|v|^2 f^N(t,x,v) \,\mathrm{d} v\,\mathrm{d}\kern0.07em x = \frac{\mathrm{d} }{\mathrm{d} t} \frac{1}{2} \sum_{p=1}^{N} w_{p} |v_{p}|^2\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} v_{p} \boldsymbol{\cdot} ( E(t,x_{p}) + v_{p} \times B(t,x_{p}) - \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p}) )\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} v_{p} \boldsymbol{\cdot} E(t,x_{p})\nonumber\\ & \quad = \sum_h \sum_{p=1}^{N} w_{p} v_{p} \boldsymbol{\cdot} \tilde E(t,x_{h}) \psi_{\eta}(x_{p} - x_{h}) h^{d_x}\nonumber\\ & \quad = \sum_h \tilde J (t,x_{h}) \boldsymbol{\cdot} \tilde E(t,x_{h}) h^{d_x}, \end{align}

where we have used theorem 2.3, the conservation of momentum at the level of the collision operator. At the continuous level, the conservation of energy identity would be found by invoking Maxwell's equations in order to relate $J\cdot E$ to the electromagnetic energy. Namely

(2.41)\begin{equation} \partial_t \left( \frac{1}{2}|E|^2 + \frac{1}{2}|B|^2 \right) = E \cdot \partial_t E + B \cdot \partial_t B ={-}J \cdot E - \boldsymbol{\nabla}_x\boldsymbol{\cdot} (E \times B), \end{equation}

a computation that requires the vector identity

(2.42)\begin{equation} \boldsymbol{\nabla}_x\boldsymbol{\cdot} (E\times B)=(\boldsymbol{\nabla}_x\times E)\boldsymbol{\cdot} B-E\boldsymbol{\cdot} (\boldsymbol{\nabla}_x \times B). \end{equation}

Remark 2.5 (Conservation of energy in the fully discrete setting)

Both the chain rule in time and the vector identity in space will not hold for a general discretisation; however, the correct approach is well known. In space, Yee's lattice (Yee Reference Yee1966; Hyman & Shashkov Reference Hyman and Shashkov1999) leads to the corresponding discrete vector identities. In time, implicit schemes (the midpoint rule, in particular) can be used to preserve the chain rule (Chen et al. Reference Chen, Chacón and Barnes2011; Markidis & Lapenta Reference Markidis and Lapenta2011; Crouseilles, Einkemmer & Faou Reference Crouseilles, Einkemmer and Faou2015; Lapenta Reference Lapenta2017; Kormann & Sonnendrücker Reference Kormann and Sonnendrücker2021). Under such conditions, we would arrive at the balance of energy:

(2.43)\begin{align} & \frac{\mathrm{d} }{\mathrm{d} t} ( \mathcal{E}_K +\mathcal{E}_E +\mathcal{E}_B ) \nonumber\\ & \quad := \frac{\mathrm{d} }{\mathrm{d} t} \left( \frac{1}{2} \sum_{p=1}^{N} w_{p} |v_{p}|^2 + \frac{1}{2}\sum_h |\tilde E(t,x_{h})|^2 h^{d_x} + \frac{1}{2}\sum_h |\tilde B(t,x_{h})|^2 h^{d_x}\right) \equiv 0. \end{align}

The numerical examples of § 3 are discretised using Yee's lattice (see § 2.6), but employ an explicit time integrator for the sake of efficiency. Therefore, energy is only conserved at the semidiscrete level (continuous in time), but generally not conserved at the fully discrete level.

Remark 2.7 (Local conservation of charge)

The local conservation of charge equation,

(2.44)\begin{equation} \partial_t\rho + \boldsymbol{\nabla}\boldsymbol{\cdot} J = 0, \end{equation}

a continuity equation for the electric charge, can be derived from Maxwell's equations and also from Vlasov's equation. This equation is intimately related to the conservation of energy; as such, it is desirable that a PIC scheme should verify a discrete version of (2.44). PIC methods which exhibit local charge conservation have been proposed in several works, including Villasenor & Buneman (Reference Villasenor and Buneman1992), Esirkepov (Reference Esirkepov2001), Chen et al. (Reference Chen, Chacón and Barnes2011) and Chen et al. (Reference Chen, Chacón, Yin, Albright, Stark and Bird2020).

2.5.2. Increment of entropy

At the continuous level, the global increment of the entropy (2.6) is a consequence of two facts: the H-theorem, and the fact that the Vlasov equation conserves entropy. Since the discrete H-theorem has already been shown in § 2.4, we must study the conservation.

The evolution of the regularised entropy is not a priori given by the expression in § 2.5.1 since the test function, $g=({\partial \mathcal {H}_{\eta,\varepsilon }}/{\partial f}) [f^N]$, also depends in time. However, a direct computation shows

(2.45)\begin{align} & \frac{\mathrm{d} }{\mathrm{d} t} {\mathcal{S}_{\eta,\varepsilon}} [f^N] = {-} \frac{\mathrm{d} }{\mathrm{d} t} \iint_{\varOmega\times{\mathbb{R}}^{d_v}} f^N \log(\tilde{f}^N) \,\mathrm{d} v\,\mathrm{d}\kern0.07em x = {-} \frac{\mathrm{d} }{\mathrm{d} t} \sum_{p=1}^{N} w_{p} \log(\tilde{f}^N(x_{p},v_{p}))\nonumber\\ & \quad = {-} \sum_{p,q=1}^{N} w_{p} w_{q} \frac{ \frac{\mathrm{d} }{\mathrm{d} t} [ \vphantom{x^2} \psi_{\eta}(x_{p} - x_{q}) \varphi_{\varepsilon}(v_{p} - v_{q}) ] }{ \tilde{f}^N(x_{p},v_{p}) }\nonumber\\ & \quad = {-} \sum_{p,q=1}^{N} w_{p} w_{q} \frac{\begin{array}{@{}l@{}} \boldsymbol{\nabla}_x \psi_{\eta}(x_{p} - x_{q}) \varphi_{\varepsilon}(v_{p} - v_{q}) \boldsymbol{\cdot} (\dot x_{p} - \dot x_{q})\\ \quad+ \psi_{\eta}(x_{p} - x_{q}) \boldsymbol{\nabla}_v \varphi_{\varepsilon}(v_{p} - v_{q}) \boldsymbol{\cdot} (\dot v_{p} - \dot v_{q}) \end{array}}{ \tilde{f}^N(x_{p},v_{p}) } . \end{align}

The first summand can be greatly simplified, noting

(2.46)\begin{align} & \sum_{p,q=1}^{N} w_{p} w_{q} \frac{ \boldsymbol{\nabla}_x \psi_{\eta}(x_{p} - x_{q}) \varphi_{\varepsilon}(v_{p} - v_{q}) \boldsymbol{\cdot} (\dot x_{p} - \dot x_{q}) }{ \tilde{f}^N(x_{p},v_{p}) }\nonumber\\ & \quad = \sum_{p,q=1}^{N} w_{p} w_{q} \dfrac{ \boldsymbol{\nabla}_x \psi_{\eta}(x_{p} - x_{q}) \varphi_{\varepsilon}(v_{p} - v_{q}) }{ \tilde{f}^N(x_{p},v_{p}) } \boldsymbol{\cdot} \dot x_{p}\nonumber\\ & \qquad + \sum_{p,q=1}^{N} w_{p} w_{q} \frac{ \boldsymbol{\nabla}_x \psi_{\eta}(x_{p} - x_{q}) \varphi_{\varepsilon}(v_{p} - v_{q}) }{ \tilde{f}^N(x_{q},v_{q}) } \boldsymbol{\cdot} \dot x_{p}\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} \left[ \frac{ \boldsymbol{\nabla}_x \tilde{f}^N(x_{p},v_{p}) }{ \tilde{f}^N(x_{p},v_{p}) } + \sum_{q=1}^{N} w_{q} \frac{ \boldsymbol{\nabla}_x \psi_{\eta}(x_{p} - x_{q}) \varphi_{\varepsilon}(v_{p} - v_{q}) }{ \tilde{f}^N(x_{q},v_{q}) } \right] \boldsymbol{\cdot} \dot x_{p}\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} \left[ \frac{ \boldsymbol{\nabla}_x \tilde{f}^N(x_{p},v_{p}) }{ \tilde{f}^N(x_{p},v_{p}) } + \left( \frac{f^N}{\tilde{f}^N} \ast_{x,v} (\psi_{\eta}\boldsymbol{\nabla}_v\varphi_{\varepsilon})\right) (x_{p},v_{p}) \right] \boldsymbol{\cdot} \dot x_{p}\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_x \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) \boldsymbol{\cdot} \dot x_{p} = \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_x \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) \boldsymbol{\cdot} v_{p}, \end{align}

where we have swapped the labels $p$ and $q$ in the second term of the second line, and have used the antisymmetry of $\boldsymbol {\nabla }_x \psi _{\eta }(x_{p} - x_{q})$. The last expression is the discrete analogue of the term $\iint _{\varOmega \times {\mathbb {R}}^{d_v}} (\log f +1 ) \boldsymbol {\nabla }_x\boldsymbol {\cdot } (fv) \,\mathrm {d} v\,\mathrm{d}\kern0.07em x$, which vanishes under reasonable assumptions. However, the discrete term is not zero in general; i.e. the PIC method does not conserve the regularised entropy exactly. We define the entropy conservation error in position,

(2.47)\begin{equation} \mathcal{C}_x := {-} \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_x \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) \boldsymbol{\cdot} v_{p}, \end{equation}

which will be shown to be numerically negligible.

The second summand in the evolution of $\mathcal {H}_{\eta,\varepsilon } [f^N]$ similarly becomes

(2.48)\begin{align} & \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] (x_{p},v_{p}) \boldsymbol{\cdot} \dot v_{p}\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] \boldsymbol{\cdot} ( E(t,x_{p}) + v_{p} \times B(t,x_{p} ) - \mathcal{U}_{\eta,\varepsilon}[f^N](x_{p},v_{p}) )\nonumber\\ & \quad = \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] \boldsymbol{\cdot} ( E(t,x_{p}) + v_{p} \times B(t,x_{p}) ) - \mathcal{D}_{\eta,\varepsilon}\nonumber\\ & \quad \leq \sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] \boldsymbol{\cdot} ( E(t,x_{p}) + v_{p} \times B(t,x_{p}) ), \end{align}

using the discrete H-theorem (theorem 2.4). The remaining quantity is the discrete analogue of the continuous term $\iint _{\varOmega \times {\mathbb {R}}^{d_v}} (\log f +1 ) \boldsymbol {\nabla }_v\boldsymbol {\cdot } (f(E+v\times B)) \,\mathrm {d} v\,\mathrm{d}\kern0.07em x$; once again, the continuous integral vanishes, but the discrete term does not, in general. We define the entropy conservation error in velocity,

(2.49)\begin{equation} \mathcal{C}_v := {-}\sum_{p=1}^{N} w_{p} \boldsymbol{\nabla}_v \frac{\delta \mathcal{H}_{\eta,\varepsilon}}{\delta f} [f^N] \boldsymbol{\cdot} ( E(t,x_{p}) + v_{p} \times B(t,x_{p}) ), \end{equation}

which will also be shown to be numerically negligible.

To summarise, we find

(2.50)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} {\mathcal{S}_{\eta,\varepsilon}} [f^N] = \mathcal{C}_x + \mathcal{C}_v {+\mathcal{D}_{\eta,\varepsilon}}, \end{equation}

where ${\mathcal {D}_{\eta,\varepsilon }\geq 0}$. Because the PIC method does not exactly conserve entropy, the C-PIC method cannot be said to exactly increase it. However, both $\mathcal {C}_x$ and $\mathcal {C}_v$ will be shown to be numerically negligible in § 3.1.2.

2.6.1. Field update

In 1D-2V, Maxwell's equations reduce to

(2.51)\begin{equation} \partial_t E_1 ={-} J_1, \quad \partial_t E_2 ={-} J_2 - \partial_x B_3, \quad \partial_t B_3 ={-} \partial_x E_2, \quad x \in \varOmega = (0,L), \end{equation}

for some $L>0$, using the convention $E=(E_1,E_2,0)$, $B=(0,0,B_3)$, $J=(J_1,J_2,0)$. We prescribe periodic boundary conditions in space.

In order to discretise (2.51), we choose a number of cells $N_x$, and compute the corresponding regularisation parameter $\eta =L/N_x$. We construct two meshes with spacing $\eta$; the spacing of the mesh must be related to the spatial regularisation parameter, see §§ 2.2.2 and 2.2.4. We define a primal mesh $\varOmega =\{ x_{i} \mid 1\leq i\leq N_x \}$ and a dual mesh $\varOmega '=\{ x_{\textrm {i}+1/2} \mid 0\leq i\leq N_x \}$, where $x_{i} := (i-1/2) \eta$.

The field update begins with the particles and the electromagnetic field at time $t^{n}$. The positions and velocities of the particles are $x_{p}^{n} = (x_{1,p}^{n})$ and $v_{p}^{n} = (v_{1,p}^{n}, v_{2,p}^{n})$ for $1\leq p\leq N$. The electric field is considered on the primal mesh; the first component is $\tilde {E}_{1,i}^{n}$, and the second component is $\tilde {E}_{2,i}^{n}$, for $1\leq i\leq N_x$. The magnetic field is considered on the dual mesh; the third component is $\tilde {B}_{3,i}^{n}$ for $0\leq i\leq N_x$.

The current is evaluated over the primal mesh $\varOmega$ through the interpolation described in § 2.2.2, from $x_{p}^{n}$ and $v_{p}^{n}$, yielding $\tilde {J}_{1,i}^{n}$ and $\tilde {J}_{2,i}^{n}$ for $1\leq i\leq N_x$. Then, the field update is

(2.52)\begin{gather} \frac{ \tilde{E}_{1,i}^{n+1} - \tilde{E}_{1,i}^{n} }{\Delta t} ={-}\tilde{J}_{1,i}^{n}, \quad\text{for } 1\leq i\leq N_x, \end{gather}

(2.53)\begin{gather}\frac{ \tilde{E}_{2,i}^{n+1} - \tilde{E}_{2,i}^{n} }{\Delta t} ={-} \tilde{J}_{2,i}^{n} - \frac{ \tilde{B}_{3,i+1/2}^{n} - \tilde{B}_{3,i-1/2}^{n} }{\Delta x}, \quad\text{for } 1\leq i\leq N_x, \end{gather}

(2.54)\begin{gather}\frac{ \tilde{B}_{3,i+1/2}^{n+1} - \tilde{B}_{3,i+1/2}^{n} }{\Delta t} ={-} \frac{ \tilde{E}_{2,i+1}^{n} - \tilde{E}_{2,i}^{n} }{\Delta x}, \quad\text{for } 1\leq i\leq N_x; \end{gather}

the periodic boundary conditions are

(2.55)\begin{equation} \tilde{B}_{3,1/2}^{n} = \tilde{B}_{3,N_x+1/2}^{n} \quad\text{and}\quad \tilde{E}_{2,N_x+1}^{n} = \tilde{E}_{2,1}^{n}. \end{equation}

This spatial discretisation is Yee's lattice (Yee Reference Yee1966; Hyman & Shashkov Reference Hyman and Shashkov1999) applied to the 1D-2V setting, a choice motivated by conservation of energy considerations, see § 2.5.1.

To conclude, the fields acting on the particles, $E_{1,p}^{n+1}$, $E_{2,p}^{n+1}$ and $B_{3,p}^{n+1}$, are computed from $\tilde {E}_{1,i}^{n+1}$, $\tilde {E}_{2,i}^{n+1}$ and $\tilde {B}_{3,i+1/2}^{n+1}$ using the extrapolation described in § 2.2.2. In the numerical experiments without magnetic field, we simply set $\tilde {B}_{3,i+1/2}\equiv 0$.

2.6.2. Particle update

Given the fields acting on the particles at time $t^{n}$, $E_{1,p}^{n}$, $E_{2,p}^{n}$ and $B_{3,p}^{n}$, as well as the collision term $\mathcal {U}_{\eta,\varepsilon }[f^N] = ( \mathcal {U}_{1,\eta,\varepsilon }[f^N], \mathcal {U}_{2,\eta,\varepsilon }[f^N] )$ described in § 2.2.3, the particle update is

(2.56)\begin{gather} \frac{ x_{1,p}^{n+1} - x_{1,p}^{n} }{\Delta t} = v_{1,p}^{n+1}, \end{gather}

(2.57)\begin{gather}\frac{ v_{1,p}^{n+1} - v_{1,p}^{n} }{\Delta t} = E_{1,p}^{n} + v_{2,p}^{n} B_{3,p}^{n} - \mathcal{U}_{1,\eta,\varepsilon}[f^N](x_{p},v_{p}), \end{gather}

(2.58)\begin{gather}\frac{ v_{2,p}^{n+1} - v_{2,p}^{n} }{\Delta t} = E_{2,p}^{n} - v_{1,p}^{n} B_{3,p}^{n} - \mathcal{U}_{2,\eta,\varepsilon}[f^N](x_{p},v_{p}). \end{gather}

2.7.1. Particle initialisation

Given an initial datum $f_0(x,v_1,v_2)$, the clear approach to particle initialisation is to construct the initial particle positions and velocities $x_{p}^0 = (x_{1,p}^0)$ and $v_{p}^0 = (v_{1,p}^0, v_{2,p}^0)$ as samples $(x_{1,p}^0,v_{1,p}^0, v_{2,p}^0)$ from the distribution $f_0$.

Whenever possible, exact sampling tools are used. For instance, the marginals of $f_0$ in velocity are Gaussian mixtures in most of the numerical experiments entertained in § 3; in those cases, sampling the velocity coordinates $(v_{1,p}^0, v_{2,p}^0)$ is trivial.

In the cases where an exact sampling is not known (in § 3.1.1, and the marginal of $f_0$ in space of the rest of experiments), we use stratified sampling. The support of the target distribution is discretised in a mesh (in our case, the spatial mesh described in § 2.6), and particles are sampled uniformly on each mesh cell. The number of particles sampled on each cell is proportional to the integral of the target distribution on that cell (which we approximate numerically).

Symmetries – it is possible to enforce certain symmetries in the particle sampling. In the test cases of § 3.2, we ensure that the initial momentum of the particle solution is zero. This is done by performing the sampling with $N/2$ particles, and then generating $N/2$ additional particles through a ${\rm \pi}$ radians rotation about the velocity origin.

2.7.2. Field initialisation

The numerical experiments of § 3 all use self-consistent initialisations for the fields. Whenever present, the fields are initialised with approximations of $\boldsymbol {\nabla }_x\boldsymbol {\cdot } E_0=\rho -\rho _{\textrm {ion}}$ and $\boldsymbol {\nabla }_x\boldsymbol {\cdot } B_0=0$.

3.1.1. Order of convergence (2V)

We begin the validation of the scheme by considering a spatially homogeneous problem. In this simplified setting, our scheme would be analogous to that of Carrillo et al. (Reference Carrillo, Hu, Wang and Wu2020); however, our ‘asymmetric’ regularisation of the entropy is precisely the choice not pursued in the reference. The validation is therefore novel. Furthermore, we are validating the random initialisation described in § 2.7.1, as well as the random batch approach described in § 2.3.2.

We shall reproduce a test from Carrillo et al. (Reference Carrillo, Hu, Wang and Wu2020), which studies the relaxation to the Maxwellian in a homogeneous setting. The initial condition is a ring-like density,

(3.1)\begin{equation} f_0(v_1,v_2) = \frac{1}{\rm \pi} \exp\,(-(v_1^2+v_2^2)) (v_1^2+v_2^2), \end{equation}

and the asymptotic steady state is the Maxwellian

(3.2)\begin{equation} f_\infty(v_1,v_2) = \frac{1}{2{\rm \pi}} \exp\left(-\frac{v_1^2+v_2^2}{2}\right), \end{equation}

for every exponent $\gamma$. However, the choice of exponent determines the evolution in time of $f$. We will study the Maxwellian ($\gamma =0$) and Coulombian ($\gamma =-d$) cases numerically.

In the Maxwellian case, an exact solution exists (a derivation can be found in the appendix of Carrillo et al. (Reference Carrillo, Hu, Wang and Wu2020)), which could also be use for validation. However, we prefer to compare both the Maxwellian and Coulombian cases on an equal footing; thus we shall study the numerical convergence of the method in both cases via relative errors.

We consider the domain $v\in (-4,4)^2$, and solve the Landau equation with datum $f_0=f(0,v)$ for $t\in (0,15)$, initialising the particles with a stratified sampling. The collision strength is $C=2^{-4}$. We choose $N_{v}\in \{ 8,16,32,64 \}$ (note $N_{v_1}=N_{v_2}=N_{v}$), $N_c\in \{ 1,2,4,8 \}$, $\Delta t\in \{ 100^{-1},200^{-1},400^{-1} \}$. We solve the problems with and without random batching (respectively, $R=1$ and $R=16$).

Figure 1 shows a typical numerical solution. Figure 2 shows the relative ${L^{2}}$ error of the solutions, defined for each value of $N_{v}$ and its corresponding solution $f_{N_{v}}$ as

(3.3)\begin{equation} \text{Error}(N_{v}) := \frac{\| f_{N_{v}} - f_{N_{v}/2} \|_{{L^{2}}}}{\| f_{N_{v}} \|_{{L^{2}}}}. \end{equation}

Interestingly, Coulomb performs better than Maxwell when the number of particles per cell is very low ($N_c=1$); in the base case, $\Delta t=100^{-1}$, the Maxwellian case appears not to converge; this can be resolved by reducing $\Delta t$. However, for $N_c\geq 4$, both cases display second-order convergence with respect to $N_{v}$, and the error is lower for the Maxwellian case. The use of random batches ($R=16$) barely affects these results.

Figure 1. Typical solution in the Maxwellian case of the validation test of § 3.1.1. Here $C=2^{-4}$; $t\in (0,15)$; $N_{v}=64$; $N_c=8$; $\Delta t=10^{-2}$; $R=16$.

Figure 2. Order of convergence for the validation test of § 3.1.1. Relative ${L^{2}}$ errors of the regularised solutions $\tilde {f}^N$. Here $C=2^{-4}$; $t\in (0,15)$.

3.1.2. Landau damping (1D-2V)

We now validate the inhomogeneous scheme by studying the Landau damping (Landau Reference Landau1936; Villani Reference Villani2013) of an electric wave. We adapt a test from Medaglia et al. (Reference Medaglia, Pareschi and Zanella2023), suitably extended to the 1D-2V setting (similar tests have been used, for instance, in Crouseilles & Filbet (Reference Crouseilles and Filbet2004), Cheng, Gamba & Morrison (Reference Cheng, Gamba and Morrison2013) and Zhang & Gamba (Reference Zhang and Gamba2017)). We verify that our numerical solution behaves consistently with respect to the known physical properties. Furthermore, we study the error in the entropy increment equality. This is a necessary check because, as reported in Touati et al. (Reference Touati, Codur, Tsung, Decyk, Mori and Silva2022), collisionless PIC simulations can in fact cause the entropy to increase. We establish here that any error due to the discretisation of the transport is several orders of magnitude smaller than the typical entropy values, and therefore the effects in our simulations are solely due to the presence of the collision operator.

The test considers a distribution in Maxwellian equilibrium, with a small spatial perturbation,

(3.4)\begin{equation} f_0(x,v_1,v_2) = \frac{1+\alpha\cos(kx)}{2{\rm \pi}} \exp\left(-\frac{v_1^2+v_2^2}{2}\right), \end{equation}

with $k=2^{-1}$, and $\alpha =10^{-1}$. We will study the Coulombian case ($\gamma =-d$). Linear theory predicts that the ${L^{2}}$ norm of the electric field will oscillate, but the amplitude of the oscillation will decay exponentially with rate

(3.5)\begin{equation} \gamma_l ={-}\frac{1}{k^3}\sqrt{\frac{\rm \pi}{8}} \exp\left(-\frac{1}{2k^2}-\frac{3}{2}\right) \end{equation}

in the collisionless case. In the presence of sufficiently weak collisions, the decay is accelerated by a linear correction to $\gamma _l+{C\gamma _{l,c}}$ (Delcroix & Bers Reference Delcroix and Bers1994; Chen Reference Chen2016), where

(3.6)\begin{equation} {\gamma_{l,c} ={-}\sqrt{\frac{2}{9{\rm \pi}}}}. \end{equation}

This is, however, not the case for stronger collisions. In the hydrodynamic limit ($C\rightarrow \infty$), in the absence of a magnetic field, the limiting behaviour of system (2.1) is described by the Euler–Poisson equations, where the corresponding initial condition will result in a standing electric wave. Therefore, the actual decay rate of the electric field can be made arbitrarily slow for sufficiently strong collisions.

We consider the domain $x\in (0,2{\rm \pi} /k)$, $v\in (-4,4)^2$, and solve the Vlasov–Ampère–Landau equation for $t\in (0,10)$, initialising the particles with a combined sampling (stratified in position, Gaussian in velocity). The collision strength is $C\in \{ 0,0.01,0.1,1 \}$. We choose $N_{x}=128$, $N_{v}=32$ (note, once again, $N_{v_1}=N_{v_2}=N_{v}$), $N_c=8$, $\Delta t=50^{-1}$ and $R=32$. Each simulation employs a total of $1\,048\,576$ particles.

Figure 3 shows the decay of the first component of the electric field $(E_1)$ in the ${L^{2}}$ norm. The collisionless case behaves as expected. The very weakly collisional case ($C=0.01$) is almost indistinguishable from collisionless; so much so, that it obscures the reference $C=0$ plot. The weakly collisional case ($C=0.1$) is within the purview of the linear analysis, and it shows an accelerated decay consistent with the theoretical rate. Finally, the strongly collisional case ($C=1.0$) breaks the linear trend and displays a decay rate slower than the collisionless one, very far from the linear prediction, which is expected.

Figure 3. Numerical Landau damping in the validation test of § 3.1.2. Here $t\in (0,10)$; $N_{x}=128$; $N_{v}=32$; $N_c=8$; $\Delta t=50^{-1}$; $R=32$. Total of $1\,048\,576$ particles. Constants $\gamma _l$ and $\gamma _{l,c}$ given in § 3.1.2.

Figure 4 shows the increment of the entropy in the collisional cases, as well as the entropy conservation errors discussed in § 2.5.2. Compared with the upcoming experiments, the distribution $f$ in this test remains close to Maxwellian equilibrium as it evolves in time. Therefore, we expect the entropy increment effects to be relatively small. Nevertheless, we observe that the entropy transport error is at least three orders of magnitude smaller than the typical entropy increment; in this scenario, the evolution of the entropy is driven solely by the action of the collision operator.

Figure 4. Entropy and entropy transport error in the Landau damping validation test of § 3.1.2. Here $t\in (0,10)$; $N_{x}=128$; $N_{v}=32$; $N_c=8$; $\Delta t=50^{-1}$; $R=32$. Total of $1\,048\,576$ particles.

3.2.1. Two-stream instability (1D-2V)

We turn to study the effects of the collisional effects in more interesting problems. Here we study the two-stream instability, a known interaction between two beams of electrons which causes a vortex to appear in phase space. We observe the smoothing of the distribution due to collisional effects, as well an improvement in the conservation of energy.

The test, adapted from Medaglia et al. (Reference Medaglia, Pareschi and Zanella2023), considers two Maxwellian travelling beams, with a very small spatial perturbation,

(3.7)\begin{align} f_0(x,v_1,v_2) = \frac{1+\alpha\cos(kx)}{2{\rm \pi}} \left[ \exp\left(-\frac{(v_1-c)^2}{2}\right) + \exp\left(-\frac{(v_1+c)^2}{2}\right) \right] \exp\left(-\frac{v_2^2}{2}\right), \end{align}

with $c=2.4$, $k=5^{-1}$ and $\alpha =200^{-1}$. The electric field is initialised self-consistently. We study the Coulombian case ($\gamma =-2$).

We consider the domain $x\in (0,2{\rm \pi} /k)$, $v\in (-6,6)^2$, and solve Vlasov–Ampère–Landau for $t\in (0,50)$, initialising the particles with a combined sampling (stratified in position, Gaussian in velocity). The collision strength is $C\in \{ 0,0.01,0.02,0.04,0.08 \}$. We choose $N_{x}=256$, $N_{}{v_x}=32$, $N_{}{v_y}=4$, $N_c=16$, $\Delta t=20^{-1}$ and $R=32$. Each simulation employs a total of $524\,288$ particles.

Figure 5 shows the vortex formation in the collisionless setting, and the smoothing of the vortex as the effects of collision are strengthened; the action of the Landau operator as a nonlinear, non-local diffusion operator is evident here. For completeness, Figure 6 shows a comparison of the final state of the vortex for each collision strength. As expected, increasing $C$ interpolates from the collisionless dynamics to (nearly) a Maxwellian equilibrium.

Figure 5. Vortex formation in the two-stream instability test of § 3.2.1 (marginals of $\tilde {f}^N$). Here $t\in (0,50)$; $N_{x}=256$; $N_{v_1}=32$; $N_{v_2}=4$; $N_c=16$; $\Delta t=20^{-1}$; $R=32$. Total of $524\,288$ particles. See Bailo et al. (Reference Bailo, Carrillo and Hu2024a,Reference Bailo, Carrillo and Hub) for animations.

Figure 6. Final vortex in the two-stream instability test of § 3.2.1 (marginals of $\tilde {f}^N$). Here $t\in (0,50)$; $N_{x}=256$; $N_{v_1}=32$; $N_{v_2}=4$; $N_c=16$; $\Delta t=20^{-1}$; $R=32$. Total of $524\,288$ particles. See Bailo et al. (Reference Bailo, Carrillo and Hu2024a,Reference Bailo, Carrillo and Hub) for animations.

Figure 7 shows the evolution of the total energy, kinetic energy, electric energy and entropy of the problem. Around time $20$, we see an exponential conversion of kinetic energy to electric energy: this corresponds to the initial appearance of the vortex. The conversion saturates around time $30$, but both energies continue to increase, thus the total energy is not exactly conserved. Due to the structure-preserving properties of our discretisation, stronger collisions lead to better conservation properties, not worse. The entropy increases throughout. In the collisionless case, the entropy is essentially constant, until the point where the vortex appears, where some variation can be seen. This could be due to larger entropy conservation errors $\mathcal {C}_x$ and $\mathcal {C}_v$, or the PIC entropy increase effects described in Touati et al. (Reference Touati, Codur, Tsung, Decyk, Mori and Silva2022).

Figure 7. Energy and entropy in the two-stream instability test of § 3.2.1. Here $t\in (0,50)$; $N_{x}=256$; $N_{v_1}=32$; $N_{v_2}=4$; $N_c=16$; $\Delta t=20^{-1}$; $R=32$. Total of $524\,288$ particles. See Bailo et al. (Reference Bailo, Carrillo and Hu2024a,Reference Bailo, Carrillo and Hub) for animations.

Figure 8 shows the growth of the $L^2$ norm of the electric field in time. In the collisionless case, the growth rate in the exponential phase is $\gamma _s = 0.2258$ (Chen Reference Chen2016; Liu & Xu Reference Liu and Xu2017; Xiao & Frank Reference Xiao and Frank2021; Medaglia et al. Reference Medaglia, Pareschi and Zanella2023). It has recently been suggested (Zhou & Bellan Reference Zhou and Bellan2023) that this growth rate should be insensitive to collisional effects; however, other sources (Sydorenko et al. Reference Sydorenko, Kaganovich, Ventzek and Chen2016) expect the rate to slow down linearly for weak collisions to $\gamma _s-C\gamma _{s,c}$, for some positive constant $\gamma _{s,c}$. Our results are roughly consistent with $\gamma _{s,c}\simeq 1$.

Figure 8. Exponential growth of the electric field in the two-stream instability test of § 3.2.1. Here $t\in (0,50)$; $N_{x}=256$; $N_{v_1}=32$; $N_{v_2}=4$; $N_c=16$; $\Delta t=20^{-1}$; $R=32$. Total of $524\,288$ particles. See Bailo et al. (Reference Bailo, Carrillo and Hu2024a,Reference Bailo, Carrillo and Hub) for animations. Constants $\gamma _s$ and $\gamma _{s,c}$ given in § 3.2.1.

Figure 12 shows the energy-conservation error for the collisionless experiment, defined simply as the absolute difference between the final and initial (discrete) energy. The collisionless simulation is repeated both with a larger number of particles (doubling $N_c$ twice) and with a smaller time step (halving $\Delta t$ three times) in order to study whether the conservation error decreases. In this test, the error is essentially insensitive to the number of particles; however, it decreases linearly with $\Delta t$. This is expected, in view of the discussion of § 2.5.1 and the fact that we employ an explicit Euler integrator. As already remarked in the §§ 1 and 2.5.1, better time integrations are available, and would lead to improved conservation properties, but their use in combination with our collision operator is too costly at present. Nevertheless, we conclude that the use of more sophisticated integrators or a finer time step would lead to much better conservation in our numerics.

3.2.2. Weibel instability (1D-2V)

To conclude, we study the Weibel instability, a problem which includes the full electromagnetic effects. Once again we observe the smoothing of the distribution due to collisional effects, as well as improved conservation properties.

This test, borrowed from Cheng, Christlieb & Zhong (Reference Cheng, Christlieb and Zhong2014), considers two Maxwellian beams

(3.8)\begin{align} f_0(x,v_1,v_2) = \frac{1}{{\rm \pi}\beta} \exp\left(-\frac{v_1^2}{\beta}\right) \left[ \exp\left(-\frac{(v_2-c)^2}{\beta}\right) + \exp\left(-\frac{(v_2+c)^2}{\beta}\right) \right], \end{align}

with $c=0.3$, and $\beta =10^{-2}$. The electric field is initialised self-consistently (to zero), and the initial magnetic field is

(3.9)\begin{equation} B_3(0,x) = \alpha \sin(kx), \end{equation}

where $k=5^{-1}$, and $\alpha =10^{-3}$. We study the Coulombian case ($\gamma =-2$).

We consider the domain $x\in (0,2{\rm \pi} /k)$, $v\in (-6,6)^2$, and solve VML for $t\in (0,125)$, initialising the particles with a combined sampling (stratified in position, Gaussian in velocity). The collision strength is $C\in \{ 0,0.0001,0.0002,0.0004,0.0008 \}$. We choose $N_{x}=32$, $N_{}{v_x}=64$, $N_{}{v_y}=64$, $N_c=8$, $\Delta t=10^{-1}$ and $R=64$. Each simulation employs a total of $1\,048\,576$ particles.

Figure 9 shows the action of the instability. Without collisional effects, the initially narrow, fast beams suddenly become wider and slow, coinciding with the exponential growth of the magnetic field. The new beams are stable, though they appear noisy in the figure due to the nature of the particle approximation. The collisional cases display similar qualitative behaviour at the beginning, but eventually the two beams begin to merge. Once again, the diffusive nature of the Landau operator is evident. For completeness, figure 10 shows a comparison of the final state of the beams for each collision strength. As before, stronger collisions drive the distribution towards a Maxwellian.

Figure 9. Electron beam collapse in the Weibel instability test of § 3.2.2 (marginals of $\tilde {f}^N$). Here $t\in (0,125)$; $N_{x}=32$; $N_{v}=64$; $N_c=8$; $\Delta t=10^{-1}$; $R=64$. Total of $1\,048\,576$ particles. See Bailo et al. (Reference Bailo, Carrillo and Hu2024a,Reference Bailo, Carrillo and Hub) for animations.

Figure 10. Final beam collapse in the Weibel instability test of § 3.2.2 (marginals of $\tilde {f}^N$). Here $t\in (0,125)$; $N_{x}=32$; $N_{v}=64$; $N_c=8$; $\Delta t=10^{-1}$; $R=64$. Total of $1\,048\,576$ particles. See Bailo et al. (Reference Bailo, Carrillo and Hu2024a,Reference Bailo, Carrillo and Hub) for animations.

Figure 11 shows the evolution of the total energy, kinetic energy, electric energy, magnetic energy and entropy of the problem. Around time $30$, we see an exponential conversion of kinetic energy to magnetic energy. The conversion saturates around time $40$, when the magnetic energy begins to oscillate around a constant value, and appears convergent in time. The electric energy behaves similarly, but its magnitude is 20 times smaller, establishing this as a mainly magnetic effect. Unlike the two-stream instability, the lack of conservation here is mainly visible in the kinetic energy, which grows steadily. Once again, stronger collisions lead to better conservation properties, not worse. Just as in the two-stream instability case, the entropy conservation error due to the transport is no longer negligible once the instability develops, leading to variations in the entropy in the collisionless case. Nevertheless, the addition of collisional terms still plays an important role in the evolution of the entropy, which increases monotonically when $C$ is large enough.

Figure 11. Energy and entropy in the Weibel instability test of § 3.2.2. $t\in (0,125)$. Here ${N_{x}=32}$; $N_{v}=64$; $N_c=8$; $\Delta t=10^{-1}$; $R=64$. Total of $1\,048\,576$ particles. See Bailo et al. (Reference Bailo, Carrillo and Hu2024a,Reference Bailo, Carrillo and Hub) for animations.

Just as in the previous section, figure 12 shows the energy-conservation error for the collisionless experiment. Again, the collisionless simulation is repeated both with a larger number of particles and with a smaller time step. In this test, the error is not insensitive to the number of particles, but it nevertheless decreases linearly with $\Delta t$. Once again, we conclude that the use of more advanced time integrators or a smaller $\Delta t$ would greatly improve the conservation errors across our experiments.

Figure 12. Collisionless conservation-of-energy error on two-stream instability (§ 3.2.1) and the Weibel instability (§ 3.2.2) as a function of $\Delta t$. Starting from the original simulation parameters, $N_c$ is doubled twice, and $\Delta t$ is halved three times; other parameters remain unchanged.