Table of contents


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


Overview: Milestone II


The topic of the second milestone is the recombination history of the universe: how did the baryons go from being ionized (consisting only of free electron and protons) to being neutral (consisting of neutral hydrogen atoms)? The final goal of this part is to compute the optical depth as a function of $x$, the logarithm of the scale factor, $a$, and its derivatives, and the visibility function, $g$, and its derivatives.

The deliverables are the following:

  • In the paper define the quantities to be computed, give a short description of the algorithms used, and add one plot of $\tau(x)$, $-\tau^{\prime}(x)$ and $\tau^{\prime\prime}(x)$, one plot of $\tilde{g}(x)$, $\tilde{g}^{\prime}(x)$ and $\tilde{g}^{\prime\prime}(x)$ (appropriately scaled so that they fit into the same plot), and one plot of $X_e(x)$.
  • Compute the time ($x_*$ and redshift $z_*$) for the last scattering surface. You can use the definition $\tau = 1$ or as the time when visibility function peaks. Compute the time ($x_{\rm rec}$ and $z_{\rm rec}$) for when recombination is half-way, i.e. $X_e = 0.5$. Compare this with the prediction from the Saha equation. If you include reionization compute these numbers with reionization turned off.
  • A transcript of the module written for the evaluations

A short introduction to the milestone can be found in this PDF.


Theoretical background


The problem of this part of the project is to compute the optical depth, $\tau(x)$, and the so-called visibility function, $g(x)$, which is needed for integrating the Boltzmann-Einstein equations later on, and computing the CMB power spectrum. Physical descriptions and derivations will not be given here, but can be found in Callin (2006) or Dodelson. Here only the required definitions will be presented.

Light that travels through a medium can be absorbed by the medium. If we have a source emiting an intensity $I_0$ then an observer a distance $x$ away from the source will observe an intensity $I(x) = I_0 e^{-\tau(x)}$. The quantity $\tau$ is called the optical depth. If $\tau \ll 1$ then the medium does nothing (we say its optically thin) and if $\tau \gg 1$ then we will see nothing (the medium is optically thick). The transition between these two regimes is when $\tau \sim 1$. In cosmology the main 'absorption' is Thompson scattering of photons of free electrons. The optical depth is defined as $$\tau(\eta) = \int_{\eta}^{\eta_0} n_e \sigma_T a d\eta'$$ and quantifies the probability for scattering a photon between some previous time and today (the number of scatterings of photons by electrons per unit time is $cn_e\sigma_T$). The components involved in this expression are: $n_e=n_e(\eta)$; the electron density (ie., the number of free electrons per cubic meter) at time $\eta$, $\sigma_T = \frac{8\pi}{3}\frac{\alpha^2\hbar^2}{m_e^2c^2}$; the Thompson cross-section (Constants.sigma_T) and the scale factor $a$. The transition between the Universe being optically thick and optically thin happens around recombination when most of the free electrons are captured by free protons to form neutral hydrogen.

The expression for $\tau$ may also be written on an differential form, such that $$\tau' = \frac{d\tau}{dx} = -\frac{n_e \sigma_T }{H}.$$ This implies that we can use existing routines for solving differential equations to compute $\tau$, if we can only somehow compute $n_e$ at any time. And that's the difficult part.

Instead of actually computing $n_e$, we rather focus on the fractional electron density, $X_e \equiv n_e / n_H$, where $n_H$ is the proton density. If we will assume that all baryons are protons (ie., no helium or heavier elements), then $$n_H = n_b \approx \frac{\rho_b}{m_H} = \frac{\Omega_{b0} \rho_c}{m_H a^3}$$ where $\rho_c \equiv \frac{3H_0^2}{8\pi G}$ is the critical density of the universe today. If we do include Helium then only a (mass) fraction $1-Y_p$ of the baryons are hydrogen so $$n_H = (1-Y_p)n_b$$ where $Y_p$ is the Helium fraction and is a new cosmological parameter (BBN gives us $Y_p \approx 0.24$ and obsrvational constraints are consistent with this - our fiducial value in this project is $Y_p = 0.24$). If you are a master student then you can take $Y_p = 0$ and ignore Helium alltogether in this project. If you are a PhD student you must include this. Note that this is the first time we actually use cosmological parameters for computing observables -- from now on, we will at least be sensitive to $\Omega_{b 0}$, $H_0$ and $Y_p$.

