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. 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 the Saha equation, the RHS of the Peebles equation and solve the ODE to get $X_e$. Solve the ODE for $\tau$ and compute the visibility function $g$ from this.
  • Plot of $X_e(x)$ from the Peebles equation compared with the results of the Saha equation. (15p).
  • Plot of $\tau(x)$, $-\tau^{\prime}(x)$ and $\tau^{\prime\prime}(x)$ (10p).
  • 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) (10p).
  • Compute the times ($x_{\rm decoupling}$, $z_{\rm decoupling}$, $t_{\rm decoupling}$) for the last scattering surface surface and for when recombination happens. You can use the definition $\tau = 1$ or as the time when visibility function peaks for when last scattering happens. You can use the definition $X_e = 0.1$ for when recombination happens. Compare these numbers with the prediction from the Saha equation. Compute the freeze-out abundance of free electrons today $X_e(x=0)$. If you include reionization compute these numbers with reionization turned off (15p).
  • Compute the sound-horizon at decoupling $r_s$ (5p).

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

Page: /

Theoretical background

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

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 the lecture notes above, in Callin (2006), Baumann 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 $$\boxed{\tau' = \frac{d\tau}{dx} = -\frac{c 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_{c0}}{m_H a^3}$$ where $\rho_{c0} \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 = (1-Y_p) \frac{\Omega_{b0} \rho_{c0}}{m_H a^3}$$ where $m_H$ is the hydrogen mass and $Y_p$ is the Helium fraction and is a new cosmological parameter (BBN gives us $Y_p \approx 0.24$ and observational constraints are consistent with this - our fiducial value in this project is $Y_p = 0.245$). 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 ignore Helium then the first one of these reads $$\boxed{\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_\gamma$, 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_\gamma = T_{\rm CMB 0} / 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, $$\boxed{\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 $$ \boxed{ \begin{align} C_r(T_b) &= \frac{\Lambda_{2s\rightarrow1s} + \Lambda_{\alpha}}{\Lambda_{2s\rightarrow1s} + \Lambda_{\alpha} + \beta^{(2)}(T_b)}, \,\,{\rm (dimensionless)},\\ H &, \,\,{\rm (dimension~1/s)}\\ \Lambda_{2s\rightarrow1s} &= 8.227 \textrm{s}^{-1}, \,\,{\rm (dimension~1/s)}\\ \Lambda_{\alpha} &= H\frac{(3\epsilon_0)^3}{(8\pi)^2 n_{1s}}, \,\,{\rm (dimension~1/s)}\\ n_{1s} &= (1-X_e)n_H, \,\,{\rm (dimension~1/m^3)}\\ n_H &= (1-Y_p)\frac{3H_0^2\Omega_{b0}}{8\pi G m_H a^3}, \,\,{\rm (dimension~1/m^3)}\\ \beta^{(2)}(T_b) &= \beta(T_b) e^{3\epsilon_0/4T_b}, \,\,{\rm (dimension~1/s)}\\ \beta(T_b) &= \alpha^{(2)}(T_b) \left(\frac{m_eT_b}{2\pi}\right)^{3/2} e^{-\epsilon_0/T_b}, \,\,{\rm (dimension~1/s)} \\ \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), \,\,{\rm (dimension~m^3/s)}\\ \phi_2(T_b) &= 0.448\ln(\epsilon_0/T_b), \,\,{\rm (dimensionless)}. \end{align} } $$ where $\alpha \simeq \frac{1}{137.0359992}$ in the expression for $\alpha^{(2)}$ is the fine-structure constant (it can also be expressed by the Thompson cross-section $\sigma_T = \frac{8\pi}{3}\frac{\alpha^2}{m_e^2}\,\,({\rm dimension~m^2})$ that is availiable in the constants we provide). 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, $$ \boxed{ \tilde{g}(x) = -\tau' e^{-\tau} } $$ which has the property that $\int_{-\infty}^{0} \tilde{g}(x)dx = 1$. This function is therefore a true probability distribution, and describes the probability density for a given photon to last time having scattered 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$. The first- and second-order derivatives of $\tilde{g}$ will be also required, and so these must also be properly splined, just as $\tau$ was.

The final thing you have to compute is the so-called sound-horizon at decoupling - the total distance a sound-wave in the photon-baryon plasma can propagate from the Big Bang and untill photons decouple. This quantity is not so relevant to theory we have gone through so far, but we will see later in this course that this length scale is imprinted in the CMB and the distribution of dark matter and galaxies in our Universe so it's important to know how big it is. The sound-speed of the coupled photon-baryon plasma is slightly lower than the sound-speed of photons ($c/\sqrt{3}$) and given by $c_s = c \sqrt{\frac{R}{3(1+R)}}$ where $R = \frac{4\Omega_{\gamma 0}}{3\Omega_{b 0} a}$ so the sound-horizon is $$\boxed{s(x) = \int_0^{a} \frac{c_s dt}{a} = \int_{-\infty}^{x} \frac{c_s dx}{\mathcal{H}} \to \frac{ds(x)}{dx} = \frac{c_s}{\mathcal{H}}\,\,\,\text{with}\,\,\,s(x_{\rm ini}) = \frac{c_s(x_{\rm ini})}{\mathcal{H}(x_{\rm ini})} }$$ This should be computed and splined (it will be useful to have in the next milestone) and evaluating it at $x = x_{\rm decoupling}$ gives us the sound-horizon at decoupling $r_s \equiv s(x_{\rm decoupling})$.

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. It's useful to know that it is possible to compute derivatives directly from a spline.

