Boltzmann Equation Solver for Thermalization

TL;DR

BEST is a Python framework that accurately solves the momentum-resolved Boltzmann equation for complex scattering processes, incorporating energy conservation, quantum statistics, and cosmological effects, validated against semi-analytical results.

Contribution

The paper introduces a novel, flexible Python code for solving the Boltzmann equation with exact conservation laws and efficient parallelization, supporting complex particle interactions and cosmological scenarios.

Findings

Validated Monte Carlo results against semi-analytical collision integrals.

Demonstrated correct energy conservation in 2 to 3 particle processes.

Showcased thermalization of a massive scalar field with proper conservation laws.

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.