Loading of Relativistic Maxwellian-type Distributions Revisited

This paper proposes a simple numerical method based on inverse transform sampling to efficiently load relativistic Maxwellian-type distributions by approximating the cumulative distribution of a shifted-Maxwellian energy distribution with an invertible function to generate random energy and momentum variates.

The Gist

Imagine you are a video game developer trying to simulate a crowd of people running through a city. To make the simulation look real, you need to tell every single "person" (or particle) in your computer how fast they should run and in which direction.

In the world of physics, specifically when dealing with particles moving near the speed of light (relativistic particles), this is a notoriously difficult puzzle. This paper, written by Takayuki Umeda, proposes a clever new way to solve that puzzle.

Here is the breakdown of the problem and the solution, using simple analogies.

The Problem: The "Unfair Dice"

In physics simulations (like Particle-in-Cell or Monte Carlo), scientists need to generate random numbers that follow a specific pattern called a Maxwell-Jüttner distribution. Think of this distribution as a very specific, complex rulebook for how fast particles should move.

For a long time, the standard way to generate these numbers was Rejection Sampling.

• The Analogy: Imagine you are trying to fill a bucket with water using a leaky cup. You scoop up water (a random number), check if it fits the bucket's shape (the distribution rule), and if it doesn't fit, you throw it away and try again.
• The Issue: If the bucket has a weird shape, you might throw away 90% of your water. This is incredibly slow and wastes computer power, which is bad for high-performance supercomputers.

The Solution: The "Magic Map"

The author suggests a different method called Inverse Transform Sampling.

• The Analogy: Instead of throwing water away, imagine you have a perfectly straight, clear slide (a uniform random number generator). You want to slide down into a curved, winding pool (the complex distribution).
• The Challenge: To slide down the right path, you need a map that tells you exactly where to start on the slide to land in the right spot in the pool. Mathematically, this map is called the "inverse cumulative distribution function."
• The Problem with the Old Map: For relativistic particles, this map is so mathematically complex that it's impossible to write down as a simple formula. It's like trying to draw a map of a cave that keeps changing shape.

The Paper's Innovation: The "Good Enough" Approximation

Umeda's paper introduces a new way to handle this. Instead of using the old, messy "Maxwell-Jüttner" rulebook, he suggests using a slightly different, but very similar, rulebook called the Relativistic Maxwellian Energy Distribution.

1. Simplifying the Shape: He shows that this new rulebook is mathematically much friendlier. It's like realizing that while the cave is complex, the entrance is actually a perfect circle.
2. The Magic Approximation: He creates a simple, easy-to-calculate formula (a "Magic Map") that approximates the complex curve with incredible accuracy.
   • He tested this map and found the error is less than 0.01% (less than one part in ten thousand).
   • Because the formula is simple, the computer can instantly calculate the inverse (the slide path) without needing to look up huge tables or guess.

How It Works in Practice

The method takes three simple random numbers (like rolling three dice) and turns them into a particle's speed and direction:

1. Die 1: Determines the energy (how fast the particle is moving).
2. Die 2 & 3: Determine the direction (which way it's pointing).

Because the math is so clean, the computer doesn't have to "reject" any numbers. Every single roll of the dice results in a valid particle. This makes the simulation run much faster and more efficiently.

Why Does This Matter?

• Speed: It's like switching from a leaky cup to a firehose. Simulations that used to take hours might now take minutes.
• Accuracy: The author proves that even though he used a "simplified" rulebook, the results are almost identical to the complex, standard rulebook.
• Versatility: This method works even if the particles are drifting in a specific direction (like a wind blowing through the crowd) or if they are moving at different temperatures.

The Bottom Line

Takayuki Umeda found a way to bypass the "leaky cup" problem of simulating fast-moving particles. By swapping a mathematically impossible map for a very accurate, easy-to-use approximation, he gave scientists a faster, cleaner tool to simulate the universe at the speed of light.

In short: He replaced a broken, slow method with a streamlined, high-speed elevator that takes particles exactly where they need to go, every single time.

Technical Summary

Here is a detailed technical summary of the paper "Loading of relativistic Maxwellian-type distributions revisited" by Takayuki Umeda.

1. Problem Statement

In particle kinetic simulations (e.g., Monte Carlo and Particle-in-Cell methods), initializing particle velocities requires generating random variates from specific distribution functions.

• Current Limitations: The standard distribution for relativistic particles is the Maxwell-Jüttner distribution. However, generating random variates from this distribution is computationally challenging because:
  • Rejection Sampling: The most common method (rejection sampling) is often inefficient for high-performance computing (HPC) due to low acceptance rates, especially if the proposal distribution is not optimally chosen.
  • Inverse Transform Sampling: While more efficient, this method requires the analytical inversion of the cumulative distribution function (CDF). For the Maxwell-Jüttner distribution, the CDF cannot be expressed analytically, necessitating the use of numerical tables and interpolation, which introduces accuracy issues and computational overhead.
• Goal: The author seeks a simple, efficient, and accurate numerical method to load relativistic particles (including drifting/shifted distributions) using inverse transform sampling without relying on rejection sampling or complex numerical tables.

2. Methodology

The paper proposes a new approach based on a Relativistic Maxwellian Energy Distribution rather than the traditional Maxwell-Jüttner distribution. The methodology involves the following steps:

A. Theoretical Formulation

• Alternative Distribution: The author introduces a relativistic shifted-Maxwellian energy distribution. This is derived by transforming the non-relativistic Maxwellian energy distribution ($f_K(K)$) into the relativistic regime by substituting kinetic energy $K$ with $mc^2(\gamma - 1)$.
• Shifted Distribution: For drifting plasmas, the distribution is boosted using a drift velocity vector $\mathbf{v}_D$ and its associated Lorentz factor $\gamma_D$. The distribution is formulated in terms of a "boosted Lorentz factor" ($\gamma_B$) and polar angles.
• Normalization: The distribution is normalized such that the integral over the momentum space equals unity.

B. Inverse Transform Sampling Implementation

To enable inverse transform sampling, the author addresses the non-invertibility of the cumulative distribution function (CDF):

1. Dimensionless Energy: The distribution is rewritten in terms of a normalized energy variable $E = \frac{mc^2}{\gamma_D T}(\gamma_B - 1)$.
2. Analytical CDF: The CDF of this normalized energy, $F_E(x)$, is derived analytically involving the error function ($\text{erf}$) and the incomplete Gamma function.
3. Approximation: Since $F_E(x)$ has no closed-form inverse, the author proposes a rational function approximation, $F_{app}(x)$, which is analytically invertible.
   • The approximation takes the form: $F_{app}(x) \approx \left[ 1 - \exp\left( -\frac{ax + bx^2}{1 + cx + dx^2} \right) \right]^{3/2}$.
   • Coefficients ($a, b, c, d$) were determined via a brute-force search to minimize the root mean square error against the true CDF for $0 < x \leq 8$. The relative error of this approximation is found to be less than$10^{-4}$.
4. Coordinate Transformation:
   • Energy/Momentum Magnitude: Uniform random variables are transformed into energy variates using the inverse of the approximated CDF.
   • Angles: The polar angle ($\theta$) and azimuthal angle ($\phi$) are generated. Crucially, for drifting distributions ($v_D \neq 0$), the polar angle distribution depends on the boosted Lorentz factor. The author derives a specific inverse CDF for the polar angle to ensure the resulting momentum vector follows the correct shifted distribution, correcting the common error of using simple isotropic angle generation for drifting frames.
   • Momentum Vector: The scalar variates ($\gamma_B, \theta, \phi$) are converted to Cartesian momentum components ($u_x, u_y, u_z$) using Jacobian transformations.

3. Key Contributions

• New Distribution Function: Introduction of the "Relativistic Maxwellian Energy Distribution" as a viable, analytically tractable alternative to the Maxwell-Jüttner distribution for simulation initialization.
• Analytic Inversion: Development of a highly accurate, invertible approximation for the CDF of the relativistic energy distribution, eliminating the need for numerical interpolation tables.
• Corrected Angular Sampling: Derivation of the correct inverse CDF for the polar angle in drifting frames, ensuring the generated momentum vectors strictly adhere to the relativistic shifted distribution (a step often overlooked in simplified methods).
• Algorithmic Efficiency: The method reduces the loading process to a simple function call: $(u_x, u_y, u_z) = \text{func}(R^{(1)}, R^{(2)}, R^{(3)})$, where $R$ are uniform random numbers. This is highly suitable for vectorization and HPC environments.

4. Results

• Validation of Energy Distribution: Numerical tests with $10^8$ samples show excellent agreement between the generated histogram and the theoretical Maxwellian energy distribution (Eq. 16) on both linear and logarithmic scales.
• Validation of Momentum/Velocity: Comparisons of the reduced momentum and velocity distributions (along the drift axis) against analytical profiles confirm the method's accuracy.
• Comparison with Maxwell-Jüttner:
  • The proposed distribution closely matches the Maxwell-Jüttner distribution in the non-relativistic limit (low temperature/high $mc^2/T$).
  • At high temperatures, the Maxwell-Jüttner distribution exhibits a more enhanced high-energy tail.
  • Physical Properties: The paper calculates and compares mean momentum, thermal energy, and kinetic energy.
    • The mean velocity for both distributions equals the drift velocity $\mathbf{v}_D$.
    • The mean momentum of the proposed Maxwellian distribution is slightly smaller than that of the Maxwell-Jüttner distribution.
    • The thermal energy of the proposed distribution is exactly $3T/2$(in the rest frame), whereas the Maxwell-Jüttner thermal energy includes relativistic corrections dependent on$mc^2/T$.
• Accuracy: The relative error between the approximate CDF and the true CDF is consistently below $10^{-4}$.

5. Significance

• Computational Efficiency: By replacing rejection sampling and numerical table lookups with a direct analytical inversion, the method significantly reduces computational cost and memory usage, making it ideal for large-scale, high-performance plasma simulations.
• Simplicity and Robustness: The algorithm is straightforward to implement as a standard library function, reducing the complexity of simulation codes.
• Versatility: The method is extended to handle anisotropic distributions and non-relativistic limits, providing a unified framework for loading various particle distributions.
• Open Source: The author provides a sample MATLAB code, facilitating immediate adoption by the scientific community.

In conclusion, Umeda's work offers a rigorous, efficient, and accurate solution to a long-standing problem in relativistic kinetic simulations, enabling faster initialization of particle ensembles without sacrificing physical fidelity.