A structure-preserving particle discretisation for the Lenard–Bernstein collision operator

Abstract

Collisions are an important dissipation mechanism in plasmas. When approximating collision operators numerically, it is important to preserve their mathematical structure in order to retain the laws of thermodynamics at the discrete level. This is particularly challenging when considering particle methods. A simple but commonly used collision operator is the Lenard–Bernstein operator, or its modified energy- and momentum-conserving counterpart. In this work, we present a macro-particle discretisation of this operator that is provably energy and momentum preserving.

1. Introduction

Structure-preserving numerical methods aim at preserving certain properties of a system of equations exactly at the discrete level. Some examples for properties of interest are symmetries and conservation laws, Lagrangian or Hamiltonian structure, or compatibility with the laws of thermodynamics. Preserving such structures is typically found to be advantageous for accuracy and robustness of numerical schemes, especially for strongly nonlinear problems and long-time simulations (Hairer & Wanner 2006). This has also been recognised in plasma physics, and the last decade has seen striking efforts towards the development of structure-preserving algorithms for problems such as magnetohydrodynamics, the Vlasov–Poisson and the Vlasov–Maxwell system (see e.g. Morrison (2017) and references therein). So far, most work has focused on dissipationless systems, with dissipative systems, such as collisional kinetic systems, being considered only more recently. However, dissipative effects, although often weak, are important for the correct simulation of physical behaviour over long simulation times. Sometimes, the neglect of dissipative effects can cause numerical problems, e.g. when small structures emerge that cannot be resolved by the computational mesh. In many cases, these structures are unphysical because dissipation would prevent their emergence in the first place. Thus, the inclusion of dissipation is important not only for physical correctness but also because it can aid numerical robustness.

Work on the structure-preserving discretisation of Vlasov-like equations has mainly focused on particle-based methods. In recent years, many authors have worked on the ideal (non-dissipative) part of the problem, including Chen, Chacón & Barnes (2011), Markidis & Lapenta (2011), Squire, Qin & Tang (2012), Evstatiev & Shadwick (2013), Qin et al. (2016), Burby (2017), Kraus et al. (2017), Zhang & Gamba (2017), Campos Pinto, Kormann & Sonnendrücker (2022).

After the discretisation of the ideal problem was well understood, focus shifted towards the structure-preserving discretisation of the collisional (dissipative) part. While early work focused on grid-based methods (see e.g. Yoon & Chang 2014; Taitano et al. 2015; Hirvijoki & Adams 2017; Kraus & Hirvijoki 2017; Shiroto & Sentoku 2019), structure-preserving discretisations for collision operators with particles have been considered more recently. Hirvijoki, Kraus & Burby (2018) considered an approach where the weights of the marker particles are varied, instead of their velocities. Carrillo et al. (2020) and Hirvijoki (2021) used finite-sized marker particles to discretise the Landau operator. Mollén et al. (2021) and Pusztay, Knepley & Adams (2022) focused on projection/interpolation techniques for computing collision operators for particles. An alternative approach is that of Tyranowski (2021), which treats the collisions as a stochastic process, effectively modelling their underlying microscopic behaviour rather than the resultant macroscopic effects modelled by various collision operators.

The aim of this work is to provide a proof-of-concept for an alternative approach to structure-preserving particle methods for collisions, specifically for a conservative version of the operator of Lenard & Bernstein (1958). We employ the deterministic particle method (Degond & Mustieles 1990; Chertock 2017) to obtain dynamical equations for the particle velocities. These are then regularised by evaluating the collisional flux on a smoothened representation of the distribution function that is obtained from a projection of the particles onto a spline basis. We show that this semi-discretisation maintains momentum and energy conservation exactly. A midpoint discretisation in time is employed in order to retain these conservation properties also at the fully discrete level.

The structure of the paper is as follows. In § 2, we detail the derivation of the conservative Lenard–Bernstein operator. In § 3, the semi-discretisation of the operator is presented and its conservation properties are verified. Section 4 describes a possible time discretisation by the midpoint rule and proves that it retains the desired conservation properties. Section 5 shows several numerical tests and examples for the one-dimensional case, including some convergence results and verification of momentum and energy conservation. Finally, the paper is concluded with a discussion of current and future work.

2. The conservative Lenard–Bernstein operator

The Lenard–Bernstein collision operator (Lenard & Bernstein 1958) is a scalar operator of the form

(2.1)\begin{equation} C[f] = \nu\frac{\partial}{\partial v} \cdot \left(\frac{\partial f}{\partial v} + vf \right), \end{equation}

