A structure-preserving collisional particle method for the Landau kinetic equation

The Landau equation is a fundamental kinetic equation that describes the evolution of the distribution of charged particles in a collisional plasma where grazing collisions are predominant [23, 30]. In the space-homogeneous case, the Landau equation governs the evolution of the density $f(t,v)$ of particles with velocity $v$ at time $t$:

The collision kernel $A$ is given by

where $\Lambda>0$ is the collision strength, $I_{d}$ is the identity matrix, and the parameter $\gamma$ can take values within the range $(-d-1,1]$. The most important case is when $d=3$ and $\gamma=-3$ that corresponds to the original Landau equation for charged particles interacting with Coulomb potentials, derived from the Boltzmann equation in the grazing collision limit.

The well-posedness of the Landau equation has been extensively studied in the literature (see e.g. [29, 12, 17, 19] and references therein), covering various aspects of existence, uniqueness, and regularity under different conditions. Recently, Guillen and Silvestre [20] achieved a significant breakthrough by proving the global well-posedness of the Landau equation with smooth initial data. Their paper also provides a detailed overview of the historical developments and key references in this area.

Various numerical methods have been developed for solving the Landau equation, including finite difference schemes [31], the Fourier–Galerkin spectral method [27], the direct simulation Monte Carlo method [2, 5, 32, 24], and particle methods [16, 8, 9]. For comprehensive reviews on this topic, we refer readers to [6, 13, 8, 4]. Among these, particle methods have seen significant advancements over the past few decades. These methods approximate the solution as a linear combination of Dirac $\delta$-functions centered at particle locations.

Particle methods are categorized into deterministic and stochastic approaches. In deterministic particle methods, particle locations and weights evolve based on ODE systems derived from the weak formulation of the target equation. Carrillo et al. [8] proposed a gradient-flow-based deterministic particle method that preserves critical physical properties, such as conservation of mass, momentum, energy, and entropy dissipation. The evaluation of the collision operator $Q(f,f)$ in this method typically requires $O(N^{2})$ operations, where $N$ is the number of velocity points. With the help of treecode summation, it can be reduced to $O(N\log N)$. More recently, the random batch particle method [9] further improved efficiency, achieving $O(N)$ complexity. However, these are primarily numerical approximations, which effectively recover density but do not fully capture individual particle trajectories.

A typical stochastic particle approximation for the Landau equation stems from its mean-field SDE formulation, which is thought to capture the inherent randomness of particle trajectories. Formally, the Landau equation corresponds to the McKean–Vlasov SDE

where $K:=\nabla\cdot A=(1-d)\Lambda|z|^{\gamma}z$, and $f$ represents both the density function of $V$ and the solution to the Landau equation. A direct particle approximation of this SDE [16] takes the form

Convergence results for this system are established in specific cases, as shown in [18, 17, 7] in some cases. The computational cost of simulating the particle system (1.3) $O(N^{2})$ per time step, which poses challenges for large-scale applications.

In this paper, we develop a stochastic particle method for the Landau equation based on our previous work [14]. The method is built upon two key principles: collisions occur pairwise, and the grazing collision between two particles can be approximate by a diffusion process. Specifically, in each collision time window (of width $\Delta t$), particles are randomly paired, and the motion of each particle is governed by the diffusion process:

where $V_{\theta(i)}$ denotes the paired particle of $V_{i}$. In [14], we introduced this particle system and proved that, under regularity conditions on the coefficients, its empirical distribution converges to the solution of the Landau equation, with the convergence rate quantified using relative entropy. This system can be formally interpreted as a random batch method for the particle system (1.3) with a batch size of $p=2$ [22], but here the random batch structure arises naturally from physical collisions rather than as a purely approximation. We remark that the random batch approximation with any batch size $p>2$ is not a consistent approximation of McKean–Vlasov SDE (1.2) unless $p\to\infty$ as $N\to\infty$. The mechanism of the approximation in (1.4) is intrinsically different from (1.3).

This paper focuses on the numerical performance of the particle system, but with the particular choice of Brownian motions

A key structure of the system with this choice of Brownian motions is that the relative velocity is a spherical Brownian motion. By exploiting the unique structure of the system and utilizing spherical Brownian motion sampling, we construct an exact numerical scheme that preserves the particle distribution, ensuring consistency with the continuous system (see Algorithm 1). Overall, the proposed method offers several notable advantages:

• •

  Structure Preservation. The method preserves important physical properties of the Landau equation, such as the conservation of mass, momentum, and energy, as well as the dissipation of entropy, even in its temporal discretization implementation (see Section 3). These properties ensure that the numerical system stays consistent with the underlying physics, making it reliable for long-time simulations. Numerical experiments later in the paper show its stability and accuracy over extended time periods.

