Boltzmann Equation Solver for Thermalization

Abstract: We present BEST (Boltzmann Equation Solver for Thermalization), a Python framework for solving the momentum-resolved Boltzmann equation for arbitrary $n_{\rm in} \to n_{\rm out}$ scattering processes. The collision integral is evaluated directly in $3(n_{\rm total}-2)$ dimensions using the VEGAS adaptive Monte Carlo algorithm with vectorized batch evaluation. Momentum conservation is enforced exactly by expressing one particle's momentum through the constraint, while energy conservation is imposed via a narrow Gaussian representation of the delta function. We identify a subtlety in the construction of the collision integral for processes with unequal initial and final multiplicities ($n_{\rm in} \neq n_{\rm out}$) involving identical particles: the full collision rate requires separate evaluation with the observed momentum pinned to each side of the reaction, weighted by the respective particle multiplicities. Failure to account for this leads to systematic violation of energy conservation. The code supports massive particles with time-dependent masses, Bose-Einstein and Fermi-Dirac quantum statistics, multiple coupled species, cosmological expansion with comoving momenta, and both Euler and Heun time integration. Parallelization is achieved by distributing independent momentum grid points across MPI ranks, yielding near-linear scaling to hundreds of cores. We validate the Monte Carlo results against a semi-analytical $2 \to 2$ collision integral with exact energy conservation, following the phase-space reduction of Ala-Mattinen et al. As a demonstration, we study thermalization of a massive scalar field through a $2 \leftrightarrow 3$ number-changing process and show that energy conservation is restored only when all identical-particle contributions are correctly summed. The code is publicly available at https://github.com/best-hep/best.

• The paper introduces Best, a Python-based framework that implements a momentum-resolved solver for general n-to-m scatterings, addressing complex phase-space integrals.
• It employs adaptive Monte Carlo integration and a rigorous decomposition of collision integrals to maintain energy-momentum conservation, crucial for accurate thermalization.
• The framework demonstrates robust thermal evolution and near-linear scalability, enabling detailed studies of non-equilibrium processes in cosmological and dark sector models.

Boltzmann Equation Solver for Thermalization: Algorithmic Framework and Implications

Introduction and Motivation

The momentum-resolved Boltzmann equation is central to modeling kinetically out-of-equilibrium dynamics in high-energy and cosmological environments, such as the early universe, baryogenesis, and the thermal history of dark sectors. Solving the Boltzmann equation for arbitrary $n_{\text{in}} \to n_{\text{out}}$ scattering processes is a computationally formidable task due to the collision integral's high-dimensional nature and intricate phase-space constraints set by energy-momentum conservation. While a variety of approaches have been developed for $2 \to 2$ processes, including analytic reductions and moment-based simplifications, these do not generalize to higher-multiplicity processes required by phenomena such as cannibal interactions, semi-annihilations, or freeze-in through multi-body decays.

The paper introduces Best (Boltzmann Equation Solver for Thermalization), a Python-based framework designed for direct, momentum-resolved solution of the Boltzmann equation for general $n_{\text{in}} \to n_{\text{out}}$ scatterings. The framework emphasizes rigorous enforcement of conservation laws, generic process specification, and supports quantum statistics, time-dependent masses, multiple species, and cosmological expansion, leveraging parallelization and adaptive Monte Carlo integration for scalability.

Collision Integral Decomposition and Treatment of Identical Particles

A critical technical advance in the framework is the formal decomposition of the collision integral for processes with unequal initial and final multiplicities involving identical particles. For reactions such as $\phi\phi \leftrightarrow \phi\phi\phi$, the full collision rate at a fixed momentum must sum over all ways the observed particle can be assigned to positions on either side of the vertex, each weighted by appropriate multiplicities.

Neglecting this decomposition introduces systematic energy nonconservation, as shown for the $2 \leftrightarrow 3$ case where omitting terms yields a $\sim 40\%$ error in total energy—far exceeding typical Monte Carlo statistical uncertainties. The framework implements the prescription

$$C_a(\mathbf{p}) = n_\alpha^a\, C_{n_\alpha}(\mathbf{p}) + n_\beta^a\, C_{n_\beta}(\mathbf{p}),$$

where $n_\alpha^a$ ($n_\beta^a$) is the count of species $a$ on each side. This decomposition is nontrivial for $2 \to 2$0, yet essential for correct thermalization and conservation law respect.

Figure 1: Top: particle labeling for $2 \to 2$1. Middle and bottom: decomposition of the collision integral; $2 \to 2$2 with the observed momentum on the two-particle side, $2 \to 2$3 on the three-particle side.

High-dimensional Monte Carlo Integration