Note that in all the equations in the lectures we are using units where $\hbar = c = k_b = 1$. You have to do dimensional analysis to restore these constants in these equations (or atleast check that the dimensions are correct). 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$. You can solve it with any initial condition and subtract from $\tau$ the value you get at $x=0$. This however likely leads to truncation errors close to $x=0$ (manifesting as noise). A better solution is to solve it from $x=0$ and going back in time. From $\tau$ compute the visibility function $\tilde{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 overflow (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 inclusion 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.

To test that the visibiility function is correct remember that this is a PDF so if you integrate it over all $x$'s from the past till today then that should give $1$, i.e. $\int_{x_{\rm min}}^0 \tilde{g}(x)\,dx = 1$.

Compare your results to the plots below and those in Callin (2006). You can also run my code online here to get results to compare to (NB: cannot guarantee there are no bugs in that code).

    Figure 1: 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: Saha equation with helium

    Above we introduced the primordial mass fraction of $^4$ Helium, $Y_p$, 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+}} = 4\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+}} = \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^2e^2$ and the helium core has $q = 2e$ instead of $q=e$ for hydrogen so the energy is therefore $4$ times as large). 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 = f_e n_b$ 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 (close to) best-fit values from Planck 2018: $z_{\rm reion} = 8$, $\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 further 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 a dimensionless expression you automatically know you have to add $k_b$ in front.

    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 as the photon temperature. This is a very good approximation for recombination, but if you want to do it more properly then we can compute this. In the absence of interactions the temperature of a plasma of non-relativistic particles goes down as $T\propto 1/a^2$ (as non-relativistic kinetic energy is $\frac{p^2}{2m}$ compared to $p$ for photons and $p \propto 1/a$), but in the early Universe interactions with photons keep it having the same temperature as photons. We can compute the temperature from the Boltzmann equation (just as we did for photons as the 0th order solution), but it also follows (see also this plus Eq. A.32-33 on page 211 for exact expressions) from the first law of thermodynamics. The result is $$\frac{dT_b}{dx} + 2T_b = +2\frac{\mu}{m_e}R\frac{d\tau}{dx}(T - T_b)$$ where the last term is energy injection from interactions with the photons and $R = \frac{4\Omega_{\gamma 0}}{3a\Omega_{b 0}}$. When $\tau^\prime \gg 1$ the equation forces $T_b \approx T$. Above $\mu$ is the mean molecular weight of baryons including free electrons and all ions of H and He as function of time (but its fine to just take $\mu = 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}}$. You will find that the baryon temperature only starts to diverge from the photon temperature at $z\lesssim 300$ (when $X_e\sim 10^{-3}$) so it will not be important for the CMB anisotropies. When solving this equation numerically then it might be useful to introduce $y = T_b/T-1$, derive an ODE for $y$, and then solve this. You might also need to adjust the accuracy parameters of the ODE solver to get it to work ($\epsilon_{\rm rel} = \epsilon_{\rm abs} = 10^{-12}$ seems to work fine).

    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. This output of this code is $X_e$ and $T_b$ which can be directly compared to your results.

    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_{\gamma 0}}{3\Omega_{b 0} a}$. When $\tau = 1$ the photons decouple and stop noticing the baryons and when $\tau_b = 1$ (we call this time $x_{\rm drag}$) the baryons stop noticing the photons. This happens a bit later since there are so many more photons than baryons the baryons. The time between photon decoupling and baryon decoupling is called the baryon-drag era (baryons gets dragged along even though photons have decoupled - we call this the "drag epoch"). If you include reionization you should compute these with the no-reionization $\tau$ (just compute both at the same time).

    The sound-horizon: $$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)}}$ is the sound-speed of the photon-baryon fluid and $R = \frac{4\Omega_{\gamma 0}}{3\Omega_{b 0} a}$. Evaluating this when $x = x_{\rm decoupling}$ ($\tau(x_{\rm decoupling}) = 1$) we get the sound-horizon at the last scattering surface $r_s$ (the scale that is imprinted in the CMB). Likewise we can evaluate this in the baryon-drag epoch ($x_{\rm drag}$) to get the scale that is imprinted in the large scale structure of the Universe (these two numbers are very similar though).