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

SymBoltz is a new Julia package for solving the linear Einstein-Boltzmann equations in cosmology. It features a symbolic-numeric interface for specifying equations, is free of approximation switching schemes, and is compatible with automatic differentiation. Cosmological models are built from replaceable physical components in a way that scales well to extended models, or alternatively written as one compact system of equations. The modeler provides their equations, and SymBoltz solves them while reducing friction in the modeling process. Symbolic knowledge enables powerful automation of tasks, such as separating computational stages (e.g., background and perturbations), generating analytical and sparse Jacobian matrices, and interpolating arbitrary variables from the solution. Implicit solvers integrate the full stiff equations at all times without approximations, which greatly simplifies the code. Performance remains comparable to existing approximation-based codes due to efficient high-order implicit methods, fast generated code, optimal handling of the Jacobian, and sparse matrix methods. Automatic differentiation gives exact derivatives of any output with respect to any input, which is important for gradient-based Markov chain Monte Carlo (MCMC) sampling in large parameter spaces, Fisher forecasting, emulator training, and sensitivity analyses. The main features form a synergy that reinforces the design of the code. Output spectra agree with established codes up to 0.1% with standard precision levels. More work is needed to implement additional features and for fast reverse-mode automatic differentiation of scalar loss functions. SymBoltz is publicly available with single-command installation and extensive documentation. We welcome all contributions to the code from the community.