Vlasov-Poisson-Landau Equations Overview

• Vlasov–Poisson–Landau equations are kinetic plasma models that couple Vlasov transport with Landau’s Coulomb collision operator to capture charge interactions in plasmas.
• The system features distinct formulations and macro–micro decompositions that enable analysis of conservation laws, global existence, and decay near Maxwellian equilibria.
• Recent numerical approaches employ structure-preserving methods to accurately simulate collisional relaxation, Landau damping, and singular limits.

Vlasov–Poisson–Landau (VPL) equations are kinetic plasma models that couple Vlasov transport under a self-consistent electrostatic field with Landau’s Coulomb-collision operator. In the normalized two-species setting, one seeks nonnegative ion and electron distributions $F_\pm(t,x,v)$ solving

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$

with

$$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$

This framework appears in whole-space, toroidal, bounded-domain, weakly collisional, and reduced-model regimes, and it supports both perturbative global theory near Maxwellians and singular-limit analysis toward fluid or dispersive equations (Strain et al., 2012).

1. Canonical formulations and conservation laws

The basic VPL system occurs in several mathematically distinct forms. In the whole-space two-species problem on $\mathbb R^3_x\times\mathbb R^3_v$, the Poisson source is the charge disparity $F_+-F_-$, and the collision operator accounts for both intra-species and inter-species Coulomb interactions. On $\mathbb T^3_x\times\mathbb R^3_v$ in the weakly collisional regime, the model is often written with an explicit small collision frequency $\nu\ll1$,

$$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$

with the electric field determined after subtraction of the spatially averaged density. In bounded domains $\Omega\subset\mathbb R^3$, one imposes specular reflection,

$$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$

on the inflow boundary. A further variant models ions in one spatial dimension with electrons prescribed by the Maxwell–Boltzmann law $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$0, leading to

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$1

These are not mere notational changes: they encode different scalings, boundary mechanisms, and asymptotic limits (Duan et al., 2023).

For the ion model under the Maxwell–Boltzmann relation, the Landau operator conserves mass, momentum, and energy and satisfies an $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$2-theorem. In the Hamiltonian part of the dynamics, the collisionless Vlasov–Poisson equations admit a non-canonical Poisson bracket, conserve the Hamiltonian

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$3

and possess Casimirs such as total mass and entropy. In metriplectic formulations, the collisional step is generated by a symmetric positive-semidefinite metric bracket for which entropy is non-decreasing while mass, momentum, and energy remain collision invariants. These structural facts are central in both analysis and discretization because they identify which parts of the dynamics are conservative, which are dissipative, and which degeneracies must be respected rather than regularized away (Finn et al., 13 Feb 2026).

2. Equilibria, linearization, and macro–micro structure

The perturbative theory is organized around Maxwellian equilibria. In the whole-space near-equilibrium setting, the global Maxwellian is

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$4

and one writes

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$5

The perturbation then satisfies

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$6

where $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$7 is the linearized Landau operator and $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$8 is the quadratic remainder. In this formulation $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$9 is nonnegative, and its null-space is spanned by $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$0. The central coercivity statement is that $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$1 controls the microscopic component $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$2, not the hydrodynamic modes; hence Poisson coupling and transport must be used to recover dissipation in the macroscopic directions (Strain et al., 2012).

For asymptotic and multiscale problems, the relevant reference state is often a local Maxwellian $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$3. Writing $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$4, one introduces the local linearized Landau operator

$$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$5

whose null-space is spanned by the five collision invariants $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$6. The associated micro–macro decomposition uses the projection $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$7 onto these invariants and $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$8, so that $$Q(G,H)(v)=\nabla_v\!\cdot\!\int_{\mathbb R^3}\Phi(v-v_*) \bigl\{G(v_*)\,\nabla_vH(v)-H(v)\,\nabla_{v_*}G(v_*)\bigr\}\,dv_*, \qquad \Phi(z)=\frac1{|z|}\Bigl(I-\frac{z\otimes z}{|z|^2}\Bigr).$$9 with $\mathbb R^3_x\times\mathbb R^3_v$0. Inverting $\mathbb R^3_x\times\mathbb R^3_v$1 on the orthogonal complement generates fluid-type equations with kinetic source terms involving $\mathbb R^3_x\times\mathbb R^3_v$2, and the leading microscopic operator in the scaled equations becomes $\mathbb R^3_x\times\mathbb R^3_v$3. This decomposition is the analytic backbone of both the KdV limit and weighted-energy closures (Duan et al., 2023).

In the whole-space two-species theory, the same macro–micro mechanism is encoded through projections $\mathbb R^3_x\times\mathbb R^3_v$4 onto the null-space of $\mathbb R^3_x\times\mathbb R^3_v$5 and dissipative norms involving the Landau collision weight $\mathbb R^3_x\times\mathbb R^3_v$6. A key estimate is

$\mathbb R^3_x\times\mathbb R^3_v$7

which is paired with weighted mixed Sobolev energies in $\mathbb R^3_x\times\mathbb R^3_v$8 and with control of $\mathbb R^3_x\times\mathbb R^3_v$9. This suggests a general principle for VPL analysis: the collision operator provides anisotropic, degenerate velocity diffusion, while transport and Poisson coupling supply the missing coercivity through hypocoercive interactions (Wang, 2012).

