Moment-preserving particle merging via non-negative least squares

This paper proposes a novel particle merging algorithm for rarefied gas dynamics simulations that conserves arbitrary velocity and spatial moments, as well as collision rates, by solving a non-negative least squares problem, thereby significantly reducing merging-induced errors in macroscopic quantities.

The Gist

Imagine you are trying to simulate a massive crowd of people moving through a city using a computer. In the world of physics, this "crowd" is a gas or plasma made of trillions of atoms. Since computers can't track every single atom, scientists use computational particles. Think of each computational particle as a "representative" or a "proxy" standing in for a huge group of real atoms.

Sometimes, these proxies get too heavy (representing too many atoms) or there are simply too many of them, slowing down the computer. To fix this, scientists need to merge them: take a group of proxies and combine them into fewer, smarter proxies.

The problem? If you just mash them together randomly, you lose important information. It's like taking a diverse group of people (some tall, some short, some fast, some slow) and replacing them with two average-looking people. You've saved space, but you've lost the shape of the crowd. The simulation might think the gas is hotter or colder than it really is, or that it's flowing in the wrong direction.

The Old Way: The "Binning" Method

Previously, scientists used a method called binning. Imagine sorting the crowd into boxes based on their speed.

• The Analogy: You put all the "fast runners" in Box A and all the "slow walkers" in Box B. Inside each box, you grab a few people and say, "Okay, you two are now one person who is the average of the two."
• The Flaw: This is a bit like averaging a marathon runner and a toddler. The result is a "medium-speed" person who doesn't actually exist. You lose the extreme details (the very fast runners and the very slow walkers), which are crucial for predicting things like heat transfer or chemical reactions.

The New Way: The "Moment-Preserving" Method

The authors of this paper, Georgii and Manuel, propose a smarter way to merge these particles using a mathematical tool called Non-Negative Least Squares (NNLS).

Here is the creative analogy for their method:

The "Perfect Portrait" Analogy
Imagine you have a complex painting made of thousands of tiny dots (the particles). You want to repaint it using fewer dots, but you must keep the exact same overall shape, color balance, and shadows (the "moments" or physical properties like energy and momentum).

1. The Goal: You don't just average the dots. Instead, you look at the original painting and ask: "Which specific dots from the original painting, if I keep them and adjust their size (weight), can recreate the exact same picture?"
2. The Math (NNLS): The computer solves a puzzle. It tries to find a combination of the original dots that, when weighted correctly, perfectly matches the mathematical "fingerprint" (the moments) of the original crowd.
3. The Constraint (Non-Negative): You can't have "negative people." In physics, a particle's weight represents a count of atoms, so it must be zero or positive. The math ensures you never create a "negative weight" particle, which would break the laws of physics.

Why This is a Big Deal

1. Keeping the "Tail" of the Distribution
In a gas, most particles move at an average speed, but a few move incredibly fast. These fast particles are the "tail" of the distribution.

• Old Method: The binning method often smears out these fast particles, making them look average.
• New Method: Because the NNLS method picks specific original particles to keep, it preserves those rare, super-fast particles. This is vital because those fast particles are often the ones causing chemical reactions or generating heat.

2. The "Reaction Rate" Upgrade
The authors also added a special feature for plasma (ionized gas). In plasma, particles collide to create new particles (ionization).

• The Problem: If you merge particles incorrectly, you might accidentally change how often they collide. It's like merging a group of aggressive drivers and a group of calm drivers into a single "average" driver, changing the likelihood of a crash.
• The Solution: They derived a way to ensure that even after merging, the rate of these collisions stays exactly the same. They even created a "shortcut" version for electrons (which are super fast) that approximates this perfectly without needing to do heavy math.

The Results: A Smoother Ride

The authors tested their new algorithm on several scenarios:

• Relaxing Gas: When gas settles down to a calm state, the new method stayed much closer to the "true" answer than the old binning method.
• Electric Fields: When simulating plasma in an electric field, the new method predicted the temperature and reaction rates with much less error.
• Heat Flow: In a simulation of heat flowing through a pipe, the new method predicted the temperature profile almost perfectly, while the old method got the walls too hot or too cold.

The Bottom Line

Think of the old method as blending a smoothie: you take all the ingredients, mix them up, and you lose the texture of the individual fruits.

The new method is like curating a playlist: you select the best, most representative songs from the original album to create a shorter playlist that still captures the exact vibe, energy, and emotion of the full album.

By using this "curated" approach (solving a specific math puzzle), the scientists can run simulations with fewer particles, saving computer power, while still getting results that are more accurate than before. It's a win-win for simulating the invisible world of gases and plasmas.

Technical Summary

1. Problem Statement

In rarefied gas dynamics and plasma simulations, stochastic particle methods like Direct Simulation Monte Carlo (DSMC) and Particle-in-Cell (PIC) are essential for modeling non-equilibrium flows. These methods use computational particles with "weights" to represent physical particles. However, specific scenarios necessitate reducing the number of computational particles:

• Large density gradients: To maintain efficiency in high-density regions.
• Axisymmetric simulations: To improve statistics near the rotational axis.
• Variable-weight collisions: Collisions between particles of unequal weights can generate new particles, leading to exponential growth in particle count.
• Resolution enhancement: Some algorithms split particles to resolve specific collisional processes, increasing the count.

