The collisional particle-in-cell method for the Vlasov–Maxwell–Landau equations

Abstract

We introduce an extension of the particle-in-cell method that captures the Landau collisional effects in the Vlasov–Maxwell–Landau equations. The method arises from a regularisation of the variational formulation of the Landau equation, leading to a discretisation of the collision operator that conserves mass, charge, momentum and energy, while increasing the (regularised) entropy. The collisional effects appear as a fully deterministic effective force, thus the method does not require any transport–collision splitting. The scheme can be used in arbitrary dimension, and for a general interaction, including the Coulomb case. We validate the scheme on scenarios such as the Landau damping, the two-stream instability and the Weibel instability, demonstrating its effectiveness in the numerical simulation of plasma.

1. Introduction

This work introduces a fully deterministic extension of the particle-in-cell (PIC) method for the Vlasov–Maxwell equations that incorporates the effects of Landau collisions. The extension is based on a regularisation of the entropic structure of the Landau operator, and is able to preserve mass, charge, momentum and energy, while mimicking the entropy increase structure of the problem. The method can be applied in arbitrary dimension and includes all particle interaction types, including the Coulomb interaction, relevant for plasmas and nuclear fusion.

The evolution of the electrons in a plasma can be modelled by the Vlasov–Landau equation,

(1.1)\begin{equation} \partial_t f + v \boldsymbol{\cdot} \boldsymbol{\nabla}_x f + a(t,x,v) \boldsymbol{\cdot} \boldsymbol{\nabla}_v f = Q[f,f], \quad x\in \varOmega \subset {\mathbb{R}}^{d_x}, v\in {\mathbb{R}}^{d_v}, \end{equation}

where $f=f(t,x,v)$ is the number distribution function of the electrons in phase space. In the most complete models, the equation is posed in three physical dimensions (${{d_x}={d_v}=3}$), and the acceleration experienced by the electrons is derived from the Lorentz force,

(1.2)\begin{equation} a(t,x,v) = \frac{q}{m} (E(t,x) + v\times B(t,x)), \end{equation}

where $m>0$ and $q<0$ are, respectively, the mass and charge of the electron.

The long-range interactions between charged particles are described by the electric and magnetic fields, $E$ and $B$, which satisfy Maxwell's equations,

(1.3)\begin{align} \varepsilon_0 \mu_0 \partial_t E = \boldsymbol{\nabla}_x\times B - \mu_0 J,\quad \partial_t B ={-}\boldsymbol{\nabla}_x\times E,\quad \varepsilon_0 \boldsymbol{\nabla}_x\boldsymbol{\cdot} E = \rho + \rho_{\text{ion}},\quad \boldsymbol{\nabla}_x\boldsymbol{\cdot} B = 0, \end{align}

where $\varepsilon _0$ and $\mu _0$ are the permittivity and permeability of free space, and are related to the speed of light by $c^2 = (\varepsilon _0 \mu _0)^{-1}$. We assume the background ion density $\rho _{\text {ion}}$ is a given positive constant. The charge density, $\rho$, and the current density, $J$, are defined as velocity moments of the distribution $f$,

(1.4a,b)\begin{equation} \rho = q\int_{{\mathbb{R}}^{d_v}} f \,\mathrm{d} v \quad\text{and}\quad J= q\int_{{\mathbb{R}}^{d_v}} vf \,\mathrm{d} v. \end{equation}

The short-range interactions between electrons are described by $Q$, the Landau collision operator,

(1.5)\begin{equation} Q[f,f](v)=\boldsymbol{\nabla}_v\boldsymbol{\cdot} \int_{{\mathbb{R}}^{d_v}} A(v-v_*) [ f(v_*) \boldsymbol{\nabla}_v f(v) - f(v) \boldsymbol{\nabla}_{v_*} f(v_*) ] \,\mathrm{d} v_*. \end{equation}

The collisional cross-section $A$ is a symmetric and positive-semidefinite matrix given by

(1.6a,b)\begin{equation} {A(z) = \mathcal{C}_{\gamma} |z|^{\gamma+2} \varPi(z), \quad \varPi(z) = \left( I_{d_v} - \frac{z\otimes z}{|z|^2} \right), } \end{equation}