3. Global existence, decay, and regularization near Maxwellians

A foundational perturbative result is the global small-data theory in $F_+-F_-$0. For integers $F_+-F_-$1 and $F_+-F_-$2, Strain–Zhu define a weighted Sobolev norm

$F_+-F_-$3

assume $F_+-F_-$4, and prove global existence and uniqueness together with

$F_+-F_-$5

They also obtain polynomial convergence to equilibrium: $F_+-F_-$6 for any small $F_+-F_-$7. The proof combines coercivity of the linearized Landau operator, a time-frequency Lyapunov functional in Fourier space, and nonlinear weighted Sobolev estimates (Strain et al., 2012).

Wang’s whole-space two-species analysis sharpens the description of decay by separating the total density from the charge disparity. Near a global Maxwellian, one obtains a unique global solution with weighted energy bounds; for low moments, the potential and the difference $F_+-F_-$8 satisfy

$F_+-F_-$9

and, under stronger exponential velocity weights, they even obey $\mathbb T^3_x\times\mathbb R^3_v$0 decay. By contrast, the sum $\mathbb T^3_x\times\mathbb R^3_v$1 decays at the optimal algebraic rate associated with the Landau equation in the whole space, assuming an additional $\mathbb T^3_x\times\mathbb R^3_v$2-moment or $\mathbb T^3_x\times\mathbb R^3_v$3 bound on the data. This decomposition clarifies that the Poisson field accelerates relaxation of charge imbalance while leaving the conserved total density in the slower Landau regime (Wang, 2012).

Regularization theory adds a different layer. For hard and soft potentials with $\mathbb T^3_x\times\mathbb R^3_v$4, the smoothing analysis in $\mathbb T^3_x\times\mathbb R^3_v$5 proves that any classical solution becomes immediately smooth with respect to all variables. The argument uses a time-weighted energy method combined with pseudo-differential calculus, with weights $\mathbb T^3_x\times\mathbb R^3_v$6 that vanish at $\mathbb T^3_x\times\mathbb R^3_v$7 for high derivatives and symbols $\mathbb T^3_x\times\mathbb R^3_v$8 designed to recover fractional gains in $\mathbb T^3_x\times\mathbb R^3_v$9 and $\nu\ll1$0. Under the global smallness assumption, the same framework yields optimal large-time decay for hard potentials,

$\nu\ll1$1

and slightly weaker high-order decay for soft potentials (Deng, 2021).

The bounded-domain theory with specular reflection shows that this near-Maxwellian stability is not restricted to translationally invariant geometries. In a general bounded domain, small high-velocity-weighted $\nu\ll1$2 and $\nu\ll1$3 perturbations generate a unique global solution, with uniform weighted energy bounds, almost-exponential time decay in both $\nu\ll1$4 and $\nu\ll1$5, and Hölder as well as velocity-derivative bounds. The key methodological shift is an $\nu\ll1$6 framework in which macro–micro energy estimates are combined with $\nu\ll1$7 theory for an ultraparabolic equation obtained after flattening and reflecting the boundary. A common misconception is that specular reflection is only a geometric nuisance; in this setting it changes the analytic architecture by replacing whole-space Fourier tools with kinetic De Giorgi–Nash–Moser-type regularity (Dong et al., 2020).

4. Weakly collisional dynamics, Landau damping, and collisionless behavior

On $\nu\ll1$8 with small collision frequency $\nu\ll1$9, the weakly collisional VPL system exhibits both Landau damping and collisional relaxation. Chaturvedi–Luk–Nguyen prove that if the initial perturbation is $$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$0 in weighted Sobolev norms involving the commuting vector fields

$$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$1

then the solution exists globally, remains $$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$2, and converges to the global Maxwellian. More precisely, for each nonzero Fourier mode,

$$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$3

while

$$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$4

The analytical mechanism couples Guo’s weighted energy method, hypocoercive cross terms, the exact commutation $$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$5 for suitable $$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$6, and pointwise resolvent estimates for the Volterra equation satisfied by the density modes (Chaturvedi et al., 2021).

A later Gevrey analysis establishes a stronger nonlinear regime. For $$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$7-independent perturbations of the global Maxwellian in Gevrey-$$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$8, one obtains uniform-in-$$\partial_t F + v\!\cdot\!\nabla_x F + E\!\cdot\!\nabla_v F=\nu\,Q(F,F),$$9 Landau damping and enhanced dissipation,

$\Omega\subset\mathbb R^3$0

for any $\Omega\subset\mathbb R^3$1, together with a quantitative collisionless limit: if $\Omega\subset\mathbb R^3$2 solves VPL and $\Omega\subset\mathbb R^3$3 solves the collisionless Vlasov–Poisson equation with the same initial data, then

$\Omega\subset\mathbb R^3$4

The proof introduces an infinite-regularity commuting vector field method, a nearly physical-side treatment of collisionless echoes, a decomposition of collisional contributions in the density Volterra equation, and a quasi-linearization procedure for the slowly evolving homogeneous mode at long times (Bedrossian et al., 24 Aug 2025).