• •

  Accurate Particle Properties. By modeling collisions directly as pairwise interactions in the grazing regime, the method provides a detailed representation of individual particle properties, complementing approaches that focus on approximating the overall particle distribution. This helps achieve a more comprehensive understanding of particle dynamics.

• •

  Computational Efficiency. The computational cost per time step is $O(N)$, making the method highly efficient for large systems. Additionally, the grouping structure of the algorithm is naturally well-suited for parallel computation, further enhancing its scalability and efficiency.

• •

  Improved Sampling. Unlike traditional direct simulation Monte Carlo (DSMC) methods, which simulate Boltzmann collisions in the grazing regime to approximate Landau collisions, our approach directly simulates the Landau collision. This eliminates an extra layer of approximation and avoids the need for rejection sampling when updating velocities, reducing computational cost and making the method simpler to implement.

This method thus provides a practical and efficient tool for solving the Landau equation, offering both physical accuracy and computational efficiency.

The rest of the paper is organized as follows. Section 2 introduces the collisional particle system and demonstrates that it preserves the fundamental physical properties of the Landau equation. Section 3 describes the construction of the exact temporal discretization scheme and discusses its distinctive features. As a comparison, the Euler–Maruyama scheme is also presented, showing that it does not conserve energy but instead results in an energy increase over time. Section 4 focuses on the numerical implementation and evaluates the performance of the proposed method through selected examples. Finally, Section 5 illustrates the applicability of our method to spatially non-homogeneous equations using the Vlasov–Poisson–Landau equation as a case study.

We first present our particle system in detail, followed by some explanations and remarks. Then, we discuss several interesting properties of this system.

The equation (2.1) is equivalent to (1.4), which is written in Itô’s form, provided that $W_{i}=-W_{\theta(i)}$ in (1.4) as well. To verify this, let $Z_{i}=V_{i}-V_{\theta(i)}$, then the dynamics of $Z_{i}$ is given by $\mathop{}\mathopen{\mathrm{missing}}{d}Z_{i}=2\sigma(Z_{i})\circ\mathop{}\mathopen{\mathrm{missing}}{d}W_{i}$. Using the relationship between Stratonovich and Itô integrals, one has

A direct calculation shows that

which along with the fact $\sigma(z)z=0$ implies $\sigma\nabla\cdot\sigma(Z_{i})\equiv 0$. Recalling that $K=\nabla\cdot A=\nabla\cdot\sigma^{2}$, we confirm that (2.1) and (1.4) are indeed equivalent. We remark that the particle system was first introduced in our previous work [14]. However, to facilitate theoretical analysis, we assumed in [14] that the Brownian motions $W_{i}$ in system (1.4) are independent, rather than satisfying $W_{i}=-W_{\theta(i)}$; in this case, (1.4) cannot be rewritten in the form of (2.1). This assumption was discussed in Remark 2.3 of [14]. We believe that under both settings, the particle system converges to the same limit equation. Establishing this convergence under more general conditions will be a focus of our future work.

From the form of equation (2.1), we observe that $Z_{i}=V_{i}-V_{\theta(i)}$ evolves as a spherical Brownian motion (SBM) within each collision window $[t_{m},t_{m+1})$. Specifically, we have

where $I_{d}-(Z_{i}\otimes Z_{i})/|Z_{i}|^{2}$ is the projection operator onto the plane orthogonal to $Z_{i}$. This implies that the magnitude $|Z_{i}|$ remains constant during the interval $[t_{m},t_{m+1})$, confirming that $Z_{i}$ is indeed an SBM within this time window. This property is crucial because it not only ensures the well-posedness of the system but also allows us to derive an exact temporal discretization scheme for Particle System (CP) in the next section.

To intuitively explain why Particle System (CP) approximates the Landau equation, consider the evolution of the velocity distribution $f^{i}$ for the $i$-th particle during $[t_{m},t_{m+1})$, which satisfies:

Here, $f^{i}$ depends on both the random pairing $\theta(i)$ and the velocity $V_{\theta(i)}$. Assuming that all particles are independently and identically distributed (i.i.d.) at $t_{m}$ (since there is particle chaos when $N$ is large), the weak correlations during $[t_{m},t_{m+1})$ imply that the expected distribution $\tilde{f}^{i}$, averaged over all pairings and velocities, approximately satisfies

which corresponds to the Landau equation (1.1). Thus, it is reasonable to expect that Particle System (CP) converges to the Landau equation.

Next, we verify that Particle System (CP) preserves several key physical properties of the Landau equation, namely the conservation of mass, momentum, and energy, as well as entropy dissipation. The conservation of mass is straightforward. For momentum and energy, we establish the following result:

To justify the dissipation of entropy, we denote by $f^{N}$ the joint law of the $N$ particles in Particle System (CP) and define the (normalized) entropy as

It follows from (2.1) that, during each time interval $[t_{m},t_{m+1})$, the law $f^{N}(t,v_{1},\dots,v_{N})$ formally satisfies the Liouville equation:

Assuming that $f^{N}$ is properly regular, one can deduce that

as $A$ is semi-positive definite. Therefore, we have

Particle System (CP) is defined in continuous time. In practical implementations, we need to discretize it in time. The simplest approach is to use the Euler–Maruyama (EM) scheme (see Algorithm 2). However, as we will show later, the EM scheme does not have the desired structure preservation property: while the conservation of mass and momentum is automatically satisfied, the EM discretization fails to preserve the total energy of the particle system.

Fortunately, the unique structure of the system allows us to construct an exact time discretization scheme (see Algorithm 1). Unlike most other methods, our time discretization algorithm introduces no additional errors. This ensures that the probability distribution of the particles at the discrete time points matches exactly with that of the original continuous-time system, thereby preserving all the system’s desirable physical properties.

The key to achieving this lies in the fact that the relative velocity $Z_{i}=V_{i}-V_{\theta(i)}$ evolves as a spherical Brownian motion. According to Stroock’s representation, a standard spherical Brownian motion $Y_{t}$ on $\mathbb{S}^{d-1}$ satisfies the following SDE:

Comparing this with (2.3), it follows that the rescaled spherical Brownian motion $|Z_{i}|Y_{kt}$, with time-scaling coefficient $k=4\Lambda|Z_{i}|^{\gamma}$, provides a (weak) solution to (2.3).

In practice, we use existing methods [21, 25] for the exact or approximate simulation of SBM. For instance, in two dimensions, SBM on the unit circle can be simulated exactly using the explicit formula $(\cos(B_{t}),\sin(B_{t}))$, where $B_{t}$ is standard Brownian motion. For higher dimensions, SBM increments can still be simulated exactly by decomposing the motion into radial and angular components, as described in [21, 25]. The angular component is uniformly distributed on $\mathbb{S}^{d-2}$, while the radial component is governed by the Wright–Fisher diffusion. The Wright–Fisher distribution can be either exactly simulated or well-approximated when the time step $\Delta t$ is small (e.g., $\Delta t\leq 0.05$).

Next, we present the Euler–Maruyama scheme for the system, and discuss the energy errors produced by this discretization.

In the EM scheme, the conservation of momentum is automatically preserved. For the energy, we have following result.

In this section, we illustrate the applicability of our method to spatially non-homogeneous equations, in particular the Vlasov–Poisson–Landau equation. Application to other models with Landau collision is similar.

In spatially non-homogeneous cases, space heterogeneity generates an electric field if the particles are charged, which affects both the velocities and positions of particles. To describe this phenomenon, the Vlasov–Poisson–Landau equation is formulated as:

where $f(x,v,t)$ is the particle distribution, $\rho(x,t)=\int_{\Omega_{v}}f(x,v,t)\mathop{}\mathopen{\mathrm{missing}}{d}v$, $\rho_{0}$ is the background charge density, $\phi$ is the potential, and $E$ is the electric field. When periodic boundary conditions are employed, the charged neutral condition must be imposed for the finiteness of the system energy so that $\int_{\Omega_{x}}(\rho-\rho_{0})\,dx=0$. For example, in plasma simulations, $f$ may represent the density of electrons and then $\rho_{0}=\frac{1}{|\Omega_{x}|}\int_{\Omega_{x}}\rho\,dx$ represents the density of ions (considering that electrons have negative charge, one may reverse the sign of $\rho$ and thus $\phi,E$ above, while $-E$ would be used in the Vlasov equation so there would be no intrinsic difference).

The total energy $\mathcal{E}_{\text{total}}$ comprises kinetic energy $\mathcal{E}_{K}$ and electric energy $\mathcal{E}_{E}$, given by:

The conservation of $\mathcal{E}_{\text{total}}$ is a critical property of the system.

To solve the Vlasov–Poisson–Landau equation, we adopt a time-splitting approach that decouples the collision step

and the advection step

Traditionally, some studies have employed the DSMC method to address the collision step in the Vlasov–Poisson–Landau equation [24, 32]. This approach involves simulating Boltzmann collisions to approximate the Landau collision. While it has been effectively combined with other schemes for solving the advection step, the method requires rejection sampling when updating velocities, which may increase computational cost and time.

For  (4.2), we use the SBM scheme to update particle velocities in the same grid, preserving the kinetic energy $\mathcal{E}_{K}$. Since particle positions remain unchanged during this step, the electric field and $\mathcal{E}_{E}$ are also unchanged.