where $f:\mathbb {R}^d \times [0,\infty ) \to \mathbb {R}$ is the single-particle distribution function, $v \in \mathbb {R}^d$ is the velocity and $\nu$ is the collision frequency, which is assumed to be constant in time. In most applications, such a collision operator is coupled to the Vlasov–Poisson or Vlasov–Maxwell equations, and so the distribution function would also depend on the position variables. Here, however, we will ignore this dependency as we study the collision operator independently of the ideal dynamics and this operator acts purely in velocity space. Specifically, we solve the following differential equation:

(2.2)\begin{equation} \partial_t f = C[f] . \end{equation}

The Lenard–Bernstein operator is applicable in velocity dimensions $d = 1,2,3$, though collision operators which describe more physics effects, such as the Landau operator, may be preferred in two and three dimensions in order to allow, for example, interchange of momentum between different components. The steady-state solution to (2.2) for this operator is a $d$-dimensional Gaussian distribution.

2.1. Construction of the conservative operator

The Lenard–Bernstein operator (2.1) preserves mass density, the zeroth moment of the distribution function. However, it does not preserve momentum density or energy density, which are the first and second moments, respectively, and whose conservation is crucial for obtaining physically correct results in numerical simulations. In order to enforce conservation of these quantities, we follow Kraus (2013) (see also Filbet & Sonnendrücker 2003) and modify the operator through an expansion, as follows:

(2.3)\begin{equation} C[f] = \nu\frac{\partial}{\partial v} \cdot \left(\frac{\partial f}{\partial v} + A_1f + A_2 vf \right) . \end{equation}

In the following, we will see that the coefficients $A_m$ are functions of the moments of the distribution function $f$. In general, preserving $k$ moments of the distribution function will require an expansion including $k$ terms in the operator. The coefficients are then computed by requiring (2.2) to obey the respective conservation laws, which here are for momentum and energy. Specifically, conservation of the $k$th moment requires the following condition to hold:

(2.4)\begin{equation} \int v^k C[f] \,\text{d}{v} = \nu \int v^k \left[ \frac{\partial}{\partial v} \cdot \left(\frac{\partial f}{\partial v} + A_1f + A_2 vf \right)\right] \text{d}{v} = 0. \end{equation}

Integrating this by parts, we obtain the following condition:

(2.5)\begin{equation} \nu \left[ v^k\left(\frac{\partial f}{\partial v} + A_1f + A_2 vf \right) \right]_{-\infty}^{+\infty} - k \nu \int v^{k-1}\left(\frac{\partial f}{\partial v} + A_1f + A_2 vf \right) \text{d}{v} = 0 .\end{equation}

Without loss of generality, we assume that $f$ and $\partial f/ \partial v$ approach zero as $|v| \to \infty$, so that the first term in (2.5) is zero and we obtain the following condition:

(2.6)\begin{equation} \int v^{k-1}\left(\frac{\partial f}{\partial v} + A_1f + A_2 vf \right) \text{d}{v} = 0, \end{equation}

for $k = 1, 2.$ Integrating the first term by parts once again, this equation becomes

(2.7)\begin{equation} \int \left[ (k -1) v^{k-2} - A_1v^{k-1} - A_2 v^k \right] f \,\text{d}{v} = 0, \end{equation}

where the assumption that $f$ approaches zero as $|v|$ tends to infinity has been utilised once more. Writing the moments as $M_m[f] = \int v^m f \,\text {d}{v}$, we obtain the following conditions:

(2.8)\begin{equation} (k - 1) M_{k-2} = A_1 M_{k-1} + A_2 M_{k}, \quad k = 1, 2 . \end{equation}

These conditions provide a linear system of equations that can be solved for the coefficients $A_1$, $A_2$:

(2.9)\begin{equation} \left.\begin{array}{c@{}} A_1M_0 + A_2M_1 = 0, \\ A_1M_1 + A_2M_2 = M_0. \end{array}\right\}\end{equation}

The solution to the system of equations in (2.9) is

(2.10)\begin{gather} A_1 = \frac{M_0M_1}{M_0M_2 - M_1^2} = \frac{-u}{\varepsilon - u^2}, \end{gather}

(2.11)\begin{gather}A_2 = \frac{-M_0^2}{M_0M_2 - M_1^2} = \frac{1}{\varepsilon - u^2}, \end{gather}

where $nu$ and $n\varepsilon$ are the momentum and energy density, respectively, and are related to the moments as follows:

(2.12a–c)\begin{equation} n = M_0 = \int f \,\text{d}{v}, \quad nu = M_1 = \int vf \,\text{d}{v}, \quad n\varepsilon = M_2 = \int v^2f \,\text{d}{v}. \end{equation}

Let us note that here, $n$, $u$ and $\varepsilon$ are just constants. However, in the general Vlasov case, these quantities have a spatial dependency. Upon inserting the expressions for $A_1, A_2$ into (2.3), we obtain the following operator:

