Table of contents

Back to the main page 1. Overview: Milestone III 2. Theoretical background 3. What you have to do 4. Testing your code 5. Appendix

Overview: Milestone III

The topic of the third milestone of this project is the evolution of structure in the universe: how did small fluctuations in the baryon-photon-dark-matter fluid grow from shortly after inflation until today? The ultimate goal of this part is to construct two-dimensional functions - of time and Fourier scale, $x$ and $k$ - for each of the main physical quantities of interest, $\Phi(x,k)$, $\Psi(x,k)$, $\delta_{\rm CDM}(x,k)$, $\delta_b(x,k)$, $v_{\rm CDM}(x,k)$, $v_b(x,k)$, $\Theta_\ell(x,k)$.

The deliverables are the following. There is 55 points availiable for the actual deliverables and 45 points for how this is presented/discussed in the report:

  • The source code used for the evaluation. This should be possible for us to easily compile and run. A common way students have chosen to do this in previous years is to make a GitHub account, upload your code there and just provide the link. This is no requirement, you can just as easily include this together with your report in Devilry.
  • A copy of your paper. This should have a section about this milestone which should contain: introduction (what we do and why this is relevant for the project), a short theory section with all relevant equations, tests showing the code gives expected results, the actual plots/results asked for and any other things you deem relevant and finally a proper description and discussion of the results tying it together with the physics.
  • Implement all the required functions in the template solve the ODE system for the perturbations.
  • For the plots below choose 3-4 values of the wavenumber $k$ that shows modes that enter the horizon in different eras. We want to see the three main regimes: large-scales (i.e., unaffected by causal physics), small-scales (i.e., early oscillations with subsequent damping), and intermediate scales (i.e., scales that have just undergone a few oscillations).
  • Plot of the density perturbations $\delta_\gamma = 4\Theta_0$, $\delta_{\rm CDM}$ and $\delta_b$ (and $\delta_\nu = 4\mathcal{N}_0$ if you include neutrinos) (20p).
  • Plot of the velocity perturbations $v_\gamma=-3\Theta_1$, $v_{\rm CDM}$ and $v_b$ (and $v_\nu=-3\mathcal{N}_1$ if you include neutrinos) (20p).
  • Plot of the photon quadrupole $\Theta_2$ (and $\mathcal{N}_2$ if you include neutrinos) (5p).
  • Plot of the gravitational potential $\Phi$ and the sum $\Phi+\Psi$ (10p).
  • Plot of the polarization multipoles $\Theta_{P0}$, $\Theta_{P1}$ and $\Theta_{P2}$ if you include polarization.

The project splits into two branches, one for Ph.D. students and one for Master students. The difference is that Ph. D. students have to consider also neutrinos and polarization, while Master students only need to take into account photons, baryons and dark matter, and only temperature fluctuations. So make sure to choose the one appropriate for you. (But of course, if you are a Master student and want to go for the more advanced problem, you are more than welcome to do so)

There might be typos in the equations here so make sure you double-check with those in Callin (2006).

A short introduction to the milestone can be found in this PDF (keynote) shown below.

Page: /

Theoretical background

See the lecture notes for the theory you should know for this milestone.

During the lectures, we have derived (or will derive) the linearized Einstein and Boltzmann equations for photons, baryons and dark matter, and their respective inflationary initial conditions. The task of the current part of the project is to solve these equations numerically. The good news is that the numerical solution of these equations follows in exactly the same path as when solving for instance the Peebles' equation or the equation for the conformal time. The bad news is that the expressions for the equations are somewhat more complicated. But if one is just a little careful about typing these in correctly, everything should work just fine.

