Symplectic particle-in-cell methods for hybrid plasma models with Boltzmann electrons and space-charge effects

This paper develops and validates geometric particle-in-cell methods that preserve the Hamiltonian structure and energy for hybrid plasma models featuring kinetic ions, Boltzmann electrons, and space-charge effects, demonstrating their effectiveness through numerical experiments on instabilities and wave phenomena.

The Gist

Imagine a giant, chaotic dance floor filled with two types of dancers: Heavy Ion Dancers and Light Electron Dancers.

In the world of plasma physics (the study of super-hot, charged gases like the sun or fusion reactors), scientists want to predict how these dancers move. Usually, you have to track every single dancer individually. But there are so many electrons (billions of billions) that tracking them all is like trying to count every grain of sand on a beach while running a marathon—it takes too long and crashes your computer.

The Problem: The "Heavy" vs. "Light" Dilemma
In this paper, the author, Yingzhe Li, tackles a specific scenario called a Hybrid Model.

• The Ions: These are the heavy, slow dancers. We track them individually because their specific moves matter a lot.
• The Electrons: These are the tiny, hyper-fast dancers. Instead of tracking each one, we treat them like a fluid or a "cloud" that instantly adjusts its shape based on where the ions are. This is the "Boltzmann" part.

However, there's a catch. Even though electrons are fast, they still push and pull on the ions (like static electricity). If you ignore this push-and-pull (called "space-charge effects"), your simulation becomes inaccurate. But if you try to calculate it perfectly, the math gets messy, and standard computer methods often create "ghosts"—fake energy or instability that makes the simulation explode or behave weirdly.

The Solution: A "Geometric" Dance Guide
The author proposes a new way to simulate this dance floor using Symplectic Particle-in-Cell (PIC) methods.

Here is the analogy:
Imagine you are choreographing this dance.

• Old Methods: You take a snapshot of the dancers, guess their next move, take another snapshot, and guess again. Over time, your guesses get slightly wrong. The dancers might slowly gain energy they didn't have, or lose energy they should have. Eventually, the dance floor looks nothing like reality.
• The New Method (Symplectic): Instead of just guessing, you use a special set of rules (a "Geometric" map) that respects the fundamental laws of physics. It's like having a dance guide who knows that energy cannot be created or destroyed. No matter how long the dance goes on, the total energy of the system stays exactly the same (or very, very close to it).

How It Works (The Magic Tricks)
The author uses two main "magic tricks" to keep the simulation stable:

1. Hamiltonian Splitting (The "Step-by-Step" Approach):
   Imagine the dance has two parts: moving forward and reacting to the music.

   • Step 1: Move all the heavy dancers forward based on their current speed.
   • Step 2: Stop them, calculate how the electron cloud pushes them, and update their speed.
   • Repeat: Do this in tiny, alternating steps. By splitting the complex math into simple, solvable chunks, the computer doesn't get confused, and the "ghosts" (instabilities) disappear.

2. Discrete Gradient (The "Energy Bank"):
   This is a stricter rule. It ensures that if the dancers gain kinetic energy (speed), the potential energy (position) drops by the exact same amount. It's like a perfect bank account where money never disappears or appears out of thin air. This method is slower to compute but guarantees the energy is 100% conserved.

Why Does This Matter?
The paper tests these methods on three tricky scenarios:

1. The Grid Instability: When the dance floor is too coarse (big tiles), standard methods make the dancers vibrate wildly and break the simulation. The new methods smooth this out, even with fewer tiles.
2. Landau Damping: This is when a wave in the plasma slowly dies out as energy transfers to the dancers. The new methods predict exactly how fast the wave dies, matching real-world physics perfectly.
3. Resonant Waves: When you push the dancers at just the right rhythm, huge waves form. The new methods capture these complex, non-linear waves without the simulation crashing.

The Bottom Line
Yingzhe Li has built a better "dance guide" for plasma simulations. By respecting the underlying geometry of the universe (conservation of energy and momentum), these new methods allow scientists to run longer, more accurate simulations of fusion reactors and space plasmas without the computer crashing or the results becoming nonsense. It's a way to make the chaotic dance of the universe predictable and stable.

Technical Summary

1. Problem Statement

The paper addresses the numerical simulation of electrostatic hybrid plasma models where ions are treated kinetically (via the Vlasov equation) while electrons are modeled as a Boltzmann fluid. A critical feature of this model is the inclusion of space-charge effects via the Poisson equation, leading to a nonlinear Poisson–Boltzmann equation.

• Challenges:
  • Standard Particle-in-Cell (PIC) methods often suffer from finite grid instability, particularly when the electron-to-ion temperature ratio ($T_e/T_i$) is large.
  • Standard methods often fail to preserve the underlying geometric structures (symplecticity, momentum conservation) and energy conservation over long-time simulations.
  • Existing structure-preserving methods often rely on the quasi-neutral limit (neglecting space-charge effects), which is insufficient for regimes where the Debye length is resolved (e.g., inertial confinement fusion).
  • There is a need for robust numerical schemes that can handle the nonlinearity of the Poisson–Boltzmann equation while maintaining stability and accuracy.

2. Methodology

The author proposes geometric structure-preserving discretizations for the Hybrid Boltzmann-Space-charge (HBS) model. The methodology follows a "variational" or "bracket" approach to ensure the discrete system inherits the properties of the continuous physical system.