(2.13)\begin{equation} C[f] = \nu\frac{\partial}{\partial v} \cdot \left(\frac{\partial f}{\partial v} + \frac{v-u}{\varepsilon - u^2} f \right),\end{equation}

which can be seen as a conservative version of the Lenard–Bernstein operator (2.1). This is the same operator as the one obtained by Filbet & Sonnendrücker (2003) and is closely related to the operator studied in Hakim et al. (2020).

2.2. The H-theorem

We follow a similar strategy to Hakim et al. (2020) to demonstrate that the continuous operator satisfies the H-theorem. Let us denote the collisional flux by

(2.14)\begin{equation} F[f]=\frac{\partial f}{\partial v} + A_1 f + A_2 vf , \end{equation}

so that

(2.15)\begin{equation} \frac{\partial f}{\partial t} = C[f] = \nu \frac{\partial}{\partial v} \cdot F[f]. \end{equation}

The change in entropy is then

(2.16)\begin{align} \frac{{\rm d}S}{{\rm d}t} & = \frac{{\rm d}}{{\rm d}t} \int f \log f \,\text{d}{v} = \int (1 + \log f ) \frac{\partial f }{\partial t} \,\text{d}{v} = \nu \int (1 + \log f ) \frac{\partial}{\partial v} \cdot F \,\text{d}{v} \nonumber\\ & ={-} \nu \int \frac{1}{f} \frac{\partial f}{\partial v} \cdot F\, \text{d}{v} ={-} \nu \int \frac{1}{f}|F|^2 \,\text{d}{v} + \int (A_1 + A_2 v) F \,\text{d}{v} , \end{align}

where integration by parts, the assumption that $F$ goes to zero as $|v| \to \infty$ and the substitution $\partial f / \partial v = F - A_1 f - A_2 v f$ were used. By (2.6), the second term in the last expression is designed to vanish, namely

(2.17)\begin{equation} \int (A_1 + A_2 v) F \,\text{d}{v} = 0. \end{equation}

Hence, the entropy is monotonically dissipated (provided that $f$ is non-negative):

(2.18)\begin{equation} \frac{{\rm d}S}{{\rm d}t} ={-} \int \frac{1}{f} F^2 \,\text{d}{v} \leq 0. \end{equation}

The entropy is stationary when $F[f_\textrm {eq}]=0$, that is when $f_\textrm {eq}$ solves the partial differential equation (PDE)

(2.19)\begin{equation} \frac{\partial f_{\rm eq}}{\partial v} ={-}\frac{v-u}{\varepsilon-u^2}f_{\rm eq}. \end{equation}

The (unique) solution with conditions $f\overset {|v|\to \infty }{\longrightarrow }0$ and $\int f \,\textrm {d}{v}=n$ is a $d$-dimensional shifted Gaussian distribution with mean $u$ and variance $\varepsilon -u^2$,

(2.20)\begin{equation} f_{\rm eq}(v) = \frac{n}{(2{\rm \pi}(\varepsilon-u^2))^{d/2}} \, {\rm e}^{-({1}/{2})(v-u)^2/(\varepsilon-u^2)} . \end{equation}

Thus, the continuous operator of (2.13) satisfies the H-theorem.

3. Semi-discrete operator

In order to discretise the conservative Lenard–Bernstein operator in velocity space, we need to introduce a second representation of the distribution function. The particle representation with Dirac delta distributions, which is usually used to solve the ideal part, is not differentiable and thus cannot be used to evaluate the collisional part. Previous works regularised the collision operator by using finite sized particles (Carrillo et al. 2020; Hirvijoki 2021). Here, we explore a different approach based on finite element or spline spaces of sufficient regularity.

The particle distribution function is given by

(3.1)\begin{equation} f_p(v,t) = \sum \limits_\alpha w_\alpha \delta(v - v_\alpha (t)), \end{equation}

where $\{v_\alpha (t)\}_{\alpha =1}^N$ are the particle velocities which evolve over time, and $N$ is the number of particles. As the particle distribution function $f_p$ is non-differentiable, we use an $L^2$ projection of $f_p$ onto a set of differentiable basis functions $\{\varphi _j\}_{j=1}^M$ for $M \ll N$ as follows:

(3.2)\begin{equation} f_s(v,t) = \sum_i \varphi_i(v) f_i(t) = \sum_{i,j} \varphi_i(v) \, \mathbb{M}_{ij}^{{-}1} \sum \limits_\alpha w_\alpha \varphi_j (v_\alpha(t)) , \end{equation}

