Getting started

Install SymBoltz.jl and load it with

using SymBoltz

1. Create the model

The first step is to define our cosmological model. This is a symbolic representation of the variables and equations that describe the various components in the universe, such as the theory of gravity (like general relativity) and particle species (like the cosmological constant, cold dark matter, photons and baryons).

To get started, we will simply load the standard ΛCDM model:

M = ΛCDM()
hierarchy(M; describe = true)

ΛCDM: Standard cosmological constant and cold dark matter cosmological model
├─ g: Spacetime FLRW metric in Newtonian gauge
├─ G: General relativity gravity
├─ γ: Photon radiation
├─ ν: Massless neutrinos
├─ c: Cold dark matter
├─ b: Baryonic matter
│  ├─ rec: Baryon-photon recombination thermodynamics (RECFAST)
│  ├─ rei1: Reionization with tanh-like (activation function) contribution to the free electron fraction
│  └─ rei2: Reionization with tanh-like (activation function) contribution to the free electron fraction
├─ h: Massive neutrino
├─ Λ: Cosmological constant
└─ I: Harrison-Zel'dovich inflation

As shown, the model is structured as a hierarchy of the physical components. Each of these components contains a block of self-contained variables and equations that are independent from the other components, making up interchangeable modules of the entire Einstein-Boltzmann system. A full model joins several such incomplete blocks into a complete set of equations for the entire Einstein-Boltzmann system.

The hierarchical structure can be inspected interactively (with TAB-completion) in the Julia REPL by evaluating M, M.G, M.G.ρ and so on. For example, to see all equations for the theory of gravity:

equations(M.G)

\[ \begin{align*} \frac{\mathrm{d} a\left( \tau \right)}{\mathrm{d}\tau} &= \left( a\left( \tau \right) \right)^{2} \sqrt{\frac{8}{3} \rho\left( \tau \right) \pi} \\ \mathtt{F_1}\left( \tau \right) &= 0 \\ \mathtt{F_2}\left( \tau \right) &= \frac{\mathrm{d}}{\mathrm{d}\tau} \frac{\mathrm{d} a\left( \tau \right)}{\mathrm{d}\tau} + \frac{ - \left( \frac{\mathrm{d} a\left( \tau \right)}{\mathrm{d}\tau} \right)^{2} \left( 1 + \frac{ - 3 P\left( \tau \right)}{\rho\left( \tau \right)} \right)}{2 a\left( \tau \right)} \\ \frac{\mathrm{d}}{\mathrm{d}\tau} \Phi\left( \tau, k \right) &= \frac{ - \frac{4}{3} \left( a\left( \tau \right) \right)^{2} \mathtt{\delta\rho}\left( \tau, k \right) \pi}{\mathscr{E}\left( \tau \right)} + \frac{ - k^{2} \Phi\left( \tau, k \right)}{3 \mathscr{E}\left( \tau \right)} - \Psi\left( \tau, k \right) \mathscr{E}\left( \tau \right) \\ k^{2} \left( \Phi\left( \tau, k \right) - \Psi\left( \tau, k \right) \right) &= 12 \left( a\left( \tau \right) \right)^{2} \Pi\left( \tau, k \right) \pi \\ \mathtt{\dot{\Psi}}\left( \tau, k \right) &= \frac{\mathrm{d}}{\mathrm{d}\tau} \Psi\left( \tau, k \right) \\ \mathtt{\dot{\Phi}}\left( \tau, k \right) &= \frac{\mathrm{d}}{\mathrm{d}\tau} \Phi\left( \tau, k \right) \end{align*} \]

2. Solve the problem

Next, we create a numerical representation of the cosmological problem we want to solve. This splits the full symbolic model into computational stages (like the background, thermodynamics and perturbations), assigns input values to parameters and defines any parameters that are solved for with the shooting method by matching conditions at the final time.

using Unitful, UnitfulAstro # for interfacing without internal code units
pars = Dict(
    M.γ.T₀ => 2.7,
    M.c.Ω₀ => 0.27,
    M.b.Ω₀ => 0.05,
    M.ν.Neff => 3.0,
    M.g.h => 0.7,
    M.b.YHe => 0.25,
    M.h.m_eV => 0.06,
    M.I.ln_As1e10 => 3.0,
    M.I.ns => 0.95
)
ks = 10 .^ range(-5, 1, length=500) / u"Mpc"
prob = CosmologyProblem(M, pars)

Cosmology problem for model ΛCDM
Stages:
  Background: 5 unknowns, 0 splines
  Perturbations: 61 unknowns, 5 splines
Parameters & initial conditions:
  c₊Ω₀ = 0.27
  b₊Ω₀ = 0.05
  γ₊T₀ = 2.7
  b₊YHe = 0.25
  h₊m_eV = 0.06
  h = 0.7
  I₊ns = 0.95
  ν₊Neff = 3.0
  I₊ln_As1e10 = 3.0

Finally, we can simply solve the problem:

sol = solve(prob, ks) # or just solve(prob) to solve only the background

Cosmology solution for model ΛCDM
Stages:
  Background: return code Terminated; solved with Rodas4P; 1304 points
  Perturbations: return codes Success; solved with KenCarp4; 79-823 points; x500 k ∈ [0.042827494000000015, 42827.49400000001] H₀/c (interpolation in-between)

