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:

  • In the paper define the quantities to be computed and give a short description of the algorithms used
  • One plot for each of the physical quantities as function of $x$, for three different $k$'s within the interval of interest. For the radiation quantities, plot only the first two multipoles (i.e. the energy density $\delta_\gamma = 4\Theta_0$ and the velocity $v_\gamma = -3\Theta_1$). Choose values of $k$ such that each of the three main regimes are shown; 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).
  • A transcript of the module written for the evaluations

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. This 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: $$ \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: $$ \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: $$ \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: $$ \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: $$ \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

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): $$ \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, \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: $$ \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: Tight coupling time

    For any given $k$ you should switch from the tight coupling equations to the full system when $\left|\frac{d\tau}{dx}\right| < 10 \cdot \text{min}(1, \frac{ck}{\mathcal{H}})$ and no later than than the start of recombination ($z \sim 2000$).

    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.