Loading non-Maxwellian Velocity Distributions in Particle Simulations

This paper presents a comprehensive set of numerical procedures and algorithms for generating various non-Maxwellian velocity distributions, including $(r,q)$, kappa, ring, shell, super-Gaussian, and filled-shell distributions, to be used in particle simulations.

The Gist

Imagine you are a chef trying to simulate a bustling kitchen (a plasma) on a computer. In a perfect, calm kitchen, the chefs (particles) move at a predictable, average speed. This is what physicists call a Maxwellian distribution—it's like a bell curve where most people are average, and very few are extremely fast or extremely slow.

But real space kitchens are chaotic! In the solar wind or around planets, particles don't behave politely. Some are stuck in a "loss cone" (escaping like chefs running out the back door), some form rings like dancers in a circle, and others have "power-law tails" where a few super-fast chefs zoom around, breaking the rules. These are non-Maxwellian distributions.

The problem? Computer simulations (Particle-in-Cell or PIC) are great at handling the "average" chefs, but they struggle to generate these weird, chaotic crowds. If you try to force a square peg into a round hole, your simulation breaks or runs incredibly slowly.

This paper is essentially a cookbook of "numerical recipes" for nine different types of chaotic particle crowds. The authors, Zenitani, Usami, and Matsukiyo, are giving simulation scientists the specific tools to generate these weird crowds accurately and efficiently.

Here is a breakdown of their "menu" using simple analogies:

1. The "Flattop" and "Kappa" Distributions (The Generalized Crowd)

• The Problem: Sometimes the crowd isn't a smooth bell curve. It might have a flat top (everyone is moving at roughly the same speed) or a long tail of super-fast runners (the Kappa distribution).
• The Solution: The authors offer two ways to generate these.
  • The "Beta-Prime" Method: Think of this as mixing two specific ingredients (random numbers) in a blender to get the perfect texture. It's fast and precise for most cases.
  • The "Piecewise Rejection" Method: Imagine you are trying to pick apples from a tree. You have a basket (the envelope). You pick an apple (a random number), and if it fits the shape of your target crowd, you keep it. If it's too big or too small, you throw it back and try again. The authors figured out how to build the perfect basket shape so you don't waste time throwing apples away.

2. The "Regularized Kappa" (The Capped Runner)

• The Problem: In some space environments, there are no infinite speeds. There is a hard ceiling (a cutoff) where particles stop. Standard math breaks down here because it expects infinite energy.
• The Solution: They created a "Post-Rejection" method. First, they generate a standard crowd, and then they act like a bouncer at a club: "If you are faster than the cutoff speed, you're out!" They also offer a "Piecewise" method that builds the crowd from the ground up to respect the speed limit, ensuring the math never crashes.

3. The "Subtracted Kappa" (The Empty Cone)

• The Problem: Imagine a cone-shaped hole in the middle of a crowd where no one is allowed to stand (a "loss cone"). This happens when particles escape along magnetic field lines.
• The Solution: The authors propose a clever trick. Instead of trying to build the hole directly, they build two crowds and subtract one from the other. It's like taking a full pie and cutting out a slice to leave a gap. This allows them to create a very narrow, precise gap in the particle speeds, which is crucial for studying Earth's magnetosphere.

4. Rings and Shells (The Dance Floor vs. The Bubble)

• The Problem: Sometimes particles form a ring (like a hula hoop) or a shell (like a hollow ball) in speed space.
• The Solution:
  • The "Gaussian" Way: They treat the ring/shell as a fuzzy cloud. You can generate this by taking a normal crowd and stretching it into a ring or shell shape.
  • The "Maxwellian" Way (The New Favorite): The authors suggest a simpler alternative. Instead of trying to force a perfect ring, imagine taking a normal crowd and spinning it around a center point (for a ring) or scattering it in all directions from a moving source (for a shell).
  • Why it's better: The "Maxwellian" versions are mathematically easier to handle. They don't have "artificial edges" that cause computer glitches. It's like using a smooth, spinning top instead of trying to balance a stack of jagged rocks.