where $\{f_i(t)\}$ are the coefficients of the projected distribution function, $f_s$, expressed in the basis $\{ \varphi _i\}$, and $\mathbb {M}_{ij} = \int \varphi _i \varphi _j \,\textrm {d}v$ are the elements of the corresponding mass matrix $\mathbb {M}$. The projected representation of the distribution function, $f_s(v)$, will be used as an auxiliary representation for the evaluation of the collision operator where differentiability is required. This type of projection also offers the benefits of smoothing the solution for appropriately chosen basis functions $\{\varphi _j\}$. We note that (3.1) remains the primary representation of the distribution function in our method, and that (3.2) is only used in order to satisfy the differentiability requirements of the collision operator.

To construct the semi-discretisation of the conservative Lenard–Bernstein operator, we return to its form in (2.3). We will discretise this equation first, and then derive the discrete coefficients $A_1$ and $A_2$ in terms of the discrete momentum and energy densities.

To discretise the conservative collisional dynamics,

(3.3)\begin{equation} \frac{\partial f}{\partial t} = \nu\frac{\partial}{\partial v} \cdot \left(\frac{\partial f}{\partial v} + A_1 f + A_2 vf \right) , \end{equation}

we apply a deterministic particle method (see Chertock (2017) for a review). Typically, deterministic particle methods are formulated for first-order transport-type problems. They can, however, be adapted to the context of diffusion problems as shown by Degond & Mustieles (1990). Following this approach, we rewrite Equation (3.3):

(3.4a,b)\begin{equation} \frac{\partial f}{\partial t} = \frac{\partial}{\partial v} (a(v, f) f ), \quad a(v,f) = \nu \left( \frac{1}{f}\frac{\partial f}{\partial v} + A_1 + A_2 v \right) . \end{equation}

This equation is approximately solved in terms of the particle distribution function (3.1), where the particle velocities $\{v_\alpha \}$ satisfy the following ordinary differential equations (ODEs):

(3.5)\begin{equation} \dot{v}_\alpha = a(v_\alpha, f) = \nu \left(\frac{ 1}{f(v_\alpha)} \frac{\partial f}{\partial v}(v_\alpha) + A_1 + A_2 v_\alpha \right). \end{equation}

Let us note that the first term on the right-hand side is not well defined if the distribution function $f$ is replaced by its particle representation $f_p$. Therefore, we will use the projection shown in (3.2) instead and replace both instances of the distribution function $f$ with the projected distribution function $f_s$ to arrive at the following equation:

(3.6)\begin{equation} \dot{v}_\alpha = a(v_\alpha, f_s) = \nu \left( \frac{ 1}{f_s(v_\alpha)}\frac{\partial f_s}{\partial v}(v_\alpha) + A_1 + A_2 v_\alpha \right). \end{equation}

The final step in obtaining the semi-discrete system of equations is to compute the coefficients $A_1$ and $A_2$. In analogy to the continuous case of § 2, $A_1$ and $A_2$ are determined by requiring conservation of the discrete momentum and energy¹:

(3.7)\begin{gather} \frac{{\rm d}}{{\rm d}t} \sum \limits_\alpha w_\alpha v_\alpha = \nu \sum \limits_\alpha w_\alpha \left[ \frac{ 1}{f_s(v_\alpha)}\frac{\partial f_s}{\partial v}(v_\alpha) + A_1 + A_2 v_\alpha \right]= 0, \end{gather}

(3.8)\begin{gather}\frac{1}{2}\frac{{\rm d}}{{\rm d}t} \sum \limits_\alpha w_\alpha v_\alpha^2 = \nu \sum \limits_\alpha w_\alpha \left[ \frac{ v_\alpha}{f_s(v_\alpha)} \frac{\partial f_s}{\partial v}(v_\alpha) + A_1 v_\alpha + A_2 v_\alpha^2 \right]= 0. \end{gather}

Upon introduction of the discrete mass, momentum and energy densities,

(3.9a–c)\begin{equation} n_h(v_\alpha) = \sum \limits_\alpha w_\alpha, \quad n_hu_h(v_\alpha) = \sum \limits_\alpha w_\alpha v_\alpha , \quad n_h\varepsilon_h(v_\alpha) = \sum \limits_\alpha w_\alpha v_\alpha^2, \end{equation}

we obtain a linear system of equations which can be solved to find the discrete $A_1$, $A_2$:

(3.10)\begin{equation} \left.\begin{array}{c@{}} \displaystyle A_1n_h + A_2n_hu_h ={-}\sum \limits_\alpha w_\alpha \dfrac{f'_s(v_\alpha)}{f_s(v_\alpha)}, \\ \displaystyle A_1n_h u_h + A_2n_h\varepsilon_h ={-}\sum \limits_\alpha w_\alpha v_\alpha \dfrac{f'_s(v_\alpha)}{f_s(v_\alpha)}. \end{array}\right\}\end{equation}