For  (4.3), which corresponds to the Vlasov–Poisson equation, the particle-in-cell (PIC) algorithm is widely employed to solve such equations [1, 28]. As a preparatory step for the algorithm, the spatial domain $\Omega_{x}$ is divided into a set of grids $\Omega_{k}$. Using the electric field $E_{k}$ at the center of each grid $\Omega_{k}$, the electric energy $\mathcal{E}_{E}$ can be approximately computed as:

where increasing the number of grids improves the accuracy of $\mathcal{E}_{E}$. Let $x_{k}$ denote the center of the $k$-th grid, $x_{i}$ the position of the $i$-th particle, $Q$ the total charge in the whole area, and $q$ the charge of each particle. The PIC algorithm proceeds in two main steps:

Step 1: Solving the Poisson equation. The charge density at the grid point $x_{k}$ is approximated by:

where $q=Q/N$, $S(x)$ is a shape function used to approximate the $\delta$-function, and $\rho_{0}:={|\Omega_{x}|^{-1}}\int_{\Omega_{x}}\rho\mathop{}\mathopen{\mathrm{missing}}{d}x$. Then the Poisson equation $-\Delta_{x}\phi=\rho-\rho_{0}$ is solved with periodic boundary conditions. The spectral scheme is utilized to obtain $\phi$, from which the electric field at the center of the $k$-th grid, $E(x_{k})$, is computed as $E(x_{k})=-\nabla_{x}\phi(x_{k})$.

Step 2: Updating velocities and positions. To update the particle velocities and positions, the electric field at the particle positions $E(x_{i})$ is needed. Using the same shape function $S(x)$, the electric field at a general location $x$ is interpolated as:

Subsequently, the leap-frog scheme is typically employed to update the particle velocities and positions, offering high accuracy [3].

However, the leap-frog scheme in this step does not conserve the total energy $\mathcal{E}_{\text{total}}$. To address this, the Vlasov–Poisson equation can be equivalently reformulated as the Vlasov–Ampère equation, and the Crank–Nicolson (CN) method [10] can be used instead. The CN method ensures strict conservation of both the total energy $\mathcal{E}_{\text{total}}$ and the total charge. The specific form is formulated as:

where $x^{n+1/2}_{i}:=\frac{1}{2}(x^{n}_{i}+x^{n+1}_{i}),\quad v^{n+1/2}_{i}:=\frac{1}{2}(v^{n}_{i}+x^{n+1}_{i}),\quad J_{\text{mean}}:=\frac{1}{n_{0}}\sum_{k}J^{n+1/2}(x_{k}).$

By combining the SBM scheme, which preserves the kinetic energy $\mathcal{E}_{K}$ in the Landau collision step, with the energy-conserving PIC algorithm, the total energy $\mathcal{E}_{\text{total}}$ can be rigorously preserved throughout the simulation. See Algorithm 3 for details.

In this section, we validate the effectiveness of our method through numerical experiments for solving the Landau equation with different values of $\gamma$ and dimensions. We compare our numerical solutions with the analytical Bobylev–Krook–Wu (BKW) solution for Maxwell molecules and the reference solution given in [9] for the Coulomb potential case. Besides, we demonstrate the reliability and accuracy of our approach for the Vlasov-Poisson-Landau equation by studying the effects of Landau collision in the phenomenon of Landau damping.

To visualize the particle solution and compare it with the exact (or reference) solution, we construct a mollified solution $f^{N}_{\epsilon}$ from the empirical measure of particle velocities as follows:

where $\psi_{\epsilon}(x)$ is the Gaussian mollifier. Compared to deterministic particle methods such as [8, 9], the parameter $\epsilon$ is only used as a post-processing parameter and can therefore be chosen more flexibly. For visualization purposes and consistency, we set $\epsilon=0.01$ uniformly across all experiments.

To measure the accuracy of our solution, we use the relative $L_{2}$ error, which is computed on a uniform mesh grid with center points $v_{l}^{c}$ of the squares as follows:

We also evaluate the entropy of the mollified solution $f^{N}_{\epsilon}$ using the following expression:

Recall that we have shown in Section 2 that the entropy of the joint distribution $H_{N}(f^{N})$ is monotonically decreasing. However, in practical computations, directly evaluating the entropy of the joint distribution of particles is infeasible because obtaining the corresponding density function is challenging. By smoothing the empirical distribution of the particles, the resulting mollified density function $f^{N}_{\epsilon}$ approximates the marginal distribution of the joint distribution. Under the assumption that all particles are approximately i.i.d., the following relationship holds approximately:

Thus, in our numerical experiments, we compute the entropy of the mollified empirical measure $f_{\epsilon}^{N}$ to verify the entropy dissipation property.