But before we write down the equations, there are a few issues that should be pointed out. First, at early times the optical depth, $\tau$, is very large. The ODE we solve for the perturbations becomes numerically unstable in this regime and we must deal with this. Large $\tau$ means that electrons at a given place only observe temperature fluctions that are very nearby. This, in turn, implies that it will only see very smooth fluctuations, since the full system is in thermodynamic equilibrium, and all gradients are efficiently washed out. The only relevant quantities in this regime are therefore 1) the monopole, $\Theta_0$, which measures the mean temperature at the position of the electron, 2) the dipole, $\Theta_1$, which is given by the velocity of the fluid due to the Doppler effect, and 3) the quadrupole, $\Theta_2$, which is the only relevant source of polarization signals. The regime where this is the case is called tight coupling.

At later times, though, the fluid becomes thinner, and the electrons start seeing further away, and then become sensitive to higher-ordered multipoles, $\Theta_l$. Fortunately, because of a very nice computational trick due to Zaldarriaga and Seljak called line of sight integration, we only need to take into account a relatively small number of these (six is enough), and so the system of relevant equations is therefore still tractable. Note that before 1996 or so, people actually included thousands of variables, to trace multipoles for the full range. Needless to say, this was slow, and other approximations were required.

A second issue is the very large value of $\tau^{\prime}$ at early times, which multiplies $(3\Theta_1 + v_b)$ in the Boltzmann equations. The latter factor is very small early on, and the product of the two is therefore numerically extremely unstable. The result is that the standard Boltzmann equation set is completely unstable if one simply implements the full expressions at early times. The solution to this problem is to use a proper approximation for $(3\Theta_1 + v_b)$ at early times. See the appendix for a derivation.

The full system

