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

TL;DR

BSFfast tackles the computational bottleneck of including bound-state formation in early-Universe annihilation by providing precomputed, tabulated $\langle \sigma v \rangle_{\rm eff}$ for a wide class of long-range interacting particles, including highly excited states up to $n=100$. It leverages exact closed-form rates for BSF, decays, and bound-to-bound transitions, and exploits mass and coupling rescaling in the frozen-coupling limit plus an approximate running-coupling scheme to cover broad parameter spaces with minimal tabulation. The authors detail model coverage, numerical implementation (up to 5050 bound-state channels, RunDec for SM running), and unitarity considerations, and demonstrate the tool with a phenomenological superWIMP scenario showing significant late-time effects on relic densities and structure formation. The resulting fast interpolation framework enables efficient integration into Boltzmann solvers and parameter scans, with public code available for the community to extend to additional gauge groups and interaction types.

Abstract

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.

Figures (6)

• Figure 1: Bound state contribution $\langle\sigma v\rangle_\text{eff,BSF}$ to the thermally averaged effective cross section for a scalar $X$ with quantum numbers identical to those of the right-handed up-type quarks ('stop'-like, see Tab. \ref{['tab:model-coverage']}) provided by BSFfast. The blue line shows the full result when taking bound state formation/ionization, bound state decays and transitions among bound states into account. The black line corresponds to the case when neglecting transitions (technically applicable to a colour-charged but electrically neutral scalar particle). The dashed line shows the contribution from the ground state only ($n=1$) for illustration. All lines correspond to $m =10^6$ GeV and use the running QCD as well as QED couplings.
• Figure 2: Comparison of QCD-S models and approximately rescaled dQCD-S models (darker lines) where the SU(3) gauge coupling at every $x$ is chosen at the potential scale $\alpha_\text{s}=\alpha_\text{s}(\sqrt{2m T})$. Upper panel: $m^2 \left< \sigma_{} v \right>_\text{eff,BSF}$ plotted over temperature for different masses (different colours). The black dashed line shows the result for setting $\alpha=0.05$ at all $x$ which approximates $\alpha_\text{QCD}(10^7\,\text{GeV})=0.0460$. Note that choosing a different $\alpha$ here merely shifts the drawn curve around but affects neither its shape nor slope. Lower panel: Relative deviation of the approximately rescaled dQCD-S to QCD-S in percent.
• Figure 3: Ratio of the summed $s$-wave BSF cross section to the corresponding unitarity bound, $(\sigma v)_\text{BSF}^{\ell=0}/(\sigma v)_\text{uni}^{\ell=0}$, for SM QCD, plotted as a function of inverse velocity. The vertical gray line indicates the virial velocity corresponding to the highest $x$ included in our tables, $v=\sqrt{6/x_\text{max}}$. The coloured region, where the curves begin to fan up at low velocities, indicates the onset of deviations between two prescriptions for the running coupling at energy scales below 1 GeV (plateau versus cutoff).
• Figure 4: Comparison of the effective thermally averaged cross section $\langle \sigma v \rangle_\text{eff,BSF}$ to estimates of the corresponding unitarity bounds. The blue solid and blue dashed curve correspond to the models with (QCD-SU) and without (QCD-S) bound-to-bound transitions, respectively, choosing $m=10^6\,$GeV. The black solid line corresponds to the most robust unitarity bound summing over all potentially relevant partial wave contributions at a given $x$ (see text for details) while the black dashed line displays the thermally averaged $s$-wave contribution.
• Figure 5: Unitarity violation in dark QCD) with constant coupling. Problematic and unproblematic regions are shown in yellow and blue, respectively, including bound states up to $n=4000$. The red and orange contours indicate where $100\%$ and $10\%$ of the $s$-wave unitarity bound $(\sigma v)^{\ell=0}_\text{uni}$ are reached. Dashed lines show an extrapolation of the low-$v$ boundary of the unitarity-violating region to $n_\text{max}\to\infty$ based on the displayed data, see text for details.