>
> This PR fixes https://github.com/SciML/SciMLSensitivity.jl/issues/1238
We need to do those two in order to make it generally stable. It’s near the top of my list. Then all MTK codes should be fine.
---
### Author:  [@marcobonici](https://discourse.julialang.org/u/marcobonici)
#### Post date: [October 4, 2025, 4:30am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/6 "2025-10-04T04:30:18Z")
@ChrisRackauckas , this is with Enzyme, right? Will you include support for other AD systems as well, like MoonCake?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 4, 2025, 4:34am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/7 "2025-10-04T04:34:07Z")
Mooncake is already supported. This is about making sure the rules, which apply to all AD front ends (Mooncake, Enzyme, ReverseDiff, Tracker, and Zygote) fully support the MTK features, which is currently missing a bit of SciMLStructures.jl support to be complete. We also have some early demonstrations of Enzyme.jl and Mooncake.jl direct adjoints on the solvers, now with integration tests that those do not regress (tests in the AD libraries that force them to keep compatibility 😅), but need to get this to support the implicit methods as well (currently only the explicit methods).
---
### Author:  [@hersle](https://discourse.julialang.org/u/hersle)
#### Post date: [October 4, 2025, 10:07am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/8 "2025-10-04T10:07:52Z")
> [@ChrisRackauckas](#):
>
> Really cool stuff, I had no idea what you were actually doing with MTK but now this might be the coolest application 😅
Thanks Chris! Yep, I remember when I first discovered it, I thought ModelingToolkit.jl + OrdinaryDiffEq.jl would be a very powerful combination for this project. I really appreciate your help and responsiveness with bug fixes and new features.
---
### Author:  [@hersle](https://discourse.julialang.org/u/hersle)
#### Post date: [October 4, 2025, 10:23am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/9 "2025-10-04T10:23:29Z")
> [@martinmestre](#):
>
> Do you think it could be applied with warm DM, like in CLASS?
Hi @martinmestre and thank you! Warm dark matter models are not implemented yet, but I think it would be very interesting to add them and to test the extensibility of the code. Do you have a particular model in mind or a reference paper with equations? I do not have much prior experience with these models, but I suppose they could be based off the implementation of massive neutrinos, which evolve momentum bins of their distribution function in time. Let me know and I will also help you!
---
### Author:  [@hersle](https://discourse.julialang.org/u/hersle)
#### Post date: [October 4, 2025, 6:30pm UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/10 "2025-10-04T18:30:01Z")
Hi and thanks, @marcobonici and @schiarenza!
> [@marcobonici](#):
>
> 1. **Backward AD**. You clearly state in the paper that this is not wokring at the moment. Do you have any idea on the time scale to have that working? with a student of mine, @schiarenza , we are working on the new version of her NonLimber code, [Blast.jl](https://github.com/sofiachiarenza/Blast.jl), which will soon include more contributions compared to the first release (e.g. RSD, Magnification bias, fNL and CMB lensing) and backward AD. It would be extremely cool if we could use your code as the cosmology engine beside one of my emulators.
> [@ChrisRackauckas](#):
>
> It’s probably more dependent on me? MTK integration with Gauss adjoint will just work out of the box, but
>
> …
>
> We need to do those two in order to make it generally stable. It’s near the top of my list. Then all MTK codes should be fine.
Like Chris says, this is mostly dependent on the underlying libraries, at least for what concerns the ODE solutions of the background and perturbations. I rely on standard ModelingToolkit.jl and OrdinaryDiffEq.jl functionality, so I hope reverse-mode support will eventually come “for free”.
I would love for you to train your emulator with SymBoltz! Right now it may not output all the quantities you want, though. Do you (plan to) use the direct output from CAMB/CLASS for RSD, fNL, etc.? Or do you compute it yourself from simpler quantities?
@ChrisRackauckas Just FYI, the ODE solving in SymBoltz relies on first solving one “background” ODE, followed by several independent “perturbation” ODEs where the “background” solution is supplied through interpolating DataInterpolations.jl splines (as time-dependent parameters). I don’t know if that makes reverse-mode more or less complicated. It would be fantastic if that works.
> [@marcobonici](#):
>
> 1. **Extensions**. Have you considered to include an extension like non clustering Dark Energy as in w0-wa? It shouldn’t be too difficult, as it is just a matter of including some contributions to the background (and would make a nice addition as this is the main extension to LCDM that is being analyzed).
Yes, w\_0w\_a dark energy [is included](https://hersle.github.io/SymBoltz.jl/stable/components/#SymBoltz.w0wa) (see also e.g. [here](https://hersle.github.io/SymBoltz.jl/stable/models/#ww%E2%82%90CDM-(CPL-parametrization)) or [here](https://hersle.github.io/SymBoltz.jl/stable/comparison/)). I am also interesting in implementing other models to properly test the extensibility of the code. Let me know if you have something in mind.
> [@marcobonici](#):
>
> 1. **CMB lensing**. Have you implemented the lensing of the CMB spectra or for the moment you just have the unlensed CMB spectra?
It is not implemented yet. For now I have focused on getting good agreement with CLASS for the background + perturbations + matter power spectrum + (unlensed) CMB power spectrum for ΛCDM, which is challenging on its own. I have not done the calculation before, but from what I can tell it looks like another line-of-sight integration with a different source function, which should be relatively easy to add to the current framework. I opened this issue, feel free to track it and drop any comments there:
>
>
> SymBoltz should be able to compute lensed CMB spectra with a few additions. See …e.g. https://arxiv.org/pdf/1312.2697#subsection.5.5. The calculation seems to be another \[line-of-sight integration\](https://hersle.github.io/SymBoltz.jl/stable/observables/#Line-of-sight-integration), only with a different source function:
>
> $$
> \\Delta\_l^\\phi(k,\\tau\_0) = \\int\_{\\tau\_i}^{\\tau\_0} d\\tau S^\\phi(k,τ) j\_l(k(\\tau\_0-\\tau))
> $$
>
> cc @marcobonici and @sofiachiarenza
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 5, 2025, 3:17am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/11 "2025-10-05T03:17:29Z")
> [@hersle](#):
>
> @ChrisRackauckas Just FYI, the ODE solving in SymBoltz relies on first solving one “background” ODE, followed by several independent “perturbation” ODEs where the “background” solution is supplied through interpolating [DataInterpolations.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/DataInterpolations) splines (as time-dependent parameters). I don’t know if that makes reverse-mode more or less complicated. It would be fantastic if that works.
That should just work. It would be better to use discrete adjoints and the interpolation of the solver if that’s what you’re trying to do though (discrete adjoints since that’s required for the solver interpolation). But DataInterpolations in this form should just work.
---
### Author:  [@Shashank](https://discourse.julialang.org/u/Shashank)
#### Post date: [October 13, 2025, 2:06pm UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/12 "2025-10-13T14:06:57Z")
I have a few questions/comments regarding getting started. I was looking at SymBoltz.jl/scripts/cmb\_source.jl, and the script does not work for me.
In line 15 I get an error at
`τs = sol1[M.τ]`
M has no member τ. Is it supposed to be M.t?
Also, on line 19:
`println("max(abs(S2-S1)) = ", maximum(abs.((S2 .- S1)[1])))`
I have to change this line for the code to run.
On line 31
`Θls = SymBoltz.los_integrate(S, ls, ks, τs, us, u′s)`
gives an error:
```julia-auto
ERROR: MethodError: no method matching los_integrate(::Tuple{…}, ::StepRange{…}, ::StepRangeLen{…}, ::Vector{…}, ::StepRangeLen{…}, ::Vector{…})
The function `los_integrate` exists, but no method is defined for this combination of argument types.
Closest candidates are:
los_integrate(::AbstractArray{T}, ::AbstractArray{T}, ::AbstractArray, ::AbstractRange, ::AbstractArray, ::AbstractRange, ::AbstractArray; integrator, verbose) where T<:Real
@ SymBoltz ~/.julia/packages/SymBoltz/kRJhS/src/spectra.jl:158
Stacktrace:
[1] top-level scope
@ REPL[36]:1
Some type information was truncated. Use `show(err)` to see complete types.
```
I am running the code using julia 1.11.0 version. Is it possible that the documentation is outdated or written with an older julia version?
I really appreciate your effort in writing this package. Thanks a lot!!!
Edit: M does have tau if I install julia version 1.12.0, but then SymBoltz does not have a member named cmb\_ks. Strange!
---
### Author:  [@hersle](https://discourse.julialang.org/u/hersle)
#### Post date: [October 15, 2025, 10:15am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/13 "2025-10-15T10:15:04Z")
> [@Shashank](#):
>
> I was looking at SymBoltz.jl/scripts/cmb\_source.jl, and the script does not work for me.
Hi @Shashank. The scripts/ directory just contains some dirty scripts I have used for debugging during development of earlier versions. It is not updated. I should probably remove it.
Please refer to the [documentation pages](https://hersle.github.io/SymBoltz.jl/stable/) for up-to-date examples. And do let me know if you want to do something not covered there 🙂
---
### Author:  [@Shashank](https://discourse.julialang.org/u/Shashank)
#### Post date: [October 15, 2025, 12:37pm UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/14 "2025-10-15T12:37:42Z")
Thanks a lot. The scripts on the documentation page work well for me, but only if I run them on version 1.12.0, but do not work on 1.11.0. The members of M are different depending on which version of julia I am using. There may be some problem with my 1.11.0 installation.
---
### Author:  [@martinmestre](https://discourse.julialang.org/u/martinmestre)
#### Post date: [October 15, 2025, 6:54pm UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/15 "2025-10-15T18:54:08Z")
Hi @hersle! Thanks! I am not expert in the field. I would need these WDM solutions for studying stellar stream sub-halo perturbations. In particular I am interested in the DM family named Relativistic Fermionic model ([https://arxiv.org/pdf/1908.10806](https://arxiv.org/pdf/1908.10806)). I will write DM to discuss further. Best.
---
### Author:  [@Shashank](https://discourse.julialang.org/u/Shashank)
#### Post date: [October 17, 2025, 11:14am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/16 "2025-10-17T11:14:35Z")
In “getting started”
```julia-auto
ls = 10:10:1000
jl = SphericalBesselCache(ls)
Dls = spectrum_cmb(:TT, prob, jl; normalization = :Dl, unit = u"μK")
plot(ls, Dls; xlabel = "l", ylabel = "l (l+1) Cₗ / 2π", label = nothing)
```
does not work. I get an error
```julia-auto
ERROR: UndefVarError: `SphericalBesselCache` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Stacktrace:
[1] top-level scope
@ REPL[3]:1
```
Is SphericalBesselCache supposed to be imported from an external module or defined by the user?
Edit: I resolved the issue. I installed SphericalFunctions.jl then uninstalled Symboltz and then reinstalled SymBoltz. Now SphericalBesselCache is there.
---
### Author:  [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)
#### Post date: [October 17, 2025, 12:46pm UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/17 "2025-10-17T12:46:26Z")
Apologies if my questions are somewhat naive, I’m out of the GR game for some years.
So the way I understand this, you solve some FLRW-style ODE for homogeneous isotropic background geometry, and then on top of that compute stability in direction of various inhomogeneities, and have some machinery to match to observational data?
Something I’d find very interesting is to look at homogeneous anisotropic background, i.e. Bianchi, especially type 8 and 9 with the fun chaotic BKL / Mixmaster attractor + particle horizons near the initial singularity / big bang. I could talk at lengths about Mixmaster, but never got to look at how its particle/causal horizons interact with inhomogenous pertubations (I’m more of an ODE than GR person and I would fail miserably if you asked me to write down equations for that).
Mixmaster may look numerically intractable on a naive glance (the heteroclinic connections can mess you up!), but is actually very tractable if you locally view it as pertubations of Bianchi 1 and 2 (and kinda 6\_0 and 7\_0 and LRS / Taub). You just need to give up “have one set of equations on one set of coordinates” and stitch together different methods in different regimes.
Any intentions on looking in that direction, ie anisotropic background geometries?
---
### Author:  [@hersle](https://discourse.julialang.org/u/hersle)
#### Post date: [October 20, 2025, 12:22pm UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/18 "2025-10-20T12:22:19Z")
> [@Shashank](#):
>
> Is SphericalBesselCache supposed to be imported from an external module or defined by the user?
>
> Edit: I resolved the issue. I installed [SphericalFunctions.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/SphericalFunctions) then uninstalled Symboltz and then reinstalled SymBoltz. Now SphericalBesselCache is there.
@Shashank I added `SphericalBesselCache` in a recent [update](https://github.com/hersle/SymBoltz.jl/releases). The reason for this error was simply that you were using the documentation for a newer version of SymBoltz than you had installed. The same is probably true what you encountered with different members of `M` above. Typing `using Pkg; Pkg.update()` would update SymBoltz and fix the problem, and reinstalling it like you did accomplishes the same. SphericalFunctions.jl has nothing to do with it, though. 🙂
---
### Author:  [@hersle](https://discourse.julialang.org/u/hersle)
#### Post date: [October 20, 2025, 12:49pm UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/19 "2025-10-20T12:49:52Z")
> [@foobar\_lv2](#):
>
> So the way I understand this, you solve some FLRW-style ODE for homogeneous isotropic background geometry, and then on top of that compute stability in direction of various inhomogeneities, and have some machinery to match to observational data?
@foobar_lv2 Yes, it solves ODEs for a linearly perturbed FLRW metric, i.e. first-order corrections to a homogeneous and isotropic universe. SymBoltz produces theoretical predictions for a given model, which can be compared to observations (observational _data_ is outside SymBoltz’ scope).
> [@foobar\_lv2](#):
>
> Something I’d find very interesting is to look at homogeneous anisotropic background, i.e. Bianchi, especially type 8 and 9 with the fun chaotic BKL / Mixmaster attractor + particle horizons near the initial singularity / big bang. I could talk at lengths about Mixmaster, but never got to look at how its particle/causal horizons interact with inhomogenous pertubations (I’m more of an ODE than GR person and I would fail miserably if you asked me to write down equations for that).
>
> Mixmaster may look numerically intractable on a naive glance (the heteroclinic connections can mess you up!), but is actually very tractable if you locally view it as pertubations of Bianchi 1 and 2 (and kinda 6\_0 and 7\_0 and LRS / Taub). You just need to give up “have one set of equations on one set of coordinates” and stitch together different methods in different regimes.
>
> Any intentions on looking in that direction, ie anisotropic background geometries?
I don’t have any hands-on familiarity with these models, but it sounds interesting. As long as they have a “background+perturbations”-structure that is not too different from (say) ΛCDM, I think they should be possible to implement in SymBoltz’ framework. Up until now my main priority has been to get the code to work and agree well with existing standard codes for close-to standard models.
---
### Author:  [@hersle](https://discourse.julialang.org/u/hersle)
#### Post date: [March 16, 2026, 11:24am UTC](https://discourse.julialang.org/t/ann-symboltz-jl-a-symbolic-numeric-approximation-free-and-differentiable-linear-einstein-boltzmann-solver-for-cosmology/132819/20 "2026-03-16T11:24:25Z")
Update: The code paper for SymBoltz has passed peer review and is now published in [Astronomy & Astrophysics](http://dx.doi.org/10.1051/0004-6361/202557450) (the [arXiv preprint](https://arxiv.org/abs/2509.24740) is updated and fully equivalent).
The code has also received some updates. In particular, with [new ModelingToolkit updates](https://github.com/SciML/ModelingToolkit.jl/pull/4037/), it now generates analytical, sparse and allocation-free Jacobians for the perturbation equations, which is very powerful when combined with the Rodas5P solver in OrdinaryDiffEq. This makes solving the expensive perturbation equations [10x faster than in CLASS](https://arxiv.org/pdf/2509.24740#section.4) when CLASS does not use approximations, and about 1x-2x (as fast) when CLASS uses all of its approximation tricks, while SymBoltz never uses approximations.
It is also updated to support ModelingToolkit v11, which has made the symbolic-to-numeric code generation faster.
It is a good time to try the package. Please let me know if you have any questions or run into any issues!