# Problems with Computing the Friedmann Equation

I am trying to write a program which uses the Friedmann equation to graph the expansion of a universe given the starting parameters $K$ (geometry of spacetime), $R$ (curvature of spacetime) and $\rho$ (average mass density). So far, I have rearranged the equation to the form:

$$\frac{dt}{da} = \frac{1}{\sqrt{\frac{8}{3}\pi G \rho a^2 - K\frac{c^2}{R^2}}}$$

I then input this equation into Wolfram Alpha and integrated it, yielding the result:

$$t = \frac{\sqrt{\frac{3}{2\pi}}\log{\left(\sqrt{2\pi\rho G}\sqrt{8\pi a^2 \rho G - \frac{3c^2K}{R^2}}+4\pi a \rho G\right)}}{2\sqrt{\rho G}} + C$$

The program then plots values of $t$ for values of $a$ ranging from $0.5$ to $1.5$, inclusive.

The problem I have is that, when I run the program, it doesn't work. When $K = 1$, all of the values come out as NaN. When $K = -1$, all of the values of $t$ came out the same and, when $K = 0$, when given the values of our universe, it doesn't look as expected.

I have a feeling that there is a problem in my rearrangement and subsequent integration of the Friedmann equation, but I don't know enough calculus to know where I have gone wrong.

Can someone guide me on what I am doing wrong?

I have been testing this with $R = 1$ and $\rho = 9 \times 10^{-27} kgm^{-3}$.

## Edit:

$$t = \frac{\sqrt{\frac{3}{2\pi}}\log{\left(\sqrt{2\pi\rho G}\sqrt{|8\pi a^2 \rho G - \frac{3c^2K}{R^2}|}+4\pi a \rho G\right)}}{2\sqrt{\rho G}} + C$$

At $K = 1$, the values no longer come out as NaN but they are now like $K = -1$ where all of the values of $t$ are the same.

• Looking at your solution for $t$ reminds me why it's good to use units where $8 \pi G = c = 1$. – Michael Seifert Aug 11 '17 at 13:36
• On a more serious note: what values are you using for $\rho$ and $R$? This may just be an issue of inconsistent units. I also note that your solution is only valid when $\rho$ is constant, which is true for vacuum energy but won't get you the matter-dominated or radiation-dominated cosmological histories of our Universe up to this point. – Michael Seifert Aug 11 '17 at 13:41
• check if the expression under square root is positive. The same applies to the argument of the logarithm, it should be better positive, and check of course if all your variables have defined values. Finally, why using Mathematica? Matlab makes better plots without requiring a long list of plot attributes for make it nice. – Frederic Thomas Aug 11 '17 at 13:45
• Then what graphics you want to plot ? I think it's easier if you have $a(t)$ instead of $t(a)$. – Cham Aug 11 '17 at 14:33
• For $\rho = \mathrm{constant}$, the general solution for $k = -1, 0, +1$ can be cast into the following form : $$a(t) = a(0) \cosh{\omega t} \pm \sqrt{\displaystyle{a^2(0) - k \, \ell_{\Lambda}^2}} \: \sinh{\omega t},$$ where $\omega = \sqrt{\frac{\Lambda}{3}}$, $\ell_{\Lambda} \equiv \frac{1}{\omega}$ and $\rho = \frac{\Lambda}{8 \pi G}$. – Cham Aug 11 '17 at 14:38

The equation you want to solve is the first equation of Friedmann-Lemaitre : $$\tag{1} \dot{a}^2 \equiv \Big( \frac{d a}{d t} \Big)^2 = \frac{8 \pi G}{3} \, \rho \, a^2 - k,$$ where $k = -1, 0, 1$ for open, euclidian and closed space, respectively. You first need to assume some state of matter : $\rho \equiv f(a)$. For example $\rho \propto a^{- 3}$ for dustlike matter or $\rho \propto a^{- 4}$ for pure radiation. If you impose $\rho = \mathrm{constant}$, then this is like imposing an empty space with a cosmological constant only. It is then preferable to write $$\tag{2} \rho = \frac{\Lambda}{8 \pi G}.$$ The solution to the following differential equation is well known : $$\tag{3} \Big( \frac{d a}{d t} \Big)^2 = \frac{\Lambda}{3} \, a^2 - k,$$ with solution $$\tag{4} a(t) = a_0 \cosh{\omega t} \pm \sqrt{a_0^2 - \frac{3 k}{\Lambda}} \, \sinh{\omega t},$$ where $$\tag{5} \omega \equiv \sqrt{\frac{\Lambda}{3}}.$$ Since the scale factor is defined up to an arbitrary constant factor, you can set $a_0 = \sqrt{\frac{3}{\Lambda}} \equiv \omega^{-1}$ if $k = 1$ or $k = 0$, and $a_0 = 0$ if $k = -1$, so that (4) becomes \begin{align} a(t) = \omega^{-1} \cosh{\omega \, t}, \qquad \text{if $k = 1$}. \tag{6}\\[12pt] a(t) = \omega^{-1} e^{\pm \, \omega \, t}, \qquad \text{if $k = 0$}. \tag{7}\\[12pt] a(t) = \omega^{-1} \sinh{\omega \, t}, \qquad \text{if $k = -1$}. \tag{8} \end{align} These solutions describe an empty deSitter universe, with a positive cosmological constant, which can be viewed as a vacuum dominated universe.