Numerical evolution of Einstein-Boltzmann equations in cosmology

[chirag1]

TL;DR Summary

How to find the cutoff variable to halt the integration of the cosmological perturbations top determine the evolution of different k-modes?

I'm trying to numerically evolve the Einstein-Boltzmann equations for cold dark matter perturbations using Runge-Kutta method of the fourth order. There are 5 standard equations:
$$
\begin{align}
\dot{\Theta}_{r,0}+k\Theta_{r,1}&=-\dot{\Phi} \\
\dot{\Theta}_{r,1}+\frac{k}{3}\Theta_{r,0} & =\frac{-k}{3}\Phi \\
\dot{\delta}+ikv &= -3\dot{\Phi} \\
\dot{v}+\frac{\dot{a}}{a}v &= ik\Phi \\
\dot{\Phi}&=\frac{1}{3\dot{a}}\frac{3H_{0}^{2}}{2}\left(\Omega_{m}\delta+4\Omega_{r}\Theta_{r,0}a^{-1}\right)-ak^{2}\Phi-\frac{\dot{a}}{a}\Phi
\end{align}
$$The problem is, we cannot integrate them all the way to the present as radiation moments are difficult to track at late times and especially so for small scale (large k) modes. The solution to this is to find a cutoff time at which we halt the integration, discard the radiation perturbations and restart the integration. I'm facing the issue of how to obtain an expression for this cutoff time here, which depends on the k-mode. I'm more surprised by the lack of presented solutions for this standard problem (this numerical integration task is given as an textbook exercise in Chapter 8 (ex. 8.2) Modern Cosmology-Dodelson 2nd edition and 1st edition also which was more than 15 years ago, but there is no solution to this textbook exercise as well!) in literature or papers.

I've tried a lot to find something but everyone is seemingly not tackling these 5 equations and taking a different approach. But for my project, I've to work on these 5 equations only. The closest I got to something was Florian Borchers' thesis: https://www.imperial.ac.uk/media/im...ations/2010/Florian-Borchers-Dissertation.pdf where they give an expression for cutoff conformal time (page 32) but give no explanation. That expression is:
$$
\eta_{\text{stop}} = \eta_{\text{today}} - \frac{2}{3}log(100k/h)
$$
They actually use stepperdopr853 method for integration and conformal time as their integration variable, while I use RK4 and scale factor. I've tried to account for it and take help of chatgpt as well and all literature that I could find but in vain. I'm very stuck.

 

[renormalize]

$$
\begin{align}
\dot{\Theta}_{r,0}+k\Theta_{r,1}&=-\dot{\Phi} \\
\dot{\Theta}_{r,1}+\frac{k}{3}\Theta_{r,0} & =\frac{-k}{3}\Phi \\
\dot{\delta}+ikv &= -3\dot{\Phi} \\
\dot{v}+\frac{\dot{a}}{a}v &= ik\Phi \\
\dot{\Phi}&=\frac{1}{3\dot{a}}\frac{3H_{0}^{2}}{2}\left(\Omega_{m}\delta+4\Omega_{r}\Theta_{r,0}a^{-1}\right)-ak^{2}\Phi-\frac{\dot{a}}{a}\Phi
\end{align}
$$

For the non-expert it would be helpful if you could define the physical meaning of each variable that appears in your equations, as well as the definition of the "dot": differentiation w.r.t. to proper time? cosmological time? conformal time? Also, at what value of the time are you starting your integration and what are the initial conditions of each variable at the start?

 

[chirag1]

For the non-expert it would be helpful if you could define the physical meaning of each variable that appears in your equations, as well as the definition of the "dot": differentiation w.r.t. to proper time? cosmological time? conformal time? Also, at what value of the time are you starting your integration and what are the initial conditions of each variable at the start?

Yes.
I don't know much of the physical meaning of them myself yet, but I'll try to explain what I know.

$\Theta_{r,0}$ is the monopole radiation term.

It corresponds to the fractional perturbation in theangle-averaged photon flux at a given position x and time t.

$\Theta_{r,1}$ is the dipole radiation term.
$\delta$ is the dark matter density perturbation and $v$ is the bulk velocity perturbation of the dark matter.
$\Phi$ is the gravitational potential which is taken as the perturbation in the metric.

$\Omega_i$ represents the density of the species $i$.
$k$ is the wavenumber and $H_0$ is the Hubble constant value today.

The dot represents differentiation with respect to conformal time $\eta$. We can change our integration variable from this $\eta$ to scale factor $a$ using
$$\frac{da}{d\eta}=a^2H$$
where $H = \frac{da/dt}{a}$ is the Hubble parameter.
Then we can further change the variable to $log_{10}a$ which is what I'm using as my integration variable.

I'm starting the integration at $a=10^{-8}$ and trying to evolve to present $a=1$.
The initial conditions are for the variables are given by inflation-induced adiabatic modes :
$$
\begin{align}
\Theta_{r,0} &= \frac{1}{2}\Phi\\
\Theta_{r,1} &= -\frac{k}{6aH}\Phi \\
\delta &= \frac{3}{2}\Phi\\
v &= \frac{ik}{2aH}\Phi\\
\end{align}
$$

For $\Phi$ itself we normalise it to 1 as the initial value, which doesn't matter much as all the variables will be scaled accordingly.

The problem in integration happens at late times specially for small scale i.e. large $k$ modes ($k \geq 0.01 Mpc^{-1}$).