SymBoltz.jl: A symbolic-numeric, approximation-free, and differentiable linear Einstein-Boltzmann solver

SymBoltz.jl is a new Julia package that provides a symbolic-numeric, approximation-free, and fully differentiable solver for linear Einstein-Boltzmann equations, offering exact derivatives and flexible model construction while maintaining performance comparable to existing approximation-based codes.

The Gist

Imagine you are trying to predict the future of the entire universe. You want to know how galaxies form, how light from the Big Bang travels to us, and how dark matter pulls everything together. To do this, scientists use complex mathematical recipes called Einstein–Boltzmann equations.

For decades, the standard tools for solving these recipes (like CAMB and CLASS) have been like old, powerful, but rigid industrial machines. They are incredibly fast and reliable, but they are built with a lot of "shortcuts" (approximations) to keep the gears turning. If you want to change the recipe slightly—say, to test a new theory of gravity—you often have to take the machine apart, weld new parts in, and hope you didn't break the engine. It's like trying to modify a car engine while it's still running at 100 mph.

Enter SymBoltz, a new tool written in the Julia programming language. Think of SymBoltz not as a rigid machine, but as a smart, magical kitchen.

Here is how it works, broken down into three simple superpowers:

1. The "Lego" Kitchen (Symbolic-Numeric Interface)

In the old way, if you wanted to swap out the "Dark Energy" ingredient for a new, weird flavor, you had to rewrite the instructions for the whole kitchen. In SymBoltz, you just write down your new recipe in plain math (like writing on a whiteboard).

• The Analogy: Imagine you are building a model city. In the old way, you had to glue every brick together permanently. If you wanted to change a building, you had to smash the whole block. In SymBoltz, you use Lego bricks. You can snap a "Dark Energy" block out and snap a "Modified Gravity" block in. The system automatically figures out how to connect them, calculates the physics, and builds the city for you. You don't need to be a master engineer; you just need to know what you want to build.

2. The "No-Shortcuts" Chef (Approximation-Free)

The old machines use shortcuts because the math is so stiff and difficult that solving it exactly would take forever. They say, "Oh, the particles are moving too fast to calculate exactly, let's just guess." This is like a chef saying, "I'll just guess the amount of salt instead of measuring it."

SymBoltz refuses to guess. It uses a special type of math solver (an implicit solver) that is like a super-precise, slow-cooking pressure cooker.

• The Analogy: Instead of rushing and guessing, this chef takes the time to solve the exact equation at every single moment. Because the software is so smart, it doesn't need to rush. It handles the "stiff" parts of the math (where things change wildly) without needing to cut corners. This means the results are cleaner and more honest, and you don't have to worry about whether your shortcut was valid for this specific universe.

3. The "X-Ray Vision" (Differentiability)

This is perhaps the coolest feature. In the old days, if you wanted to know how changing one tiny ingredient (like the amount of dark matter) would change the final taste of the universe, you had to bake the cake 100 times, changing the ingredient slightly each time, and compare the results. It was slow and clumsy.

SymBoltz has X-ray vision. Because it knows the exact math behind every step, it can instantly tell you: "If you add a pinch more salt, the cake will be 0.05% sweeter."

• The Analogy: Imagine you are driving a car. The old way is to drive to a destination, turn around, drive back, change the route slightly, and drive again to see if it's faster. SymBoltz is like having a GPS that instantly calculates the exact effect of every possible turn on your arrival time. This allows scientists to use powerful AI-like methods (called MCMC) to find the "perfect" universe model much faster than ever before.

Why Does This Matter?

We are living in a time where our telescopes are getting incredibly sharp. We are finding cracks in our current understanding of the universe (tensions between what we see and what we think). To fix these cracks, we need to test wild, new ideas.

• Old Tools: Great for standard recipes, but hard to tweak.
• SymBoltz: Built for experimentation. It lets scientists write down their wildest theories in a few lines of code, run them without guessing, and instantly see how they change the universe.

The Bottom Line

SymBoltz is a modern, flexible, and precise simulator for the universe. It trades the "brute force" speed of old machines for elegance, accuracy, and ease of use. It's like upgrading from a hammer and chisel to a 3D printer: you can still build the same things, but now you can design and build anything you can imagine, and the machine handles the heavy lifting for you.

The paper shows that SymBoltz is already as accurate as the best tools we have, and it opens the door to a future where testing new theories of the cosmos is as easy as writing a few equations.

Technical Summary

1. Problem Statement

Cosmological modeling is facing a crisis of precision. Next-generation surveys (e.g., Euclid, Rubin Observatory, DESI) are revealing tensions with the standard $\Lambda$CDM model, necessitating the exploration of complex, extended cosmological models. However, existing tools for solving linear Einstein–Boltzmann equations (the "Boltzmann codes" like CAMB and CLASS) face significant limitations:

• Approximation Schemes: To handle the extreme stiffness of the equations (caused by widely varying timescales, e.g., tightly coupled baryon-photon fluids), traditional codes rely on approximation schemes (e.g., Tight-Coupling Approximation, Ultrarelativistic Fluid Approximation). These require switching between different sets of equations based on thresholds, complicating the code and introducing potential errors when modifying models.
• Monolithic Structure: Traditional codes are written in low-level languages (C/Fortran) with rigid, pipeline-like structures. Modifying them to include new physics (e.g., new particle species or modified gravity) requires scattering changes across multiple files (input, background, perturbations, output), making the codebase fragmented and difficult to maintain.
• Lack of Differentiability: Modern statistical inference (e.g., Hamiltonian Monte Carlo, Fisher forecasting) and machine learning (emulator training) require exact gradients of outputs with respect to inputs. Traditional codes lack native differentiability, forcing users to rely on slow, approximate finite differences or train separate neural network emulators.
• The "Two-Language" Problem: Data analysis often occurs in high-level languages (Python), while solvers are in low-level languages, creating friction in the workflow.

