Observable quantities

This page shows observable quantities that can be derived from solutions of the Einstein-Boltzmann system, such as power spectra and distances.

Primordial power spectra

SymBoltz.spectrum_primordial — Function

spectrum_primordial(k, h, As, ns=1.0; kp = k_dimensionless(0.05 / u"Mpc", h))

Compute the primordial power spectrum

\[P₀(k) = 2π² Aₛ (k/kₚ)^{nₛ-1} / k³\]

with spectral amplitude As, spectral index ns and pivot scale wavenumber kp at the wavenumber(s) k.

source

Example

using SymBoltz, Unitful, UnitfulAstro, Plots
M = SymBoltz.ΛCDM()
pars = Dict(M.g.h => 0.7, M.I.As => 2e-9, M.I.ns => 0.95)
ks = 10 .^ range(-5, +1, length=100) / u"Mpc"
Ps = spectrum_primordial(ks, M, pars)
plot(log10.(ks*u"Mpc"), log10.(Ps/u"Mpc^3"); xlabel = "log10(k/Mpc⁻¹)", ylabel = "log10(P/Mpc³)")

Matter power spectra

SymBoltz.spectrum_matter — Function

spectrum_matter(sol::CosmologySolution, k, τ = sol[τ][end]; species = [:c, :b, :h])

Compute the power spectrum

\[P(k,τ) = P₀(k) |Δ(k,τ)|²\]

of the total gauge-invariant overdensity

\[Δ = δ + (3ℰ/k²) θ = (∑ₛδρ)/(∑ₛρₛ) + (3ℰ/k²) (∑ₛ(ρₛ+Pₛ)θₛ) / (∑ₛ(ρₛ+Pₛ))\]

for the given species at wavenumber(s) k and conformal time(s) tau (final, if omitted) from the solution sol. By default, the species are cold dark matter, baryons and massive neutrinos, which are matter-like at late times in the ΛCDM model.

source

spectrum_matter(prob::CosmologyProblem, k, τ = nothing; species = [:c, :b, :h], kwargs...)

Solve the problem prob with exact wavenumber(s) k, and then compute the power spectrum with the solution sol.

source

SymBoltz.spectrum_matter_nonlinear — Function

spectrum_matter_nonlinear(sol::CosmologySolution, k)

Compute the nonlinear matter power spectrum from the cosmology solution sol at wavenumber(s) k using halofit implemented in MatterPower.jl.

source

Example

using SymBoltz, Unitful, UnitfulAstro, Plots
M = SymBoltz.ΛCDM()
pars = SymBoltz.parameters_Planck18(M)
prob = CosmologyProblem(M, pars)
ks = 10 .^ range(-5, +2, length=200) / u"Mpc"
sol = solve(prob, ks)

# Linear power spectrum
Ps = spectrum_matter(sol, ks)
plot(log10.(ks*u"Mpc"), log10.(Ps/u"Mpc^3"); xlabel = "log10(k/Mpc⁻¹)", ylabel = "log10(P/Mpc³)", label = "linear (SymBoltz)")

# Nonlinear power spectrum (from halofit)
Ps = spectrum_matter_nonlinear(sol, ks)
plot!(log10.(ks*u"Mpc"), log10.(Ps/u"Mpc^3"); label = "non-linear (halofit)", legend_position = :bottomleft)

CMB power spectra

SymBoltz.spectrum_cmb — Function

spectrum_cmb(ΘlAs::AbstractMatrix, ΘlBs::AbstractMatrix, P0s::AbstractVector, ls::AbstractVector, ks::AbstractVector; integrator = TrapezoidalRule(), normalization = :Cl)

Compute the angular power spectrum

\[Cₗᴬᴮ = (2/π) ∫\mathrm{d}k \, k² P₀(k) Θₗᴬ(k,τ₀) Θₗᴮ(k,τ₀)\]

for the given ls. If normalization == :Dl, compute $Dₗ = Cₗ l (l+1) / 2π$ instead.

source

spectrum_cmb(modes::AbstractVector, prob::CosmologyProblem, ls::AbstractVector; normalization = :Cl, unit = nothing, Δkτ0 = 2π/2, Δkτ0_S = 8.0, kτ0min = 0.1*ls[begin], kτ0max = 3*ls[end], u = (τ->tanh(τ)), u⁻¹ = (u->atanh(u)), Nlos = 768, integrator = TrapezoidalRule(), bgopts = (alg = Rodas4P(), reltol = 1e-9, abstol = 1e-9), ptopts = (alg = KenCarp4(), reltol = 1e-8, abstol = 1e-8), thread = true, kwargs...)

Compute the CMB power spectra modes (:TT, :EE, :TE or an array thereof) $C_l^{AB}$'s at angular wavenumbers ls from the cosmological solution sol. If unit is nothing the spectra are of dimensionless temperature fluctuations relative to the present photon temperature; while if unit is a temperature unit the spectra are of dimensionful temperature fluctuations.

source

spectrum_cmb(modes::AbstractVector, prob::CosmologyProblem, ls::AbstractVector, ls_interpolate_from::AbstractVector; kwargs...)

Same, but compute the spectrum properly only for ls_interpolate_from and then interpolate the results to all ls.

source

Example

using SymBoltz, Unitful, UnitfulAstro, Plots
M = SymBoltz.ΛCDM()
pars = SymBoltz.parameters_Planck18(M)
prob = CosmologyProblem(M, pars)
ls = 10:5:1500