where $\mathcal {C}_{\gamma }>0$ is the collision strength. The matrix $\varPi (z)$ is the projection matrix onto the perpendicular of $z$, and $I$ is the ${d_v}\times {d_v}$ identity matrix. The exponent $\gamma$ determines the type of interaction, and is chosen in the range $-{d_v}-1\leq \gamma \leq 1$ so that the expression in (1.5) is integrable. The most physically relevant choice for plasma is $\gamma =-{d_v}=-3$, the Coulomb interaction; in this case, $\mathcal {C}_{-3} = |\log \delta | 8^{-1} {\rm \pi}^{-1} \varepsilon _0^{-2} m^{-2} q^4$, where $\log \delta$ is the so-called Coulomb logarithm.

The Landau collision operator is sometimes referred to as the Landau–Fokker–Planck operator, since it may be rewritten as a nonlinear and non-local Fokker–Planck operator,

(1.7)\begin{equation} Q[f,f] = \boldsymbol{\nabla}_v\boldsymbol{\cdot} (A_f \boldsymbol{\nabla}_v f - f a_f), \end{equation}

where the diffusion matrix $A_f$ and the drift $a_f$ are given by

(1.8a,b)\begin{align} A_f(t,v) = \int_{{\mathbb{R}}^{d_v}} A(v-v_*) f(v_*) \,\mathrm{d} v_* \quad\text{and}\quad a_f(t,v) = \int_{{\mathbb{R}}^{d_v}} A(v-v_*) \boldsymbol{\nabla}_{v_*} f(v_*) \,\mathrm{d} v_*. \end{align}

In this form, sometimes also known as the Rosenbluth–Fokker–Planck operator, we see that the Landau operator is akin to a singular, non-local and nonlinear diffusion operator.

The numerical methods for collisionless plasma (the Vlasov–Maxwell equations) can be categorised as PIC methods, finite difference/volume/element methods, and semi-Lagrangian methods (see Sonnendrücker (2013) for a review). Particle-in-cell methods (Hockney & Eastwood 1988; Birdsall & Langdon 2018) have often been favoured because they generally scale better, in view of the high dimensionality of the problem. Modern PIC methods can be designed to preserve certain structural properties; for example, GEMPIC (Kraus et al. 2017) captures well the long-time behaviour of the equation by employing a symplectic integrator; on the other hand, the exact conservation of energy can be achieved using implicit time stepping (Chen, Chacón & Barnes 2011; Markidis & Lapenta 2011; Lapenta 2017; Ricketson & Chacón 2020; Kormann & Sonnendrücker 2021), often exploiting a specific discretisation of Maxwell's equations (Yee 1966; Hyman & Shashkov 1999). Particle-in-cell and grid-based methods exist with simultaneous energy and local charge conservation properties (Chen et al. 2011; Taitano & Chacón 2015; Taitano, Knoll & Chacón 2015b; Chen et al. 2020).

The main difficulty in the numerical simulation of collisionless plasma is the handling of the small-scale effects that arise from the electromagnetic fields, as they often lead to fine structure filamentation in phase space. Adding collisional effects will ease this difficulty while introducing a bigger one: the curse of dimensionality. Ultimately, the discretisation of the diffusive terms in six dimensions is very challenging from a computational efficiency perspective.

Several numerical methods have been proposed to discretise the Landau collision operator (1.5), both in the spatially homogeneous and inhomogeneous settings. We classify them here as deterministic and stochastic methods. Among the deterministic class, in the homogeneous setting, we highlight the classical entropy schemes (Degond & Lucquin-Desreux 1994; Buet & Cordier 1998; Crouseilles & Filbet 2004), which preserve the conservation and dissipation properties of the Landau equation, and thus lead to the correct Maxwellian stationary states. Methods that discretise the Landau operator in the form (1.7) are given in Taitano et al. (2015a), Taitano, Chacón & Simakov (2016) and Taitano et al. (2021). Implicit and asymptotic-preserving methods based on these entropy schemes have been developed in Lemou & Mieussens (2005) and Jin & Yan (2011). We refer the reader to the review Dimarco & Pareschi (2014).

