Boltzmann Equation Solver for Thermalization

This paper introduces BEST, a parallelized Python framework that solves the momentum-resolved Boltzmann equation for arbitrary scattering processes using VEGAS Monte Carlo integration, while addressing a critical subtlety in handling identical particles to ensure exact energy conservation and validate the method through thermalization studies of massive scalar fields.

The Gist

Imagine the universe as a giant, bustling dance floor. In the very early days, particles (the dancers) were moving wildly, bumping into each other, swapping partners, and changing the music's tempo. Over time, they need to settle down into a smooth, synchronized rhythm called "thermal equilibrium."

The Boltzmann Equation is the rulebook that predicts how this dance evolves. It tells us how the crowd's energy and movement change every second. However, calculating this rulebook is incredibly hard, especially when the dancers aren't just swapping partners (2 people $\to$ 2 people) but are doing complex group moves like splitting into three or merging from three (2 $\leftrightarrow$ 3).

This paper introduces Best, a new computer program designed to solve this dance-floor math for any number of dancers, not just pairs. Here is the breakdown in simple terms:

1. The Problem: The "Impossible" Math

For simple pair-swaps (2 $\to$ 2), physicists have had good ways to calculate the outcome for decades. But when particles do complex things—like two particles colliding to create three new ones ($2 \to 3$)—the math explodes.

• The Analogy: Imagine trying to predict the outcome of a game of billiards where, instead of two balls hitting, you have to track the trajectory of 5 balls simultaneously, all while they are changing weight and speed. The number of possibilities becomes so huge that standard calculators crash.
• The Old Way: Scientists used to simplify the problem by ignoring the complex details or assuming the particles were already dancing in a perfect rhythm. But this fails when the universe is chaotic and the particles are "out of sync."

2. The Solution: "Best" (The Super-Computer Dance Coach)

The author, Jong-Hyun Yoon, built Best (Boltzmann Equation Solver for Thermalization). It's a Python tool that acts like a super-observant coach watching every single dancer.

• How it works: Instead of trying to solve the whole dance floor at once, it picks one dancer (a specific momentum) and asks, "What happens to you when you bump into everyone else?"
• The Monte Carlo Method: Since there are too many possibilities to count, Best uses a trick called "Monte Carlo." Imagine throwing thousands of darts at a dartboard to guess the average score. Best throws millions of "virtual collisions" to figure out the average outcome. It uses a smart algorithm (Vegas) that learns where the "action" is happening and throws more darts there.

3. The Big Discovery: The "Missing Step"

The most important part of this paper is a subtle mistake everyone else was making.

• The Scenario: Imagine a reaction where 2 particles turn into 3.
• The Mistake: Previous calculators treated the "2" side and the "3" side as if they were the same. They thought, "If I watch a particle on the '2' side, it's the same as watching one on the '3' side."
• The Reality: They are not the same.
  • The Analogy: Think of a family reunion. If you are one of two parents, your role is different from being one of three children. If you only count the "children" side of the equation, you miss half the story.
• The Consequence: If you ignore this difference, your math says energy is being created or destroyed out of thin air. The paper shows that without fixing this, the simulation loses about 40% of its energy, which is a disaster for physics.
• The Fix: Best automatically counts every possible position a particle could be in (on the "2" side and the "3" side) and adds them up correctly. This ensures energy is perfectly conserved.

4. Why This Matters

This tool allows scientists to study scenarios that were previously impossible to simulate:

• Dark Matter: Some theories suggest dark matter particles eat each other (3 $\to$ 2) or split apart. Best can simulate this to see if it explains what we see in the universe.
• The Early Universe: It helps us understand how the universe cooled down and settled into its current state.
• Mass Changes: It can handle particles that change their "weight" (mass) as the universe expands, like a dancer changing costumes mid-performance.

5. The Power of Parallelism

To make these calculations fast, Best uses a "crowdsourcing" approach.

• The Analogy: If one person has to calculate the dance moves for 1,000 dancers, it takes forever. But if you have 1,000 people, each calculating the moves for just one dancer simultaneously, it's instant.
• Best splits the work across hundreds of computer cores (processors), making it incredibly fast and scalable.