2. Methodology

The author introduces SymBoltz.jl, a new Julia package designed to address these issues through three core pillars:

A. Symbolic-Numeric Interface

• Component-Based Modeling: Instead of a rigid pipeline, SymBoltz structures models around physical components (metric, gravity, photons, baryons, dark matter, etc.). Users define equations symbolically using ModelingToolkit.jl and Symbolics.jl.
• Automatic Compilation: The symbolic equations are compiled into fast, low-level numerical code. The system automatically:
  • Separates equations into computational stages (background, perturbations, line-of-sight integration).
  • Distinguishes between "unknown" variables (integrated via ODEs) and "observed" variables (algebraically derived), allowing users to reference any variable without manual re-derivation.
  • Generates analytical Jacobian matrices and detects sparsity patterns automatically.
• Modularity: New physics can be added by defining a self-contained submodel (e.g., a specific dark energy equation of state) and plugging it into the main system, avoiding code duplication.

B. Approximation-Freeness via Implicit Solvers

• Implicit Integration: SymBoltz avoids all approximation schemes (TCA, UFA, etc.) by solving the full, stiff Einstein–Boltzmann equations using implicit ODE solvers (from OrdinaryDiffEq.jl).
• Optimization: To overcome the historical performance penalty of implicit solvers, SymBoltz utilizes:
  • High-order methods: Specifically Rosenbrock methods (e.g., Rodas5P) and ESDIRK methods.
  • Analytical Jacobians: Generated symbolically to ensure exactness and speed.
  • Sparse Linear Algebra: The Jacobian is stored in a sparse format with a precomputed sparsity pattern, allowing efficient LU factorization (using LinearSolve.jl).
• Stage Separation: The system automatically splines background variables to serve as inputs for perturbation stages, maintaining the efficiency of sequential solving without manual intervention.

C. Differentiability

• Automatic Differentiation (AD): SymBoltz is fully compatible with ForwardDiff.jl. It can compute exact derivatives of any output (e.g., power spectra $P(k)$, $C_\ell$) with respect to any input parameter.
• Current Limitation: Currently supports forward-mode AD (efficient for many outputs, few inputs). Reverse-mode AD (essential for MCMC with scalar loss functions) is a target for future work.

3. Key Contributions

1. First Approximation-Free Julia Solver: SymBoltz is the first solver to integrate the full stiff Einstein–Boltzmann equations without switching approximation schemes, relying instead on modern implicit solvers and symbolic code generation.
2. Symbolic-Abstraction Layer: It inverts the traditional Boltzmann code architecture, allowing modelers to write equations in a high-level, physical-component format that is automatically compiled into optimized numerical code.
3. Native Differentiability: It provides a differentiable pipeline for cosmological observables, enabling gradient-based inference and Fisher forecasting without external emulation or finite differences.
4. Automated Jacobian & Sparsity Handling: The system automatically generates analytical, sparse Jacobian matrices, a critical requirement for the performance of implicit solvers, removing the burden of manual derivation from the user.

4. Results

• Accuracy: SymBoltz outputs agree with CLASS (a gold-standard code) to within 0.1% for matter power spectra ($P(k)$) and CMB angular power spectra ($C_\ell$) using standard precision settings.
• Performance:
  • In approximation-free mode, SymBoltz is >10x faster than CLASS using its implicit solver (ndf15) at comparable precision.
  • When CLASS uses its full suite of approximations and explicit solvers, SymBoltz achieves comparable speed while remaining approximation-free.
  • The speed is achieved through high-order implicit methods, analytical Jacobians, and sparse matrix optimizations.
• Differentiability Validation:
  • Fisher Forecasting: SymBoltz successfully computed Fisher matrices for CMB constraints using automatic differentiation, matching results from finite-difference methods in CLASS but without the need for step-size tuning.
  • MCMC Sampling: A preliminary demonstration using the NUTS sampler (via Turing.jl) on Type Ia supernova data showed the pipeline works, though reverse-mode gradients are not yet optimized for full perturbation-based MCMCs.

5. Significance

SymBoltz represents a paradigm shift in cosmological modeling software:

• Accessibility: It lowers the barrier to entry for developing new cosmological models, allowing researchers to prototype complex theories (e.g., modified gravity, interacting dark sectors) in a few lines of high-level code.
• Robustness: By eliminating approximation switching, it removes a major source of systematic error and complexity, ensuring that the physics is solved as written.
• Future-Proofing: Its differentiability makes it a critical tool for the era of "precision cosmology," enabling efficient gradient-based sampling in high-dimensional parameter spaces and seamless integration with machine learning workflows.
• Ecosystem: As a Julia package, it aims to foster a modular ecosystem where Boltzmann solvers, N-body simulations, and inference tools can interoperate seamlessly, potentially solving the "two-language problem" in astrophysics.

In conclusion, SymBoltz.jl demonstrates that approximation-free, differentiable, and modular Einstein–Boltzmann solvers are not only feasible but can be highly performant, offering a modern alternative to the decades-old, approximation-heavy codes currently in use.