5. Super-Gaussian and Filled-Shells (The Extreme and The Full)

• The Problem: Some distributions are "sharper" than a bell curve (Super-Gaussian), and others are just a solid block of particles up to a certain speed (Filled-Shell).
• The Solution: They provide direct formulas. For the filled-shell, it's as simple as rolling a die to pick a speed, ensuring no one is faster than the limit. For the Super-Gaussian, they use a special "gamma" ingredient to sharpen the peak.

Why Does This Matter?

Before this paper, if a scientist wanted to simulate a specific type of weird particle behavior, they might have to invent a new, complex math trick from scratch, or their computer would run so slowly it would take years to get an answer.

This paper says: "Don't reinvent the wheel. Here are the blueprints for nine different wheels."

By providing these "numerical recipes" (algorithms), the authors make it much easier for scientists to model:

• How solar wind hits Earth's magnetic shield.
• How particles get accelerated in space storms.
• How plasma behaves in fusion reactors.

In short, they turned a difficult, messy math problem into a set of clear, copy-paste instructions, making the simulation of space plasma faster, more accurate, and accessible to everyone.

Technical Summary

1. Problem Statement

In heliophysics and plasma physics, Particle-in-Cell (PIC) simulations are essential for modeling kinetic processes such as wave-particle interactions and particle acceleration. While generating Maxwellian (thermal) velocity distributions is straightforward using standard algorithms (e.g., Box-Muller), generating non-Maxwellian Velocity Distribution Functions (VDFs) is significantly more challenging.

Many observed plasma environments (solar wind, magnetospheres, bow shocks) exhibit complex VDFs, including:

• Power-law tails (Kappa distributions).
• Flattop distributions (common in magnetospheres).
• Loss-cone distributions (due to particle escape along magnetic field lines).
• Ring and shell distributions (from pickup ions).
• Super-Gaussian and filled-shell distributions.

General-purpose sampling methods like Acceptance-Rejection or Inverse Transform Sampling often suffer from inefficiency (low acceptance rates in 3D) or lack of analytical inverses for complex distributions. Furthermore, many specific non-Maxwellian models lack dedicated, robust numerical recipes for implementation in simulation codes.

2. Methodology

The authors propose a comprehensive set of "numerical recipes" to generate nine distinct non-Maxwellian VDFs. The methodology relies on transforming basic random number generators (Uniform, Normal, and Gamma variates) into the target distributions.

The core approach involves:

• Analytical Derivation: Expressing target VDFs as transformations of known statistical distributions (e.g., Gamma, Beta-prime, Student's t).
• Rejection Sampling: Developing piecewise envelope functions to ensure high acceptance efficiency for distributions that cannot be directly inverted.
• Hybrid Methods: Combining standard distributions with specific rejection steps to handle features like loss cones or high-energy cutoffs.

The paper provides explicit algorithms (presented as pseudocode in Tables I–VII) for the following distributions:

A. Generalized (r, q) Distribution

• Context: Generalizes Kappa and flattop distributions.
• Methods:
  1. Beta-prime Method: Uses the ratio of two Gamma variates. Efficient for Kappa-like regimes but risks division by zero if shape parameters $< 1$.
  2. Piecewise Rejection: Splits the distribution into low-velocity and high-velocity envelopes. Recommended for flattop-like regimes or when the Beta-prime method fails.

B. Regularized Kappa Distribution (RKD)

• Context: Kappa distribution with a high-energy cutoff to prevent divergence of higher-order moments (allowing $\kappa < 3/2$).
• Methods:
  1. Post-Rejection: Generate a standard Kappa (Student's t) distribution, then reject particles based on an exponential cutoff factor. Efficient for $\kappa > 1/2$.
  2. Piecewise Rejection: Uses inverse transform on piecewise envelope functions. Works for the entire range $\kappa \ge 0$, including the difficult $\kappa \le 1/2$ regime.

C. Subtracted Kappa Distribution

• Context: A recent model for narrow loss-cone distributions, extending the subtracted Maxwellian.
• Method: A compact algorithm using a weighted sum of exponential variates to model the loss-cone filling parameter ($\Delta$) and opening angle ($\beta$). It avoids complex power-law pitch-angle dependencies.

D. Ring and Shell Distributions (Gaussian Width)

• Context: Models for pickup ions (PUIs) with finite thermal width.
• Methods:
  1. Inverse Transform: Numerical inversion of the Cumulative Distribution Function (CDF).
  2. Piecewise Rejection: Uses log-concave properties to define exponential envelopes, offering high efficiency ($\sim 86-89\%$).

E. Ring and Shell Maxwellians

• Context: Proposed as simpler, physically motivated alternatives to the Gaussian-width distributions. They represent the scattering of a seed Maxwellian population (axisymmetric for rings, spherical for shells).
• Method: Generate a standard Maxwellian and rotate the velocity vector. This is computationally simpler and avoids artificial boundary issues at $v=0$ found in conventional ring/shell distributions.

F. Super-Gaussian and Filled-Shell Distributions

• Context: Super-Gaussian (generalized exponential power) and power-law filled shells.
• Method:
  • Super-Gaussian: Transforms a Gamma variate ($v \propto X^{1/p}$).
  • Filled-Shell: Uses a simple power-law transformation of a uniform variate ($v \propto U^{1/(3+p)}$).

3. Key Contributions

1. Comprehensive Algorithmic Library: The paper provides ready-to-implement algorithms for nine distinct non-Maxwellian VDFs, filling a gap in the literature where such recipes were previously scattered or unknown.
2. Efficiency Optimization: The authors rigorously compare different sampling strategies (e.g., Beta-prime vs. Rejection for (r,q); Post-rejection vs. Piecewise for RKD) and define the optimal parameter space for each method to maximize acceptance efficiency.
3. Introduction of Maxwellian Alternatives: The proposal of "Ring Maxwellian" and "Shell Maxwellian" distributions offers a physically robust and computationally cheaper alternative to traditional ring/shell distributions, particularly when the drift velocity is comparable to the thermal width.
4. Handling Singularities: The methods address numerical singularities (e.g., division by zero in Gamma variates with shape $<1$) that often plague standard implementations.
5. Validation: All algorithms are validated against analytical solutions using $10^6$ particles and Kullback-Leibler (KL) divergence tests, confirming high numerical accuracy.

4. Results

• Accuracy: Monte Carlo simulations for all nine distributions show excellent agreement with theoretical phase-space densities. The KL divergence approaches zero as the particle number increases, confirming the statistical validity of the methods.
• Efficiency:
  • The Piecewise Rejection method for (r,q) and RKD distributions achieves high acceptance rates (often $>95\%$ for specific parameters), significantly outperforming naive rejection methods in 3D.
  • The Ring/Shell Maxwellians are shown to be mathematically simpler to implement and calculate moments for (pressure, energy) compared to their Gaussian-width counterparts.
• Parameter Sensitivity: The study identifies specific regimes where certain methods fail (e.g., Beta-prime method failing for low $\kappa$ in RKD) and provides the correct alternative (Piecewise Rejection).

5. Significance

• Enabling Advanced Kinetic Modeling: By providing robust, efficient, and validated numerical recipes, this work lowers the barrier for simulating complex plasma phenomena (e.g., magnetic reconnection with pickup ions, electron heating in magnetotails) that rely on non-Maxwellian VDFs.
• Standardization: The paper standardizes the generation of these distributions, reducing the risk of errors in custom implementations by simulation scientists.
• GPU Considerations: The authors note that while the algorithms are efficient on CPUs, rejection-based methods may face challenges on GPUs due to thread divergence. They suggest that the proposed methods (especially the non-rejection ones like Beta-prime or Maxwellian rotation) are better suited for parallel architectures.
• Future-Proofing: The inclusion of very recent models (Subtracted Kappa, Regularized Kappa) ensures that the community can immediately utilize the latest theoretical developments in observational space physics.

In summary, this paper serves as a definitive "cookbook" for loading non-Maxwellian velocity distributions in particle simulations, bridging the gap between theoretical plasma models and practical numerical implementation.