Numerically solving for cosmological scale factor in mathematica

Question

I would like to solve the scale factor $a(t)$ in a FRW metric.

The Friedmann equation for the scale factor is:

$H^2 \equiv \left( \frac{\dot{a}}{a} \right)^2 = H_0^2\left[\Omega_R \left(\frac{a(t)}{a(t_0)}\right)^{-4}+\Omega_M \left(\frac{a(t)}{a(t_0)}\right)^{-3} +\Omega_\Lambda \right]\\ $

$ \frac{d^2a(t)}{dt^2} = H_0^2\left[ \Omega_R \frac{a(t_0)^4}{a(t)^2} + {\Omega_M} \frac{a(t_0)^3}{a(t)} +\Omega_\Lambda \right]$

I believe $a(t_0\equiv t=0)=1$ which represents the Universe today. $H_0 \simeq 1.44\times 10^{-42} \ { \text{GeV}} $ is the expansion rate measured today, $\Omega_R\simeq 9.2\times 10^{-5}$ for radiation, $\Omega_M \simeq 0.31$ for matter, and $\Omega_\Lambda\simeq 0.69$ for dark energy. In the above equation I take $a(t_0) = 1$ for the reason I outline above.

I am trying to solve this numerically in mathematica but I am unsure what to do for the boundary conditions as I am getting an error when I solve the above double differential equation. I take $a(0)=1$ i.e. the Universe today but what do I take for the other boundary condition. I am solving for earlier times so is this a negative $t$? i.e. is the other boundary condition $a(-100000)=0$.

I suggest this as at the very start of the Universe $a(t_{\text{early}})= 0$.

Answer 1

I think you are making things complicated by writing the Friedmann equation in terms of of$\frac{a(t)}{a(t_0)}$.

The general form is to write it in this way

$$\frac{H^2} {H_0^2} = \Omega_ra(t)^{-4}+\Omega_ma(t)^{-3}+\Omega_{\Lambda}$$

So this becomes

$$H^2= H_0^2[\Omega_ra(t)^{-4}+\Omega_ma(t)^{-3}+\Omega_{\Lambda}]$$

$$\frac{\dot{a}(t)^2} {a(t)^2} =H_0^2[\Omega_ra(t)^{-4}+\Omega_ma(t)^{-3}+\Omega_{\Lambda}]$$

$$\dot{a}(t)^2=H_0^2[\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2]$$

let's take the square root

$$\dot{a}(t)=H_0 \sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}$$

From now on we can write

$$\frac{da}{dt}=H_0\sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}$$

Or

$$\frac{da}{H_0\sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}} = dt$$

Let us call $$F(a) = \frac{1}{H_0\sqrt{\Omega_ra(t)^{-1}+\Omega_ma(t)^{-2}+\Omega_{\Lambda}a(t)^2}}$$

So you have $$F(a)da = dt$$ which can be solved numerically.

From here you can easily get $t(a)$ (time as a function of the scale factor). But since we want $a(t)$ you just have to plot the data points as $x=t$ and $y=a$.

For the boundry conditions as many other poeple pointed out $a(t)$ should start from $0$. While doing the computation $0$ will cause an error since in the $F(a)$ we have $1/a^2$ etc. However you can avoid this by simply starting from $a = 0.0000...1$. For the upper limit you can put any number (non-negative) you want.

$H_0 \simeq 1.44\times 10^{-42} \ { \text{GeV}} $ is the expansion rate measured today,

You should use $H_0$ in terms of its normal units, $[H] = kms^{-1}Mpc^{-1}$. Or simply $H_0 \approx 70 km/s/Mpc$.

Answer 2

G. Smith is right on several points.

One is that when you multiplied by $a(t)^2$ you forgot the third term in the right-hand-side of your second equation. It should be $a(t)^2 \Omega_\Lambda$

But the worst is the left-hand-side.

It should be $(da/dt)^2$, not $d^2a/dt^2$.

It is a first order equation (albeit of the second degree), not a second order equation. So you need only one boundary condition, namely $a(t)=1$ now, that is $a(t_0)=1$

You solve it by taking the positive square root of the right-hand side, and this is the appropriate value of $da/dt$.

How you solve this numerically, starting from $a(t_0)=1$ as single initial condition for a first order equation is your problem.

In fact, I did it recently, using the most primitive explicit integration scheme, forwards in time, but only for two or three billion years, to check a value I had found in the literature, and it checked perfectly. I forgot on which wiki page that value was, I am sorry to say. (Of course I completely neglected $\Omega_R$)

As for the units, $H_0$ has the dimension of the inverse of a time (velocity divided by distance). Therefore $H_0$ multiplied by the Planck constant would indeed be an energy. I did not check it, but a reasonable value for $H_0$ might be $1.44\ 10^{-42}$ GeV/$\hbar$, but you need to divide by ℏ to get a meaningful result.

Why on earth would anyone use a GeV/$\hbar$ as a unit for "inverse time" is hard to imagine. But at least it is meaningful. Just GeV as unit of "inverse time" is meaningless..