Rapid and accurate numerical evolution of linear cosmological perturbations with noncold relics

We describe the implementation of a new approach to the numerical evaluation of the effects of noncold relics on the evolution of cosmological perturbations. The Boltzmann hierarchies used to compute the contributions of these relics to the stress-energy tensor are replaced with a set of integral equations. These integral equations take the form of convolutions and are solved iteratively with the rest of the system. We develop efficient algorithms for evaluating these convolutions using nonuniform fast Fourier transforms (NUFFTs). This approach enables efficient and accurate evaluation of the cosmic microwave background anisotropies and matter power spectra, all the way through the history of the Universe, without relying on semianalytic approximations at late times. We implement this method in the Boltzmann solver class, resulting in a new code called classier (for class Integral Equation Revision), and apply it to massive-neutrino perturbations as a demonstration. The implementation is optimized to accurately capture the distinct behaviors of perturbations in both super-/near-horizon and subhorizon regimes. Our results match the accuracy of a fully converged Boltzmann hierarchy solution while avoiding numerical artifacts from truncation of the Boltzmann hierarchy at finite multipole and offering substantial speedups depending on the required precision and the range of scales of interest. This new framework provides a practical and robust alternative for the truncated Boltzmann hierarchy approach, especially for studying beyond $\mathrm{\ensuremath{\Lambda}}\mathrm{CDM}$ noncold relics with signatures on small scales. classier is publicly available.