Dls = spectrum_cmb([:TT, :EE, :TE], prob, ls; normalization = :Dl, unit = u"μK")
pTT = plot(ls, Dls[:, 1]; ylabel = "Dₗᵀᵀ")
pEE = plot(ls, Dls[:, 2]; ylabel = "Dₗᴱᴱ")
pTE = plot(ls, Dls[:, 3]; ylabel = "Dₗᵀᴱ", xlabel = "l")
plot(pTT, pEE, pTE, layout = (3, 1), size = (600, 700), legend = nothing)

Two-point correlation function

SymBoltz.correlation_function — Function

correlation_function(sol::CosmologySolution; N = 2048, spline = true)

Compute the two-point correlation function in real space by Fourier transforming the matter power spectrum of sol with N points the FFTLog algorithm implemented in TwoFAST. Returns N radii and correlation function values (e.g. r, ξ).

source

Example

using SymBoltz, Unitful, UnitfulAstro, Plots
M = SymBoltz.ΛCDM()
pars = SymBoltz.parameters_Planck18(M)
prob = CosmologyProblem(M, pars)
ks = 10 .^ range(-5, +3, length=300) / u"Mpc"
sol = solve(prob, ks)
rs, ξs = correlation_function(sol)
rs = rs / (SymBoltz.k0*sol.bg.ps[:h]) * u"Mpc" # TODO: auto units
plot(rs, @. ξs * rs^2; xlims = (0, 200), xlabel = "r", ylabel = "r² ξ")

Matter density fluctuations

SymBoltz.variance_matter — Function

variance_matter(sol::CosmologySolution, R)

Compute the variance $⟨δ²⟩$ of the linear matter density field with a top-hat filter with radius R. Wraps the implementation in MatterPower.jl.

source

SymBoltz.stddev_matter — Function

stddev_matter(sol::CosmologySolution, R)

Compute the standard deviation $√(⟨δ²⟩)$ of the linear matter density field with a top-hat filter with radius R.

source

using SymBoltz, Unitful, UnitfulAstro, Plots
M = SymBoltz.ΛCDM()
pars = SymBoltz.parameters_Planck18(M)
prob = CosmologyProblem(M, pars)
ks = 10 .^ range(-5, +3, length=300) / u"Mpc"
sol = solve(prob, ks)

h = sol[M.g.h]
Rs = 10 .^ range(0, 2, length=100) * u"Mpc"
σs = stddev_matter.(sol, Rs)
plot(log10.(Rs/(u"Mpc"/h)), log10.(σs); xlabel = "lg(R / (Mpc/h))", ylabel = "lg(σ)", label = nothing)

R8 = 8 * u"Mpc"/h
σ8 = stddev_matter(sol, R8)
scatter!((log10(R8/(u"Mpc"/h)), log10(σ8)), series_annotation = text("  σ₈ = $(round(σ8; digits=3))", :left), label = nothing)

Luminosity distance

SymBoltz.distance_luminosity — Function

distance_luminosity(sol::CosmologySolution, ivs = sol.bg.t, τ0 = sol[sol.prob.M.τ0])

Compute luminosity distances

\[d_L = \frac{r}{a} = \chi \, \mathrm{sinc} (\sqrt{K} (τ₀-τ)),\]

at the independent variable values ivs relative to the (present) time τ0.

source

using SymBoltz, Plots
M = SymBoltz.RMΛ(K = SymBoltz.curvature(SymBoltz.metric()))
pars = Dict(
    M.r.Ω₀ => 5e-5,
    M.m.Ω₀ => 0.3,
    M.K.Ω₀ => 0.1,
    M.r.T₀ => NaN,
    M.g.h => 0.7
)
prob = CosmologyProblem(M, pars)
sol = solve(prob)

zs = 0.0:1.0:10.0
τs = SymBoltz.timeseries(sol, M.g.z, zs) # times at given redshifts
dLs = SymBoltz.distance_luminosity(sol, τs) / SymBoltz.Gpc
plot(zs, dLs; marker=:dot, xlabel="z", ylabel="dL / Gpc", label=nothing)

Sound horizon (BAO scale)

SymBoltz.sound_horizon — Function

sound_horizon(sol::CosmologySolution)

Cumulatively integrate the sound horizon

\[ rₛ(τ) = ∫_0^τ dτ cₛ = ∫_0^τ \frac{dτ}{√(3(1+3ρ_b/4ρ_γ))}\]

to the time steps of the solution sol.

source

using SymBoltz, Plots
M = SymBoltz.ΛCDM()
pars = SymBoltz.parameters_Planck18(M)
prob = CosmologyProblem(M, pars)
sol = solve(prob)
τs = sol[M.τ]
rs = sound_horizon(sol)
plot(τs, rs; xlabel = "τ / H₀⁻¹", ylabel = "rₛ / (c/H₀)")

Line-of-sight integration

SymBoltz.los_integrate — Function

los_integrate(Ss::AbstractMatrix{T}, ls::AbstractVector, τs::AbstractVector, ks::AbstractVector, Rl::Function; integrator = TrapezoidalRule(), verbose = false) where {T <: Real}

For the given ls and ks, compute the line-of-sight-integrals

\[Iₗ(k) = ∫dτ S(k,τ) Rₗ(k(τ₀-τ))\]

over the source function values Ss against the radial functions Rl (e.g. the spherical Bessel functions $jₗ(x)$). The element Ss[i,j] holds the source function value $S(kᵢ, τⱼ)$.

source