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 Reference Hairer and Wanner2006). 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 (Reference Morrison2017) 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 (Reference Chen, Chacón and Barnes2011), Markidis & Lapenta (Reference Markidis and Lapenta2011), Squire, Qin & Tang (Reference Squire, Qin and Tang2012), Evstatiev & Shadwick (Reference Evstatiev and Shadwick2013), Qin et al. (Reference Qin, Liu, Xiao, Zhang, He, Wang, Sun, Burby, Ellison and Zhou2016), Burby (Reference Burby2017), Kraus et al. (Reference Kraus, Kormann, Morrison and Sonnendrücker2017), Zhang & Gamba (Reference Zhang and Gamba2017), Campos Pinto, Kormann & Sonnendrücker (Reference Campos Pinto, Kormann and Sonnendrücker2022).
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 Reference Yoon and Chang2014; Taitano et al. Reference Taitano, Chacón, Simakov and Molvig2015; Hirvijoki & Adams Reference Hirvijoki and Adams2017; Kraus & Hirvijoki Reference Kraus and Hirvijoki2017; Shiroto & Sentoku Reference Shiroto and Sentoku2019), structure-preserving discretisations for collision operators with particles have been considered more recently. Hirvijoki, Kraus & Burby (Reference Hirvijoki, Kraus and Burby2018) considered an approach where the weights of the marker particles are varied, instead of their velocities. Carrillo et al. (Reference Carrillo, Hu, Wang and Wu2020) and Hirvijoki (Reference Hirvijoki2021) used finite-sized marker particles to discretise the Landau operator. Mollén et al. (Reference Mollén, Adams, Knepley, Hager and Chang2021) and Pusztay, Knepley & Adams (Reference Pusztay, Knepley and Adams2022) focused on projection/interpolation techniques for computing collision operators for particles. An alternative approach is that of Tyranowski (Reference Tyranowski2021), 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 (Reference Lenard and Bernstein1958). We employ the deterministic particle method (Degond & Mustieles Reference Degond and Mustieles1990; Chertock Reference Chertock2017) 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 Reference Lenard and Bernstein1958) is a scalar operator of the form
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:
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 (Reference Kraus2013) (see also Filbet & Sonnendrücker Reference Filbet and Sonnendrücker2003) 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 (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:
Integrating this by parts, we obtain the following condition:
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:
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] = \int v^m f \,\text {d}{v}$
, 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 (2.9) is
where $nu$
and $n\varepsilon$
are the momentum and energy density, respectively, and are related to the moments as follows:
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:
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 (Reference Filbet and Sonnendrücker2003) and is closely related to the operator studied in Hakim et al. (Reference Hakim, Francisquez, Juno and Hammett2020).
2.2. The H-theorem
We follow a similar strategy to Hakim et al. (Reference Hakim, Francisquez, Juno and Hammett2020) to demonstrate that the continuous operator satisfies the H-theorem. Let us denote the collisional flux by
so that
The change in entropy is then
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
Hence, the entropy is monotonically dissipated (provided that $f$
is non-negative):
The entropy is stationary when $F[f_\textrm {eq}]=0$
, that is when $f_\textrm {eq}$
solves the partial differential equation (PDE)
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$
,
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. Reference Carrillo, Hu, Wang and Wu2020; Hirvijoki Reference Hirvijoki2021). 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_\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:
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,
we apply a deterministic particle method (see Chertock (Reference Chertock2017) 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 (Reference Degond and Mustieles1990). Following this approach, we rewrite Equation (3.3):
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):
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:
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 energyFootnote 1:
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_\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
To show that this holds, we begin from the condition used to construct the projection,
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
which by linearity of the integral becomes
Since we have that $1 = \sum _j c_j \varphi _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 $\{\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
With this, the implicit midpoint scheme for the system (3.6) 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 time step,
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
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:
The last term in this equation cancels the second term on the right-hand side of (4.4), and we are left with
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. Reference Bezanson, Edelman, Karpinski and Shah2017), and the package is available publicly (Kraus, Blickhan & Jeyakumar Reference Kraus, Blickhan and Jeyakumar2023).
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:
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:
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:
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:
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:
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.
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. Reference Kraus, Blickhan and Jeyakumar2023). 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.
Acknowledgements
Editor P. Helander thanks the referees for their advice in evaluating this article.
Declaration of interest
The authors report no conflict of interest.
Appendix 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 \sim N(\mu, \sigma ^2)$
is given by
where the cumulants are the coefficients of the Taylor expansion in $s$
, $\kappa _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.Footnote 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 $\kappa _3(X)$
and $\kappa _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,
With this, moments of the distribution function are the Taylor coefficients of the moment-generating function $M_k(t)=\partial _s^{(k)}M(0;t)=\mathbb {E}(X_t^k)$
. If we write the Lenard–Bernstein collision operator as
then, after integrating (and neglecting boundary terms originating from integration by parts), the relaxation equation $\partial _t f=C[f]$
becomes a PDE for $M(s;t)$
, namely
where the substitution $\int \textrm {e}^{sv}v^k f \, \textrm {d} v=\partial _s^{(k)}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.
(i) For arbitrary $p$
, the steady-state moment-generating function $M_\infty (s)=\lim _{t\to \infty } M(s;t)$ can be inferred from the ODE
(A7)\begin{equation} \sum_{k=0}^p A_k M _{\infty}^{(k)} = s M_{\infty} . \end{equation}In particular if $p=1$
, we have a first-order ODE
(A8)\begin{equation} M'_\infty = \frac{s-A_0}{A_1} M_\infty \iff K_\infty'=\frac{s-A_0}{A_1} , \end{equation}together with the requirement that $K_\infty (0)=\ln 1 =0$
. The solution is the cumulant-generating function of a Gaussian with mean $\mu =-A_0/A_1$ and variance $\sigma ^2=1/A_1$:
(A9)\begin{equation} K_\infty(s) = \mu s + \tfrac{1}{2}\sigma ^2 s^2. \end{equation}
(ii) In the case where $p=1$
, normalising time to the collision frequency multiplied by the variance $t\mapsto t/(\sigma ^2\nu )$, the relaxation equation in terms of the cumulant-generating function is
(A10)\begin{equation} \partial_t K + s\partial_s K = \mu s + \sigma^2s^2 , \end{equation}which is a linear non-homogeneous first-order PDE. We apply the method of characteristics. Let the curve $s(\tau )$
and $t(\tau )$ and $\kappa (\tau )=K(s(\tau ),t(\tau ))$ be such that
(A11)\begin{equation} \left.\begin{array}{c@{}} \displaystyle \dfrac{{\rm d} t}{{\rm d} \tau} = 1, \\ \displaystyle \dfrac{{\rm d} s}{{\rm d} \tau} = s(\tau),\\ \displaystyle \dfrac{{\rm d} \kappa}{{\rm d} \tau} = \dfrac{{\rm d} }{{\rm d} \tau} K(s(\tau),t(\tau)) = \partial_s K \dfrac{{\rm d}s}{{\rm d} \tau} + \partial_t K\dfrac{{\rm d}t}{{\rm d} \tau} = \mu s(\tau) + \sigma^2 s(\tau)^2 , \end{array}\right\} \end{equation}with initial conditions $s(0)=s_0$
, $t(0)=0$ and $\kappa (0)=K(s_0,0)=\mu s_0 + \frac {1}{2}\sigma ^2 s_0^2+ R(s_0)$. The solutions are
(A12)\begin{equation} \left.\begin{array}{c@{}} t(\tau)=\tau,\\ s(\tau;s_0)=s_0 \, {\rm e}^\tau,\\ \kappa(\tau;s_0) = \mu s_0 \, {\rm e}^\tau + \tfrac{1}{2}\sigma^2 s_0^2 \,{\rm e}^{2\tau}+R(s_0). \end{array}\right\} \end{equation}Inverting $s(\tau ;s_0)$
and $t(\tau ;s_0)$, we obtain the time-evolution of the cumulant-generating function as
(A13)\begin{equation} \left.\begin{array}{c@{}} \tau = t , \\ s_0(s,t) = s \,{\rm e}^{{-}t} , \\ K(s,t) = \mu s + \tfrac{1}{2} \sigma^2 s^2 + R(s\, {\rm e}^{{-}t}) = K_\infty(s) + R(s\, {\rm e}^{{-}t}). \end{array}\right\} \end{equation}As time tends to infinity the initial ‘excess’ cumulants are lost,
(A14)\begin{equation} \lim_{t\to\infty}R(s\, {\rm e}^{{-}t})=R(0)=0 . \end{equation}We are now in a position to determine the time-evolution of the individual cumulants by differentiating with respect to $s$
,
(A15)\begin{equation} \left.\begin{array}{c@{}} \partial_sK(s,t)=\mu+\sigma^2s+R'(s\, {\rm e}^{{-}t})\, {\rm e}^{{-}t} , \\ \partial_s^2 K(s,t)=\sigma^2 + R''(s\, {\rm e}^{{-}t})\, {\rm e}^{{-}2t} , \\ \partial_s^{(k)}K(s,t) =R^{(k)}(s\, {\rm e}^{{-}t})\, {\rm e}^{{-}kt}, \quad k>2 . \end{array}\right\} \end{equation}By evaluating at $s=0$
, we have
(A16)\begin{equation} \left.\begin{array}{c@{}} K_1(t) = \mu + R_{1} \, {\rm e}^{{-}t} , \\ K_2(t) = \sigma^2 + R_{2} \, {\rm e}^{{-}2t} , \\ K_k(t) = R_{k}\, {\rm e}^{{-}kt}, \quad k>2 , \end{array}\right\} \end{equation}where $R_1=\partial _s K(0,0)-\mu$
, $R_2=\partial _s^2K(0,0) -\sigma ^2$ and $R_{k}=\partial _s^{(k)}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 (3.6), instead of solving the system of equations shown in (3.10). 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 time step of $\Delta t = 5 \times 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 11. Decay of the fourth- and fifth-order cumulants, $\kappa _4$
and $\kappa _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.
Open access