Yet, the efficient approximation of the Landau operator remains a major challenge, even in the spatially homogeneous setting. The non-local nature of the collision operator leads to a quadratic complexity $O(N^2)$, where $N$ is the number of discretisation elements. Faster algorithms that attempt to reduce this cost to $O(N\log N)$ have been proposed, including multigrid algorithms (Buet et al. 1997), finite-volume methods (Taitano et al. 2015a, 2016), fast multipole expansions (Lemou 1998, 2004) and Fourier spectral methods (Pareschi, Russo & Toscani 2000); the latter deserves special attention, since it reduces the complexity to $O(N\log N)$ through the fast Fourier transform by exploiting the convolutional properties of the Landau operator. The spectral method has been coupled with discretisation methods for the transport part to treat the inhomogeneous problems (Filbet & Pareschi 2002; Dimarco et al. 2015; Zhang & Gamba 2017; Hu, Jin & Shu 2018; Li, Ren & Wang 2021), see Gamba (2017) for a review. We also refer to semi-Lagrangian techniques (Ku et al. 2016).

Nevertheless, stochastic or Monte Carlo approaches remain the most frequently used tools among practitioners to simulate Coulomb collisions in an inhomogeneous setting, and we highlight two major classes. The first is based on binary collisions (Takizuka & Abe 1977; Nanbu 1997), in the spirit of the direct simulation Monte Carlo method developed by Bird (Bird 1994) for the Boltzmann equation, see Bobylev & Nanbu (2000) Caflisch et al. (2008) and Dimarco, Caflisch & Pareschi (2010) for generalisations and Medaglia, Pareschi & Zanella (2023) for a recent application in uncertainty quantification. The second is based on the drift-diffusion formulation of the Landau operator given in (1.7) (Manheimer, Lampe & Joyce 1997; Dimits et al. 2013); see also Rosin et al. (2014) for a multilevel Monte Carlo extension. The advantages and drawbacks of these stochastic approaches are well known: they are physically motivated and easy to implement, but their convergence is slow, they require many realisations due to statistical noise, and they cannot preserve the entropy structure of the problem.

The design of more efficient discretisation techniques for the Landau operator is a matter of utmost practical relevance, as many of the methods described above are prohibitively costly and cannot be used directly in real-world applications, such as the design and on-line control of modern nuclear fusion reactors. Mature plasma simulation tools such as NESO (Threlfall et al. 2023) or XGC1 (Chang & Ku 2008; Ku, Chang & Diamond 2009) typically include collisional effects through Monte Carlo approaches or mesh-based methods. The use of surrogate models trained on synthetic data for the collisional operator has recently been proposed (Miller et al. 2021).

This work introduces the collisional particle-in-cell (C-PIC) method; a deterministic extension of PIC methods also able to approximate the Landau collisional effects in the Vlasov–Maxwell–Landau (VML) equations. To construct the method, we generalise a recent particle approximation of the Landau operator introduced by Carrillo et al. (2020) in the homogeneous setting. The method exploits the variational properties of the Landau operator to propose a regularised entropic structure which is naturally discretised by particles, in a way that is fully compatible with any PIC method. The increment of the (regularised) entropy is preserved at the discrete level, as are the conservations of mass, charge, momentum and energy.

This approach has proven robust and flexible in the homogeneous setting, and has already been used for uncertainty quantification in the Landau equation (Bailo et al. 2023) and for the multispecies Landau equation (Zonta, Pusztay & Hirvijoki 2022; Carrillo, Hu & Van Fleet 2023b). The analysis of the homogeneous method and its convergence has been presented in Carrillo, Delgadino & Wu (2022a, 2023a) and Carrillo et al. (2024). Moreover, a random batch technique was introduced in Carrillo, Jin & Tang (2022b) that can considerably reduce the computational cost of the method, while retaining all structural properties.

The generalisation to the inhomogeneous setting presented in this work presents a crucial departure from the previous works: the entropic structure is now regularised at a global level, and the resulting regularised Landau operator delocalises in space. Since the particle approximation of the operator is performed in phase space, no splitting of transport and collisions is required; the collisional effects appear simply as an effective force, alongside the Lorentz force of the PIC approach. There is therefore no stochasticity and no splitting in our scheme. Moreover, the spatial structure, combined with the random batch approach, can be leveraged for an efficient implementation of the method that performs comparably to classical PIC methods.

