# Numerically solving for cosmological scale factor in mathematica

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$$.

• $t_0$ doesn’t mean $t=0$. – G. Smith Dec 1 '19 at 23:47
• $t=0$ is the Big Bang; $t=t_0=13.8$ Gy is now. – G. Smith Dec 2 '19 at 0:04
• Integrate forward and backward from now. – G. Smith Dec 2 '19 at 0:42
• Why do you want to express the Hubble constant in GeV? I suggest expressing it in inverse gigayears. – G. Smith Dec 2 '19 at 3:06
• If all you did was multiply through by $a(t)^2$, why is $\Omega_\Lambda$ not multiplied by that? And how did the square of a first derivative turn into a second derivative? – G. Smith Dec 2 '19 at 16:55

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..

• Expressing $H$ in $\mathrm{GeV}$ is perfectly meaningful, it's called using natural units. – Javier Dec 3 '19 at 17:28
• Natural units would be $H$ in terms of the Planck mass, not GeVs ! – Alfred Dec 3 '19 at 18:12
• $H$ would be about 1.2 $10^{-61}$ (dimensionless) in "natural units" – Alfred Dec 3 '19 at 18:25

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$$.