KineticXGPU: A Tensorized Collision Operator for Dark-Sector Self-Scattering

Technical Summary: KineticXGPU: A Tensorized Collision Operator for Dark-Sector Self-Scattering

Problem Statement
The microscopic nature of dark matter (DM) remains unknown, with alternative production histories like freeze-in gaining traction due to the lack of direct detection signals. In scenarios where DM is produced via multiple channels at different times (e.g., early thermal decays and late condensate decays), the resulting phase-space distribution can be non-thermal and bimodal. While self-interactions can relax these distributions toward thermal equilibrium (Local Thermal Equilibrium, LTE), the efficiency of this process depends on the elastic $2 \to 2$ scattering rate.

Solving the full Boltzmann equation (fBE) to track the momentum-dependent evolution of such distributions is numerically demanding. The elastic collision operator couples different regions of phase space, requiring a summation over the momentum grid for every bin. Existing public tools for dark matter evolution (e.g., DRAKE, BEST) often rely on CPU-based implementations or fluid approximations (coupled Boltzmann equations, cBE) that assume a thermal shape. Fluid approximations fail to capture non-thermal features, while CPU-based momentum-resolved solvers struggle with the computational cost of the high-dimensional integrals required for the collision term.

Methodology
The authors present KineticXGPU, a PyTorch-based implementation of the $2 \to 2$ elastic self-collision operator designed for GPU acceleration.

1. Discretization and Tensorization:

   • The Boltzmann equation is solved for the phase-space distribution $f(t, p)$ in an expanding universe. The authors switch to comoving momentum $q = a(t)p$ to remove the Hubble redshifting term from the left-hand side.
   • The momentum space is discretized on a logarithmic grid. The collision operator, which involves a four-dimensional integral over momenta constrained by energy-momentum conservation, is reduced to a sum over grid points.
   • Crucially, the discretized operator is reformulated as a series of tensor contractions. The gain-minus-loss structure ($f_n f_m - f_i f_{\tilde{j}}$) is evaluated by summing over internal momentum indices ($n, m$) for each external bin ($i$). This structure allows the operations to be highly parallelized, mapping naturally onto GPU architectures.

2. Numerical Implementation:

   • Integration: The time evolution uses the scale factor $a$ (or $u = \log a$) as the independent variable. The integration strategy is adaptive:
     • At early times where self-scattering is negligible ($\Gamma_{2\to2} \ll H$), a fixed-step fourth-order Runge-Kutta method is used with the collision term set to zero.
     • Once the scattering rate becomes significant ($\Gamma_{2\to2}/H > 10^{-2}$), the solver switches to an adaptive Heun method to handle the nonlinearity of the collision term.
   • Conservation: To ensure strict conservation of particle number and energy (zeroth and first moments), the authors implement a projection technique. The elastic contribution is projected onto the subspace with vanishing moments to correct for numerical errors arising from finite grid ranges, interpolation, and quadrature.
   • Validation: The code is validated against the semi-analytical evaluator in the BEST framework. Comparisons show agreement in the collision operator values, with residuals primarily occurring near non-thermal features where interpolation sensitivity is highest.

3. Physical Scenario:
   The code is applied to a "two-source freeze-in" scenario involving a scalar DM particle $\chi$.

   • Source 1: Higgs portal decays ($h \to \chi\chi$) occurring early in the thermal bath.
   • Source 2: Decay of a cold scalar condensate ($\Phi \to \chi\chi$) occurring later, injecting a monoenergetic population with a distinct momentum.
   • Interaction: A quartic contact interaction ($\lambda \chi^4$) drives elastic self-scattering.

Key Results

• Performance: The tensorized approach demonstrates a significant computational advantage. Runtime benchmarks show that the GPU implementation (using PyTorch) outperforms the CPU implementation and the semi-analytical BEST evaluator, particularly as the number of momentum bins ($N_{grid}$) increases. The GPU runtime scales efficiently due to the parallel nature of the tensor contractions.
• Evolution of Bimodal Distributions:
  • In the absence of self-interactions, the final DM distribution retains a clear bimodal shape, reflecting the two distinct production epochs.
  • As the self-coupling strength ($\lambda$) increases, elastic scattering progressively erases the separation between the two momentum components.
  • For sufficiently strong interactions, the distribution relaxes toward a Maxwell–Boltzmann shape, consistent with LTE.
• Comparison with Fluid Approximations: The authors compare the full phase-space solution (fBE) with a coupled fluid treatment (cBE) that assumes a thermal shape with time-dependent temperature and chemical potential.
  • In the weakly interacting regime, the fBE solution remains visibly non-thermal, while the cBE assumes a thermal shape.
  • The fluid approximation fails to capture the persistence of the non-thermal features and the specific relaxation dynamics observed in the full kinetic treatment.

Significance and Claims
The paper claims that KineticXGPU provides a robust, publicly available tool for studying dark sector evolution where momentum-dependent effects are critical. Its primary significance lies in:

1. Enabling Kinetic Studies: It allows for the study of relaxation processes in scenarios where the hierarchy of rates does not guarantee LTE, specifically when multiple production channels create non-thermal distributions.
2. Computational Efficiency: By leveraging GPU acceleration through tensor contractions, it makes solving the full Boltzmann equation with elastic self-scattering computationally feasible, overcoming the bottlenecks of previous CPU-based approaches.
3. Physical Insight: The application to the two-source scenario demonstrates that self-interactions can significantly alter the small-scale matter power spectrum by smoothing or erasing non-thermal features, effects that would be missed by fluid approximations.

The authors conclude that while the code focuses on the unperturbed distribution and the self-collision operator, the computational framework is general and can be extended to other dark sector models or perturbation calculations. The code is made available on GitHub to facilitate reproducibility and further development in the field.