The solution to this linear system is as follows:

(3.11)\begin{equation} \begin{array}{c@{}} \displaystyle A_1 = \dfrac{1}{n_h\varepsilon_h - n_hu_h^2} \sum \limits_\alpha w_\alpha(u_hv_\alpha - \varepsilon_h)\dfrac{f'_s(v_\alpha)}{f_s(v_\alpha)}, \\ \displaystyle A_2 = \dfrac{1}{n_h\varepsilon_h - n_hu_h^2} \sum \limits_\alpha w_\alpha (u_h - v_\alpha)\dfrac{f'_s(v_\alpha)}{f_s(v_\alpha)}. \end{array}\end{equation}

We note that these quantities implicitly depend on time through their dependence on the particle velocities $\{v_\alpha (t)\}$, and so must be recomputed at every time step.

We also note that in order for the projection to preserve the moments of the distribution function, it is sufficient that the functions $\{1, v, v^2\}$ are contained in $span \{\varphi _j\}$. This can be demonstrated by considering the example of conservation of mass, which requires

(3.12)\begin{equation} \int 1 f_p \,\text{d}{v} = \int 1 f_s \,\text{d}{v}. \end{equation}

To show that this holds, we begin from the condition used to construct the projection,

(3.13)\begin{equation} \int \varphi_j f_p \,\text{d}{v} = \int \varphi_j f_s \,\text{d}{v}, \quad \forall\, j \in 1, \ldots, M. \end{equation}

Let $1 \in span \{ \varphi _j \}$. Then, there exists a set of coefficients $\{c_j\}$ such that $1 = \sum _j c_j \varphi _j$. Multiplying both sides of (3.13) with the corresponding $c_j$ and summing over $j$, we have

(3.14)\begin{equation} \sum_j c_j \int \varphi_j f_p \,\text{d}{v} = \sum_j c_j \int \varphi_j f_s \,\text{d}{v}, \end{equation}

which by linearity of the integral becomes

(3.15)\begin{equation} \int \sum_j c_j \varphi_j f_p \,\text{d}{v} = \int \sum_j c_j \varphi_j f_s \,\text{d}{v}. \end{equation}

Since we have that $1 = \sum _j c_j \varphi _j$, we then have that

(3.16)\begin{equation} \int 1 f_p \,\text{d}{v} = \int 1 f_s \,\text{d}{v}, \end{equation}

which is the required condition for mass conservation. Momentum and energy conservation follow similarly from requiring the functions $v$ and $v^2$ to be in the span of basis $\{\varphi _j\}$. We note that the spline basis chosen in the numerical implementation, as detailed in § 5, satisfies this property since it is a cubic basis.

We remark that at this time, a proof of monotonic entropy dissipation for the semi-discrete Lenard–Bernstein operator remains elusive. We have, however, numerically demonstrated that our discretisation maintains this property in § 5.

4. Time discretisation

In this chapter, we describe the temporal discretisation of the system of ODEs (3.6) by the implicit midpoint scheme, which is a possible choice for a method that maintains the discrete conservation laws exactly (up to machine precision). Let us start by introducing some notation. Denote by $h$ a single time step, by $t_n = t_0 + nh$ the time after the $n$th time step, and by $v_\alpha ^n \approx v_\alpha (t_n)$ the corresponding particle velocity. Further, let $v_\alpha ^{n + 1/2} = (v_\alpha ^n + v_\alpha ^{n + 1})/2$ be the particle velocity at the midpoint and, following (3.2), denote the spline representation of the distribution function at the midpoint by

(4.1)\begin{equation} f_s^{n+1/2}(v) = \sum_{i,j} \varphi_i (v) \mathbb{M}^{{-}1}_{ij}\sum \limits_\alpha w_\alpha \varphi_j(v_\alpha^{n + 1/2}) . \end{equation}

With this, the implicit midpoint scheme for the system (3.6) can be written as follows:

(4.2)\begin{equation} v_\alpha^{n+1} = v_\alpha^n + h a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ). \end{equation}

In order to verify conservation of the total momentum, let us consider the particle momentum at the $(n+1)$th time step,

(4.3)\begin{align} \sum \limits_\alpha w_\alpha v_\alpha^{n+1} & = \sum \limits_\alpha w_\alpha [ v_\alpha^n + h a( v_\alpha^{n + 1/2}, f_s^{n+1/2})] \nonumber\\ & = \sum \limits_\alpha w_\alpha v_\alpha^n + h \sum \limits_\alpha w_\alpha a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) \nonumber\\ & = \sum \limits_\alpha w_\alpha v_\alpha^n + h \sum \limits_\alpha w_\alpha \left[ \frac{ 1}{f_s^{n+1/2}}\frac{\partial f_s^{n+1/2}}{\partial v}(v_\alpha^{n + 1/2}) + A_1 + A_2 v_\alpha^{n + 1/2} \right]. \end{align}

