🔭 astrophysics

Chebyshev interpolation in Einstein-Boltzmann codes

This paper demonstrates that replacing traditional cubic splines with Chebyshev polynomial interpolation in Einstein-Boltzmann codes significantly reduces interpolation error and accelerates computations for cosmological spectra by achieving higher precision with fewer evaluation points, while also presenting a practical method to handle non-integer multipole nodes.

Original authors: Herman Sletmoen

Published 2026-08-26

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Herman Sletmoen

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

To understand the universe, cosmologists build elaborate mathematical models that describe how matter and light have behaved since the Big Bang. These models rely on a set of complex equations that track the evolution of the cosmos over time, across different scales of distance, and in every direction of the sky. To turn these equations into predictions that can be compared with real telescope data, scientists use powerful computer programs known as Einstein–Boltzmann codes. These programs act as virtual laboratories, simulating the birth of the cosmic microwave background—the faint afterglow of the Big Bang—and the distribution of galaxies. However, running these simulations is incredibly difficult because the equations are stiff and change rapidly at certain moments, such as when the universe cooled enough for atoms to form. To make the calculations manageable, the programs traditionally solve the equations only at a few specific points and then fill in the gaps between them using a method called interpolation, which is essentially drawing a smooth curve through known data points to guess the values in between. The accuracy of the entire simulation depends heavily on how well this guessing game works.

A researcher at the University of Oslo, Herman Sletmoen, has found a way to make this guessing game significantly more accurate and faster. In a new study, Sletmoen demonstrates that by changing the mathematical technique used to fill in those gaps, the computer can achieve much higher precision with far fewer calculations. The traditional method, which uses cubic splines—a technique that connects points with small, smooth curves—works well but requires a large number of data points to reach high levels of accuracy. Sletmoen replaced this with a different approach based on Chebyshev polynomials, a mathematical tool that is exceptionally good at fitting smooth curves when the data points are placed at specific, non-uniform intervals. By solving the difficult physics equations only at these special points and then using the new method to fill in the rest, the computer can produce results with interpolation errors up to ten thousand times lower than the old method for the same amount of work.

The core of the problem lies in the nature of the data these codes generate. The universe evolves smoothly over time and space, except for a few brief moments of rapid change. When a computer tries to simulate this, it must calculate the behavior of light and matter for billions of different distances and angles. Doing this for every single possibility is impossible, so the code calculates a few thousand examples and then interpolates the rest. For decades, the standard tool for this job has been cubic splines. While reliable, this method is somewhat blunt; to get a very precise answer, it needs to be fed a dense grid of data points, which slows down the simulation. Sletmoen's work shows that because the underlying physics is smooth, a more sophisticated mathematical approach can do the job with a fraction of the data. By placing the calculation points where the new method needs them most, rather than spreading them out evenly, the computer can capture the fine details of the universe's evolution much more efficiently.

In the study, Sletmoen tested this new approach on the two most critical variables in the simulation: the scale of the waves in the early universe and the angle at which we observe the sky. When the code used the new method to fill in the gaps for the wave scales, the error in the final prediction dropped dramatically. While the old method needed around two hundred data points to reach a certain level of precision, the new method reached the same or better precision with only fifty to eighty points. This reduction in required calculations translates directly to speed. For a typical simulation of the cosmic microwave background temperature, the new method made the entire process two and a half to four times faster while maintaining the same high standard of accuracy. This is a significant gain, as these simulations are often run thousands of times to test different theories about the universe.

The challenge was not just in the math, but in how the computer handles the data. The new method works best when the points are placed at specific, uneven intervals, but the physical quantities being measured, such as the angle on the sky, only exist as whole numbers. You cannot have a measurement at half a degree in the same way you can have a measurement at half a meter. To solve this, Sletmoen developed a clever workaround. Instead of forcing the computer to calculate values at impossible, non-integer angles, the method rounds the special calculation points to the nearest whole number. This small adjustment allows the powerful new math to be used without breaking the rules of the physical simulation. The study found that this "rounded" approach performed just as well as the ideal mathematical version, proving that the technique is robust enough for real-world cosmology.

The results suggest a shift in how these complex simulations are built. Older computer codes often used shortcuts and approximations to speed things up, switching between different equations depending on the situation. While fast, these shortcuts can introduce small bumps or irregularities in the data that make the traditional interpolation methods struggle. The new approach shines brightest when used with a newer generation of codes that avoid these shortcuts and solve the equations continuously, keeping the data perfectly smooth. By pairing these smooth, approximation-free codes with the new interpolation method, scientists can get the best of both worlds: the speed of fewer calculations and the precision of a highly accurate model. This combination allows researchers to explore the universe with a clarity that was previously too expensive to compute, opening the door to more detailed tests of cosmological theories without needing to wait for faster computers.