The rest of this work is organised as follows. In § 2, we recall the physical properties and variational structure of the VML equations, introduce the C-PIC method and discuss its properties. In § 3, we perform a range of numerical simulations (${d_x}=1$ and ${d_v}=2$) to validate the method and demonstrate its effectiveness, including explorations of the collisional effects on the Landau damping, the two-stream instability and the Weibel instability. We conclude in § 4, where we also present the outlook of this work.

2. The collisional particle-in-cell method

This section introduces the C-PIC method and discusses its properties.

2.1. The VML equations

The VML equations ((1.1), (1.3) and (1.5)) may, after non-dimensionalisation, be written as

(2.1a)\begin{gather} \partial_t f + v \boldsymbol{\cdot} \boldsymbol{\nabla}_x f + (E+v\times B) \boldsymbol{\cdot} \boldsymbol{\nabla}_v f = Q[f,f], \quad x \in \varOmega\subseteq {\mathbb{R}}^{d_x}, v \in {\mathbb{R}}^{d_v}, \end{gather}

(2.1b)\begin{gather}{ \partial_t E = \boldsymbol{\nabla}_x \times B - J, } \quad { \partial_t B={-}\boldsymbol{\nabla}_x\times E, } \end{gather}

(2.1c)\begin{gather}{ \boldsymbol{\nabla}_x \boldsymbol{\cdot} E = \rho-\rho_{\text{ion}}, } \quad { \boldsymbol{\nabla}_x \boldsymbol{\cdot} B=0, } \end{gather}

(2.1d)\begin{gather}Q[f,f] = \boldsymbol{\nabla}_v\boldsymbol{\cdot}(f\mathcal{U}[f]), \end{gather}

(2.1e)\begin{gather}\mathcal{U}[f](x,v) = \int_{{\mathbb{R}}^{d_v}} A (v - v_*) b [f] (x,v,v_*) f (x,v_*) \,\mathrm{d} v_*, \end{gather}

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

(2.1g)\begin{gather}{H [f](x) ={-}S [f](x),} \end{gather}

(2.1h)\begin{gather}{S [f](x) ={-}\int_{{\mathbb{R}}^{d_v}} f \log f \,\mathrm{d} v.} \end{gather}

Details of the non-dimensionalisation are given in the Appendix (A). Note that some terms have changed sign because the negative fundamental charge $q$ has been normalised. We will only consider quadrangular domains $\varOmega$ equipped with periodic boundary conditions. The cross-section matrix $A$ is the non-dimensionalised version of (1.6a,b),

(2.1i)\begin{equation} A(z) = C |z|^{\gamma+2} \varPi(z), \quad \varPi(z) = \left( I_{d_v} - \frac{z\otimes z}{|z|^2} \right), \end{equation}

where $C>0$ is the dimensionless collision strength, $\varPi (z)$ is the projection matrix onto $z^\perp$, and $I_{d_v}$ is the identity matrix in ${d_v}$ dimensions. By construction, the matrix $A$ is positive semidefinite.

Equation (2.1) possesses a wealth of physical properties: it conserves the mass, charge, momentum, and total energy of the system. It also causes the thermodynamic entropy (2.1h) to increase globally, a property often called the H-Theorem. The name ‘H-Theorem’ refers to (2.1g), the information-theoretical entropy $H$, which differs from $S$ by a sign. In the sequel, when we refer to entropy, we refer to $S$. We shall briefly recall these classical results for the sake of a complete exposition.

The Landau collision operator (2.1d) can be written in weak form as

(2.2)\begin{align} \int_{{\mathbb{R}}^{d_v}} g(v) Q[f,f] \,\mathrm{d} v & ={-} \frac{1}{2} \iint_{{\mathbb{R}}^{d_v} \times {\mathbb{R}}^{d_v}} ( \boldsymbol{\nabla}_v g(v) - \boldsymbol{\nabla}_{v_*} g(v_*) )\nonumber\\ & \quad \boldsymbol{\cdot} A (v - v_*) b [f] (v,v_*) f f_* \,\mathrm{d} v \,\mathrm{d} v_*; \end{align}

