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.

Figures (4)

• Figure 1: Top: particle labeling for $\phi\phi \leftrightarrow \phi\phi\phi$. Middle and bottom: decomposition of the collision integral. The red leg denotes the observed momentum $\mathbf{p}$; black legs are integrated over. Each term is the net rate (gain minus loss) at the observed momentum. Middle: $C_2(\mathbf{p})$, with $\mathbf{p}$ on the 2-particle side. Bottom: $C_3(\mathbf{p})$, with $\mathbf{p}$ on the 3-particle side. The full collision integral is $C_\phi(\mathbf{p}) = 2\,C_2 + 3\,C_3$.
• Figure 2: Comparison of the collision rate $|C[f](p)|$ computed by Vegas Monte Carlo (points) and the semi-analytical method (solid line), for a non-thermal initial distribution with $\lambda = 1$ and $H = 0$. Left panels: massless ($m_\phi = 0$). Right panels: massive ($m_\phi = 1$). Upper panels show the absolute rate on a log-log scale; lower panels show the ratio. Both methods include the identical-particle multiplicity factor of 2. For both cases, agreement is within a few percent where the collision rate is large; increased scatter appears near the zero crossings of $C[f]$ where the Monte Carlo signal-to-noise ratio is low.
• Figure 3: Thermalization via $2 \to 2$ elastic scattering with massive $\phi$ ($m_\phi = 1$, $\lambda = 1$, $H = 0$). Left panel: evolution of $f(p)$ from a non-thermal initial condition (red) toward the Bose-Einstein distribution (black dashed), with $T_{\rm eq}$ and $\mu$ determined by the conserved energy and particle number. Right panel: conservation of particle number $N/N_0$ (solid) and energy $E/E_0$ (dashed).
• Figure 4: $2 \leftrightarrow 3$ number-changing process with massive $\phi$ ($m_\phi = 1$, $\lambda_5 = 1$, $H = 0$). Left panels: evolution of $p^2 f(p)$ from a non-thermal initial condition (red) toward the equilibrium Bose-Einstein distribution with $\mu = 0$ (black dashed). Right panels: conservation of particle number $N/N_0$ (solid) and energy $E/E_0$ (dashed). Top: incorrect calculation using only $C_3$ (shown for illustration), leading to $\sim 40\%$ energy loss. Bottom: correct decomposition $2\,C_2 + 3\,C_3$, which restores energy conservation.