Photon temperature multipoles: $$ \boxed{ \begin{align} \Theta^\prime_0 &= -\frac{ck}{\mathcal{H}} \Theta_1 - \Phi^\prime, \\ \Theta^\prime_1 &= \frac{ck}{3\mathcal{H}} \Theta_0 - \frac{2ck}{3\mathcal{H}}\Theta_2 + \frac{ck}{3\mathcal{H}}\Psi + \tau^\prime\left[\Theta_1 + \frac{1}{3}v_b\right], \\ \Theta^\prime_\ell &= \frac{\ell ck}{(2\ell+1)\mathcal{H}}\Theta_{\ell-1} - \frac{(\ell+1)ck}{(2\ell+1)\mathcal{H}} \Theta_{\ell+1} + \tau^\prime\left[\Theta_\ell - \frac{1}{10}\Pi \delta_{\ell,2}\right], \quad\quad 2 \le \ell \lt \ell_{\textrm{max}} \\ \Theta_{\ell}^\prime &= \frac{ck}{\mathcal{H}} \Theta_{\ell-1}-c\frac{\ell+1}{\mathcal{H}\eta(x)}\Theta_\ell+\tau^\prime\Theta_\ell, \quad\quad \ell = \ell_{\textrm{max}}\\ \end{align} } $$ Photon polarization multipoles: $$ \boxed{ \begin{align} \Theta^\prime_{P0} &= -\frac{ck}{\mathcal{H}}\Theta_{P1} + \tau^\prime\left[\Theta_{P0} - \frac{1}{2}\Pi \right]\\ \Theta_{P\ell}^\prime &= \frac{\ell ck}{(2\ell+1)\mathcal{H}} \Theta_{\ell-1}^P - \frac{(\ell+1)ck}{(2\ell+1)\mathcal{H}} \Theta_{\ell+1}^P + \tau^\prime\left[\Theta_\ell^P - \frac{1}{10}\Pi\delta_{\ell,2}\right],\quad\quad 1 \le \ell \lt \ell_{\textrm{max}} \\ \Theta_{P,\ell}^\prime &= \frac{ck}{\mathcal{H}} \Theta_{\ell-1}^P-c\frac{\ell+1}{\mathcal{H}\eta(x)}\Theta_\ell^P+\tau^\prime\Theta_\ell^P, \quad\quad \ell = \ell_{\textrm{max}}\\ \end{align} } $$ Neutrino multipoles: $$ \boxed{ \begin{align} \mathcal{N}^\prime_0 &= -\frac{ck}{\mathcal{H}} \mathcal{N}_1 - \Phi^\prime, \\ \mathcal{N}^\prime_1 &= \frac{ck}{3\mathcal{H}} \mathcal{N}_0 - \frac{2ck}{3\mathcal{H}}\mathcal{N}_2 + \frac{ck}{3\mathcal{H}}\Psi \\ \mathcal{N}^\prime_\ell &= \frac{\ell ck}{(2\ell+1)\mathcal{H}} \mathcal{N}_{\ell-1} - \frac{(\ell+1)ck}{(2\ell+1)\mathcal{H}}\mathcal{N}_{\ell+1},\quad\quad 2 \le \ell \lt \ell_{\textrm{max},\nu} \\ \mathcal{N}_{\ell}^\prime &= \frac{ck}{\mathcal{H}} \mathcal{N}_{\ell-1}-c\frac{\ell+1}{\mathcal{H}\eta(x)}\mathcal{N}_{\ell}, \quad\quad \ell = \ell_{\textrm{max},\nu}\\ \end{align} } $$
Cold dark matter and baryons: $$ \boxed{ \begin{align} \delta_{\rm CDM}^\prime &= \frac{ck}{\mathcal{H}} v_{\rm CDM} - 3\Phi^\prime \\ v_{\rm CDM}^\prime &= -v_{\rm CDM} -\frac{ck}{\mathcal{H}} \Psi \\ \delta_b^\prime &= \frac{ck}{\mathcal{H}}v_b -3\Phi^\prime \\ v_b^\prime &= -v_b - \frac{ck}{\mathcal{H}}\Psi + \tau^\prime R(3\Theta_1 + v_b) \\ \end{align} } $$ Metric perturbations: $$ \boxed{ \begin{align} \Phi^\prime &= \Psi - \frac{c^2k^2}{3\mathcal{H}^2} \Phi + \frac{H_0^2}{2\mathcal{H}^2} \left[\Omega_{\rm CDM 0} a^{-1} \delta_{\rm CDM} + \Omega_{b 0} a^{-1} \delta_b + 4\Omega_{\gamma 0} a^{-2}\Theta_0 + 4\Omega_{\nu 0}a^{-2}\mathcal{N}_0\right] \\ \Psi &= -\Phi - \frac{12H_0^2}{c^2k^2a^2}\left[\Omega_{\gamma 0}\Theta_2 + \Omega_{\nu 0}\mathcal{N}_2\right] \\ \end{align} } $$ In the equations above $\Pi = \Theta_2 + \Theta_0^P + \Theta_2^P$ and $R = \frac{4\Omega_{\gamma 0}}{3\Omega_{b 0} a}$ (note that our $R$ is $1/R$ in Dodelson). Note that only one of the potentials are dynamical - $\Psi$ follows from $\Phi$ so you don't have to solve this. If you are a master student you don't have to implement neutrinos and polarization. In that case just ignore the neutrino and polarization equations and put $\mathcal{N}_\ell = 0$ and $\Theta^P_\ell = 0$ in all the other equations above.

The tight coupling regime

