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

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.


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 vivid 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 focused on dissipation-less 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.
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 et al. (2018) consider an approach where the weights of the marker particles are varied, instead of their velocities.Carrillo et al. (2020) and Hirvijoki (2021) use finite-sized marker particles to discretise the Landau operator.Mollén et al. (2021) and Pusztay et al. (2022) focus 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 (Chertock (2017); Degond & Mustieles (1990)) 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 Section 2, we detail the derivation of the conservative Lenard-Bernstein operator.In Section 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.

The conservative Lenard-Bernstein operator
The Lenard-Bernstein collision operator (Lenard & Bernstein (1958)) is a scalar operator of the form where f : R d × [0, ∞) → R is the single-particle distribution function, v ∈ R d is the velocity and ν 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: 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) for this operator is an d-dimensional Gaussian distribution.

Construction of the conservative operator
The Lenard-Bernstein operator (1) preserves mass density, the zeroth moment of the distribution function.However, it does not preserve momentum density nor 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, 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 Equation (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: Integrating this by parts, we obtain the following condition: Without loss of generality, we assume that f and ∂f /∂v approach zero as |v| → ∞, so that the first term in Equation ( 5) is zero and we obtain the following condition: for k = 1, 2. Integrating the first term by parts once again, this equation becomes: where the assumption that f approaches zero as |v| tends to infinity has been utilised once more.Writing the moments as M m [f ] = v m f dv, we obtain the following conditions: These conditions provide a linear system of equations that can be solved for the coefficients A 1 , A 2 : The solution to the system of equations in (9) is: where nu and nε are the momentum and energy density, respectively, and are related to the moments as follows: Let us note that here, n, u, and ε 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 (3), we obtain the following operator: which can be seen as a conservative version of the Lenard-Bernstein operator (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).

H-theorem
We follow a similar strategy to Hakim et al. (2020) to demonstrate that the continuous operator satisfies a H-theorem.Let us denote the collisional flux by The change in entropy is then where integration by parts, the assumption that F goes to zero as |v| → ∞ and the substitution ∂f /∂v = F − A 1 f − A 2 vf were used.By Equation ( 6), the second term in the last expression is designed to vanish, namely Hence, the entropy is monotonically dissipated (provided that f is non-negative): The entropy is stationary when F [f eq ] = 0, that is when f eq solves the PDE The (unique) solution with conditions Thus, the continuous operator of ( 13) satisfies the H-theorem.

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 particlerepresentation 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 where {v α (t)} N α=1 are the particle velocities which evolve over time, where 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 {φ j } M j=1 for M ≪ N as follows: where {f i (t)} are the coefficients of the projected distribution function, f s , expressed in the basis {φ i }, and M ij = φ i φ j dx are the elements of the corresponding mass matrix 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 {φ j }.We note that Equation ( 20) remains the primary representation of the distribution function in our method, and that Equation ( 21) 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 Equation ( 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, 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 ( 22): This equation is approximately solved in terms of the particle distribution function (20), where the particle velocities {v α } satisfy the following ordinary differential equations: 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 Equation ( 21) instead and replace both instances of the distribution function f with the projected distribution function f s to arrive at the following 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 Section 2, A 1 and A 2 are determined by requiring conservation of the discrete momentum and energy1 : Upon introduction of the discrete mass, momentum, and energy densities, we obtain a linear system of equations which can be solved to find the discrete A 1 , A 2 : . ( The solution to this linear system is as follows: We note that these quantities implicitly depend on time through their dependence on the particle velocities {v α (t)}, and so must be recomputed at every timestep.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{φ j }.This can be demonstrated by considering the example of conservation of mass, which requires: To show this holds, we begin from the condition used to construct the projection: Let 1 ∈ span{φ j }.Then, there exists a set of coefficients {c j } such that 1 = j c j φ j .
Multiplying both sides of (32) with the corresponding c j and summing over j, we have: which by linearity of the integral becomes: Since we have that 1 = j c j φ j , we then have that 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 {φ j }.We note that the spline basis chosen in the numerical implementation, as detailed in Section 5, satisfies this property since it is a cubic basis.
We remark that at this time, a proof of monotonic entropy dissipation for the semidiscrete Lenard-Bernstein operator remains elusive.We have, however, numerically demonstrated that our discretisation maintains this property in Section 5.

Time discretisation
In this chapter, we describe the temporal discretisation of the system of ODEs ( 25) 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 timestep, by t n = t 0 + nh the time after the n-th timestep, and by v n α ≈ v α (t n ) the corresponding particle velocity.Further, let v n+1/2 α = (v n α +v n+1 α )/2 be the particle velocity at the midpoint and, following Equation ( 21), denote the spline representation of the distribution function at the midpoint by With this, the implicit midpoint scheme for the system (25) can be written as follows: In order to verify conservation of the total momentum, let us consider the particle momentum at the n + 1-th timestep: We observe that the sum in the second term of this equation is exactly given by the lefthand side of Equation ( 26), evaluated at v n+1/2 α , and so is equal to zero as ( 26) 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 timestep, we have: Now, we split the last term and replace one v n α by using the implicit midpoint rule, i.e.
to obtain the following: The last term in this equation cancels the second term on the right-hand side of ( 39), and we are left with: The second term on the right-hand side of (41) equals the expression in (27) evaluated at the midpoint v n+1/2 α , and therefore is zero.Thus, the particle energy is preserved exactly in time.

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 et al., 2023).

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: Figure 1: Convergence of the L 2 norm of the particle velocity gradient, ∥ vα ∥ 2 , when computed using a true Maxwellian f , against the number of splines.The dashed line shows the reference curve y = x −4 .
where f j are the coefficients of the spline for which we are solving, and ) is the Maxwellian of mean µ = 0 and variance σ 2 = 1.Rearranging this expression, we obtain the following spline coefficients for the projected Maxwellian: where M ij = φ i φ j dv are the elements of the mass matrix M as before.The splineprojected representation of the Maxwellian distribution f M (v) is then given by f s,M (v) = j f j φ j (v) with the coefficients f j from (43).We use this projected Maxwellian to compute the right-hand side of Equation ( 25), as follows where the coefficients A 1 and A 2 are computed using the projected Maxwellian distribution f s,M in (30).The L 2 norm of the time derivative of the particle velocities, ∥ vα ∥ 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 Equation ( 44) and never for projection.Figure 1 shows the convergence of ∥ vα ∥ 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).
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. v2 α /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 Equations ( 21) and ( 25).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.
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.

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 µ = 2, obtaining a distribution with mean µ = 2 and variance σ 2 = 1 as shown in the following equation: The particles are then initialised by independently and identically (iid) 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 Equation ( 21).A cubic-spline basis of 41 elements was used, with equally spaced knots on the velocity domain v ∈ [−10, 10].The time integration is performed using the implicit midpoint scheme, which is of second order.A time-step of h = 8 × 10 −4 and a collision frequency of ν = 1 was used.The simulation was run until a final normalised time of t = 1, corresponding to 1250 timesteps.
In the initial condition, shown on the left of Figure 3, 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).
Figure 3: The initial and final distributions when the initial condition is chosen to be a normal distribution of mean µ = 2 and variance σ 2 = 1, for a sample of N = 1000 particles.
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 = f s log f s dv, 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.
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.

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 = ±2 as per the following equation: and the sampled distribution is shown in Figure 7 on the left.The sample for the particle distribution function is again of size N = 1000 particles.We use the identical numerical setup as the previous example, in particular the same time-step of h = 8 × 10 −4 and the same collision frequency of ν = 1.The final distribution obtained after time integration Figure 7: The initial and final distribution functions when the initial condition is chosen to be a bi-Maxwellian, in both particle and spline bases, for N = 1000 particles.
until t = 5 is shown in Figure 7 on the right.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 on 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.

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 ∈ [−2, 2], i.e.: A sample of N = 200 particles is taken from this distribution.A timestep of h = 1 × 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 on the order of 10 −14 (which is equivalent to momentum being preserved at this level).The energy is preserved to a relative error on 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 Figure 8: Evolution of the normalised entropy over the simulation, where the initial condition is a bi-Maxwellian.
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.

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 (25) becoming undefined due to division by zero.

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 language with the respective repository available at (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

A Time-evolution of cumulants
One useful fact about the steady-state solution of the Lenard-Bernstein operator being a normal distribution is the fact that its cumulants of order three and above are all zero.The cumulant is a closely related quantity to a moment, being defined through a cumulant generating function which is obtained by taking the natural logarithm of the moment generating function of the distribution.The cumulant generating function for a normally-distributed random variable X ∼ N (µ, σ2 ) is given by where the cumulants are the coefficients of the Taylor expansion in s, κ n = K (n) (0).In this instance, the cumulant generating function has no terms at order three and above, implying that the corresponding cumulants of the normal distribution are zero.We can also see that the first and second cumulants are simply the mean and variance, respectively 2 .For ease of computation, the third and fourth cumulants can be related to central moments of the random variable (those centred around the mean) through the following relations: where κ 3 (X) and κ 4 (X) are the third and fourth cumulants, respectively.In the discrete setting, the behaviour of the discretised cumulants can act as a quantitative check of how close the solution is to the known solution.
In fact, the time-evolution of the cumulants can be solved analytically in the case where the coefficients A k of the Lenard-Bernstein collision operator are held fixed.Let the moment-generating function be the Wick-rotation of the Fourier transform of the distribution: (s;t) .
With this, moments of the distribution function are the Taylor coefficients of the moment- If we write the Lenard-Bernstein collision operator as then, after integrating (and neglecting boundary terms originating from integration by parts), the relaxation equation where the substitution e sv v k f dv = ∂ (k) s M was exploited.By setting s = 0 we see that the zeroth-moment is conserved, which can be attributed to the collision operator being a divergence.The construction above is fairly general in the sense that the number of terms in the collision operator is capped by an arbitrary p.There are two distinct questions to address analytically.
As time tends to infinity the initial "excess" cumulants are lost: We are now in a position to determine the time-evolution of the individual cumulants by differentiating with respect to s: ∂ s K(s, t) = µ + σ 2 s + R ′ (se −t )e −t , ∂ 2 s K(s, t) = σ 2 + R ′′ (se −t )e −2t , ∂ (k)  s K(s, t) = R (k) (se −t )e −kt , k > 2.
By evaluating at s = 0, we have where R 1 = ∂ s K(0, 0) − µ, R 2 = ∂ 2 s K(0, 0) − σ 2 and R k = ∂ (k) s K(0, 0) = R (k) (0) for k > 2 are the initial residual cumulants.This shows that the decay rate of the k th cumulant is proportional to k, namely that the decay rate scales linearly with the order of the cumulant.
We check this scaling numerically in our method by fixing the A 1 and A 2 coefficients in equation ( 25), instead of solving the system of equations shown in equations ( 29).Initialising the simulation with a double Maxwellian distribution, as shown in the second example, for N = 1000 particles, M = 41 splines of order 4, and a timestep of ∆t = 5 × 10 −3 , we fix the coefficients to be A 0 = 0 and A 1 = 1.The evolution of the higher order cumulants is shown in Figure 11.We observe that the cumulants scale at the predicted rate, until they reach the level of accuracy supported by the chosen resolution (approximately 10 −2 for a particle resolution of N = 1000.The saturation point of the cumulants corresponds to the solution reaching equilibrium.

Figure 2 :
Figure 2: Convergence of the sample variance of the particle velocity time derivative, v2 α /N , with the number of particles, shown on a logarithmic scale.The dashed line represents the reference curve y = 1/N .

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

Figure 9 :
Figure 9: The initial and final distribution functions for the case of a uniform initial condition, with the initial condition shown on the left and the final result shown on the right.A sample of N = 200 particles is used.

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

Figure 11 :
Figure 11: Decay of the fourth-and fifth-order cumulants, κ 4 and κ 5 , over the course of a simulation where the initial distribution is taken as a sample of N = 1000 particles from a double Maxwellian distribution.The cumulants are normalised to their initial value and analytic scaling laws are shown as dashed lines.