We observe that the sum in the second term of this equation is exactly given by the left-hand side of (3.7), evaluated at $v_\alpha ^{n+1/2}$, and so is equal to zero as (3.7) is satisfied for all times in the scheme by construction. Thus, the particle momentum is conserved exactly by the implicit midpoint scheme.

A similar result can be demonstrated for the total particle energy. At the $(n+1)$th time step, we have

(4.4)\begin{align} \sum \limits_\alpha w_\alpha ( v_\alpha^{n+1} )^2 & = \sum \limits_\alpha w_\alpha [ v_\alpha^n + h a ( v_\alpha^{n + 1/2}, f_s^{n+1/2} )]^2 \nonumber\\ & = \sum \limits_\alpha w_\alpha (v_\alpha^{n})^2 + h^2 \sum \limits_\alpha w_\alpha a( v_\alpha^{n + 1/2}, f_s^{n+1/2} )^2 \nonumber\\ & \quad + 2h \sum \limits_\alpha w_\alpha v_\alpha^n a( v_\alpha^{n + 1/2}, f_s^{n+1/2}). \end{align}

Now, we split the last term and replace one $v_\alpha ^n$ by using the implicit midpoint rule, i.e. $v_\alpha ^n = v_\alpha ^{n+1} - ha( v_\alpha ^{n + 1/2}, f_s^{n+1/2} )$ to obtain the following:

(4.5)\begin{align} & 2h \sum \limits_\alpha w_\alpha v_\alpha^n a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) \nonumber\\ & \quad = h \sum \limits_\alpha w_\alpha v_\alpha^n a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) + h \sum \limits_\alpha w_\alpha v_\alpha^n a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) \nonumber\\ & \quad = h \sum \limits_\alpha w_\alpha v_\alpha^n a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) \nonumber\\ & \qquad + h\sum \limits_\alpha w_\alpha [ v_\alpha^{n+1} - ha( v_\alpha^{n + 1/2}, f_s^{n+1/2} )]a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) \nonumber\\ & \quad = h\sum \limits_\alpha w_\alpha (v_\alpha^n + v_\alpha^{n+1}) a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) - h^2 \sum \limits_\alpha w_\alpha a( v_\alpha^{n + 1/2}, f_s^{n+1/2})^2. \end{align}

The last term in this equation cancels the second term on the right-hand side of (4.4), and we are left with

(4.6)\begin{equation} \sum \limits_\alpha w_\alpha ( v_\alpha^{n+1} )^2 = \sum \limits_\alpha w_\alpha (v_\alpha^{n})^2 + 2h\sum \limits_\alpha w_\alpha v_\alpha^{n + 1/2}a( v_\alpha^{n + 1/2}, f_s^{n+1/2} ) . \end{equation}

The second term on the right-hand side of (4.6) equals the expression in (3.8) evaluated at the midpoint $v_\alpha ^{n + 1/2}$, and therefore is zero. Thus, the particle energy is preserved exactly in time.

5. Numerical results

In this chapter, we present numerical experiments for the one-dimensional conservative Lenard–Bernstein operator using B-spline basis functions of arbitrary order for projecting the particle distribution function. The implementation is based on the Julia programming language (Bezanson et al. 2017), and the package is available publicly (Kraus, Blickhan & Jeyakumar 2023).

5.1. Convergence of the semi-discrete operator

We demonstrate convergence properties of the semi-discretisation under projection onto a B-spline basis and particle sampling. First, we verify that the semi-discretisation indeed preserves the Maxwellian as an equilibrium solution under projection. To this end, we compute the projection of an exact Maxwellian onto the spline basis as follows:

(5.1)\begin{equation} \int \sum_j f_j \varphi_j(v) \varphi_i(v) \,\text{d}{v} = \int f_M(v) \varphi_i(v) \,\text{d}{v} , \end{equation}

where ${f_j}$ are the coefficients of the spline for which we are solving, and $f_M(v) = 1/\sqrt {2{\rm \pi} } \exp (-v^2/2)$ is the Maxwellian of mean $\mu = 0$ and variance $\sigma ^2 = 1$. Rearranging this expression, we obtain the following spline coefficients for the projected Maxwellian:

(5.2)\begin{equation} f_j = \mathbb{M}_{ij}^{{-}1} \int f_M(v) \varphi_i(v) \,\text{d}{v}, \end{equation}