Processes with $2 \to 2$4 particles entail a $2 \to 2$5-dimensional integration for each momentum grid point and time step. The code enforces exact momentum conservation by eliminating one momentum variable analytically, and approximates energy conservation via a narrow Gaussian delta function. Adaptive Monte Carlo (Vegas) integration, with separate learning curves for gain and loss terms, efficiently samples the high-variance regions of phase space, and integrator reuse between grid points and time steps enhances efficiency.

The method is benchmarked by comparing collision rates computed via Monte Carlo against a semi-analytical reduction for $2 \to 2$6 scatterings, showing agreement at the few percent level throughout the spectrum, with increased variance only near phase-space regions where the rate crosses zero.

Figure 2: Comparison of the collision rate $2 \to 2$7 from Vegas integration (points) versus semi-analytical calculation (solid line) for non-thermal initial distributions, in both massless (left) and massive (right) cases. Ratios indicate agreement within a few percent at high rates and increased scatter near the zero crossings due to Monte Carlo noise.

Thermalization Dynamics and Conservation Laws

Applying the framework to non-thermal initial conditions, the system exhibits robust evolution toward equilibrium Bose-Einstein distributions. In $2 \to 2$8 elastic scattering, both particle number and energy are conserved at the percent level, and the chemical potential and temperature of the final distribution are dictated by initial conserved quantities.

Figure 3: Thermalization via $2 \to 2$9 elastic scattering; left: distribution $n_{\text{in}} \to n_{\text{out}}$0 relaxes toward equilibrium. Right: energy and particle number conservation throughout dynamical evolution.

For $n_{\text{in}} \to n_{\text{out}}$1 number-changing processes, the numerical experiments illustrate that including only a subset of collision terms leads to spurious energy loss and incorrect equilibrium states. Incorporating the correct identical-particle decomposition exactly restores energy conservation and ensures that chemical potentials relax appropriately (to zero for $n_{\text{in}} \to n_{\text{out}}$2).

Figure 4: $n_{\text{in}} \to n_{\text{out}}$3 number-changing process; top: using only $n_{\text{in}} \to n_{\text{out}}$4 yields $n_{\text{in}} \to n_{\text{out}}$5 energy loss and incorrect equilibrium. Bottom: full decomposition recovers correct energy conservation and thermalization to Bose-Einstein equilibrium with $n_{\text{in}} \to n_{\text{out}}$6.

Parallelization and Computational Performance

The framework achieves near-linear scaling up to one processor per momentum grid point using MPI. Vegas's adaptive mapping and importance sampling are reused across grid points, minimizing redundant adaptation and integration cost. Benchmark wall times are $n_{\text{in}} \to n_{\text{out}}$7(minutes) per time step for realistic grid sizes and integration parameters, enabling practical computations on high-performance clusters for high-dimensional processes.

Implications and Prospects for Future Research

Practical Implications:

• The ability to directly evolve full, non-equilibrium momentum spectra for arbitrary processes enables precision studies of scenarios where kinetic and chemical equilibrium are simultaneously lost, such as cannibal and strongly interacting massive particle (SIMP) dark matter.
• The implementation is agnostic to model assumptions about cross section structure, quantum statistics, or cosmological expansion, suggesting broad utility across particle and cosmological applications involving number-changing and freeze-in/out processes.

Theoretical Implications:

• The correct treatment of identical particles in phase-space resolved collision integrals fills a notable gap in the literature. Previous solvers handled $n_{\text{in}} \to n_{\text{out}}$8 and $n_{\text{in}} \to n_{\text{out}}$9 processes where such subtleties are masked by symmetry, but higher-multiplicity interactions cannot ignore these terms without encountering unphysical results.
• The demonstrated technical feasibility of Monte Carlo evaluation in high dimensions under rigorous conservation constraints opens the door to exploring a much wider range of thermalization phenomena and initial conditions that were previously computationally inaccessible.

Future Directions:

• Extension to momentum-dependent, non-contact matrix elements, as arise in inelastic or multiparticle dark sector models.
• Self-consistent treatment of time-dependent masses (e.g., across cosmological phase transitions), which impact both kinematics and occupation numbers.
• Studies of freeze-in production, cascade decay chains, and scenarios where multi-particle interactions set the relic abundance.
• Rigorous uncertainty quantification in the tails of the distribution where indirect detection or structure formation may be sensitive.

Conclusion

The Best framework constitutes a comprehensive, extensible platform for simulating the non-equilibrium thermal dynamics of particle species in the early universe, particularly when kinetic and chemical equilibrium may break down and multi-particle interactions are dominant. Its core advances—treatment of the collision integral for arbitrary process multiplicities, careful attention to identical-particle contributions, and scalable parallel computation—provide a foundation for robust and theoretically consistent modeling of out-of-equilibrium phenomena in cosmology and high-energy physics.