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

  • 1
    $\begingroup$ $t_0$ doesn’t mean $t=0$. $\endgroup$
    – G. Smith
    Commented Dec 1, 2019 at 23:47
  • 1
    $\begingroup$ $t=0$ is the Big Bang; $t=t_0=13.8$ Gy is now. $\endgroup$
    – G. Smith
    Commented Dec 2, 2019 at 0:04
  • $\begingroup$ Integrate forward and backward from now. $\endgroup$
    – G. Smith
    Commented Dec 2, 2019 at 0:42
  • $\begingroup$ Why do you want to express the Hubble constant in GeV? I suggest expressing it in inverse gigayears. $\endgroup$
    – G. Smith
    Commented Dec 2, 2019 at 3:06
  • 2
    $\begingroup$ 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? $\endgroup$
    – G. Smith
    Commented Dec 2, 2019 at 16:55

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

  • $\begingroup$ Expressing $H$ in $\mathrm{GeV}$ is perfectly meaningful, it's called using natural units. $\endgroup$
    – Javier
    Commented Dec 3, 2019 at 17:28
  • $\begingroup$ Natural units would be $H$ in terms of the Planck mass, not GeVs ! $\endgroup$
    – Alfred
    Commented Dec 3, 2019 at 18:12
  • $\begingroup$ $H$ would be about 1.2 $10^{-61}$ (dimensionless) in "natural units" $\endgroup$
    – Alfred
    Commented Dec 3, 2019 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}]$$


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.