These results refute a simplistic dichotomy between “collisionless damping” and “collisional relaxation.” In the small-$\Omega\subset\mathbb R^3$5 regime, phase mixing persists on $\Omega\subset\mathbb R^3$6 time scales, enhanced dissipation emerges on $\Omega\subset\mathbb R^3$7, and the nonlinear collisional dynamics remain close to the collisionless prediction for $\Omega\subset\mathbb R^3$8. The VPL system therefore furnishes a mathematically controlled interpolation between Vlasov–Poisson and collisional kinetic behavior, rather than a discontinuous transition from one to the other.

5. Singular limits and reduced plasma models

One important reduction is the KdV limit for ions under the Maxwell–Boltzmann electron relation. In one space dimension, after the Gardner–Morikawa scaling

$\Omega\subset\mathbb R^3$9

with

$$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$0

the rescaled VPL system admits a Chapman–Enskog expansion around local Maxwellians. Matching powers of $$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$1 yields

$$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$2

and the leading density profile satisfies

$$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$3

The rigorous theorem constructs unique global smooth solutions on any fixed interval $$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$4, proves the energy bound $$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$5, and shows

$$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$6

The method is a weighted energy analysis of the micro–macro system around local Maxwellians (Duan et al., 2023).

A different singular reduction is the massless electron limit of the two-species VPL system. Writing $$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$7 and rescaling the electron velocity produces multiple velocity scales. In the formal $$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$8 limit, the electron equation loses its time derivative at leading order, and entropy–dissipation forces the electron distribution to become a local Maxwellian,

$$F(t,x,v)=F(t,x,R_xv),\qquad R_xv=v-2(v\!\cdot\!n(x))\,n(x),$$9

A crucial structural identity shows that the ion–electron Landau term drops out identically, so the limiting ion equation closes as

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$00

with an additional ODE for $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$01 from energy conservation. The rigorous convergence theorem gives

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$02

and a faster decay of a lower-weight electron error norm. The same work shows that the analogous vanishing mechanism fails for the hard-sphere Boltzmann cross-collision operator, whose leading $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$03 term survives on Maxwellians. This is a genuine model distinction, not a technical artifact (Flynn et al., 2023).

These asymptotic limits illustrate two different reduction paradigms. The KdV limit derives a dispersive scalar equation from a kinetic plasma model under long-wave, small-amplitude, collisional scaling. The massless-electron limit instead removes an entire transport equation and replaces electrons by a Maxwell–Boltzmann background, thereby clarifying why a direct two-fluid hydrodynamic limit is obstructed by ion–electron collisions in the full VPL system but not in the limiting Landau structure identified by the scaling.

6. Structure-preserving computation

The VPL system is unusually demanding numerically because its transport part is Hamiltonian, its collisional part is dissipative but invariant-preserving, and the Coulomb kernel is both nonlocal and singular. One modern response is a particle-in-cell discretization combined with discrete-gradient time integrators. In this framework, the collisionless Vlasov–Poisson step is advanced by a discrete-gradient Poisson integrator,

$$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$04

which exactly preserves the discrete Hamiltonian, while the collisional Landau step is advanced by a discrete-gradient-dependent integrator acting on particle velocities. The scheme guarantees conservation of mass, momentum, and energy and preserves monotonicity of entropy production in both the time-continuous and discrete systems. Implemented in PETSc using TS “discgrad,” DMPlex, DMParticle, and SNES, it reproduces benchmark behaviors such as Landau damping and electron–positron equilibration while maintaining bounded energy drift and essentially monotone discrete entropy up to solver tolerance (Finn et al., 13 Feb 2026).

An earlier deterministic conservative solver adopts a different architecture based on splitting. Over each time step, one solves a Vlasov–Poisson subproblem with a Runge–Kutta discontinuous Galerkin method and a homogeneous Landau subproblem with a conservative spectral method. Conservation across the two grids is enforced by a dedicated routine in Fourier space that projects the computed collision term onto the subspace satisfying the moment constraints for $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$05. The resulting method preserves mass, momentum, and energy to discretization precision, supports hybrid MPI+OpenMP parallelization, and has been tested on linear and nonlinear Landau damping, temperature relaxation, and two-stream instability; in the two-stream example, total energy remains conserved to $$\partial_tF_\pm + v\!\cdot\!\nabla_x F_\pm \pm E\!\cdot\!\nabla_vF_\pm = Q(F_\pm,F_\pm)+Q(F_\mp,F_\pm), \qquad -\Delta_x\phi=\int_{\mathbb R^3}(F_+-F_-)\,dv, \quad E=-\nabla_x\phi,$$06 (Zhang et al., 2016).

These numerical approaches reflect a broader shift in VPL computation from mere consistency toward geometric and thermodynamic compatibility. For this equation class, conservation of moments and non-decreasing entropy are not optional diagnostics; they are discrete surrogates for the same structural constraints that underlie coercivity estimates, long-time relaxation, and asymptotic limits in the analytical theory.