## 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,

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,

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,

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$,

*a,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,

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

*a,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,

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

*a,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 (Reference Sonnendrücker2013) for a review). Particle-in-cell methods (Hockney & Eastwood Reference Hockney and Eastwood1988; Birdsall & Langdon Reference Birdsall and Langdon2018) 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.* Reference Kraus, Kormann, Morrison and Sonnendrücker2017) 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 Reference Chen, Chacón and Barnes2011; Markidis & Lapenta Reference Markidis and Lapenta2011; Lapenta Reference Lapenta2017; Ricketson & Chacón Reference Ricketson and Chacón2020; Kormann & Sonnendrücker Reference Kormann and Sonnendrücker2021), often exploiting a specific discretisation of Maxwell's equations (Yee Reference Yee1966; Hyman & Shashkov Reference Hyman and Shashkov1999). Particle-in-cell and grid-based methods exist with simultaneous energy and local charge conservation properties (Chen *et al.* Reference Chen, Chacón and Barnes2011; Taitano & Chacón Reference Taitano and Chacón2015; Taitano, Knoll & Chacón Reference Taitano, Knoll and Chacón2015*b*; Chen *et al.* Reference Chen, Chacón, Yin, Albright, Stark and Bird2020).

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 Reference Degond and Lucquin-Desreux1994; Buet & Cordier Reference Buet and Cordier1998; Crouseilles & Filbet Reference Crouseilles and Filbet2004), 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.* (Reference Taitano, Chacón, Simakov and Molvig2015*a*), Taitano, Chacón & Simakov (Reference Taitano, Chacón and Simakov2016) and Taitano *et al.* (Reference Taitano, Chacón, Simakov and Anderson2021). Implicit and asymptotic-preserving methods based on these entropy schemes have been developed in Lemou & Mieussens (Reference Lemou and Mieussens2005) and Jin & Yan (Reference Jin and Yan2011). We refer the reader to the review Dimarco & Pareschi (Reference Dimarco and Pareschi2014).

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.* Reference Buet, Cordier, Degond and Lemou1997), finite-volume methods (Taitano *et al.* Reference Taitano, Chacón, Simakov and Molvig2015*a*, Reference Taitano, Chacón and Simakov2016), fast multipole expansions (Lemou Reference Lemou1998, Reference Lemou2004) and Fourier spectral methods (Pareschi, Russo & Toscani Reference Pareschi, Russo and Toscani2000); 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 Reference Filbet and Pareschi2002; Dimarco *et al.* Reference Dimarco, Li, Pareschi and Yan2015; Zhang & Gamba Reference Zhang and Gamba2017; Hu, Jin & Shu Reference Hu, Jin and Shu2018; Li, Ren & Wang Reference Li, Ren and Wang2021), see Gamba (Reference Gamba2017) for a review. We also refer to semi-Lagrangian techniques (Ku *et al.* Reference Ku, Hager, Chang, Kwon and Parker2016).

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 Reference Takizuka and Abe1977; Nanbu Reference Nanbu1997), in the spirit of the direct simulation Monte Carlo method developed by Bird (Bird Reference Bird1994) for the Boltzmann equation, see Bobylev & Nanbu (Reference Bobylev and Nanbu2000) Caflisch *et al.* (Reference Caflisch, Wang, Dimarco, Cohen and Dimits2008) and Dimarco, Caflisch & Pareschi (Reference Dimarco, Caflisch and Pareschi2010) for generalisations and Medaglia, Pareschi & Zanella (Reference Medaglia, Pareschi and Zanella2023) 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 Reference Manheimer, Lampe and Joyce1997; Dimits *et al.* Reference Dimits, Cohen, Caflisch, Rosin and Ricketson2013); see also Rosin *et al.* (Reference Rosin, Ricketson, Dimits, Caflisch and Cohen2014) 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.* Reference Threlfall, Akers, Arter, Barnes, Barton, Cantwell, Challenor, Cook, Coveney and Dodwell2023) or XGC1 (Chang & Ku Reference Chang and Ku2008; Ku, Chang & Diamond Reference Ku, Chang and Diamond2009) 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.* Reference Miller, Churchill, Dener, Chang, Munson and Hager2021).

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.* (Reference Carrillo, Hu, Wang and Wu2020) 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.* Reference Bailo, Carrillo, Medaglia and Zanella2023) and for the multispecies Landau equation (Zonta, Pusztay & Hirvijoki Reference Zonta, Pusztay and Hirvijoki2022; Carrillo, Hu & Van Fleet Reference Carrillo, Hu and Van Fleet2023*b*). The analysis of the homogeneous method and its convergence has been presented in Carrillo, Delgadino & Wu (Reference Carrillo, Delgadino and Wu2022*a*, Reference Carrillo, Delgadino and Wu2023*a*) and Carrillo *et al.* (Reference Carrillo, Delgadino, Desvillettes and Wu2024). Moreover, a random batch technique was introduced in Carrillo, Jin & Tang (Reference Carrillo, Jin and Tang2022*b*) 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