$f$ and $f_*$ stand, respectively, for $f(v)$ and $f(v_*)$, and the local dependence in $x$ has been omitted for simplicity. The expression vanishes for $g=1$, $g=v$ and $g=\tfrac {1}{2}|v|^2$; these correspond, respectively, to the conservation of mass (and charge), momentum and kinetic energy. The result is immediate for mass and momentum; for kinetic energy, one uses the fact that $A$ projects $b$ onto the perpendicular to $v-v_*$. Summarising,

(2.3a–c)\begin{equation} \int_{{\mathbb{R}}^{d_v}} Q[f,f] \,\mathrm{d} v = 0,\quad \int_{{\mathbb{R}}^{d_v}} v Q[f,f] \,\mathrm{d} v = 0,\quad \int_{{\mathbb{R}}^{d_v}} \frac{v^2}{2} Q[f,f] \,\mathrm{d} v = 0. \end{equation}

The weak form is also useful to prove the H-theorem: one chooses the test function $g={\partial H}/{\partial f}=\log f+1$ in order to arrive at

(2.4)\begin{align} & \int_{{\mathbb{R}}^{d_v}} (\log f+1) Q[f,f] \,\mathrm{d} v \nonumber\\ & \quad ={-} \frac{1}{2} \iint_{{\mathbb{R}}^{d_v} \times {\mathbb{R}}^{d_v}} b [f] (v,v_*) \boldsymbol{\cdot} A (v - v_*) b [f] (v,v_*) f f_* \,\mathrm{d} v \,\mathrm{d} v_* \leq 0, \end{align}

where the sign follows from the positive semidefinite property of the matrix $A$.

A corollary to these properties is that the kernel of $Q$ is precisely the set of Maxwellian distributions parametrised by density, momentum and energy, as is the case for the Boltzmann equation (Villani 1998; Gualdani & Zamponi 2017).

At the level of the full VML system, the increment of the total entropy,

(2.5)\begin{equation} {\mathcal{S}[f] := \int_\varOmega S[f](x) \,\mathrm{d}\kern0.07em x ={-}\iint_{\varOmega \times {\mathbb{R}}^{d_v}} f \log f \,\mathrm{d} v \,\mathrm{d}\kern0.07em x,} \end{equation}

is shown by multiplying the Vlasov equation (2.1a) by the test function $g={\partial \mathcal {H}}/{\partial f}=\log f+1$, integrating over the phase space and using the boundary conditions/behaviour at infinity of $f$ to arrive at

(2.6)\begin{equation} {\frac{\mathrm{d} }{\mathrm{d} t} \mathcal{S}[f] ={-}\iint_{\varOmega\times{\mathbb{R}}^{d_v}} (\log f+1) Q[f,f] \,\mathrm{d} v \,\mathrm{d}\kern0.07em x \geq 0,} \end{equation}

using the H-theorem.

The global conservation properties are similarly derived by multiplying (2.1a) by the correct test function and integrating. Mass ($g=1$) immediately leads to

(2.7)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \iint_{\varOmega\times{\mathbb{R}}^{d_v}} f \,\mathrm{d} v \,\mathrm{d}\kern0.07em x = 0. \end{equation}

Momentum ($g=v$) results in

(2.8)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \iint_{\varOmega\times{\mathbb{R}}^{d_v}} v f \,\mathrm{d} v \,\mathrm{d}\kern0.07em x = \int_{\varOmega} (\rho E + J\times B) \,\mathrm{d}\kern0.07em x; \end{equation}

a classical computation involving Maxwell's equations leads to

(2.9)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \left( \iint_{\varOmega\times{\mathbb{R}}^{d_v}} v f \,\mathrm{d} v \,\mathrm{d}\kern0.07em x + \int_\varOmega E\times B \,\mathrm{d}\kern0.07em x \right) = \int_{\varOmega} \rho_{\text{ion}} E \,\mathrm{d}\kern0.07em x; \end{equation}

