Source-linked AI summary
A discontinuous Galerkin method for the Vlasov-Poisson system
R. E. Heath, I. M. Gamba, P. J. Morrison, C. Michler
TL;DR
The Vlasov-Poisson system is challenging to simulate because filamentation produces stiff gradients. This paper evaluates a discontinuous Galerkin approach on linear and nonlinear problems, reporting comparisons with theory and established methods that demonstrate its viability.
Problem
Filamentation produces stiff gradients in the Vlasov-Poisson system, making numerical simulation challenging.
Source: §1 · p. 2Method
The paper applies discontinuous Galerkin approximations on phase-space meshes, including a piecewise-constant scheme for its plasma simulations.
Source: §3.1 · p. 10, §3.5 · p. 20Results
The computed examples are benchmarked against theory and established methods, and Maxwellian linearized-system results compare very well with theoretical Landau-damping calculations.
Source: §1 · p. 4, §5 · p. 55Takeaways & Limitations
The paper concludes that the DG method with upwinding is viable for solving the Vlasov-Poisson system.
Source: §5 · p. 55The piecewise-constant scheme preserves positivity, but finding a scheme that also preserves higher-order moments such as momentum and energy remains a challenge.
Source: §3.5 · p. 20
Abstract
from arXiv · showhide
A discontinuous Galerkin method for approximating the Vlasov-Poisson system of equations describing the time evolution of a collisionless plasma is proposed. The method is mass conservative and, in the case that piecewise constant functions are used as a basis, the method preserves the positivity of the electron distribution function and weakly enforces continuity of the electric field through mesh interfaces and boundary conditions. The performance of the method is investigated by computing several examples and error estimates associated system's approximation are stated. In particular, computed results are benchmarked against established theoretical results for linear advection and the phenomenon of linear Landau damping for both the Maxwell and Lorentz distributions. Moreover, two nonlinear problems are considered: nonlinear Landau damping and a version of the two-stream instability are computed. For the latter, fine scale details of the resulting long-time BGK-like state are presented. Conservation laws are examined and various comparisons to theory are made. The results obtained demonstrate that the discontinuous Galerkin method is a viable option for integrating the Vlasov-Poisson system.
1 Introduction
Filamentation and Landau damping make the Vlasov-Poisson system challenging to simulate, while existing approaches have limitations in some regimes. The paper proposes a coupled Upwind Penalty Galerkin method and evaluates it on linear and nonlinear problems.
1 Introduction
- Filamentation creates stiff gradients, and Landau damping describes exponential collisionless decay of electric-field modes, complicating numerical simulation.Filamentation arises as characteristics wrap through phase space, leaving disparate values of the distribution along different characteristics.Source: §1 · p. 2
- Particle-in-Cell methods can suffer numerical noise proportional to 1/√n when the distribution tail is not negligible and many particles are needed.Fourier-transform approaches with filamentation filtering may be inadequate when the physics of interest depends on filamentary structure, as in Landau damping.Source: §1 · p. 3
- The coupled Upwind Penalty Galerkin method uses an upwind Galerkin scheme for the Vlasov equation and a DG interior-penalty scheme for Poisson’s equation.The paper describes the method for general polynomial bases; it is locally and globally mass conservative.Source: §1 · p. 3
- Piecewise-constant basis functions preserve positivity, while the method’s local and discontinuous character supports mesh refinement around filamentation and stiff gradients.Weakly imposed boundary conditions also accommodate varied boundary conditions.Source: §1 · pp. 3–4
- The method is evaluated through accuracy, convergence, and conservation tests on linear and nonlinear problems, benchmarked against theory and established methods.The paper includes advection, Landau damping, and two-stream instability examples.Source: §1 · pp. 4–5
- The introduction situates the proposal among DG work on kinetic equations, including a positivity-preserving Vlasov-Boltzmann scheme that also supports higher-order polynomial approximations.The paper notes that a higher-order UPG approximation for the Vlasov system is under development.Source: §1 · p. 4
2 The Vlasov-Poisson System
The section sets out the nondimensionalized Vlasov-Poisson model, its boundary conditions, and its conservation properties. It also introduces linearization about an equilibrium to study Landau damping.
2 The Vlasov-Poisson System
- The model is nondimensionalized using characteristic plasma time, Debye length, and thermal velocity scales.Time is scaled by the inverse plasma frequency.Source: §2 · p. 5
- The Vlasov-Poisson model couples transport of the electron distribution in phase space to an electrostatic field and potential.The transport equation uses the divergence-free phase-space vector α(x, v, t) = (v, −E(v, t)).Source: §2 · pp. 5–6
- For bounded spatial domains, the distribution has an inflow-boundary condition; this work assumes periodic spatial boundaries and decay in velocity.The Poisson equation may use Dirichlet, Neumann, Robin, or periodic spatial boundary conditions; under homogeneous or periodic conditions with no Dirichlet boundary, neutrality is required.Source: §2 · p. 6
- The system conserves particle number, momentum, energy, and Casimir invariants, and the paper checks these invariants in its nonlinear computations.The Casimir family includes enstrophy, entropy, and particle number.Source: §2 · pp. 6–7
- Linearization around an equilibrium distribution yields a system used to study the exponential long-time decay of electric-field modes known as Landau damping.The stated linearization assumes the electric-field perturbation term is small and drops it from the perturbed Vlasov equation.Source: §2 · p. 7
3 Method Formulation
The section formulates the UPG method by coupling an upwind discontinuous Galerkin approximation of the Vlasov equation with an interior-penalty Galerkin Poisson solver, using product meshes and discontinuous elementwise polynomial approximations on a bounded velocity domain. It describes the iteration between the electric field and distribution function, an error estimate, and a positivity-preserving piecewise-constant scheme.
3 Method Formulation
- The phase-space mesh is the Cartesian product of independently refined spatial and velocity meshes, and the approximation is polynomial on each element without requiring continuity across elements.The Poisson spatial mesh must match the Vlasov spatial mesh, although their polynomial degrees may differ.Source: §3 · p. 8, §3.3 · p. 14
- Coupling the two discretizations yields a nonlinear semi-discrete first-order ODE system for the distribution approximation, with the electric field and potential computed through the Poisson approximation.The system can be advanced with an explicit conservative time integrator such as Runge–Kutta.Source: §3 · p. 9
3.1 Preliminaries
- The computational velocity domain is bounded, and the approximation is assumed to vanish on its boundary; the resulting error depends on the true solution’s density outside that domain.For spatially periodic problems, the authors consider this reasonable when the velocity domain is sufficiently large to contain most of the approximate density.Source: §3.1 · p. 9
- The approximation spaces use tensor-product Q polynomials on product meshes, while P spaces are also suitable and may reduce degrees of freedom in higher dimensions.For n-dimensional calculations, the stated basis counts are (r + 1)^n for Qr(K) and r^(n+1) for Pr(K).Source: §3.1 · p. 10
3.2 Upwind Galerkin Approximation of the Vlasov Equation
- The Vlasov discretization replaces the face value with an upwind numerical flux, selecting the trace from the element that the flow field α enters across the face.This general formulation covers phase spaces of dimension 2n for n = 1, 2, or 3, with inflow boundary conditions and arbitrary-degree piecewise polynomials.Source: §3.2 · pp. 11–12
3.3 Interior Penalty Approximations of the Poisson Equation
- The Poisson approximation uses SIPG, IIPG, or NIPG schemes, which weakly enforce approximate continuity across interior faces and Dirichlet boundary regions.The three schemes differ in the parameter c_s, set to -1, 0, and 1 for SIPG, IIPG, and NIPG, respectively.Source: §3.3 · pp. 14–16
- The Poisson bilinear form yields a uniquely solvable matrix system, subject to the stated boundary and compatibility conditions.For periodic boundaries, the conservative UPG scheme satisfies the compatibility condition at the numerical level; the resulting sparse matrix is non-block-diagonal.Source: §3.3 · p. 17
3.4 Discontinuous Galerkin Approximation of the Vlasov-Poisson System
- The UPG method computes the Poisson potential and its local gradients to form a discontinuous electric field, then uses that field in the Vlasov ODE system.This combines the UG approximation of the Vlasov equation with an interior-penalty approximation of Poisson.Source: §3.4 · pp. 17–18
- The coupled ODE system can be advanced by an explicit integrator using a Gummel-map iteration; the simulations use a conservative high-order Runge-Kutta integrator.At each step, the field-dependent update is computed before advancing the distribution function.Source: §3.4 · p. 18
- The semi-discrete error estimate for the potential, electric field, and distribution function is optimal under the stated smoothness conditions.The estimate concerns suitable norms of Φ−Φh, E−Eh, and f−fh.Source: §3.4 · p. 19
3.5 Two-Dimensional Phase Space with Piecewise Constant Approximation
- Piecewise-constant approximations in the two-dimensional phase space give a positivity-preserving scheme used for the plasma simulations.The approach extends to higher dimensions, but preserving positivity together with higher-order moments such as momentum and energy remains a challenge.Source: §3.5 · p. 20
- With a given electric field, the piecewise-constant approximation reduces to a linear ODE system whose interior-face terms are integrated over element interfaces.The basis is formed from products of piecewise-constant spatial and velocity functions.Source: §3.5 · pp. 20–21
4 Numerical Results
The numerical results assess the method across linear advection and linear and nonlinear Landau damping, including Maxwellian and Lorentzian equilibria where specified. They also examine nonlinear mode evolution, entropy dissipation, phase-space smoothing, and the two-stream run’s evolution toward a BGK-like state, alongside conservation behavior and fine-scale features.
4 Numerical Results
- The section presents six tests: four for linear dynamics and phase-space filamentation, and two for nonlinear evolution; piecewise constants approximate f, quadratic polynomials approximate Φ, and time uses a conservative fourth-order Runge-Kutta method.The tests include nonlinear Landau damping and a two-stream instability calculation; this part reports the nonlinear Landau damping results.Source: §4 · pp. 21–22, §4.2 · p. 34
4.1.1 Advection Results
- For Maxwellian linear advection, the computed net charge density converges under mesh refinement toward the analytic decay max_x |ρtot(x, t)| = 0.1 e^−t²/8.The results agree with earlier works at early times but show no later recurrence, which the authors attribute to the fine mesh and the DG method’s monotonic dissipative error.Source: §4.1.1 · pp. 25–26, Figure 1 · p. 26
- For Lorentzian linear advection, computations match the theoretical decay at early times but diverge later, with larger wavenumbers harder to resolve.The Lorentz equilibrium’s heavy velocity tail requires a larger velocity domain and provides a more stringent numerical test; its phase-mixing decay also shares analytic structure with Landau damping.Source: §4.1 · p. 23, §4.1.1 · p. 26
4.1.2 Linear Landau Damping Results
- For Lorentzian linear Landau damping, numerical rates and oscillation frequencies match theory to two decimal digits across k = 1/8, 1/6, 1/4, and 1/2.The k = 1/2 case damps for less time because velocity-space filamentation develops more rapidly; the heavy-tailed equilibrium also makes this a stringent truncation test.Source: §4.1.2 · pp. 29–30, Figure 5 · p. 33
- Under mesh refinement, Maxwellian linear Landau damping converges to the theoretical decay rate −0.153 and frequency 1.415, each to three decimal digits.The decay rate remains sustained longer on finer meshes before leveling off; the computed distribution also develops increasing filamentation over time.Source: §4.1.2 · p. 28, Figure 4 · p. 32
4.2.1 Nonlinear Landau Damping
- In nonlinear Landau damping, mode-one reaches a minimum near t ≈ 15, all four modes peak near t ≈ 40, and the bounce time is about 20.The damping rate substantially exceeds the linear rate of γ = −0.153 because higher modes are excited early and energy transfers from the first mode through nonlinear interactions.Source: §4.2.1 · p. 34
- The first-mode damping rate is γ = −0.287, consistent with −0.281 and −0.243; differing fits, finer-mesh coupling, and short-time entropy dissipation may explain deviations.Over the run to T = 120, total energy is conserved within a few percent and entropy is well conserved later in the run.Source: §4.2.1 · pp. 34–35, Figure 7 · p. 37
- The spatial average develops an early plateau near the wave phase velocity, then smooths around t ≈ 40 into a nearly constant state with a persistent electron hole.A small periodic dimpling at the maximum persists at late times.Source: §4.2.1 · p. 35, Figure 9 · p. 45
4.2.2 Nonlinear Two-Stream Instability
- The two-stream run progresses from exponential linear growth through trapping and mode saturation to an apparent BGK-like electron-hole state; early growth agrees with linear theory.The first four modes peak near t ≈ 18, with growth rates on the order of a couple of tenths of ωp.Source: §4.2.2 · pp. 38–39, Figure 12 · p. 48, Figure 13 · p. 48
- Over 100 time units, particle-number error stays ≤0.01% (nonlinear iteration contributes about 10^-4), momentum is exact, enstrophy is within about 10%, and total-energy relative error is below 4%.Kinetic and electrostatic energy oscillations cancel in their sum; particle-number and enstrophy errors increase monotonically.Source: §4.2.2 · pp. 39–40, Figure 14 · p. 49, Figure 15 · p. 50
- At t = 100, f100 plotted against particle energy appears as a graph across all 9 million uniformly meshed phase-space pairs; piecewise-constant DG keeps f100 > 0, and the authors identify a saturated BGK state.Source: §4.2.2 · pp. 41–42, Figure 16 · p. 50
- The fitted two-stream model gives a reasonable fit to the t = 100 distribution and potential, but underrepresents charge near the distribution’s maximum.The fit uses ΦM = 1.06, E∗ = 1.59, β = 1, and a = 0.1148; the authors describe it as a rough self-consistency test.Source: §4.2.2 · p. 43, Figure 17 · p. 51
- Fine-scale f100 features include a cusp at the E = 0 trapping boundary, steepening near minimum energy, and direction-independent splitting from E ≈ 1.3 that suggests incomplete saturation.The cusp and steepening may reflect either the dynamics or the numerical scheme; the authors leave their interpretation for future study.Source: §4.2.2 · p. 42, Figure 18 · p. 52, Figure 19 · p. 53
- E2 appears as a graph over Φ even in unsaturated states, and the model provides a reasonable fit to the t = 100 E2-versus-Φ data.This comparison is presented as a fair indication of self-consistency, not as an optimized fit.Source: §4.2.2 · p. 44, Figure 20 · p. 54, Figure 21 · p. 54
5 Conclusions
The paper demonstrates the viability of an upwind discontinuous Galerkin method for the Vlasov-Poisson system through benchmark and nonlinear computations. It also identifies future work in higher-order bases, adaptive and parallel implementations, and additional plasma applications.
5 Conclusions
- The DG method with upwinding is demonstrated as viable for Vlasov-Poisson problems, with results compared against theory across linear and nonlinear examples.The computations include linear advection, Maxwellian and Lorentzian Landau damping, nonlinear Landau damping, and the symmetric two-stream instability.Source: §5 · p. 55
- Lorentzian-equilibrium simulations verify the wave-number dependence of Landau damping, apparently for the first time in a simulation.For the Maxwellian equilibrium, computed results also compare very well with theoretical Landau-damping calculations.Source: §5 · p. 55
- In nonlinear Landau damping and two-stream calculations, monitored constants of motion have monotonic error attributed to nonlinear coupling; the two-stream run reaches a long-time BGK state.The authors discuss reasons their results lack numerical recurrence and display fine-scale distribution features of the final BGK state.Source: §5 · p. 55
- Higher-order polynomial bases, an improved temporal integrator, a positivity-preserving limiter, and recurrence studies are identified as ongoing or upcoming work.For piecewise-constant functions, recurrence can occur, but its amplitude decays.Source: §5 · p. 55
- Nonuniform-mesh operation has been seen in unreported work as a first step toward an adaptive code; parallel implementation and applications with physical boundary conditions or collision operators remain prospective.The authors characterize the DG method as an attractive alternative for plasma computations, noting that Vlasov-Poisson serves as a test bed for more sophisticated kinetic theories.Source: §5 · pp. 55–56
A Appendix: Dispersion Relations
The appendix analyzes dispersion relations for Lorentzian and two-stream equilibria. It derives the Lorentzian damping behavior and characterizes the two-stream equilibrium as unstable.
Dispersion relation setup
- The analysis uses perturbations proportional to exp(ikx−iωt); for Landau modes, ω = ωR+iγ gives the asymptotic electric-field decay rate γ and oscillation frequency ωR.The stated interpretation applies after analytic continuation into the lower half of the ω-plane.Source: A Appendix: Dispersion Relations · p. 56
- For bounded domains, the real positive wave numbers are k = 2πn/L, with integer n; the appendix numerically searches for roots to produce Fig. 13.For the stated system size L = 4π, the authors evaluate ϵ with k replaced by k/2.Source: A Appendix: Dispersion Relations · pp. 56, 58
Dispersion relation derivation
- The Lorentzian dispersion function is obtained by evaluating its integral using Cauchy's theorem and summing residues at v = i and v = u.For the two-stream case, the plasma Z-function is related to the Hilbert transform and error function, and its lower-half-plane values are obtained by analytic continuation.Source: A Appendix: Dispersion Relations · pp. 57–58
Lorentzian equilibrium
- For the Lorentz equilibrium, analytic continuation gives Landau damping with γ(k) = −k and |ωR(k)| = 1, while there are no discrete modes in the upper half-plane.The passage relates the absence of discrete growing or damped modes to the monotonicity of the equilibrium in v^2.Source: A Appendix: Dispersion Relations · p. 57
Two-stream equilibrium
- For the two-stream equilibrium, the dispersion relation has a root in the upper half-plane, indicating instability.The appendix expresses the relation using the plasma Z-function for computational convenience.Source: A Appendix: Dispersion Relations · p. 57