Summary

Best is a new, open-source software that lets physicists simulate the chaotic dance of particles in the early universe with unprecedented accuracy. It fixes a hidden error in how we count particle interactions, ensuring that the laws of physics (like energy conservation) are never broken, even when particles are doing complex group moves. It's like upgrading from a blurry, low-resolution map of the universe to a crystal-clear, high-definition simulation.

Technical Summary

1. Problem Statement

The Boltzmann equation governs the evolution of particle phase-space distributions in cosmological scenarios, such as dark matter freeze-out, freeze-in, and dark sector thermalization. While the equation is conceptually simple, solving the collision integral for momentum-resolved distributions is computationally formidable.

• Dimensionality Challenge: For a general scattering process with $n_{in}$ initial and $n_{out}$ final particles, the collision integral requires integration over a phase space of dimension $3(n_{total} - 2)$ (where $n_{total} = n_{in} + n_{out}$).
  • Standard $2 \to 2$ processes require 6 dimensions (reducible to 2 via analytical methods).
  • Higher-multiplicity processes (e.g., $2 \to 3$, $3 \to 2$, $2 \to 4$) require 9, 12, or more dimensions, rendering analytical reduction impossible.
• Identical Particle Subtlety: For processes where $n_{in} \neq n_{out}$ involving identical particles, a critical subtlety exists in constructing the collision integral. Standard textbook treatments often integrate over momentum to find number densities, absorbing symmetry factors. However, in momentum-resolved calculations where the observed momentum $p$ is fixed, contributions from the "gain" and "loss" terms differ structurally depending on which side of the reaction the observed particle resides. Previous solvers often neglected this distinction for asymmetric processes, leading to systematic violations of energy conservation.
• Limitations of Existing Solvers: Existing codes (e.g., Drake, methods by Ala-Mattinen et al.) rely on analytical reductions valid only for $2 \to 2$ (or $1 \to 2$) processes. They cannot handle arbitrary $n_{in} \to n_{out}$ processes or the specific identical-particle decomposition required for asymmetric reactions.

2. Methodology: The "Best" Framework

The author presents Best (Boltzmann Equation Solver for Thermalization), a Python-based framework designed to solve the momentum-resolved Boltzmann equation for arbitrary scattering processes using direct numerical evaluation.

A. Core Algorithm

• Direct Monte Carlo Integration: Instead of analytical reduction, the collision integral is evaluated directly in $3(n_{total}-2)$ dimensions using the VEGAS adaptive Monte Carlo algorithm.
• Vectorization & Parallelization:
  • Vectorized Batch Evaluation: The integrand is evaluated in batches to maximize CPU/GPU efficiency.
  • MPI Parallelization: The momentum grid points are distributed across MPI ranks. Each rank computes the collision integral for its assigned grid points independently, yielding near-linear scaling to hundreds of cores.
  • Integrator Reuse: The VEGAS adaptive importance sampling map is preserved across time steps and grid points, as the integrand structure varies slowly, significantly reducing the number of evaluations needed per step.

B. Handling Conservation Laws

• Momentum Conservation: Enforced exactly by expressing the momentum of one "conserved" particle as a function of the others: $\vec{p}_{cons} = \sum \vec{p}_{in} - \sum \vec{p}_{out}$.
• Energy Conservation: Imposed via a narrow Gaussian representation of the Dirac delta function:
  $$\delta(E_{in} - E_{out}) \approx \frac{1}{\sigma_E \sqrt{2\pi}} \exp\left(-\frac{(E_{in} - E_{out})^2}{2\sigma_E^2}\right)$$
  The width $\sigma_E$ is adaptive; if the Monte Carlo error is high, the width increases to improve sampling, then narrows to sharpen the constraint.
• Optimization of Conserved Particle: For asymmetric processes ($n_{in} \neq n_{out}$), the conserved particle is placed on the side with fewer particles. This ensures the reconstructed momentum remains physical (on-shell) over a broader region of the integration domain.

C. Identical-Particle Decomposition (Key Innovation)