*a*)\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}

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

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

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

*e*)\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}

*f*)\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}

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

*h*)\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.6*a*,*b*),

*i*)\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.1*h*) to increase globally, a property often called the H-Theorem. The name ‘H-Theorem’ refers to (2.1*g*), 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.1*d*) can be written in weak form as

$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,

*a–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

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 Reference Villani1998; Gualdani & Zamponi Reference Gualdani and Zamponi2017).

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

is shown by multiplying the Vlasov equation (2.1*a*) 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

using the H-theorem.

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

Momentum ($g=v$) results in

a classical computation involving Maxwell's equations leads to

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

once again, Maxwell's equations lead to

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

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

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:

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

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

(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

and define the splines

*a,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}

*a,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

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$:

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,

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,

*a,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 Reference Yee1966; Hyman & Shashkov Reference Hyman and Shashkov1999) 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$,

*a*)\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}

*b*)\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}

*c*)\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.22*c*) can be rewritten as

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

and

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

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

(2.26*a*)\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.26(Step II) Compute the values of $\boldsymbol {\nabla }_v ({\partial \mathcal {H}_{\eta,\varepsilon }}/{\partial f}) [f^N]$ at each particle,*b*)\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}(2.26(Step III) Compute the velocity field,*c*)\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}(2.26where*d*)\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}(2.26*e*)\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

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

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

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

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

(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)$.

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

(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 Tang2022*b*) 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

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

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 Tang2022*b*) 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. 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

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

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

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.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)$:

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

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

a computation that requires the vector identity

##### 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:

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,

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

The first summand can be greatly simplified, noting

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,

which will be shown to be numerically negligible.

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

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,

which will also be shown to be numerically negligible.

To summarise, we find

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

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

the periodic boundary conditions are

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.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 (Reference Bailo, Carrillo and Hu2024*a*). Videos of the simulations can be found in the permanent repository in Bailo, Carrillo & Hu (Reference Bailo, Carrillo and Hu2024*b*).

### 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.* (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,

and the asymptotic steady state is the Maxwellian

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

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.

#### 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,

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

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

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

### 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.* (Reference Medaglia, Pareschi and Zanella2023), considers two Maxwellian travelling beams, with a very small spatial perturbation,

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 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 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 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

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

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

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.

## 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.* Reference Threlfall, Akers, Arter, Barnes, Barton, Cantwell, Challenor, Cook, Coveney and Dodwell2023). 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.* Reference Ku, Hager, Chang, Kwon and Parker2016), 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.* (Reference Carrillo, Hu and Van Fleet2023*b*) did for the homogeneous scheme of Carrillo *et al.* (Reference Carrillo, Hu, Wang and Wu2020), 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 (Reference Bailo and Rey2022) for a different class of schemes and collision operators. Following Bailo *et al.* (Reference Bailo, Carrillo, Medaglia and Zanella2023), stochastic Galerkin expansions could be used to perform uncertainty quantification at the inhomogeneous level, as was done in Medaglia *et al.* (Reference Medaglia, Pareschi and Zanella2023) 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.* (Reference Miller, Churchill, Dener, Chang, Munson and Hager2021) performs a similar task at the gyrokinetic level with data from XGC1 (Chang & Ku Reference Chang and Ku2008; Ku *et al.* Reference Ku, Chang and Diamond2009).