Ideally we would just solve the original ODE system for all times, however this does not work as thew system is very stiff and numerically unstable when $\tau$ is very large. We therefore have to solve some slightly different equations in the tight coupling regime. In the tight coupling regime, the only differences are 1) that one should only include $\ell = 0$ and 1 for $\Theta_\ell$ (and none for polarization) - all higher moments are given by those, and 2) that the expressions for $\Theta'_1$ and $v_b'$ are quite a bit more involved (see the appendix for how to derive this): $$ \boxed{ \begin{align} q &= \frac{-[(1-R)\tau^\prime + (1+R)\tau^{\prime\prime}](3\Theta_1+v_b) - \frac{ck}{\mathcal{H}}\Psi + (1-\frac{\mathcal{H}^\prime}{\mathcal{H}})\frac{ck}{\mathcal{H}}(-\Theta_0 + 2\Theta_2) - \frac{ck}{\mathcal{H}}\Theta_0^\prime}{(1+R)\tau^\prime + \frac{\mathcal{H}^\prime}{\mathcal{H}} - 1}\\ v_b^\prime &= \frac{1}{1+R} \left[-v_b - \frac{ck}{\mathcal{H}}\Psi + R(q + \frac{ck}{\mathcal{H}}(-\Theta_0 + 2\Theta_2) - \frac{ck}{\mathcal{H}}\Psi)\right]\\ \Theta^\prime_1 &= \frac{1}{3} (q - v_b^\prime) \end{align} } $$ In the tight coupling regime, we get the same expressions for the higher-ordered photon moments as given by the initial conditions, $$ \boxed{ \begin{align} \Theta_2 &= \left\{ \begin{array}{l} -\frac{8ck}{15\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(with polarization)} \\ -\frac{20ck}{45\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(without polarization)} \end{array}\right. \\ \Theta_\ell &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau'} \Theta_{\ell-1}, \quad\quad \ell \gt 2\\ \Theta_0^P &= \frac{5}{4}\Theta_2 \\ \Theta_1^P &= -\frac{ck}{4\mathcal{H}\tau'}\Theta_2 \\ \Theta_2^P &= \frac{1}{4}\Theta_2 \\ \Theta_\ell^P &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau'} \Theta_{\ell-1}^P, \quad\quad \ell \gt 2 \end{align} } $$

Initial conditions

The initial conditions are given by: $$ \boxed{ \begin{align} \Psi &= -\frac{1}{\frac{3}{2} + \frac{2f_\nu}{5}}\\ \Phi &= -(1+\frac{2f_\nu}{5})\Psi \\ \delta_{\rm CDM} &= \delta_b = -\frac{3}{2} \Psi \\ v_{\rm CDM} &= v_b = -\frac{ck}{2\mathcal{H}} \Psi\\ &\text{Photons:}\\ \Theta_0 &= -\frac{1}{2} \Psi \\ \Theta_1 &= +\frac{ck}{6\mathcal{H}}\Psi \\ \Theta_2 &= \left\{ \begin{array}{l} -\frac{8ck}{15\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(with polarization)} \\ -\frac{20ck}{45\mathcal{H}\tau^\prime} \Theta_1, \quad\quad \textrm{(without polarization)} \end{array}\right. \\ \Theta_\ell &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau^\prime} \Theta_{\ell-1}\\ &\text{Photon Polarization:}\\ \Theta_0^P &= \frac{5}{4} \Theta_2 \\ \Theta_1^P &= -\frac{ck}{4\mathcal{H}\tau'} \Theta_2 \\ \Theta_2^P &= \frac{1}{4}\Theta_2 \\ \Theta_\ell^P &= -\frac{\ell}{2\ell+1} \frac{ck}{\mathcal{H}\tau^\prime} \Theta_{\ell-1}^P \\ &\text{Neutrinos:}\\ \mathcal{N}_0 &= -\frac{1}{2} \Psi \\ \mathcal{N}_1 &= +\frac{ck}{6\mathcal{H}}\Psi \\ \mathcal{N}_2 &= -\frac{c^2k^2 a^2 (\Phi+\Psi)}{12H_0^2\Omega_{\nu 0}}\\ \mathcal{N}_\ell &= \frac{ck}{(2\ell+1)\mathcal{H}} \mathcal{N}_{\ell-1}, \quad\quad \ell \ge 3 \end{align} } $$ where $f_{\nu} = \frac{\Omega_{\nu 0}}{\Omega_{\gamma 0} + \Omega_{\nu 0}}$. If you don't include neutrinos then set $f_\nu = 0$. Note that $\Psi$ is not a dynamical variable in the code and is only used here to set the rest of the variables. Since the equation system is linear we are free to choose the normalization of $\Psi$ as we want when we solve it (the normalization can be done in the end). The particular normalization we use here is such that when we in the next milestone are to compute power-spectra then $\Psi_{\rm true}^2 = \Psi_{\rm ours}^2 \mathcal{P}_\mathcal{R}$ where $\mathcal{P}_\mathcal{R}$ is the usual curvature perturbation power-spectum, the perturbations set up by inflation, $\mathcal{P}_\mathcal{R}(k) = A_s(k/k_{\rm pivot})^{n_s-1}$ with $A_s,n_s,k_{\rm pivot}$ (primoridal amplitude, spectral index and the pivot scale) are the same parameters as is standard in the litterature (and used in codes like CAMB and CLASS).

What you have to do

Implement a class/module (Perturbations.h if you use the C++ template) that takes in a BackgroundCosmology and RecombinationHistory object - the ones you created in the previous two milestones - and use this to evolve the perturbations of the Universe.

"All" you have to do is to set up the initial conditions, make the function that sets the right hand side of the coupled ODE system and then solve it and store the solution. However one thing that complicates it is that when integrating the perturbations you can't just integrate the full system directly - it is unstable due to the tight coupling between photons and baryons in the very early Universe. You therefore have to start off solving the tight coupling system: here we only have to include two photon multipoles $\Theta_0$ and $\Theta_1$ (and no polarization multipoles). Once tight coupling ends you have to switch to the full system. When doing this you will have to set the initial conditions from the tight coupling solution + you have to give a value to the multipoles we have in the full system, but don't include in the tight-coupling regime (the value of these are given by the same relations as we have in the initial conditions, e.g. $\Theta_2$ is given by the value of $\Theta_1$). This means you will need to make two functions to set initial conditions to your ODE vector - one at the start and one after tight coupling end - and you will have to implement two functions that sets the right hand side of the ODE system - one for tight coupling and one for the full system. Once you have this you will have to make a vector of $k$-values for which we will solve the system and integrate the perturbations from the start till today for all these $k$-values. You have to extract and store the results of the ODE and spline $f(x,k)$ for all of the quantities ($\delta_b, v_b, \Phi,\Theta_0$ etc.) that are required in the line-of-sight integrals.

The vector of the perturbations that we integral will have between $10$ and $30$ or so elements depening on how many multipoles we include (and this vector will be different in the two regimes). This means its very important to keep track over which index in the vector corresponds to which quantity. For example if we don't have polarization or neutrinos then one possible way to do this is to place (for the full system): $\delta_{\rm CDM}$ is $[0]$, $v_{\rm CDM}$ is $[1]$, $\delta_b$ is $[2]$, $v_b$ is $[3]$, $\Phi$, is $[4]$ and $\Theta_\ell$ is $[5+\ell]$ for $\ell = 0,1,\ldots,\ell_{\rm max}$ for a total of $6+\ell_{\rm max}$ components. And in the tight coupling regime we just have $7$. You can choose to do this "by hand", however its easy to make a mistake. A much better way, less prone to making index mistakes (plus it would make it much easier to add new components if you needed to do that), is to precompute indices for the different components (in both regimes) and have variables that tells us what is the index of each component. Thus to access $\Phi$ we write something like $y[$index_Phi$]$. In the code template I have addeed such a system with examples for how to use it, but do what you think is best (and understands!).

Testing your code

The most useful thing for you is probably going to be comparing to plots shown below, but we can also do some more quantitative checks by comparing to analytical solutions. We don't have exact analytical solutions to the full equation set - it's way to complicated for that - however we can derive some analytical approximation for some of the quantities in certain regimes that we can use to check that we are computing things correctly. Here are just some simple examples:

We can check that the temperature multipoles are integrated correctly by comparing it with the analytical approximation we derived in the lectures. In its very simplest form we have (before recombination) $\Theta_0 + \Psi \propto \cos(k\eta/\sqrt{3})$. This is not perfect, but you should find similar looking oscillations. A full analytical approximation is derived in Hu & Sugiyama (1995) which is good to $\sim 10-20\%$.

We can check that the gravitational potential is solved correctly by comparing this to analytical expections (see Dodelson or Baumann for expressions). For example one simple and concrete test: if we make a run with just matter and radiation then on super-horizon scales (e.g. $k \lesssim 0.001/$Mpc) the gravitational potential in the radiation dominated regime is related to the gravitational potential in the matter dominated regime via $\Phi_{\rm matter-era} = \frac{9}{10}\Phi_{\rm radiation-era}$. We also have $\Phi \approx \Phi_{\rm ini} \frac{\sin(y) -y \cos(y)}{y^3/3}$ with $y = k\eta/\sqrt{3}$ for small scale modes that enter the horizon in the radiation era ($k\gtrsim 0.1/Mpc$) and we expect $\Phi \approx $ constant in the matter era.

We can check that the dark matter perturbations are integrated correctly by comparing to analytical approximations (the Meszaros equation, see Dodelson or Baumann for equations). For example on subhorizon scales $k\gg \mathcal{H}$ we should have $\delta \propto a = e^x$ in the matter era and $\delta\propto \log(a)$ (so basically frozen) in the radiation era. Another useful analytical approximation is that the growth rate of matter perturbations satisfy $f \equiv \frac{1}{\delta}\frac{d\delta}{dx} \approx \Omega_M(x)^{0.55}$ in the matter era.




    Figure 1: Here are some plots to compare your results to. We show the evolution of the perturbation for a toy-cosmology with parameters $\Omega_{b 0} = 0.05$, $\Omega_{\rm CDM 0} = 0.45$, $\Omega_{\Lambda 0} = 0.5$, $\Omega_{\nu 0} = 0$, $Y_p = 0$ and $h=0.7$. Polarization and neutrions are not included. You can also run my code online here to get results to compare to (NB: cannot guarantee there are no bugs in that code). For the baryon overdensity and velocity we plot the absolute value above (as it can be negative).

    Appendix

    Appendix: Tight coupling time

    Tight coupling should be used when all the following conditions holds: 1) $\left|\frac{d\tau}{dx}\right| > 10$, 2) $\left|\frac{d\tau}{dx}\right| > 10\cdot \frac{ck}{\mathcal{H}}$ and 3) we are not past the onset of recombination which we here just take to be $x=-8.3$ ($z=4000$). You should make a routine that, for a given $k$, starts at the initial time and then searches for when tight coupling ends. This would be the time when $\left|\frac{d\tau}{dx}\right| = 10 \cdot \text{min}(1, \frac{ck}{\mathcal{H}})$ and if the time you find has $x > -8.3$ then set the end time of tight coupling to be $-8.3$. It will only be for the small scales modes $k \gtrsim 0.15 /$ Mpc that you will find that tight coupling ends earlier than this.

    Appendix: Tight coupling equations

    The aim here is to derive an (numerically stable) equation for $[3\Theta_1 + v_b]$ in the tight coupling regime. The equations for $\Theta_1$ and $v_b$ are given by $$3\Theta^\prime_1 = \frac{ck}{\mathcal{H}}\left( \Theta_0 - 2\Theta_2 + \Psi \right) + \tau^\prime\left[3\Theta_1 + v_b\right], \\ v_b^\prime = -v_b - \frac{ck}{\mathcal{H}}\Psi + \tau^\prime R[3\Theta_1 + v_b]$$ Summing these gives us $$[3\Theta_1 + v_b]^\prime = \frac{ck}{\mathcal{H}}( \Theta_0 - 2\Theta_2) + \tau^\prime(1+R)\left[3\Theta_1 + v_b\right] - v_b\tag{A}$$ Taking the derivative of this expression, using $R^\prime = - R$, gives us $$[3\Theta_1 + v_b]^{\prime\prime} = -\frac{ck}{\mathcal{H}}\frac{\mathcal{H}^\prime}{\mathcal{H}}( \Theta_0 - 2\Theta_2) + \frac{ck}{\mathcal{H}}( \Theta_0 - 2\Theta_2)^\prime + (\tau^{\prime\prime}(1+R) - R\tau^\prime)\left[3\Theta_1 + v_b\right] + (\tau^\prime(1+R) - 1)\left[3\Theta_1 + v_b\right]^\prime + 3\Theta_1^\prime$$ where we have used $-v_b^\prime = -[3\Theta_1 + v_b]^\prime + 3\Theta_1^\prime$ to get ridd of the $-v_b^\prime$ term. Substituting in the equation for $\Theta_1^\prime$ in this last term we arrive at $$[3\Theta_1 + v_b]^{\prime\prime} = (\tau^{\prime\prime}(1+R) + (1-R)\tau^\prime)\left[3\Theta_1 + v_b\right] + (\tau^\prime(1+R) - 1)\left[3\Theta_1 + v_b\right]^\prime + \frac{ck}{\mathcal{H}}\left( \Theta_0 - 2\Theta_2 + \Psi + \Theta_0' - 2\Theta_2' - \frac{\mathcal{H}^\prime}{\mathcal{H}}(\Theta_0 - 2\Theta_2)\right)$$ Now comes the approximation. In the tight coupling regime one can show that to a good approximation $(3\Theta_1 + v_b) \propto \frac{1}{\tau'} \propto \eta$ since $\tau' \propto \frac{1}{a}$ and $\eta \propto a$ in a radiation dominated Universe. This means that $\frac{d^2}{d\eta^2}(3\Theta_1 + v_b) \approx 0$ or in terms of $x$ $$[3\Theta_1 + v_b]^{\prime\prime} \approx -\frac{\mathcal{H}^\prime}{\mathcal{H}}[3\Theta_1 + v_b]^{\prime}$$ Using this approximation in the expression above and solving for $q\equiv \left[3\Theta_1 + v_b\right]^\prime \to \Theta_1' = \frac{q-v_b'}{3}$ we arrive at $$q = \frac{-(\tau^{\prime\prime}(1+R) + (1-R)\tau^\prime)\left[3\Theta_1 + v_b\right] - \frac{ck}{\mathcal{H}}\left(\Psi + \Theta_0' - 2\Theta_2' + (1- \frac{\mathcal{H}^\prime}{\mathcal{H}})(\Theta_0 - 2\Theta_2)\right)}{\tau^\prime(1+R) + \frac{\mathcal{H}^\prime}{\mathcal{H}} - 1}$$ Further we can show that $\Theta_2^\prime \sim 3\Theta_2 \ll \Theta_0$ so this term can be ignored. Finally using (A) to solve for $\tau^\prime(1+R)\left[3\Theta_1 + v_b\right]$ we get $$\tau^\prime(1+R)\left[3\Theta_1 + v_b\right] = \frac{q - \frac{ck}{\mathcal{H}}( \Theta_0 - 2\Theta_2) + v_b}{1+R}$$ which can be substitute into the equation for $v_b^\prime$ to give us $$v_b^\prime = \frac{1}{1+R} \left[-v_b - \frac{ck}{\mathcal{H}}\Psi + R(q + \frac{ck}{\mathcal{H}}(-\Theta_0 + 2\Theta_2) - \frac{ck}{\mathcal{H}}\Psi)\right]$$

    Appendix: Parallelizing the integration

    The integration is the thing that will take most of the time. If the integration is slow and you want to make it go faster here is one very 'easy' way to speed it up by parallelizing it. Most computers these days have more than one core that can do work so we can take advantage of that. The loop over $k$ where we integrate can be done in parallel and adding this to the code is fairly easy: add the compiler flag -fopenmp and before the loop add a OpenMP tag telling the compiler to do it parallel: 

    #pragma omp parallel for schedule(dynamic, 1)
for(int i = 0; i < n_k_values; i++){
  //... do some work
  //... do some work
  //... do some work
}

    Thats it. But... you will have to be very careful with what you do inside that loop or things can go very wrong. You have to make sure that the multiple threads you start don't write to the same place in memory - if that happens things will go wrong. If you want to try this now is the time to go and read an introduction to OpenMP so that you understand what is going on. And for debugging (does not work on a Mac) you can add the compiler flag -fsanitize=thread and it will check for errors. In any case: if you want to try this then only do this when you have made it work. That way you have something to compare to and you won't waste time you don't have.