Beta Decay - Derivation of the Friedmann Equation

Starting with the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric in geometrised units (where $c=1$ ) :

$\mathrm{d}s^2 = -\mathrm{d}t^2 + a^2(t)\left[\frac{1}{1-kr^2}\mathrm{d}r^2 + r^2 \left( \mathrm{d}\theta^2 + \sin^2\theta \mathrm{d}\phi^2 \right) \right]$

Where $a(t)$ is the scale factor of the universe. When $a(t) = 1$ , $k$ is the Gaussian curvature.

We can calculate the $00$ component of the Einstein Tensor to be:

$\mathbf{G}_{00} = \frac{3(k + \dot{a}^2)}{a^2}$

This was calculated using the Python code here.

We also know that the $00$ component of the Energy-Momentum Tensor is the energy density:

$\mathbf{T}_{00} = \rho$

We can now substitute values into the Einstein Field Equation:

$\mathbf{G}_{\mu\nu} = 8\pi G \mathbf{T}_{\mu\nu}$

To get:

$\mathbf{G}_{00} = 8\pi G \mathbf{T}_{00}$

$\frac{3(k + \dot{a}^2)}{a^2} = 8\pi G\rho$

$\frac{k}{a^2}+\left(\frac{\dot{a}}{a}\right)^2 = \frac{8\pi G}{3}\rho$

$H^2=\frac{8\pi G}{3}\rho - \frac{k}{a^2}$

Where $H$ is the Hubble rate. In the present day, the Hubble rate is simply the Hubble constant and represents the rate at which the universe is expanding.

Converting from geometrised units we get:

$H^2 = \frac{8\pi G}{3}\rho-\frac{kc^2}{a^2}$

Which is known as the first Friedmann Equation.

To get the equation used in the simulation, first we need to expand $\rho$ , the energy density:

$\rho = \frac{\rho_{0,R}}{a^4} + \frac{\rho_{0,M}}{a^3} + \rho_{0,\Lambda}$

Where $\rho_{0,R}$ is the current radiation density, $\rho_{0,M}$ is the current mass density and $\rho_{0,\Lambda}$ is the current vacuum density.

If we divide both sides of the Friedmann Equation by the Hubble constant squared, $H_0^2$ , we get:

$\frac{H^2}{H_0^2} = \frac{8\pi G}{3H_0^2}\left(\frac{\rho_{0,R}}{a^4} + \frac{\rho_{0,M}}{a^3} + \rho_{0,\Lambda}\right)-\frac{kc^2}{H_0^2a^2}$

Which can be simplified by defining coefficients called "Density Parameters":

$\frac{H^2}{H_0^2} = \Omega_{0,R}a^{-4} + \Omega_{0,M}a^{-3} + \Omega_{0,\Lambda}+\Omega_{0,k}a^{-2}$

Where $\Omega_0 = \Omega_{0,R}+\Omega_{0,M}+\Omega_{0,\Lambda}$ and $\Omega_{0,k} = 1-\Omega_0$ .

We can then rearrange the equation in terms of $a$ and $\dot{a}$ (because $H = \frac{\dot{a}}{a}$ ):

$\dot{a} = H_0\sqrt{ \Omega_{0,R}a^{-2} + \Omega_{0,M}a^{-1} + \Omega_{0,\Lambda}a^2+\Omega_{0,k}a}$

Now, the Euler Method can be used to numerically solve the above equation:

$\begin{align} a(t+h) &\approx a(t) + h\times\dot{a}(t)\\ &\approx a(t) + h\times H_0\sqrt{ \frac{\Omega_{0,R}}{a^2} + \frac{\Omega_{0,M}}{a} + \Omega_{0,\Lambda}a^2+\Omega_{0,k}a} \end{align}$

Here, $h = \frac{24}{7\times10^4} \approx 0.00037142857$ .

This is the equation used in the simulation. Check it out!