The Challenge: Reducing particle count ("merging") inevitably discards information. Traditional merging methods (e.g., binning-based approaches) often only conserve low-order moments (mass, momentum, energy). This leads to significant errors in higher-order macroscopic quantities (stress tensor, heat flux) and distorts the velocity distribution function (VDF), particularly in the high-energy tails. This distortion cascades into errors in collision operator evaluations and reaction rates, degrading simulation accuracy.

2. Methodology

The authors propose a novel Moment-Preserving Particle Merging algorithm based on solving a Non-Negative Least Squares (NNLS) problem.

Core Formulation

The merging process is formulated as finding a sparse weight vector $\mathbf{w}$ for a subset of pre-merge particles such that specific moments are conserved.

• Objective: Find new weights $\mathbf{w} \geq 0$ for a subset of $N$ pre-merge particles (reducing to $M$ particles) such that:
  $$\mathbf{V}\mathbf{w} = \mathbf{m}$$
  Where $\mathbf{V}$ is a matrix containing velocity and spatial monomials of the pre-merge particles, and $\mathbf{m}$ is the vector of target moments to be conserved.
• Sparsity: The constraint that $\mathbf{w}$ must have only $M$ non-zero entries (where $M < N$) naturally arises from the underdetermined nature of the system when the number of moment constraints is less than the number of particles.
• Algorithm: The authors utilize the Lawson-Hanson algorithm to solve the NNLS system. To ensure numerical stability (since Vandermonde matrices are often ill-conditioned), they introduce a scaling strategy:
  1. Shift velocities and positions relative to the center of mass.
  2. Scale rows and columns based on reference variances and norms.
  3. Solve the scaled system and map the solution back to physical weights.

Extensions

• Rate-Preserving Merging: To address the impact of merging on chemical reaction rates, the authors extend the NNLS system to include constraints for exact collision rates. This involves adding rows to the matrix $\mathbf{V}$ and entries to the vector $\mathbf{m}$ representing the collision integrals.
• Approximate Rate-Preserving Merging: For electron-neutral collisions in plasmas (where electron velocities $\gg$ neutral velocities), an approximate rate constraint is derived. This simplifies the calculation by decoupling the electron rate from the specific neutral velocities, significantly reducing computational cost while maintaining accuracy for high-speed electron processes.

3. Key Contributions

1. Generalized Moment Conservation: The method generalizes previous moment-preserving approaches to spatially inhomogeneous problems, conserving arbitrary orders of velocity and spatial moments.
2. NNLS Framework: It establishes a robust framework using NNLS to find sparse weight solutions, guaranteeing non-negative weights and preserving the physical bounds of particle positions and velocities (by reusing pre-merge particle coordinates).
3. Numerical Stability: The introduction of specific scaling techniques (centering and normalization) significantly improves the conditioning of the linear system, making the algorithm stable for high-order moment conservation.
4. Rate Conservation Extensions: The derivation of exact and approximate rate-preserving constraints allows the method to be applied to plasma simulations where reaction rates are critical.

4. Numerical Results

The algorithm was implemented in the open-source code Merzbild.jl and compared against a state-of-the-art adaptive octree binning method.

• Equilibrium Distribution Merging:

  • When merging particles from a Maxwell-Boltzmann distribution, the NNLS method preserved the high-speed tails of the distribution significantly better than octree binning.
  • Trade-off: NNLS produced a wider distribution of particle weights (higher variance) compared to octree, which could increase collision evaluation costs in some scenarios, though this was less pronounced when pre-merge weights were already non-uniform.

• 0D BKW Relaxation (Boltzmann Equation):

  • In the relaxation of a pseudo-Maxwell gas, the NNLS method showed significantly lower bias in the evolution of 4th and 8th order moments compared to octree merging.
  • The octree method exhibited jumps in moment values immediately after merging, whereas NNLS maintained smooth evolution consistent with the analytical solution.

• Ionization in Plasma (Constant Electric Field):

  • Simulations of Ar/Ar+/e- plasma showed that NNLS merging (especially with rate preservation) yielded ionization rate coefficients and electron temperatures closer to reference values (Bolsig+ and high-particle-count DSMC) than octree merging.
  • Immediate Post-Merge Error: In the first 50 timesteps after merging, NNLS demonstrated significantly lower bias in ionization rates and electron temperature compared to octree, indicating faster recovery from merging-induced perturbations.

• 1D Fourier Flow (Inhomogeneous):

  • In a temperature-gradient flow between two plates, NNLS merging produced temperature and number density profiles much closer to a fixed-weight reference solution.
  • Octree merging overestimated temperatures and showed large deviations in number density near the walls.
  • Surface Quantities: NNLS reduced the error in incident number density flux, surface pressure, and kinetic energy flux by a factor of 10 or more compared to octree at low particle counts (e.g., kinetic energy flux error dropped from ~10% to <1%).

5. Significance

This work provides a mathematically rigorous and numerically stable approach to particle reduction in kinetic simulations. By leveraging NNLS, the authors demonstrate that conserving higher-order moments is not just a theoretical improvement but a practical necessity for accuracy in:

• Rarefied Gas Dynamics: Accurately capturing heat flux and stress tensors in non-equilibrium flows.
• Plasma Physics: Correctly modeling ionization rates and electron temperatures, which are highly sensitive to the high-energy tail of the distribution function.

The proposed method allows for simulations with lower average particle counts while maintaining accuracy comparable to high-resolution simulations, potentially offering significant computational savings for complex, multi-scale problems. The extension to rate preservation makes it particularly valuable for reactive plasma simulations where standard merging often fails to capture correct chemistry.