where $\mathbb {M}_{ij} = \int \varphi _i \varphi _j \,\textrm {d}{v}$ are the elements of the mass matrix $\mathbb {M}$ as before. The spline-projected representation of the Maxwellian distribution $f_M(v)$ is then given by $f_{s,M}(v) = \sum _j f_j \varphi _j(v)$ with the coefficients $f_j$ from (5.2). We use this projected Maxwellian to compute the right-hand side of (3.6), as follows:

(5.3)\begin{equation} \dot{v}_\alpha = \nu \left( \frac{ 1}{f_{s,M}(v_\alpha)}\frac{\partial f_{s,M}}{\partial v}(v_\alpha) + A_1 + A_2 v_\alpha \right), \end{equation}

where the coefficients $A_1$ and $A_2$ are computed using the projected Maxwellian distribution $f_{s,M}$ in (3.11). The $L^2$ norm of the time derivative of the particle velocities, $\| \dot {v}_\alpha \|_2$, can then be used to check if the Maxwellian is an equilibrium solution under projection, as this norm should approach zero with increasing spline resolution in such case. We compute this quantity for a range of spline resolutions, using a sample of $N = 100\,000$ particles from a normal distribution where the sample is strictly used for evaluation of (5.3) and never for projection. Figure 1 shows the convergence of $\| \dot {v}_\alpha \|_2$ with an increasing number of splines, computed using cubic spline basis functions. As expected, the norm of the time derivative approaches zero with an increasing number of splines at a rate which corresponds to the order of the splines used (cubic splines have order $k = 4$).

Figure 1. Convergence of the $L^2$ norm of the particle velocity gradient, $\| \dot {v}_\alpha \|_2$, when computed using a true Maxwellian $f$, against the number of splines. The dashed line shows the reference curve $y = x^{-4}$.

Secondly, we demonstrate convergence of the semi-discretisation under particle sampling. Here, we instead compute the sample variance of the particle velocity time derivatives, i.e. $\sum \dot {v}_\alpha ^2/N$, keeping the spline resolution fixed and varying the number of particles in the sample. Here, we directly project the particles to compute the spline representation of $f$, as per (3.2) and (3.6). The results of this are shown in figure 2, and we observe that the sample variance converges at a rate slightly above $1/N$, which is the expected convergence rate for the sample variance generated by statistical sampling methods.

Figure 2. Convergence of the sample variance of the particle velocity time derivative, $\sum \dot {v}_\alpha ^2/N$, with the number of particles, shown on a logarithmic scale. The dashed line represents the reference curve $y = 1/N$.

These two results demonstrate that the Maxwellian distribution remains an equilibrium solution under semi-discretisation, and the method demonstrates the expected convergence properties with both number of splines and particles.

5.2. Relaxation of a shifted normal distribution

In the first example, we initialise the distribution function with a standard normal distribution shifted to the right by $\mu = 2$, obtaining a distribution with mean $\mu = 2$ and variance $\sigma ^2 = 1$ as shown in the following equation:

(5.4)\begin{equation} f(v,t=0) = \frac{1}{\sqrt{2{\rm \pi}}}\,{\rm e}^{- ({1}/{2})(v - 2)^2}. \end{equation}

The particles are then initialised by independently and identically sampling from this distribution, for $N = 1000$ particles. The spline distribution is initialised by $L^2$ projecting the initial particle distribution onto the spline basis, with the spline coefficients computed as per (3.2). A cubic-spline basis of 41 elements was used, with equally spaced knots on the velocity domain $v \in [-10, 10]$. The time integration is performed using the implicit midpoint scheme, which is of second order. A time step of $h = 8 \times 10^{-4}$ and a collision frequency of $\nu = 1$ was used. The simulation was run until a final normalised time of $t = 1$, corresponding to $1250$ time steps.

In the initial condition, shown in figure 3(a), we observe slight variations in the spline-projected distribution function due to the sampling error of the particles. In the final distribution, we observe the expected behaviour of the projected distribution approaching an exact normal with the initial variations smoothed out, and due to the momentum conservation, the mean of the distribution stays at the same value ($v = 2$). Here, the particle momentum and the energy are conserved up to machine precision. The evolution of the energy and momentum error are shown in figure 4. The evolution of the entropy, $S = \int f_s \log {f_s} \,\textrm {d}v$, is shown in figure 5 (normalised by the magnitude of its initial value), and we observe the monotonic dissipation of the entropy over the course of the simulation as expected.

Figure 3. The (a) initial and (b) final distributions when the initial condition is chosen to be a normal distribution of mean $\mu = 2$ and variance $\sigma ^2 = 1$, for a sample of $N=1000$ particles.

Figure 4. Energy and momentum conservation during the simulation, for the initial condition of a shifted normal distribution.

