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

SymBoltz is a new Julia package that solves the linear Einstein-Boltzmann equations. 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 in model space, or alternatively written as one compact system of equations. The modeler should simply write down their equations, and SymBoltz solves them and eliminates friction in the modeling process. Symbolic knowledge enables powerful automation of tasks, such as separating computational stages like the background and perturbations, generating the analytical Jacobian matrix and its sparsity pattern, 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 as good as in existing approximation-based codes due to high-order implicit methods that take long time steps, fast generated code, optimal handling of the Jacobian and efficient 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 methods in large parameter spaces, training of emulators, Fisher forecasting and sensitivity analysis. The main features form a synergy that reinforces the design of the code. Results agree with established codes to 0.1% with standard precision. More work is needed to implement additional features and for fast reverse-mode automatic differentiation of scalar loss functions. SymBoltz is available at https://github.com/hersle/SymBoltz.jl with single-command installation and extensive documentation, and welcomes all contributions.