## Acknowledgements

*Editor L.O. Silva thanks the referees for their advice in evaluating this article.*

## Funding

R.B. and J.A.C. were supported by the Advanced Grant Nonlocal-CPD (Nonlocal PDEs for Complex Particle Dynamics: Phase Transitions, Patterns and Synchronization) of the European Research Council Executive Agency (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 883363) and by the EPSRC grant EP/T022132/1 ‘Spectral element methods for fractional differential equations, with applications in applied analysis and medical imaging.’ J.H. was partially supported by AFOSR grant FA9550-21-1-0358 and DOE grant DE-SC0023164.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Dimensionless VML equations for the Coulomb case in dimension three

For the sake of completeness, we describe here the non-dimensionalisation that leads to (2.1), which roughly follows (Degond Reference Degond2007). Throughout, we assume we are in the Coulomb setting, $\gamma =-{d_v}=-3$.

To begin, we assume the typical density is $n_0$ and the typical temperature is $T_0$. The thermal speed is defined as

and the typical distribution magnitude is defined as $f_0 = n_0 v_0^{-3}$ for dimensional consistency, which results in the scaling $\,f = f_0 \bar {f}$.

The velocity coordinate will be scaled by $v = v_0 \bar {v}$. To define scales for time and space, we make use of the Debye length,

and the plasma frequency,

We can now define $t_0 = \omega _p^{-1}$ and $x_0 = \lambda _D$, and scale $t = t_0 \bar {t}$ and $x = x_0 \bar {x}$.

The typical force has magnitude $F_0:= k_B T_0 x_0^{-1}$; the fields will therefore be scaled by $E = E_0 \bar {E}$ and $B = B_0 \bar {B}$, with scales $E_0 = F_0q^{-1}$ and $B_0 = E_0 v_0^{-1}$.

We assume the Coulomb collisional cross-section normalises with a factor

see Degond (Reference Degond2007) for details. We define the collision frequency $\nu _0 := \sigma _0 n_0 v_0$, and scale the collision operator by $Q = Q_0 \bar {Q}$, where $Q_0:= \nu _0 f_0$.

We first scale Maxwell's equations. Ampère's equation $\varepsilon _0 \mu _0 \partial _t E = \boldsymbol {\nabla }_x\times B - \mu _0 J$ becomes $\partial _{\bar {t}} \bar {E} = ({c^2}/{v_0^2})\boldsymbol {\nabla }_{\bar {x}} \times \bar {B} - \bar {J}$ in the new scale, where $J = J_0 \bar {J}$ and $J_0 = q v_0 n_0$ for consistency. We make here the non-relativistic assumption $c=v_0$ to arrive at $\partial _{\bar {t}} \bar {E} = \boldsymbol {\nabla }_{\bar {x}} \times \bar {B} - \bar {J}$. Poisson's equation $\varepsilon _0 \boldsymbol {\nabla }_x\boldsymbol {\cdot } E = \rho + \rho _{\textrm {ion}}$ becomes $\boldsymbol {\nabla }_{\bar {x}} \boldsymbol {\cdot } \bar {E} = \bar {\rho } - \bar {\rho }_{\textrm {ion}}$, where $\rho = \rho _0 \bar {\rho }$, $\rho _0 = q n_0$ and $\rho _{\textrm {ion}} = - \rho _0 \bar {\rho }_{\textrm {ion}}$. Faraday's and Gauss’ equations remain unchanged in the new variables.

We now scale the Vlasov-Landau equation, $\partial _t f + v \boldsymbol {\cdot } \boldsymbol {\nabla }_x f + \frac {q}{m} (E + v\times B) \boldsymbol {\cdot } \boldsymbol {\nabla }_v f = Q[f,f]$. Under our assumptions, the equation becomes

where

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

where $C=\nu _0 \omega _p^{-1}$ is the dimensionless collision strength, obtained as the quotient of the collisional frequency and the plasma frequency.