Figure 5. Evolution of the normalised entropy, i.e. $S/|S(t=0)|$ where $S = \int f_s \log {f_s} \,\textrm {d}v$, during the simulation.

We also observe that the method works well even for small numbers of particles. Results for the same simulation using a sample of $N=200$ particles instead are shown in figure 6. The particle energy and momentum are again conserved up to machine precision in relative error, with similar behaviour as shown in figure 4.

Figure 6. The (a) initial and (b) final distributions when the initial condition is chosen to be a normal distribution of mean $\mu = 2$ and variance $\sigma ^2 = 1$, for a sample of $N = 200$ particles.

5.3. Relaxation of a bi-Maxwellian distribution

Next, we consider a double Maxwellian distribution as the initial condition. Each peak is a standard normal distribution which has been shifted from the origin by $v = \pm 2$ as per the following equation:

(5.5)\begin{equation} f(v,t=0) = \frac{1}{\sqrt{2{\rm \pi}}} ( {\rm e}^{- ({1}/{2})(v - 2)^2} + {\rm e}^{- ({1}/{2})(v + 2)^2}), \end{equation}

and the sampled distribution is shown in figure 7(a). The sample for the particle distribution function is again of size $N = 1000$ particles. We use the identical numerical set-up as the previous example, in particular the same time step of $h = 8 \times 10^{-4}$ and the same collision frequency of $\nu = 1$. The final distribution obtained after time integration until $t = 5$ is shown in figure 7(b). Again, the equilibrated distribution is a Gaussian with equal mean and variance to the initial condition. The particle energy is conserved up to machine precision in relative error. The particle momentum is conserved up to a relative error of the order of $10^{-14}$. The behaviour of the two quantities over time is oscillatory and similar to figure 4. The normalised entropy also decreases monotonically in time, and this is illustrated in figure 8.

Figure 7. The (a) initial and (b) final distribution functions when the initial condition is chosen to be a bi-Maxwellian, in both particle and spline bases, for $N=1000$ particles.

Figure 8. Evolution of the normalised entropy over the simulation, where the initial condition is a bi-Maxwellian.

5.4. Relaxation of a uniform distribution

In the last example, we initialise the distribution function with a uniform distribution shifted and scaled to the interval $v \in [-2, 2]$, i.e.

(5.6)\begin{equation} f(v, t=0) = \begin{cases} \frac{1}{4}, & \text{if } v \in [{-}2, 2],\\ 0, & \text{else.} \end{cases} \end{equation}

A sample of $N = 200$ particles is taken from this distribution. A time step of $h = 1 \times 10^{-4}$ is used and the final integration time is $t = 1$. All other parameters are kept the same. Figure 9 shows the initial and final distributions obtained in the simulation, demonstrating again the expected result that the initial distribution relaxes to a normal distribution. The resultant particle distribution function retains the same mean up to a relative error of the order of $10^{-14}$ (which is equivalent to momentum being preserved at this level). The energy is preserved to a relative error of the order of $10^{-15}$. The entropy decreases monotonically as expected and is illustrated in figure 10. We note that the method performs as well here as it does in the other examples despite this being a more challenging case, as the uniform distribution is discontinuous and therefore not amenable to being represented using B-spline basis functions.

Figure 9. The (a) initial and (a) final distribution functions for the case of a uniform initial condition. A sample of $N=200$ particles is used.

Figure 10. Evolution of the normalised entropy over the simulation, where the initial condition is a uniform distribution.

5.5. Remarks

It is important to ensure that the chosen velocity domain for a simulation is sufficiently large such that no particle leaves this domain at any time. There is no sensible method for returning a particle to the simulation domain, as the true velocity space domain for this problem is infinite. Practically, a particle leaving the domain will return zero when the spline projected distribution is evaluated on the particle's velocity, which leads to the evaluation of (3.6) becoming undefined due to division by zero.

6. Conclusion

In this work, we have outlined the development of structure-preserving particle-based algorithms for the simulation of collision operators. While the approach itself is general, it has been specifically applied to a conservative version of the Lenard–Bernstein operator. We have derived an energy- and momentum-preserving particle discretisation for this operator, and the implicit midpoint method is shown to exactly preserve these quantities in time as well. We have demonstrated the convergence properties of the semi-discretisation under the projection used, as well as under particle sampling. Numerical examples for the one-dimensional case demonstrated the viability of the method and verified its conservation properties. The method is implemented in the Julia programming language with the respective repository available (Kraus et al. 2023). The proposed method can be coupled to any Vlasov–Poisson or Vlasov–Maxwell particle solver, and future research will detail such a coupling and its benefits.

Currently, we are adapting the approach presented here for the Landau collision operator using the metriplectic formulation, and the results will be reported in a follow-up paper.