in the case where the background ion density $\rho _{\text {ion}}$ is constant and there is no magnetic field, this quantity is exactly zero. Similarly, energy ($g=\tfrac {1}{2}|v|^2$) leads to

(2.10)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \frac{1}{2} \iint_{\varOmega\times{\mathbb{R}}^{d_v}} |v|^2 f \,\mathrm{d} v \,\mathrm{d}\kern0.07em x = \int_{\varOmega} J\cdot E \,\mathrm{d}\kern0.07em x; \end{equation}

once again, Maxwell's equations lead to

(2.11)\begin{equation} \frac{\mathrm{d} }{\mathrm{d} t} \left( \frac{1}{2} \iint_{\varOmega\times{\mathbb{R}}^{d_v}} |v|^2 f \,\mathrm{d} v \,\mathrm{d}\kern0.07em x + \frac{1}{2} \int_\varOmega (|E|^2 + |B|^2) \,\mathrm{d}\kern0.07em x \right) = 0. \end{equation}

This computation is described in more detail in § 2.5.

2.2. Description of the method

To devise the C-PIC method for (2.1), we shall seek an approximate particle solution of the form

(2.12)\begin{equation} f^N(t,x,v) = \sum_{p=1}^{N} w_{p} \delta(x - x_{p}(t)) \delta(v - v_{p}(t)), \end{equation}

where $N$ is the number of particles, and where $w_{p}$, $x_{p}(t)$ and $v_{p}(t)$ are, respectively, the $p\text {th}$ particle's weight, position and velocity. Following the characteristics, we require $x_{p}(t)$ and $v_{p}(t)$ to solve

(2.13)\begin{equation} \left\{ \begin{aligned} & \displaystyle \frac{\mathrm{d}\kern0.07em x_{p}}{\mathrm{d} t} = v_{p},\\ & \displaystyle \frac{\mathrm{d} v_{p}}{\mathrm{d} t} = E(t,x_{p}) + v_{p} \times B(t,x_{p}) - \mathcal{U}[f^N](x_{p},v_{p}), \end{aligned} \right. \end{equation}

where $E(t,x_{p})$ and $B(t,x_{p})$ are the electromagnetic fields acting on the particle, and $\mathcal {U}[f^N](x_{p},v_{p})$ is the effective force arising from the collision term. In the absence of collisions, system (2.13) will reduce to the classical PIC method for the Vlasov–Maxwell system ((1.1),(1.3)).

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.

Remark 2.1 (Staggered grids)

While the presentation in this section places $\tilde J$, $\tilde E$ and $\tilde B$ on the same mesh for simplicity, this is not required. In particular, staggered grids such as Yee's lattice (Yee 1966; Hyman & Shashkov 1999) can be used instead. This is exactly what is used here, see § 2.6 for the full implementation details.

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 2006).

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. (2020) 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. Computational optimisation of the method

In order to maximise the computational performance of the method, we describe two optimisations that we employ in our implementation of the scheme. The first, a cell list, does not alter the method. The second, the random batch method, does alter the scheme, but does not significantly affect its accuracy. In both cases, the structural properties of the scheme persist.

2.3.1. Cell list optimisation

Cell lists (Allen & Tildesley 1990, 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 2020; Carrillo et al. 2022b) 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. (2022b) for details, where the random batch method was used to accelerate the particle method for the homogeneous Landau equation (Carrillo et al. 2020). Furthermore, the use of batches does not significantly affect the accuracy of the method, see § 3.1.

2.4. Properties of the regularised collision operator

This section discusses the structural properties of the discrete collision operator (2.26).

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. Global properties of the C-PIC method

This section discusses the properties of the full C-PIC method (2.13). By construction, the discretisation of the collision operator preserves the conservation properties and H-theorem, as discussed in § 2.4. Therefore, C-PIC inherits the properties of the underlying PIC implementation.

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 1966; Hyman & Shashkov 1999) 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. 2011; Markidis & Lapenta 2011; Crouseilles, Einkemmer & Faou 2015; Lapenta 2017; Kormann & Sonnendrücker 2021). 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.6 (Conservation of momentum)

