BSFfast: Rapid computation of bound-state effects on annihilation in the early Universe

Bound-state formation (BSF) can have a large impact on annihilation of new physics particles with long-range interactions in the early Universe. In particular, the inclusion of excited bound states has been found to strongly reduce the dark matter abundance and qualitatively modify the associated freeze-out dynamics. While these effects can be captured by an effective annihilation cross section, its explicit computation is numerically expensive and therefore impractical for repeated use in Boltzmann solvers or parameter scans. In this work we present BSFfast, a lightweight numerical tool that provides precomputed, tabulated effective BSF cross sections for a wide class of phenomenologically relevant models, including highly excited bound states and, where applicable, the full network of radiative bound-to-bound transitions. We exploit rescaling relations of the cross section to efficiently cover models with additional free parameters and provide fast interpolation routines in Mathematica, python and C for use in Boltzmann solvers. As an illustration, we apply BSFfast to a superWIMP scenario with a colored mediator, demonstrating that the tool enables phenomenological studies that would otherwise be computationally prohibitive. The code is publicly available on GitHub.