Observable quantities

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

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.

spectrum_matter([species,] sol::CosmologySolution, k, τ = sol[τ][end])

Compute the power spectrum

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

of the total gauge-invariant overdensity

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

for one or more species (default: m for matter) at wavenumber(s) k and conformal time(s) τ (final, if omitted) from the solution sol.

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

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

spectrum_matter([species,] prob::CosmologyProblem, k::NTuple{2, Number}, τs = nothing; sourceopts = (atol = 4.0, rtol = 4e-3), coarse_length = 9, kwargs...)

Compute the matter power spectrum on the interval $k$ with adaptively chosen wavenumbers. Returns wavenumbers and power spectrum values.

The interval is first divided into a grid with coarse_length logarithmically spaced wavenumbers. It is then adaptively refined with the absolute and relative tolerances in sourceopts.

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.

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

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.

spectrum_cmb(modes::AbstractVector{<:Symbol}, prob::CosmologyProblem, jl::SphericalBesselCache; normalization = :Cl, unit = nothing, kτ0s = 0.1*jl.l[begin]:2π/2:10*jl.l[end], xs = 0.0:0.0008:1.0, l_limber = 50, integrator = TrapezoidalRule(), bgopts = (alg = Rodas4P(), reltol = 1e-9, abstol = 1e-9), ptopts = (alg = KenCarp4(), reltol = 1e-8, abstol = 1e-8), sourceopts = (rtol = 1e-3, atol = 0.9), coarse_length = 9, thread = true, verbose = false, kwargs...)

Compute angular CMB power spectra $Cₗᴬᴮ$ at angular wavenumbers ls from the cosmological problem prob. The requested modes are specified as a vector of symbols in the form :AB, where A and B are T (temperature), E (E-mode polarization) or ψ (lensing). 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. Returns a matrix of $Cₗ$ if normalization is :Cl, or $Dₗ = l(l+1)/2π$ if normalization is :Dl.

The lensing line-of-sight integral uses the Limber approximation for l ≥ l_limber.

Source functions are computed on a $k$-grid that is adaptively refined from an initial grid with size coarse_length. The refinement criterion is controlled with sourceopts.

Examples

using SymBoltz, Unitful
M = ΛCDM()
pars = parameters_Planck18(M)
prob = CosmologyProblem(M, pars)

ls = 10:10:1000
jl = SphericalBesselCache(ls)
modes = [:TT, :TE, :ψψ, :ψT]
Dls = spectrum_cmb(modes, prob, jl; normalization = :Dl, unit = u"μK")

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

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

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, ξ).

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.

stddev_matter(sol::CosmologySolution, R)

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

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.

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_grid(sol::CosmologySolution, Ss::AbstractVector, τs; thread = true)

Evaluate and return source functions $S(τ,k)$ from the solution sol. The source functions are given by symbolic expressions Ss, and evaluated on a grid with conformal times τs and the same wavenumbers as sol is solved for.

source_grid(Ss_coarse::AbstractArray, ks_coarse, ks_fine; ktransform = identity, thread = true)

Interpolate values Ss_coarse of source functions $S(τ,k)$ from a coarse wavenumber grid ks_coarse to a fine grid ks_fine. The interpolation is linear in ktransform(k) (e.g. identity for interpolation in $k$ or log for interpolation in $\ln k$. Conformal times are unchanged.

source_grid(prob::CosmologyProblem, Ss::AbstractArray, τs, ks; bgopts = (), ptopts = (), thread = true, verbose = false)

Compute and evaluate source functions $S(τ,k)$ with symbolic expressions Ss on a grid with conformal times τs and wavenumbers ks from the problem prob.

The options bgopts and ptopts are passed to the background and perturbation solves.

source_grid_adaptive(prob::CosmologyProblem, Ss::AbstractVector, τs, ks[, bgsol]; bgopts = (), ptopts = (), ktransform = (identity, identity), sort = true, thread = true, verbose = false, kwargs...)

Adaptively compute and evaluate source functions $S(τ,k)$ with symbolic expressions Ss on a grid with fixed conformal times τs, but adaptively refined grid of wavenumbers from the problem prob. The source functions are first evaluated on the (coarse) initial grid ks. Each subinterval $(k₁, k₂)$ of ks is then adaptively refined until the linear interpolation $Sᵢ = (S(k₁)+S(k₂))/2$ to the midpoint $k=(k₁+k₂)/2$ approximates the actual value $S(k)$ there within some tolerance. The comparison $Sᵢ ≈ S$ is done with isapprox(Sᵢ, S; kwargs...), where S and Sᵢ are vectors with the (conformal) timeseries of the source function for that wavenumber. It receives the keyword arguments kwargs passed to this function, so atol, rtol and/or norm can be specified to tune the tolerance.

If τs is nothing, the source function is evaluated at the final time only (today).

Returns the refined wavenumbers sorted in ascending order and a grid with the corresponding source function values. If not sort, the wavenumbers and source function values are instead left in the order in which they were inserted in the refinement process.

The options bgopts and ptopts are passed to the background and perturbation solves.

los_integrate(Ss::AbstractArray{T, 3}, ls::AbstractVector, τs::AbstractVector, ks::AbstractVector, jl::SphericalBesselCache; l_limber = typemax(Int), integrator = TrapezoidalRule(), thread = true, verbose = false) where {T <: Real}

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

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

over the source function values Ss against the spherical Bessel functions $jₗ(x)$ cached in jl. The element Ss[i,j,k] holds the source function value $Sᵢ(kⱼ, τₖ)$. The limber approximation

\[Iₗ ≈ √(π/(2l+1)) S(k,τ₀-(l+1/2)/k)\]

is used for l ≥ l_limber.