Technical Summary: Chebyshev Interpolation in Einstein–Boltzmann Codes

Problem Statement
Einstein–Boltzmann codes (e.g., CLASS, CAMB, SymBoltz) compute theoretical cosmological predictions by solving stiff ordinary differential equations (ODEs) and line-of-sight integrals. To manage computational cost, these codes evolve perturbations on coarse grids of independent variables—conformal time ( $\tau$ ), wavenumber ( $k$ ), and angular multipole ( $\ell$ )—and interpolate to finer grids. While interpolation is a standard optimization, traditional methods like cubic splines introduce errors that scale algebraically ( $n^{-4}$ ). Furthermore, "approximation-based" codes often switch equations at specific thresholds to speed up calculations, introducing non-smooth features ("kinks") that degrade the performance of global interpolation methods. Conversely, "approximation-free" codes (e.g., SymBoltz, DISCO-EB) maintain fully smooth solutions but require efficient interpolation to offset their higher per-point computational cost.

Methodology
The paper proposes replacing traditional cubic splines with Chebyshev polynomial interpolation for the independent variables $k$ and $\ell$ . The methodology relies on the following principles:

Node Selection: Instead of equispaced points, which suffer from Runge's phenomenon, the method utilizes Chebyshev nodes (specifically the second kind, including endpoints) where the density of points scales as $(1-x^2)^{-1/2}$ .
Evaluation: Interpolation is performed using the barycentric formula, which is numerically stable and requires only $O(n)$ operations per evaluation point.
Variable-Specific Implementation:
- Time ( $\tau$ ): The author explicitly does not apply Chebyshev interpolation here. The adaptive nature of ODE solvers and the presence of local physical features (e.g., recombination) make rigid node placement inefficient. The paper notes that specialized dense output methods within ODE solvers are superior for time interpolation.
- Wavenumber ( $k$ ): Perturbation ODEs are solved explicitly only at Chebyshev nodes within the interval $[k_{min}, k_{max}]$ . Source functions $S(\tau, k)$ are then interpolated to arbitrary $k$ values using the barycentric formula.
- Multipole ( $\ell$ ): A challenge arises because $\ell$ $ℓ$ is physically defined only as an integer, while Chebyshev nodes are non-integer. The paper presents two solutions:
  - Generalized: Extend spherical Bessel functions $j_\ell(x)$ to non-integer $\ell$ (solving the ODE continuously) to allow standard Chebyshev interpolation.
  - Integer-Rounded: Compute $j_\ell(x)$ only at integer $\ell$ but round the Chebyshev nodes to the nearest integers. This "near-Chebyshev" approach maintains stability and avoids the need for non-integer Bessel function implementations.
- Piecewise Handling: For the lensing spectrum ( $C^{LL}_\ell$ ), where the Limber approximation introduces a non-smooth transition, the author employs a piecewise Chebyshev interpolant (separate polynomials for $\ell < 10$ and $\ell \ge 10$ ) to restore rapid convergence.

Key Results
The author implemented these methods in the approximation-free code SymBoltz and compared performance against cubic splines on matter and CMB spectra:

Convergence Rate: Chebyshev interpolation exhibits geometric (exponential) convergence for smooth functions, whereas cubic splines converge algebraically.
Error Reduction: For source functions in $k$ , Chebyshev interpolation achieves errors four orders of magnitude lower than cubic splines using the same number of points (e.g., at $N=80$ ).
Efficiency: To reach a relative error of $10^{-4}$ to $10^{-5}$ in the CMB temperature spectrum ( $C^{TT}_\ell$ ), Chebyshev interpolation requires only 50–80 points per variable ( $k$ and $\ell$ ). In contrast, cubic splines require approximately 200 points to reach similar precision.
Speedup: This reduction in required sampling points translates to a 2.5× to 4× speedup in the full computation for a fixed high-precision target.
Robustness: The integer-rounded Chebyshev method for $\ell$ performs comparably to the generalized non-integer method, suggesting it is a viable strategy for curved universes where hyperspherical Bessel functions are typically computed via integer-based recurrences.

Significance and Claims
The paper argues that Chebyshev interpolation creates a "synergy" with approximation-free Einstein–Boltzmann codes. While approximation-based codes rely on equation switching to gain speed (at the cost of smoothness), approximation-free codes generate smooth solutions that are ideally suited for the rapid convergence of Chebyshev polynomials. By leveraging this smoothness, approximation-free codes can recover performance advantages without sacrificing precision.

The author positions this work as a step toward simpler numerical algorithms that exploit the physical properties of the equations (smoothness) rather than relying on complex approximation schemes. They acknowledge that the improvement is most dramatic at high precision levels and that the global nature of Chebyshev interpolation requires careful handling of non-smooth features (like the Limber switch), which can be mitigated via piecewise interpolation. The implementation is available in the SymBoltz package.

More like this