using MKL

3. Use the solution

You are now free to do whatever you want with the solution object. For example, to get the time points used by the solver and corresponding values of the scale factor $a(τ)$:

τs = sol[M.τ]
as = sol(M.g.a, τs)

1304-element Vector{Float64}:
 1.0e-8
 1.0027207449885391e-8
 1.0037274056317198e-8
 1.0056767150657421e-8
 1.0076260245653028e-8
 1.0098167994331648e-8
 1.0121203640482273e-8
 1.0146012166380299e-8
 1.0172211415325421e-8
 1.019994871640612e-8
 ⋮
 0.6786180270048812
 0.7124705016246922
 0.7482788734312426
 0.7862769150744753
 0.8267384408965704
 0.8699856552095218
 0.9163997004561728
 0.9664341162138457
 1.0

Similarly, to get $\Phi(k,τ)$ for the 500 wavenumbers we solved for at the same times:

Φs = sol(M.g.Φ, τs, ks)

1304×500 Matrix{Float64}:
 0.704214  0.704214  0.704214  …  0.704214    0.704214    0.704214
 0.704214  0.704214  0.704214     0.704213    0.704213    0.704213
 0.704214  0.704214  0.704214     0.704212    0.704212    0.704212
 0.704214  0.704214  0.704214     0.704211    0.704211    0.704211
 0.704214  0.704214  0.704214     0.704211    0.704211    0.70421
 0.704214  0.704214  0.704214  …  0.70421     0.70421     0.704209
 0.704214  0.704214  0.704214     0.704209    0.704209    0.704208
 0.704214  0.704214  0.704214     0.704208    0.704208    0.704207
 0.704214  0.704214  0.704214     0.704207    0.704206    0.704206
 0.704214  0.704214  0.704214     0.704206    0.704205    0.704205
 ⋮                             ⋱                          
 0.544375  0.544375  0.544374  …  3.07457e-5  2.92565e-5  2.78386e-5
 0.537644  0.537644  0.537644     3.03629e-5  2.88922e-5  2.7492e-5
 0.530279  0.530279  0.530279     2.99444e-5  2.8494e-5   2.71131e-5
 0.522239  0.522239  0.522239     2.9488e-5   2.80597e-5  2.66998e-5
 0.513482  0.513482  0.513482     2.89912e-5  2.7587e-5   2.625e-5
 0.503968  0.503968  0.503968  …  2.84518e-5  2.70737e-5  2.57616e-5
 0.493655  0.493655  0.493655     2.78674e-5  2.65176e-5  2.52325e-5
 0.482503  0.482503  0.482502     2.72358e-5  2.59166e-5  2.46606e-5
 0.475039  0.475039  0.475039     2.68133e-5  2.55146e-5  2.4278e-5

You could plot this with using Plots; plot(log10.(as), transpose(Φs)), but this is more convenient with the included plot recipe:

using Plots
ks_plot = [1e-3, 1e-2, 1e-1, 1e-0] / u"Mpc"
plot(sol, log10(M.g.a), M.g.Φ, ks_plot) # lg(a) vs. Φ for 4 wavenumbers

We can also calculate the matter power spectrum:

Ps = spectrum_matter(sol, ks)
plot(log10.(ks/u"1/Mpc"), log10.(Ps/u"Mpc^3"); xlabel = "lg(k/Mpc⁻¹)", ylabel = "lg(P/Mpc³)", label = nothing)

Similarly, we can calculate the angular CMB (TT) power spectrum:

ls = 10:10:1000
Dls = spectrum_cmb(:TT, prob, ls; normalization = :Dl, unit = u"μK")
plot(ls, Dls; xlabel = "l", ylabel = "l (l+1) Cₗ / 2π", label = nothing)

And here is a condensed plot with several quantities:

p = plot(layout=(3, 3), size=(900, 700), tickfontsize=6, labelfontsize=6, legendfontsize=5)
plot!(p[1], sol, log10(M.g.a), [M.b.ρ, M.c.ρ, M.γ.ρ, M.ν.ρ, M.h.ρ, M.Λ.ρ] ./ M.G.ρ)
plot!(p[2], sol, log10(M.g.a), [M.b.w, M.c.w, M.γ.w, M.ν.w, M.h.w, M.Λ.w])
plot!(p[3], sol, log10(M.g.a), log10(M.g.E))
plot!(p[4], sol, log10(M.g.a), [M.b.rec.XHe⁺⁺, M.b.rec.XHe⁺, M.b.rec.XH⁺, M.b.Xe])
plot!(p[5], sol, log10(M.g.a), log10.([M.γ.T, M.b.T] ./ M.γ.T₀))
plot!(p[6], sol, log10(M.g.a), log10(abs(M.b.κ)))
plot!(p[7], sol, log10(M.g.a), [M.g.Φ, M.g.Ψ], ks_plot)
plot!(p[8], sol, log10(M.g.a), log10.(abs.([M.b.δ, M.c.δ, M.γ.δ, M.ν.δ, M.h.δ])), ks_plot; klabel = false)
plot!(p[9], sol, log10(M.g.a), log10.(abs.([M.b.θ, M.c.θ, M.γ.θ, M.ν.θ, M.h.θ])), ks_plot; klabel = false)