Now, we have two different equations available for $X_e$ as a function of temperature and density, namely the so-called Saha and Peebles' equations. If we ignor Helium then the first one of these reads $$\frac{X_e^2}{1-X_e} = \frac{1}{n_b} \left(\frac{m_e T_b}{2\pi}\right)^{3/2} e^{-\epsilon_0/T_b},$$ where $T_b$ is the baryon temperature of the universe, and $\epsilon = 13.6$ eV is the ionization energy of hydrogen (ie., the energy a photon needs in order to rip an electron away from a proton). In principle, one would have to solve separately for both $T_b$ and the photon temperature, $T_r$, but in practice it is an excellent approximation to set these equal (see the appendix for how the baryon temperature evolves). We can therefore assume that $T_b = T_r = T_{\rm CMB 0} / a = 2.725 \textrm{K} / a$. With the above information, Saha's equation reduces to a standard second-order equation in $X_e$, and can be solved directly using the normal formula, $y = (-b \pm \sqrt{b^2 - 4ac})/2a$. If we include helium ($Y_p \not = 0$) then it becomes a bit more complicated, see the appendix for more details.

Saha's equation is an excellent approximation when $X_e \approx 1$. When $X_e$ is noticeably smaller than one, better approximations are required, and one such approximation is the Peebles' equation, $$\frac{dX_e}{dx} = \frac{C_r(T_b)}{H} \left[\beta(T_b)(1-X_e) - n_H \alpha^{(2)}(T_b)X_e^2\right],$$ where $$ \begin{align} C_r(T_b) &= \frac{\Lambda_{2s\rightarrow1s} + \Lambda_{\alpha}}{\Lambda_{2s\rightarrow1s} + \Lambda_{\alpha} + \beta^{(2)}(T_b)}, \\ \Lambda_{2s\rightarrow1s} &= 8.227 \textrm{s}^{-1}\\ \Lambda_{\alpha} &= H\frac{(3\epsilon_0)^3}{(8\pi)^2 n_{1s}}\\ n_{1s} &= (1-X_e)n_H \\ \beta^{(2)}(T_b) &= \beta(T_b) e^{3\epsilon_0/4T_b} \\ \beta(T_b) &= \alpha^{(2)}(T_b) \left(\frac{m_e T_b}{2\pi}\right)^{3/2} e^{-\epsilon_0/T_b} \\ \alpha^{(2)}(T_b) &= \frac{64\pi}{\sqrt{27\pi}} \frac{\alpha^2}{m_e^2}\sqrt{\frac{\epsilon_0}{T_b}}\phi_2(T_b) \\ \phi_2(T_b) &= 0.448\ln(\epsilon_0/T_b). \end{align} $$ This looks a bit scary, but it's not too bad, really. First of all, the various constants are simply physical constants describing simple atomic physics; Peebles' equation takes into account the transition rates between the ground-state ($1s$) and the first excited state ($2s$) of the hydrogen atom. For even higher-accurate work, many more states should be included, and also other atoms, most notably the helium atom. But here we will be satisfied with the Peebles' equation. Second, the Peebles' equation is simply yet another linear first-order differential equation -- which by now is something we know how to solve. The fact that the right-hand side is slightly more complicated than previous examples doesn't really make things that much harder.

Next, we must decide when to switch from Saha's equation to Peebles' equation. For simplicity, we simply say that when $X_e \gt 0.99$, we use Saha; when $X_e \lt 0.99$, we use Peebles.

When we have computed $n_e$, the next step is to spline it, such that you can evaluate it at arbitrary values of $x$. One thing: Because $n_e$ varies over many, many orders of magnitude, it is useful to spline $\log n_e$ rather than $n_e$ itself; this function varies much more slowly with $x$, and is therefore easier to interpolate. But do remember to exponentiate the resulting splined function value, after interpolating to the desired $x$.

We are now ready to compute the optical depth, by solving the ODE for $\tau$ mentioned above. For initial conditions, remember that the optical depth at our place in the universe today is precisely zero.