A similar analysis can be performed for the conservation of momentum. Similarly, a specific chain rule in time and vector identity in space are required to hold in order to derive the correct conservation identity. In general, these will not hold, even in the discretisations that do conserve the energy. We are not aware of any discretisation capable of exactly conserving energy and momentum simultaneously.

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 (1992), Esirkepov (2001), Chen et al. (2011) and Chen et al. (2020).

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. Full discretisation and time stepping

This section describes the full discretisation of the field solver and our time stepping scheme. We limit the implementation to the 1D-2V setting (one dimension in space, two dimensions in velocity; ${d_x}=1$, ${d_v}=2$) for the sake of computation, although the method as described in § 2.2 also applies in higher dimensions.

We solve the VML equations over a time interval $(0,T)$, where $T>0$, and discretise it with a uniform step $\Delta t$. We describe here the update from time $t^{n}$ to $t^{n+1}$, where $t^{n}:= n\Delta t$.

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 1966; Hyman & Shashkov 1999) 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. Initialisation from continuous data

This section discusses the initialisation of the C-PIC method (2.13). We again limit the discussion to the 1D-2V setting, though our approach can be generalised to higher dimensions.

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. Numerical experiments

In this section we validate the C-PIC method and demonstrate its effectiveness in dealing with a range of collisional plasma simulations. Interactive versions of the simulations presented in this section are available online in Bailo, Carrillo & Hu (2024a). Videos of the simulations can be found in the permanent repository in Bailo, Carrillo & Hu (2024b).

3.1. Validation

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. (2020); 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. (2020), 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. (2020)), 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 1936; Villani 2013) of an electric wave. We adapt a test from Medaglia et al. (2023), suitably extended to the 1D-2V setting (similar tests have been used, for instance, in Crouseilles & Filbet (2004), Cheng, Gamba & Morrison (2013) and Zhang & Gamba (2017)). 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. (2022), 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 1994; Chen 2016), 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. Study of collisional effects

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. (2023), 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. (2024a,b) 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. (2024a,b) 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. (2022).

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. (2024a,b) 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 2016; Liu & Xu 2017; Xiao & Frank 2021; Medaglia et al. 2023). It has recently been suggested (Zhou & Bellan 2023) that this growth rate should be insensitive to collisional effects; however, other sources (Sydorenko et al. 2016) 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. (2024a,b) 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 (2014), 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. (2024a,b) 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. (2024a,b) 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. (2024a,b) 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.

4. Conclusion and outlook

In this work, we have introduced a generalisation of the PIC method able to simulate the Landau collisions in a spatially inhomogeneous plasma. We have derived the method from the gradient-flow structure of the Landau equation, preserving a variational structure, and retaining all conservation properties for the collisional update. We have validated the numerical method, and presented collisional simulations of phenomena such as the Landau damping, the two-stream instability, and the Weibel instability, which demonstrate its ability to simulate a collisional plasma.

Several extensions of the work could be considered. A clear future task is to perform 3D-3V (three dimensions in space, three dimensions in velocity) simulations of a collisional plasma with full electromagnetic effects using C-PIC; this remains a tough computational challenge, and will likely involve the implementation of C-PIC within a mature plasma simulation ecosystem, such as NESO (Threlfall et al. 2023). It would also be of interest to study whether the techniques employed in this work could be applied to the gyrokinetic Landau operator (Ku et al. 2016), as a way to reduce the overall computational cost. Furthermore, a natural objective is the extension of the method to a two species setting, as Carrillo et al. (2023b) did for the homogeneous scheme of Carrillo et al. (2020), though the large ion–electron mass ratio will likely lead to challenging restrictions on the time integrator; a similar challenge was addressed in Bailo & Rey (2022) for a different class of schemes and collision operators. Following Bailo et al. (2023), stochastic Galerkin expansions could be used to perform uncertainty quantification at the inhomogeneous level, as was done in Medaglia et al. (2023) for the Vlasov–Poisson–BGK equations. The use of C-PIC as a source of synthetic data to train a surrogate model is not out of the question; Miller et al. (2021) performs a similar task at the gyrokinetic level with data from XGC1 (Chang & Ku 2008; Ku et al. 2009).