The paper identifies and solves a critical issue for processes with $n_{in} \neq n_{out}$ and identical particles (e.g., $\phi\phi \leftrightarrow \phi\phi\phi$).

• The Problem: The observed momentum $p$ can appear on the initial side (2-particle side) or the final side (3-particle side). These correspond to distinct integrals, $C_2(p)$ and $C_3(p)$, which are not equal.
• The Solution: The full collision rate must be a weighted sum:
  $$C_{total}(p) = n_{\alpha} C_{n_{\alpha}}(p) + n_{\beta} C_{n_{\beta}}(p)$$
  where $n_{\alpha}$ and $n_{\beta}$ are the multiplicities of the species on each side.
• Consequence: Failing to include both terms (e.g., calculating only $C_3$) leads to a systematic violation of energy conservation (up to ~40% in tested cases), as the gain/loss balance is incorrect.

D. Additional Features

• Quantum Statistics: Supports Bose-Einstein (stimulated emission) and Fermi-Dirac (Pauli blocking) statistics exactly in both gain and loss terms.
• Time-Dependent Masses: Supports massive particles with masses that change over time (crucial for phase transitions).
• Cosmological Expansion: Uses comoving momenta $q = ap$ to handle Hubble expansion, converting to physical momenta $p = q/a(t)$ within the collision integral.
• Interpolation/Extrapolation: Uses cubic splines within the grid and physically motivated exponential extrapolation outside the grid to handle distribution tails.

3. Key Results

The code was validated and demonstrated through several benchmarks:

1. $2 \to 2$ Benchmark (Validation):

   • Compared against a semi-analytical method (exact energy conservation, reduced to 2D) for both massless and massive scalar fields.
   • Result: The Monte Carlo results agreed with the semi-analytical solution within a few percent across the momentum range. Discrepancies only appeared near zero-crossings of the collision rate where Monte Carlo signal-to-noise is low.

2. Thermalization via $2 \to 2$:

   • Simulated a massive scalar field relaxing from a non-thermal distribution to a Bose-Einstein equilibrium.
   • Result: The system correctly conserved particle number and energy ($\sim 1\%$ error) and converged to the expected equilibrium temperature and chemical potential.

3. $2 \leftrightarrow 3$ Number-Changing Process (Critical Demonstration):

   • Simulated the $\phi\phi \leftrightarrow \phi\phi\phi$ process.
   • Incorrect Approach: Using only the term where $p$ is on the 3-particle side ($C_3$) resulted in a ~40% loss of energy over time.
   • Correct Approach: Using the full decomposition ($2C_2 + 3C_3$) restored energy conservation to within $\sim 1\%$ (consistent with statistical noise).
   • Equilibrium: The system correctly thermalized to a Bose-Einstein distribution with chemical potential $\mu = 0$ (required by $3\mu = 2\mu$ for number-changing equilibrium).

4. Performance:

   • Achieved near-linear scaling on the KISTI Nurion supercomputer.
   • A single time step for a $2 \to 2$ process (272 grid points, $10^6$ evaluations) took ~50s on 272 MPI ranks.
   • A $2 \to 3$ process (9 dimensions) took ~300s on 544 ranks.

4. Significance and Impact

• Solving the "High-Dimension" Barrier: Best enables the study of cosmological processes previously considered computationally intractable due to high-dimensional phase spaces (e.g., $3 \to 2$ cannibal dark matter, semi-annihilation, multi-body production).
• Correcting Systematic Errors: The paper highlights a previously unaddressed error in momentum-resolved solvers for asymmetric processes. The proposed decomposition is essential for energy conservation in models involving number-changing interactions.
• Versatility: The framework is general-purpose, supporting arbitrary $n \to m$ processes, multiple coupled species, time-dependent masses, and full quantum statistics without relying on relaxation time approximations or specific distribution ansatzes.
• Open Source: The code is publicly available, providing the community with a robust tool for investigating complex dark sector dynamics, freeze-in/out scenarios, and thermalization across cosmological phase transitions.

In summary, Best represents a significant advancement in computational cosmology, moving beyond the limitations of $2 \to 2$ approximations to provide a rigorous, scalable, and accurate solver for the full momentum-resolved Boltzmann equation.