The final item to compute is the so-called visibility function, $$ \tilde{g}(x) = -\tau' e^{-\tau},\,\,\,\,\,\,\int_{-\infty}^{0} \tilde{g}(x)dx = 1. $$ This function is actually a true probability distribution, and describes the probability density for a given photon to scatter at time $x$. As you will see, this function is sharply peaked around $x=-7$, corresponding to redshifts of $z\sim 1100$. The fact that this is so sharply peaked is the reason why we often refer to the recombination period as the surface of last scattering: This process happened during a very thin shell around us, at a redshift of $z \sim 1100$.

Finally, first- and second-order derivatives of $\tilde{g}$ will be also required, and so these must also be properly splined, just as $\tau$ was.


What you have to do


Implement a class/module (RecombinationHistory.h if you use the C++ template) that takes in recombination parameters (the helium abundance $Y_p$ - if you are a master student you won't need this so just put it to $0$ or remove it completely from the code) and a BackgroundCosmology object - the one you created in the previous milestone - and use this to solve the recombination history of the Universe.

I have added a few of the splines you will create to RecombinationHistory.h See this header for the names. If you need more splines then you can add them there.

Note that in all the equations we are using units where $\hbar = c = k_b = 1$. You have to do dimensional analysis to restore these constants in these equations. For example in $e^{-\epsilon_0/T_b}$ the expression in the exponential has to be dimensionless ($e^x= 1 + x + x^2 \ldots$ so taking the exponential of something with units gives nonsense as $x$ has different units than $x^2$ and so on). We have that $\epsilon_0 \sim 13.6056$ eV is an energy. To get this in Joules we can use that $1$ eV $ = 1.60217653\cdot 10^{-19}$ J. The only constant that is related to temperature is $k_b$ which has units of Joules per Kelvin. Thus $k_b T_b$ when $T_b$ is given in Kelvin is an energy so the correct expression is $e^{-\epsilon_0 / (k_b T_b)}$ with $\epsilon_0$ expressed in Joules and $T_b$ in Kelvin. All the physical constants you will need are placed in the Constants struct in Utils.h (for example write Constants.k_b to get the Boltzmann constant). Note that all of these constants are in SI units. See the appendix for more info about this.

Start by solving for the electron density $n_e$ ($X_e$) by solving the Saha and Peebles equations. The former is easy: it's just a quadratic equation for $X_e$ so it's easy to solve (but be careful: the analytical solution at early times will be on the form "huge" - "huge" which should give approximately $1$ so you might need to use the approximation $\sqrt{1+x} \approx 1 + \frac{x}{2}$ when $|x|\ll 1$). The Saha equation is only valid as long as $X_e \approx 1$ so you will have to switch to the Peebles equation when $X_e \lt 0.99$. Spline the result for $X_e$ and $n_e$ as this will be needed when computing $\tau$.

Once you have computed $n_e$ you can integrate the ODE for $\tau$ (again remember units). Remember to fix the integration constant by ensuring $\tau(x=0) = 0$ (e.g. you can solve it with any initial condition and subtract from $\tau$ the value you get at $x=0$). From $\tau$ compute the visibility function $g$. Spline both of these functions.

One thing to be careful about: at late times $T_b$ becomes small so $\epsilon_0/T_b$ becomes huge and $e^{\epsilon_0/T_b}$ becomes enormous and can overlow (the maximum double is $\sim 10^{300}$). For $e^{-\epsilon_0/T_b}$ this is not an issue as this will underlow to $0$ which is fine, but if you compute it as the equations are given (i.e. first compute $\beta$ and then use this to compute $\beta_2$) then you must be careful to avoid overflow for $\beta_2$. Either write out the full expression for $\beta_2$ and combine the exponentials or check if $\epsilon_0/T_b$ is larger than say $200$ or so and if it is then simply put $\beta_2 = 0$ instead of computing the exponential.

If you are a PhD student you will have to include helium and reionization. The only difference for the latter is that you must correct the result for $X_e$ that you get from the Peebles equation before solving for $\tau$. It is important to have enough points in your $x$-arrays used to store $\tau$ and $g$ before splining them as to be able to capture the reionization process ($X_e$ goes from $\sim 0$ up to $\sim 1$ in a very narrow $z$-range).


Testing your code


In this module you will have to work with the root of all evil - physical units. The only nice thing here is that if you get any of the equations wrong then the result will be completely wrong and the code will likely crash (its hard to do a mistake and get a small error which is good). But make sure you tripple check your equations before implementing them.

If you include helium then I suggest you implement it without helium first (so just the one Saha equation) before including helium (the three Saha equations) and keep them as seperate methods. This gives you something to compare to and you should check that you get the same result in your helium solver when $Y_p = 0$ as your hydrogen-only solver. For the Peebles equation the includsion of Helium is just a change to $n_H$.

If you include reionization be sure to plot the results for $\tau$ and $g$ (and all the derivatives we need) in the critical regions around $z_{\rm recombination}$ and $z_{\rm reionization}$. Make sure that these curves are nice and smooth here. If they are not then you are likely not including enough points in your arrays to capture the reionization process.

Compare your results to the plots below and those in Callin (2006)





Here are some plots to compare your results to. We show the free electron fraction $X_e(x)$, the optical debth $\tau$ and the visibility function $\tilde{g}$ and its derivatives with and without including reionization 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$. For reionization we used $z_{\rm reion} = 11$ and $\Delta z_{\rm reion} = 0.5$. For the $X_e$ plot with Helium recombination we have $Y_p = 0.24$ and $z_{\rm reion He} = 3.5$ and $\Delta z_{\rm reion He} = 0.5$.


Appendix



Appendix: Saha equation with helium


Above we introduced the primordial mass fraction of $^4$ Helium, Yp, and noted that this is defined as the ratio between the total mass of helium and the total baryon mass. If we define a baryon as a proton or a neutron with mass $m_b = m_H$ - the mass a of hydrogen atom then since each helium atom has $4$ baryons we have $$Y_p = \frac{m_{\rm He}n_{\rm He}}{m_bn_b} = \frac{4n_{\rm He}}{n_b}$$ This gives us $\frac{n_H}{n_b} = \frac{n_b - 4n_{\rm He}}{n_b} = 1 - Y_p$ which is the relation we quoted in the main text. We can compute the recombination history of helium in the same way of for hydrogen using the Saha equation. Defining $$x_{\rm He+} = \frac{n_{\rm He+}}{n_{\rm He}},\,\,\, x_{\rm He++} = \frac{n_{\rm He++}}{n_{\rm He}},\,\,\, x_{\rm H+} = \frac{n_{\rm H+}}{n_{\rm H}},\,\,\,$$ we can derive three Saha equations $$n_e\frac{x_{\rm He+}}{1-x_{\rm He+}-x_{\rm He++}} = 2\left(\frac{m_eT_b}{2\pi}\right)^{3/2}e^{-\chi_0/T_b}$$ $$n_e\frac{x_{\rm He++}}{x_{\rm He+}} = 2\left(\frac{m_eT_b}{2\pi}\right)^{3/2}e^{-\chi_1/T_b}$$ $$n_e\frac{x_{\rm H+}}{1-x_{\rm H+}} = 2\left(\frac{m_eT_b}{2\pi}\right)^{3/2}e^{-\epsilon_0/T_b}$$ where $m_e$ is the electron mass, $\chi_0 = 24.5874$eV is the ionization energy of neutral helium and $\chi_1 = 4\epsilon_0 = 54.42279$ eV the ionization energy of $He^+$ (this latter is given by the same expression as for hydrogen - the one one that follows from Bohrs model $E\propto Q_e^2 Q_p^2$ - but with twice the change in the core). The free electron number density is $n_e = 2n_{\rm He++} + n_{\rm He+} + n_{\rm H+}$ or $$\frac{n_e}{n_b} = (2x_{\rm He++} + x_{\rm He+})\frac{Y_p}{4} + x_{\rm H+}(1-Y_p) \equiv f_e \implies X_e = \frac{f_e}{1-Y_p}$$ This is a bit more complicated to solve than for just hydrogen, but not much. One simple algorithm is to do it iteratively. We pick a value for $f_e \sim 1$ which gives us $n_e$ and allows us to solve algebraically for $x_{\rm He++},x_{\rm He+},x_{\rm H+}$ using the three Saha equations. Then in the end we compute what value for $f_e$ this corresponds to from the definition of $f_e$ above and use this as the starting value for the next iteration. The solution quickly converges in a few steps (use a stopping criterion like for example $|f_e - f_{e,old}| \lt 10^{-10}$).

Note that the Peebles equation with helium is just as we quoted it in the main text, but with using $n_H = (1-Y_p)n_b$ (in the equation and in the expression for $n_{1s}$) instead of simply $n_H = n_b$ as we did for hydrogen only.


Appendix: Reionization


Around redshift of $z\sim 10$ stars are made in large numbers and high energy radiation from these stars hits hydrogen and helium atoms and kicks loose the electrons. This is a complicated process and we will not get into details about it here. What we are mainly interested in is how this affects the propagation of CMB photons and how that again affects the CMB spectrum. For this purpose we can put up a very simple model that turns out to work pretty very well: what happens is that free electron fraction $X_e$ quickly rises from its Peebles's equation prediction ($X_e \sim 10^{-4}$ at late times) - the one we solved for above - to a value close to $1$. A simple model that described this is (the same as used in CAMB): $$X_e = X_e^{\rm Peebles} + \frac{1+f_{\rm He}}{2}\left(1+\tanh\frac{y_{\rm reion}-y}{\Delta y_{\rm reion}}\right)$$ where $y = (1+z)^{3/2} = e^{-3x/2}$, $f_{\rm He} = \frac{Y_p}{4(1-Y_p)}$, $y_{\rm reion} = (1+z_{\rm reion})^{3/2}$, $\Delta y_{\rm reion} = \frac{3}{2}\sqrt{1+z_{\rm reion}}\Delta z_{\rm reion}$ (to be complete we should add another parameter telling us what fraction of the atoms gets reionized, but we will assume that this value is $1$). This then has two free parameters: $z_{\rm reion}$ and $\Delta z_{\rm reion}$ that depends on when and how fast this process happens that we don't know a priori but can be fitted to the CMB data (as we will see later on one clear signal of this is a bump in the polarization $E$ mode spectrum at low $\ell$). In addition to this Helium probably gets doubly reionized at late reshifts which we can model by adding a term (again just using CAMB's model) $$\frac{f_{\rm He}}{2}\left(1 + \tanh \frac{z_{\rm He reion}-z}{\Delta z_{\rm He reion}}\right)$$ to the equation above. The fiducial values we'll use for the free parameters here are $z_{\rm reion} = 11$, $\Delta z_{\rm reion} = 0.5$, $z_{\rm He reion} = 3.5$ and $\Delta z_{\rm He reion} = 0.5$.

See the plots below for how $X_e$ and $\tau$ are affected by reionization. Notice in particular that we get a plateau in the graph for $\tau$ and the value of $\tau$ here, $\tau(z = z_{\rm reion})$, is what is most commonly used to describe reionization (i.e. instead of talking about when reionization happened, $z_{\rm reion}$, we instead talk about the value of the optical debth at this plateau).

If you include reionization its useful (but not required) to compute $\tau$ with and without reionization as $\tau$ as we compute it defined from our perspective ($\tau = 0$ today), but the early Universe don't know about reionization. If we for example want to figure out when the Universe becomes opaque (define this to be $\tau = 1$) then we should use $\tau_{\rm noreion}$ for this. Computing this can be done in the same ODE as that for $\tau$, just add another equation to the system.


Appendix: Dimensional analysis


In this course we often use natural units for which $c = \hbar = k_b = 1$. When we need to solve our equations numerically we will sometimes need to go back to normal units. Here is how you can do this.

Note that $[c] = m / s$, $[\hbar] = Js = kg\cdot m^2/s$ and $[k_b] = J / K$. If you have an equation like $$\frac{d\tau}{dx} = -\frac{n_e \sigma_T}{H}$$ and you want to introduce the missing constants then since $\tau$ and $x$ is dimensionless the right hand side must be dimensionless. $[\sigma_T] = m^2$, $[n_e] = 1/m^3$ and $[H] = 1/s$ so the right hand side currently has units of $m^{-3} \cdot m^2 \cdot s = s/m$. You can do this mechanically by putting up the equation $$1 = [\hbar]^N [c]^M [k_b]^P \cdot [n_e \sigma_T / H] = (kg m^2 s^{-1})^N (m s^{-1})^M (kg m^2 s^{-2} K^{-1})^P \cdot \frac{s}{m}$$ And then solving the linear system that arises when requiring each of the units to match up. In this example $$0 = 2N + M + 2P-1 \,\,\,(\text{the power of 'm' on the RHS})$$ $$0 = -N - M - 2P+1 \,\,\,(\text{the power of 's' on the RHS})$$ $$0 = N + P \,\,\,(\text{the power of 'kg' on the RHS})$$ $$0 = -P \,\,\,(\text{the power of 'K' on the RHS})$$ The solution is $N=0$, $M=1$ and $P=0$ so the equation with the constants reintroduced is simply $\tau' = -\frac{cn_e\sigma_T}{H}$. This is very tedious, but it will always work if you do it correctly. However in this simple example its quite easy to see directly that we need to multiply it with $c$ as this automatically cancels the $s/m$ unit we already have.

There are some simplifications one can do. For example temperature only enters through the Boltzmann constant $k_b$ (Joules per Kelvin) so whenever you have a $T$ in an expression you know you have to add $k_b$ in front, i.e. make the replacement $T\to k_bT$.

Note that we cannot apply a (non-linear) function to something with units so whenever you have something like $e^{-\frac{\epsilon}{T}}$ you know that the argument of the exponential must be dimensionless so $[\epsilon/T] = 1$ and since $\epsilon$ is an energy (and so is $k_bT)$ the expression with units introduced is simply $e^{-\frac{\epsilon}{k_bT}}$.


Appendix: More things one can compute...


These things are optional, but here are some extra things one can compute.

Baryon temperature: We have made the approximation that the baryon temperature is always the same. This is a decent approximation (given all the other approximations we use), but if you want to do it more properly then we can compute this. The first law of thermodynamics gives us $$\frac{dT_b}{dx} = -2T_b - 2\frac{\mu}{m_e}\frac{d\tau}{dx}(T - T_b)$$ where the last term is energy injection from interactions with the photons. When $\tau^\prime \gg 1$ the equation forces $T_b \approx T$. Above $\mu \approx (1-Y_p)m_H$ is the mean molecular weight including free electrons and all ions of H and He as function of time (but its good enough to just take $\mu = (1-Y_p)m_H$). This equation has to be solved together with the Peebles equation (which depends on $T_b$) with the initial condition $T_b(x_{\rm ini}) = T(x_{\rm ini}) = T_{\rm CMB 0} e^{-x_{\rm ini}}$.

Improving recombination accuracy: One way to improve the accuracy of the recombination calculation is to include a fudge factor in the Peebles equation (this is what Recfast does). The value of this has been found by comparing to more accurate calculations (which take way too long to be used in a Boltzmann code). This amounts to replacing $\alpha_2$ with $f\alpha_2$ where $f$ is the so-called fudge factor. The best value for this is around $1.14$ for Recfast (and probably a bit higher for us). Including this and the point above you should be able to match the Recfast results for $X_e$ and $T_b$ very well. There will still be some discrepancies though around the time when Helium becomes neutral since our treatment of Helium is done with the approximative Saha equations. If you want to compare with Recfast you can download and run it from this page.

Baryon optical depth: In addition to the optical depth for the photons we can likewise compute an optical depth for the baryons: $$\frac{d\tau_b}{dx} = \frac{cR\sigma_Tn_e}{H} = R\frac{d\tau}{dx}$$ where $R = \frac{4\Omega_{r0}}{3\Omega_{b 0} a}$. When $\tau = 1$ the photons decouple and stop noticing the baryons and when $\tau_b = 1$ the baryons stop noticing the photons. If you include reionization you should compute these with the no-reionization $\tau$ (just compute both at the same time).

The sound horizon: $$r_s(x) = \int_0^a \frac{c_s dt}{a} = \int_{-\infty}^x \frac{c_s dx}{\mathcal{H}}$$ where $c_s = c \sqrt{\frac{R}{3(1+R)}}$. Evaluating this when $x = x_*$ ($\tau(x_*) = 1$) we get the sound horizon at the last scattering surface (the scale that is imprinted in the CMB and in the large scale structure of the Universe).