A. Mathematical Formulation

• Model: The HBS model consists of the Vlasov equation for ions, $E = -\nabla \phi$, and the Poisson–Boltzmann equation:
  $$-\Delta \phi = Z \int f dv - n_0 \exp\left(\frac{\phi - \phi_0}{T_e}\right)$$
• Hamiltonian Structure: The model is formulated as a Hamiltonian system using:
  1. Action Integral: Derived by combining Low's action principle with the Boltzmann electron term.
  2. Poisson Bracket: A canonical bracket for the distribution function $f$ and a bracket for the potential $\phi$.
• Discretization Strategy:
  • Phase Space (Particles): The distribution function $f$ is discretized using a Particle-in-Cell (PIC) approach with $N_p$ macro-particles and shape functions (B-splines).
  • Field (Potential): The electrostatic potential $\phi$ is discretized using Finite Difference Methods (FDM) on a uniform grid.
  • Equivalence: The paper proves that discretizing the Action Integral and discretizing the Poisson Bracket yield equivalent finite-dimensional Hamiltonian systems.

B. Time Integration Schemes

Two specific time-stepping methods are employed to preserve different geometric properties:

1. Hamiltonian Splitting Method (Symplectic):
   • The Hamiltonian is split into kinetic ($H_1$) and potential/field ($H_2$) parts.
   • Subsystems are solved exactly (or via fixed-point iteration for the nonlinear Poisson–Boltzmann equation).
   • Composed using Lie splitting (1st order) or Strang splitting (2nd order).
   • Property: Preserves the symplectic structure and exhibits excellent long-term stability, though energy is conserved only approximately (bounded error).
2. Discrete Gradient Method (Energy Conserving):
   • An implicit method using a discrete gradient operator.
   • Property: Preserves the total energy exactly (up to machine precision) by construction.

C. Electromagnetic Extension

The framework is extended to an electromagnetic hybrid model (involving a high-frequency laser field) proposed by Vu.

• A new Poisson bracket is constructed as the sum of the HBS bracket and the canonical bracket of the Schrödinger equation (for the vector potential envelope).
• A corresponding Hamiltonian splitting scheme is derived for this coupled system.

3. Key Contributions

1. Structure-Preserving HBS Scheme: Development of the first symplectic and energy-conserving PIC methods specifically for the HBS model with full space-charge effects (non-quasi-neutral).
2. Proof of Equivalence: Rigorous proof that the discretization via Action Integral and Poisson Bracket leads to the same discrete Hamiltonian system.
3. Conservation Laws: Demonstration that the global neutrality condition is conserved under periodic or zero Neumann boundary conditions.
4. Asymptotic Preserving: The schemes are shown to be asymptotic preserving, meaning they correctly reduce to the structure-preserving schemes for the quasi-neutral limit ($\lambda_D \to 0$) and the large electron temperature limit.
5. Suppression of Instabilities: The geometric nature of the methods is shown to effectively suppress finite grid instability, a common failure mode in standard PIC codes for hybrid models.

4. Numerical Results

The methods were validated using the Python package STRUPHY through three benchmark experiments:

• Finite Grid Instability:

  • Setup: Equilibrium simulation with high $T_e/T_i$.
  • Result: Both splitting and discrete gradient methods successfully suppressed the rapid linear growth of ion kinetic energy seen in standard PIC. The ion kinetic energy oscillated stably around the initial value.
  • Accuracy: Momentum errors were kept low ($\sim 10^{-4}$), and energy errors were minimal ($\sim 10^{-5}$ for splitting, $\sim 10^{-13}$ for discrete gradient).

• Landau Damping:

  • Linear: Simulated linear ion Landau damping. The numerical damping rates matched the theoretical dispersion relation derived from the HBS model.
  • Nonlinear: Simulated nonlinear Landau damping with large $T_e$. The methods correctly captured the initial exponential decay, the oscillatory phase, and the subsequent exponential growth (instability), matching known Vlasov-Poisson results.
  • Energy Conservation: The discrete gradient method maintained total energy with errors around $10^{-12}$, while the splitting method maintained bounded errors around $10^{-2}$ to $10^{-4}$ depending on parameters.

• Resonantly Excited Nonlinear Ion Waves:

  • Setup: Simulation with a time-dependent ponderomotive driving term ($\phi_0$).
  • Result: The methods reproduced the resonant response peak ($R(t) \approx 5$) and the correct frequency modulation (fast and slow modes) observed in previous literature.
  • Phase Space: Contour plots showed the formation of 5 vortices, consistent with the driving mode $k=5$.

5. Significance

• Long-Term Reliability: By preserving geometric structures (symplecticity and energy), these methods are superior for long-time plasma simulations where standard methods drift or become unstable.
• Physical Fidelity: The ability to resolve space-charge effects without sacrificing stability allows for accurate modeling of regimes (like inertial confinement fusion) where the Debye length is critical.
• Robustness: The suppression of finite grid instability makes these methods highly suitable for high-temperature plasma simulations where traditional PIC codes often fail.
• Extensibility: The framework provides a blueprint for extending structure-preserving methods to more complex electromagnetic hybrid models involving laser-plasma interactions.

In conclusion, this work provides a rigorous mathematical foundation and practical numerical algorithms for simulating hybrid plasmas with space-charge effects, offering a significant improvement in stability and conservation properties over existing non-geometric approaches.