Chebyshev interpolation in Einstein-Boltzmann codes

Einstein-Boltzmann codes compute theoretical predictions of cosmological models and rely heavily on interpolation in their independent variables: time τ, wavenumber k and multipole. We give a practical summary of interpolation with Chebyshev polynomials, which converges rapidly for smooth functions and thus pairs naturally with approximation-free Einstein-Boltzmann codes. By solving the perturbations and line-of-sight integrals at Chebyshev nodes in k and, we show that Chebyshev polynomials interpolate to higher precision than traditional cubic splines from fewer explicit solutions. On a set of example spectra for matter and the cosmic microwave background (CMB), we find up to four orders of magnitude lower interpolation error using the same number of points. For a typical CMB temperature spectrum computed with interpolation in both k and, Chebyshev polynomials converge to 10-4 - 10-5 relative error with only 50-80 points per variable, while cubic splines approach 10-4 error with 200 points, translating to a 2.5 times - 4 times speedup. The exact improvement depends on the target function and is generally more dramatic at high precision levels. Standard Chebyshev-interpolation needs line-of-sight integrals generalized to non-integer, but we show a way to avoid this by rounding the nodes to integers. Chebyshev interpolation is implemented in SymBoltz, which is available at https://github.com/hersle/SymBoltz.jl.